Coverage for /usr/lib/python3/dist-packages/sympy/polys/polyroots.py: 6%
654 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Algorithms for computing symbolic roots of polynomials. """
4import math
5from functools import reduce
7from sympy.core import S, I, pi
8from sympy.core.exprtools import factor_terms
9from sympy.core.function import _mexpand
10from sympy.core.logic import fuzzy_not
11from sympy.core.mul import expand_2arg, Mul
12from sympy.core.numbers import Rational, igcd, comp
13from sympy.core.power import Pow
14from sympy.core.relational import Eq
15from sympy.core.sorting import ordered
16from sympy.core.symbol import Dummy, Symbol, symbols
17from sympy.core.sympify import sympify
18from sympy.functions import exp, im, cos, acos, Piecewise
19from sympy.functions.elementary.miscellaneous import root, sqrt
20from sympy.ntheory import divisors, isprime, nextprime
21from sympy.polys.domains import EX
22from sympy.polys.polyerrors import (PolynomialError, GeneratorsNeeded,
23 DomainError, UnsolvableFactorError)
24from sympy.polys.polyquinticconst import PolyQuintic
25from sympy.polys.polytools import Poly, cancel, factor, gcd_list, discriminant
26from sympy.polys.rationaltools import together
27from sympy.polys.specialpolys import cyclotomic_poly
28from sympy.utilities import public
29from sympy.utilities.misc import filldedent
33z = Symbol('z') # importing from abc cause O to be lost as clashing symbol
36def roots_linear(f):
37 """Returns a list of roots of a linear polynomial."""
38 r = -f.nth(0)/f.nth(1)
39 dom = f.get_domain()
41 if not dom.is_Numerical:
42 if dom.is_Composite:
43 r = factor(r)
44 else:
45 from sympy.simplify.simplify import simplify
46 r = simplify(r)
48 return [r]
51def roots_quadratic(f):
52 """Returns a list of roots of a quadratic polynomial. If the domain is ZZ
53 then the roots will be sorted with negatives coming before positives.
54 The ordering will be the same for any numerical coefficients as long as
55 the assumptions tested are correct, otherwise the ordering will not be
56 sorted (but will be canonical).
57 """
59 a, b, c = f.all_coeffs()
60 dom = f.get_domain()
62 def _sqrt(d):
63 # remove squares from square root since both will be represented
64 # in the results; a similar thing is happening in roots() but
65 # must be duplicated here because not all quadratics are binomials
66 co = []
67 other = []
68 for di in Mul.make_args(d):
69 if di.is_Pow and di.exp.is_Integer and di.exp % 2 == 0:
70 co.append(Pow(di.base, di.exp//2))
71 else:
72 other.append(di)
73 if co:
74 d = Mul(*other)
75 co = Mul(*co)
76 return co*sqrt(d)
77 return sqrt(d)
79 def _simplify(expr):
80 if dom.is_Composite:
81 return factor(expr)
82 else:
83 from sympy.simplify.simplify import simplify
84 return simplify(expr)
86 if c is S.Zero:
87 r0, r1 = S.Zero, -b/a
89 if not dom.is_Numerical:
90 r1 = _simplify(r1)
91 elif r1.is_negative:
92 r0, r1 = r1, r0
93 elif b is S.Zero:
94 r = -c/a
95 if not dom.is_Numerical:
96 r = _simplify(r)
98 R = _sqrt(r)
99 r0 = -R
100 r1 = R
101 else:
102 d = b**2 - 4*a*c
103 A = 2*a
104 B = -b/A
106 if not dom.is_Numerical:
107 d = _simplify(d)
108 B = _simplify(B)
110 D = factor_terms(_sqrt(d)/A)
111 r0 = B - D
112 r1 = B + D
113 if a.is_negative:
114 r0, r1 = r1, r0
115 elif not dom.is_Numerical:
116 r0, r1 = [expand_2arg(i) for i in (r0, r1)]
118 return [r0, r1]
121def roots_cubic(f, trig=False):
122 """Returns a list of roots of a cubic polynomial.
124 References
125 ==========
126 [1] https://en.wikipedia.org/wiki/Cubic_function, General formula for roots,
127 (accessed November 17, 2014).
128 """
129 if trig:
130 a, b, c, d = f.all_coeffs()
131 p = (3*a*c - b**2)/(3*a**2)
132 q = (2*b**3 - 9*a*b*c + 27*a**2*d)/(27*a**3)
133 D = 18*a*b*c*d - 4*b**3*d + b**2*c**2 - 4*a*c**3 - 27*a**2*d**2
134 if (D > 0) == True:
135 rv = []
136 for k in range(3):
137 rv.append(2*sqrt(-p/3)*cos(acos(q/p*sqrt(-3/p)*Rational(3, 2))/3 - k*pi*Rational(2, 3)))
138 return [i - b/3/a for i in rv]
140 # a*x**3 + b*x**2 + c*x + d -> x**3 + a*x**2 + b*x + c
141 _, a, b, c = f.monic().all_coeffs()
143 if c is S.Zero:
144 x1, x2 = roots([1, a, b], multiple=True)
145 return [x1, S.Zero, x2]
147 # x**3 + a*x**2 + b*x + c -> u**3 + p*u + q
148 p = b - a**2/3
149 q = c - a*b/3 + 2*a**3/27
151 pon3 = p/3
152 aon3 = a/3
154 u1 = None
155 if p is S.Zero:
156 if q is S.Zero:
157 return [-aon3]*3
158 u1 = -root(q, 3) if q.is_positive else root(-q, 3)
159 elif q is S.Zero:
160 y1, y2 = roots([1, 0, p], multiple=True)
161 return [tmp - aon3 for tmp in [y1, S.Zero, y2]]
162 elif q.is_real and q.is_negative:
163 u1 = -root(-q/2 + sqrt(q**2/4 + pon3**3), 3)
165 coeff = I*sqrt(3)/2
166 if u1 is None:
167 u1 = S.One
168 u2 = Rational(-1, 2) + coeff
169 u3 = Rational(-1, 2) - coeff
170 b, c, d = a, b, c # a, b, c, d = S.One, a, b, c
171 D0 = b**2 - 3*c # b**2 - 3*a*c
172 D1 = 2*b**3 - 9*b*c + 27*d # 2*b**3 - 9*a*b*c + 27*a**2*d
173 C = root((D1 + sqrt(D1**2 - 4*D0**3))/2, 3)
174 return [-(b + uk*C + D0/C/uk)/3 for uk in [u1, u2, u3]] # -(b + uk*C + D0/C/uk)/3/a
176 u2 = u1*(Rational(-1, 2) + coeff)
177 u3 = u1*(Rational(-1, 2) - coeff)
179 if p is S.Zero:
180 return [u1 - aon3, u2 - aon3, u3 - aon3]
182 soln = [
183 -u1 + pon3/u1 - aon3,
184 -u2 + pon3/u2 - aon3,
185 -u3 + pon3/u3 - aon3
186 ]
188 return soln
190def _roots_quartic_euler(p, q, r, a):
191 """
192 Descartes-Euler solution of the quartic equation
194 Parameters
195 ==========
197 p, q, r: coefficients of ``x**4 + p*x**2 + q*x + r``
198 a: shift of the roots
200 Notes
201 =====
203 This is a helper function for ``roots_quartic``.
205 Look for solutions of the form ::
207 ``x1 = sqrt(R) - sqrt(A + B*sqrt(R))``
208 ``x2 = -sqrt(R) - sqrt(A - B*sqrt(R))``
209 ``x3 = -sqrt(R) + sqrt(A - B*sqrt(R))``
210 ``x4 = sqrt(R) + sqrt(A + B*sqrt(R))``
212 To satisfy the quartic equation one must have
213 ``p = -2*(R + A); q = -4*B*R; r = (R - A)**2 - B**2*R``
214 so that ``R`` must satisfy the Descartes-Euler resolvent equation
215 ``64*R**3 + 32*p*R**2 + (4*p**2 - 16*r)*R - q**2 = 0``
217 If the resolvent does not have a rational solution, return None;
218 in that case it is likely that the Ferrari method gives a simpler
219 solution.
221 Examples
222 ========
224 >>> from sympy import S
225 >>> from sympy.polys.polyroots import _roots_quartic_euler
226 >>> p, q, r = -S(64)/5, -S(512)/125, -S(1024)/3125
227 >>> _roots_quartic_euler(p, q, r, S(0))[0]
228 -sqrt(32*sqrt(5)/125 + 16/5) + 4*sqrt(5)/5
229 """
230 # solve the resolvent equation
231 x = Dummy('x')
232 eq = 64*x**3 + 32*p*x**2 + (4*p**2 - 16*r)*x - q**2
233 xsols = list(roots(Poly(eq, x), cubics=False).keys())
234 xsols = [sol for sol in xsols if sol.is_rational and sol.is_nonzero]
235 if not xsols:
236 return None
237 R = max(xsols)
238 c1 = sqrt(R)
239 B = -q*c1/(4*R)
240 A = -R - p/2
241 c2 = sqrt(A + B)
242 c3 = sqrt(A - B)
243 return [c1 - c2 - a, -c1 - c3 - a, -c1 + c3 - a, c1 + c2 - a]
246def roots_quartic(f):
247 r"""
248 Returns a list of roots of a quartic polynomial.
250 There are many references for solving quartic expressions available [1-5].
251 This reviewer has found that many of them require one to select from among
252 2 or more possible sets of solutions and that some solutions work when one
253 is searching for real roots but do not work when searching for complex roots
254 (though this is not always stated clearly). The following routine has been
255 tested and found to be correct for 0, 2 or 4 complex roots.
257 The quasisymmetric case solution [6] looks for quartics that have the form
258 `x**4 + A*x**3 + B*x**2 + C*x + D = 0` where `(C/A)**2 = D`.
260 Although no general solution that is always applicable for all
261 coefficients is known to this reviewer, certain conditions are tested
262 to determine the simplest 4 expressions that can be returned:
264 1) `f = c + a*(a**2/8 - b/2) == 0`
265 2) `g = d - a*(a*(3*a**2/256 - b/16) + c/4) = 0`
266 3) if `f != 0` and `g != 0` and `p = -d + a*c/4 - b**2/12` then
267 a) `p == 0`
268 b) `p != 0`
270 Examples
271 ========
273 >>> from sympy import Poly
274 >>> from sympy.polys.polyroots import roots_quartic
276 >>> r = roots_quartic(Poly('x**4-6*x**3+17*x**2-26*x+20'))
278 >>> # 4 complex roots: 1+-I*sqrt(3), 2+-I
279 >>> sorted(str(tmp.evalf(n=2)) for tmp in r)
280 ['1.0 + 1.7*I', '1.0 - 1.7*I', '2.0 + 1.0*I', '2.0 - 1.0*I']
282 References
283 ==========
285 1. http://mathforum.org/dr.math/faq/faq.cubic.equations.html
286 2. https://en.wikipedia.org/wiki/Quartic_function#Summary_of_Ferrari.27s_method
287 3. https://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormula.html
288 4. https://people.bath.ac.uk/masjhd/JHD-CA.pdf
289 5. http://www.albmath.org/files/Math_5713.pdf
290 6. https://web.archive.org/web/20171002081448/http://www.statemaster.com/encyclopedia/Quartic-equation
291 7. https://eqworld.ipmnet.ru/en/solutions/ae/ae0108.pdf
292 """
293 _, a, b, c, d = f.monic().all_coeffs()
295 if not d:
296 return [S.Zero] + roots([1, a, b, c], multiple=True)
297 elif (c/a)**2 == d:
298 x, m = f.gen, c/a
300 g = Poly(x**2 + a*x + b - 2*m, x)
302 z1, z2 = roots_quadratic(g)
304 h1 = Poly(x**2 - z1*x + m, x)
305 h2 = Poly(x**2 - z2*x + m, x)
307 r1 = roots_quadratic(h1)
308 r2 = roots_quadratic(h2)
310 return r1 + r2
311 else:
312 a2 = a**2
313 e = b - 3*a2/8
314 f = _mexpand(c + a*(a2/8 - b/2))
315 aon4 = a/4
316 g = _mexpand(d - aon4*(a*(3*a2/64 - b/4) + c))
318 if f.is_zero:
319 y1, y2 = [sqrt(tmp) for tmp in
320 roots([1, e, g], multiple=True)]
321 return [tmp - aon4 for tmp in [-y1, -y2, y1, y2]]
322 if g.is_zero:
323 y = [S.Zero] + roots([1, 0, e, f], multiple=True)
324 return [tmp - aon4 for tmp in y]
325 else:
326 # Descartes-Euler method, see [7]
327 sols = _roots_quartic_euler(e, f, g, aon4)
328 if sols:
329 return sols
330 # Ferrari method, see [1, 2]
331 p = -e**2/12 - g
332 q = -e**3/108 + e*g/3 - f**2/8
333 TH = Rational(1, 3)
335 def _ans(y):
336 w = sqrt(e + 2*y)
337 arg1 = 3*e + 2*y
338 arg2 = 2*f/w
339 ans = []
340 for s in [-1, 1]:
341 root = sqrt(-(arg1 + s*arg2))
342 for t in [-1, 1]:
343 ans.append((s*w - t*root)/2 - aon4)
344 return ans
346 # whether a Piecewise is returned or not
347 # depends on knowing p, so try to put
348 # in a simple form
349 p = _mexpand(p)
352 # p == 0 case
353 y1 = e*Rational(-5, 6) - q**TH
354 if p.is_zero:
355 return _ans(y1)
357 # if p != 0 then u below is not 0
358 root = sqrt(q**2/4 + p**3/27)
359 r = -q/2 + root # or -q/2 - root
360 u = r**TH # primary root of solve(x**3 - r, x)
361 y2 = e*Rational(-5, 6) + u - p/u/3
362 if fuzzy_not(p.is_zero):
363 return _ans(y2)
365 # sort it out once they know the values of the coefficients
366 return [Piecewise((a1, Eq(p, 0)), (a2, True))
367 for a1, a2 in zip(_ans(y1), _ans(y2))]
370def roots_binomial(f):
371 """Returns a list of roots of a binomial polynomial. If the domain is ZZ
372 then the roots will be sorted with negatives coming before positives.
373 The ordering will be the same for any numerical coefficients as long as
374 the assumptions tested are correct, otherwise the ordering will not be
375 sorted (but will be canonical).
376 """
377 n = f.degree()
379 a, b = f.nth(n), f.nth(0)
380 base = -cancel(b/a)
381 alpha = root(base, n)
383 if alpha.is_number:
384 alpha = alpha.expand(complex=True)
386 # define some parameters that will allow us to order the roots.
387 # If the domain is ZZ this is guaranteed to return roots sorted
388 # with reals before non-real roots and non-real sorted according
389 # to real part and imaginary part, e.g. -1, 1, -1 + I, 2 - I
390 neg = base.is_negative
391 even = n % 2 == 0
392 if neg:
393 if even == True and (base + 1).is_positive:
394 big = True
395 else:
396 big = False
398 # get the indices in the right order so the computed
399 # roots will be sorted when the domain is ZZ
400 ks = []
401 imax = n//2
402 if even:
403 ks.append(imax)
404 imax -= 1
405 if not neg:
406 ks.append(0)
407 for i in range(imax, 0, -1):
408 if neg:
409 ks.extend([i, -i])
410 else:
411 ks.extend([-i, i])
412 if neg:
413 ks.append(0)
414 if big:
415 for i in range(0, len(ks), 2):
416 pair = ks[i: i + 2]
417 pair = list(reversed(pair))
419 # compute the roots
420 roots, d = [], 2*I*pi/n
421 for k in ks:
422 zeta = exp(k*d).expand(complex=True)
423 roots.append((alpha*zeta).expand(power_base=False))
425 return roots
428def _inv_totient_estimate(m):
429 """
430 Find ``(L, U)`` such that ``L <= phi^-1(m) <= U``.
432 Examples
433 ========
435 >>> from sympy.polys.polyroots import _inv_totient_estimate
437 >>> _inv_totient_estimate(192)
438 (192, 840)
439 >>> _inv_totient_estimate(400)
440 (400, 1750)
442 """
443 primes = [ d + 1 for d in divisors(m) if isprime(d + 1) ]
445 a, b = 1, 1
447 for p in primes:
448 a *= p
449 b *= p - 1
451 L = m
452 U = int(math.ceil(m*(float(a)/b)))
454 P = p = 2
455 primes = []
457 while P <= U:
458 p = nextprime(p)
459 primes.append(p)
460 P *= p
462 P //= p
463 b = 1
465 for p in primes[:-1]:
466 b *= p - 1
468 U = int(math.ceil(m*(float(P)/b)))
470 return L, U
473def roots_cyclotomic(f, factor=False):
474 """Compute roots of cyclotomic polynomials. """
475 L, U = _inv_totient_estimate(f.degree())
477 for n in range(L, U + 1):
478 g = cyclotomic_poly(n, f.gen, polys=True)
480 if f.expr == g.expr:
481 break
482 else: # pragma: no cover
483 raise RuntimeError("failed to find index of a cyclotomic polynomial")
485 roots = []
487 if not factor:
488 # get the indices in the right order so the computed
489 # roots will be sorted
490 h = n//2
491 ks = [i for i in range(1, n + 1) if igcd(i, n) == 1]
492 ks.sort(key=lambda x: (x, -1) if x <= h else (abs(x - n), 1))
493 d = 2*I*pi/n
494 for k in reversed(ks):
495 roots.append(exp(k*d).expand(complex=True))
496 else:
497 g = Poly(f, extension=root(-1, n))
499 for h, _ in ordered(g.factor_list()[1]):
500 roots.append(-h.TC())
502 return roots
505def roots_quintic(f):
506 """
507 Calculate exact roots of a solvable irreducible quintic with rational coefficients.
508 Return an empty list if the quintic is reducible or not solvable.
509 """
510 result = []
512 coeff_5, coeff_4, p_, q_, r_, s_ = f.all_coeffs()
514 if not all(coeff.is_Rational for coeff in (coeff_5, coeff_4, p_, q_, r_, s_)):
515 return result
517 if coeff_5 != 1:
518 f = Poly(f / coeff_5)
519 _, coeff_4, p_, q_, r_, s_ = f.all_coeffs()
521 # Cancel coeff_4 to form x^5 + px^3 + qx^2 + rx + s
522 if coeff_4:
523 p = p_ - 2*coeff_4*coeff_4/5
524 q = q_ - 3*coeff_4*p_/5 + 4*coeff_4**3/25
525 r = r_ - 2*coeff_4*q_/5 + 3*coeff_4**2*p_/25 - 3*coeff_4**4/125
526 s = s_ - coeff_4*r_/5 + coeff_4**2*q_/25 - coeff_4**3*p_/125 + 4*coeff_4**5/3125
527 x = f.gen
528 f = Poly(x**5 + p*x**3 + q*x**2 + r*x + s)
529 else:
530 p, q, r, s = p_, q_, r_, s_
532 quintic = PolyQuintic(f)
534 # Eqn standardized. Algo for solving starts here
535 if not f.is_irreducible:
536 return result
537 f20 = quintic.f20
538 # Check if f20 has linear factors over domain Z
539 if f20.is_irreducible:
540 return result
541 # Now, we know that f is solvable
542 for _factor in f20.factor_list()[1]:
543 if _factor[0].is_linear:
544 theta = _factor[0].root(0)
545 break
546 d = discriminant(f)
547 delta = sqrt(d)
548 # zeta = a fifth root of unity
549 zeta1, zeta2, zeta3, zeta4 = quintic.zeta
550 T = quintic.T(theta, d)
551 tol = S(1e-10)
552 alpha = T[1] + T[2]*delta
553 alpha_bar = T[1] - T[2]*delta
554 beta = T[3] + T[4]*delta
555 beta_bar = T[3] - T[4]*delta
557 disc = alpha**2 - 4*beta
558 disc_bar = alpha_bar**2 - 4*beta_bar
560 l0 = quintic.l0(theta)
561 Stwo = S(2)
562 l1 = _quintic_simplify((-alpha + sqrt(disc)) / Stwo)
563 l4 = _quintic_simplify((-alpha - sqrt(disc)) / Stwo)
565 l2 = _quintic_simplify((-alpha_bar + sqrt(disc_bar)) / Stwo)
566 l3 = _quintic_simplify((-alpha_bar - sqrt(disc_bar)) / Stwo)
568 order = quintic.order(theta, d)
569 test = (order*delta.n()) - ( (l1.n() - l4.n())*(l2.n() - l3.n()) )
570 # Comparing floats
571 if not comp(test, 0, tol):
572 l2, l3 = l3, l2
574 # Now we have correct order of l's
575 R1 = l0 + l1*zeta1 + l2*zeta2 + l3*zeta3 + l4*zeta4
576 R2 = l0 + l3*zeta1 + l1*zeta2 + l4*zeta3 + l2*zeta4
577 R3 = l0 + l2*zeta1 + l4*zeta2 + l1*zeta3 + l3*zeta4
578 R4 = l0 + l4*zeta1 + l3*zeta2 + l2*zeta3 + l1*zeta4
580 Res = [None, [None]*5, [None]*5, [None]*5, [None]*5]
581 Res_n = [None, [None]*5, [None]*5, [None]*5, [None]*5]
583 # Simplifying improves performance a lot for exact expressions
584 R1 = _quintic_simplify(R1)
585 R2 = _quintic_simplify(R2)
586 R3 = _quintic_simplify(R3)
587 R4 = _quintic_simplify(R4)
589 # hard-coded results for [factor(i) for i in _vsolve(x**5 - a - I*b, x)]
590 x0 = z**(S(1)/5)
591 x1 = sqrt(2)
592 x2 = sqrt(5)
593 x3 = sqrt(5 - x2)
594 x4 = I*x2
595 x5 = x4 + I
596 x6 = I*x0/4
597 x7 = x1*sqrt(x2 + 5)
598 sol = [x0, -x6*(x1*x3 - x5), x6*(x1*x3 + x5), -x6*(x4 + x7 - I), x6*(-x4 + x7 + I)]
600 R1 = R1.as_real_imag()
601 R2 = R2.as_real_imag()
602 R3 = R3.as_real_imag()
603 R4 = R4.as_real_imag()
605 for i, s in enumerate(sol):
606 Res[1][i] = _quintic_simplify(s.xreplace({z: R1[0] + I*R1[1]}))
607 Res[2][i] = _quintic_simplify(s.xreplace({z: R2[0] + I*R2[1]}))
608 Res[3][i] = _quintic_simplify(s.xreplace({z: R3[0] + I*R3[1]}))
609 Res[4][i] = _quintic_simplify(s.xreplace({z: R4[0] + I*R4[1]}))
611 for i in range(1, 5):
612 for j in range(5):
613 Res_n[i][j] = Res[i][j].n()
614 Res[i][j] = _quintic_simplify(Res[i][j])
615 r1 = Res[1][0]
616 r1_n = Res_n[1][0]
618 for i in range(5):
619 if comp(im(r1_n*Res_n[4][i]), 0, tol):
620 r4 = Res[4][i]
621 break
623 # Now we have various Res values. Each will be a list of five
624 # values. We have to pick one r value from those five for each Res
625 u, v = quintic.uv(theta, d)
626 testplus = (u + v*delta*sqrt(5)).n()
627 testminus = (u - v*delta*sqrt(5)).n()
629 # Evaluated numbers suffixed with _n
630 # We will use evaluated numbers for calculation. Much faster.
631 r4_n = r4.n()
632 r2 = r3 = None
634 for i in range(5):
635 r2temp_n = Res_n[2][i]
636 for j in range(5):
637 # Again storing away the exact number and using
638 # evaluated numbers in computations
639 r3temp_n = Res_n[3][j]
640 if (comp((r1_n*r2temp_n**2 + r4_n*r3temp_n**2 - testplus).n(), 0, tol) and
641 comp((r3temp_n*r1_n**2 + r2temp_n*r4_n**2 - testminus).n(), 0, tol)):
642 r2 = Res[2][i]
643 r3 = Res[3][j]
644 break
645 if r2 is not None:
646 break
647 else:
648 return [] # fall back to normal solve
650 # Now, we have r's so we can get roots
651 x1 = (r1 + r2 + r3 + r4)/5
652 x2 = (r1*zeta4 + r2*zeta3 + r3*zeta2 + r4*zeta1)/5
653 x3 = (r1*zeta3 + r2*zeta1 + r3*zeta4 + r4*zeta2)/5
654 x4 = (r1*zeta2 + r2*zeta4 + r3*zeta1 + r4*zeta3)/5
655 x5 = (r1*zeta1 + r2*zeta2 + r3*zeta3 + r4*zeta4)/5
656 result = [x1, x2, x3, x4, x5]
658 # Now check if solutions are distinct
660 saw = set()
661 for r in result:
662 r = r.n(2)
663 if r in saw:
664 # Roots were identical. Abort, return []
665 # and fall back to usual solve
666 return []
667 saw.add(r)
669 # Restore to original equation where coeff_4 is nonzero
670 if coeff_4:
671 result = [x - coeff_4 / 5 for x in result]
672 return result
675def _quintic_simplify(expr):
676 from sympy.simplify.simplify import powsimp
677 expr = powsimp(expr)
678 expr = cancel(expr)
679 return together(expr)
682def _integer_basis(poly):
683 """Compute coefficient basis for a polynomial over integers.
685 Returns the integer ``div`` such that substituting ``x = div*y``
686 ``p(x) = m*q(y)`` where the coefficients of ``q`` are smaller
687 than those of ``p``.
689 For example ``x**5 + 512*x + 1024 = 0``
690 with ``div = 4`` becomes ``y**5 + 2*y + 1 = 0``
692 Returns the integer ``div`` or ``None`` if there is no possible scaling.
694 Examples
695 ========
697 >>> from sympy.polys import Poly
698 >>> from sympy.abc import x
699 >>> from sympy.polys.polyroots import _integer_basis
700 >>> p = Poly(x**5 + 512*x + 1024, x, domain='ZZ')
701 >>> _integer_basis(p)
702 4
703 """
704 monoms, coeffs = list(zip(*poly.terms()))
706 monoms, = list(zip(*monoms))
707 coeffs = list(map(abs, coeffs))
709 if coeffs[0] < coeffs[-1]:
710 coeffs = list(reversed(coeffs))
711 n = monoms[0]
712 monoms = [n - i for i in reversed(monoms)]
713 else:
714 return None
716 monoms = monoms[:-1]
717 coeffs = coeffs[:-1]
719 # Special case for two-term polynominals
720 if len(monoms) == 1:
721 r = Pow(coeffs[0], S.One/monoms[0])
722 if r.is_Integer:
723 return int(r)
724 else:
725 return None
727 divs = reversed(divisors(gcd_list(coeffs))[1:])
729 try:
730 div = next(divs)
731 except StopIteration:
732 return None
734 while True:
735 for monom, coeff in zip(monoms, coeffs):
736 if coeff % div**monom != 0:
737 try:
738 div = next(divs)
739 except StopIteration:
740 return None
741 else:
742 break
743 else:
744 return div
747def preprocess_roots(poly):
748 """Try to get rid of symbolic coefficients from ``poly``. """
749 coeff = S.One
751 poly_func = poly.func
752 try:
753 _, poly = poly.clear_denoms(convert=True)
754 except DomainError:
755 return coeff, poly
757 poly = poly.primitive()[1]
758 poly = poly.retract()
760 # TODO: This is fragile. Figure out how to make this independent of construct_domain().
761 if poly.get_domain().is_Poly and all(c.is_term for c in poly.rep.coeffs()):
762 poly = poly.inject()
764 strips = list(zip(*poly.monoms()))
765 gens = list(poly.gens[1:])
767 base, strips = strips[0], strips[1:]
769 for gen, strip in zip(list(gens), strips):
770 reverse = False
772 if strip[0] < strip[-1]:
773 strip = reversed(strip)
774 reverse = True
776 ratio = None
778 for a, b in zip(base, strip):
779 if not a and not b:
780 continue
781 elif not a or not b:
782 break
783 elif b % a != 0:
784 break
785 else:
786 _ratio = b // a
788 if ratio is None:
789 ratio = _ratio
790 elif ratio != _ratio:
791 break
792 else:
793 if reverse:
794 ratio = -ratio
796 poly = poly.eval(gen, 1)
797 coeff *= gen**(-ratio)
798 gens.remove(gen)
800 if gens:
801 poly = poly.eject(*gens)
803 if poly.is_univariate and poly.get_domain().is_ZZ:
804 basis = _integer_basis(poly)
806 if basis is not None:
807 n = poly.degree()
809 def func(k, coeff):
810 return coeff//basis**(n - k[0])
812 poly = poly.termwise(func)
813 coeff *= basis
815 if not isinstance(poly, poly_func):
816 poly = poly_func(poly)
817 return coeff, poly
820@public
821def roots(f, *gens,
822 auto=True,
823 cubics=True,
824 trig=False,
825 quartics=True,
826 quintics=False,
827 multiple=False,
828 filter=None,
829 predicate=None,
830 strict=False,
831 **flags):
832 """
833 Computes symbolic roots of a univariate polynomial.
835 Given a univariate polynomial f with symbolic coefficients (or
836 a list of the polynomial's coefficients), returns a dictionary
837 with its roots and their multiplicities.
839 Only roots expressible via radicals will be returned. To get
840 a complete set of roots use RootOf class or numerical methods
841 instead. By default cubic and quartic formulas are used in
842 the algorithm. To disable them because of unreadable output
843 set ``cubics=False`` or ``quartics=False`` respectively. If cubic
844 roots are real but are expressed in terms of complex numbers
845 (casus irreducibilis [1]) the ``trig`` flag can be set to True to
846 have the solutions returned in terms of cosine and inverse cosine
847 functions.
849 To get roots from a specific domain set the ``filter`` flag with
850 one of the following specifiers: Z, Q, R, I, C. By default all
851 roots are returned (this is equivalent to setting ``filter='C'``).
853 By default a dictionary is returned giving a compact result in
854 case of multiple roots. However to get a list containing all
855 those roots set the ``multiple`` flag to True; the list will
856 have identical roots appearing next to each other in the result.
857 (For a given Poly, the all_roots method will give the roots in
858 sorted numerical order.)
860 If the ``strict`` flag is True, ``UnsolvableFactorError`` will be
861 raised if the roots found are known to be incomplete (because
862 some roots are not expressible in radicals).
864 Examples
865 ========
867 >>> from sympy import Poly, roots, degree
868 >>> from sympy.abc import x, y
870 >>> roots(x**2 - 1, x)
871 {-1: 1, 1: 1}
873 >>> p = Poly(x**2-1, x)
874 >>> roots(p)
875 {-1: 1, 1: 1}
877 >>> p = Poly(x**2-y, x, y)
879 >>> roots(Poly(p, x))
880 {-sqrt(y): 1, sqrt(y): 1}
882 >>> roots(x**2 - y, x)
883 {-sqrt(y): 1, sqrt(y): 1}
885 >>> roots([1, 0, -1])
886 {-1: 1, 1: 1}
888 ``roots`` will only return roots expressible in radicals. If
889 the given polynomial has some or all of its roots inexpressible in
890 radicals, the result of ``roots`` will be incomplete or empty
891 respectively.
893 Example where result is incomplete:
895 >>> roots((x-1)*(x**5-x+1), x)
896 {1: 1}
898 In this case, the polynomial has an unsolvable quintic factor
899 whose roots cannot be expressed by radicals. The polynomial has a
900 rational root (due to the factor `(x-1)`), which is returned since
901 ``roots`` always finds all rational roots.
903 Example where result is empty:
905 >>> roots(x**7-3*x**2+1, x)
906 {}
908 Here, the polynomial has no roots expressible in radicals, so
909 ``roots`` returns an empty dictionary.
911 The result produced by ``roots`` is complete if and only if the
912 sum of the multiplicity of each root is equal to the degree of
913 the polynomial. If strict=True, UnsolvableFactorError will be
914 raised if the result is incomplete.
916 The result can be be checked for completeness as follows:
918 >>> f = x**3-2*x**2+1
919 >>> sum(roots(f, x).values()) == degree(f, x)
920 True
921 >>> f = (x-1)*(x**5-x+1)
922 >>> sum(roots(f, x).values()) == degree(f, x)
923 False
926 References
927 ==========
929 .. [1] https://en.wikipedia.org/wiki/Cubic_equation#Trigonometric_and_hyperbolic_solutions
931 """
932 from sympy.polys.polytools import to_rational_coeffs
933 flags = dict(flags)
935 if isinstance(f, list):
936 if gens:
937 raise ValueError('redundant generators given')
939 x = Dummy('x')
941 poly, i = {}, len(f) - 1
943 for coeff in f:
944 poly[i], i = sympify(coeff), i - 1
946 f = Poly(poly, x, field=True)
947 else:
948 try:
949 F = Poly(f, *gens, **flags)
950 if not isinstance(f, Poly) and not F.gen.is_Symbol:
951 raise PolynomialError("generator must be a Symbol")
952 f = F
953 except GeneratorsNeeded:
954 if multiple:
955 return []
956 else:
957 return {}
958 else:
959 n = f.degree()
960 if f.length() == 2 and n > 2:
961 # check for foo**n in constant if dep is c*gen**m
962 con, dep = f.as_expr().as_independent(*f.gens)
963 fcon = -(-con).factor()
964 if fcon != con:
965 con = fcon
966 bases = []
967 for i in Mul.make_args(con):
968 if i.is_Pow:
969 b, e = i.as_base_exp()
970 if e.is_Integer and b.is_Add:
971 bases.append((b, Dummy(positive=True)))
972 if bases:
973 rv = roots(Poly((dep + con).xreplace(dict(bases)),
974 *f.gens), *F.gens,
975 auto=auto,
976 cubics=cubics,
977 trig=trig,
978 quartics=quartics,
979 quintics=quintics,
980 multiple=multiple,
981 filter=filter,
982 predicate=predicate,
983 **flags)
984 return {factor_terms(k.xreplace(
985 {v: k for k, v in bases})
986 ): v for k, v in rv.items()}
988 if f.is_multivariate:
989 raise PolynomialError('multivariate polynomials are not supported')
991 def _update_dict(result, zeros, currentroot, k):
992 if currentroot == S.Zero:
993 if S.Zero in zeros:
994 zeros[S.Zero] += k
995 else:
996 zeros[S.Zero] = k
997 if currentroot in result:
998 result[currentroot] += k
999 else:
1000 result[currentroot] = k
1002 def _try_decompose(f):
1003 """Find roots using functional decomposition. """
1004 factors, roots = f.decompose(), []
1006 for currentroot in _try_heuristics(factors[0]):
1007 roots.append(currentroot)
1009 for currentfactor in factors[1:]:
1010 previous, roots = list(roots), []
1012 for currentroot in previous:
1013 g = currentfactor - Poly(currentroot, f.gen)
1015 for currentroot in _try_heuristics(g):
1016 roots.append(currentroot)
1018 return roots
1020 def _try_heuristics(f):
1021 """Find roots using formulas and some tricks. """
1022 if f.is_ground:
1023 return []
1024 if f.is_monomial:
1025 return [S.Zero]*f.degree()
1027 if f.length() == 2:
1028 if f.degree() == 1:
1029 return list(map(cancel, roots_linear(f)))
1030 else:
1031 return roots_binomial(f)
1033 result = []
1035 for i in [-1, 1]:
1036 if not f.eval(i):
1037 f = f.quo(Poly(f.gen - i, f.gen))
1038 result.append(i)
1039 break
1041 n = f.degree()
1043 if n == 1:
1044 result += list(map(cancel, roots_linear(f)))
1045 elif n == 2:
1046 result += list(map(cancel, roots_quadratic(f)))
1047 elif f.is_cyclotomic:
1048 result += roots_cyclotomic(f)
1049 elif n == 3 and cubics:
1050 result += roots_cubic(f, trig=trig)
1051 elif n == 4 and quartics:
1052 result += roots_quartic(f)
1053 elif n == 5 and quintics:
1054 result += roots_quintic(f)
1056 return result
1058 # Convert the generators to symbols
1059 dumgens = symbols('x:%d' % len(f.gens), cls=Dummy)
1060 f = f.per(f.rep, dumgens)
1062 (k,), f = f.terms_gcd()
1064 if not k:
1065 zeros = {}
1066 else:
1067 zeros = {S.Zero: k}
1069 coeff, f = preprocess_roots(f)
1071 if auto and f.get_domain().is_Ring:
1072 f = f.to_field()
1074 # Use EX instead of ZZ_I or QQ_I
1075 if f.get_domain().is_QQ_I:
1076 f = f.per(f.rep.convert(EX))
1078 rescale_x = None
1079 translate_x = None
1081 result = {}
1083 if not f.is_ground:
1084 dom = f.get_domain()
1085 if not dom.is_Exact and dom.is_Numerical:
1086 for r in f.nroots():
1087 _update_dict(result, zeros, r, 1)
1088 elif f.degree() == 1:
1089 _update_dict(result, zeros, roots_linear(f)[0], 1)
1090 elif f.length() == 2:
1091 roots_fun = roots_quadratic if f.degree() == 2 else roots_binomial
1092 for r in roots_fun(f):
1093 _update_dict(result, zeros, r, 1)
1094 else:
1095 _, factors = Poly(f.as_expr()).factor_list()
1096 if len(factors) == 1 and f.degree() == 2:
1097 for r in roots_quadratic(f):
1098 _update_dict(result, zeros, r, 1)
1099 else:
1100 if len(factors) == 1 and factors[0][1] == 1:
1101 if f.get_domain().is_EX:
1102 res = to_rational_coeffs(f)
1103 if res:
1104 if res[0] is None:
1105 translate_x, f = res[2:]
1106 else:
1107 rescale_x, f = res[1], res[-1]
1108 result = roots(f)
1109 if not result:
1110 for currentroot in _try_decompose(f):
1111 _update_dict(result, zeros, currentroot, 1)
1112 else:
1113 for r in _try_heuristics(f):
1114 _update_dict(result, zeros, r, 1)
1115 else:
1116 for currentroot in _try_decompose(f):
1117 _update_dict(result, zeros, currentroot, 1)
1118 else:
1119 for currentfactor, k in factors:
1120 for r in _try_heuristics(Poly(currentfactor, f.gen, field=True)):
1121 _update_dict(result, zeros, r, k)
1123 if coeff is not S.One:
1124 _result, result, = result, {}
1126 for currentroot, k in _result.items():
1127 result[coeff*currentroot] = k
1129 if filter not in [None, 'C']:
1130 handlers = {
1131 'Z': lambda r: r.is_Integer,
1132 'Q': lambda r: r.is_Rational,
1133 'R': lambda r: all(a.is_real for a in r.as_numer_denom()),
1134 'I': lambda r: r.is_imaginary,
1135 }
1137 try:
1138 query = handlers[filter]
1139 except KeyError:
1140 raise ValueError("Invalid filter: %s" % filter)
1142 for zero in dict(result).keys():
1143 if not query(zero):
1144 del result[zero]
1146 if predicate is not None:
1147 for zero in dict(result).keys():
1148 if not predicate(zero):
1149 del result[zero]
1150 if rescale_x:
1151 result1 = {}
1152 for k, v in result.items():
1153 result1[k*rescale_x] = v
1154 result = result1
1155 if translate_x:
1156 result1 = {}
1157 for k, v in result.items():
1158 result1[k + translate_x] = v
1159 result = result1
1161 # adding zero roots after non-trivial roots have been translated
1162 result.update(zeros)
1164 if strict and sum(result.values()) < f.degree():
1165 raise UnsolvableFactorError(filldedent('''
1166 Strict mode: some factors cannot be solved in radicals, so
1167 a complete list of solutions cannot be returned. Call
1168 roots with strict=False to get solutions expressible in
1169 radicals (if there are any).
1170 '''))
1172 if not multiple:
1173 return result
1174 else:
1175 zeros = []
1177 for zero in ordered(result):
1178 zeros.extend([zero]*result[zero])
1180 return zeros
1183def root_factors(f, *gens, filter=None, **args):
1184 """
1185 Returns all factors of a univariate polynomial.
1187 Examples
1188 ========
1190 >>> from sympy.abc import x, y
1191 >>> from sympy.polys.polyroots import root_factors
1193 >>> root_factors(x**2 - y, x)
1194 [x - sqrt(y), x + sqrt(y)]
1196 """
1197 args = dict(args)
1199 F = Poly(f, *gens, **args)
1201 if not F.is_Poly:
1202 return [f]
1204 if F.is_multivariate:
1205 raise ValueError('multivariate polynomials are not supported')
1207 x = F.gens[0]
1209 zeros = roots(F, filter=filter)
1211 if not zeros:
1212 factors = [F]
1213 else:
1214 factors, N = [], 0
1216 for r, n in ordered(zeros.items()):
1217 factors, N = factors + [Poly(x - r, x)]*n, N + n
1219 if N < F.degree():
1220 G = reduce(lambda p, q: p*q, factors)
1221 factors.append(F.quo(G))
1223 if not isinstance(f, Poly):
1224 factors = [ f.as_expr() for f in factors ]
1226 return factors