Thiele's interpolation formula: Difference between revisions

Content added Content deleted

Inline

Revision as of 01:56, 1 December 2010

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 (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 }}}}}}

ρ 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 = π

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

Translation of: ALGOL 68

- Currying, array slicing and exception handling removed, and uses some GCC extensions.

<lang c>#include <stdlib.h>

include <stdio.h>
include <math.h>

/* The MODE of lx and ly here should really be a UNION of "something double", "something COMPLex", and "something SYMBOLIC" ... */

double thiele(int upb_lx, double lx[], double ly[], double x) {

 int upb_xx=upb_lx+1;
 double *xx=lx-1, *yy=ly-1; /* shift base index to 1 */
 int n=upb_xx;

/* Assuming that the values of xx are distinct ... */

 double p [/*0:*/n-1 +1][/*1:*/n +1];
 int i, j;
 for(i=1; i<=upb_xx; i++) p[0][i]=yy[i];
 for(i=1; i<=n-1; i++)p[1][i]=(xx[i]-xx[1+i])/(p[0][i]-p[0][1+i]);
 for(i=2; i<=n-1; i++){
   for(j=1; j<=n-i; j++){
     p[i][j]=(xx[j]-xx[j+i])/(p[i-1][j]-p[i-1][j+1])+p[i-2][j+1];
   }
 }
 double a=0;
 for(i=n-1; i>=2; i--){ a=(x-xx[i])/(p[i][1]-p[i-2][1]+a); }
 return yy[1]+(x-xx[1])/(p[1][1]+a);

}

main(){

 double lwb_x=0, upb_x=1.55, delta_x = 0.05; 

 int upb_trig_table = ((upb_x-lwb_x)/delta_x);
 typedef double trig_table_t [/*0:*/upb_trig_table +1];
 trig_table_t x_OF_trig_table, sin_x_OF_trig_table, cos_x_OF_trig_table, tan_x_OF_trig_table;

 void init_trig_table(){
   int i;
   for(i = 0; i<= upb_trig_table; i++){
     double x = i*delta_x; 
     x_OF_trig_table[i]=x;
     sin_x_OF_trig_table[i]=sin(x);
     cos_x_OF_trig_table[i]=cos(x);
     tan_x_OF_trig_table[i]=tan(x);
   }
 }

 init_trig_table();

 int upb_row = 4;
 double inv_sin(double x){ return thiele(upb_trig_table, sin_x_OF_trig_table, x_OF_trig_table, x);}
 double inv_cos(double x){ return thiele(upb_trig_table, cos_x_OF_trig_table, x_OF_trig_table, x);}
 double inv_tan(double x){ return thiele(upb_trig_table, tan_x_OF_trig_table, x_OF_trig_table, x);}

 char *result_fmt = "pi estimate using %s interpolation: %0.13f\n";
 printf(result_fmt, "sin", 6*inv_sin(1.0/2));
 printf(result_fmt, "cos", 3*inv_cos(1.0/2));
 printf(result_fmt, "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

D

<lang d>import std.stdio ; import std.math ;

U[] myMap(U,V) (U delegate(V) f, V[] a) {

   V[] r = new V[](a.length) ;
   foreach(i, v ; a)
       r[i] = f(v) ;
   return r ;

}

T[] myRng(T)(T start, T end, T step) {

   T[] r ;
   for(; start < end ; start += step)
       r ~= start ;
   return r ;

}

alias real delegate(real) RealFun ; const real End = 1.55L , Step = 0.05L , Start = 0.0L ;

class Thiele {

   private real[] F, X, RhoY, RhoX ;
   const int NN ;
   
   this(real[] f, real[] x) {
       F   = f.dup ;
       X   = x.dup ;
       NN  = X.length ;
       assert(NN > 2, "at leat 3 values") ;
       assert(X.length == F.length, "input arrays not of same size") ;
       RhoY = rhoN(F, X) ;
       RhoX = rhoN(X, F) ;
   }

   this(RealFun fun, real start = Start, real end = End, real step = Step) {
       this(myMap(fun, myRng(start, end, step)), 
            myRng(start, end, step)) ;
   }

   private static real[] rhoN(real[] f, real[] x) {
       int N = x.length;
       real[][] p = new real[][] (N, N) ;
       for(int i = 0; i < N ; i++)
           p[i][0] = f[i] ;
       for(int i = 0; i < N - 1 ; i++)
           p[i][1] = (x[i] - x[i+1]) / (p[i][0] - p[i+1][0]) ;
       for(int j = 2; j < N - 1; j++)
           for(int i = 0; i < N - j - 1 ; i++)
               p[i][j] = p[i+1][j-2] + (x[i] - x[i+j]) /  
                         (p[i][j-1] - p[i+1][j-1]) ;                           ;
       return p[1].dup ;    
   }

   real evalY(real x) {
       real a = 0.0L ;
       for(int i = NN - 3; i >= 2 ; i--)
           a = (x - X[i]) / (RhoY[i] - RhoY[i-2] + a) ;
       return F[1] + (x - X[1]) / (RhoY[1] + a) ;        
   }
   real evalX(real y) { // inverse
       real a = 0.0L ;
       for(int i = NN - 3; i >= 2 ; i--)
           a = (y - F[i]) / (RhoX[i] - RhoX[i-2] + a) ;
       return X[1] + (y - F[1]) / (RhoX[1] + a) ;        
   }

}

void main() {

   auto fun = (real x) { return std.math.sin(x) ; } ;
   auto t = new Thiele(fun) ;
   writefln(" %d interpolating points", t.NN) ;
   writefln("std.math.sin(0.5) : %20.18f", std.math.sin(0.5L)) ;
   writefln("  Thiele sin(0.5) : %20.18f", t.evalY(0.5L)) ;    
   writefln("*%20.18f library constant", std.math.PI) ;
   writefln(" %20.18f 6*sin`(0.5)", t.evalX(0.5L) * 6.0L) ;
   t = new Thiele((real x) {return std.math.cos(x) ; }) ;
   writefln(" %20.18f 3*cos`(0.5)", t.evalX(0.5L) * 3.0L) ;
   t = new Thiele((real x) {return std.math.tan(x) ; }) ;
   writefln(" %20.18f 4*tan`(0.5)", t.evalX(1.0L) * 4.0L) ;

}</lang>

output:

 32 interpolating points
std.math.sin(0.5) : 0.479425538604203000
  Thiele sin(0.5) : 0.479425538604203000
*3.141592653589793239 library constant
 3.141592653589793238 6*sin`(0.5)
 3.141592653589793238 3*cos`(0.5)
 3.141592653589793238 4*tan`(0.5)

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

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>

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

@@ Line 372: / Line 372: @@
 #uncomment to see the function
 #"{$asin}"
-*$asin.invoke(.5)[0]
+*$asin.InvokeReturnAsIs(.5)
 $acos = (Thiele-Interpolation $acost)
 #uncomment to see the function
 #"{$acos}"
-*$acos.invoke(.5)[0]
+*$acos.InvokeReturnAsIs(.5)
 $atan = (Thiele-Interpolation $atant)
 #uncomment to see the function
 #"{$atan}"
-*$atan.invoke(1)[0]</lang>
+*$atan.InvokeReturnAsIs(1)</lang>
 =={{header|Tcl}}==