"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "mpmath/calculus/polynomials.py" between
mpmath-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, error=False) def polyroots(ctx, coeffs, maxsteps=50, cleanup=True, extraprec=10,
: error=False):
""" """
Computes all roots (real or complex) of a given polynomial. The roots are Computes all roots (real or complex) of a given polynomial.
returned as a sorted list, where real roots appear first followed by
complex conjugate roots as adjacent elements. The polynomial should be
given as a list of coefficients, in the format used by :func:`~mpmath.polyva
l`.
The leading coefficient must be nonzero.
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 typicall The roots are computed to the current working precision accuracy. If this
y 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
answer with too low `extraprec`.
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 number. 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 real 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 []
orig = ctx.prec tol = +ctx.eps
weps = +ctx.eps with ctx.extraprec(extraprec):
try:
ctx.prec += 10
tol = ctx.eps * 128
deg = len(coeffs) - 1 deg = len(coeffs) - 1
# Must be monic # Must be monic
lead = ctx.convert(coeffs[0]) lead = ctx.convert(coeffs[0])
if lead == 1: if lead == 1:
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):
if not abs(err[i]) < tol: 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])) < weps: 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])) < weps: 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)))
finally:
ctx.prec = orig
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

Home  |  About  |  All  |  Newest  |  Fossies Dox  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTPS