# Find the roots of a complicated polynomial

For the polynomial

$$p(x) = (x-1)(x-2) \cdots (x-20) - 2^{-23}x^{19}\, ,$$

I tried to use fzero() in MATLAB, and I set the interval to be $$[0.5,1.5]\cdots [19.5,20.5]$$. It did work for $$[0.5, 1.5]$$ to about $$[7.5, 8.5]$$. However, when it goes to larger numbers it seems not to work, due to the sign.

Does anyone know how to fix that? Or some other methods?

• Your polynomial is very poorly behaved, try plotting it to see it for yourself. There is usually no good way to work around that. – Kirill Feb 5 '19 at 10:17
• If you rewrite your polynomial with a Horner form for evaluation and represent it via chebyshev functions in chebfun you can get something which is close to the real roots. – Bort Feb 5 '19 at 17:19
• This is Wilkinson's polynomial, an example of a polynomial where root finding is notoriousy difficult due to its bad conditioning (a small perturbation in the coefficient of $x^{19}$ leads to a "dispersion" of the roots $(1,2,3, \ldots, 30)$ to complex conjugate pairs. There is no solution in floating point arithmetic except increasing the precision (but then a smaller perturbation will cause the same issue). – GertVdE Feb 6 '19 at 19:25
• Just an addendum on conditioning of a problem and stability of an algorithm: if your problem is ill-conditioned, a good and stable algorithm will perform up to a certain level but you will always run into trouble, no matter what. Ill-conditioning is inherent to your problem, no algorithm can rescue you from that. – GertVdE Feb 7 '19 at 8:41

One technique is to use a library with arbitrary large integers and use fixed point arithmetic. Just scale up x by some integer factor, then compute your polynomial in fixed-point, keeping only the sign.

From two values, that you know positive and negative, you can then perform a binary search. At the very end, divide back by the fixed point factor and get the value you are looking for.

The below is still work-in-progress, but should give you the right idea. (Using vanilla python, because it supports large integers)

It finds the following roots. It has not been debugged thoroughly, but hopefully, this gives you an idea of a way to do it.

1.0
2.0
2.999999999999805
4.000000000261023
4.999999927551538
6.000006943952296
6.999697233936014
8.007267603450376
8.917250248517071
20.846908101482256


Finding the complex roots is left as an exercise to the reader :-)

# an arbitrary fixed-point precision factor
factor = 2 ** 80

def poly(x):
# this function computes 2^23*(x−1)(x−2)...(x−20) − x^19
# it has the same sign as (x−1)(x−2)...(x−20) − x^19 / 2^23

pos = 1
for i in range(20):
pos = pos * (x - (i + 1) * factor)
pos = (2 ** 23) * pos
# pos is (2 ** 23) * (x−1)(x−2)...(x−20) (factor^20)

neg = (x ** 19) * factor
# neg is x^19 (factor^19) (factor)

#print(pos,neg)

ret = pos - neg
# ret is now 2^23*(x−1)(x−2)...(x−20) (factor^20) - x^19 (factor^19) (factor)
# which is 2^23*(x−1)(x−2)...(x−20) (factor^20) - x^19 (factor^20)
# which is (factor^20)( 2^23*(x−1)(x−2)...(x−20) - x^19)
# which has the same sign as (2^23*(x−1)(x−2)...(x−20) - x^19)

return ret

def search(x,y):

# print("Searching between {} and {}".format(x/factor,y/factor))
while((y - x) > 1):
z = (x + y) // 2

if (poly(z) < 0) == (poly(x) < 0):
x = z
else:
y = z

print(y / factor) # here, we go back to floating point and incur loss of precision

def main():

poly(2 * factor)

last_v = None
last_sign = None

steps = 100000 # number of steps
lower = -10   # range from
upper = 30    # range to

for i in range(steps):
v = (((steps - i) * lower) + (i * upper)) * factor // steps

sign = poly(v) > 0

if (last_sign is not None) and (sign != last_sign):
search(last_v,v)

last_v = v
last_sign = sign

main()