Thiele's interpolation formula: Difference between revisions

Updated D entry
(Updated D entry)
Line 355:
 
struct Domain {
const real b, e, s;
 
auto range() const /*pure nothrow*/ {
// return iota(b, e + s, s);
return iota(cast()b, e + s, s);
}
}
Line 365 ⟶ 364:
real eval0(alias RY, alias X, alias Y)(in real x) pure nothrow {
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*/
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; // Not nothrow.
this.X = x.idup;
rhoY = rhoN(Y, X);
Line 385 ⟶ 384:
 
this(in real delegate(real) f,
Domain d=Domain(0.0L, 1.55L, 0.05L)) {
immutable /*pure nothrow*/ {
auto xrng = d.range().array();
this(array(map!f(xrng)),auto xrng) = d.range.array;
auto this(xrng = d.range()map!f.array(, xrng);
}
 
auto rhoN(immutable(real)[] y, immutable(real)[] x) /*pure*/ {
pure /*nothrow*/ {
int N = x.length;
real[][]immutable pint N = new real[][](N, N)x.length;
auto p = new real[0][](N, = y[]N); // array/vector operation follow
p[10][0..$-1] = (xy[0 .. $-1] - x[1 .. $]) /;
p[1][0 .. $ - 1] = (p[0]x[0 .. $-1] - p[0]x[1 .. $]); /
foreach (int j; 2 (p[0][0 .. N$-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().idup;
}
 
alias eval = eval0!(rhoY, X, Y) eval;
alias inverse = eval0!(rhoX, Y, X) inverse;
}
 
void main() {
// can't pass sin,cos,tan directly
autoimmutable tsin = immutable(Thiele)((real x) => x.sin(x));
autoimmutable tcos = immutable(Thiele)((real x) => x.cos(x));
autoimmutable ttan = immutable(Thiele)((real x) => x.tan(x));
 
writefln(" %d interpolating points\n", tsin.X.length);
writefln("std.math.sin(0.5): %20.18f", sin(0.5L).sin);
writefln(" Thiele sin(0.5): %20.18f\n", tsin.eval(0.5L));