Thiele's interpolation formula: Difference between revisions

{{task|Arithmetic operations}}

{{Wikipedia|Thiele's interpolation formula}}

 

 

 

 

 

 

Demonstrate Thiele's interpolation function by:

# Using columns from this table define an inverse - using Thiele's interpolation - for each trig function;

# Finally: demonstrate the following well known trigonometric identities:

 

=={{header|Ada}}==

thiele.ads:

 

generic

X, Y, RhoX : Real_Array (1 .. Length);

end record;

 

thiele.adb:

use type Real_Array;

 

end Inverse;

 

 

example:

with Ada.Numerics.Generic_Elementary_Functions;

with Thiele;

Long_Float'Image (4.0 * Inverse (Tan, 1.0)));

end;

 

output:

{{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.}}

 

# The MODE of lx and ly here should really be a UNION of "something REAL",

"tan", 4*inv tan(1)

))

Output:

<pre>

 

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.

#include <string.h>

#include <math.h>

printf("%16.14f\n", 4 * i_tan(1.));

return 0;

3.14159265358979

 

=={{header|Common Lisp}}==

Using the notations from above the C code instead of task desc.

(defparameter *thiele-length* 256)

(defparameter *rho-cache* (make-hash-table :test #'equal))

(format t "~f~%" (* 6 (inv-sin .5)))

(format t "~f~%" (* 3 (inv-cos .5)))

3.1415926535885172

 

=={{header|D}}==

 

struct Domain {

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);

{{out}}

<pre> 32 interpolating points

3.1415926535897932382 3 * inv_cos(0.5)

3.1415926535897932382 4 * inv_tan(1.0)</pre>

 

=={{header|Go}}==

{{trans|ALGOL 68}}

 

import (

return y[0] + (xin-x[0])/(ρ0[1]+a)

}

Output:

<pre>

=={{header|Haskell}}==

Caching of rho is automatic due to lazy lists.

 

 

{{out}}

 

=={{header|J}}==

span =: {. - {: NB. head - tail

spans =: span\ NB. apply span to successive infixes

 

<pre>

</pre>

 

NB. abscissae_of_knots coef ordinates_of_knots

NB. returns the interpolation coefficients for eval

(p{~>:i)+(x-i{xx)%(p{~i+2)+a

)

 

<pre>

</pre>

 

thiele =: 2 : 0

p =. _2 _{.,:n

(p{~>:i)+(y-i{m)%a+p{~i+2

)

 

<pre>

100*(_1&o. %~ _1&o. - asin) 1r5*i.6 NB. % error arcsin

0 1.4052 4.50319 9.32495 16.9438 39.321

</pre>

 

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.

 

 

let rec rdiff a l r =

Printf.printf "4*arctan(1.0) = %.15f\n" (4.0*.(interpolate 1.0 tan_tab n));;

 

Output:

<pre>Degree 8 interpolation:

4*arctan(1.0) = 3.141592653589793</pre>

 

 

}

 

<pre>

</pre>

 

=={{header|PicoLisp}}==

{{trans|C}}

(load "@lib/math.l")

 

 

(de iTan (X)

Test:

(prinl (round (* 3 (iCos 0.5)) 15))

Output:

<pre>3.141592653589793

 

=={{header|PowerShell}}==

{

$rho=@()

#uncomment to see the function

#"{$atan}"

 

=={{header|Racket}}==

#lang racket

(define xs (for/vector ([x (in-range 0.0 1.6 0.05)]) x))

(* 3 (icos 0.5))

(* 4 (itan 1.))

Output:

3.141592653589793

3.1415926535897936

3.1415926535897953

 

=={{header|Tcl}}==

{{works with|Tcl|8.5}}

{{trans|D}}

### Create a thiele-interpretation function with the given name that interpolates

### off the given table.

expr {$f1 + ($x - [lindex $X 1]) / ([lindex $rho 1] + $a)}

}} $X [lindex $p 1] [lindex $F 1]

Demonstration code:

for {set i 0} {$i < 32} {incr i} {

lappend trigTable(x) [set x [expr {0.05 * $i}]]

puts "pi estimate using sin interpolation: [expr {6 * [invSin 0.5]}]"

puts "pi estimate using cos interpolation: [expr {3 * [invCos 0.5]}]"

Output:

<pre>

pi estimate using cos interpolation: 3.141592653589793

pi estimate using tan interpolation: 3.141592653589794

</pre>

 

{{trans|C}}

Please see the C example for the comments I've removed (this is an as pure-as-I-make-it translation).

 

fcn rho(xs,ys,rs, i,n){

print("%16.14f\n".fmt( 6.0 * i_sin(0.5)));

print("%16.14f\n".fmt( 3.0 * i_cos(0.5)));

{{out}}

<pre>