Thiele's interpolation formula: Difference between revisions

#* &nbsp; <big><big> 4 &times; tan^-1 1 = <math>\pi</math></big></big>

<br><br>

 

=={{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;

{{out}}

<pre>3.14159265358979

=={{header|C++}}==

{{trans|C}}

#include <iostream>

#include <iomanip>

 

return 0;

{{out}}

<pre>3.141592653589793115997963

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

thiele xs ys = f rho1 (tail xs)

where

[(sin, (2, 31)), (cos, (2, 100)), (tan, (4, 1000))]

-- N points taken uniformly from 0 to Pi/d

{{out}}

<pre>3.141592653589795

 

=={{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>

=={{header|Java}}==

{{trans|C}}

 

public class Test {

System.out.printf("%16.14f%n", 4 * thiele(t_tan, xval, r_tan, 1.0, 0));

}

<pre>3.14159265358979

3.14159265358979

Accuracy improves with a larger table and smaller step size.

{{trans|C}}

const N2 = N * div(N - 1, 2)

const step = 0.01

 

thiele_tables()

3.1415926535898335

3.141592653589818

=={{header|Kotlin}}==

{{trans|C}}

 

const val N = 32

println("%16.14f".format(3 * thiele(tcos, xval, rcos, 0.5, 0)))

println("%16.14f".format(4 * thiele(ttan, xval, rtan, 1.0, 0)))

 

{{out}}

3.14159265358980

</pre>

 

=={{header|Nim}}==

{{trans|Java}}

import math

 

echo fmt"{6 * thiele(tsin, xval, rsin, 0.5, 0):16.14f}"

echo fmt"{3 * thiele(tcos, xval, rcos, 0.5, 0):16.14f}"

{{out}}

<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:

=={{header|Perl}}==

{{trans|Sidef}}

use warnings;

use feature 'say';

say 6 * &$sin_inverse(0.5);

say 3 * &$cos_inverse(0.5);

{{out}}

<pre>3.14159265358979

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

{{out}}

(64 bit, obviously 3 fewer digits on 32 bit)

=={{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|Python}}==

{{trans|Go}}

 

import math

print('{:16.14f}'.format(6*iSin(.5)))

print('{:16.14f}'.format(3*iCos(.5)))

{{out}}

<pre>

 

=={{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|Raku}}==

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

 

multi sub ρ(&f, @x where * < 1) { 0 } # Identity

multi sub ρ(&f, @x where * == 1) { &f(@x[0]) }

say "cos interpolation: $cos_pi";

say "tan interpolation: $tan_pi";

 

Output:

cos interpolation: 3.1387286696692

tan interpolation: 3.14159090545243</pre>

 

=={{header|Sidef}}==

{{trans|Python}}

var ρ = {|i| [y[i]]*(y.len-i) }.map(^y)

 

say 6*iSin(0.5)

say 3*iCos(0.5)

{{out}}

<pre>

{{trans|Kotlin}}

 

let N2 = N * (N - 1) / 2

let step = 0.05

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

 

{{out}}

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

{{trans|C}}

{{libheader|Wren-fmt}}

 

var N = 32

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

 

{{out}}

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