barycentric_lagrange.py

from math import isinf

def product(xs):
    p = 1
    for x in xs:
        p *= x
    return p

# O(n) time per interpolation once the weights are computed for the sample points.
def interp(xs, ws, fs, x):
    y, l = 0, 0
    for xi, wi, fi in zip(xs, ws, fs):
        if x == xi: return fi
        li = wi/(x - xi)
        yi = li*fi
        if isinf(yi): return fi # Only check that should be needed, but Python raises div-by-zero FP exceptions.
        y += yi
        l += li
    return y/l

# O(n^2) time, but only depends on sample points, not on function values.
def weights(xs):
    return [1/product(xi - xj for j, xj in enumerate(xs) if i != j) for i, xi in enumerate(xs)]

xs = [-1, 0, 1]
ws = weights(xs)
r = 1e200
fs = [r, 0, r]
print([interp(xs, ws, fs, x) for x in (-3, -2, 1, 0, 1, 2, 3)])