A sympy-based newton polynomial constructor.

newtonPolynomial.py

"""
A sympy-based newton polynomial constructor. 

Given a set of function inputs and outputs, the newtonPolynomial function will construct an 
expression that for every input gives the corresponding output. For intermediate values,
the polynomial interpolates (giving varying results based on the shape of your input). 

This is useful when the result needs to be used outside of Python, because the 
expression can easily be copied. To convert the expression to a python function object,
use sympy.lambdify.
"""

from sympy import symbols, expand, lambdify
from sympy.mpmath import tan, pi

import math

from operator import mul
from functools import reduce

# sympy symbols
x = symbols('x')

# convenience functions
product = lambda *args: reduce(mul, *(list(args) + [1]))

# test data
labels = [(-3/2), (-3/4), 0, 3/4, 3/2]
points = [math.tan(v) for v in labels]

def apxnewtongenc(maxdiff, x, y):
	# taken from http://adorio-research.org/wordpress/?p=11165
    c = y[:]
    n = len(x) 
    for i in range(1, maxdiff+1):
       for j in range(n-1, i-1, -1): 
         c[j] = (c[j] - c[j-1])/(x[j]-x[j-i])
    return c

def newtonPolynomial(xs, ys):
	# based on https://en.wikipedia.org/wiki/Newton_polynomial#Example

	# make table, take heads
	k = len(xs)
	heads = apxnewtongenc(k - 1, xs, ys)

	# 
	coeffs = [product((x - v for v in xs[:i])) for i in range(k)]

	# the generated sympy expression
	eq = sum(expand(h * c) for h, c in zip(heads, coeffs))

	return eq
 
if __name__ == "__main__":
    func = newtonPolynomial(labels, points)
    print(func)
    pyfunc = lambdify(x, func)
    for a, b in zip(labels, points):
    	assert(pyfunc(a) - b < 1e-6)