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.)

References