

 This page uses content from Wikipedia. The original article was at Thiele's interpolation formula. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)


Thiele's interpolation formula is an interpolation formula for a function f(•) of a single variable.   It is expressed as a continued fraction:

 



f
(
x
)
=
f
(

x

1


)
+








x
−

x

1












ρ

1


(

x

1


,

x

2


)
+








x
−

x

2












ρ

2


(

x

1


,

x

2


,

x

3


)
−
f
(

x

1


)
+








x
−

x

3












ρ

3


(

x

1


,

x

2


,

x

3


,

x

4


)
−

ρ

1


(

x

1


,

x

2


)
+
⋯

















{\displaystyle f(x)=f(x_{1})+{\cfrac {x-x_{1}}{\rho _{1}(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 _{1}(x_{1},x_{2})+\cdots }}}}}}}






ρ


{\displaystyle \rho }

   represents the   reciprocal difference,   demonstrated here for reference:

 




ρ

1


(

x

0


,

x

1


)
=




x

0


−

x

1




f
(

x

0


)
−
f
(

x

1


)





{\displaystyle \rho _{1}(x_{0},x_{1})={\frac {x_{0}-x_{1}}{f(x_{0})-f(x_{1})}}}


 




ρ

2


(

x

0


,

x

1


,

x

2


)
=




x

0


−

x

2





ρ

1


(

x

0


,

x

1


)
−

ρ

1


(

x

1


,

x

2


)



+
f
(

x

1


)


{\displaystyle \rho _{2}(x_{0},x_{1},x_{2})={\frac {x_{0}-x_{2}}{\rho _{1}(x_{0},x_{1})-\rho _{1}(x_{1},x_{2})}}+f(x_{1})}


 




ρ

n


(

x

0


,

x

1


,
…
,

x

n


)
=




x

0


−

x

n





ρ

n
−
1


(

x

0


,

x

1


,
…
,

x

n
−
1


)
−

ρ

n
−
1


(

x

1


,

x

2


,
…
,

x

n


)



+

ρ

n
−
2


(

x

1


,
…
,

x

n
−
1


)


{\displaystyle \rho _{n}(x_{0},x_{1},\ldots ,x_{n})={\frac {x_{0}-x_{n}}{\rho _{n-1}(x_{0},x_{1},\ldots ,x_{n-1})-\rho _{n-1}(x_{1},x_{2},\ldots ,x_{n})}}+\rho _{n-2}(x_{1},\ldots ,x_{n-1})}


Demonstrate Thiele's interpolation function by:

 Building a   32   row trig table of values   for   



x


{\displaystyle x}

   from   0   by   0.05   to   1.55   of the trig functions:
   sin
   cos
   tan
 Using columns from this table define an inverse - using Thiele's interpolation - for each trig function;
 Finally: demonstrate the following well known trigonometric identities:
    6 × sin-1 ½ = 



π


{\displaystyle \pi }


    3 × cos-1 ½ = 



π


{\displaystyle \pi }


    4 × tan-1 1 = 



π


{\displaystyle \pi }

 


PROC raise exception = ([]STRING msg)VOID: ( putf(stand error,("Exception:", $" "g$, msg, $l$)); stop );
 
# 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
    raise exception(("Unequal length arrays supplied: ",whole(UPB xx,0)," NE ",whole(UPB yy,0))); SKIP
  FI
END;
 
test:(
  FORMAT real fmt = $g(0,real width-2)$;
 
  REAL lwb x=0, upb x=1.55, delta x = 0.05; 
 
  [0:ENTIER ((upb x-lwb x)/delta x)]STRUCT(REAL x, sin x, cos x, tan x) trig table;
 
  PROC init trig table = VOID: 
    FOR i FROM LWB trig table TO UPB trig table DO 
      REAL x = lwb x+i*delta x; 
      trig table[i]:=(x, sin(x), cos(x), tan(x))
    OD;
 
  init trig table;
 
# Curry the thiele function to create matching inverse trigonometric functions #
  PROC (REAL)REAL inv sin = thiele(sin x OF trig table, x OF trig table,),
                  inv cos = thiele(cos x OF trig table, x OF trig table,),
                  inv tan = thiele(tan x OF trig table, x OF trig table,);
 
  printf(($"pi estimate using "g" interpolation: "f(real fmt)l$, 
    "sin", 6*inv sin(1/2),
    "cos", 3*inv cos(1/2),
    "tan", 4*inv tan(1)
  ))
)