Jump to content

Thiele's interpolation formula: Difference between revisions

Updated D entry
m (Link to continued fraction.)
(Updated D entry)
Line 352:
 
=={{header|D}}==
{{works with|D|2.051}}
<lang d>import std.stdio, std.range, std.array, std.algorithm, std.math;
 
struct Domain {
real b, e, s;
 
auto range() const /*pure nothrow*/ {
// return iota(b, e + s, s);
return iota(cast()b, e + s, s);
}
}
 
pure nothrow real eval0(alias RY, alias X, alias Y)(in real x) pure nothrow {
real a = 0.0L;
foreach_reverse (i; 2 .. X.length - 3)
Line 372 ⟶ 373:
immutable real[] Y, X, rhoY, rhoX;
 
this(real[] y, real[] x) /*pure*/ nothrow
in {
assert(x.length > 2, "at leat 3 values");
Line 383 ⟶ 384:
}
 
this(in real delegate(real) f, Domain d=Domain(0.0L, 1.55L, 0.05L)) {
auto xrngDomain d= arrayDomain(d0.range(0L, 1.55L, 0.05L)); {
auto xrng = d.range().array();
this(array(map!f(xrng)), xrng);
}
 
auto rhoN(immutable(real)[] y, immutable(real)[] x) /*pure*/ {
int N = x.length;
real[][] p = new real[][] (N, N);
p[0][] = y[]; // array/vector operation follow
p[1][0..$-1] = (x[0 .. $-1] - x[1 .. $]) /
(p[0][0 .. $-1] - p[0][1 .. $]);
foreach (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 array(p.map!q{a[1]}(p).array().idup;
}
 
Line 407 ⟶ 409:
 
void main() {
// can't pass sin,cos,tan directly :-(
auto tsin = Thiele((real x){ return=> sin(x); });
auto tcos = Thiele((real x){ return=> cos(x); });
auto ttan = Thiele((real x){ return=> tan(x); });
 
writefln(" %d interpolating points\n", tsin.X.length);
Line 421 ⟶ 423:
writefln(" %20.19f 4 * inv_tan(1.0)", ttan.inverse(1.0L) * 4.0L);
}</lang>
{{out}}
 
output:
<pre> 32 interpolating points
 
Line 432 ⟶ 433:
3.1415926535897932382 3 * inv_cos(0.5)
3.1415926535897932382 4 * inv_tan(1.0)</pre>
 
=={{header|Go}}==
{{trans|ALGOL 68}}
Cookies help us deliver our services. By using our services, you agree to our use of cookies.