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 |