def test_roots_preprocessed():
    E, F, J, L = symbols("E,F,J,L")

    f = -21601054687500000000*E**8*J**8/L**16 + \
        508232812500000000*F*x*E**7*J**7/L**14 - \
        4269543750000000*E**6*F**2*J**6*x**2/L**12 + \
        16194716250000*E**5*F**3*J**5*x**3/L**10 - \
        27633173750*E**4*F**4*J**4*x**4/L**8 + \
        14840215*E**3*F**5*J**3*x**5/L**6 + \
        54794*E**2*F**6*J**2*x**6/(5*L**4) - \
        1153*E*J*F**7*x**7/(80*L**2) + \
        633*F**8*x**8/160000

    assert roots(f, x) == {}

    R1 = roots(f.evalf(), x, multiple=True)
    R2 = [-1304.88375606366, 97.1168816800648, 186.946430171876, 245.526792947065,
          503.441004174773, 791.549343830097, 1273.16678129348, 1850.10650616851]

    w = Wild('w')
    p = w*E*J/(F*L**2)

    assert len(R1) == len(R2)

    for r1, r2 in zip(R1, R2):
        match = r1.match(p)
        assert match is not None and abs(match[w] - r2) < 1e-10

def test_roots_inexact():
    R1 = sorted([ r.evalf() for r in roots(x**2 + x + 1,   x) ])
    R2 = sorted([ r         for r in roots(x**2 + x + 1.0, x) ])

    for r1, r2 in zip(R1, R2):
        assert abs(r1 - r2) < 1e-12

    f = x**4 + 3.0*sqrt(2.0)*x**3 - (78.0 + 24.0*sqrt(3.0))*x**2 + 144.0*(2*sqrt(3.0) + 9.0)

    R1 = sorted(roots(f, multiple=True))
    R2 = sorted([-12.7530479110482, -3.85012393732929, 4.89897948556636, 7.46155167569183])

    for r1, r2 in zip(R1, R2):
        assert abs(r1 - r2) < 1e-10

def test_roots_binomial():
    assert roots_binomial(Poly(5*x, x)) == [0]
    assert roots_binomial(Poly(5*x**4, x)) == [0, 0, 0, 0]
    assert roots_binomial(Poly(5*x + 2, x)) == [-Rational(2, 5)]

    A = 10**Rational(3, 4)/10

    assert roots_binomial(Poly(5*x**4 + 2, x)) == \
        [-A - A*I, -A + A*I, A - A*I, A + A*I]

    a1 = Symbol('a1', nonnegative=True)
    b1 = Symbol('b1', nonnegative=True)

    r0 = roots_quadratic(Poly(a1*x**2 + b1, x))
    r1 = roots_binomial(Poly(a1*x**2 + b1, x))

    assert powsimp(r0[0]) == powsimp(r1[0])
    assert powsimp(r0[1]) == powsimp(r1[1])
    for a, b, s, n in cartes((1, 2), (1, 2), (-1, 1), (2, 3, 4, 5)):
        if a == b and a != 1:  # a == b == 1 is sufficient
            continue
        p = Poly(a*x**n + s*b)
        ans = roots_binomial(p)
        assert ans == _nsort(ans)

    # issue 8813
    assert roots(Poly(2*x**3 - 16*y**3, x)) == {
        2*y*(-S(1)/2 - sqrt(3)*I/2): 1,
        2*y: 1,
        2*y*(-S(1)/2 + sqrt(3)*I/2): 1}

    def _indicial(self):
        list_coeff = self.annihilator.listofpoly
        R = self.annihilator.parent.base
        x = self.x
        s = R.zero
        y = R.one

        def _pole_degree(poly):
            root_all = roots(poly.rep, filter='Z')
            if 0 in root_all.keys():
                return root_all[0]
            else:
                return 0

        degree = [j.degree() for j in list_coeff]
        degree = max(degree)
        inf = 10 * (max(1, degree) + max(1, self.annihilator.order))

        deg = lambda q: inf if q.is_zero else _pole_degree(q)
        b = deg(list_coeff[0])
        print (b)

        for j in range(1, len(list_coeff)):
            b = min(b, deg(list_coeff[j]) - j)
            print(b)

        for i, j in enumerate(list_coeff):
            listofdmp = j.all_coeffs()
            degree = len(listofdmp) - 1
            if - i - b <= 0:
                s = s + listofdmp[degree - i - b] * y
            y *= x - i
        return roots(s.rep, filter='R').keys()

def test_roots2():
    """Just test that calculating these roots does not hang
    (final result is not checked)
    """
    a, b, c, d, x = symbols("a,b,c,d,x")

    f1 = x**2*c + (a/b) + x*c*d - a
    f2 = x**2*(a + b*(c-d)*a) + x*a*b*c/(b*d-d) + (a*d-c/d)

    assert roots(f1, x).values() == [1, 1]
    assert roots(f2, x).values() == [1, 1]

    (zz, yy, xx, zy, zx, yx, k) = symbols("zz,yy,xx,zy,zx,yx,k")

    e1 = (zz-k)*(yy-k)*(xx-k) + zy*yx*zx + zx-zy-yx
    e2 = (zz-k)*yx*yx + zx*(yy-k)*zx + zy*zy*(xx-k)

    assert roots(e1 - e2, k).values() == [1, 1, 1]

    def doit(self, **hints):
        if not hints.get('roots', True):
            return self

        _roots = roots(self.poly, multiple=True)

        if len(_roots) < self.poly.degree():
            return self
        else:
            return Add(*[self.fun(r) for r in _roots])

    def simplify(self, x):
        """simplify(self, x)

           Compute a simplified representation of the function using
           property number 4.

           x can be:

           - a symbol

           Examples
           ========

           >>> from sympy import DiracDelta
           >>> from sympy.abc import x, y

           >>> DiracDelta(x*y).simplify(x)
           DiracDelta(x)/Abs(y)
           >>> DiracDelta(x*y).simplify(y)
           DiracDelta(y)/Abs(x)

           >>> DiracDelta(x**2 + x - 2).simplify(x)
           DiracDelta(x - 1)/3 + DiracDelta(x + 2)/3

           See Also
           ========

           is_simple, Directdelta

        """
        from sympy.polys.polyroots import roots

        if not self.args[0].has(x) or (len(self.args) > 1 and self.args[1] != 0 ):
            return self
        try:
            argroots = roots(self.args[0], x)
            result = 0
            valid = True
            darg = abs(diff(self.args[0], x))
            for r, m in argroots.items():
                if r.is_real is not False and m == 1:
                    result += self.func(x - r)/darg.subs(x, r)
                else:
                    # don't handle non-real and if m != 1 then
                    # a polynomial will have a zero in the derivative (darg)
                    # at r
                    valid = False
                    break
            if valid:
                return result
        except PolynomialError:
            pass
        return self

def test_roots_slow():
    """Just test that calculating these roots does not hang. """
    a, b, c, d, x = symbols("a,b,c,d,x")

    f1 = x**2*c + (a/b) + x*c*d - a
    f2 = x**2*(a + b*(c - d)*a) + x*a*b*c/(b*d - d) + (a*d - c/d)

    assert list(roots(f1, x).values()) == [1, 1]
    assert list(roots(f2, x).values()) == [1, 1]

    (zz, yy, xx, zy, zx, yx, k) = symbols("zz,yy,xx,zy,zx,yx,k")

    e1 = (zz - k)*(yy - k)*(xx - k) + zy*yx*zx + zx - zy - yx
    e2 = (zz - k)*yx*yx + zx*(yy - k)*zx + zy*zy*(xx - k)

    assert list(roots(e1 - e2, k).values()) == [1, 1, 1]

    f = x**3 + 2*x**2 + 8
    R = list(roots(f).keys())

    assert not any(i for i in [f.subs(x, ri).n(chop=True) for ri in R])

    def simplify(self, x):
        """simplify(self, x)

           Compute a simplified representation of the function using
           property number 4.

           x can be:

           - a symbol

           Examples
           ========

           >>> from sympy import DiracDelta
           >>> from sympy.abc import x, y

           >>> DiracDelta(x*y).simplify(x)
           DiracDelta(x)/Abs(y)
           >>> DiracDelta(x*y).simplify(y)
           DiracDelta(y)/Abs(x)

           >>> DiracDelta(x**2 + x - 2).simplify(x)
           DiracDelta(x - 1)/3 + DiracDelta(x + 2)/3

           See Also
           ========

           is_simple, Directdelta

        """
        from sympy.polys.polyroots import roots

        if not self.args[0].has(x) or (len(self.args)>1 and self.args[1] != 0 ):
            return self
        try:
            argroots = roots(self.args[0], x, \
                                                     multiple=True)
            result = 0
            valid = True
            darg = diff(self.args[0], x)
            for r in argroots:
                #should I care about multiplicities of roots?
                if r.is_real and not darg.subs(x,r).is_zero:
                    result = result + DiracDelta(x - r)/abs(darg.subs(x,r))
                else:
                    valid = False
                    break
            if valid:
                return result
        except PolynomialError:
            pass
        return self

def test_roots_slow():
    """Just test that calculating these roots does not hang. """
    a, b, c, d, x = symbols("a,b,c,d,x")

    f1 = x**2*c + (a/b) + x*c*d - a
    f2 = x**2*(a + b*(c-d)*a) + x*a*b*c/(b*d-d) + (a*d-c/d)

    assert roots(f1, x).values() == [1, 1]
    assert roots(f2, x).values() == [1, 1]

    (zz, yy, xx, zy, zx, yx, k) = symbols("zz,yy,xx,zy,zx,yx,k")

    e1 = (zz-k)*(yy-k)*(xx-k) + zy*yx*zx + zx-zy-yx
    e2 = (zz-k)*yx*yx + zx*(yy-k)*zx + zy*zy*(xx-k)

    assert roots(e1 - e2, k).values() == [1, 1, 1]

    f = x**3 + 2*x**2 + 8
    R = roots(f).keys()

    assert f.subs(x, R[0]).simplify() == 0
    assert f.subs(x, R[1]).simplify() == 0
    assert f.subs(x, R[2]).simplify() == 0

def test_roots_cubic():
    assert roots_cubic(Poly(2*x**3, x)) == [0, 0, 0]
    assert roots_cubic(Poly(x**3 - 3*x**2 + 3*x - 1, x)) == [1, 1, 1]

    assert roots_cubic(Poly(x**3 + 1, x)) == \
        [-1, S.Half - I*sqrt(3)/2, S.Half + I*sqrt(3)/2]
    assert roots_cubic(Poly(2*x**3 - 3*x**2 - 3*x - 1, x))[0] == \
         S.Half + 3**Rational(1, 3)/2 + 3**Rational(2, 3)/2
    eq = -x**3 + 2*x**2 + 3*x - 2
    assert roots(eq, trig=True, multiple=True) == \
           roots_cubic(Poly(eq, x), trig=True) == [
        S(2)/3 + 2*sqrt(13)*cos(acos(8*sqrt(13)/169)/3)/3,
        -2*sqrt(13)*sin(-acos(8*sqrt(13)/169)/3 + pi/6)/3 + S(2)/3,
        -2*sqrt(13)*cos(-acos(8*sqrt(13)/169)/3 + pi/3)/3 + S(2)/3,
        ]

    def evalf(self, points, method='RK4'):
        """
        Finds numerical value of a holonomic function using numerical methods.
        (RK4 by default). A set of points (real or complex) must be provided
        which will be the path for the numerical integration.

        The path should be given as a list [x1, x2, ... xn]. The numerical
        values will be computed at each point in this order x1 --> x2 --> x3
        ... --> xn.

        Returns values of the function at x1, x2, ... xn in a list.

        Examples
        =======

        >>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
        >>> from sympy.polys.domains import ZZ, QQ
        >>> from sympy import symbols
        >>> x = symbols('x')
        >>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')
        >>> # a straight line on the real axis from (0 to 1)
        >>> r = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]

        # using Runge-Kutta 4th order on e^x from 0.1 to 1.
        # exact solution at 1 is 2.71828182845905
        >>> HolonomicFunction(Dx - 1, x, 0, [1]).evalf(r)
        [1.10517083333333, 1.22140257085069, 1.34985849706254, 1.49182424008069,
        1.64872063859684, 1.82211796209193, 2.01375162659678, 2.22553956329232,
        2.45960141378007, 2.71827974413517]

        # using Euler's method for the same
        >>> HolonomicFunction(Dx - 1, x, 0, [1]).evalf(r, method='Euler')
        [1.1, 1.21, 1.331, 1.4641, 1.61051, 1.771561, 1.9487171, 2.14358881,
        2.357947691, 2.5937424601]

        One can also observe that the value obtained using Runge-Kutta 4th order
        is much more accurate than Euler's method.
        """

        from sympy.holonomic.numerical import _evalf
        for i in roots(self.annihilator.listofpoly[-1].rep):
            if i == self.x0 or i in points:
                raise TypeError("Provided path contains a singularity")
        return _evalf(self, points, method=method)

def test_roots_mixed():
    f = -1936 - 5056*x - 7592*x**2 + 2704*x**3 - 49*x**4

    _re, _im = intervals(f, all=True)
    _nroots = nroots(f)
    _sroots = roots(f, multiple=True)

    _re = [ Interval(a, b) for (a, b), _ in _re ]
    _im = [ Interval(re(a), re(b))*Interval(im(a), im(b)) for (a, b), _ in _im ]

    _intervals = _re + _im
    _sroots = [ r.evalf() for r in _sroots ]

    _nroots = sorted(_nroots, key=lambda x: x.sort_key())
    _sroots = sorted(_sroots, key=lambda x: x.sort_key())

    for _roots in (_nroots, _sroots):
        for i, r in zip(_intervals, _roots):
            if r.is_real:
                assert r in i
            else:
                assert (re(r), im(r)) in i

def test_roots():
    assert roots(1, x) == {}
    assert roots(x, x) == {S.Zero: 1}
    assert roots(x**9, x) == {S.Zero: 9}
    assert roots(((x - 2)*(x + 3)*(x - 4)).expand(), x) == {-S(3): 1, S(2): 1, S(4): 1}

    assert roots(2*x + 1, x) == {-S.Half: 1}
    assert roots((2*x + 1)**2, x) == {-S.Half: 2}
    assert roots((2*x + 1)**5, x) == {-S.Half: 5}
    assert roots((2*x + 1)**10, x) == {-S.Half: 10}

    assert roots(x**4 - 1, x) == {I: 1, S.One: 1, -S.One: 1, -I: 1}
    assert roots((x**4 - 1)**2, x) == {I: 2, S.One: 2, -S.One: 2, -I: 2}

    assert roots(((2*x - 3)**2).expand(), x) == { Rational(3, 2): 2}
    assert roots(((2*x + 3)**2).expand(), x) == {-Rational(3, 2): 2}

    assert roots(((2*x - 3)**3).expand(), x) == { Rational(3, 2): 3}
    assert roots(((2*x + 3)**3).expand(), x) == {-Rational(3, 2): 3}

    assert roots(((2*x - 3)**5).expand(), x) == { Rational(3, 2): 5}
    assert roots(((2*x + 3)**5).expand(), x) == {-Rational(3, 2): 5}

    assert roots(((a*x - b)**5).expand(), x) == { b/a: 5}
    assert roots(((a*x + b)**5).expand(), x) == {-b/a: 5}

    assert roots(x**2 + (-a - 1)*x + a, x) == {a: 1, S.One: 1}

    assert roots(x**4 - 2*x**2 + 1, x) == {S.One: 2, -S.One: 2}

    assert roots(x**6 - 4*x**4 + 4*x**3 - x**2, x) == \
        {S.One: 2, -1 - sqrt(2): 1, S.Zero: 2, -1 + sqrt(2): 1}

    assert roots(x**8 - 1, x) == {
        sqrt(2)/2 + I*sqrt(2)/2: 1,
        sqrt(2)/2 - I*sqrt(2)/2: 1,
        -sqrt(2)/2 + I*sqrt(2)/2: 1,
        -sqrt(2)/2 - I*sqrt(2)/2: 1,
        S.One: 1, -S.One: 1, I: 1, -I: 1
    }

    f = -2016*x**2 - 5616*x**3 - 2056*x**4 + 3324*x**5 + 2176*x**6 - \
        224*x**7 - 384*x**8 - 64*x**9

    assert roots(f) == {S(0): 2, -S(2): 2, S(2): 1, -S(7)/2: 1, -S(3)/2: 1, -S(1)/2: 1, S(3)/2: 1}

    assert roots((a + b + c)*x - (a + b + c + d), x) == {(a + b + c + d)/(a + b + c): 1}

    assert roots(x**3 + x**2 - x + 1, x, cubics=False) == {}
    assert roots(((x - 2)*(
        x + 3)*(x - 4)).expand(), x, cubics=False) == {-S(3): 1, S(2): 1, S(4): 1}
    assert roots(((x - 2)*(x + 3)*(x - 4)*(x - 5)).expand(), x, cubics=False) == \
        {-S(3): 1, S(2): 1, S(4): 1, S(5): 1}
    assert roots(x**3 + 2*x**2 + 4*x + 8, x) == {-S(2): 1, -2*I: 1, 2*I: 1}
    assert roots(x**3 + 2*x**2 + 4*x + 8, x, cubics=True) == \
        {-2*I: 1, 2*I: 1, -S(2): 1}
    assert roots((x**2 - x)*(x**3 + 2*x**2 + 4*x + 8), x ) == \
        {S(1): 1, S(0): 1, -S(2): 1, -2*I: 1, 2*I: 1}

    r1_2, r1_3, r1_9, r4_9, r19_27 = [ Rational(*r)
        for r in ((1, 2), (1, 3), (1, 9), (4, 9), (19, 27)) ]

    U = -r1_2 - r1_2*I*3**r1_2
    V = -r1_2 + r1_2*I*3**r1_2
    W = (r19_27 + r1_9*33**r1_2)**r1_3

    assert roots(x**3 + x**2 - x + 1, x, cubics=True) == {
        -r1_3 - U*W - r4_9*(U*W)**(-1): 1,
        -r1_3 - V*W - r4_9*(V*W)**(-1): 1,
        -r1_3 - W - r4_9*(  W)**(-1): 1,
    }

    f = (x**2 + 2*x + 3).subs(x, 2*x**2 + 3*x).subs(x, 5*x - 4)

    r13_20, r1_20 = [ Rational(*r)
        for r in ((13, 20), (1, 20)) ]

    s2 = sqrt(2)
    assert roots(f, x) == {
        r13_20 + r1_20*sqrt(1 - 8*I*s2): 1,
        r13_20 - r1_20*sqrt(1 - 8*I*s2): 1,
        r13_20 + r1_20*sqrt(1 + 8*I*s2): 1,
        r13_20 - r1_20*sqrt(1 + 8*I*s2): 1,
    }

    f = x**4 + x**3 + x**2 + x + 1

    r1_4, r1_8, r5_8 = [ Rational(*r) for r in ((1, 4), (1, 8), (5, 8)) ]

    assert roots(f, x) == {
        -r1_4 + r1_4*5**r1_2 + I*(r5_8 + r1_8*5**r1_2)**r1_2: 1,
        -r1_4 + r1_4*5**r1_2 - I*(r5_8 + r1_8*5**r1_2)**r1_2: 1,
        -r1_4 - r1_4*5**r1_2 + I*(r5_8 - r1_8*5**r1_2)**r1_2: 1,
        -r1_4 - r1_4*5**r1_2 - I*(r5_8 - r1_8*5**r1_2)**r1_2: 1,
    }

    f = z**3 + (-2 - y)*z**2 + (1 + 2*y - 2*x**2)*z - y + 2*x**2

    assert roots(f, z) == {
        S.One: 1,
#.........这里部分代码省略.........

    def series(self, n=6, coefficient=False, order=True):
        """
        Finds the power series expansion of given holonomic function.

        Examples
        ========

        >>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
        >>> from sympy.polys.domains import ZZ, QQ
        >>> from sympy import symbols
        >>> x = symbols('x')
        >>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')

        >>> HolonomicFunction(Dx - 1, x, 0, [1]).series()  # e^x
        1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + O(x**6)

        >>> HolonomicFunction(Dx**2 + 1, x, 0, [0, 1]).series(n=8)  # sin(x)
        x - x**3/6 + x**5/120 - x**7/5040 + O(x**8)

        See Also
        ========

        HolonomicFunction.to_sequence
        """

        recurrence = self.to_sequence()
        l = len(recurrence.u0) - 1
        k = recurrence.recurrence.order
        x = self.x
        seq_dmp = recurrence.recurrence.listofpoly
        R = recurrence.recurrence.parent.base
        K = R.get_field()
        seq = []

        if 0 in roots(seq_dmp[-1].rep, filter='Z').keys():
            singular = True
        else:
            singular = False

        for i, j in enumerate(seq_dmp):
            seq.append(K.new(j.rep))

        sub = [-seq[i] / seq[k] for i in range(k)]
        sol = [i for i in recurrence.u0]

        if l + 1 >= n:
            pass
        else:
            # use the initial conditions to find the next term
            for i in range(l + 1 - k, n - k):
                coeff = S(0)
                for j in range(k):
                    if i + j >= 0:
                        coeff += DMFsubs(sub[j], i) * sol[i + j]
                sol.append(coeff)

        if coefficient:
            return sol

        ser = S(0)
        for i, j in enumerate(sol):
            ser += x**i * j
        if order:
            return ser + Order(x**n, x)
        else:
            return ser

    def to_sequence(self):
        """
        Finds the recurrence relation in power series expansion
        of the function.

        Examples
        ========

        >>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
        >>> from sympy.polys.domains import ZZ, QQ
        >>> from sympy import symbols
        >>> x = symbols('x')
        >>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')

        >>> HolonomicFunction(Dx - 1, x, 0, [1]).to_sequence()
        HolonomicSequence((-1) + (n + 1)Sn, n), u(0) = 1

        See Also
        ========

        HolonomicFunction.series

        References
        ==========

        hal.inria.fr/inria-00070025/document
        """

        dict1 = {}
        n = symbols('n', integer=True)
        dom = self.annihilator.parent.base.dom
        R, _ = RecurrenceOperators(dom.old_poly_ring(n), 'Sn')

        for i, j in enumerate(self.annihilator.listofpoly):
            listofdmp = j.all_coeffs()
            degree = len(listofdmp) - 1
            for k in range(degree + 1):
                coeff = listofdmp[degree - k]
                if coeff == 0:
                    continue
                if i - k in dict1:
                    dict1[i - k] += (coeff * rf(n - k + 1, i))
                else:
                    dict1[i - k] = (coeff * rf(n - k + 1, i))

        sol = []
        lower = min(dict1.keys())
        upper = max(dict1.keys())

        for j in range(lower, upper + 1):
            if j in dict1.keys():
                sol.append(dict1[j].subs(n, n - lower))
            else:
                sol.append(S(0))
        # recurrence relation
        sol = RecurrenceOperator(sol, R)

        if not self._have_init_cond:
            return HolonomicSequence(sol)
        if self.x0 != 0:
            return HolonomicSequence(sol)
        # computing the initial conditions for recurrence
        order = sol.order
        all_roots = roots(sol.listofpoly[-1].rep, filter='Z')
        all_roots = all_roots.keys()

        if all_roots:
            max_root = max(all_roots)
            if max_root >= 0:
                order += max_root + 1

        y0 = _extend_y0(self, order)
        u0 = []
        # u(n) = y^n(0)/factorial(n)
        for i, j in enumerate(y0):
            u0.append(j / factorial(i))

        return HolonomicSequence(sol, u0)

 def _pole_degree(poly):
     root_all = roots(poly.rep, filter='Z')
     if 0 in root_all.keys():
         return root_all[0]
     else:
         return 0

def test_roots():
    assert roots(1, x) == {}
    assert roots(x, x) == {S.Zero: 1}
    assert roots(x**9, x) == {S.Zero: 9}
    assert roots(((x-2)*(x+3)*(x-4)).expand(), x) == {-S(3): 1, S(2): 1, S(4): 1}

    assert roots(2*x+1, x) == {-S.Half: 1}
    assert roots((2*x+1)**2, x) == {-S.Half: 2}
    assert roots((2*x+1)**5, x) == {-S.Half: 5}
    assert roots((2*x+1)**10, x) == {-S.Half: 10}

    assert roots(x**4 - 1, x) == {I: 1, S.One: 1, -S.One: 1, -I: 1}
    assert roots((x**4 - 1)**2, x) == {I: 2, S.One: 2, -S.One: 2, -I: 2}

    assert roots(((2*x-3)**2).expand(), x) == { Rational(3,2): 2}
    assert roots(((2*x+3)**2).expand(), x) == {-Rational(3,2): 2}

    assert roots(((2*x-3)**3).expand(), x) == { Rational(3,2): 3}
    assert roots(((2*x+3)**3).expand(), x) == {-Rational(3,2): 3}

    assert roots(((2*x-3)**5).expand(), x) == { Rational(3,2): 5}
    assert roots(((2*x+3)**5).expand(), x) == {-Rational(3,2): 5}

    assert roots(((a*x-b)**5).expand(), x) == { b/a: 5}
    assert roots(((a*x+b)**5).expand(), x) == {-b/a: 5}

    assert roots(x**4-2*x**2+1, x) == {S.One: 2, -S.One: 2}

    assert roots(x**6-4*x**4+4*x**3-x**2, x) == \
        {S.One: 2, -1 - sqrt(2): 1, S.Zero: 2, -1 + sqrt(2): 1}

    assert roots(x**8-1, x) == {
         2**S.Half/2 + I*2**S.Half/2: 1,
         2**S.Half/2 - I*2**S.Half/2: 1,
        -2**S.Half/2 + I*2**S.Half/2: 1,
        -2**S.Half/2 - I*2**S.Half/2: 1,
        S.One: 1, -S.One: 1, I: 1, -I: 1
    }

    f = -2016*x**2 - 5616*x**3 - 2056*x**4 + 3324*x**5 + 2176*x**6 - 224*x**7 - 384*x**8 - 64*x**9

    assert roots(f) == {S(0): 2, -S(2): 2, S(2): 1, -S(7)/2: 1, -S(3)/2: 1, -S(1)/2: 1, S(3)/2: 1}

    assert roots((a+b+c)*x - (a+b+c+d), x) == {(a+b+c+d)/(a+b+c): 1}

    assert roots(x**3+x**2-x+1, x, cubics=False) == {}
    assert roots(((x-2)*(x+3)*(x-4)).expand(), x, cubics=False) == {-S(3): 1, S(2): 1, S(4): 1}
    assert roots(((x-2)*(x+3)*(x-4)*(x-5)).expand(), x, cubics=False) == \
            {-S(3): 1, S(2): 1, S(4): 1, S(5): 1}
    assert roots(x**3 + 2*x**2 + 4*x + 8, x) == {-S(2): 1, -2*I: 1, 2*I: 1}
    assert roots(x**3 + 2*x**2 + 4*x + 8, x, cubics=True) == \
                {-2*I: 1, 2*I: 1, -S(2): 1}
    assert roots((x**2 - x)*(x**3 + 2*x**2 + 4*x + 8), x ) == \
                {S(1): 1, S(0): 1, -S(2): 1, -2*I: 1, 2*I: 1}

    r1_2, r1_3, r1_9, r4_9, r19_27 = [ Rational(*r) \
        for r in ((1,2), (1,3), (1,9), (4,9), (19,27)) ]

    U = -r1_2 - r1_2*I*3**r1_2
    V = -r1_2 + r1_2*I*3**r1_2
    W = (r19_27 + r1_9*33**r1_2)**r1_3

    assert roots(x**3+x**2-x+1, x, cubics=True) == {
        -r1_3 - U*W - r4_9*(U*W)**(-1): 1,
        -r1_3 - V*W - r4_9*(V*W)**(-1): 1,
        -r1_3 -   W - r4_9*(  W)**(-1): 1,
    }

    f = (x**2+2*x+3).subs(x, 2*x**2 + 3*x).subs(x, 5*x-4)

    r1_2, r13_20, r1_100 = [ Rational(*r) \
        for r in ((1,2), (13,20), (1,100)) ]

    assert roots(f, x) == {
        r13_20 + r1_100*(25 - 200*I*2**r1_2)**r1_2: 1,
        r13_20 - r1_100*(25 - 200*I*2**r1_2)**r1_2: 1,
        r13_20 + r1_100*(25 + 200*I*2**r1_2)**r1_2: 1,
        r13_20 - r1_100*(25 + 200*I*2**r1_2)**r1_2: 1,
    }

    f = x**4 + x**3 + x**2 + x + 1

    r1_4, r1_8, r5_8 = [ Rational(*r) for r in ((1,4), (1,8), (5,8)) ]

    assert roots(f, x) == {
        -r1_4 + r1_4*5**r1_2 + I*(r5_8 + r1_8*5**r1_2)**r1_2: 1,
        -r1_4 + r1_4*5**r1_2 - I*(r5_8 + r1_8*5**r1_2)**r1_2: 1,
        -r1_4 - r1_4*5**r1_2 + I*(r5_8 - r1_8*5**r1_2)**r1_2: 1,
        -r1_4 - r1_4*5**r1_2 - I*(r5_8 - r1_8*5**r1_2)**r1_2: 1,
    }

    f = z**3 + (-2 - y)*z**2 + (1 + 2*y - 2*x**2)*z - y + 2*x**2

    assert roots(f, z) == {
        S.One: 1,
        S.Half + S.Half*y + S.Half*(1 - 2*y + y**2 + 8*x**2)**S.Half: 1,
        S.Half + S.Half*y - S.Half*(1 - 2*y + y**2 + 8*x**2)**S.Half: 1,
    }

    assert roots(a*b*c*x**3 + 2*x**2 + 4*x + 8, x, cubics=False) == {}
#.........这里部分代码省略.........

def test_issue_14522():
    eq = Poly(x**4 + x**3*(16 + 32*I) + x**2*(-285 + 386*I) + x*(-2824 - 448*I) - 2058 - 6053*I, x)
    roots_eq = roots(eq)
    assert all(eq(r) == 0 for r in roots_eq)

def test_roots_composite():
    assert len(roots(Poly(y**3 + y**2*sqrt(x) + y + x, y, composite=True))) == 3