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

61,094 bytes added ,  3 months ago
Added FreeBASIC
(Added language PowerShell)
(Added FreeBASIC)
 
(74 intermediate revisions by 30 users not shown)
Line 1:
{{task|Arithmetic operations}}
{{Wikipedia|Thiele's interpolation formula}}
'''[[wp:Thiele's_interpolation_formula|Thiele's interpolation formula]]''' is an interpolation formula for a function ''f''(•) of a single variable. It is expressed as a continued fraction:
 
<br>
:<math> f(x) = f(x_1) + \cfrac{x-x_1}{\rho(x_1,x_2) + \cfrac{x-x_2}{\rho_2(x_1,x_2,x_3) - f(x_1) + \cfrac{x-x_3}{\rho_3(x_1,x_2,x_3,x_4) - \rho(x_1,x_2) + \cdots}}} </math>
'''[[wp:Thiele's_interpolation_formula|Thiele's interpolation formula]]''' is an interpolation formula for a function ''f''(•) of a single variable. &nbsp; It is expressed as a [[continued fraction]]:
 
:: <big><big><math> f(x) = f(x_1) + \cfrac{x-x_1}{\rho_1(x_1,x_2) + \cfrac{x-x_2}{\rho_2(x_1,x_2,x_3) - f(x_1) + \cfrac{x-x_3}{\rho_3(x_1,x_2,x_3,x_4) - \rho_1(x_1,x_2) + \cdots}}} </math></big></big>
ρ represents the [[wp:reciprocal difference|reciprocal difference]], demonstrated here for reference:
 
<big><big><math>\rho</math></big></big> &nbsp; represents the &nbsp; [[wp:reciprocal difference|reciprocal difference]], &nbsp; demonstrated here for reference:
:<math>\rho_1(x_0, x_1) = \frac{x_0 - x_1}{f(x_0) - f(x_1)}</math>
 
:: <big><big><math>\rho_2rho_1(x_0, x_1, x_2) = \frac{x_0 - x_2x_1}{\rho_1f(x_0, x_1) - \rho_1f(x_1, x_2)} + f(x_1)</math></big></big>
 
:: <big><big><math>\rho_nrho_2(x_0, x_1,\ldots,x_n x_2) = \frac{x_0 -x_n x_2}{\rho_{n-1}rho_1(x_0, x_1,\ldots,x_{n-1}) - \rho_{n-1}rho_1(x_1, x_2,\ldots,x_n)} +\rho_{n-2} f(x_1,\ldots,x_{n-1})</math></big></big>
 
:: <big><big><math>\rho_n(x_0,x_1,\ldots,x_n)=\frac{x_0-x_n}{\rho_{n-1}(x_0,x_1,\ldots,x_{n-1})-\rho_{n-1}(x_1,x_2,\ldots,x_n)}+\rho_{n-2}(x_1,\ldots,x_{n-1})</math></big></big>
 
Demonstrate Thiele's interpolation function by:
# Building a &nbsp; '''32''' &nbsp; row ''trig table'' of values of&nbsp; thefor trig&nbsp; functions<big><big><math> ''sin'',x ''cos''</math></big></big> and&nbsp; ''tan''.from e.g.&nbsp; '''for0''' x&nbsp; by &nbsp; '''from0.05''' 0&nbsp; to &nbsp; '''by1.55''' 0.05&nbsp; '''to'''of the trig functions: 1.55...
#* &nbsp; '''sin'''
#* &nbsp; '''cos'''
#* &nbsp; '''tan'''
# Using columns from this table define an inverse - using Thiele's interpolation - for each trig function;
# Finally: demonstrate the following well known trigonometric identities:
#* &nbsp; <big><big> 6 &times; sin<sup>-1</sup> &frac12; = &<math>\pi;</math></big></big>
#* &nbsp; <big><big> 3 &times; cos<sup>-1</sup> &frac12; = &<math>\pi;</math></big></big>
#* &nbsp; <big><big> 4 &times; tan<sup>-1</sup> 1 = &<math>\pi;</math></big></big>
<br><br>
 
=={{header|11l}}==
{{trans|Python}}
 
<syntaxhighlight lang="11l">F thieleInterpolator(x, y)
V ρ = enumerate(y).map((i, yi) -> [yi] * (@y.len - i))
L(i) 0 .< ρ.len - 1
ρ[i][1] = (x[i] - x[i + 1]) / (ρ[i][0] - ρ[i + 1][0])
L(i) 2 .< ρ.len
L(j) 0 .< ρ.len - i
ρ[j][i] = (x[j] - x[j + i]) / (ρ[j][i - 1] - ρ[j + 1][i - 1]) + ρ[j + 1][i - 2]
V ρ0 = ρ[0]
F t(xin)
V a = 0.0
L(i) (@=ρ0.len - 1 .< 1).step(-1)
a = (xin - @=x[i - 1]) / (@=ρ0[i] - @=ρ0[i - 2] + a)
R @=y[0] + (xin - @=x[0]) / (@=ρ0[1] + a)
R t
 
V xVal = (0.<32).map(i -> i * 0.05)
V tSin = xVal.map(x -> sin(x))
V tCos = xVal.map(x -> cos(x))
V tTan = xVal.map(x -> tan(x))
V iSin = thieleInterpolator(tSin, xVal)
V iCos = thieleInterpolator(tCos, xVal)
V iTan = thieleInterpolator(tTan, xVal)
print(‘#.14’.format(6 * iSin(0.5)))
print(‘#.14’.format(3 * iCos(0.5)))
print(‘#.14’.format(4 * iTan(1)))</syntaxhighlight>
 
{{out}}
<pre>
3.14159265358979
3.14159265358979
3.14159265358980
</pre>
 
=={{header|Ada}}==
thiele.ads:
<syntaxhighlight 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;</syntaxhighlight>
 
thiele.adb:
<syntaxhighlight 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;</syntaxhighlight>
 
example:
<syntaxhighlight 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;</syntaxhighlight>
 
output:
<pre>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</pre>
 
=={{header|ALGOL 68}}==
Line 26 ⟶ 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 78 ⟶ 249:
"tan", 4*inv tan(1)
))
)</langsyntaxhighlight>
Output:
<pre>
Line 85 ⟶ 256:
pi estimate using tan interpolation: 3.1415926535898
</pre>
 
=={{header|C}}==
The recursive relations of <math>\rho</math>s can be made clearer: Given <math>N+1</math> sampled points <math>x_0, x_1, \cdots x_N</math>, rewrite the symbol <math>\rho</math> as
{{trans|ALGOL 68}} - Currying, array slicing and exception handling removed, and uses some GCC extensions.
 
<lang c>#include <stdlib.h>
:<math>\rho_{n,i} = \rho_n(x_i, x_{i+1}, \cdots x_{i+n})</math> where <math>1 \leq n \leq N</math>,
#include <stdio.h>
 
with suplements
 
:<math>\rho_{0, i} = f(x_i)\qquad\text{and}\qquad\rho_{n, i} = 0\quad\text{for}\quad n < 0</math>.
 
Now the recursive relation is simply
 
:<math>\rho_{n,i} = \displaystyle{x_i - x_{i + n} \over \rho_{n-1,i} - \rho_{n-1,i+1}} + \rho_{n-2,i+1}, \ 0\leq i+n \leq N</math>
 
Also note how <math>f(x_1)</math> in the interpolation formula can be replaced by <math>\rho_{0,1}</math>; define Thiele interpolation at step <math>n</math> as
 
:<math>\displaystyle{F_n(x) = \rho_{n,1} - \rho_{n - 2, 1} + { x - x_{n+1}\over F_{n+1}(x)}}</math>
 
with the termination <math>F_N(x) = 1</math>, and the interpolation formula is now <math>f(x) = F_0(x)</math>, easily implemented as a recursive function.
 
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.
<syntaxhighlight lang="c">#include <stdio.h>
#include <string.h>
#include <math.h>
 
#define N 32
/* The MODE of lx and ly here should really be a UNION of "something double",
#define N2 (N * (N - 1) / 2)
"something COMPLex", and "something SYMBOLIC" ... */
#define STEP .05
 
double thiele(int upb_lx, double lx[], double ly[], double x)
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;
int upb_xx=upb_lx+1;
if (!n) return y[i];
double *xx=lx-1, *yy=ly-1; /* shift base index to 1 */
 
int n=upb_xx;
int idx = (N - 1 - n) * (N - n) / 2 + i;
/* Assuming that the values of xx are distinct ... */
if (r[idx] != r[idx]) /* only happens if value not computed yet */
double p [/*0:*/n-1 +1][/*1:*/n +1];
r[idx] = (x[i] - x[i + n])
int i, j;
/ (rho(x, y, r, i, n - 1) - rho(x, y, r, i + 1, n - 1))
for(i=1; i<=upb_xx; i++) p[0][i]=yy[i];
+ rho(x, y, r, i + 1, n - 2);
for(i=1; i<=n-1; i++)p[1][i]=(xx[i]-xx[1+i])/(p[0][i]-p[0][1+i]);
return r[idx];
for(i=2; i<=n-1; i++){
for(j=1; j<=n-i; j++){
p[i][j]=(xx[j]-xx[j+i])/(p[i-1][j]-p[i-1][j+1])+p[i-2][j+1];
}
}
double a=0;
for(i=n-1; i>=2; i--){ a=(x-xx[i])/(p[i][1]-p[i-2][1]+a); }
return yy[1]+(x-xx[1])/(p[1][1]+a);
}
main(){
 
double thiele(double *x, double lwb_x=0*y, upb_x=1.55double *r, delta_xdouble =xin, 0.05;int n)
{
if (n > N - 1) return 1;
int upb_trig_table = ((upb_x-lwb_x)/delta_x);
return rho(x, y, r, 0, n) - rho(x, y, r, 0, n - 2)
typedef double trig_table_t [/*0:*/upb_trig_table +1];
+ (xin - x[n]) / thiele(x, y, r, xin, n + 1);
trig_table_t x_OF_trig_table, sin_x_OF_trig_table, cos_x_OF_trig_table, tan_x_OF_trig_table;
}
void init_trig_table(){
int i;
for(i = 0; i<= upb_trig_table; i++){
double x = i*delta_x;
x_OF_trig_table[i]=x;
sin_x_OF_trig_table[i]=sin(x);
cos_x_OF_trig_table[i]=cos(x);
tan_x_OF_trig_table[i]=tan(x);
}
}
init_trig_table();
int upb_row = 4;
double inv_sin(double x){ return thiele(upb_trig_table, sin_x_OF_trig_table, x_OF_trig_table, x);}
double inv_cos(double x){ return thiele(upb_trig_table, cos_x_OF_trig_table, x_OF_trig_table, x);}
double inv_tan(double x){ return thiele(upb_trig_table, tan_x_OF_trig_table, x_OF_trig_table, x);}
char *result_fmt = "pi estimate using %s interpolation: %0.13f\n";
printf(result_fmt, "sin", 6*inv_sin(1.0/2));
printf(result_fmt, "cos", 3*inv_cos(1.0/2));
printf(result_fmt, "tan", 4*inv_tan(1));
}</lang>
Output:
<pre>
pi estimate using sin interpolation: 3.1415926535898
pi estimate using cos interpolation: 3.1415926535898
pi estimate using tan interpolation: 3.1415926535898
</pre>
 
#define i_sin(x) thiele(t_sin, xval, r_sin, x, 0)
=={{header|D}}==
#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;
}</syntaxhighlight>
{{out}}
<pre>3.14159265358979
3.14159265358979
3.14159265358979</pre>
 
=={{header|C++}}==
{{trans|C}}
{{works with|C++14}}
<syntaxhighlight lang="cpp">#include <cmath>
#include <iostream>
#include <iomanip>
#include <string.h>
 
constexpr unsigned int N = 32u;
double xval[N], t_sin[N], t_cos[N], t_tan[N];
 
constexpr unsigned int N2 = N * (N - 1u) / 2u;
double r_sin[N2], r_cos[N2], r_tan[N2];
 
double ρ(double *x, double *y, double *r, int i, int n) {
<lang d>import std.stdio ;
if (n < 0)
import std.math ;
return 0;
if (!n)
return y[i];
 
unsigned int idx = (N - 1 - n) * (N - n) / 2 + i;
U[] myMap(U,V) (U delegate(V) f, V[] a) {
Vif (r[idx] r != new Vr[idx](a.length) ;
r[idx] = (x[i] - x[i + n]) / (ρ(x, y, r, i, n - 1) - ρ(x, y, r, i + 1, n - 1)) + ρ(x, y, r, i + 1, n - 2);
foreach(i, v ; a)
return r[iidx] = f(v) ;
return r ;
}
 
double thiele(double *x, double *y, double *r, double xin, unsigned int n) {
T[] myRng(T)(T start, T end, T step) {
return n > N - 1 ? 1. : ρ(x, y, r, 0, n) - ρ(x, y, r, 0, n - 2) + (xin - x[n]) / thiele(x, y, r, xin, n + 1);
T[] r ;
for(; start < end ; start += step)
r ~= start ;
return r ;
}
 
inline auto i_sin(double x) { return thiele(t_sin, xval, r_sin, x, 0); }
alias real delegate(real) RealFun ;
inline auto i_cos(double x) { return thiele(t_cos, xval, r_cos, x, 0); }
const real End = 1.55L , Step = 0.05L , Start = 0.0L ;
inline auto i_tan(double x) { return thiele(t_tan, xval, r_tan, x, 0); }
 
classint Thielemain() {
constexpr double step = .05;
for (auto i = 0u; 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 (auto i = 0u; i < N2; i++)
r_sin[i] = r_cos[i] = r_tan[i] = NAN;
 
std::cout << std::setw(16) << std::setprecision(25)
private real[] F, X, RhoY, RhoX ;
<< 6 * i_sin(.5) << std::endl
const int NN ;
<< 3 * i_cos(.5) << std::endl
<< 4 * i_tan(1.) << std::endl;
this(real[] f, real[] x) {
 
F = f.dup ;
return X = x.dup 0;
}</syntaxhighlight>
NN = X.length ;
{{out}}
assert(NN > 2, "at leat 3 values") ;
<pre>3.141592653589793115997963
assert(X.length == F.length, "input arrays not of same size") ;
3.141592653589793560087173
RhoY = rhoN(F, X) ;
3.141592653589794892354803</pre>
RhoX = rhoN(X, F) ;
 
=={{header|Common Lisp}}==
Using the notations from above the C code instead of task desc.
<syntaxhighlight lang="lisp">;; 256 is heavy overkill, but hey, we memoized
(defparameter *thiele-length* 256)
(defparameter *rho-cache* (make-hash-table :test #'equal))
 
(defmacro make-thele-func (f name xx0 xx1)
(let ((xv (gensym)) (yv (gensym))
(x0 (gensym)) (x1 (gensym)))
`(let* ((,xv (make-array (1+ *thiele-length*)))
(,yv (make-array (1+ *thiele-length*)))
(,x0 ,xx0)
(,x1 ,xx1))
(loop for i to *thiele-length* with x do
(setf x (+ ,x0 (* (/ (- ,x1 ,x0) *thiele-length*) i))
(aref ,yv i) x
(aref ,xv i) (funcall ,f x)))
(defun ,name (x) (thiele x ,yv ,xv, 0)))))
 
(defun rho (yv xv n i)
(let (hit (key (list yv xv n i)))
(if (setf hit (gethash key *rho-cache*))
hit
(setf (gethash key *rho-cache*)
(cond ((zerop n) (aref yv i))
((minusp n) 0)
(t (+ (rho yv xv (- n 2) (1+ i))
(/ (- (aref xv i)
(aref xv (+ i n)))
(- (rho yv xv (1- n) i)
(rho yv xv (1- n) (1+ i)))))))))))
 
(defun thiele (x yv xv n)
(if (= n *thiele-length*)
1
(+ (- (rho yv xv n 1) (rho yv xv (- n 2) 1))
(/ (- x (aref xv (1+ n)))
(thiele x yv xv (1+ n))))))
 
(make-thele-func #'sin inv-sin 0 (/ pi 2))
(make-thele-func #'cos inv-cos 0 (/ pi 2))
(make-thele-func #'tan inv-tan 0 (/ pi 2.1)) ; tan(pi/2) is INF
 
(format t "~f~%" (* 6 (inv-sin .5)))
(format t "~f~%" (* 3 (inv-cos .5)))
(format t "~f~%" (* 4 (inv-tan 1)))</syntaxhighlight>output (SBCL):<syntaxhighlight lang="text">3.141592653589793
3.1415926535885172
3.141592653589819</syntaxhighlight>
 
=={{header|D}}==
<syntaxhighlight lang="d">import std.stdio, std.range, std.array, std.algorithm, std.math;
 
struct Domain {
const real b, e, s;
 
auto range() const pure /*nothrow*/ @safe /*@nogc*/ {
return iota(b, e + s, s);
}
}
 
real thiseval0(RealFunalias funRY, realalias startX, =alias Start,Y)(in real endx) =pure End,nothrow real@safe step = Step)@nogc {
real a = 0.0L;
this(myMap(fun, myRng(start, end, step)),
foreach_reverse (immutable i; 2 .. X.length - 3)
myRng(start, end, step)) ;
a = (x - X[i]) / (RY[i] - RY[i-2] + a);
return Y[1] + (x - X[1]) / (RY[1] + a);
}
 
immutable struct Thiele {
immutable real[] Y, X, rhoY, rhoX;
 
this(real[] y, real[] x) immutable pure nothrow /*@safe*/
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);
}
 
private staticthis(in real[] rhoNfunction(real[]) f,pure real[]nothrow x)@safe {@nogc f,
int NDomain d = xDomain(0.0L, 1.55L, 0.length;05L))
immutable pure /*nothrow @safe*/ {
real[][] p = new real[][] (N, N) ;
for(intauto ixrng = 0d.range.array; i < N ; i++)
p[i][0] = this(xrng.map!f[i].array, xrng);
for(int i = 0; i < N - 1 ; i++)
p[i][1] = (x[i] - x[i+1]) / (p[i][0] - p[i+1][0]) ;
for(int j = 2; j < N - 1; j++)
for(int i = 0; i < N - j - 1 ; i++)
p[i][j] = p[i+1][j-2] + (x[i] - x[i+j]) /
(p[i][j-1] - p[i+1][j-1]) ; ;
return p[1].dup ;
}
 
realauto evalYrhoN(immutable real[] y, immutable real[] x) {
pure nothrow @safe {
real a = 0.0L ;
for(immutable int iN = NN - 3x.length; i >= 2 ; i--)
auto ap = (xnew - Xreal[i]) / (RhoY[i](N, - RhoY[i-2] + aN) ;
return Fp[10] + (x - X[1]) /= (RhoYy[1] + a) ;
p[1][0 .. $ - 1] = (x[0 .. $-1] - x[1 .. $]) /
}
(p[0][0 .. $-1] - p[0][1 .. $]);
real evalX(real y) { // inverse
realforeach a(immutable =int j; 2 0.0L. ;N - 1) {
for(int i = NN -immutable 3;M i >= 2N ;- j i--) 1;
ap[j][0..M] = (y p[j- F2][i1..M+1]) /+ (RhoXx[i0..M] - RhoXx[i-2] j..M+ aj]) ;/
return X[1] + (y - F[1]) / (RhoX[1] + a) ; (p[j-1][0 .. M] - p[j-1][1 .. M+1]);
}
return p.map!q{ a[1] }.array;
}
 
alias eval = eval0!(rhoY, X, Y);
alias inverse = eval0!(rhoX, Y, X);
}
 
void main() {
// Can't pass sin, cos and tan directly.
auto fun = (real x) { return std.math.sin(x) ; } ;
autoimmutable ttsin = new Thiele(fun)x => x.sin);
writefln("immutable %dtcos interpolating= points",Thiele(x t=> x.NNcos) ;
immutable ttan = Thiele(x => x.tan);
writefln("std.math.sin(0.5) : %20.18f", std.math.sin(0.5L)) ;
writefln(" Thiele sin(0.5) : %20.18f", t.evalY(0.5L)) ;
writefln("*%20.18f library constant", std.math.PI) ;
writefln(" %20.18f 6*sin`(0.5)", t.evalX(0.5L) * 6.0L) ;
t = new Thiele((real x) {return std.math.cos(x) ; }) ;
writefln(" %20.18f 3*cos`(0.5)", t.evalX(0.5L) * 3.0L) ;
t = new Thiele((real x) {return std.math.tan(x) ; }) ;
writefln(" %20.18f 4*tan`(0.5)", t.evalX(1.0L) * 4.0L) ;
}</lang>
 
writefln(" %d interpolating points\n", tsin.X.length);
output:
writefln("std.math.sin(0.5): %20.18f", 0.5L.sin);
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);
}</syntaxhighlight>
{{out}}
<pre> 32 interpolating points
std.math.sin(0.5) : 0.479425538604203000
Thiele sin(0.5) : 0.479425538604203000
*3.141592653589793239 library constant
3.141592653589793238 6*sin`(0.5)
3.141592653589793238 3*cos`(0.5)
3.141592653589793238 4*tan`(0.5)</pre>
 
std.math.sin(0.5): 0.479425538604203000
=={{header|Perl 6}}==
Thiele sin(0.5): 0.479425538604203000
{{works with|Rakudo|2010.09-32}}<br>
Implemented to parallel the (generalized) formula. (i.e. clearer, but naive and very slow.)
<lang perl6>use v6;
 
*3.1415926535897932385 library constant
# reciprocal difference:
3.1415926535897932380 6 * inv_sin(0.5)
multi sub rho($f, @x where { +@x < 1 }) { 0 } # Identity
3.1415926535897932382 3 * inv_cos(0.5)
multi sub rho($f, @x where { +@x == 1 }) { $f(@x[0]) }
3.1415926535897932382 4 * inv_tan(1.0)</pre>
multi sub rho($f, @x where { +@x > 1 }) {
 
my $ord = +@x;
=={{header|FreeBASIC}}==
{{trans|Phix}}
return
<syntaxhighlight lang="vbnet">Const As Integer n1 = 32
( @x[0] - @x[* -1] ) # ( x - x[n] )
Const As Integer n2 = (n1 * (n1 - 1) / 2)
/ ( rho($f, @x[^($ord -1)]) # / ( rho[n-1](x[0], ..., x[n-1])
Const As Double paso = 0.05
- rho($f, @x[1..^($ord)]) ) # - rho[n-1](x[1], ..., x[n]) )
Const As Double INF = 1e308
+ rho($f, @x[1..^($ord -1)]); # + rho[n-2](x[1], ..., x[n-1])
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}}
<syntaxhighlight lang="go">package main
 
import (
"fmt"
"math"
)
 
func main() {
// task 1: build 32 row trig table
const nn = 32
const step = .05
xVal := make([]float64, nn)
tSin := make([]float64, nn)
tCos := make([]float64, nn)
tTan := make([]float64, nn)
for i := range xVal {
xVal[i] = float64(i) * step
tSin[i], tCos[i] = math.Sincos(xVal[i])
tTan[i] = tSin[i] / tCos[i]
}
// task 2: define inverses
iSin := thieleInterpolator(tSin, xVal)
iCos := thieleInterpolator(tCos, xVal)
iTan := thieleInterpolator(tTan, xVal)
// task 3: demonstrate identities
fmt.Printf("%16.14f\n", 6*iSin(.5))
fmt.Printf("%16.14f\n", 3*iCos(.5))
fmt.Printf("%16.14f\n", 4*iTan(1))
}
 
func thieleInterpolator(x, y []float64) func(float64) float64 {
# Thiele:
n := len(x)
multi sub thiele($x, %f, $ord where { $ord == +%f }) { 1 } # Identity
ρ := make([][]float64, n)
multi sub thiele($x, %f, $ord) {
for i := range ρ {
my $f = {%f{$^a}}; # f(x) as a table lookup
ρ[i] = make([]float64, n-i)
ρ[i][0] = y[i]
# Caveat: depends on the fact that Rakudo maintains key order within hashes
}
my $a = rho($f, %f.keys[^($ord +1)]);
for i := 0; i < n-1; i++ {
my $b = rho($f, %f.keys[^($ord -1)]);
ρ[i][1] = (x[i] - x[i+1]) / (ρ[i][0] - ρ[i+1][0])
}
my $num = $x - %f.keys[$ord];
for i := 2; i < n; i++ {
my $cont = thiele($x, %f, $ord +1);
for j := 0; j < n-i; j++ {
ρ[j][i] = (x[j]-x[j+i])/(ρ[j][i-1]-ρ[j+1][i-1]) + ρ[j+1][i-2]
# Thiele always takes this form:
return $a - $b + ( $num / $cont );}
}
// ρ0 used in closure. the rest of ρ becomes garbage.
ρ0 := ρ[0]
return func(xin float64) float64 {
var a float64
for i := n - 1; i > 1; i-- {
a = (xin - x[i-1]) / (ρ0[i] - ρ0[i-2] + a)
}
return y[0] + (xin-x[0])/(ρ0[1]+a)
}
}</syntaxhighlight>
Output:
<pre>
3.14159265358979
3.14159265358979
3.14159265358980
</pre>
 
=={{header|Haskell}}==
Caching of rho is automatic due to lazy lists.
<syntaxhighlight lang="haskell">thiele :: [Double] -> [Double] -> Double -> Double
thiele xs ys = f rho1 (tail xs)
where
f _ [] _ = 1
f r@(r0:r1:r2:rs) (x:xs) v = r2 - r0 + (v - x) / f (tail r) xs v
rho1 = (!! 1) . (++ [0]) <$> rho
rho = repeat 0 : repeat 0 : ys : rnext (tail rho) xs (tail xs)
where
rnext _ _ [] = []
rnext r@(r0:r1:rs) x xn =
let z_ = zipWith
in z_ (+) (tail r0) (z_ (/) (z_ (-) x xn) (z_ (-) r1 (tail r1))) :
rnext (tail r) x (tail xn)
 
-- Inverted interpolation function of f
invInterp :: (Double -> Double) -> [Double] -> Double -> Double
invInterp f xs = thiele (map f xs) xs
 
main :: IO ()
main =
mapM_
print
[ 3.21 * inv_sin (sin (pi / 3.21))
, pi / 1.2345 * inv_cos (cos 1.2345)
, 7 * inv_tan (tan (pi / 7))
]
where
[inv_sin, inv_cos, inv_tan] =
uncurry ((. div_pi) . invInterp) <$>
[(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]</syntaxhighlight>
{{out}}
<pre>3.141592653589795
3.141592653589791
3.141592653587783</pre>
 
=={{header|J}}==
<syntaxhighlight lang="j">
span =: {. - {: NB. head - tail
spans =: span\ NB. apply span to successive infixes
</syntaxhighlight>
 
<pre>
span 12 888 6 4 8 3
9
4 spans 12 888 6 4 8 3
8 880 3
</pre>
 
<syntaxhighlight 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
)
</syntaxhighlight>
 
<pre>
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│
└─────────┴────────────┘
</pre>
 
<syntaxhighlight 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
)
</syntaxhighlight>
 
<pre>
'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
</pre>
 
=={{header|Java}}==
{{trans|C}}
<syntaxhighlight lang="java">import static java.lang.Math.*;
 
public class Test {
final static int N = 32;
final static int N2 = (N * (N - 1) / 2);
final static double STEP = 0.05;
 
static double[] xval = new double[N];
static double[] t_sin = new double[N];
static double[] t_cos = new double[N];
static double[] t_tan = new double[N];
 
static double[] r_sin = new double[N2];
static double[] r_cos = new double[N2];
static double[] r_tan = new double[N2];
 
static double rho(double[] x, double[] y, double[] r, int i, int n) {
if (n < 0)
return 0;
 
if (n == 0)
return y[i];
 
int idx = (N - 1 - n) * (N - n) / 2 + i;
if (r[idx] != r[idx])
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];
}
 
static 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);
}
 
public static void main(String[] args) {
for (int 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 (int i = 0; i < N2; i++)
r_sin[i] = r_cos[i] = r_tan[i] = Double.NaN;
 
System.out.printf("%16.14f%n", 6 * thiele(t_sin, xval, r_sin, 0.5, 0));
System.out.printf("%16.14f%n", 3 * thiele(t_cos, xval, r_cos, 0.5, 0));
System.out.printf("%16.14f%n", 4 * thiele(t_tan, xval, r_tan, 1.0, 0));
}
}</syntaxhighlight>
<pre>3.14159265358979
3.14159265358979
3.14159265358980</pre>
 
=={{header|Julia}}==
Accuracy improves with a larger table and smaller step size.
{{trans|C}}
<syntaxhighlight lang="julia">const N = 256
const N2 = N * div(N - 1, 2)
const step = 0.01
const xval_table = zeros(Float64, N)
const tsin_table = zeros(Float64, N)
const tcos_table = zeros(Float64, N)
const ttan_table = zeros(Float64, N)
const rsin_cache = Dict{Float64, Float64}()
const rcos_cache = Dict{Float64, Float64}()
const rtan_cache = Dict{Float64, Float64}()
 
function rho(x, y, rhocache, i, n)
if n < 0
return 0.0
elseif n == 0
return y[i+1]
end
idx = (N - 1 - n) * div(N - n, 2) + i
if !haskey(rhocache, idx)
rhocache[idx] = (x[i+1] - x[i + n+1]) / (rho(x, y, rhocache, i, n - 1) -
rho(x, y, rhocache, i + 1, n - 1)) + rho(x, y, rhocache, i + 1, n - 2)
end
rhocache[idx]
end
 
function thiele(x, y, r, xin, n)
if n > N - 1
return 1.0
end
rho(x, y, r, 0, n) - rho(x, y, r, 0, n - 2) + (xin - x[n + 1]) / thiele(x, y, r, xin, n + 1)
end
 
function thiele_tables()
for i in 1:N
xval_table[i] = (i-1) * step
tsin_table[i] = sin(xval_table[i])
tcos_table[i] = cos(xval_table[i])
ttan_table[i] = tsin_table[i] / tcos_table[i]
end
println(6 * thiele(tsin_table, xval_table, rsin_cache, 0.5, 0))
println(3 * thiele(tcos_table, xval_table, rcos_cache, 0.5, 0))
println(4 * thiele(ttan_table, xval_table, rtan_cache, 1.0, 0))
end
 
thiele_tables()
</syntaxhighlight>{{output}}<pre>
3.1415926535898335
3.141592653589818
3.141592653589824
</pre>
 
=={{header|Kotlin}}==
{{trans|C}}
<syntaxhighlight lang="scala">// version 1.1.2
 
const val N = 32
const val N2 = N * (N - 1) / 2
const val STEP = 0.05
 
val xval = DoubleArray(N)
val tsin = DoubleArray(N)
val tcos = DoubleArray(N)
val ttan = DoubleArray(N)
val rsin = DoubleArray(N2) { Double.NaN }
val rcos = DoubleArray(N2) { Double.NaN }
val rtan = DoubleArray(N2) { Double.NaN }
 
fun rho(x: DoubleArray, y: DoubleArray, r: DoubleArray, i: Int, n: Int): Double {
if (n < 0) return 0.0
if (n == 0) return y[i]
val idx = (N - 1 - n) * (N - n) / 2 + i
if (r[idx].isNaN()) {
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]
}
 
fun thiele(x: DoubleArray, y: DoubleArray, r: DoubleArray, xin: Double, n: Int): Double {
## Demo
if (n > N - 1) return 1.0
sub mk-inv($fn, $d, $lim) {
return rho(x, y, r, 0, n) - rho(x, y, r, 0, n - 2) +
my %h;
for 0..$lim { %h{ $fn ($_xin *- $dx[n]) }/ =thiele(x, $_y, *r, $dxin, }n + 1)
return %h;
}
 
fun main(args: Array<String>) {
sub MAIN($tblsz) {
my %invsin =for mk-inv(&sin,i in 0.05, $tblszuntil N); {
xval[i] = i * STEP
my %invcos = mk-inv(&cos, 0.05, $tblsz);
tsin[i] = Math.sin(xval[i])
my %invtan = mk-inv(&tan, 0.05, $tblsz);
tcos[i] = Math.cos(xval[i])
ttan[i] = tsin[i] / tcos[i]
}
println("%16.14f".format(6 * thiele(tsin, xval, rsin, 0.5, 0)))
println("%16.14f".format(3 * thiele(tcos, xval, rcos, 0.5, 0)))
println("%16.14f".format(4 * thiele(ttan, xval, rtan, 1.0, 0)))
}</syntaxhighlight>
 
{{out}}
<pre>
3.14159265358979
3.14159265358979
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}}
<syntaxhighlight lang="nim">import strformat
import math
 
const N = 32
const N2 = N * (N - 1) div 2
const STEP = 0.05
 
var xval = newSeq[float](N)
var tsin = newSeq[float](N)
var tcos = newSeq[float](N)
var ttan = newSeq[float](N)
var rsin = newSeq[float](N2)
var rcos = newSeq[float](N2)
var rtan = newSeq[float](N2)
 
proc rho(x, y: openArray[float], r: var openArray[float], i, n: int): float =
if n < 0:
return 0
if n == 0:
return y[i]
mylet $sin_piidx = 6(N - 1 - n) * thiele(0.5,N %invsin,- 0n); div 2 + i
if r[idx] != r[idx]:
my $cos_pi = 3 * thiele(0.5, %invcos, 0);
my $tan_pi r[idx] = 4(x[i] *- thiele(1.0,x[i %invtan,+ 0n]); /
(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]
 
proc thiele(x, y: openArray[float], r: var openArray[float], xin: float, n: int): float =
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)
 
for i in 0..<N:
xval[i] = float(i) * STEP
tsin[i] = sin(xval[i])
tcos[i] = cos(xval[i])
ttan[i] = tsin[i] / tcos[i]
 
for i in 0..<N2:
rsin[i] = NaN
rcos[i] = NaN
rtan[i] = NaN
echo fmt"{6 * thiele(tsin, xval, rsin, 0.5, 0):16.14f}"
say "pi = {pi}";
echo fmt"{3 * thiele(tcos, xval, rcos, 0.5, 0):16.14f}"
say "estimations using a table of $tblsz elements:";
echo fmt"{4 * thiele(ttan, xval, rtan, 1.0, 0):16.14f}"</syntaxhighlight>
say "sin interpolation: $sin_pi";
{{out}}
say "cos interpolation: $cos_pi";
<pre>
say "tan interpolation: $tan_pi";
3.14159265358979
}</lang>
3.14159265358979
3.14159265358979
</pre>
 
=={{header|OCaml}}==
Output (table size of 6 for want of resources):
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.
 
<syntaxhighlight lang="ocaml">let xv, fv = fst, snd
 
let rec rdiff a l r =
if l > r then 0.0 else
if l = r then fv a.(l) else
if l+1 = r then (xv a.(l) -. xv a.(r)) /. (fv a.(l) -. fv a.(r)) else
(xv a.(l) -. xv a.(r)) /. (rdiff a l (r-1) -. rdiff a (l+1) r) +. rdiff a (l+1) (r-1)
 
let rec thiele x a a0 k n =
if k = n then 1.0 else
rdiff a a0 (a0+k) -. rdiff a a0 (a0+k-2) +. (x -. xv a.(a0+k)) /. thiele x a a0 (k+1) n
 
let interpolate x a n =
let m = Array.length a in
let dist i = abs_float (x -. xv a.(i)) in
let nearer i j = if dist j < dist i then j else i in
let rec closest i j = if j = m then i else closest (nearer i j) (j+1) in
let c = closest 0 1 in
let c' = if c < n/2 then 0 else if c > m-n then m-n else c-(n/2) in
thiele x a c' 0 n
 
let table a b n f =
let g i =
let x = a +. (b-.a)*.(float i)/.(float (n-1)) in
(f x, x) in
Array.init n g
 
let [sin_tab; cos_tab; tan_tab] = List.map (table 0.0 1.55 32) [sin; cos; tan]
 
let test n =
Printf.printf "\nDegree %d interpolation:\n" n;
Printf.printf "6*arcsin(0.5) = %.15f\n" (6.0*.(interpolate 0.5 sin_tab n));
Printf.printf "3*arccos(0.5) = %.15f\n" (3.0*.(interpolate 0.5 cos_tab n));
Printf.printf "4*arctan(1.0) = %.15f\n" (4.0*.(interpolate 1.0 tan_tab n));;
 
List.iter test [8; 12; 16]</syntaxhighlight>
Output:
<pre>Degree 8 interpolation:
6*arcsin(0.5) = 3.141592654456238
3*arccos(0.5) = 3.141592653520809
4*arctan(1.0) = 3.141592653437432
 
Degree 12 interpolation:
6*arcsin(0.5) = 3.141592653587590
3*arccos(0.5) = 3.141592653562618
4*arctan(1.0) = 3.141592653589756
 
Degree 16 interpolation:
6*arcsin(0.5) = 3.141592653589793
3*arccos(0.5) = 3.141592653589793
4*arctan(1.0) = 3.141592653589793</pre>
 
=={{header|Perl}}==
{{trans|Sidef}}
<syntaxhighlight lang="perl">use strict;
use warnings;
use feature 'say';
use Math::Trig;
use utf8;
 
sub thiele {
my($x, $y) = @_;
 
my @ρ;
push @ρ, [($$y[$_]) x (@$y-$_)] for 0 .. @$y-1;
for my $i (0 .. @ρ - 2) {
$ρ[$i][1] = (($$x[$i] - $$x[$i+1]) / ($ρ[$i][0] - $ρ[$i+1][0]))
}
for my $i (2 .. @ρ - 2) {
for my $j (0 .. (@ρ - 2) - $i) {
$ρ[$j][$i] = ((($$x[$j]-$$x[$j+$i]) / ($ρ[$j][$i-1]-$ρ[$j+1][$i-1])) + $ρ[$j+1][$i-2])
}
}
my @ρ0 = @{$ρ[0]};
 
return sub {
my($xin) = @_;
 
my $a = 0;
for my $i (reverse 2 .. @ρ0 - 2) {
$a = (($xin - $$x[$i-1]) / ($ρ0[$i] - $ρ0[$i-2] + $a))
}
$$y[0] + (($xin - $$x[0]) / ($ρ0[1] + $a))
}
}
 
my(@x,@sin_table,@cos_table,@tan_table);
push @x, .05 * $_ for 0..31;
push @sin_table, sin($_) for @x;
push @cos_table, cos($_) for @x;
push @tan_table, tan($_) for @x;
 
my $sin_inverse = thiele(\@sin_table, \@x);
my $cos_inverse = thiele(\@cos_table, \@x);
my $tan_inverse = thiele(\@tan_table, \@x);
 
say 6 * &$sin_inverse(0.5);
say 3 * &$cos_inverse(0.5);
say 4 * &$tan_inverse(1.0);</syntaxhighlight>
{{out}}
<pre>3.14159265358979
3.14159265358979
3.1415926535898</pre>
 
=={{header|Phix}}==
{{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).
<!--<syntaxhighlight lang="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>
<span style="color: #000000;">STEP</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0.05</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">inf</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1e300</span><span style="color: #0000FF;">*</span><span style="color: #000000;">1e300</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">nan</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-(</span><span style="color: #000000;">inf</span><span style="color: #0000FF;">/</span><span style="color: #000000;">inf</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">xval</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">t_sin</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">t_cos</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">t_tan</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">N</span><span style="color: #0000FF;">),</span><span style="color: #000000;">4</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">N</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">xval</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">i</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;">STEP</span>
<span style="color: #000000;">t_sin</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sin</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xval</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span>
<span style="color: #000000;">t_cos</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">cos</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xval</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span>
<span style="color: #000000;">t_tan</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">t_sin</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">t_cos</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">enum</span> <span style="color: #000000;">R_SIN</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">R_COS</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">R_TAN</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">R_TRIG</span><span style="color: #0000FF;">=$</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">rhot</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">nan</span><span style="color: #0000FF;">,</span><span style="color: #000000;">N2</span><span style="color: #0000FF;">),</span><span style="color: #000000;">R_TRIG</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">rho</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">rdx</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">int</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">int</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">n</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">idx</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: #000000;">n</span><span style="color: #0000FF;">)</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;">n</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">2</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">i</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">;</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">rhot</span><span style="color: #0000FF;">[</span><span style="color: #000000;">rdx</span><span style="color: #0000FF;">][</span><span style="color: #000000;">idx</span><span style="color: #0000FF;">]=</span><span style="color: #000000;">nan</span> <span style="color: #008080;">then</span> <span style="color: #000080;font-style:italic;">-- value not computed yet</span>
<span style="color: #000000;">rhot</span><span style="color: #0000FF;">[</span><span style="color: #000000;">rdx</span><span style="color: #0000FF;">][</span><span style="color: #000000;">idx</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</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;">x</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</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: #0000FF;">(</span><span style="color: #000000;">rho</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rdx</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">i</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;">rho</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rdx</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</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;">rho</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rdx</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">rhot</span><span style="color: #0000FF;">[</span><span style="color: #000000;">rdx</span><span style="color: #0000FF;">][</span><span style="color: #000000;">idx</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">thiele</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">rdx</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">xin</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">n</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: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">rho</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rdx</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</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;">rho</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rdx</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
<span style="color: #0000FF;">+</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">xin</span><span style="color: #0000FF;">-</span><span style="color: #000000;">x</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;">thiele</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rdx</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">xin</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: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">fmt</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">machine_bits</span><span style="color: #0000FF;">()=</span><span style="color: #000000;">32</span><span style="color: #0000FF;">?</span><span style="color: #008000;">"%32s : %.14f\n"</span>
<span style="color: #0000FF;">:</span><span style="color: #008000;">"%32s : %.17f\n"</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;">"PI"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">PI</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;">"6*arcsin(0.5)"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">6</span><span style="color: #0000FF;">*</span><span style="color: #7060A8;">arcsin</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0.5</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;">"3*arccos(0.5)"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">*</span><span style="color: #7060A8;">arccos</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0.5</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*arctan(1)"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">4</span><span style="color: #0000FF;">*</span><span style="color: #7060A8;">arctan</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</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;">"6*thiele(t_sin,xval,R_SIN,0.5,0)"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">6</span><span style="color: #0000FF;">*</span><span style="color: #000000;">thiele</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t_sin</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xval</span><span style="color: #0000FF;">,</span><span style="color: #000000;">R_SIN</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;">"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>
<!--</syntaxhighlight>-->
{{out}}
(64 bit, obviously 3 fewer digits on 32 bit)
<pre>
PI : 3.14159265358979324
pi = 3.14159265358979
6*arcsin(0.5) : 3.14159265358979324
estimations using a table of 6 elements:
3*arccos(0.5) : 3.14159265358979324
sin interpolation: 3.14153363985515
4*arctan(1) : 3.14159265358979324
cos interpolation: 1.68779321655997
6*thiele(t_sin,xval,R_SIN,0.5,0) : 3.14159265358979324
tan interpolation: 3.14826236377727
3*thiele(t_cos,xval,R_COS,0.5,0) : 3.14159265358979324
4*thiele(t_tan,xval,R_TAN,1,0) : 3.14159265358979324
</pre>
 
=={{header|PicoLisp}}==
{{trans|C}}
<syntaxhighlight lang="picolisp">(scl 17)
(load "@lib/math.l")
 
(setq
*X-Table (range 0.0 1.55 0.05)
*SinTable (mapcar sin *X-Table)
*CosTable (mapcar cos *X-Table)
*TanTable (mapcar tan *X-Table)
*TrigRows (length *X-Table) )
 
(let N2 (>> 1 (* *TrigRows (dec *TrigRows)))
(setq
*InvSinTable (need N2)
*InvCosTable (need N2)
*InvTanTable (need N2) ) )
 
(de rho (Tbl Inv I N)
(cond
((lt0 N) 0)
((=0 N) (get *X-Table I))
(T
(let Idx (+ I (>> 1 (* (- *TrigRows 1 N) (- *TrigRows N))))
(or
(get Inv Idx)
(set (nth Inv Idx) # only happens if value not computed yet
(+
(rho Tbl Inv (inc I) (- N 2))
(*/
(- (get Tbl I) (get Tbl (+ I N)))
1.0
(-
(rho Tbl Inv I (dec N))
(rho Tbl Inv (inc I) (dec N)) ) ) ) ) ) ) ) ) )
 
(de thiele (Tbl Inv X N)
(if (> N *TrigRows)
1.0
(+
(-
(rho Tbl Inv 1 (dec N))
(rho Tbl Inv 1 (- N 3)) )
(*/
(- X (get Tbl N))
1.0
(thiele Tbl Inv X (inc N)) ) ) ) )
 
(de iSin (X)
(thiele *SinTable *InvSinTable X 1) )
 
(de iCos (X)
(thiele *CosTable *InvCosTable X 1) )
 
(de iTan (X)
(thiele *TanTable *InvTanTable 1.0 1) )</syntaxhighlight>
Test:
<syntaxhighlight lang="picolisp">(prinl (round (* 6 (iSin 0.5)) 15))
(prinl (round (* 3 (iCos 0.5)) 15))
(prinl (round (* 4 (iTan 1.0)) 15))</syntaxhighlight>
Output:
<pre>3.141592653589793
3.141592653589793
3.141592653589793</pre>
 
=={{header|PowerShell}}==
<langsyntaxhighlight PowerShelllang="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
}
 
Line 372 ⟶ 1,357:
#uncomment to see the function
#"{$asin}"
6*$asin.invokeInvokeReturnAsIs(.5)[0]
$acos = (Thiele-Interpolation $acost)
#uncomment to see the function
#"{$acos}"
3*$acos.invokeInvokeReturnAsIs(.5)[0]
$atan = (Thiele-Interpolation $atant)
#uncomment to see the function
#"{$atan}"
4*$atan.invokeInvokeReturnAsIs(1)[0]</langsyntaxhighlight>
 
=={{header|Python}}==
{{trans|Go}}
<syntaxhighlight lang="python">#!/usr/bin/env python3
 
import math
 
def thieleInterpolator(x, y):
ρ = [[yi]*(len(y)-i) for i, yi in enumerate(y)]
for i in range(len(ρ)-1):
ρ[i][1] = (x[i] - x[i+1]) / (ρ[i][0] - ρ[i+1][0])
for i in range(2, len(ρ)):
for j in range(len(ρ)-i):
ρ[j][i] = (x[j]-x[j+i]) / (ρ[j][i-1]-ρ[j+1][i-1]) + ρ[j+1][i-2]
ρ0 = ρ[0]
def t(xin):
a = 0
for i in range(len(ρ0)-1, 1, -1):
a = (xin - x[i-1]) / (ρ0[i] - ρ0[i-2] + a)
return y[0] + (xin-x[0]) / (ρ0[1]+a)
return t
 
# task 1: build 32 row trig table
xVal = [i*.05 for i in range(32)]
tSin = [math.sin(x) for x in xVal]
tCos = [math.cos(x) for x in xVal]
tTan = [math.tan(x) for x in xVal]
# task 2: define inverses
iSin = thieleInterpolator(tSin, xVal)
iCos = thieleInterpolator(tCos, xVal)
iTan = thieleInterpolator(tTan, xVal)
# task 3: demonstrate identities
print('{:16.14f}'.format(6*iSin(.5)))
print('{:16.14f}'.format(3*iCos(.5)))
print('{:16.14f}'.format(4*iTan(1)))</syntaxhighlight>
{{out}}
<pre>
3.14159265358979
3.14159265358979
3.14159265358980
</pre>
 
=={{header|Racket}}==
<syntaxhighlight lang="racket">
#lang racket
(define xs (for/vector ([x (in-range 0.0 1.6 0.05)]) x))
(define (x i) (vector-ref xs i))
 
(define-syntax define-table
(syntax-rules ()
[(_ f tf rf if)
(begin (define tab (for/vector ([x xs]) (f x)))
(define (tf n) (vector-ref tab n))
(define cache (make-vector (/ (* 32 31) 2) #f))
(define (rf n thunk)
(or (vector-ref cache n)
(let ([v (thunk)])
(vector-set! cache n v)
v)))
(define (if t) (thiele tf x rf t 0)))]))
 
(define-table sin tsin rsin isin)
(define-table cos tcos rcos icos)
(define-table tan ttan rtan itan)
 
(define (rho x y r i n)
(cond
[(< n 0) 0]
[(= n 0) (y i)]
[else (r (+ (/ (* (- 32 1 n) (- 32 n)) 2) i)
(λ() (+ (/ (- (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)))))]))
(define (thiele x y r xin n)
(cond
[(> n 31) 1]
[(+ (rho x y r 0 n) (- (rho x y r 0 (- n 2)))
(/ (- xin (x n)) (thiele x y r xin (+ n 1))))]))
 
(* 6 (isin 0.5))
(* 3 (icos 0.5))
(* 4 (itan 1.))
</syntaxhighlight>
Output:
<syntaxhighlight lang="racket">
3.141592653589793
3.1415926535897936
3.1415926535897953
</syntaxhighlight>
 
=={{header|Raku}}==
(formerly Perl 6)
{{Works with|rakudo|2018.09}}<br>
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" line># reciprocal difference:
multi sub ρ(&f, @x where * < 1) { 0 } # Identity
multi sub ρ(&f, @x where * == 1) { &f(@x[0]) }
multi sub ρ(&f, @x where * > 1) {
( @x[0] - @x[* - 1] ) # ( x - x[n] )
/ (ρ(&f, @x[^(@x - 1)]) # / ( ρ[n-1](x[0], ..., x[n-1])
- ρ(&f, @x[1..^@x]) ) # - ρ[n-1](x[1], ..., x[n]) )
+ ρ(&f, @x[1..^(@x - 1)]); # + ρ[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
# must sort hash keys to maintain order between invocations
my $a = ρ(&f, %f.keys.sort[^($ord +1)]);
my $b = ρ(&f, %f.keys.sort[^($ord -1)]);
my $num = $x - %f.keys.sort[$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 = 12) {
 
my ($sin_pi, $cos_pi, $tan_pi);
my $p1 = Promise.start( { my %invsin = mk-inv(&sin, 0.05, $tblsz); $sin_pi = 6 * thiele(0.5, %invsin, 0) } );
my $p2 = Promise.start( { my %invcos = mk-inv(&cos, 0.05, $tblsz); $cos_pi = 3 * thiele(0.5, %invcos, 0) } );
my $p3 = Promise.start( { my %invtan = mk-inv(&tan, 0.05, $tblsz); $tan_pi = 4 * thiele(1.0, %invtan, 0) } );
await $p1, $p2, $p3;
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";
}</syntaxhighlight>
 
Output:
 
<pre>pi = 3.14159265358979
estimations using a table of 12 elements:
sin interpolation: 3.14159265358961
cos interpolation: 3.1387286696692
tan interpolation: 3.14159090545243</pre>
 
=={{header|Rust}}==
<syntaxhighlight lang="rust">
const N: usize = 32;
const STEP: f64 = 0.05;
 
fn main() {
let x: Vec<f64> = (0..N).map(|i| i as f64 * STEP).collect();
let sin = x.iter().map(|x| x.sin()).collect::<Vec<_>>();
let cos = x.iter().map(|x| x.cos()).collect::<Vec<_>>();
let tan = x.iter().map(|x| x.tan()).collect::<Vec<_>>();
 
println!(
"{}\n{}\n{}",
6. * thiele(&sin, &x, 0.5),
3. * thiele(&cos, &x, 0.5),
4. * thiele(&tan, &x, 1.)
);
}
 
fn thiele(x: &[f64], y: &[f64], xin: f64) -> f64 {
let mut p: Vec<Vec<f64>> = (0..N).map(|i| (i..N).map(|_| 0.0).collect()).collect();
 
(0..N).for_each(|i| p[i][0] = y[i]);
 
(0..N - 1).for_each(|i| p[i][1] = (x[i] - x[i + 1]) / (p[i][0] - p[i + 1][0]));
 
(2..N).for_each(|i| {
(0..N - i).for_each(|j| {
p[j][i] = (x[j] - x[j + i]) / (p[j][i - 1] - p[j + 1][i - 1]) + p[j + 1][i - 2];
})
});
 
let mut a = 0.;
(2..N).rev().for_each(|i| {
a = (xin - x[i - 1]) / (p[0][i] - p[0][i - 2] + a);
});
y[0] + (xin - x[0]) / (p[0][1] + a)
}
 
</syntaxhighlight>
{{out}}
<pre>
3.141592653589793
3.1415926535897936
3.1415926535897953
</pre>
 
=={{header|Sidef}}==
{{trans|Python}}
<syntaxhighlight lang="ruby">func thiele(x, y) {
var ρ = {|i| [y[i]]*(y.len-i) }.map(^y)
 
for i in ^(ρ.end) {
ρ[i][1] = ((x[i] - x[i+1]) / (ρ[i][0] - ρ[i+1][0]))
}
for i (2 .. ρ.end) {
for j (0 .. ρ.end-i) {
ρ[j][i] = (((x[j]-x[j+i]) / (ρ[j][i-1]-ρ[j+1][i-1])) + ρ[j+1][i-2])
}
}
 
var ρ0 = ρ[0]
 
func t(xin) {
var a = 0
for i (ρ0.len ^.. 2) {
a = ((xin - x[i-1]) / (ρ0[i] - ρ0[i-2] + a))
}
y[0] + ((xin-x[0]) / (ρ0[1]+a))
}
return t
}
 
# task 1: build 32 row trig table
var xVal = {|k| k * 0.05 }.map(^32)
var tSin = xVal.map { .sin }
var tCos = xVal.map { .cos }
var tTan = xVal.map { .tan }
 
# task 2: define inverses
var iSin = thiele(tSin, xVal)
var iCos = thiele(tCos, xVal)
var iTan = thiele(tTan, xVal)
 
# task 3: demonstrate identities
say 6*iSin(0.5)
say 3*iCos(0.5)
say 4*iTan(1)</syntaxhighlight>
{{out}}
<pre>
3.14159265358979323846438729976818601771260734312
3.14159265358979323846157620314930763214337987744
3.14159265358979323846264318595256260456200366896
</pre>
 
=={{header|Swift}}==
 
{{trans|Kotlin}}
 
<syntaxhighlight lang="swift">let N = 32
let N2 = N * (N - 1) / 2
let step = 0.05
 
var xval = [Double](repeating: 0, count: N)
var tsin = [Double](repeating: 0, count: N)
var tcos = [Double](repeating: 0, count: N)
var ttan = [Double](repeating: 0, count: N)
var rsin = [Double](repeating: .nan, count: N2)
var rcos = [Double](repeating: .nan, count: N2)
var rtan = [Double](repeating: .nan, count: N2)
 
func rho(_ x: [Double], _ y: [Double], _ r: inout [Double], _ i: Int, _ n: Int) -> Double {
guard n >= 0 else {
return 0
}
 
guard n != 0 else {
return y[i]
}
 
let idx = (N - 1 - n) * (N - n) / 2 + i
 
if r[idx] != r[idx] {
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]
}
 
func thiele(_ x: [Double], _ y: [Double], _ r: inout [Double], _ xin: Double, _ n: Int) -> Double {
guard n <= N - 1 else {
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)
}
 
for i in 0..<N {
xval[i] = Double(i) * step
tsin[i] = sin(xval[i])
tcos[i] = cos(xval[i])
ttan[i] = tsin[i] / tcos[i]
}
 
print(String(format: "%16.14f", 6 * thiele(tsin, xval, &rsin, 0.5, 0)))
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>
 
{{out}}
 
<pre>3.14159265358979
3.14159265358979
3.14159265358980</pre>
 
=={{header|Tcl}}==
{{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 428 ⟶ 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 445 ⟶ 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 452 ⟶ 1,744:
pi estimate using tan interpolation: 3.141592653589794
</pre>
 
=={{header|Wren}}==
{{trans|C}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "./fmt" for Fmt
 
var N = 32
var N2 = N * (N - 1) / 2
var STEP = 0.05
 
var xval = List.filled(N, 0.0)
var tsin = List.filled(N, 0.0)
var tcos = List.filled(N, 0.0)
var ttan = List.filled(N, 0.0)
var rsin = List.filled(N2, 0/0)
var rcos = List.filled(N2, 0/0)
var rtan = List.filled(N2, 0/0)
 
var rho
rho = Fn.new { |x, y, r, i, n|
if (n < 0) return 0
if (n == 0) return y[i]
var idx = (N - 1 - n) * (N - n) / 2 + i
if (r[idx].isNan) {
r[idx] = (x[i] - x[i + n]) /
(rho.call(x, y, r, i, n - 1) - rho.call(x, y, r, i + 1, n - 1)) +
rho.call(x, y, r, i + 1, n - 2)
}
return r[idx]
}
 
var thiele
thiele = Fn.new { |x, y, r, xin, n|
if (n > N - 1) return 1
return rho.call(x, y, r, 0, n) - rho.call(x, y, r, 0, n -2) +
(xin - x[n]) / thiele.call(x, y, r, xin, n + 1)
}
 
for (i in 0...N) {
xval[i] = i * STEP
tsin[i] = xval[i].sin
tcos[i] = xval[i].cos
ttan[i] = tsin[i] / tcos[i]
}
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))</syntaxhighlight>
 
{{out}}
<pre>
3.14159265358979
3.14159265358979
3.14159265358979
</pre>
 
=={{header|zkl}}==
{{trans|C}}
Please see the C example for the comments I've removed (this is an as pure-as-I-make-it translation).
<syntaxhighlight lang="zkl">const N=32, N2=(N * (N - 1) / 2), STEP=0.05;
 
fcn rho(xs,ys,rs, i,n){
if (n < 0) return(0.0);
if (not n) return(ys[i]);
idx := (N - 1 - n) * (N - n) / 2 + i;
if (Void==rs[idx])
rs[idx] = (xs[i] - xs[i + n])
/ (rho(xs, ys, rs, i, n - 1) - rho(xs, ys, rs, i + 1, n - 1))
+ rho(xs, ys, rs, i + 1, n - 2);
return(rs[idx]);
}
fcn thiele(xs,ys,rs, xin, n){
if (n > N - 1) return(1.0);
rho(xs, ys, rs, 0, n) - rho(xs, ys, rs, 0, n - 2)
+ (xin - xs[n]) / thiele(xs, ys, rs, xin, n + 1);
}
 
///////////
reg t_sin=L(), t_cos=L(), t_tan=L(),
r_sin=L(), r_cos=L(), r_tan=L(), xval=L();
i_sin := thiele.fpM("11101",t_sin, xval, r_sin, 0);
i_cos := thiele.fpM("11101",t_cos, xval, r_cos, 0);
i_tan := thiele.fpM("11101",t_tan, xval, r_tan, 0);
foreach i in (N){
xval.append(x:=STEP*i);
t_sin.append(x.sin());
t_cos.append(x.cos());
t_tan.append(t_sin[i] / t_cos[i]);
}
foreach i in (N2){ r_sin+Void; r_cos+Void; r_tan+Void; }
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)));</syntaxhighlight>
{{out}}
<pre>
3.14159265358979
3.14159265358979
3.14159265358979
</pre>
 
 
{{omit from|GUISS}}
Jjuanhdez
2,169

edits

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