• specified slope
• specified curvature
• "not-a-knot" smooth ends

a= (xa[khi]-x)/(h=xa[khi]-xa[klo]);
b= (x-xa[klo] )/h;
*y= a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]
   + b*b*b-b)*y2a[khi])*h*h/6;

b= x-xa[klo];
a= xa[khi]-x;
h= xa[khi]-xa[klo];
*y= ((ya[klo]-y2a[klo]/6*b*(a+h))*a+(ya[khi]
   -y2a[khi]/6*a*(b+h))*b)/h;

The spline fit function in Numerical Recipes In C has a problem with a potential aliasing of array elements which can force the compiler to generate less efficient code. The problem can be alleviated simply by changing the order of lines in the code. The prototype reads

void spline(float x[]
   float y[],int n,float y2[])

Intermediate results can be saved where they will be used again on the next step. Since the algorithm is recursive anyway, the additional recurrences will not pose an obstacle to any known compilers, as they might in parallelizable code. The improvement in speed is not as large as might be expected; sequential processors do not have enough registers and must use memory storage of these temporaries, and pipeline processors can handle many of the duplicate floating point operations by increasing parallelism.

Integration of the spline fit to find the area under a curve gives slightly over half the error of Simpson's rule integration using the same data points. This improvement is hardly worth the extra effort unless it is impossible to get the data spaced appropriately for Simpson's rule.

A cubic spline curve exaggerates errors or uncertainty in the data points. Usual ways of dealing with this are to use the "spline in tension" method, or to put "springs" between the curve and the knots and allow the curve to adjust to greater smoothness by deviating slightly from the input data. Tension or stiffness of the springs is an arbitrary parameter, requiring human interaction to get a visually satisfactory result. Repeated smoothing keeps on changing the data; unlike the Chebyshef smoothing, the adjusted splines never produce a curve which doesn't want to be adjusted further.

In practice, there is a fairly low limit (around 15) to the number of points which can be spline fit reliably without provision for tension or data adjustment. Demanding geometric problems, like fitting the contour of an airplane wing, can be handled by fitting splines in several segments.

A tactic which may be useful for geometric curve fitting by splines or any other method is to parameterize x, y, and z as a function of stair step distance. Ideally, the arc length along the curve would be a parameter; the stairstep distance is much easier to obtain. Fitting two or three functions of one independent variable is not much more time-consuming.

The theory is covered well in the textbooks, and a reliable set of C functions are presented in Numerical Recipes In C. I will concentrate on ways to use Chebyshef fits, and some implementation issues which should be addressed for best results.

Determination of Chebyshef coefficients requires the function values to be supplied at points dictated by the method, where the polynomial will match the specified values. Using specific points may seem awkward when the control points consist of arbitrary data at arbitrary locations, but Chebyshef fitting has been applied successfully by using linear interpolated values. Placement of the input points is less critical than for splines, although, since the function evaluations are concentrated toward the ends of the interval, the data are weighted more heavily there. The number of function evaluations needed for Chebyshef fitting is the minimum required to specify a given number of polynomial terms, so this method is as economical as any for functions which are expensive to evaluate. A dozen Chebyshef terms are likely to be sufficient for close approximation to any smooth curve.

The authors of Numerical Recipes in C advise against converting the Chebyshef polynomial into an ordinary polynomial before evaluation. This advice makes sense when there will be few evaluations of the polynomial, or when little is known about its behavior. Existing software often violates this advice, possibly because Clenshaw's algorithm was not well known until recently. Clenshaw's algorithm is 30 percent faster than evaluation by more obvious methods on sequential processors; even with careful coding for best performance, the difference drops to 20 percent on typical pipeline processors, but the small advantage in accuracy remains.

On co-processor architectures which have a built-in cosine function, you should avoid going through subroutine calls that narrow arguments and results to double precision. The greatest accuracy improvement comes from use of a long double value for PI. If the system does not support true long double constants, you may be able to generate a long double value which can be obtained by tricks such as

#define PI
(sqrt(3.125*3.125)+.01659265358979323846)

Pipelined processors should have a pipelined version of cos(), since it should be possible to calculate a group of four cosines as fast as two individual cosines. An intermediate level of performance may be obtained by supplying source code for cos() and using an inlining compiler option. While cosf() can be handled by a macro, the greater number of terms in cos() makes such an approach impractical. The subject of group math library functions could fill another article.

In order to get a series of sufficient accuracy to use in a library function, the whole Chebyshef fitting process has to be carried out in a higher precision. long double precision should be sufficient to calculate coefficients for double precision function approximations.

If a Chebyshef series is fitted directly to a function which has zeros, there will be errors near the zeros which are large relative to the local value of the function, even though they are less than the errors at the end of the interval. One way of getting around this is to subtract off the linear term of a series expansion around the zero, and make the Chebyshef fit to the remainder. Similar tactics might be applied where the function is known to be approximated by a more complicated function.

A better method for the standard library functions is to divide out the zero, so that the Chebyshef polynomial is fit to a function such as sin(x)/x. The Chebyshef polynomial will have only even terms, with the size of the odd terms calculated being an indication of numerical error. The end result will be an approximation which has less relative error over most of the interval than a minimax polynomial obtained by laborious iteration to minimize the maximum relative errors. The general strategy is to approximate a function by a known dominant term times a Chebyshef polynomial which takes care of the fine adjustments.

Most Chebyshef economized polynomials can be determined by straightforward programming using a system with long double arithmetic or a multiple precision package. In some cases, as with log((1+x)/(1-x)), it is better to derive a Taylor series which matches the form of the planned approximation and fit the Chebyshef polynomial to this series. Working out the cancelling terms algebraically rather than numerically avoids much of the roundoff error.

Normally, Chebyshef economized polynomials would be used only as a software development tool, to provide a best fit, most economical polynomial for quick evaluation of a function. If the number of evaluations of a curve is not great, or the convergence properties of a polynomial are unknown, it is better to use the Clenshaw algorithms for evaluation. Still, an economized polynomial might well be better than a least squares fit regression, and one might want to find such a polynomial for comparison.

exp(x) = 1/exp(-x)

tan(PI/2-x) = 1/tan(x)

tan(x) = -tan(-x)

(p(x)+q(x))/(p(x)-q(x))

q(2x)/p(2x)

You could arrive at similar results by deriving continued fraction approximations and converting to rational polynomials, a process which might be educational but wouldn't help in terms of practical results. The coefficients obtained analytically for rational polynomials (for unbounded intervals) are likely to be far from optimum for a closed interval.

People who play with these approximations have their favorite iterative numerical methods for adjusting a reasonable guess at the coefficients to greater accuracy. These calculations take almost twice as many bits of precision as the desired final accuracy; you would be lucky to get an approximation accurate to 12 decimal places on an 8087. Even IBM 3081 quad precision is marginal for obtaining a rational approximation for tan(). The author has had satisfactory results from a method acquired from Ruzinsky, which runs 10 to 20 seconds on a VAX 8550. With a multiple precision software package, higher accuracy fits should be feasible but time-consuming.

References

Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., Numerical Recipes In C, Cambridge University Press, 1988.

Ruzinsky, Steven A., "A Simple Minimax Algorithm," Dr. Dobb's Journal, #93, pp. 84-91, July 1984.