

 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:

 



f
(
x
)
=
f
(

x

1


)
+








x
−

x

1












ρ

1


(

x

1


,

x

2


)
+








x
−

x

2












ρ

2


(

x

1


,

x

2


,

x

3


)
−
f
(

x

1


)
+








x
−

x

3












ρ

3


(

x

1


,

x

2


,

x

3


,

x

4


)
−

ρ

1


(

x

1


,

x

2


)
+
⋯

















{\displaystyle 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 }}}}}}}






ρ


{\displaystyle \rho }

   represents the   reciprocal difference,   demonstrated here for reference:

 




ρ

1


(

x

0


,

x

1


)
=




x

0


−

x

1




f
(

x

0


)
−
f
(

x

1


)





{\displaystyle \rho _{1}(x_{0},x_{1})={\frac {x_{0}-x_{1}}{f(x_{0})-f(x_{1})}}}


 




ρ

2


(

x

0


,

x

1


,

x

2


)
=




x

0


−

x

2





ρ

1


(

x

0


,

x

1


)
−

ρ

1


(

x

1


,

x

2


)



+
f
(

x

1


)


{\displaystyle \rho _{2}(x_{0},x_{1},x_{2})={\frac {x_{0}-x_{2}}{\rho _{1}(x_{0},x_{1})-\rho _{1}(x_{1},x_{2})}}+f(x_{1})}


 




ρ

n


(

x

0


,

x

1


,
…
,

x

n


)
=




x

0


−

x

n





ρ

n
−
1


(

x

0


,

x

1


,
…
,

x

n
−
1


)
−

ρ

n
−
1


(

x

1


,

x

2


,
…
,

x

n


)



+

ρ

n
−
2


(

x

1


,
…
,

x

n
−
1


)


{\displaystyle \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})}


Demonstrate Thiele's interpolation function by:

 Building a   32   row trig table of values   for   



x


{\displaystyle x}

   from   0   by   0.05   to   1.55   of the trig functions:
   sin
   cos
   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:
    6 × sin-1 ½ = 



π


{\displaystyle \pi }


    3 × cos-1 ½ = 



π


{\displaystyle \pi }


    4 × tan-1 1 = 



π


{\displaystyle \pi }

 


#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;
}