Thiele's interpolation formula: Difference between revisions

Content added Content deleted
m (Link to continued fraction.)
(Updated D entry)
Line 352: Line 352:


=={{header|D}}==
=={{header|D}}==
{{works with|D|2.051}}
<lang d>import std.stdio, std.range, std.array, std.algorithm, std.math;
<lang d>import std.stdio, std.range, std.array, std.algorithm, std.math;


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

auto range() {
auto range() const /*pure nothrow*/ {
return iota(b, e + s, s);
// 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) {
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 (i; 2 .. X.length - 3)
Line 372: Line 373:
immutable real[] Y, X, rhoY, rhoX;
immutable real[] Y, X, rhoY, rhoX;


this(real[] y, real[] x)
this(real[] y, real[] x) /*pure*/ nothrow
in {
in {
assert(x.length > 2, "at leat 3 values");
assert(x.length > 2, "at leat 3 values");
Line 383: Line 384:
}
}


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


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


Line 407: Line 409:


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


writefln(" %d interpolating points\n", tsin.X.length);
writefln(" %d interpolating points\n", tsin.X.length);
Line 421: Line 423:
writefln(" %20.19f 4 * inv_tan(1.0)", ttan.inverse(1.0L) * 4.0L);
writefln(" %20.19f 4 * inv_tan(1.0)", ttan.inverse(1.0L) * 4.0L);
}</lang>
}</lang>
{{out}}

output:
<pre> 32 interpolating points
<pre> 32 interpolating points


Line 432: Line 433:
3.1415926535897932382 3 * inv_cos(0.5)
3.1415926535897932382 3 * inv_cos(0.5)
3.1415926535897932382 4 * inv_tan(1.0)</pre>
3.1415926535897932382 4 * inv_tan(1.0)</pre>

=={{header|Go}}==
=={{header|Go}}==
{{trans|ALGOL 68}}
{{trans|ALGOL 68}}
Retrieved from "https://rosettacode.org/wiki/Thiele%27s_interpolation_formula"