This implements
http://rosettacode.org/wiki/Thiele's_interpolation_formula
I've marked with NB. ???? spots I identify for possible improvement
(and for my personal edification).
span =: {. - {: NB. head - tail
spans =: span\ NB. apply span to successive infixes
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. NB. ???? Replace explicit loop with ^:
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. NB. ???? Continued fraction
a =. (x-i{xx)%-/(p{~i+2),(i{p),a
end.
(p{~>:i)+(x-i{xx)%(p{~i+2)+a
)
trig_table =: 1 2 3 o./ angles =: 5r100*i.32
(;:'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│
└─────────┴────────────┘
NB. ????
NB. The conjunction combines coef with eval.
NB. Wish to have J pre-compute the coefficients.
NB. Maybe M. is useful here?
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
)
's c t' =: trig_table
asin =: s thiele angles
6*asin 0.5
3.14159
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
