NB. Horadam.ijs Script
NB. Exact Horadam Numbers with a Chebyshev Accent
NB. Written by Cliff Reiter, May 3, 1999
NB. (p,q,r,s) horadam n
NB. exp0 through exp5 list various experiments
NB. Updated to J6.02 December 2010

mp=:+/ . *
case=: 2&| + +:@(1&<)
w0=:(=i.2)"_
w1=:0 1&,:@x:@[
we=:[ mp~@:w -:@]
wo=:w1 mp [ w <:@]
w=: w0`w1`we`wo@.(case@])
horadam=:(2&{.@[ mp _2&{.@[ {.@w <:@])"1 0
chebyshev=: (, <:@+:@*: , _1: , +:)@]"0 horadam [

wrsc=: 1!:2&2
sxfmt=:(10&{. , '...('"_ , ":@(#-'_'&=@{.) , ')...'"_ , _10&{.)@":"0
sxmfmt=: ([,'    '"_ ,])/"2 @:sxfmt

horadamb=:4 : 0"1 0
'p q r s'=.x: x
k0=.((%r),0)"_
k1=.0 1x"_
ke=.((r,1)&mp@:*: , {: * ((+:r),s)&mp)@:k@-:
ko=.(0 1,:r,s)&mp@k@<:
k=.k0`k1`ke`ko@.case
((p*r),q) mp k <:y
)

exp0=:0 : 0                NB. warm up experiments
2 3 w0 ''
]W=.2 3 w1 ''
W mp W
(,:case)i.8
2 3 w 2
2 3 w 50
1 0 2 3 horadam 51
,.2 3 _3 3 horadam 249 250 251
sxfmt 1 0 2 3 horadam 51
,.2 3 _3 3 horadam 249 250 251
1 1 1 1 horadam >:i.10
1 1 1 1 horadam 100
)

exp1=:0 : 0                NB. timing experiemnts/results
   5000 (6!:2) 'F1=.1 1 1 1 horadam 10'
0.0011124
   1000 (6!:2) 'F1=:1 1 1 1 horadam 10'
0.001125
   500 (6!:2) 'F2=:1 1 1 1 horadam 100'
0.002374
   50 (6!:2) 'F3=:1 1 1 1 horadam 1000'
0.005
   10 (6!:2) 'F4=:1 1 1 1 horadam 10000'
0.0438
   (6!:2) 'F5=:1 1 1 1 horadam 100000'
3.86
   (6!:2) 'F6=:1 1 1 1 horadam 1000000'
398.172
   (6!:2) 'F7=:1 1 1 1 horadam 10000000'
62125.9
   sxfmt F4
3364476487...(2090)...9947366875
   sxfmt F5
2597406934...(20899)...3428746875
   sxfmt F6
1953282128...(208988)...8242546875
   sxfmt F7
1129834378...(2089877)...6380546875
)

exp2=:0 : 0   NB. special families
1 3 _1 2 horadam 1 2 3 4 5 6 7 8  NB. arthi starting 1 3
3 8 _1 2 horadam 1 2 3 4 5 6 7 8  NB. arthi starting 3 8
1 3 0 3 horadam >:i.10  NB. geom with ratio 3
1 1 1 1 horadam >:i.10  NB. fib
1 3 1 1 horadam >:i.10  NB. lucas
1 2 1 2 horadam >:i.10  NB. pell
1 3 1 2 horadam >:i.10  NB. pell-lucas
1 1 2 1 horadam >:i.10  NB. Jacobsthal
1 5 2 1 horadam >:i.10  NB. Jacobsthal-lucas
)

exp3=:0 : 0   NB. to plot Chebyshev polys
load 'plot'
$Y=:1 2 3 4 5 chebyshev"0 1 X=:_1++:(i.%<:)100
plot X;Y
)

exp4=: 0 : 0  NB. floating parameter example
x: _0.6 0.98 _2.7225 _3.3
float=: x:^:_1
float x: _0.6 0.98 _2.7225 _3.3
float (x: _0.6 0.98 _2.7225 _3.30) horadam 100 101 102
)

exp5=:0 : 0     NB. display recursive calls
w=: w0`w1`we`wo@.(case@wrsc@])
1 1 1 1 horadam 100
)