function"/"(Left, Right : Real_Array)return Real_Array is Result : Real_Array (Left'Range);beginif Left'Length /= Right'Length thenraise Constraint_Error with"arrays not same size";endif;for I in Result'Rangeloop Result (I) := Left (I) / Right (I);endloop;return Result;end"/";

with the termination FN(x)=1{\displaystyle F_{N}(x)=1}, and the interpolation formula is now f(x)=F0(x){\displaystyle f(x)=F_{0}(x)}, easily implemented as a recursive function.

Note that each ρn{\displaystyle \rho _{n}} needs to look up ρn−1{\displaystyle \rho _{n-1}} twice, so the total look ups go up as O(2N){\displaystyle O(2^{N})} while there are only O(N2){\displaystyle O(N^{2})} values. This is a text book situation for memoization.

This example shows how the accuracy changes with the degree of interpolation. The table 'columns' are only constructed implicitly during the recursive calculation of rdiff and thiele, 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 x for accuracy.

Implemented to parallel the generalized formula, making for clearer, but slower, code. Offsetting that, the use of Promise allows concurrent calculations, so running all three types of interpolation should not take any longer than running just one (presuming available cores).

$asin=(Thiele-Interpolation $asint)#uncomment to see the function#"{$asin}"6*$asin.InvokeReturnAsIs(.5)$acos=(Thiele-Interpolation $acost)#uncomment to see the function#"{$acos}"3*$acos.InvokeReturnAsIs(.5)$atan=(Thiele-Interpolation $atant)#uncomment to see the function#"{$atan}"4*$atan.InvokeReturnAsIs(1)