Spline path computation using integer addition

A while back, I posted to the RepRap orgs with an idea for using spline paths instead of linear paths on an extruder.

Whether this is a good idea or not, I wanted to put down the algorithm for performing the integer spline computations.

Given:
X0, X1, dXn, dXn as the endpoints, and derivatives of a function x(t) for t=0 and t=n.

Desired:
A set of functions to compute x(t+1) give x(t), using only addition.

Solution:
Using algebra, and solving for the various functions needed:

Cubic spline basis:
x(t) = at^3 + bt^2 + ct + d
x'(t) = 3at^2 + 2bt + c

Endpoints and derivitives at endpoints yields the following:

a = ( n * dX0 + 2 * X0 - 2 * X1 + n * dXn ) / ( n^3 )
b = -( 2 * n * dX0 + 3 * X0 - 3 * X1 + n * dX1 ) / ( n^2 )
c = dX0
d = X0

To solve using only addition, functions p(t), q(t), and r(t) are defined as solutions to:

x(t+1) = x(t) + p(t)
p(t+1) = p(t) + q(t)
q(t+1) = q(t) + r(t)
r(t+1) = r(t)

The solution to this system is:

p(t) = 3at^2 + 3at + 2bt + a + b + c
q(t) = 6at + 6a + 2b
r(t) = 6a

So, to get the points x(0), x(1), x(2), ..., x(n-1), x(n) using only integer addition:

x(0) = X0
p(0) = a + b + c
q(0) = 6a + 2b
r(0) = 6a

x(t+1) = x(t) + p(t)
p(t+1) = p(t) + q(t)
q(t+1) = q(t) + r(t)
r(t+1) = r(t)

To use only integer addition, all numbers above are coded using a fixed point integer:

POS(t) = x(t) / 65536

The accuracy is relatively good for 64 steps. The cumulative errors in the calculation require greater accuracy for more steps than this.

For N dimensions, you need 4 values per dimension per spline path. Note that for continuous spline cubics, the x(0) point does not need to be computed for any except the first because it will be the x(n) of the previous cubic.