Thiele's interpolation formula

In mathematics, Thiele's interpolation formula is a formula that defines a rational function $f(x)$ from a finite set of inputs $x_i$ and their function values $f(x_i)$. The problem of generating a function whose graph passes through a given set of function values is called interpolation. This interpolation formula is named after the Danish mathematician Thorvald N. Thiele. It is expressed as a continued fraction, where ρ represents the reciprocal difference:

$f(x) = f(x_1) + \cfrac{x-x_1}{\rho(x_1,x_2) + \cfrac{x-x_2}{\rho_2(x_1,x_2,x_3) - f(x_1) + \cfrac{x-x_3}{\rho_3(x_1,x_2,x_3,x_4) - \rho(x_1,x_2) + \cdots}}}$

Implementation

Here is a version of Thiele interpolation in Algol 68:

 ¢ The MODE of lx and ly here should really be a UNION of "something REAL", "something COMPLex",
and "something SYMBOLIC" ... ¢
PROC thiele:=([]REAL lx,ly, REAL x) REAL:
BEGIN
[]REAL xx=lx[@1],yy=ly[@1];
INT n=UPB xx;
IF UPB yy=n THEN
¢ Assuming that the values of xx are distinct ... ¢
[0:n-1,1:n]REAL p;
p[0,]:=yy[];
FOR i TO n-1 DO p[1,i]:=(xx[i]-xx[1+i])/(p[0,i]-p[0,1+i]) OD;
FOR i FROM 2 TO n-1 DO
FOR j TO n-i DO
p[i,j]:=(xx[j]-xx[j+i])/(p[i-1,j]-p[i-1,j+1])+p[i-2,j+1]
OD
OD;
REAL a:=0;
FOR i FROM n-1 BY -1 TO 2 DO a:=(x-xx[i])/(p[i,1]-p[i-2,1]+a) OD;
yy[1]+(x-xx[1])/(p[1,1]+a)
ELSE
error; ~ ¢ Unequal length arrays supplied ¢
FI
END;


Although it works in most cases, it is sensitive to input values, especially those due to equally spaced abscissæ. (Essentially, in such a case, it may reduce to the ratio of two polynomials, which may have factors in common - thus yielding a 0/0 form.)