Code generation

Suppose you have worked on formula which you used Maple to derive, and you'd like to now include this in your C or C++ program. The worst you could do is simply try to copy it over, for as you know, copying by hand is one of the most excruciating and error prone operations. Instead, you can just ask Maple to do the work for you:

Find the interpolating cubic polynomial passing through the points $(t_n - 4h, f_{n - 4})$, $(t_n - 3h, f_{n - 3})$, $(t_n - 2h, f_{n - 2})$ and $(t_n, f_{n})$. We are assuming that the data-point at time $t_{n - 1} = t_n - h$ was lost, but we would still like to estimate the rate of change, the concavity, and the area underneath the curve of the data described by $f$.

Find the interpolating polynomial with respect to the variable $\tau$, and then differentiate with respect to $\tau$ either once or twice, and then evaluate the resulting polynomial at $\tau = t$, or integrate the polynomial with $\tau$ going from $t - h$ to $t$:

[> poly := CurveFitting['PolynomialInterpolation']( [t - 4*h, t - 3*h, t - 2*h, t],
                   [f[n - 4], f[n - 3], f[n - 2], f[n]], tau ):
[> simplify( eval( diff( poly, tau ), tau = t ), 'size' );

$\frac{13 f_{n}-9 f_{n -4}+32 f_{n -3}-36 f_{n -2}}{12 h}$

[> CodeGeneration[C]( %, 'optimize' );

t1 = (double) (1 / h * (13 * f[n - 1] - 9 * f[n - 5] + 32 * f[n - 4] - 36 * f[n - 3])) / 0.12e2;

[> simplify( eval( diff( poly, tau, tau ), tau = t ), 'size' );

$\frac{3 f_{n}-5 f_{n -4}+16 f_{n -3}-14 f_{n -2}}{4 h^{2}}$

[> CodeGeneration[C]( %, 'optimize' );

t2 = (int) (h * h);
t3 = (double) ((3 * f[n - 1] - 5 * f[n - 5] + 16 * f[n - 4] - 14 * f[n - 3]) / t2) / 0.4e1;

[> simplify( int( poly, tau = t - h..t ), 'size' );

$\frac{h \left(19 f_{n -4}+55 f_{n}-72 f_{n -3}+94 f_{n -2}\right)}{96}$

[> CodeGeneration[C]( %, 'optimize' );

t4 = (double) (h * (19 * f[n - 5] + 55 * f[n - 1] - 72 * f[n - 4] + 94 * f[n - 3])) / 0.96e2;

You'll notice that Maple may have tried to be too clever: Maples arrays are indexed starting at 1, and C or C++ arrays are indexed starting at 0, so Maple shifted the indices to the left by one. Whether or not this is desirable is up to you.

Note: I personally would continue to make additional changes to the code, using regular expressions, to make it more readable for humans (the compiler doesn't care if it 0.96e2 or 96.0, but honestly, a human trying to interpret the code would certainly prefer the latter), and we note that the optimizer assumed that $h$ was an integer, so I removed that optimization:

t1 =   (13.0*f[n] -  9.0*f[n - 4] + 32.0*f[n - 3] - 36.0*f[n - 2])/(12.0*h);
t3 =   ( 3.0*f[n] -  5.0*f[n - 4] + 16.0*f[n - 3] - 14.0*f[n - 2])/(4.0*h*h);
t4 = h*(55.0*f[n] + 19.0*f[n - 4] - 72.0*f[n - 3] + 94.0*f[n - 2])/96.0;