user
|  |
| --- |
| 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 {\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 {\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 {\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 {\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 {\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:

1. Building a   **32**   row *trig table* of values   for   

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

     π
     {\displaystyle {\displaystyle \pi }}
   * 3 × cos-1 ½ = 

     π
     {\displaystyle {\displaystyle \pi }}
   * 4 × tan-1 1 = 

     π
     {\displaystyle {\displaystyle \pi }}
assistant
(formerly Perl 6)

**Works with**: rakudo version 2018.09

Implemented to parallel the generalized formula, making for clearer, but slower, code. Offsetting that, the use of `Promise` allows concurrent calculations, so running all three types of interpolation should not take any longer than running just one (presuming available cores).

```
# reciprocal difference:
multi sub ρ(&f, @x where * < 1) { 0 } # Identity
multi sub ρ(&f, @x where * == 1) { &f(@x[0]) }
multi sub ρ(&f, @x where * > 1) {
    ( @x[0] - @x[* - 1] )       # ( x - x[n] )
    / (ρ(&f, @x[^(@x - 1)])     # / ( ρ[n-1](x[0], ..., x[n-1])
    - ρ(&f, @x[1..^@x]) )       # - ρ[n-1](x[1], ..., x[n]) )
    + ρ(&f, @x[1..^(@x - 1)]);  # + ρ[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
 
  # must sort hash keys to maintain order between invocations
  my $a = ρ(&f, %f.keys.sort[^($ord +1)]);
  my $b = ρ(&f, %f.keys.sort[^($ord -1)]);
 
  my $num = $x - %f.keys.sort[$ord];
  my $cont = thiele($x, %f, $ord +1);
 
  # Thiele always takes this form:
  return $a - $b + ( $num / $cont );
}
 
## Demo
sub mk-inv(&fn, $d, $lim) {
  my %h;
  for 0..$lim { %h{ &fn($_ * $d) } = $_ * $d }
  return %h;
}
 
sub MAIN($tblsz = 12) {

  my ($sin_pi, $cos_pi, $tan_pi);
  my $p1 = Promise.start( { my %invsin = mk-inv(&sin, 0.05, $tblsz); $sin_pi = 6 * thiele(0.5, %invsin, 0) } );
  my $p2 = Promise.start( { my %invcos = mk-inv(&cos, 0.05, $tblsz); $cos_pi = 3 * thiele(0.5, %invcos, 0) } );
  my $p3 = Promise.start( { my %invtan = mk-inv(&tan, 0.05, $tblsz); $tan_pi = 4 * thiele(1.0, %invtan, 0) } );
  await $p1, $p2, $p3;
 
  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";
}

```

Output:

```
pi = 3.14159265358979
estimations using a table of 12 elements:
sin interpolation: 3.14159265358961
cos interpolation: 3.1387286696692
tan interpolation: 3.14159090545243
```