1
$\begingroup$

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?

$\endgroup$
  • 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 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 at 17:19
  • 6
    $\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 at 19:25
  • $\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 at 8:41
4
$\begingroup$

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()
$\endgroup$

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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