example.py

from decimal import Decimal, getcontext
from xmath import fma

def horner_double(polynomial, x):
    s = 0.0
    for c in polynomial:
        s = s*x + c
    return s

def horner_fma(polynomial, x):
    s = 0.0
    for c in polynomial:
        s = fma(s, x, c)
    return s

def _two_sum(a, b):
    x = a + b
    t = x - a
    y = (a - (x - t)) + (b - t)
    return x, y

def _two_product_fma(a, b):
    x = a*b
    y = fma(a, b, -x)
    return x, y

def horner_compensated(polynomial, x):
    # Based on "Compensated Horner Scheme"
    # by S. Graillat, Ph. Langlois, N. Louvet, 2005
    s = error = 0.0
    for c in polynomial:
        p, pi = _two_product_fma(s, x)
        s, sigma = _two_sum(p, c)
        error = error*x + (pi + sigma)
    return s + error

def horner_decimal(polynomial, x):
    polynomial, x = map(Decimal, polynomial), Decimal(x)
    s = Decimal(0.0)
    for c in polynomial:
        s = s*x + c
    return float(s)

getcontext().prec = 100

polynomial, x = (1.0, -78.0, 2717.0, -55770.0, 749463.0, -6926634.0,
    44990231.0, -206070150.0, 657206836.0, -1414014888.0, 1931559552.0,
    -1486442880.0, 479001600.0), 12.001 # Wilkinson's polynomial

for f in (horner_double, horner_fma, horner_compensated, horner_decimal):
    print(repr(f(polynomial, x)), '\t', f.__name__)