

 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 }

 


(scl 17)
(load "@lib/math.l")
 
(setq
   *X-Table (range 0.0 1.55 0.05)
   *SinTable (mapcar sin *X-Table)
   *CosTable (mapcar cos *X-Table)
   *TanTable (mapcar tan *X-Table)
   *TrigRows (length *X-Table) )
 
(let N2 (>> 1 (* *TrigRows (dec *TrigRows)))
   (setq
      *InvSinTable (need N2)
      *InvCosTable (need N2)
      *InvTanTable (need N2) ) )
 
(de rho (Tbl Inv I N)
   (cond
      ((lt0 N) 0)
      ((=0 N) (get *X-Table I))
      (T
         (let Idx (+ I (>> 1 (* (- *TrigRows 1 N) (- *TrigRows N))))
            (or
               (get Inv Idx)
               (set (nth Inv Idx)  # only happens if value not computed yet
                  (+
                     (rho Tbl Inv (inc I) (- N 2))
                     (*/
                        (- (get Tbl I) (get Tbl (+ I N)))
                        1.0
                        (-
                           (rho Tbl Inv I (dec N))
                           (rho Tbl Inv (inc I) (dec N)) ) ) ) ) ) ) ) ) )
 
(de thiele (Tbl Inv X N)
   (if (> N *TrigRows)
      1.0
      (+
         (-
            (rho Tbl Inv 1 (dec N))
            (rho Tbl Inv 1 (- N 3)) )
         (*/
            (- X (get Tbl N))
            1.0
            (thiele Tbl Inv X (inc N)) ) ) ) )
 
(de iSin (X)
   (thiele *SinTable *InvSinTable X 1) )
 
(de iCos (X)
   (thiele *CosTable *InvCosTable X 1) )
 
(de iTan (X)
   (thiele *TanTable *InvTanTable 1.0 1) )