Thiele's interpolation formula

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})+{\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 }}}}}}

ρ represents the reciprocal difference, demonstrated here for reference:

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

\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})

\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 of the trig functions sin, cos and tan. e.g. for x from 0 by 0.05 to 1.55...
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 ½ = π
- 3 × cos^-1 ½ = π
- 4 × tan^-1 1 = π

Ada

thiele.ads: <lang Ada>with Ada.Numerics.Generic_Real_Arrays;

generic

  type Real is digits <>;

package Thiele is

  package Real_Arrays is new Ada.Numerics.Generic_Real_Arrays (Real);
  subtype Real_Array is Real_Arrays.Real_Vector;

  type Thiele_Interpolation (Length : Natural) is private;

  function Create (X, Y : Real_Array) return Thiele_Interpolation;
  function Inverse (T : Thiele_Interpolation; X : Real) return Real;

private

  type Thiele_Interpolation (Length : Natural) is record
     X, Y, RhoX : Real_Array (1 .. Length);
  end record;

end Thiele;</lang>

thiele.adb: <lang Ada>package body Thiele is

  use type Real_Array;

  function "/" (Left, Right : Real_Array) return Real_Array is
     Result : Real_Array (Left'Range);
  begin
     if Left'Length /= Right'Length then
        raise Constraint_Error with "arrays not same size";
     end if;
     for I in Result'Range loop
        Result (I) := Left (I) / Right (I);
     end loop;
     return Result;
  end "/";

  function Rho (X, Y : Real_Array) return Real_Array is
     N      : constant Natural                      := X'Length;
     P      : array (1 .. N) of Real_Array (1 .. N) :=
       (others => (others => 9.9));
     Result : Real_Array (1 .. N);
  begin
     P (1) (1 .. N)      := Y (1 .. N);
     P (2) (1 .. N - 1)  := (X (1 .. N - 1) - X (2 .. N)) /
       (P (1) (1 .. N - 1) - P (1) (2 .. N));
     for I in 3 .. N loop
        P (I) (1 .. N - I + 1)  := P (I - 2) (2 .. N - I + 2) +
          (X (1 .. N - I + 1) - X (I .. N)) /
          (P (I - 1) (1 .. N - I + 1) - P (I - 1) (2 .. N - I + 2));
     end loop;
     for I in X'Range loop
        Result (I) := P (I) (1);
     end loop;
     return Result;
  end Rho;

  function Create (X, Y : Real_Array) return Thiele_Interpolation is
  begin
     if X'Length < 3 then
        raise Constraint_Error with "at least 3 values";
     end if;
     if X'Length /= Y'Length then
        raise Constraint_Error with "input arrays not same size";
     end if;
     return (Length => X'Length, X => X, Y => Y, RhoX => Rho (X, Y));
  end Create;

  function Inverse (T : Thiele_Interpolation; X : Real) return Real is
     A : Real := 0.0;
  begin
     for I in reverse 3 .. T.Length loop
        A := (X - T.X (I - 1)) / (T.RhoX (I) - T.RhoX (I - 2) + A);
     end loop;
     return T.Y (1) + (X - T.X (1)) / (T.RhoX (2) + A);
  end Inverse;

end Thiele;</lang>

example: <lang Ada>with Ada.Text_IO; with Ada.Numerics.Generic_Elementary_Functions; with Thiele;

procedure Main is

  package Math is new Ada.Numerics.Generic_Elementary_Functions
    (Long_Float);
  package Float_Thiele is new Thiele (Long_Float);
  use Float_Thiele;

  Row_Count : Natural := 32;

  X_Values   : Real_Array (1 .. Row_Count);
  Sin_Values : Real_Array (1 .. Row_Count);
  Cos_Values : Real_Array (1 .. Row_Count);
  Tan_Values : Real_Array (1 .. Row_Count);

begin

  -- build table
  for I in 1 .. Row_Count loop
     X_Values (I)   := Long_Float (I) * 0.05 - 0.05;
     Sin_Values (I) := Math.Sin (X_Values (I));
     Cos_Values (I) := Math.Cos (X_Values (I));
     Tan_Values (I) := Math.Tan (X_Values (I));
  end loop;
  declare
     Sin : Thiele_Interpolation := Create (Sin_Values, X_Values);
     Cos : Thiele_Interpolation := Create (Cos_Values, X_Values);
     Tan : Thiele_Interpolation := Create (Tan_Values, X_Values);
  begin
     Ada.Text_IO.Put_Line
       ("Internal Math.Pi:    " &
        Long_Float'Image (Ada.Numerics.Pi));
     Ada.Text_IO.Put_Line
       ("Thiele 6*InvSin(0.5):" &
        Long_Float'Image (6.0 * Inverse (Sin, 0.5)));
     Ada.Text_IO.Put_Line
       ("Thiele 3*InvCos(0.5):" &
        Long_Float'Image (3.0 * Inverse (Cos, 0.5)));
     Ada.Text_IO.Put_Line
       ("Thiele 4*InvTan(1):  " &
        Long_Float'Image (4.0 * Inverse (Tan, 1.0)));
  end;

end Main;</lang>

output:

Internal Math.Pi:     3.14159265358979E+00
Thiele 6*InvSin(0.5): 3.14159265358979E+00
Thiele 3*InvCos(0.5): 3.14159265358979E+00
Thiele 4*InvTan(1):   3.14159265358979E+00

ALGOL 68

Works with: ALGOL 68 version Revision 1 - except the Currying (aka partial parametrisation) in test block is a proposal for ALGOL 68 Rev2

Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny - Currying is supported.

<lang algol68>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)
 ))

)</lang> Output:

pi estimate using sin interpolation: 3.1415926535898
pi estimate using cos interpolation: 3.1415926535898
pi estimate using tan interpolation: 3.1415926535898

C

The recursive relations of $\rho$ s can be made clearer: Given $N+1$ sampled points $x_{0},x_{1},\cdots x_{N}$ , rewrite the symbol $\rho$ as

\rho _{n,i}=\rho _{n}(x_{i},x_{i+1},\cdots x_{i+n})

where

1\leq n\leq N

,

with suplements

\rho _{0,i}=f(x_{i})\qquad {\text{and}}\qquad \rho _{n,i}=0\quad {\text{for}}\quad n<0

.

Now the recursive relation is simply

\rho _{n,i}=\displaystyle {x_{i}-x_{i+n} \over \rho _{n-1,i}-\rho _{n-1,i+1}}+\rho _{n-2,i+1},\ 0\leq i+n\leq N

Also note how $f(x_{1})$ in the interpolation formula can be replaced by $\rho _{0,1}$ ; define Thiele interpolation at step $n$ as

\displaystyle {F_{n}(x)=\rho _{n,1}-\rho _{n-2,1}+{x-x_{n+1} \over F_{n+1}(x)}}

with the termination $F_{N}(x)=1$ , and the interpolation formula is now $f(x)=F_{0}(x)$ , easily implemented as a recursive function.

Note that each $\rho _{n}$ needs to look up $\rho _{n-1}$ twice, so the total look ups go up as $O(2^{N})$ while there are only $O(N^{2})$ values. This is a text book situation for memoization. <lang c>#include <stdio.h>

include <string.h>
include <math.h>

define N 32
define N2 (N * (N - 1) / 2)
define STEP .05

double xval[N], t_sin[N], t_cos[N], t_tan[N];

/* rho tables, layout: rho_{n-1}(x0) rho_{n-2}(x0), rho_{n-1}(x1), .... rho_0(x0), rho_0(x1), ... rho_0(x_{n-1})

  rho_i row starts at index (n - 1 - i) * (n - i) / 2  	*/

double r_sin[N2], r_cos[N2], r_tan[N2];

/* both rho and thiele functions recursively resolve values as decribed by

  formulas.  rho is cached, thiele is not. */

/* rho_n(x_i, x_{i+1}, ..., x_{i + n}) */ double rho(double *x, double *y, double *r, int i, int n) { if (n < 0) return 0; if (!n) return y[i];

int idx = (N - 1 - n) * (N - n) / 2 + i; if (r[idx] != r[idx]) /* only happens if value not computed yet */ r[idx] = (x[i] - x[i + n]) / (rho(x, y, r, i, n - 1) - rho(x, y, r, i + 1, n - 1)) + rho(x, y, r, i + 1, n - 2); return r[idx]; }

double thiele(double *x, double *y, double *r, double xin, int n) { if (n > N - 1) return 1; return rho(x, y, r, 0, n) - rho(x, y, r, 0, n - 2) + (xin - x[n]) / thiele(x, y, r, xin, n + 1); }

define i_sin(x) thiele(t_sin, xval, r_sin, x, 0)
define i_cos(x) thiele(t_cos, xval, r_cos, x, 0)
define i_tan(x) thiele(t_tan, xval, r_tan, x, 0)

int main() { int i; for (i = 0; i < N; i++) { xval[i] = i * STEP; t_sin[i] = sin(xval[i]); t_cos[i] = cos(xval[i]); t_tan[i] = t_sin[i] / t_cos[i]; } for (i = 0; i < N2; i++) /* init rho tables to NaN */ r_sin[i] = r_cos[i] = r_tan[i] = 0/0.;

printf("%16.14f\n", 6 * i_sin(.5)); printf("%16.14f\n", 3 * i_cos(.5)); printf("%16.14f\n", 4 * i_tan(1.)); return 0; }</lang>output<lang>3.14159265358979 3.14159265358979 3.14159265358979</lang>

Common Lisp

Using the notations from above the C code instead of task desc. <lang lisp>;; 256 is heavy overkill, but hey, we memoized (defparameter *thiele-length* 256) (defparameter *rho-cache* (make-hash-table :test #'equal))

(defmacro make-thele-func (f name xx0 xx1)

 (let ((xv (gensym)) (yv (gensym))

(x0 (gensym)) (x1 (gensym)))

   `(let* ((,xv (make-array (1+ *thiele-length*)))

(,yv (make-array (1+ *thiele-length*))) (,x0 ,xx0) (,x1 ,xx1))

      (loop for i to *thiele-length* with x do

(setf x (+ ,x0 (* (/ (- ,x1 ,x0) *thiele-length*) i)) (aref ,yv i) x (aref ,xv i) (funcall ,f x)))

      (defun ,name (x) (thiele x ,yv ,xv, 0)))))

(defun rho (yv xv n i)

 (let (hit (key (list yv xv n i)))
   (if (setf hit (gethash key *rho-cache*))
     hit
     (setf (gethash key *rho-cache*)

(cond ((zerop n) (aref yv i)) ((minusp n) 0) (t (+ (rho yv xv (- n 2) (1+ i)) (/ (- (aref xv i) (aref xv (+ i n))) (- (rho yv xv (1- n) i) (rho yv xv (1- n) (1+ i)))))))))))

(defun thiele (x yv xv n)

 (if (= n *thiele-length*)
   1
   (+ (- (rho yv xv n 1) (rho yv xv (- n 2) 1))
      (/ (- x (aref xv (1+ n)))

(thiele x yv xv (1+ n))))))

(make-thele-func #'sin inv-sin 0 (/ pi 2)) (make-thele-func #'cos inv-cos 0 (/ pi 2)) (make-thele-func #'tan inv-tan 0 (/ pi 2.1)) ; tan(pi/2) is INF

(format t "~f~%" (* 6 (inv-sin .5))) (format t "~f~%" (* 3 (inv-cos .5))) (format t "~f~%" (* 4 (inv-tan 1)))</lang>output (SBCL):<lang>3.141592653589793 3.1415926535885172 3.141592653589819</lang>

D

<lang d>import std.stdio, std.range, std.array, std.algorithm, std.math;

struct Domain {

   const real b, e, s;

   auto range() const pure /*nothrow*/ @safe /*@nogc*/ {
       return iota(b, e + s, s);
   }

}

real eval0(alias RY, alias X, alias Y)(in real x) pure nothrow @safe @nogc {

   real a = 0.0L;
   foreach_reverse (immutable i; 2 .. X.length - 3)
       a = (x - X[i]) / (RY[i] - RY[i-2] + a);
   return Y[1] + (x - X[1]) / (RY[1] + a);

}

immutable struct Thiele {

   immutable real[] Y, X, rhoY, rhoX;

   this(real[] y, real[] x) immutable pure nothrow /*@safe*/
   in {
       assert(x.length > 2, "at leat 3 values");
       assert(x.length == y.length, "input arrays not of same size");
   } body {
       this.Y = y.idup;
       this.X = x.idup;
       rhoY = rhoN(Y, X);
       rhoX = rhoN(X, Y);
   }

   this(in real function(real) pure nothrow @safe @nogc f,
        Domain d = Domain(0.0L, 1.55L, 0.05L))
   immutable pure /*nothrow @safe*/ {
       auto xrng = d.range.array;
       this(xrng.map!f.array, xrng);
   }

   auto rhoN(immutable real[] y, immutable real[] x)
   pure nothrow @safe {
       immutable int N = x.length;
       auto p = new real[][](N, N);
       p[0][] = y[];
       p[1][0 .. $ - 1] = (x[0 .. $-1] - x[1 .. $]) /
                          (p[0][0 .. $-1] - p[0][1 .. $]);
       foreach (immutable int j; 2 .. N - 1) {
           immutable M = N - j - 1;
           p[j][0..M] = p[j-2][1..M+1] + (x[0..M] - x[j..M+j]) /
                        (p[j-1][0 .. M] - p[j-1][1 .. M+1]);
       }
       return p.map!q{ a[1] }.array;
   }

   alias eval = eval0!(rhoY, X, Y);
   alias inverse = eval0!(rhoX, Y, X);

}

void main() {

   // Can't pass sin, cos and tan directly.
   immutable tsin = Thiele(x => x.sin);
   immutable tcos = Thiele(x => x.cos);
   immutable ttan = Thiele(x => x.tan);

   writefln(" %d interpolating points\n", tsin.X.length);
   writefln("std.math.sin(0.5): %20.18f", 0.5L.sin);
   writefln("  Thiele sin(0.5): %20.18f\n", tsin.eval(0.5L));

   writefln("*%20.19f library constant", PI);
   writefln(" %20.19f 6 * inv_sin(0.5)", tsin.inverse(0.5L) * 6.0L);
   writefln(" %20.19f 3 * inv_cos(0.5)", tcos.inverse(0.5L) * 3.0L);
   writefln(" %20.19f 4 * inv_tan(1.0)", ttan.inverse(1.0L) * 4.0L);

}</lang>

Output:

 32 interpolating points

std.math.sin(0.5): 0.479425538604203000
  Thiele sin(0.5): 0.479425538604203000

*3.1415926535897932385 library constant
 3.1415926535897932380 6 * inv_sin(0.5)
 3.1415926535897932382 3 * inv_cos(0.5)
 3.1415926535897932382 4 * inv_tan(1.0)

Go

Translation of: ALGOL 68

<lang go>package main

import (

   "fmt"
   "math"

)

func main() {

   // task 1: build 32 row trig table
   const nn = 32
   const step = .05
   xVal := make([]float64, nn)
   tSin := make([]float64, nn)
   tCos := make([]float64, nn)
   tTan := make([]float64, nn)
   for i := range xVal {
       xVal[i] = float64(i) * step
       tSin[i], tCos[i] = math.Sincos(xVal[i])
       tTan[i] = tSin[i] / tCos[i]
   }
   // task 2: define inverses
   iSin := thieleInterpolator(tSin, xVal)
   iCos := thieleInterpolator(tCos, xVal)
   iTan := thieleInterpolator(tTan, xVal)
   // task 3: demonstrate identities
   fmt.Printf("%16.14f\n", 6*iSin(.5))
   fmt.Printf("%16.14f\n", 3*iCos(.5))
   fmt.Printf("%16.14f\n", 4*iTan(1))

}

func thieleInterpolator(x, y []float64) func(float64) float64 {

   n := len(x)
   ρ := make([][]float64, n)
   for i := range ρ {
       ρ[i] = make([]float64, n-i)
       ρ[i][0] = y[i]
   }
   for i := 0; i < n-1; i++ {
       ρ[i][1] = (x[i] - x[i+1]) / (ρ[i][0] - ρ[i+1][0])
   }
   for i := 2; i < n; i++ {
       for j := 0; j < n-i; j++ {
           ρ[j][i] = (x[j]-x[j+i])/(ρ[j][i-1]-ρ[j+1][i-1]) + ρ[j+1][i-2]
       }
   }
   // ρ0 used in closure.  the rest of ρ becomes garbage.
   ρ0 := ρ[0]
   return func(xin float64) float64 {
       var a float64
       for i := n - 1; i > 1; i-- {
           a = (xin - x[i-1]) / (ρ0[i] - ρ0[i-2] + a)
       }
       return y[0] + (xin-x[0])/(ρ0[1]+a)
   }

}</lang> Output:

3.14159265358979
3.14159265358979
3.14159265358980

Haskell

Caching of rho is automatic due to lazy lists. <lang haskell>thiele xs ys = f rho1 (tail xs) where f _ [] _ = 1 f r@(r0:r1:r2:rs) (x:xs) v = r2 - r0 + (v-x) / f (tail r) xs v

rho1 = (map ((!!1).(++[0])) rho)

rho = [0,0..] : [0,0..] : ys : rnext (tail rho) xs (tail xs) where rnext _ _ [] = [] rnext r@(r0:r1:rs) x xn = let z_ = zipWith in (z_ (+) (tail r0) (z_ (/) (z_ (-) x xn) (z_ (-) r1 (tail r1)))) : rnext (tail r) x (tail xn)

-- inverted interpolation function of f inv_interp f xs = thiele (map f xs) xs

main = do print $ 3.21 * inv_sin (sin (pi / 3.21)) print $ pi/1.2345 * inv_cos (cos (1.2345)) print $ 7 * inv_tan (tan (pi / 7)) where inv_sin = inv_interp sin $ div_pi 2 31 inv_cos = inv_interp cos $ div_pi 2 100 inv_tan = inv_interp tan $ div_pi 4 1000 -- because we can -- uniformly take n points from 0 to Pi/d div_pi d n = map (* (pi / (d * n))) [0..n]</lang>

Output:

3.141592653589795
3.1415926535897802
3.1415926535835275

J

<lang j> span =: {. - {: NB. head - tail spans =: span\ NB. apply span to successive infixes </lang>

   span 12 888 6 4 8 3
9
   4 spans 12 888 6 4 8 3
8 880 3

<lang j> NB. abscissae_of_knots coef ordinates_of_knots NB. returns the interpolation coefficients for eval coef =: 4 : 0

p =. _2 _{.,:y
for_i. i. # x do.
  p =. (p , ([: }. - }. p {~ _2:) + (x spans~ 2+]) % 2 spans - }. [: {: p"_) i
end.
x; , _ 1 {. p

)

NB. unknown_abscissae eval coefficients eval =: 4 : 0

'xx p' =. y
a =. 0
i =. <: # xx
while. 0 < i=.<:i do.
  a =. (x-i{xx)%-/(p{~i+2),(i{p),a
end.
(p{~>:i)+(x-i{xx)%(p{~i+2)+a

) </lang>

   trig_table =: 1 2 3 o./ angles =: 5r100*i.32

   0 _1 }. ": (;:'angle sin cos tan'),.<"1] 8j4": _ 5{.angles,trig_table
┌─────┬────────────────────────────────────────
│angle│  0.0000  0.0500  0.1000  0.1500  0.2000
├─────┼────────────────────────────────────────
│sin  │  0.0000  0.0500  0.0998  0.1494  0.1987
├─────┼────────────────────────────────────────
│cos  │  1.0000  0.9988  0.9950  0.9888  0.9801
├─────┼────────────────────────────────────────
│tan  │  0.0000  0.0500  0.1003  0.1511  0.2027
└─────┴────────────────────────────────────────


   ('Thiele pi';'error'),;/"1(,. 1p1&-)6 3 4 * 1r2 1r2 1 eval"0 1 trig_table coef"1 angles
┌─────────┬────────────┐
│Thiele pi│error       │
├─────────┼────────────┤
│3.14159  │_4.44089e_15│
├─────────┼────────────┤
│3.14159  │_4.44089e_16│
├─────────┼────────────┤
│3.14159  │_7.10543e_15│
└─────────┴────────────┘

<lang j> thiele =: 2 : 0

p =. _2 _{.,:n
for_i. i.#m do.
  p =. (p , ([: }. - }. p {~ _2:) + (m spans~ 2+]) % 2 spans - }. [: {: p"_) i
end.
p =. , _ 1 {. p
a =. 0
i =. <:#m
while. 0 < i=.<:i do.
  a =. (y-i{m)%-/(p{~i+2),(i{p),a
end.
(p{~>:i)+(y-i{m)%a+p{~i+2

) </lang>

   's c t' =: trig_table
   asin =: s thiele angles

   6*asin 0.5
3.14159

   1r5 * i.6
0 1r5 2r5 3r5 4r5 1
   100*(_1&o. %~ _1&o. - asin) 1r5*i.6   NB. % error arcsin
0 1.4052 4.50319 9.32495 16.9438 39.321

OCaml

This example shows how the accuracy changes with the degree of interpolation. The table 'columns' are only constructed implicitly during the recursive calculation of rdiff and thiele, but (as mentioned in the C code example) using memoization or explicit tabulation would speed up the calculation. The interpolation uses the nearest points around x for accuracy.

<lang OCaml>let xv, fv = fst, snd

let rec rdiff a l r =

  if l > r then 0.0 else
  if l = r then fv a.(l) else
  if l+1 = r then (xv a.(l) -. xv a.(r)) /. (fv a.(l) -. fv a.(r)) else
  (xv a.(l) -. xv a.(r)) /. (rdiff a l (r-1) -. rdiff a (l+1) r) +. rdiff a (l+1) (r-1)

let rec thiele x a a0 k n =

  if k = n then 1.0 else
  rdiff a a0 (a0+k) -. rdiff a a0 (a0+k-2) +. (x -. xv a.(a0+k)) /. thiele x a a0 (k+1) n

let interpolate x a n =

  let m = Array.length a in
  let dist i = abs_float (x -. xv a.(i)) in
  let nearer i j = if dist j < dist i then j else i in
  let rec closest i j = if j = m then i else closest (nearer i j) (j+1) in
  let c = closest 0 1 in
  let c' = if c < n/2 then 0 else if c > m-n then m-n else c-(n/2) in
  thiele x a c' 0 n

let table a b n f =

  let g i =
     let x = a +. (b-.a)*.(float i)/.(float (n-1)) in
     (f x, x) in
  Array.init n g

let [sin_tab; cos_tab; tan_tab] = List.map (table 0.0 1.55 32) [sin; cos; tan]

let test n =

  Printf.printf "\nDegree %d interpolation:\n" n;
  Printf.printf "6*arcsin(0.5) = %.15f\n" (6.0*.(interpolate 0.5 sin_tab n));
  Printf.printf "3*arccos(0.5) = %.15f\n" (3.0*.(interpolate 0.5 cos_tab n));
  Printf.printf "4*arctan(1.0) = %.15f\n" (4.0*.(interpolate 1.0 tan_tab n));;

List.iter test [8; 12; 16]</lang> Output:

Degree 8 interpolation:
6*arcsin(0.5) = 3.141592654456238
3*arccos(0.5) = 3.141592653520809
4*arctan(1.0) = 3.141592653437432

Degree 12 interpolation:
6*arcsin(0.5) = 3.141592653587590
3*arccos(0.5) = 3.141592653562618
4*arctan(1.0) = 3.141592653589756

Degree 16 interpolation:
6*arcsin(0.5) = 3.141592653589793
3*arccos(0.5) = 3.141592653589793
4*arctan(1.0) = 3.141592653589793

Perl 6

Works with: Rakudo version 2010.09-32

Implemented to parallel the (generalized) formula. (i.e. clearer, but naive and very slow.) <lang perl6>use v6;

reciprocal difference:

multi sub rho($f, @x where { +@x < 1 }) { 0 } # Identity multi sub rho($f, @x where { +@x == 1 }) { $f(@x[0]) } multi sub rho($f, @x where { +@x > 1 }) {

 my $ord = +@x;
 
 return
   ( @x[0] - @x[* -1] )            # ( x - x[n] )
   / ( rho($f, @x[^($ord -1)])     # / ( rho[n-1](x[0], ..., x[n-1])
       - rho($f, @x[1..^($ord)]) ) # - rho[n-1](x[1], ..., x[n]) )
   + rho($f, @x[1..^($ord -1)]);   # + rho[n-2](x[1], ..., x[n-1])

}

Thiele:

multi sub thiele($x, %f, $ord where { $ord == +%f }) { 1 } # Identity multi sub thiele($x, %f, $ord) {

 my $f = {%f{$^a}};                # f(x) as a table lookup
 
 # Caveat: depends on the fact that Rakudo maintains key order within hashes
 my $a = rho($f, %f.keys[^($ord +1)]);
 my $b = rho($f, %f.keys[^($ord -1)]);
 
 my $num = $x - %f.keys[$ord];
 my $cont = thiele($x, %f, $ord +1);
 
 # Thiele always takes this form:
 return $a - $b + ( $num / $cont );

}

1. Demo

sub mk-inv($fn, $d, $lim) {

 my %h;
 for 0..$lim { %h{ $fn($_ * $d) } = $_ * $d }
 return %h;

}

sub MAIN($tblsz) {

 my %invsin = mk-inv(&sin, 0.05, $tblsz);
 my %invcos = mk-inv(&cos, 0.05, $tblsz);
 my %invtan = mk-inv(&tan, 0.05, $tblsz);
 
 my $sin_pi = 6 * thiele(0.5, %invsin, 0);
 my $cos_pi = 3 * thiele(0.5, %invcos, 0);
 my $tan_pi = 4 * thiele(1.0, %invtan, 0);
 
 say "pi = {pi}";
 say "estimations using a table of $tblsz elements:";
 say "sin interpolation: $sin_pi";
 say "cos interpolation: $cos_pi";
 say "tan interpolation: $tan_pi";

}</lang>

Output (table size of 6 for want of resources):

pi = 3.14159265358979
estimations using a table of 6 elements:
sin interpolation: 3.14153363985515
cos interpolation: 1.68779321655997
tan interpolation: 3.14826236377727

PicoLisp

Translation of: C

<lang PicoLisp>(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) )</lang>

Test: <lang PicoLisp>(prinl (round (* 6 (iSin 0.5)) 15)) (prinl (round (* 3 (iCos 0.5)) 15)) (prinl (round (* 4 (iTan 1.0)) 15))</lang> Output:

3.141592653589793
3.141592653589793
3.141592653589793

PowerShell

<lang PowerShell>Function Reciprocal-Difference( [Double[][]] $function ) {

   $rho=@()
   $rho+=0
   $funcl = $function.length
   if( $funcl -gt 0 )
   {
       -2..($funcl-1) | ForEach-Object {
           $i=$_
           #Write-Host "$($i+1) - $($rho[$i+1]) - $($rho[$i+1].GetType())"
           $rho[$i+2] = $( 0..($funcl-$i-1) | Where-Object {$_ -lt $funcl} | ForEach-Object {
               $j=$_
               switch ($i) {
                   {$_ -lt 0 } { 0 }
                   {$_ -eq 0 } { $function[$j][1] }
                   {$_ -gt 0 } { ( $function[$j][0] - $function[$j+$i][0] ) / ( $rho[$i+1][$j] - $rho[$i+1][$j+1] ) + $rho[$i][$j+1] }
               }
           if( $_ -lt $funcl )
           {
               $rho += 0
           }
       })
       }
   }
   $rho

}

Function Thiele-Interpolation ( [Double[][]] $function ) {

   $funcl = $function.length
   $invoke = "{`n`tparam([Double] `$x)`n"
   if($funcl -gt 1)
   {
       $rho = Reciprocal-Difference $function
       ($funcl-1)..0 | ForEach-Object { 
           $invoke += "`t"
           $invoke += '$x{0} = {1} - {2}' -f $_, @($rho[$_+2])[0], @($rho[$_])[0]
           if($_ -lt ($funcl-1))
           {
               $invoke += ' + ( $x - {0} ) / $x{1} ' -f $function[$_][0], ($_+1)
           }
           $invoke += "`n"
       }
       $invoke+="`t`$x0`n}"
   } else {
       $invoke += "`t`$x`n}"
   }
   invoke-expression $invoke

}

$sint=@{}; 0..31 | ForEach-Object { $_ * 0.05 } | ForEach-Object { $sint[$_] = [Math]::sin($_) } $cost=@{}; 0..31 | ForEach-Object { $_ * 0.05 } | ForEach-Object { $cost[$_] = [Math]::cos($_) } $tant=@{}; 0..31 | ForEach-Object { $_ * 0.05 } | ForEach-Object { $tant[$_] = [Math]::tan($_) } $asint=New-Object 'Double[][]' 32,2; $sint.GetEnumerator() | Sort-Object Value | ForEach-Object {$i=0}{ $asint[$i][0] = $_.Value; $asint[$i][1] = $_.Name; $i++ } $acost=New-Object 'Double[][]' 32,2; $cost.GetEnumerator() | Sort-Object Value | ForEach-Object { $i=0 }{ $acost[$i][0] = $_.Value; $acost[$i][1] = $_.Name; $i++ } $atant=New-Object 'Double[][]' 32,2; $tant.GetEnumerator() | Sort-Object Value | ForEach-Object {$i=0}{ $atant[$i][0] = $_.Value; $atant[$i][1] = $_.Name; $i++ }

$asin = (Thiele-Interpolation $asint)

uncomment to see the function
"{$asin}"

6*$asin.InvokeReturnAsIs(.5) $acos = (Thiele-Interpolation $acost)

uncomment to see the function
"{$acos}"

3*$acos.InvokeReturnAsIs(.5) $atan = (Thiele-Interpolation $atant)

uncomment to see the function
"{$atan}"

4*$atan.InvokeReturnAsIs(1)</lang>

Python

Translation of: Go

<lang python>#!/usr/bin/env python3

import math

def thieleInterpolator(x, y):

   ρ = [[yi]*(len(y)-i) for i, yi in enumerate(y)]
   for i in range(len(ρ)-1):
       ρ[i][1] = (x[i] - x[i+1]) / (ρ[i][0] - ρ[i+1][0])
   for i in range(2, len(ρ)):
       for j in range(len(ρ)-i):
           ρ[j][i] = (x[j]-x[j+i]) / (ρ[j][i-1]-ρ[j+1][i-1]) + ρ[j+1][i-2]
   ρ0 = ρ[0]
   def t(xin):
       a = 0
       for i in range(len(ρ0)-1, 1, -1):
           a = (xin - x[i-1]) / (ρ0[i] - ρ0[i-2] + a)
       return y[0] + (xin-x[0]) / (ρ0[1]+a)
   return t

task 1: build 32 row trig table

xVal = [i*.05 for i in range(32)] tSin = [math.sin(x) for x in xVal] tCos = [math.cos(x) for x in xVal] tTan = [math.tan(x) for x in xVal]

task 2: define inverses

iSin = thieleInterpolator(tSin, xVal) iCos = thieleInterpolator(tCos, xVal) iTan = thieleInterpolator(tTan, xVal)

task 3: demonstrate identities

print('{:16.14f}'.format(6*iSin(.5))) print('{:16.14f}'.format(3*iCos(.5))) print('{:16.14f}'.format(4*iTan(1)))</lang>

Output:

3.14159265358979
3.14159265358979
3.14159265358980

Racket

lang racket

(define xs (for/vector ([x (in-range 0.0 1.6 0.05)]) x)) (define (x i) (vector-ref xs i))

(define-syntax define-table

 (syntax-rules ()
   [(_ f tf rf if) 
    (begin (define tab (for/vector ([x xs]) (f x)))
           (define (tf n) (vector-ref tab n))
           (define cache (make-vector (/ (* 32 31) 2) #f))
           (define (rf n thunk)
             (or (vector-ref cache n)
                 (let ([v (thunk)])
                   (vector-set! cache n v)
                   v)))
           (define (if t) (thiele tf x rf t 0)))]))

(define-table sin tsin rsin isin) (define-table cos tcos rcos icos) (define-table tan ttan rtan itan)

(define (rho x y r i n)

 (cond
   [(< n 0) 0]
   [(= n 0) (y i)]
   [else (r (+ (/ (* (- 32 1 n) (- 32 n)) 2) i)
            (λ() (+ (/ (- (x i) (x (+ i n)))
                       (- (rho x y r i (- n 1)) (rho x y r (+ i 1) (- n 1))))
                    (rho x y r (+ i 1) (- n 2)))))]))

(define (thiele x y r xin n)

 (cond
   [(> n 31) 1]
   [(+ (rho x y r 0 n) (- (rho x y r 0 (- n 2))) 
       (/ (- xin (x n)) (thiele x y r xin (+ n 1))))]))

(* 6 (isin 0.5)) (* 3 (icos 0.5)) (* 4 (itan 1.)) </lang> Output: <lang racket> 3.141592653589793 3.1415926535897936 3.1415926535897953 </lang>

Tcl

Works with: Tcl version 8.5

Translation of: D

<lang tcl>#

1. 1. Create a thiele-interpretation function with the given name that interpolates
  2. off the given table.

proc thiele {name : X -> F} {

   # Sanity check
   if {[llength $X] != [llength $F]} {

error "unequal length lists supplied: [llength $X] != [llength $F]"

   #
   ### Compute the table of reciprocal differences
   #
   set p [lrepeat [llength $X] [lrepeat [llength $X] 0.0]]
   set i 0
   foreach x0 [lrange $X 0 end-1] x1 [lrange $X 1 end] \

f0 [lrange $F 0 end-1] f1 [lrange $F 1 end] { lset p $i 0 $f0 lset p $i 1 [expr {($x0 - $x1) / ($f0 - $f1)}] lset p [incr i] 0 $f1

   }
   for {set j 2} {$j<[llength $X]-1} {incr j} {

for {set i 0} {$i<[llength $X]-$j} {incr i} { lset p $i $j [expr { [lindex $p $i+1 $j-2] + ([lindex $X $i] - [lindex $X $i+$j]) / ([lindex $p $i $j-1] - [lindex $p $i+1 $j-1]) }] }

   #
   ### Make pseudo-curried function that actually evaluates Thiele's formula
   #
   interp alias {} $name {} apply {{X rho f1 x} {

set a 0.0 foreach Xi [lreverse [lrange $X 2 end]] \ Ri [lreverse [lrange $rho 2 end]] \ Ri2 [lreverse [lrange $rho 0 end-2]] { set a [expr {($x - $Xi) / ($Ri - $Ri2 + $a)}] } expr {$f1 + ($x - [lindex $X 1]) / ([lindex $rho 1] + $a)}

   }} $X [lindex $p 1] [lindex $F 1]

}</lang> Demonstration code: <lang tcl>proc initThieleTest {} {

   for {set i 0} {$i < 32} {incr i} {

lappend trigTable(x) [set x [expr {0.05 * $i}]] lappend trigTable(sin) [expr {sin($x)}] lappend trigTable(cos) [expr {cos($x)}] lappend trigTable(tan) [expr {tan($x)}]

   thiele invSin : $trigTable(sin) -> $trigTable(x)
   thiele invCos : $trigTable(cos) -> $trigTable(x)
   thiele invTan : $trigTable(tan) -> $trigTable(x)

} initThieleTest puts "pi estimate using sin interpolation: [expr {6 * [invSin 0.5]}]" puts "pi estimate using cos interpolation: [expr {3 * [invCos 0.5]}]" puts "pi estimate using tan interpolation: [expr {4 * [invTan 1.0]}]"</lang> Output:

pi estimate using sin interpolation: 3.1415926535897936
pi estimate using cos interpolation: 3.141592653589793
pi estimate using tan interpolation: 3.141592653589794

zkl

Translation of: C

Please see the C example for the comments I've removed (this is an as pure-as-I-make-it translation). <lang zkl>const N=32, N2=(N * (N - 1) / 2), STEP=0.05;

fcn rho(xs,ys,rs, i,n){

  if (n < 0) return(0.0);
  if (not n) return(ys[i]);

  idx := (N - 1 - n) * (N - n) / 2 + i;
  if (Void==rs[idx])
     rs[idx] = (xs[i] - xs[i + n])

/ (rho(xs, ys, rs, i, n - 1) - rho(xs, ys, rs, i + 1, n - 1)) + rho(xs, ys, rs, i + 1, n - 2);

  return(rs[idx]);

}

fcn thiele(xs,ys,rs, xin, n){

  if (n > N - 1) return(1.0);
  rho(xs, ys, rs, 0, n) - rho(xs, ys, rs, 0, n - 2)
     + (xin - xs[n]) / thiele(xs, ys, rs, xin, n + 1);

}

///////////

reg t_sin=L(), t_cos=L(), t_tan=L(),

   r_sin=L(), r_cos=L(), r_tan=L(),  xval=L();

i_sin := thiele.fpM("11101",t_sin, xval, r_sin, 0); i_cos := thiele.fpM("11101",t_cos, xval, r_cos, 0); i_tan := thiele.fpM("11101",t_tan, xval, r_tan, 0);

foreach i in (N){

  xval.append(x:=STEP*i);
  t_sin.append(x.sin());
  t_cos.append(x.cos());
  t_tan.append(t_sin[i] / t_cos[i]);

} foreach i in (N2){ r_sin+Void; r_cos+Void; r_tan+Void; }

print("%16.14f\n".fmt( 6.0 * i_sin(0.5))); print("%16.14f\n".fmt( 3.0 * i_cos(0.5))); print("%16.14f\n".fmt( 4.0 * i_tan(1.0)));</lang>

Output:

3.14159265358979
3.14159265358979
3.14159265358979