Rosetta Code

Thiele's interpolation formula: Difference between revisions

Page Discussion

Thiele's interpolation formula (view source)

Revision as of 21:27, 9 March 2024

4,118 bytes added ,  2 months ago
Added FreeBASIC
(Added FreeBASIC)
 
(3 intermediate revisions by 3 users not shown)
Line 30:
{{trans|Python}}
 
<langsyntaxhighlight lang="11l">F thieleInterpolator(x, y)
V ρ = enumerate(y).map((i, yi) -> [yi] * (@y.len - i))
L(i) 0 .< ρ.len - 1
Line 54:
print(‘#.14’.format(6 * iSin(0.5)))
print(‘#.14’.format(3 * iCos(0.5)))
print(‘#.14’.format(4 * iTan(1)))</langsyntaxhighlight>
 
{{out}}
Line 65:
=={{header|Ada}}==
thiele.ads:
<langsyntaxhighlight Adalang="ada">with Ada.Numerics.Generic_Real_Arrays;
 
generic
Line 81:
X, Y, RhoX : Real_Array (1 .. Length);
end record;
end Thiele;</langsyntaxhighlight>
 
thiele.adb:
<langsyntaxhighlight Adalang="ada">package body Thiele is
use type Real_Array;
 
Line 139:
end Inverse;
 
end Thiele;</langsyntaxhighlight>
 
example:
<langsyntaxhighlight Adalang="ada">with Ada.Text_IO;
with Ada.Numerics.Generic_Elementary_Functions;
with Thiele;
Line 184:
Long_Float'Image (4.0 * Inverse (Tan, 1.0)));
end;
end Main;</langsyntaxhighlight>
 
output:
Line 197:
{{works with|ALGOL 68G|Any - tested with release [http://sourceforge.net/projects/algol68/files/algol68g/algol68g-1.18.0/algol68g-1.18.0-9h.tiny.el5.centos.fc11.i386.rpm/download 1.18.0-9h.tiny] - Currying is supported.}}
{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - due to extensive use of '''format'''[ted] ''transput'' - Also slicing a '''struct''' table and currying unimplemented.}}
<langsyntaxhighlight 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",
Line 249:
"tan", 4*inv tan(1)
))
)</langsyntaxhighlight>
Output:
<pre>
Line 277:
 
Note that each <math>\rho_n</math> needs to look up <math>\rho_{n-1}</math> twice, so the total look ups go up as <math>O(2^N)</math> while there are only <math>O(N^2)</math> values. This is a text book situation for memoization.
<langsyntaxhighlight lang="c">#include <stdio.h>
#include <string.h>
#include <math.h>
Line 340:
printf("%16.14f\n", 4 * i_tan(1.));
return 0;
}</langsyntaxhighlight>
{{out}}
<pre>3.14159265358979
Line 349:
{{trans|C}}
{{works with|C++14}}
<langsyntaxhighlight lang="cpp">#include <cmath>
#include <iostream>
#include <iomanip>
Line 397:
 
return 0;
}</langsyntaxhighlight>
{{out}}
<pre>3.141592653589793115997963
Line 405:
=={{header|Common Lisp}}==
Using the notations from above the C code instead of task desc.
<langsyntaxhighlight lang="lisp">;; 256 is heavy overkill, but hey, we memoized
(defparameter *thiele-length* 256)
(defparameter *rho-cache* (make-hash-table :test #'equal))
Line 448:
(format t "~f~%" (* 6 (inv-sin .5)))
(format t "~f~%" (* 3 (inv-cos .5)))
(format t "~f~%" (* 4 (inv-tan 1)))</langsyntaxhighlight>output (SBCL):<syntaxhighlight lang="text">3.141592653589793
3.1415926535885172
3.141592653589819</langsyntaxhighlight>
 
=={{header|D}}==
<langsyntaxhighlight lang="d">import std.stdio, std.range, std.array, std.algorithm, std.math;
 
struct Domain {
Line 524:
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);
}</langsyntaxhighlight>
{{out}}
<pre> 32 interpolating points
Line 535:
3.1415926535897932382 3 * inv_cos(0.5)
3.1415926535897932382 4 * inv_tan(1.0)</pre>
 
=={{header|FreeBASIC}}==
{{trans|Phix}}
<syntaxhighlight lang="vbnet">Const As Integer n1 = 32
Const As Integer n2 = (n1 * (n1 - 1) / 2)
Const As Double paso = 0.05
Const As Double INF = 1e308
Const As Double NaN = -(INF/INF)
 
Dim As Double xVal(n1), tSin(n1), tCos(n1), tTan(n1)
For i As Integer = 1 To n1
xVal(i) = (i-1) * paso
tSin(i) = Sin(xVal(i))
tCos(i) = Cos(xVal(i))
tTan(i) = Tan(xVal(i))
Next i
 
Dim As Integer rSin, rCos, rTan, rTrig
 
Dim Shared rhot(rTrig, n2) As Double
For i As Integer = 0 To rTrig
For j As Integer = 0 To n2
rhot(i, j) = NaN
Next j
Next i
 
Function rho(x() As Double, y() As Double, Byval rdx As Integer, _
Byval i As Integer, Byval n As Integer) As Double
If n < 0 Then Return 0
If n = 0 Then Return y(i+1)
Dim As Integer idx = (n1 - 1 - n) * (n1 - n) / 2 + i + 1
If rhot(rdx, idx) = NaN Then 'valor aún no calculado
rhot(rdx, idx) = (x(i+1) - x(i+1 + n)) _
/ (rho(x(), y(), rdx, i, n-1) - rho(x(), y(), rdx, i+1, n-1)) _
+ rho(x(), y(), rdx, i+1, n-2)
End If
Return rhot(rdx, idx)
End Function
 
Function thieleInterpolator(x() As Double, y() As Double, _
Byval rdx As Integer, Byval xin As Double, Byval n As Integer) As Double
If n > n1-1 Then Return 1
Return rho(x(), y(), rdx, 0, n) - rho(x(), y(), rdx, 0, n-2) _
+ (xin-x(n+1)) / thieleInterpolator(x(), y(), rdx, xin, n+1)
End Function
 
Print " PI : 3.141592653589793"
Print " 6*arcsin(0.5) : "; 6 * Asin(0.5)
Print " 3*arccos(0.5) : "; 3 * Acos(0.5)
Print " 4*arctan(1.0) : "; 4 * Atn(1.0)
 
Print "6*thiele(tSin,xVal,rSin,0.5,0) : "; 6 * thieleInterpolator(tSin(), xVal(), rSin, 0.5, 0)
Print "3*thiele(tCos,xVal,rCos,0.5,0) : "; 3 * thieleInterpolator(tCos(), xVal(), rCos, 0.5, 0)
Print "4*thiele(tTan,xVal,rTan,1.0,0) : "; 4 * thieleInterpolator(tTan(), xVal(), rTan, 1.0, 0)
 
Sleep</syntaxhighlight>
 
=={{header|Go}}==
{{trans|ALGOL 68}}
<langsyntaxhighlight lang="go">package main
 
import (
Line 592 ⟶ 649:
return y[0] + (xin-x[0])/(ρ0[1]+a)
}
}</langsyntaxhighlight>
Output:
<pre>
Line 602 ⟶ 659:
=={{header|Haskell}}==
Caching of rho is automatic due to lazy lists.
<langsyntaxhighlight lang="haskell">thiele :: [Double] -> [Double] -> Double -> Double
thiele xs ys = f rho1 (tail xs)
where
Line 633 ⟶ 690:
[(sin, (2, 31)), (cos, (2, 100)), (tan, (4, 1000))]
-- N points taken uniformly from 0 to Pi/d
div_pi (d, n) = (* (pi / (d * n))) <$> [0 .. n]</langsyntaxhighlight>
{{out}}
<pre>3.141592653589795
Line 640 ⟶ 697:
 
=={{header|J}}==
<syntaxhighlight lang="j">
<lang j>
span =: {. - {: NB. head - tail
spans =: span\ NB. apply span to successive infixes
</syntaxhighlight>
</lang>
 
<pre>
Line 652 ⟶ 709:
</pre>
 
<syntaxhighlight lang="j">
<lang j>
NB. abscissae_of_knots coef ordinates_of_knots
NB. returns the interpolation coefficients for eval
Line 673 ⟶ 730:
(p{~>:i)+(x-i{xx)%(p{~i+2)+a
)
</syntaxhighlight>
</lang>
 
<pre>
Line 702 ⟶ 759:
</pre>
 
<syntaxhighlight lang="j">
<lang j>
thiele =: 2 : 0
p =. _2 _{.,:n
Line 716 ⟶ 773:
(p{~>:i)+(y-i{m)%a+p{~i+2
)
</syntaxhighlight>
</lang>
 
<pre>
Line 733 ⟶ 790:
=={{header|Java}}==
{{trans|C}}
<langsyntaxhighlight lang="java">import static java.lang.Math.*;
 
public class Test {
Line 787 ⟶ 844:
System.out.printf("%16.14f%n", 4 * thiele(t_tan, xval, r_tan, 1.0, 0));
}
}</langsyntaxhighlight>
<pre>3.14159265358979
3.14159265358979
Line 795 ⟶ 852:
Accuracy improves with a larger table and smaller step size.
{{trans|C}}
<langsyntaxhighlight lang="julia">const N = 256
const N2 = N * div(N - 1, 2)
const step = 0.01
Line 840 ⟶ 897:
 
thiele_tables()
</langsyntaxhighlight>{{output}}<pre>
3.1415926535898335
3.141592653589818
Line 848 ⟶ 905:
=={{header|Kotlin}}==
{{trans|C}}
<langsyntaxhighlight lang="scala">// version 1.1.2
 
const val N = 32
Line 890 ⟶ 947:
println("%16.14f".format(3 * thiele(tcos, xval, rcos, 0.5, 0)))
println("%16.14f".format(4 * thiele(ttan, xval, rtan, 1.0, 0)))
}</langsyntaxhighlight>
 
{{out}}
Line 898 ⟶ 955:
3.14159265358980
</pre>
 
=={{header|Mathematica}} / {{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">num = 32;
num2 = num (num - 1)/2;
step = 0.05;
ClearAll[\[Rho], Thiele]
\[Rho][x_List, y_List, i_Integer, n_Integer] := Module[{idx},
If[n < 0,
0
,
If[n == 0,
y[[i + 1]]
,
idx = (num - 1 - n) (num - n)/2 + i + 1;
If[r[[idx]] === Null,
r[[idx]] = (x[[1 + i]] -
x[[1 + i + n]])/(\[Rho][x, y, i, n - 1] - \[Rho][x, y,
i + 1, n - 1]) + \[Rho][x, y, i + 1, n - 2];
];
r[[idx]]
]
]
]
Thiele[x_List, y_List, xin_, n_Integer] := Module[{},
If[n > num - 1,
1
,
\[Rho][x, y, 0, n] - \[Rho][x, y, 0, n - 2] + (xin - x[[n + 1]])/
Thiele[x, y, xin, n + 1]
]
]
xval = Range[0, num - 1] step;
funcvals = Sin[xval];
r = ConstantArray[Null, num2];
6 Thiele[funcvals, xval, 0.5, 0]
funcvals = Cos[xval];
r = ConstantArray[Null, num2];
3 Thiele[funcvals, xval, 0.5, 0]
funcvals = Tan[xval];
r = ConstantArray[Null, num2];
4 Thiele[funcvals, xval, 1.0, 0]</syntaxhighlight>
{{out}}
<pre>3.14159
3.14159
3.14159</pre>
 
=={{header|Nim}}==
{{trans|Java}}
<langsyntaxhighlight Nimlang="nim">import strformat
import math
 
Line 948 ⟶ 1,050:
echo fmt"{6 * thiele(tsin, xval, rsin, 0.5, 0):16.14f}"
echo fmt"{3 * thiele(tcos, xval, rcos, 0.5, 0):16.14f}"
echo fmt"{4 * thiele(ttan, xval, rtan, 1.0, 0):16.14f}"</langsyntaxhighlight>
{{out}}
<pre>
Line 959 ⟶ 1,061:
This example shows how the accuracy changes with the degree of interpolation. The table 'columns' are only constructed implicitly during the recursive calculation of <em>rdiff</em> and <em>thiele</em>, but (as mentioned in the C code example) using memoization or explicit tabulation would speed up the calculation. The interpolation uses the nearest points around <em>x</em> for accuracy.
 
<langsyntaxhighlight OCamllang="ocaml">let xv, fv = fst, snd
 
let rec rdiff a l r =
Line 994 ⟶ 1,096:
Printf.printf "4*arctan(1.0) = %.15f\n" (4.0*.(interpolate 1.0 tan_tab n));;
 
List.iter test [8; 12; 16]</langsyntaxhighlight>
Output:
<pre>Degree 8 interpolation:
Line 1,013 ⟶ 1,115:
=={{header|Perl}}==
{{trans|Sidef}}
<langsyntaxhighlight lang="perl">use strict;
use warnings;
use feature 'say';
Line 1,057 ⟶ 1,159:
say 6 * &$sin_inverse(0.5);
say 3 * &$cos_inverse(0.5);
say 4 * &$tan_inverse(1.0);</langsyntaxhighlight>
{{out}}
<pre>3.14159265358979
Line 1,066 ⟶ 1,168:
{{trans|C}}
To be honest I was slightly wary of this, what with tables being passed by reference and fairly heavy use of closures in other languages, but in the end all it took was a simple enum (R_SIN..R_TRIG).
<!--<langsyntaxhighlight Phixlang="phix">-->
<span style="color: #008080;">constant</span> <span style="color: #000000;">N</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">32</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">N2</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">N</span> <span style="color: #0000FF;">*</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">N</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">),</span>
Line 1,116 ⟶ 1,218:
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">,{</span><span style="color: #008000;">"3*thiele(t_cos,xval,R_COS,0.5,0)"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">*</span><span style="color: #000000;">thiele</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t_cos</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xval</span><span style="color: #0000FF;">,</span><span style="color: #000000;">R_COS</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">,{</span><span style="color: #008000;">"4*thiele(t_tan,xval,R_TAN,1,0)"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">4</span><span style="color: #0000FF;">*</span><span style="color: #000000;">thiele</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t_tan</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xval</span><span style="color: #0000FF;">,</span><span style="color: #000000;">R_TAN</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)})</span>
<!--</langsyntaxhighlight>-->
{{out}}
(64 bit, obviously 3 fewer digits on 32 bit)
Line 1,131 ⟶ 1,233:
=={{header|PicoLisp}}==
{{trans|C}}
<langsyntaxhighlight PicoLisplang="picolisp">(scl 17)
(load "@lib/math.l")
 
Line 1,184 ⟶ 1,286:
 
(de iTan (X)
(thiele *TanTable *InvTanTable 1.0 1) )</langsyntaxhighlight>
Test:
<langsyntaxhighlight PicoLisplang="picolisp">(prinl (round (* 6 (iSin 0.5)) 15))
(prinl (round (* 3 (iCos 0.5)) 15))
(prinl (round (* 4 (iTan 1.0)) 15))</langsyntaxhighlight>
Output:
<pre>3.141592653589793
Line 1,195 ⟶ 1,297:
 
=={{header|PowerShell}}==
<langsyntaxhighlight PowerShelllang="powershell">Function Reciprocal-Difference( [Double[][]] $function )
{
$rho=@()
Line 1,263 ⟶ 1,365:
#uncomment to see the function
#"{$atan}"
4*$atan.InvokeReturnAsIs(1)</langsyntaxhighlight>
 
=={{header|Python}}==
{{trans|Go}}
<langsyntaxhighlight lang="python">#!/usr/bin/env python3
 
import math
Line 1,298 ⟶ 1,400:
print('{:16.14f}'.format(6*iSin(.5)))
print('{:16.14f}'.format(3*iCos(.5)))
print('{:16.14f}'.format(4*iTan(1)))</langsyntaxhighlight>
{{out}}
<pre>
Line 1,307 ⟶ 1,409:
 
=={{header|Racket}}==
<langsyntaxhighlight lang="racket">
#lang racket
(define xs (for/vector ([x (in-range 0.0 1.6 0.05)]) x))
Line 1,347 ⟶ 1,449:
(* 3 (icos 0.5))
(* 4 (itan 1.))
</syntaxhighlight>
</lang>
Output:
<langsyntaxhighlight lang="racket">
3.141592653589793
3.1415926535897936
3.1415926535897953
</syntaxhighlight>
</lang>
 
=={{header|Raku}}==
Line 1,360 ⟶ 1,462:
Implemented to parallel the generalized formula, making for clearer, but slower, code. Offsetting that, the use of <code>Promise</code> allows concurrent calculations, so running all three types of interpolation should not take any longer than running just one (presuming available cores).
 
<syntaxhighlight lang="raku" perl6line># reciprocal difference:
multi sub ρ(&f, @x where * < 1) { 0 } # Identity
multi sub ρ(&f, @x where * == 1) { &f(@x[0]) }
Line 1,406 ⟶ 1,508:
say "cos interpolation: $cos_pi";
say "tan interpolation: $tan_pi";
}</langsyntaxhighlight>
 
Output:
Line 1,417 ⟶ 1,519:
 
=={{header|Rust}}==
<langsyntaxhighlight lang="rust">
const N: usize = 32;
const STEP: f64 = 0.05;
Line 1,455 ⟶ 1,557:
}
 
</syntaxhighlight>
</lang>
{{out}}
<pre>
Line 1,465 ⟶ 1,567:
=={{header|Sidef}}==
{{trans|Python}}
<langsyntaxhighlight lang="ruby">func thiele(x, y) {
var ρ = {|i| [y[i]]*(y.len-i) }.map(^y)
 
Line 1,503 ⟶ 1,605:
say 6*iSin(0.5)
say 3*iCos(0.5)
say 4*iTan(1)</langsyntaxhighlight>
{{out}}
<pre>
Line 1,515 ⟶ 1,617:
{{trans|Kotlin}}
 
<langsyntaxhighlight lang="swift">let N = 32
let N2 = N * (N - 1) / 2
let step = 0.05
Line 1,564 ⟶ 1,666:
print(String(format: "%16.14f", 3 * thiele(tcos, xval, &rcos, 0.5, 0)))
print(String(format: "%16.14f", 4 * thiele(ttan, xval, &rtan, 1.0, 0)))
</syntaxhighlight>
</lang>
 
{{out}}
Line 1,575 ⟶ 1,677:
{{works with|Tcl|8.5}}
{{trans|D}}
<langsyntaxhighlight lang="tcl">#
### Create a thiele-interpretation function with the given name that interpolates
### off the given table.
Line 1,618 ⟶ 1,720:
expr {$f1 + ($x - [lindex $X 1]) / ([lindex $rho 1] + $a)}
}} $X [lindex $p 1] [lindex $F 1]
}</langsyntaxhighlight>
Demonstration code:
<langsyntaxhighlight lang="tcl">proc initThieleTest {} {
for {set i 0} {$i < 32} {incr i} {
lappend trigTable(x) [set x [expr {0.05 * $i}]]
Line 1,635 ⟶ 1,737:
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]}]"</langsyntaxhighlight>
Output:
<pre>
Line 1,646 ⟶ 1,748:
{{trans|C}}
{{libheader|Wren-fmt}}
<langsyntaxhighlight ecmascriptlang="wren">import "./fmt" for Fmt
 
var N = 32
Line 1,688 ⟶ 1,790:
Fmt.print("$16.14f", 6 * thiele.call(tsin, xval, rsin, 0.5, 0))
Fmt.print("$16.14f", 3 * thiele.call(tcos, xval, rcos, 0.5, 0))
Fmt.print("$16.14f", 4 * thiele.call(ttan, xval, rtan, 1.0, 0))</langsyntaxhighlight>
 
{{out}}
Line 1,700 ⟶ 1,802:
{{trans|C}}
Please see the C example for the comments I've removed (this is an as pure-as-I-make-it translation).
<langsyntaxhighlight lang="zkl">const N=32, N2=(N * (N - 1) / 2), STEP=0.05;
 
fcn rho(xs,ys,rs, i,n){
Line 1,739 ⟶ 1,841:
print("%16.14f\n".fmt( 6.0 * i_sin(0.5)));
print("%16.14f\n".fmt( 3.0 * i_cos(0.5)));
print("%16.14f\n".fmt( 4.0 * i_tan(1.0)));</langsyntaxhighlight>
{{out}}
<pre>
Jjuanhdez
2,130

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/158239...361527"