Thiele's interpolation formula: Difference between revisions
(→{{header|C}}: replaced with readable code) |
m (→{{header|C}}: comments) |
||
Line 244: | Line 244: | ||
int idx = (N - 1 - n) * (N - n) / 2 + i; |
int idx = (N - 1 - n) * (N - n) / 2 + i; |
||
if (r[idx] != r[idx]) |
if (r[idx] != r[idx]) /* only happens if value not computed yet */ |
||
r[idx] = (x[i] - x[i + n]) |
r[idx] = (x[i] - x[i + n]) |
||
/ (rho(x, y, r, i, n - 1) - rho(x, y, r, i + 1, n - 1)) |
/ (rho(x, y, r, i, n - 1) - rho(x, y, r, i + 1, n - 1)) |
||
Line 272: | Line 272: | ||
} |
} |
||
for (i = 0; i < N2; i++) |
for (i = 0; i < N2; i++) |
||
/* init rho tables to NaN */ |
|||
r_sin[i] = r_cos[i] = r_tan[i] = 0/0.; |
r_sin[i] = r_cos[i] = r_tan[i] = 0/0.; |
||
Revision as of 07:43, 27 July 2011
You are encouraged to solve this task according to the task description, using any language you may know.
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:
ρ represents the reciprocal difference, demonstrated here for reference:
Demonstrate Thiele's interpolation function by:
- Building a 32 row trig table of values of the trig functions sin, cos and tan. e.g. for x from 0 by 0.05 to 1.55...
- 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 ½ = π
- 3 × cos-1 ½ = π
- 4 × tan-1 1 = π
Ada
thiele.ads: <lang Ada>with Ada.Numerics.Generic_Real_Arrays;
generic
type Real is digits <>;
package Thiele is
package Real_Arrays is new Ada.Numerics.Generic_Real_Arrays (Real); subtype Real_Array is Real_Arrays.Real_Vector;
type Thiele_Interpolation (Length : Natural) is private;
function Create (X, Y : Real_Array) return Thiele_Interpolation; function Inverse (T : Thiele_Interpolation; X : Real) return Real;
private
type Thiele_Interpolation (Length : Natural) is record X, Y, RhoX : Real_Array (1 .. Length); end record;
end Thiele;</lang>
thiele.adb: <lang Ada>package body Thiele is
use type Real_Array;
function "/" (Left, Right : Real_Array) return Real_Array is Result : Real_Array (Left'Range); begin if Left'Length /= Right'Length then raise Constraint_Error with "arrays not same size"; end if; for I in Result'Range loop Result (I) := Left (I) / Right (I); end loop; return Result; end "/";
function Rho (X, Y : Real_Array) return Real_Array is N : constant Natural := X'Length; P : array (1 .. N) of Real_Array (1 .. N) := (others => (others => 9.9)); Result : Real_Array (1 .. N); begin P (1) (1 .. N) := Y (1 .. N); P (2) (1 .. N - 1) := (X (1 .. N - 1) - X (2 .. N)) / (P (1) (1 .. N - 1) - P (1) (2 .. N)); for I in 3 .. N loop P (I) (1 .. N - I + 1) := P (I - 2) (2 .. N - I + 2) + (X (1 .. N - I + 1) - X (I .. N)) / (P (I - 1) (1 .. N - I + 1) - P (I - 1) (2 .. N - I + 2)); end loop; for I in X'Range loop Result (I) := P (I) (1); end loop; return Result; end Rho;
function Create (X, Y : Real_Array) return Thiele_Interpolation is begin if X'Length < 3 then raise Constraint_Error with "at least 3 values"; end if; if X'Length /= Y'Length then raise Constraint_Error with "input arrays not same size"; end if; return (Length => X'Length, X => X, Y => Y, RhoX => Rho (X, Y)); end Create;
function Inverse (T : Thiele_Interpolation; X : Real) return Real is A : Real := 0.0; begin for I in reverse 3 .. T.Length loop A := (X - T.X (I - 1)) / (T.RhoX (I) - T.RhoX (I - 2) + A); end loop; return T.Y (1) + (X - T.X (1)) / (T.RhoX (2) + A); end Inverse;
end Thiele;</lang>
example: <lang Ada>with Ada.Text_IO; with Ada.Numerics.Generic_Elementary_Functions; with Thiele;
procedure Main is
package Math is new Ada.Numerics.Generic_Elementary_Functions (Long_Float); package Float_Thiele is new Thiele (Long_Float); use Float_Thiele;
Row_Count : Natural := 32;
X_Values : Real_Array (1 .. Row_Count); Sin_Values : Real_Array (1 .. Row_Count); Cos_Values : Real_Array (1 .. Row_Count); Tan_Values : Real_Array (1 .. Row_Count);
begin
-- build table for I in 1 .. Row_Count loop X_Values (I) := Long_Float (I) * 0.05 - 0.05; Sin_Values (I) := Math.Sin (X_Values (I)); Cos_Values (I) := Math.Cos (X_Values (I)); Tan_Values (I) := Math.Tan (X_Values (I)); end loop; declare Sin : Thiele_Interpolation := Create (Sin_Values, X_Values); Cos : Thiele_Interpolation := Create (Cos_Values, X_Values); Tan : Thiele_Interpolation := Create (Tan_Values, X_Values); begin Ada.Text_IO.Put_Line ("Internal Math.Pi: " & Long_Float'Image (Ada.Numerics.Pi)); Ada.Text_IO.Put_Line ("Thiele 6*InvSin(0.5):" & Long_Float'Image (6.0 * Inverse (Sin, 0.5))); Ada.Text_IO.Put_Line ("Thiele 3*InvCos(0.5):" & Long_Float'Image (3.0 * Inverse (Cos, 0.5))); Ada.Text_IO.Put_Line ("Thiele 4*InvTan(1): " & Long_Float'Image (4.0 * Inverse (Tan, 1.0))); end;
end Main;</lang>
output:
Internal Math.Pi: 3.14159265358979E+00 Thiele 6*InvSin(0.5): 3.14159265358979E+00 Thiele 3*InvCos(0.5): 3.14159265358979E+00 Thiele 4*InvTan(1): 3.14159265358979E+00
ALGOL 68
<lang algol68>PROC raise exception = ([]STRING msg)VOID: ( putf(stand error,("Exception:", $" "g$, msg, $l$)); stop );
- The MODE of lx and ly here should really be a UNION of "something REAL",
"something COMPLex", and "something SYMBOLIC" ... #
PROC thiele=([]REAL lx,ly, REAL x) REAL: BEGIN
[]REAL xx=lx[@1],yy=ly[@1]; INT n=UPB xx; IF UPB yy=n THEN
- Assuming that the values of xx are distinct ... #
[0:n-1,1:n]REAL p; p[0,]:=yy[]; FOR i TO n-1 DO p[1,i]:=(xx[i]-xx[1+i])/(p[0,i]-p[0,1+i]) OD; FOR i FROM 2 TO n-1 DO FOR j TO n-i DO p[i,j]:=(xx[j]-xx[j+i])/(p[i-1,j]-p[i-1,j+1])+p[i-2,j+1] OD OD; REAL a:=0; FOR i FROM n-1 BY -1 TO 2 DO a:=(x-xx[i])/(p[i,1]-p[i-2,1]+a) OD; yy[1]+(x-xx[1])/(p[1,1]+a) ELSE raise exception(("Unequal length arrays supplied: ",whole(UPB xx,0)," NE ",whole(UPB yy,0))); SKIP FI
END;
test:(
FORMAT real fmt = $g(0,real width-2)$;
REAL lwb x=0, upb x=1.55, delta x = 0.05;
[0:ENTIER ((upb x-lwb x)/delta x)]STRUCT(REAL x, sin x, cos x, tan x) trig table;
PROC init trig table = VOID: FOR i FROM LWB trig table TO UPB trig table DO REAL x = lwb x+i*delta x; trig table[i]:=(x, sin(x), cos(x), tan(x)) OD;
init trig table;
- Curry the thiele function to create matching inverse trigonometric functions #
PROC (REAL)REAL inv sin = thiele(sin x OF trig table, x OF trig table,), inv cos = thiele(cos x OF trig table, x OF trig table,), inv tan = thiele(tan x OF trig table, x OF trig table,);
printf(($"pi estimate using "g" interpolation: "f(real fmt)l$, "sin", 6*inv sin(1/2), "cos", 3*inv cos(1/2), "tan", 4*inv tan(1) ))
)</lang> Output:
pi estimate using sin interpolation: 3.1415926535898 pi estimate using cos interpolation: 3.1415926535898 pi estimate using tan interpolation: 3.1415926535898
C
<lang c>#include <stdio.h>
- include <string.h>
- include <math.h>
- define N 32
- define N2 (N * (N + 1) / 2)
- define STEP .05
double xval[N], t_sin[N], t_cos[N], t_tan[N];
/* rho tables, layout: rho_{n-1}(x0) rho_{n-2}(x0), rho_{n-1}(x1), .... rho_0(x0), rho_0(x1), ... rho_0(x_{n-1})
rho_i row starts at index (n - 1 - i) * (n - i) / 2 */
double r_sin[N2], r_cos[N2], r_tan[N2];
/* both rho and thiele functions recursively resolve values as decribed by
formulas. rho is cached, thiele is not. */
/* rho_n(x_i, x_{i+1}, ..., x_{i + n}) */ double rho(double *x, double *y, double *r, int i, int n) { if (n < 0) return 0; if (!n) return y[i];
int idx = (N - 1 - n) * (N - n) / 2 + i; if (r[idx] != r[idx]) /* only happens if value not computed yet */ r[idx] = (x[i] - x[i + n]) / (rho(x, y, r, i, n - 1) - rho(x, y, r, i + 1, n - 1)) + rho(x, y, r, i + 1, n - 2); return r[idx]; }
double thiele(double *x, double *y, double *r, double xin, int n) { if (n > N - 1) return 1; return rho(x, y, r, 0, n) - rho(x, y, r, 0, n - 2) + (xin - x[n]) / thiele(x, y, r, xin, n + 1); }
- define i_sin(x) thiele(t_sin, xval, r_sin, x, 0)
- define i_cos(x) thiele(t_cos, xval, r_cos, x, 0)
- define i_tan(x) thiele(t_tan, xval, r_tan, x, 0)
int main() { int i; for (i = 0; i < N; i++) { xval[i] = i * STEP; t_sin[i] = sin(xval[i]); t_cos[i] = cos(xval[i]); t_tan[i] = t_sin[i] / t_cos[i]; } for (i = 0; i < N2; i++) /* init rho tables to NaN */ r_sin[i] = r_cos[i] = r_tan[i] = 0/0.;
printf("%16.14f\n", 6 * i_sin(.5)); printf("%16.14f\n", 3 * i_cos(.5)); printf("%16.14f\n", 4 * i_tan(1.)); return 0; }</lang>output<lang>3.14159265358979 3.14159265358979 3.14159265358979</lang>
D
<lang d>import std.stdio, std.range, std.array, std.algorithm, std.math;
struct Domain {
real b, e, s; auto range() { return iota(b, e + s, s); }
}
pure nothrow real eval0(alias RY, alias X, alias Y)(in real x) {
real a = 0.0L; foreach_reverse (i; 2 .. X.length - 3) a = (x - X[i]) / (RY[i] - RY[i-2] + a); return Y[1] + (x - X[1]) / (RY[1] + a);
}
struct Thiele {
immutable real[] Y, X, rhoY, rhoX;
this(real[] y, real[] x) 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(real delegate(real) f, Domain d=Domain(0.0L, 1.55L, 0.05L)) { auto xrng = array(d.range()); this(array(map!f(xrng)), xrng); }
auto rhoN(immutable(real)[] y, immutable(real)[] x) { 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(map!q{a[1]}(p)).idup; }
alias eval0!(rhoY, X, Y) eval; alias eval0!(rhoX, Y, X) inverse;
}
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); writefln("std.math.sin(0.5): %20.18f", sin(0.5L)); 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);
}</lang>
output:
32 interpolating points std.math.sin(0.5): 0.479425538604203000 Thiele sin(0.5): 0.479425538604203000 *3.1415926535897932385 library constant 3.1415926535897932380 6 * inv_sin(0.5) 3.1415926535897932382 3 * inv_cos(0.5) 3.1415926535897932382 4 * inv_tan(1.0)
J
<lang j> span =: {. - {: NB. head - tail spans =: span\ NB. apply span to successive infixes </lang>
span 12 888 6 4 8 3 9 4 spans 12 888 6 4 8 3 8 880 3
<lang j> NB. abscissae_of_knots coef ordinates_of_knots NB. returns the interpolation coefficients for eval coef =: 4 : 0
p =. _2 _{.,:y for_i. i. # x do. p =. (p , ([: }. - }. p {~ _2:) + (x spans~ 2+]) % 2 spans - }. [: {: p"_) i end. x; , _ 1 {. p
)
NB. unknown_abscissae eval coefficients eval =: 4 : 0
'xx p' =. y a =. 0 i =. <: # xx while. 0 < i=.<:i do. a =. (x-i{xx)%-/(p{~i+2),(i{p),a end. (p{~>:i)+(x-i{xx)%(p{~i+2)+a
) </lang>
trig_table =: 1 2 3 o./ angles =: 5r100*i.32 0 _1 }. ": (;:'angle sin cos tan'),.<"1] 8j4": _ 5{.angles,trig_table ┌─────┬──────────────────────────────────────── │angle│ 0.0000 0.0500 0.1000 0.1500 0.2000 ├─────┼──────────────────────────────────────── │sin │ 0.0000 0.0500 0.0998 0.1494 0.1987 ├─────┼──────────────────────────────────────── │cos │ 1.0000 0.9988 0.9950 0.9888 0.9801 ├─────┼──────────────────────────────────────── │tan │ 0.0000 0.0500 0.1003 0.1511 0.2027 └─────┴──────────────────────────────────────── ('Thiele pi';'error'),;/"1(,. 1p1&-)6 3 4 * 1r2 1r2 1 eval"0 1 trig_table coef"1 angles ┌─────────┬────────────┐ │Thiele pi│error │ ├─────────┼────────────┤ │3.14159 │_4.44089e_15│ ├─────────┼────────────┤ │3.14159 │_4.44089e_16│ ├─────────┼────────────┤ │3.14159 │_7.10543e_15│ └─────────┴────────────┘
<lang j> thiele =: 2 : 0
p =. _2 _{.,:n for_i. i.#m do. p =. (p , ([: }. - }. p {~ _2:) + (m spans~ 2+]) % 2 spans - }. [: {: p"_) i end. p =. , _ 1 {. p a =. 0 i =. <:#m while. 0 < i=.<:i do. a =. (y-i{m)%-/(p{~i+2),(i{p),a end. (p{~>:i)+(y-i{m)%a+p{~i+2
) </lang>
's c t' =: trig_table asin =: s thiele angles 6*asin 0.5 3.14159 1r5 * i.6 0 1r5 2r5 3r5 4r5 1 100*(_1&o. %~ _1&o. - asin) 1r5*i.6 NB. % error arcsin 0 1.4052 4.50319 9.32495 16.9438 39.321
Perl 6
Implemented to parallel the (generalized) formula. (i.e. clearer, but naive and very slow.) <lang perl6>use v6;
- reciprocal difference:
multi sub rho($f, @x where { +@x < 1 }) { 0 } # Identity multi sub rho($f, @x where { +@x == 1 }) { $f(@x[0]) } multi sub rho($f, @x where { +@x > 1 }) {
my $ord = +@x; return ( @x[0] - @x[* -1] ) # ( x - x[n] ) / ( rho($f, @x[^($ord -1)]) # / ( rho[n-1](x[0], ..., x[n-1]) - rho($f, @x[1..^($ord)]) ) # - rho[n-1](x[1], ..., x[n]) ) + rho($f, @x[1..^($ord -1)]); # + rho[n-2](x[1], ..., x[n-1])
}
- Thiele:
multi sub thiele($x, %f, $ord where { $ord == +%f }) { 1 } # Identity multi sub thiele($x, %f, $ord) {
my $f = {%f{$^a}}; # f(x) as a table lookup # Caveat: depends on the fact that Rakudo maintains key order within hashes my $a = rho($f, %f.keys[^($ord +1)]); my $b = rho($f, %f.keys[^($ord -1)]); my $num = $x - %f.keys[$ord]; my $cont = thiele($x, %f, $ord +1); # Thiele always takes this form: return $a - $b + ( $num / $cont );
}
- Demo
sub mk-inv($fn, $d, $lim) {
my %h; for 0..$lim { %h{ $fn($_ * $d) } = $_ * $d } return %h;
}
sub MAIN($tblsz) {
my %invsin = mk-inv(&sin, 0.05, $tblsz); my %invcos = mk-inv(&cos, 0.05, $tblsz); my %invtan = mk-inv(&tan, 0.05, $tblsz); my $sin_pi = 6 * thiele(0.5, %invsin, 0); my $cos_pi = 3 * thiele(0.5, %invcos, 0); my $tan_pi = 4 * thiele(1.0, %invtan, 0); say "pi = {pi}"; say "estimations using a table of $tblsz elements:"; say "sin interpolation: $sin_pi"; say "cos interpolation: $cos_pi"; say "tan interpolation: $tan_pi";
}</lang>
Output (table size of 6 for want of resources):
pi = 3.14159265358979 estimations using a table of 6 elements: sin interpolation: 3.14153363985515 cos interpolation: 1.68779321655997 tan interpolation: 3.14826236377727
PowerShell
<lang PowerShell>Function Reciprocal-Difference( [Double[][]] $function ) {
$rho=@() $rho+=0 $funcl = $function.length if( $funcl -gt 0 ) { -2..($funcl-1) | ForEach-Object { $i=$_ #Write-Host "$($i+1) - $($rho[$i+1]) - $($rho[$i+1].GetType())" $rho[$i+2] = $( 0..($funcl-$i-1) | Where-Object {$_ -lt $funcl} | ForEach-Object { $j=$_ switch ($i) { {$_ -lt 0 } { 0 } {$_ -eq 0 } { $function[$j][1] } {$_ -gt 0 } { ( $function[$j][0] - $function[$j+$i][0] ) / ( $rho[$i+1][$j] - $rho[$i+1][$j+1] ) + $rho[$i][$j+1] } } if( $_ -lt $funcl ) { $rho += 0 } }) } } $rho
}
Function Thiele-Interpolation ( [Double[][]] $function ) {
$funcl = $function.length $invoke = "{`n`tparam([Double] `$x)`n" if($funcl -gt 1) { $rho = Reciprocal-Difference $function ($funcl-1)..0 | ForEach-Object { $invoke += "`t" $invoke += '$x{0} = {1} - {2}' -f $_, @($rho[$_+2])[0], @($rho[$_])[0] if($_ -lt ($funcl-1)) { $invoke += ' + ( $x - {0} ) / $x{1} ' -f $function[$_][0], ($_+1) } $invoke += "`n" } $invoke+="`t`$x0`n}" } else { $invoke += "`t`$x`n}" } invoke-expression $invoke
}
$sint=@{}; 0..31 | ForEach-Object { $_ * 0.05 } | ForEach-Object { $sint[$_] = [Math]::sin($_) } $cost=@{}; 0..31 | ForEach-Object { $_ * 0.05 } | ForEach-Object { $cost[$_] = [Math]::cos($_) } $tant=@{}; 0..31 | ForEach-Object { $_ * 0.05 } | ForEach-Object { $tant[$_] = [Math]::tan($_) } $asint=New-Object 'Double[][]' 32,2; $sint.GetEnumerator() | Sort-Object Value | ForEach-Object {$i=0}{ $asint[$i][0] = $_.Value; $asint[$i][1] = $_.Name; $i++ } $acost=New-Object 'Double[][]' 32,2; $cost.GetEnumerator() | Sort-Object Value | ForEach-Object { $i=0 }{ $acost[$i][0] = $_.Value; $acost[$i][1] = $_.Name; $i++ } $atant=New-Object 'Double[][]' 32,2; $tant.GetEnumerator() | Sort-Object Value | ForEach-Object {$i=0}{ $atant[$i][0] = $_.Value; $atant[$i][1] = $_.Name; $i++ }
$asin = (Thiele-Interpolation $asint)
- uncomment to see the function
- "{$asin}"
6*$asin.InvokeReturnAsIs(.5) $acos = (Thiele-Interpolation $acost)
- uncomment to see the function
- "{$acos}"
3*$acos.InvokeReturnAsIs(.5) $atan = (Thiele-Interpolation $atant)
- uncomment to see the function
- "{$atan}"
4*$atan.InvokeReturnAsIs(1)</lang>
Tcl
<lang tcl>#
- Create a thiele-interpretation function with the given name that interpolates
- off the given table.
proc thiele {name : X -> F} {
# Sanity check if {[llength $X] != [llength $F]} {
error "unequal length lists supplied: [llength $X] != [llength $F]"
}
# ### Compute the table of reciprocal differences # set p [lrepeat [llength $X] [lrepeat [llength $X] 0.0]] set i 0 foreach x0 [lrange $X 0 end-1] x1 [lrange $X 1 end] \
f0 [lrange $F 0 end-1] f1 [lrange $F 1 end] { lset p $i 0 $f0 lset p $i 1 [expr {($x0 - $x1) / ($f0 - $f1)}] lset p [incr i] 0 $f1
} for {set j 2} {$j<[llength $X]-1} {incr j} {
for {set i 0} {$i<[llength $X]-$j} {incr i} { lset p $i $j [expr { [lindex $p $i+1 $j-2] + ([lindex $X $i] - [lindex $X $i+$j]) / ([lindex $p $i $j-1] - [lindex $p $i+1 $j-1]) }] }
}
# ### Make pseudo-curried function that actually evaluates Thiele's formula # interp alias {} $name {} apply {{X rho f1 x} {
set a 0.0 foreach Xi [lreverse [lrange $X 2 end]] \ Ri [lreverse [lrange $rho 2 end]] \ Ri2 [lreverse [lrange $rho 0 end-2]] { set a [expr {($x - $Xi) / ($Ri - $Ri2 + $a)}] } expr {$f1 + ($x - [lindex $X 1]) / ([lindex $rho 1] + $a)}
}} $X [lindex $p 1] [lindex $F 1]
}</lang> Demonstration code: <lang tcl>proc initThieleTest {} {
for {set i 0} {$i < 32} {incr i} {
lappend trigTable(x) [set x [expr {0.05 * $i}]] lappend trigTable(sin) [expr {sin($x)}] lappend trigTable(cos) [expr {cos($x)}] lappend trigTable(tan) [expr {tan($x)}]
}
thiele invSin : $trigTable(sin) -> $trigTable(x) thiele invCos : $trigTable(cos) -> $trigTable(x) thiele invTan : $trigTable(tan) -> $trigTable(x)
} initThieleTest puts "pi estimate using sin interpolation: [expr {6 * [invSin 0.5]}]" puts "pi estimate using cos interpolation: [expr {3 * [invCos 0.5]}]" puts "pi estimate using tan interpolation: [expr {4 * [invTan 1.0]}]"</lang> Output:
pi estimate using sin interpolation: 3.1415926535897936 pi estimate using cos interpolation: 3.141592653589793 pi estimate using tan interpolation: 3.141592653589794