

 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 }

 


constant N = 32,
         N2 = (N * (N - 1) / 2),
         STEP = 0.05
 
constant inf = 1e300*1e300,
         nan = -(inf/inf)
 
sequence {xval, t_sin, t_cos, t_tan} = repeat(repeat(0,N),4)
 
for i=1 to N do
    xval[i] = (i-1) * STEP
    t_sin[i] = sin(xval[i])
    t_cos[i] = cos(xval[i])
    t_tan[i] = t_sin[i] / t_cos[i]
end for
 
enum R_SIN, R_COS, R_TAN, R_TRIG=$
 
sequence rhot = repeat(repeat(nan,N2),R_TRIG)
 
function rho(sequence x, y, integer rdx, int i, int n)
    if n<0 then return 0 end if
    if n=0 then return y[i+1] end if
 
    integer idx = (N - 1 - n) * (N - n) / 2 + i + 1;
    if rhot[rdx][idx]=nan then -- value not computed yet
        rhot[rdx][idx] = (x[i+1] - x[i+1 + n])
                        / (rho(x, y, rdx, i, n-1) - rho(x, y, rdx, i+1, n-1))
                        + rho(x, y, rdx, i+1, n-2)
    end if
    return rhot[rdx][idx]
end function
 
function thiele(sequence x, y, integer rdx, atom xin, integer n)
    if n>N-1 then return 1 end if
    return rho(x, y, rdx, 0, n) - rho(x, y, rdx, 0, n-2)
            + (xin-x[n+1]) / thiele(x, y, rdx, xin, n+1)
end function
 
constant fmt = iff(machine_bits()=32?"%32s : %.14f\n"
                                    :"%32s : %.17f\n")
printf(1,fmt,{"PI",PI})
printf(1,fmt,{"6*arcsin(0.5)",6*arcsin(0.5)})
printf(1,fmt,{"3*arccos(0.5)",3*arccos(0.5)})
printf(1,fmt,{"4*arctan(1)",4*arctan(1)})
 
printf(1,fmt,{"6*thiele(t_sin,xval,R_SIN,0.5,0)",6*thiele(t_sin,xval,R_SIN,0.5,0)})
printf(1,fmt,{"3*thiele(t_cos,xval,R_COS,0.5,0)",3*thiele(t_cos,xval,R_COS,0.5,0)})
printf(1,fmt,{"4*thiele(t_tan,xval,R_TAN,1,0)",4*thiele(t_tan,xval,R_TAN,1,0)})
