

 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 }

 


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