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?

  • 8
    $\begingroup$ Your polynomial is very poorly behaved, try plotting it to see it for yourself. There is usually no good way to work around that. $\endgroup$
    – Kirill
    Feb 5 '19 at 10:17
  • 1
    $\begingroup$ 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. $\endgroup$
    – Bort
    Feb 5 '19 at 17:19
  • 7
    $\begingroup$ 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). $\endgroup$
    – GertVdE
    Feb 6 '19 at 19:25
  • 1
    $\begingroup$ 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. $\endgroup$
    – 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.


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)


    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
            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):

        last_v = v
        last_sign = sign


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.