## "Fossies" - the Fresh Open Source Software Archive

### Source code changes of the file "mpmath/calculus/polynomials.py" betweenmpmath-0.18.tar.gz and mpmath-0.19.tar.gz

About: mpmath is a Python library for arbitrary-precision floating-point arithmetic.

polynomials.py  (mpmath-0.18):polynomials.py  (mpmath-0.19)
skipping to change at line 47 skipping to change at line 47
for c in coeffs[1:]: for c in coeffs[1:]:
if derivative: if derivative:
q = p + x*q q = p + x*q
p = c + x*p p = c + x*p
if derivative: if derivative:
return p, q return p, q
else: else:
return p return p
@defun @defun
def polyroots(ctx, coeffs, maxsteps=50, cleanup=True, extraprec=10, def polyroots(ctx, coeffs, maxsteps=50, cleanup=True, extraprec=10,
error=False):
""" """
Computes all roots (real or complex) of a given polynomial. Computes all roots (real or complex) of a given polynomial.
With *error=True*, :func:`~mpmath.polyroots` returns a tuple *(roots, err)* The roots are returned as a sorted list, where real roots appear first
where followed by complex conjugate roots as adjacent elements. The polynomial
*err* is an estimate of the maximum error among the computed roots. should be given as a list of coefficients, in the format used by
:func:`~mpmath.polyval`. The leading coefficient must be nonzero.
With *error=True*, :func:`~mpmath.polyroots` returns a tuple *(roots, err)*
where *err* is an estimate of the maximum error among the computed roots.
**Examples** **Examples**
Finding the three real roots of `x^3 - x^2 - 14x + 24`:: Finding the three real roots of `x^3 - x^2 - 14x + 24`::
>>> from mpmath import * >>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True >>> mp.dps = 15; mp.pretty = True
>>> nprint(polyroots([1,-1,-14,24]), 4) >>> nprint(polyroots([1,-1,-14,24]), 4)
[-4.0, 2.0, 3.0] [-4.0, 2.0, 3.0]
skipping to change at line 100 skipping to change at line 102
... print(r) ... print(r)
... ...
1.0 1.0
(-0.8090169943749474241 + 0.58778525229247312917j) (-0.8090169943749474241 + 0.58778525229247312917j)
(-0.8090169943749474241 - 0.58778525229247312917j) (-0.8090169943749474241 - 0.58778525229247312917j)
(0.3090169943749474241 + 0.95105651629515357212j) (0.3090169943749474241 + 0.95105651629515357212j)
(0.3090169943749474241 - 0.95105651629515357212j) (0.3090169943749474241 - 0.95105651629515357212j)
**Precision and conditioning** **Precision and conditioning**
Provided there are no repeated roots, :func:`~mpmath.polyroots` can The roots are computed to the current working precision accuracy. If this
accuracy cannot be achieved in `maxsteps` steps, then a `NoConvergence`
compute all roots of an arbitrary polynomial to high precision:: exception is raised. The algorithm internally is using the current working
precision extended by `extraprec`. If `NoConvergence` was raised, that is
caused either by not having enough extra precision to achieve convergence
(in which case increasing `extraprec` should fix the problem) or too low
`maxsteps` (in which case increasing `maxsteps` should fix the problem), or
a combination of both.
The user should always do a convergence study with regards to `extraprec`
to ensure accurate results. It is possible to get convergence to a wrong
Provided there are no repeated roots, :func:`~mpmath.polyroots` can
typically compute all roots of an arbitrary polynomial to high precision::
>>> mp.dps = 60 >>> mp.dps = 60
>>> for r in polyroots([1, 0, -10, 0, 1]): >>> for r in polyroots([1, 0, -10, 0, 1]):
... print r ... print r
... ...
-3.14626436994197234232913506571557044551247712918732870123249 -3.14626436994197234232913506571557044551247712918732870123249
-0.317837245195782244725757617296174288373133378433432554879127 -0.317837245195782244725757617296174288373133378433432554879127
0.317837245195782244725757617296174288373133378433432554879127 0.317837245195782244725757617296174288373133378433432554879127
3.14626436994197234232913506571557044551247712918732870123249 3.14626436994197234232913506571557044551247712918732870123249
>>> >>>
skipping to change at line 126 skipping to change at line 141
**Algorithm** **Algorithm**
:func:`~mpmath.polyroots` implements the Durand-Kerner method [1], which :func:`~mpmath.polyroots` implements the Durand-Kerner method [1], which
uses complex arithmetic to locate all roots simultaneously. uses complex arithmetic to locate all roots simultaneously.
The Durand-Kerner method can be viewed as approximately performing The Durand-Kerner method can be viewed as approximately performing
simultaneous Newton iteration for all the roots. In particular, simultaneous Newton iteration for all the roots. In particular,
the convergence to simple roots is quadratic, just like Newton's the convergence to simple roots is quadratic, just like Newton's
method. method.
Although all roots are internally calculated using complex arithmetic, Although all roots are internally calculated using complex arithmetic, any
any root found to have an imaginary part smaller than the estimated root found to have an imaginary part smaller than the estimated numerical
numerical error is truncated to a real Real roots are placed error is truncated to a real number (small real parts are also chopped).
first in the returned list, sorted by value. The remaining complex Real roots are placed first in the returned list, sorted by value. The
roots are sorted by their parts so that conjugate roots end up remaining complex roots are sorted by their real parts so that conjugate
next to each other. roots end up next to each other.
**References** **References**
1. http://en.wikipedia.org/wiki/Durand-Kerner_method 1. http://en.wikipedia.org/wiki/Durand-Kerner_method
""" """
if len(coeffs) <= 1: if len(coeffs) <= 1:
if not coeffs or not coeffs[0]: if not coeffs or not coeffs[0]:
raise ValueError("Input to polyroots must not be the zero polynomial ") raise ValueError("Input to polyroots must not be the zero polynomial ")
# Constant polynomial with no roots # Constant polynomial with no roots
return [] return []
tol = +ctx.eps
with ctx.extraprec(extraprec):
tol =
deg = len(coeffs) - 1 deg = len(coeffs) - 1
# Must be monic # Must be monic
coeffs = [ctx.convert(c) for c in coeffs] coeffs = [ctx.convert(c) for c in coeffs]
else: else:
coeffs = [c/lead for c in coeffs] coeffs = [c/lead for c in coeffs]
f = lambda x: ctx.polyval(coeffs, x) f = lambda x: ctx.polyval(coeffs, x)
roots = [ctx.mpc((0.4+0.9j)**n) for n in xrange(deg)] roots = [ctx.mpc((0.4+0.9j)**n) for n in xrange(deg)]
err = [ctx.one for n in xrange(deg)] err = [ctx.one for n in xrange(deg)]
# Durand-Kerner iteration until convergence # Durand-Kerner iteration until convergence
for step in xrange(maxsteps): for step in xrange(maxsteps):
if abs(max(err)) < tol: if abs(max(err)) < tol:
break break
for i in xrange(deg): for i in xrange(deg):
p = roots[i]
p = roots[i] x = f(p)
x = f(p) for j in range(deg):
for j in range(deg): if i != j:
if i != j: try:
try: x /= (p-roots[j])
x /= (p-roots[j]) except ZeroDivisionError:
except ZeroDivisionError: continue
continue roots[i] = p - x
roots[i] = p - x err[i] = abs(x)
err[i] = abs(x) if abs(max(err)) >= tol:
# Remove small imaginary parts raise ctx.NoConvergence("Didn't converge in maxsteps=%d steps." \
% maxsteps)
# Remove small real or imaginary parts
if cleanup: if cleanup:
for i in xrange(deg): for i in xrange(deg):
if abs(ctx._im(roots[i])) < if abs(roots[i]) < tol:
roots[i] = ctx.zero
elif abs(ctx._im(roots[i])) < tol:
roots[i] = roots[i].real roots[i] = roots[i].real
elif abs(ctx._re(roots[i])) < : elif abs(ctx._re(roots[i])) < tol:
roots[i] = roots[i].imag * 1j roots[i] = roots[i].imag * 1j
roots.sort(key=lambda x: (abs(ctx._im(x)), ctx._re(x))) roots.sort(key=lambda x: (abs(ctx._im(x)), ctx._re(x)))
if error: if error:
err = max(err) err = max(err)
err = max(err, ctx.ldexp(1, -orig+1)) err = max(err, ctx.ldexp(1, -orig+1))
return [+r for r in roots], +err return [+r for r in roots], +err
else: else:
return [+r for r in roots] return [+r for r in roots]
End of changes. 10 change blocks.
41 lines changed or deleted 51 lines changed or added