Finding the roots of a simple polynomial with mpmath

Question

I'm trying to use mpmath.polyroots to find the roots of a simple polynomial with integer coefficients x*(x-4)**3, which when expanded has the coefficient vector of [1, -12, 48, 64, 0]. The following code fails:

import mpmath
p = [  1, -12,  48, -64,   0]
print mpmath.polyroots(p,maxsteps=2000)

with the error:

Traceback (most recent call last):
  File "poly.py", line 3, in <module>
    print mpmath.polyroots(p,maxsteps=2000)
  File "/usr/local/lib/python2.7/dist-packages/mpmath/calculus/polynomials.py", line 188, in polyroots
    % maxsteps)
mpmath.libmp.libhyper.NoConvergence: Didn't converge in maxsteps=2000 steps.

Increasing the number of steps does not help. The expected answer is obviously [0,4,4,4].

Is mpmath unable to find the roots of a polynomial if there is some multiplicity? How can I resolve this?

Answer 1

The documentation mentions that increasing extraprec may be necessary to achieve convergence. This is the case for your example, due to root multiplicity.

>>> mpmath.polyroots(p, maxsteps=100, extraprec=110)
[mpf('0.0'), mpf('4.0'), mpf('4.0'), mpf('4.0')]

The extraprec parameter means the extra number of digits to be used in the process of calculation, compared the number of digits needed in the result (which is 15 by default, set globally in mpmath.mp.dps). Its default value is 10; increasing it to 110 achieves convergence in 100 steps. This actually takes less time than trying (and failing) to find the roots with maxsteps=2000 at the default extraprec level.