Isolate Complex Roots of Polynomials¶

AUTHOR:

Carl Witty (2007-11-18): initial version

This is an implementation of complex root isolation. That is, given a polynomial with exact complex coefficients, we compute isolating intervals for the complex roots of the polynomial. (Polynomials with integer, rational, Gaussian rational, or algebraic coefficients are supported.)

We use a simple algorithm. First, we compute a squarefree decomposition of the input polynomial; the resulting polynomials have no multiple roots. Then, we find the roots numerically, using NumPy (at low precision) or Pari (at high precision). Then, we verify the roots using interval arithmetic.

EXAMPLES:

Sage

sage: x = polygen(ZZ)
sage: (x^5 - x - 1).roots(ring=CIF)
[(1.167303978261419?, 1),
 (-0.764884433600585? - 0.352471546031727?*I, 1),
 (-0.764884433600585? + 0.352471546031727?*I, 1),
 (0.181232444469876? - 1.083954101317711?*I, 1),
 (0.181232444469876? + 1.083954101317711?*I, 1)]

Python

>>> from sage.all import *
>>> x = polygen(ZZ)
>>> (x**Integer(5) - x - Integer(1)).roots(ring=CIF)
[(1.167303978261419?, 1),
 (-0.764884433600585? - 0.352471546031727?*I, 1),
 (-0.764884433600585? + 0.352471546031727?*I, 1),
 (0.181232444469876? - 1.083954101317711?*I, 1),
 (0.181232444469876? + 1.083954101317711?*I, 1)]

Sage Live

x = polygen(ZZ)
(x^5 - x - 1).roots(ring=CIF)

sage.rings.polynomial.complex_roots.complex_roots(p, skip_squarefree=False, retval='interval', min_prec=0)[source]¶

Compute the complex roots of a given polynomial with exact coefficients (integer, rational, Gaussian rational, and algebraic coefficients are supported). Returns a list of pairs of a root and its multiplicity.

Roots are returned as a ComplexIntervalFieldElement; each interval includes exactly one root, and the intervals are disjoint.

By default, the algorithm will do a squarefree decomposition to get squarefree polynomials. The skip_squarefree parameter lets you skip this step. (If this step is skipped, and the polynomial has a repeated root, then the algorithm will loop forever!)

You can specify retval=’interval’ (the default) to get roots as complex intervals. The other options are retval=’algebraic’ to get elements of QQbar, or retval=’algebraic_real’ to get only the real roots, and to get them as elements of AA.

EXAMPLES:

Sage

sage: from sage.rings.polynomial.complex_roots import complex_roots
sage: x = polygen(ZZ)
sage: complex_roots(x^5 - x - 1)
[(1.167303978261419?, 1),
 (-0.764884433600585? - 0.352471546031727?*I, 1),
 (-0.764884433600585? + 0.352471546031727?*I, 1),
 (0.181232444469876? - 1.083954101317711?*I, 1),
 (0.181232444469876? + 1.083954101317711?*I, 1)]
sage: v = complex_roots(x^2 + 27*x + 181)

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.complex_roots import complex_roots
>>> x = polygen(ZZ)
>>> complex_roots(x**Integer(5) - x - Integer(1))
[(1.167303978261419?, 1),
 (-0.764884433600585? - 0.352471546031727?*I, 1),
 (-0.764884433600585? + 0.352471546031727?*I, 1),
 (0.181232444469876? - 1.083954101317711?*I, 1),
 (0.181232444469876? + 1.083954101317711?*I, 1)]
>>> v = complex_roots(x**Integer(2) + Integer(27)*x + Integer(181))

Sage Live

from sage.rings.polynomial.complex_roots import complex_roots
x = polygen(ZZ)
complex_roots(x^5 - x - 1)
v = complex_roots(x^2 + 27*x + 181)

Unfortunately due to numerical noise there can be a small imaginary part to each root depending on CPU, compiler, etc, and that affects the printing order. So we verify the real part of each root and check that the imaginary part is small in both cases:

Sage

sage: v  # random
[(-14.61803398874990?..., 1), (-12.3819660112501...? + 0.?e-27*I, 1)]
sage: sorted((v[0][0].real(),v[1][0].real()))
[-14.61803398874989?, -12.3819660112501...?]
sage: v[0][0].imag().upper() < 1e25
True
sage: v[1][0].imag().upper() < 1e25
True

sage: K.<im> = QuadraticField(-1)
sage: eps = 1/2^100
sage: x = polygen(K)
sage: p = (x-1)*(x-1-eps)*(x-1+eps)*(x-1-eps*im)*(x-1+eps*im)

Python

>>> from sage.all import *
>>> v  # random
[(-14.61803398874990?..., 1), (-12.3819660112501...? + 0.?e-27*I, 1)]
>>> sorted((v[Integer(0)][Integer(0)].real(),v[Integer(1)][Integer(0)].real()))
[-14.61803398874989?, -12.3819660112501...?]
>>> v[Integer(0)][Integer(0)].imag().upper() < RealNumber('1e25')
True
>>> v[Integer(1)][Integer(0)].imag().upper() < RealNumber('1e25')
True

>>> K = QuadraticField(-Integer(1), names=('im',)); (im,) = K._first_ngens(1)
>>> eps = Integer(1)/Integer(2)**Integer(100)
>>> x = polygen(K)
>>> p = (x-Integer(1))*(x-Integer(1)-eps)*(x-Integer(1)+eps)*(x-Integer(1)-eps*im)*(x-Integer(1)+eps*im)

Sage Live

v  # random
sorted((v[0][0].real(),v[1][0].real()))
v[0][0].imag().upper() < 1e25
v[1][0].imag().upper() < 1e25
K.<im> = QuadraticField(-1)
eps = 1/2^100
x = polygen(K)
p = (x-1)*(x-1-eps)*(x-1+eps)*(x-1-eps*im)*(x-1+eps*im)

This polynomial actually has all-real coefficients, and is very, very close to (x-1)^5:

Sage

sage: [RR(QQ(a)) for a in list(p - (x-1)^5)]
[3.87259191484932e-121, -3.87259191484932e-121]
sage: rts = complex_roots(p)
sage: [ComplexIntervalField(10)(rt[0] - 1) for rt in rts]
[-7.8887?e-31, 0, 7.8887?e-31, -7.8887?e-31*I, 7.8887?e-31*I]

Python

>>> from sage.all import *
>>> [RR(QQ(a)) for a in list(p - (x-Integer(1))**Integer(5))]
[3.87259191484932e-121, -3.87259191484932e-121]
>>> rts = complex_roots(p)
>>> [ComplexIntervalField(Integer(10))(rt[Integer(0)] - Integer(1)) for rt in rts]
[-7.8887?e-31, 0, 7.8887?e-31, -7.8887?e-31*I, 7.8887?e-31*I]

Sage Live

[RR(QQ(a)) for a in list(p - (x-1)^5)]
rts = complex_roots(p)
[ComplexIntervalField(10)(rt[0] - 1) for rt in rts]

We can get roots either as intervals, or as elements of QQbar or AA.

Sage

sage: p = (x^2 + x - 1)
sage: p = p * p(x*im)
sage: p
-x^4 + (im - 1)*x^3 + im*x^2 + (-im - 1)*x + 1

Python

>>> from sage.all import *
>>> p = (x**Integer(2) + x - Integer(1))
>>> p = p * p(x*im)
>>> p
-x^4 + (im - 1)*x^3 + im*x^2 + (-im - 1)*x + 1

Sage Live

p = (x^2 + x - 1)
p = p * p(x*im)
p

Two of the roots have a zero real component; two have a zero imaginary component. These zero components will be found slightly inaccurately, and the exact values returned are very sensitive to the (non-portable) results of NumPy. So we post-process the roots for printing, to get predictable doctest results.

Sage

sage: def tiny(x):
....:     return x.contains_zero() and x.absolute_diameter() <  1e-14
sage: def smash(x):
....:     x = CIF(x[0]) # discard multiplicity
....:     if tiny(x.imag()): return x.real()
....:     if tiny(x.real()): return CIF(0, x.imag())
sage: rts = complex_roots(p); type(rts[0][0]), sorted(map(smash, rts))
(<class 'sage.rings.complex_interval.ComplexIntervalFieldElement'>,
 [-1.618033988749895?, -0.618033988749895?*I,
  1.618033988749895?*I, 0.618033988749895?])
sage: rts = complex_roots(p, retval='algebraic'); type(rts[0][0]), sorted(map(smash, rts))
(<class 'sage.rings.qqbar.AlgebraicNumber'>,
 [-1.618033988749895?, -0.618033988749895?*I,
  1.618033988749895?*I, 0.618033988749895?])
sage: rts = complex_roots(p, retval='algebraic_real'); type(rts[0][0]), rts
(<class 'sage.rings.qqbar.AlgebraicReal'>,
 [(-1.618033988749895?, 1), (0.618033988749895?, 1)])

Python

>>> from sage.all import *
>>> def tiny(x):
...     return x.contains_zero() and x.absolute_diameter() <  RealNumber('1e-14')
>>> def smash(x):
...     x = CIF(x[Integer(0)]) # discard multiplicity
...     if tiny(x.imag()): return x.real()
...     if tiny(x.real()): return CIF(Integer(0), x.imag())
>>> rts = complex_roots(p); type(rts[Integer(0)][Integer(0)]), sorted(map(smash, rts))
(<class 'sage.rings.complex_interval.ComplexIntervalFieldElement'>,
 [-1.618033988749895?, -0.618033988749895?*I,
  1.618033988749895?*I, 0.618033988749895?])
>>> rts = complex_roots(p, retval='algebraic'); type(rts[Integer(0)][Integer(0)]), sorted(map(smash, rts))
(<class 'sage.rings.qqbar.AlgebraicNumber'>,
 [-1.618033988749895?, -0.618033988749895?*I,
  1.618033988749895?*I, 0.618033988749895?])
>>> rts = complex_roots(p, retval='algebraic_real'); type(rts[Integer(0)][Integer(0)]), rts
(<class 'sage.rings.qqbar.AlgebraicReal'>,
 [(-1.618033988749895?, 1), (0.618033988749895?, 1)])

Sage Live

def tiny(x):
    return x.contains_zero() and x.absolute_diameter() <  1e-14
def smash(x):
    x = CIF(x[0]) # discard multiplicity
    if tiny(x.imag()): return x.real()
    if tiny(x.real()): return CIF(0, x.imag())
rts = complex_roots(p); type(rts[0][0]), sorted(map(smash, rts))
rts = complex_roots(p, retval='algebraic'); type(rts[0][0]), sorted(map(smash, rts))
rts = complex_roots(p, retval='algebraic_real'); type(rts[0][0]), rts

sage.rings.polynomial.complex_roots.interval_roots(p, rts, prec)[source]¶

We are given a squarefree polynomial p, a list of estimated roots, and a precision.

We attempt to verify that the estimated roots are in fact distinct roots of the polynomial, using interval arithmetic of precision prec. If we succeed, we return a list of intervals bounding the roots; if we fail, we return None.

EXAMPLES:

Sage

sage: x = polygen(ZZ)
sage: p = x^3 - 1
sage: rts = [CC.zeta(3)^i for i in range(0, 3)]
sage: from sage.rings.polynomial.complex_roots import interval_roots
sage: interval_roots(p, rts, 53)
[1, -0.500000000000000? + 0.866025403784439?*I,
    -0.500000000000000? - 0.866025403784439?*I]
sage: interval_roots(p, rts, 200)
[1, -0.500000000000000000000000000000000000000000000000000000000000?
      + 0.866025403784438646763723170752936183471402626905190314027904?*I,
    -0.500000000000000000000000000000000000000000000000000000000000?
      - 0.866025403784438646763723170752936183471402626905190314027904?*I]

Python

>>> from sage.all import *
>>> x = polygen(ZZ)
>>> p = x**Integer(3) - Integer(1)
>>> rts = [CC.zeta(Integer(3))**i for i in range(Integer(0), Integer(3))]
>>> from sage.rings.polynomial.complex_roots import interval_roots
>>> interval_roots(p, rts, Integer(53))
[1, -0.500000000000000? + 0.866025403784439?*I,
    -0.500000000000000? - 0.866025403784439?*I]
>>> interval_roots(p, rts, Integer(200))
[1, -0.500000000000000000000000000000000000000000000000000000000000?
      + 0.866025403784438646763723170752936183471402626905190314027904?*I,
    -0.500000000000000000000000000000000000000000000000000000000000?
      - 0.866025403784438646763723170752936183471402626905190314027904?*I]

Sage Live

x = polygen(ZZ)
p = x^3 - 1
rts = [CC.zeta(3)^i for i in range(0, 3)]
from sage.rings.polynomial.complex_roots import interval_roots
interval_roots(p, rts, 53)
interval_roots(p, rts, 200)

sage.rings.polynomial.complex_roots.intervals_disjoint(intvs)[source]¶

Given a list of complex intervals, check whether they are pairwise disjoint.

EXAMPLES:

Sage

sage: from sage.rings.polynomial.complex_roots import intervals_disjoint
sage: a = CIF(RIF(0, 3), 0)
sage: b = CIF(0, RIF(1, 3))
sage: c = CIF(RIF(1, 2), RIF(1, 2))
sage: d = CIF(RIF(2, 3), RIF(2, 3))
sage: intervals_disjoint([a,b,c,d])
False
sage: d2 = CIF(RIF(2, 3), RIF(2.001, 3))
sage: intervals_disjoint([a,b,c,d2])
True

Python

>>> from sage.all import *
>>> from sage.rings.polynomial.complex_roots import intervals_disjoint
>>> a = CIF(RIF(Integer(0), Integer(3)), Integer(0))
>>> b = CIF(Integer(0), RIF(Integer(1), Integer(3)))
>>> c = CIF(RIF(Integer(1), Integer(2)), RIF(Integer(1), Integer(2)))
>>> d = CIF(RIF(Integer(2), Integer(3)), RIF(Integer(2), Integer(3)))
>>> intervals_disjoint([a,b,c,d])
False
>>> d2 = CIF(RIF(Integer(2), Integer(3)), RIF(RealNumber('2.001'), Integer(3)))
>>> intervals_disjoint([a,b,c,d2])
True

Sage Live

from sage.rings.polynomial.complex_roots import intervals_disjoint
a = CIF(RIF(0, 3), 0)
b = CIF(0, RIF(1, 3))
c = CIF(RIF(1, 2), RIF(1, 2))
d = CIF(RIF(2, 3), RIF(2, 3))
intervals_disjoint([a,b,c,d])
d2 = CIF(RIF(2, 3), RIF(2.001, 3))
intervals_disjoint([a,b,c,d2])