#!/usr/bin/python
# -*- coding: utf-8 -*-
"""Calculate the values needed to perform the method of differences.
If you have some polynomial like ax^4 + bx^3 + cx^2 + dx + e,
you can calculate its value at any point by taking the dot
product of the vector [a b c d e] with the vector [x^4 x^3 x^2
x 1]. You can calculate its value at a sequence of points by
multiplying [a b c d e] through a matrix consisting of such
vectors. Thus, its value at such a sequence of points is a
linear function of the coefficients, and you can calculate the
coefficients from the samples by inverting that matrix if
there are enough samples.
That’s probably obvious to some people and impenetrable to
others.
The method of finite differences uses a simpler method to calculate a
sequence of evenly-spaced samples. It uses a state machine with n
finite differences d[0:n]. It works like this:
"""
import itertools, sys
def tabulate(d):
while True:
yield d[0]
for jj in range(len(d)-1): d[jj] += d[jj+1]
"""
tabulate([1, 3, 2]) produces the sequence of perfect
squares. This is a sort of discrete analogue to Taylor series
approximation: d[0] is the current value of the function, d[1]
is sort of its first derivative, d[2] is sort of its second
derivative, and so on. Specifically d[1] is the mean value of
the first derivative from here to the next timestep, and d[2]
is the mean value of the derivative of d[1] over the same
distance.
That’s not a particularly helpful way to look at it, so here’s
another one. The output values are linear combinations of the
differences:
d[0]
d[0] + d[1]
d[0] + 2d[1] + d[2]
d[0] + 3d[1] + 3d[2] + d[3]
So actually they’re the differences multiplied through a
triangular matrix consisting of Pascal’s triangle, the
binomial coefficients. So now we have three different
representations of the same polynomial, all linear functions
of each other: its values at a set of sample points, its
coefficients, and its finite differences.
So what we’d really like is a way to transform the
coefficients into the finite differences. That can be done
just by inverting the Pascal’s triangle matrix and multiplying
the result with the matrix of [x^4 x^3 ...] vectors, but I’m a
little vague on how to invert matrices, really, and so I’ll
start with a more direct approach.
d[0] = e
d[0] + d[1] = a + b + c + d + e, so
d[1] = a + b + c + d
d[0] + 2d[1] + d[2] = 16a + 8b + 4c + 2d + e, so
2d[1] + d[2] = 16a + 8b + 4c + 2d, so
d[2] = 14a + 6b + 2c
d[0] + 3d[1] + 3d[2] + d[3] = 81a + 27b + 9c + 3d + e, so
3d[1] + 3d[2] + d[3] = 81a + 27b + 9c + 3d, so
3d[2] + d[3] = 78a + 24b + 6c, so
d[3] = 36a + 6b
God, I am loving this lower triangular matrix shit.
Anyway, it’s pretty clear how to proceed in this way. You have
two matrices with corresponding rows, representing these big
old linear equations, each corresponding to a point on the
polynomial; you do corresponding operations on them to turn
the left one into the identity matrix, and then the right one
gives you the matrix for converting coefficients into
differences. As derived above, the matrix for up to cubics is:
1 0 0 0
0 1 1 1
0 0 2 6
0 0 0 6
"""
cubicmatrix = [[1, 0, 0, 0],
[0, 1, 1, 1],
[0, 0, 2, 6],
[0, 0, 0, 6]]
def matrix_vector_multiply(matrix, vector):
return [sum(a*b for a, b in zip(vector, row))
for row in matrix]
def tabulate_cubic(coefficients, k):
diffs = matrix_vector_multiply(cubicmatrix, coefficients)
return list(itertools.islice(tabulate(diffs), k))
def tabulate_polynomial(coefficients, k):
return [evaluate_polynomial(coefficients, ii) for ii in range(k)]
def evaluate_polynomial(coefficients, xx):
return sum(c*xx**n for n, c in enumerate(coefficients))
def ok(a, b): assert a == b, (a, b)
ok(tabulate_cubic( [1, 0, 0, 0], 10),
tabulate_polynomial([1, 0, 0, 0], 10))
ok(tabulate_cubic( [0, 1, 0, 0], 10),
tabulate_polynomial([0, 1, 0, 0], 10))
ok(tabulate_cubic( [0, 0, 1, 0], 10),
tabulate_polynomial([0, 0, 1, 0], 10))
ok(tabulate_cubic( [0, 0, 0, 1], 10),
tabulate_polynomial([0, 0, 0, 1], 10))
ok(tabulate_cubic( [30, -2, 1, -5], 10),
tabulate_polynomial([30, -2, 1, -5], 10))
"""
So to calculate that matrix, first we need the two matrices
whose rows represent the linear functions being equated.
"""
def pascaltri(n):
row = [1] + [0] * (n-1)
for ii in range(n):
yield row
row = [1] + [row[ii] + row[ii+1]
for ii in range(len(row)-1)]
def powermatrix(n):
return [[ii**jj for jj in range(n)]
for ii in range(n)]
"""
And then we need a function to do the row operations to solve
the simultaneous linear equations for the pascaltri side,
which happens to be lower triangular.
It also has the special property that its diagonal is already
all-ones, which means we can avoid any division. In turn, this ensures
that the result will be entirely integers if the inputs are entirely
integers, since integers are closed under addition, subtraction, and
multiplication.
"""
def solve_lower_triangular(lowertri, other):
n = len(lowertri)
assert n == len(other)
assert all(len(row) == n for row in lowertri)
assert all(len(row) == n for row in other)
for rownum in range(n):
assert lowertri[rownum][rownum] == 1
for colnum in range(rownum):
# The number to eliminate:
a = lowertri[rownum][colnum]
# Do the necessary row operations:
lowertri[rownum] = [lowertri[rownum][ii] -
a * lowertri[colnum][ii]
for ii in range(n)]
other[rownum] = [other[rownum][ii] -
a * other[colnum][ii]
for ii in range(n)]
assert lowertri[rownum][colnum] == 0
return other
def coefficients_to_diffs_matrix(n):
differences = list(pascaltri(n))
coefficients = powermatrix(n)
return solve_lower_triangular(differences, coefficients)
assert coefficients_to_diffs_matrix(4) == cubicmatrix
"""
Here’s the resulting matrix for 11th-order polynomials:
1 0 0 0 0 0 0
0 0 0 0 0
0 1 1 1 1 1 1
1 1 1 1 1
0 0 2 6 14 30 62
126 254 510 1022 2046
0 0 0 6 36 150 540
1806 5796 18150 55980 171006
0 0 0 0 24 240 1560
8400 40824 186480 818520 3498000
0 0 0 0 0 120 1800
16800 126000 834120 5103000 29607600
0 0 0 0 0 0 720
15120 191520 1905120 16435440 129230640
0 0 0 0 0 0 0
5040 141120 2328480 29635200 322494480
0 0 0 0 0 0 0
0 40320 1451520 30240000 479001600
0 0 0 0 0 0 0
0 0 362880 16329600 419126400
0 0 0 0 0 0 0
0 0 0 3628800 199584000
0 0 0 0 0 0 0
0 0 0 0 39916800
And then it remains only to compute the differences and apply
them.
(The real usefulness here, if any, is not in applying them ---
after all, that merely tabulates the values of the polynomial,
which you can do with a one-liner --- but in calculating them
so that the polynomial can be applied elsewhere.)
"""
def show_polynomial(coefficients):
rv = []
first = True
terms = list(enumerate(coefficients))
for ii, coefficient in reversed(terms):
if coefficient == 0:
continue
if not first:
rv.append("+" if coefficient > 0 else "-")
coefficient = abs(coefficient)
first = False
if ii == 0:
power = ""
elif ii == 1:
power = "x"
else:
power = "x^%d" % ii
if coefficient == 1 and power != '':
num = ""
elif coefficient == -1 and power != '':
num = "-"
else:
num = str(coefficient)
rv.append("%s%s" % (num, power))
if first: rv.append("0")
return ' '.join(rv)
ok(show_polynomial([0, 3]), "3x")
ok(show_polynomial([0]), "0")
ok(show_polynomial([2, 3]), "3x + 2")
ok(show_polynomial([2, 1]), "x + 2")
ok(show_polynomial([1, 3]), "3x + 1")
ok(show_polynomial([0, 0, 1]), "x^2")
ok(show_polynomial([1, 0, 1]), "x^2 + 1")
ok(show_polynomial([-1, 0, 1]), "x^2 - 1")
ok(show_polynomial([-1, 0, -1]), "-x^2 - 1")
ok(show_polynomial([-1, -2, 1]), "x^2 - 2x - 1")
ok(show_polynomial([0, 0, -1, 1]), "x^3 - x^2")
ok(show_polynomial([-1]), "-1")
def demo(coefficients, n):
print "f(x) =", show_polynomial(coefficients)
matrix = coefficients_to_diffs_matrix(len(coefficients))
print "Matrix to calculate differences:"
for row in matrix:
for value in row: print "%10d" % value,
print
differences = matrix_vector_multiply(matrix, coefficients)
print "Differences:",
for ii in differences: print ii,
print
values = tabulate(differences)
for ii in range(n):
print ("f(%d) =" % ii), values.next()
usage_message = """Usage: %s n a b c d e
n is the number of values to tabulate.
a b c d e are the coefficients for a polynomial:
ax^4 + bx^3 + cx^2 + dx + e
but you can specify any number of coefficients, giving a
polynomial of any degree.
"""
if __name__ == '__main__':
if len(sys.argv) == 1:
sys.stderr.write(usage_message % sys.argv[0])
sys.exit(1)
coefficients = [int(arg) for arg in sys.argv[2:]]
coefficients.reverse()
demo(coefficients, int(sys.argv[1]))
--
To unsubscribe: http://lists.canonical.org/mailman/listinfo/kragen-hacks
