Precision error in Python

Question

So I am ultimately trying to use Horner's rule (http://mathworld.wolfram.com/HornersRule.html) to evaluate polynomials, and creating a function to evaluate the polynomial. Whatever, so my problem is with how I wrote the function; it works for easy polynomials like 3x^2 + 2x^1 + 5 and so on. But once you get to evaluating a polynomial with a floating point number (something crazy like 1.8953343e-20, etc. ) it loses it's precision. Because I am using this function to evaluate roots of a polynomial using Newton's Method (http://tutorial.math.lamar.edu/Classes/CalcI/NewtonsMethod.aspx), I need this to be precise, so it doesn't lose it's value through a small rounding error and whatnot.

I have already troubleshooted with two other people that the problem lies within the evaluatePoly() function, and not my other functions that evaluates Newton's Method. Also, I originally evaluated the polynomial normally (multiplying x to the degree, multiplying by constant, etc.) and it pulled out the correct answer. However, the assignment requires one to use Horner's rule for easier calculation.

This is my following code:

def evaluatePoly(poly, x_):
    """Evaluates the polynomial at x = x_ and returns the result as a floating
point number using Horner's rule"""
    #http://mathworld.wolfram.com/HornersRule.html  
    total = 0.
    polyRev = poly[::-1]
    for nn in polyRev:
        total = total * x_
        total = total + nn
    return total

Note: I have already tried setting nn, x_, (total * x_) as floats using float().

This is the output I am receiving:

Polynomial: 5040x^0 + 1602x^1 + 1127x^2 - 214x^3 - 75x^4 + 4x^5 + 1x^6

Derivative: 1602x^0 + 2254x^1 - 642x^2 - 300x^3 + 20x^4 + 6x^5

(6.9027369297630505, False)

Starting at 100.00, no root found, last estimate was 6.90, giving value f(6.90) = -6.366463e-12
(6.9027369297630505, False)

Starting at 10.00, no root found, last estimate was 6.90, giving value f(6.90) = -6.366463e-12

(-2.6575456505038764, False)

Starting at 0.00, no root found, last estimate was -2.66, giving value f(-2.66) = 8.839758e+03

(-8.106973924480215, False)

Starting at -10.00, no root found, last estimate was -8.11, giving value f(-8.11) = -1.364242e-11

(-8.106973924480215, False)

Starting at -100.00, no root found, last estimate was -8.11, giving value f(-8.11) = -1.364242e-11

This is the output I need:

Polynomial: 5040x^0 + 1602x^1 + 1127x^2 - 214x^3 - 75x^4 + 4x^5 + 1x^6

Derivative: 1602x^0 + 2254x^1 - 642x^2 - 300x^3 + 20x^4 + 6x^5

(6.9027369297630505, False)

Starting at 100.00, no root found, last estimate was 6.90,giving value f(6.90) = -2.91038e-11

(6.9027369297630505, False)

Starting at 10.00, no root found, last estimate was 6.90,giving value f(6.90) = -2.91038e-11

(-2.657545650503874, False)

Starting at 0.00, no root found, last estimate was -2.66,giving value f(-2.66) = 8.83976e+03

(-8.106973924480215, True)

Starting at -10.00, root found at x = -8.11, giving value f(-8.11)= 0.000000e+00

(-8.106973924480215, True)

Starting at -100.00, root found at x = -8.11, giving value f(-8.11)= 0.000000e+00

Note: Please ignore the tuples of the errored output; That is the result of my newton's method, where the first result is the root and the second result is indicating whether it is a root or not.

Answer 1

Evaluation of a polynomial close to a root requires, by construction of the problem, that large terms cancel to yield a small result. The size of the intermediate terms can be bounded in a worst-case sense by the evaluation of the polynomial that has all its coefficients set to the absolute values of the coefficients of the original polynomial. For the first root that gives

In [6]: x0=6.9027369297630505
In [7]: evaluatePoly(poly,x0)                                     
Out[7]: -6.366462912410498e-12

In [8]: evaluatePoly([abs(c) for c in poly],abs(x0))
Out[8]: 481315.82997756737

This value is a first estimate for the magnification factor of the floating point errors, multiplied with the machine epsilon 2.22e-16 this gives a bound on the accumulated floating point errors of any evaluation method of 1.07e-10 and indeed both evaluation methods give values comfortably below this bound, indicating that the root was found within the capabilities of the floating point format.

Looking at the graph of the evaluation around x0 one sees that the basic assumption of a smooth curve fails at that magnification, and the x-axis is crossed in a jump so that no better value for x0 can be found:

Answer 2

Try this:

def evaluatePoly(poly, x_):
    '''Evaluate the polynomial poly at x = x_ and return the result as a
    floating-point number using Horner's rule'''
    total= 0
    degree =0 
    for coef in poly:
        total += (x_**degree) * coef
        degree += 1