Thiele's interpolation formula: Difference between revisions

Content added Content deleted
(Updated D entry)
Line 355: Line 355:


struct Domain {
struct Domain {
real b, e, s;
const real b, e, s;


auto range() const /*pure nothrow*/ {
auto range() const /*pure nothrow*/ {
// return iota(b, e + s, s);
return iota(b, e + s, s);
return iota(cast()b, e + s, s);
}
}
}
}
Line 365: Line 364:
real eval0(alias RY, alias X, alias Y)(in real x) pure nothrow {
real eval0(alias RY, alias X, alias Y)(in real x) pure nothrow {
real a = 0.0L;
real a = 0.0L;
foreach_reverse (i; 2 .. X.length - 3)
foreach_reverse (immutable i; 2 .. X.length - 3)
a = (x - X[i]) / (RY[i] - RY[i-2] + a);
a = (x - X[i]) / (RY[i] - RY[i-2] + a);
return Y[1] + (x - X[1]) / (RY[1] + a);
return Y[1] + (x - X[1]) / (RY[1] + a);
}
}


struct Thiele {
immutable struct Thiele {
immutable real[] Y, X, rhoY, rhoX;
immutable real[] Y, X, rhoY, rhoX;


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


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


auto rhoN(immutable(real)[] y, immutable(real)[] x) /*pure*/ {
auto rhoN(immutable(real)[] y, immutable(real)[] x)
pure /*nothrow*/ {
int N = x.length;
real[][] p = new real[][](N, N);
immutable int N = x.length;
p[0][] = y[]; // array/vector operation follow
auto p = new real[][](N, N);
p[1][0..$-1] = (x[0 .. $-1] - x[1 .. $]) /
p[0][] = y[];
(p[0][0 .. $-1] - p[0][1 .. $]);
p[1][0 .. $ - 1] = (x[0 .. $-1] - x[1 .. $]) /
foreach (int j; 2 .. N - 1) {
(p[0][0 .. $-1] - p[0][1 .. $]);
foreach (immutable int j; 2 .. N - 1) {
immutable M = N - j - 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][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]);
(p[j-1][0 .. M] - p[j-1][1 .. M+1]);
}
}
return p.map!q{a[1]}().array().idup;
return p.map!q{ a[1] }.array;
}
}


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


void main() {
void main() {
// can't pass sin,cos,tan directly
// can't pass sin,cos,tan directly
auto tsin = Thiele((real x) => sin(x));
immutable tsin = immutable(Thiele)(x => x.sin);
auto tcos = Thiele((real x) => cos(x));
immutable tcos = immutable(Thiele)(x => x.cos);
auto ttan = Thiele((real x) => tan(x));
immutable ttan = immutable(Thiele)(x => x.tan);


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


Retrieved from "https://rosettacode.org/wiki/Thiele%27s_interpolation_formula"