aboutsummaryrefslogtreecommitdiff
path: root/src/interp/sfsfun.boot
blob: 1ccba154657e9470997a2f62ef2a95d006f96aca (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
-- Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
-- All rights reserved.
-- Copyright (C) 2007-2012, Gabriel Dos Reis.
-- All rights reserved.
--
-- Redistribution and use in source and binary forms, with or without
-- modification, are permitted provided that the following conditions are
-- met:
--
--     - Redistributions of source code must retain the above copyright
--       notice, this list of conditions and the following disclaimer.
--
--     - Redistributions in binary form must reproduce the above copyright
--       notice, this list of conditions and the following disclaimer in
--       the documentation and/or other materials provided with the
--       distribution.
--
--     - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--       names of its contributors may be used to endorse or promote products
--       derived from this software without specific prior written permission.
--
-- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
-- IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
-- TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
-- PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
-- OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
-- EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
-- PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
-- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
-- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
-- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


-- NOTEfrom TTT: at least BesselJAsymptOrder needs work

-- 1. This file contains the contents of BWC's original files:
--       floaterrors.boot
--       floatutils.boot
--       rgamma.boot
--       cgamma.boot
--       rpsi.boot
--       cpsi.boot
--       f01.boot
--       chebf01cmake.boot
--       chebevalsf.boot
--       besselIJ.boot

-- 2. All declarations have been commented out with "--@@"
--    since the boot translator is generating bad lisp code from them.

-- 3. The functions PsiAsymptotic, PsiEps and PsiAsymptoticOrder
--    had inconpatible definitions in rpsi.boot and cpsi.boot --
--    the local variables were declared float in one file and COMPLEX in
--    the other.  The type declarations have been commented out and the
--    duplicate definitions have been deleted.

-- 4. BesselIJ was not compiling.  I have modified the code from that
--    file to make it compile.  It should be checked for correctness.

-- SMW June 25, 1991

-- "Fixes" to BesselJ, B. Char June 14, 1992.  Needs extensive testing and
--   further fixes to BesselI and BesselJ.
-- More fixes to BesselJ, T. Tsikas 24 Feb, 1995.




import macros
namespace BOOT

FloatError(formatstring,arg) ==
--        ERROR(formatstring,arg)
        ERROR formatToString(formatstring,arg)

nangenericcomplex () ==
        COMPLEX NaNQ()



fracpart(x) ==
        second(MULTIPLE_-VALUE_-LIST(FLOOR(x)))

intpart(x) ==
        first(MULTIPLE_-VALUE_-LIST(FLOOR(x)))

negintp(x) ==
        if ZEROP IMAGPART(x) and x<0.0 and ZEROP fracpart(x)
        then
                true
        else
                false

--- Small float implementation of Gamma function.  Valid for
--- real arguments.  Maximum error (relative) seems to be
--- 2-4 ulps for real x 2<x<9, and up to ten ulps for larger x
--- up to overflow.  See Hart & Cheney.
--- Bruce Char, April, 1990.

horner(l,x) ==
        result := 0
        for el in l repeat
                result := result *x + el
        return result

rgamma (x) ==
        if COMPLEXP(x) then FloatError('"Gamma not implemented for complex value ~D",x)
        ZEROP (x-1.0) => 1.0
        if x>20 then gammaStirling(x) else gammaRatapprox(x)

lnrgamma (x) ==
        if x>20 then lnrgammaRatapprox(x) else LOG(gammaRatapprox(x))

cbeta(z,w) ==
        cgamma(z)*cgamma(w)/(cgamma(z+w))

gammaStirling(x) ==
       EXP(lnrgamma(x))

lnrgammaRatapprox(x) ==
       (x-.5)*LOG(x) - x + LOG(SQRT(2.0*PI)) + phiRatapprox(x)

phiRatapprox(x) ==
        arg := 1/(x**2)
        p := horner([.0666629070402007526,_
                     .6450730291289920251389,_
                     .670827838343321349614617,_
                     .12398282342474941538685913],arg);
        q := horner([1.0,7.996691123663644194772,_
                      8.09952718948975574728214,_
                      1.48779388109699298468156],arg);
        result := p/(x*q)
        result

gammaRatapprox (x) ==
        if (x>=2 and x<=3)
        then
                result := gammaRatkernel(x)
        else
                if x>3
                then
                     n := FLOOR(x)-2
                     a := x-n-2
                     reducedarg := 2+a
                     prod := */[reducedarg+i for i in 0..n-1]
                     result := prod* gammaRatapprox(reducedarg)
                else
                   if (2>x and x>0)
                   then
                     n := 2-FLOOR(x)
                     a := x-FLOOR(x)
                     reducedarg := 2+a
                     prod := */[x+i for i in 0..n-1]
                     result := gammaRatapprox(reducedarg)/prod
                 else
                        Pi := PI
                        lx := MULTIPLE_-VALUE_-LIST(FLOOR(x))
                        intpartx := first(lx)+1
                        restx := second(lx)
                        if ZEROP restx  -- case of negative non-integer value
                        then
                          FloatError ('"Gamma undefined for non-positive integers: ~D",x)
                          result := nangenericcomplex ()
                        else
                          result := Pi/(gammaRatapprox(1.0-x)*(-1.0)**(intpartx+1)*SIN(restx*Pi))
        result

gammaRatkernel(x) ==
           p := horner(reverse([3786.01050348257245475108,_
                        2077.45979389418732098416,_
                        893.58180452374981423868,_
                        222.1123961680117948396,_
                        48.95434622790993805232,_
                        6.12606745033608429879,_
                        .778079585613300575867]),x-2)
           q := horner(reverse([3786.01050348257187258861,_
                        476.79386050368791516095,_
                        -867.23098753110299445707,_
                        83.55005866791976957459,_
                        50.78847532889540973716,_
                        -13.40041478578134826274,_
                        1]),x-2.0)
           p/q

-- cgamma(z) Gamma function for complex arguments.
--    Bruce Char    April-May, 1990.
--
-- Our text for complex gamma is H. Kuki's paper Complex Gamma
-- Function with Error Control", CACM vol. 15, no. 4, ppp. 262-267.
-- (April, 1972.)  It uses the reflection formula and the basic
-- z+1 recurrence to transform the argument into something that
-- Stirling's asymptotic formula can handle.
--
-- However along the way it does a few tricky things to reduce
-- problems due to roundoff/cancellation error for particular values.

-- cgammat is auxiliary "t" function (see p. 263 Kuki)
cgammat(x) ==
        MAX(0.1, MIN(10.0, 10.0*SQRT(2.0) - abs(x)))

cgamma (z) ==
        z2 := IMAGPART(z)
        z1 := REALPART(z)       --- call real valued gamma if z is real
        if ZEROP z2
        then    result := rgamma(z1)
        else
                result := clngamma(z1,z2,z)
                result := EXP(result)
        result

lncgamma(z) ==
   clngamma(REALPART z, IMAGPART z, z)

clngamma(z1,z2,z) ==
        --- conjugate of gamma is gamma of conjugate.  map 2nd and 4th quads
        --- to first and third quadrants
        if z1<0.0
        then if z2 > 0.0
                then result := CONJUGATE(clngammacase1(z1,-z2,COMPLEX(z1,-z2)))
                else result := clngammacase1(z1,z2,z)
        else if z2 < 0.0
                then result := CONJUGATE(clngammacase23(z1,-z2,_
                                COMPLEX(z1,-z2)))
                else result := clngammacase23(z1,z2,z)
        result

clngammacase1(z1,z2,z) ==
        result1 := PiMinusLogSinPi(z1,z2,z)
        result2 := clngamma(1.0-z1,-z2,1.0-z)
        result1-result2

PiMinusLogSinPi(z1,z2,z) ==
        cgammaG(z1,z2)  - logH(z1,z2,z)

cgammaG(z1,z2) ==
        LOG(2*PI) + PI*z2 - COMPLEX(0.0,1.0)*PI*(z1-.5)

logH(z1,z2,z) ==
        z1bar := second(MULTIPLE_-VALUE_-LIST(FLOOR(z1))) ---frac part of z1
        piz1bar := PI*z1bar
        piz2 := PI*z2
        twopiz2 := 2.0*piz2
        i := COMPLEX(0.0,1.0)
        part2 := EXP(twopiz2)*(2.0*SIN(piz1bar)**2 + SIN(2.0*piz1bar)*i)
        part1 := -TANH(piz2)*(1.0+EXP(twopiz2))
--- part1 is another way of saying 1 - exp(2*Pi*z1bar)
        LOG(part1+part2)

clngammacase23(z1,z2,z) ==
        tz2 := cgammat(z2)
        if (z1 < tz2)
        then result:= clngammacase2(z1,z2,tz2,z)
        else result:= clngammacase3(z)
        result

clngammacase2(z1,z2,tz2,z) ==
        n := float(CEILING(tz2-z1))
        zpn := z+n
        (z-.5)*LOG(zpn) - (zpn) + cgammaBernsum(zpn) - cgammaAdjust(logS(z1,z2,z,n,zpn))

logS(z1,z2,z,n,zpn) ==
        sum := 0.0
        for k in 0..(n-1) repeat
                if z1+k < 5.0 - 0.6*z2
                then sum := sum + LOG((z+k)/zpn)
                else sum := sum + LOG(1.0 - (n-k)/zpn)
        sum

--- on p. 265, Kuki, logS result should have its imaginary part
--- adjusted by 2 Pi if it is negative.
cgammaAdjust(z) ==
        if IMAGPART(z)<0.0
        then result := z + COMPLEX(0.0, 2.0*PI)
        else result := z
        result

clngammacase3(z) ==
        (z- .5)*LOG(z) - z + cgammaBernsum(z)

cgammaBernsum (z) ==
        sum := LOG(2.0*PI)/2.0
        zterm := z
        zsquaredinv := 1.0/(z*z)
        l:= [.083333333333333333333, -.0027777777777777777778,_
                .00079365079365079365079,  -.00059523809523809523810,_
                .00084175084175084175084, -.0019175269175269175269,_
                .0064102564102564102564]
        for m in 1..7 for el in l repeat
                zterm := zterm*zsquaredinv
                sum := sum + el*zterm
        sum




--- nth derivatives of ln gamma for real x, n = 0,1,....
--- requires files floatutils, rgamma
$PsiAsymptoticBern := VECTOR(0.0, 0.1666666666666667, -0.03333333333333333, 0.02380952380952381,_
              -0.03333333333333333, 0.07575757575757576, -0.2531135531135531, 1.166666666666667,_
              -7.092156862745098, 54.97117794486216, -529.1242424242424, 6192.123188405797,_
              -86580.25311355311, 1425517.166666667, -27298231.06781609, 601580873.9006424,_
              -15116315767.09216, 429614643061.1667, -13711655205088.33, 488332318973593.2,_
              -19296579341940070.0,  841693047573682600.0, -40338071854059460000.0)


PsiAsymptotic(n,x) ==
        xn := x**n
        xnp1 := xn*x
        xsq := x*x
        xterm := xsq
        factterm := rgamma(n+2)/2.0/rgamma(float(n+1))
        --- initialize to 1/n!
        sum := AREF($PsiAsymptoticBern,1)*factterm/xterm
        for k in 2..22 repeat
                xterm := xterm * xsq
                if n=0 then factterm := 1.0/float(2*k)
                else if n=1 then factterm := 1
                else factterm := factterm * float(2*k+n-1)*float(2*k+n-2)/(float(2*k)*float(2*k-1))
                sum := sum + AREF($PsiAsymptoticBern,k)*factterm/xterm
        PsiEps(n,x) + 1.0/(2.0*xnp1) + 1.0/xn * sum


PsiEps(n,x) ==
        if n = 0
        then
                result := -LOG(x)
        else
                result :=  1.0/(float(n)*(x**n))
        result

PsiAsymptoticOrder(n,x,nterms) ==
        sum := 0
        xterm := 1.0
        np1 := n+1
        for k in 0..nterms repeat
                xterm := (x+float(k))**np1
                sum := sum + 1.0/xterm
        sum


rPsi(n,x) ==
        if x<=0.0
        then
                if ZEROP fracpart(x)
                then FloatError('"singularity encountered at ~D",x)
                else
                        m := MOD(n,2)
                        sign := (-1)**m
                        if fracpart(x)=.5
                        then
                                skipit := 1
                        else
                                skipit := 0
                        sign*((PI**(n+1))*cotdiffeval(n,PI*(-x),skipit) + rPsi(n,1.0-x))
        else if n=0
        then
                - rPsiW(n,x)
        else
                rgamma(float(n+1))*rPsiW(n,x)*(-1)**MOD(n+1,2)

---Amos' w function, with w(0,x) picked to be -psi(x) for x>0
rPsiW(n,x) ==
        if (x <=0 or n < 0)
        then
                HardError('"rPsiW not implemented for negative n or non-positive x")
        nd := 6         -- magic number for number of digits in a word?
        alpha := 3.5 + .40*nd
        beta := 0.21 + (.008677e-3)*(nd-3) + (.0006038e-4)*(nd-3)**2
        xmin := float(FLOOR(alpha + beta*n) + 1)
        if n>0
        then
                a := MIN(0,1.0/float(n)*LOG($DoubleFloatPrecision/MIN(1.0,x)))
                c := EXP(a)
                if abs(a) >= 0.001
                then
                        fln := x/c*(1.0-c)
                else
                        fln := -x*a/c
                bign := FLOOR(fln) + 1
--- Amos says to use alternative series for large order if ordinary
--- backwards recurrence too expensive
                if (bign < 15) and (xmin > 7.0+x)
                then
                        return PsiAsymptoticOrder(n,x,bign)
        if x>= xmin
        then
                return PsiAsymptotic(n,x)
---ordinary case -- use backwards recursion
        PsiBack(n,x,xmin)

PsiBack(n,x,xmin) ==
        xintpart := PsiIntpart(x)
        x0 := x-xintpart                ---frac part of x
        result := PsiAsymptotic(n,x0+xmin+1.0)
        for k in xmin..xintpart by -1 repeat
--- Why not decrement from x?   See Amos p. 498
                result := result + 1.0/((x0 + float(k))**(n+1))
        result


PsiIntpart(x) ==
        if x<0
        then
                result :=  -PsiInpart(-x)
        else
                result := FLOOR(x)
        return result


---Code for computation of derivatives of cot(z), necessary for
--- polygamma reflection formula.  If you want to compute n-th derivatives of
---cot(Pi*x), you have to multiply the result of cotdiffeval by Pi**n.

-- MCD: This is defined at the Lisp Level.
-- COT(z) ==
--         1.0/TAN(z)

cotdiffeval(n,z,skipit) ==
---skip=1 if arg z is known to be an exact multiple of Pi/2
        a := MAKE_-ARRAY(n+2)
        AREF(a,0) := 0.0
        AREF(a,1) := 1.0
        for i in 2..n repeat
          AREF(a,i) := 0.0
        for l in 1..n repeat
                m := MOD(l+1,2)
                for k in m..l+1 by 2 repeat
                        if k<1
                        then
                                t1 := 0
                        else
                                t1 := -AREF(a,k-1)*(k-1)
                        if k>l
                        then
                                t2 := 0
                        else
                                t2 := -AREF(a,k+1)*(k+1)
                        AREF(a,k) := t1+t2
        --- evaluate d^N/dX^N cot(z) via Horner-like rule
        v := COT(z)
        sq := v**2
        s := AREF(a,n+1)
        for i in (n-1)..0 by -2 repeat
                s := s*sq + AREF(a,i)
        m := MOD(n,2)
        if m=0
        then
                s := s*v
        if skipit=1
        then
                if m=0
                then
                        return 0
                else
                        return AREF(a,0)
        else
                return s
--- nth derivatives of ln gamma for complex z, n=0,1,...
--- requires files rpsi (and dependents), floaterrors
--- currently defined only in right half plane until reflection formula
--- works

--- B. Char, June, 1990.

cPsi(n,z) ==
        x := REALPART(z)
        y := IMAGPART(z)
        if ZEROP y 
        then    --- call real function if real
                return rPsi(n,x)
        if y<0.0
        then    -- if imagpart(z) negative, take conjugate of conjugate
                conjresult := cPsi(n,COMPLEX(x,-y))
                return COMPLEX(REALPART(conjresult),-IMAGPART(conjresult))
        nterms := 22
        bound := 10.0
        if x<0.0 --- and abs(z)>bound and abs(y)<bound
        then
                FloatError('"Psi implementation can't compute at ~S ",[n,z])
---             return cpsireflect(n,x,y,z)
        else if (x>0.0 and abs(z)>bound ) --- or (x<0.0 and abs(y)>bound)
        then
                return PsiXotic(n,PsiAsymptotic(n,z))
        else            --- use recursion formula
                m := CEILING(SQRT(bound*bound - y*y) - x + 1.0)
                result := COMPLEX(0.0,0.0)
                for k in 0..(m-1) repeat
                        result := result + 1.0/((z + float(k))**(n+1))
                return PsiXotic(n,result+PsiAsymptotic(n,z+m))

PsiXotic(n,result) ==
        rgamma(float(n+1))*(-1)**MOD(n+1,2)*result

cpsireflect(n,z) ==
        m := MOD(n,2)
        sign := (-1)**m
        sign*PI**(n+1)*cotdiffeval(n,PI*z,0) + cPsi(n,1.0-z)

--- c parameter to 0F1, possibly complex
--- z argument to 0F1
--- Depends on files floaterror, floatutils

--- Program transcribed from Fortran,, p. 80 Luke 1977

chebf01 (c,z) ==
--- w scale factor so that 0<z/w<1
--- n    n+2 coefficients will be produced stored in an array
---  indexed from 0 to n+1.
--- See Luke's books for further explanation
        n := 75 --- ad hoc decision
---     if abs(z)/abs(c) > 200.0 and abs(z)>10000.0
---     then
---             FloatError('"cheb0F1 not implemented for ~S < 1",[c,z])
        w := 2.0*z
--- arr will be used to store the Cheb. series coefficients
        four:= 4.0
        start := EXPT(10.0, -200)
        n1 := n+1
        n2 := n+2
        a3 := 0.0
        a2 := 0.0
        a1 := start     -- arbitrary starting value
        z1 := four/w
        ncount := n1
        arr := MAKE_-ARRAY(n2)
        AREF(arr,ncount) := start  -- start off
        x1 := n2
        c1 := 1.0 - c
        for ncount in n..0 by -1 repeat
                divfac := 1.0/x1
                x1 := x1 -1.0
                AREF(arr,ncount) := x1*((divfac+z1*(x1-c1))*a1 +_
                        (1.0/x1 + z1*(x1+c1+1.0))*a2-divfac*a3)
                a3 := a2
                a2 := a1
                a1 := AREF(arr,ncount)
        AREF(arr,0) := AREF(arr,0)/2.0
--  compute scale factor
        rho := AREF(arr,0)
        sum := rho
        p := 1.0
        for i in 1..n1 repeat
                rho := rho - p*AREF(arr,i)
                sum := sum+AREF(arr,i)
                p := -p
        for l in 0..n1 repeat
          AREF(arr,l) := AREF(arr,l)/rho
        sum := sum/rho
---     Now evaluate array at argument
        b := 0.0
        temp := 0.0
        for i in (n+1)..0 by -1 repeat
                cc := b
                b := temp
                temp := -cc + AREF(arr,i)
        temp


brutef01(c,z) ==
--  Use series definition.  Won't work well if cancellation occurs
        term := 1.0
        sum := term
        n := 0.0
        oldsum := 0.0
        maxnterms := 10000
        for i in 0..maxnterms until oldsum=sum repeat
                oldsum := sum
                term := term*z/(c+n)/(n+1.0)
                sum := sum + term
                n := n+1.0
        sum

f01(c,z)==
        if negintp(c)
        then
                FloatError('"0F1 not defined for negative integer parameter value ~S",c)
-- conditions when we'll permit the computation
        else if abs(c)<1000.0 and abs(z)<1000.0
        then
                brutef01(c,z)
        else if ZEROP IMAGPART(z) and ZEROP IMAGPART(c) and z>=0.0 and c>=0.0
        then
                brutef01(c,z)
---     else
---             t := SQRT(-z)
---             c1 := c-1.0
---             p := PHASE(c)
---             if abs(c)>10.0*abs(t) and p>=0.0 and PHASE(c)<.90*PI
---             then BesselJAsymptOrder(c1,2*t)*cgamma(c/(t**(c1)))
---             else if abs(t)>10.0*abs(c) and abs(t)>50.0
---             then BesselJAsympt(c1,2*t)*cgamma(c/(t**(c1)))
---             else
---                     FloatError('"0F1 not implemented for ~S",[c,z])
        else if (10.0*abs(c)>abs(z)) and abs(c)<1.0E4 and abs(z)<1.0E4
        then
                brutef01(c,z)
        else
                FloatError('"0F1 not implemented for ~S",[c,z])

--- c parameter to 0F1
--- w scale factor so that 0<z/w<1
--- n    n+2 coefficients will be produced stored in an array
---  indexed from 0 to n+1.
--- See Luke's books for further explanation

--- Program transcribed from Fortran,, p. 80 Luke 1977
chebf01coefmake (c,w,n) ==
--- arr will be used to store the Cheb. series coefficients
        four:= 4.0
        start := EXPT(10.0, -200)
        n1 := n+1
        n2 := n+2
        a3 := 0.0
        a2 := 0.0
        a1 := start     -- arbitrary starting value
        z1 := four/w
        ncount := n1
        arr := MAKE_-ARRAY(n2)
        AREF(arr,ncount) := start  -- start off
        x1 := n2
        c1 := 1.0 - c
        for ncount in n..0 by -1 repeat
                divfac := 1.0/x1
                x1 := x1 -1.0
                AREF(arr,ncount) := x1*((divfac+z1*(x1-c1))*a1 +_
                        (1.0/x1 + z1*(x1+c1+1.0))*a2-divfac*a3)
                a3 := a2
                a2 := a1
                a1 := AREF(arr,ncount)
        AREF(arr,0) := AREF(arr,0)/2.0
--  compute scale factor
        rho := AREF(arr,0)
        sum := rho
        p := 1.0
        for i in 1..n1 repeat
                rho := rho - p*AREF(arr,i)
                sum := sum+AREF(arr,i)
                p := -p
        for l in 0..n1 repeat
          AREF(arr,l) := AREF(arr,l)/rho
        sum := sum/rho
        return([sum,arr])




---evaluation of Chebychev series of degree n at x, where the series's
---coefficients are given by the list in descending order (coef. of highest
---power first)

---May be numerically unstable for certain lists of coefficients;
--- could possibly reverse sequence of coefficients

--- Cheney and Hart p. 15.

--- B. Char, March 1990

chebeval (coeflist,x) ==
        b := 0;
        temp := 0;
        y := 2*x;

        for el in coeflist repeat
                c := b;
                b := temp
                temp := y*b -c + el
        (temp -c)/2


chebevalarr(coefarr,x,n) ==
        b := 0;
        temp := 0;
        y := 2*x;

        for i in 1..n repeat
                c := b;
                b := temp
                temp := y*b -c + coefarr.i
        (temp -c)/2

--- If plist is a list of coefficients for the Chebychev approximation
--- of a function f(x), then chebderiveval computes the Chebychev approximation
--- of f'(x).  See Luke, "Special Functions and their approximations, vol. 1
--- Academic Press 1969., p. 329 (from Clenshaw and Cooper)

--- < definition to be supplied>

--- chebstareval(plist,n) computes a Chebychev approximation from a
--- coefficient list, using shifted Chebychev polynomials of the first kind
--- The defining relation is that T*(n,x) = T(n,2*x-1).  Thus the interval
--- [0,1] of T*n is the interval [-1,1] of Tn.

chebstareval(coeflist,x) ==          -- evaluation of T*(n,x)
        b := 0
        temp := 0
        y := 2*(2*x-1)

        for el in coeflist repeat
                c := b;
                b := temp
                temp := y*b -c + el
        temp - y*b/2


chebstarevalarr(coefarr,x,n) ==          -- evaluation of sum(C(n)*T*(n,x))

        b := 0
        temp := 0
        y := 2*(2*x-1)

        for i in (n+1)..0 by -1 repeat
                c := b
                b := temp
                temp := y*b -c + AREF(coefarr,i)
        temp - y*b/2

--Float definitions for Bessel functions I and J.
--External references:  cgamma, rgamma, chebf01coefmake, chebevalstarsf
-- floatutils

---BesselJ works for complex and real values of v and z
BesselJ(v,z) ==
---Ad hoc boundaries for approximation
        B1:= 10
        B2:= 10
        n := 50         --- number of terms in Chebychev series.
        --- tests for negative integer order
        (float?(v) and ZEROP fracpart(v) and (v<0)) or (COMPLEXP(v) and ZEROP IMAGPART(v) and ZEROP fracpart(REALPART(v)) and REALPART(v)<0.0) =>
             --- odd or even according to v (9.1.5 A&S)
             --- $J_{-n}(z)=(-1)^{n} J_{n}(z)$
             BesselJ(-v,z)*EXPT(-1.0,v)
        (float?(z) and  (z<0)) or (COMPLEXP(z) and REALPART(z)<0.0) =>
          --- negative argument (9.1.35 A&S) 
          --- $J_{\nu}(z e^{m \pi i}) = e^{m \nu \pi i} J_{\nu}(z)$
             BesselJ(v,-z)*EXPT(-1.0,v)
        ZEROP z and ((float?(v) and (v>=0.0)) or (COMPLEXP(v) and 
           ZEROP IMAGPART(v) and REALPART(v)>=0.0)) =>  --- zero arg, pos. real order
            ZEROP v => 1.0  --- J(0,0)=1
            0.0  --- J(v,0)=0 for real v>0
        rv := abs(v)
        rz := abs(z)
        (rz>B1) and (rz > B2*rv) =>  --- asymptotic argument
            BesselJAsympt(v,z)
        (rv>B1) and (rv > B2*rz) => --- asymptotic order
            BesselJAsymptOrder(v,z)
        (rz< B1) and (rv<B1) =>       --- small order and argument
                 arg := -(z*z)/4.0
                 w := 2.0*arg
                 vp1 := v+1.0
                 [sum,arr] := chebf01coefmake(vp1,w,n)
                 ---if we get NaNs then half n
                 while not _=(sum,sum) repeat
                        n:=FLOOR(n/2)
                        [sum,arr] := chebf01coefmake(vp1,w,n)
                 ---now n is safe, can we increase it (we know that 2*n is bad)?
                 chebstarevalarr(arr,arg/w,n)/cgamma(vp1)*EXPT(z/2.0,v)
        true => BesselJRecur(v,z)
        FloatError('"BesselJ not implemented for ~S", [v,z])

BesselJRecur(v,z) ==
        -- boost order
        --Numerical.Recipes. suggest so:=v+sqrt(n.s.f.^2*v)
        so:=15.0*z
        -- reduce order until non-zero
        while ZEROP abs(BesselJAsymptOrder(so,z)) repeat so:=so/2.0 
        if abs(so)<abs(z) then so:=v+18.*SQRT(v)
        m:= FLOOR(abs(so-v))+1
        w := newVector m
        AREF(w,m-1) := BesselJAsymptOrder(v+m-1,z)
        AREF(w,m-2) := BesselJAsymptOrder(v+m-2,z)
        for i in m-3 .. 0 by -1 repeat
          AREF(w,i) := 2.0 * (v+i+1.0) * AREF(w,i+1) /z -AREF(w,i+2)
        AREF(w,0)

BesselI(v,z) ==
        B1 := 15.0
        B2 := 10.0
        ZEROP(z) and float?(v) and (v>=0.0) =>  --- zero arg, pos. real order
            ZEROP(v) => 1.0  --- I(0,0)=1
            0.0             --- I(v,0)=0 for real v>0
--- Transformations for negative integer orders
        float?(v) and ZEROP(fracpart(v)) and (v<0) => BesselI(-v,z)
--- Halfplane transformations for Re(z)<0
        REALPART(z)<0.0 => BesselI(v,-z)*EXPT(-1.0,v)
--- Conjugation for complex order and real argument
        REALPART(v)<0.0 and not ZEROP IMAGPART(v) and float?(z) =>
              CONJUGATE(BesselI(CONJUGATE(v),z))
---We now know that Re(z)>= 0.0
        abs(z) > B1 =>    --- asymptotic argument case
                                FloatError('"BesselI not implemented for ~S",[v,z])
        abs(v) > B1 =>
                                FloatError('"BesselI not implemented for ~S",[v,z])
---     case of small argument and order
        REALPART(v)>= 0.0 =>  besselIback(v,z)
        REALPART(v)< 0.0 =>
                        chebterms := 50
                        besselIcheb(z,v,chebterms)
        FloatError('"BesselI not implemented for ~S",[v,z])

--- Compute n terms of the chebychev series for f01
besselIcheb(z,v,n) ==
        arg := (z*z)/4.0
        w := 2.0*arg;
        vp1 := v+1.0;
        [sum,arr] := chebf01coefmake(vp1,w,n)
        result := chebstarevalarr(arr,arg/w,n)/cgamma(vp1)*EXPT(z/2.0,v)

besselIback(v,z) ==
        ipv := IMAGPART(v)
        rpv := REALPART(v)
        lm := MULTIPLE_-VALUE_-LIST(FLOOR(rpv))
        m := first(lm)    --- floor of real part of v
        n := 2*MAX(20,m+10)  --- how large the back recurrence should be
        tv := second(lm)+(v-rpv) ---  fractional part of real part of v
                        --- plus imaginary part of v
        vp1 := tv+1.0;
        result := BesselIBackRecur(v,m,tv,z,'"I",n)
        result := result/cgamma(vp1)*EXPT(z/2.0,tv)

--- Backward recurrence for Bessel functions.  Luke (1975), p. 247.
--- works for -Pi< arg z <= Pi and  -Pi < arg v <= Pi
BesselIBackRecur(largev,argm,v,z,type,n) ==
--- v + m = largev
        one := 1.0
        two := 2.0
        zero := 0.0
        start := EXPT(10.0,-40)
        z2 := two/z
        m2 := n+3
        w:=MAKE_-ARRAY(m2+1)
        AREF(w,m2) := zero --- start off
        if type = '"I"
        then
                val := one
        else
                val := -one
        m1 := n+2
        AREF(w,m1) := start
        m := n+1
        xm := float(m)
        ct1 := z2*(xm+v)
        --- initialize
        for m in (n+1)..1 by -1 repeat
                AREF(w,m) := AREF(w,m+1)*ct1 + val*AREF(w,m+2)
                ct1 := ct1 - z2
        m := 1 + FLOOR(n/2)
        m2 := m + m -1
        if (v=0)
        then
                pn := AREF(w, m2 + 2)
                for m2 in (2*m-1)..3 by -2 repeat
                        pn := AREF(w, m2) - val *pn
                pn := AREF(w,1) - val*(pn+pn)
        else
                v1 := v-one
                xm := float(m)
                ct1 := v + xm + xm
                pn := ct1*AREF(w, m2 + 2)
                for m2 in (m+m -1)..3 by -2 repeat
                        ct1 := ct1 - two
                        pn := ct1*AREF(w,m2) - val*pn/xm*(v1+xm)
                        xm := xm - one
                pn := AREF(w,1) - val * pn
        m1 := n+2
        for m in 1..m1 repeat
                AREF(w,m) := AREF(w,m)/pn
        AREF(w,argm+1)




---Asymptotic functions for large values of z.  See p. 204 Luke 1969 vol. 1.

--- mu is 4*v**2
--- zsqr is z**2
--- zfth is z**4

BesselasymptA(mu,zsqr,zfth) ==
        (mu -1)/(16.0*zsqr) * (1 + (mu - 13.0)/(8.0*zsqr) + _
                (mu**2 - 53.0*mu + 412.0)/(48.0*zfth))

BesselasymptB(mu,z,zsqr,zfth) ==
        musqr := mu*mu
        z + (mu-1.0)/(8.0*z) *(1.0 + (mu - 25.0)/(48.0*zsqr) + _
                (musqr - 114.0*mu + 1073.0)/(640.0*zfth) +_
                (5.0*mu*musqr - 1535.0*musqr + 54703.0*mu - 375733.0)/(128.0*zsqr*zfth))

--- Asymptotic series only works when |v| < |z|.

BesselJAsympt (v,z) ==
        pi := PI
        mu := 4.0*v*v
        zsqr := z*z
        zfth := zsqr*zsqr
        SQRT(2.0/(pi*z))*EXP(BesselasymptA(mu,zsqr,zfth))*_
                COS(BesselasymptB(mu,z,zsqr,zfth) - pi*v/2.0 - pi/4.0)


---Asymptotic series for I.  See Whittaker, p. 373.
--- valid for -3/2 Pi < arg z < 1/2 Pi

BesselIAsympt(v,z,n) ==
        i := COMPLEX(0.0, 1.0)
        if (REALPART(z) = 0.0)
        then return EXPT(i,v)*BesselJ(v,-IMAGPART(z))
        sum1 := 0.0
        sum2 := 0.0
        fourvsq := 4.0*v**2
        two := 2.0
        eight := 8.0
        term1 := 1.0
---             sum1, sum2, fourvsq,two,i,eight,term1])
        for r in 1..n repeat
                term1 := -term1 *(fourvsq-(two*float(r)-1.0)**2)/_
                        (float(r)*eight*z)
                sum1 := sum1 + term1
                sum2 := sum2 + abs(term1)
        sqrttwopiz := SQRT(two*PI*z)
        EXP(z)/sqrttwopiz*(1.0 + sum1 ) +_
                EXP(-(float(n)+.5)*PI*i)*EXP(-z)/sqrttwopiz*(1.0+ sum2)


---Asymptotic formula for BesselJ when order is large comes from
---Debye (1909).  See Olver, Asymptotics and Special Functions, p. 134.
---Expansion good for 0<=phase(v)<Pi
---A&S recommend "uniform expansion" with complicated coefficients and Airy function.
---Debye's Formula is in 9.3.7,9.3.9,9.3.10 of A&S
---AXIOM recurrence for u_{k} 
---f(0)==1::EXPR INT
---f(n)== (t^2)*(1-t^2)*D(f(n-1),t)/2 + (1/8)*integrate( (1-5*t^2)*f(n-1),t)
BesselJAsymptOrder(v,z) ==
        sechalpha := z/v
        alpha := ACOSH(1.0/sechalpha)
        tanhalpha := SQRT(1.0-(sechalpha*sechalpha))
    --  cothalpha := 1.0/tanhalpha
        ca := 1.0/tanhalpha

        Pi := PI
        ca2:=ca*ca
        ca4:=ca2*ca2
        ca8:=ca4*ca4
        EXP(-v*(alpha-tanhalpha))/SQRT(2.0*Pi*v*tanhalpha)*_
        (1.0+_
        horner([              -5.0,                3.0],_
                                                                ca2)*ca/(v*24.0)+_
        horner([             385.0,             -462.0,              81.0],_
                                                                ca2)*ca2/(1152.0*v*v)+_
        horner([         -425425.0,           765765.0,         -369603.0,             30375.0],_
                                                                ca2)*ca2*ca/(414720.0*v*v*v)+_
        horner([       185910725.0,       -446185740.0,       349922430.0,        -94121676.0,         4465125.0],_
                                                                ca2)*ca4/(39813120.0*v*v*v*v)+_
        horner([   -188699385875.0,     566098157625.0,   -614135872350.0,     284499769554.0,    -49286948607.0,      1519035525.0],_
                                                                ca2)*ca4*ca/(6688604160.0*v*v*v*v*v)+_
        horner([1023694168371875.0,-3685299006138750.0,5104696716244125.0,-3369032068261860.0,1050760774457901.0,-127577298354750.0,2757049477875.0],_
                                                                ca2)*ca4*ca2/(4815794995200.0*v*v*v*v*v*v))
                

---  See Olver, p. 376-382.
BesselIAsymptOrder(v,vz) ==
        z := vz/v
        Pi := PI
---     Use reflection formula (Atlas, p. 492)  if v not in right half plane;  Is this always accurate?
        if REALPART(v)<0.0
        then return BesselIAsymptOrder(-v,vz) + 2.0/Pi*SIN(-v*Pi)*BesselKAsymptOrder(-v,vz)
---     Use the reflection formula (Atlas, p. 496) if z not in right half plane;
        if REALPART(vz) < 0.0
        then return EXPT(-1.0,v)*BesselIAsymptOrder(v,-vz)
        vinv := 1.0/v
        opzsqroh := SQRT(1.0+z*z)
        eta := opzsqroh + LOG(z/(1.0+opzsqroh))
        p := 1.0/opzsqroh
        p2 := p*p
        p4 := p2*p2
        u0p := 1.
        u1p := 1.0/8.0*p-5.0/24.0*p*p2
        u2p := (9.0/128.0+(-77.0/192.0+385.0/1152.0*p2)*p2)*p2
        u3p := (75.0/1024.0+(-4563.0/5120.0+(17017.0/9216.0-85085.0/82944.0*p2)_
                *p2)*p2)*p2*p
        u4p := (3675.0/32768.0+(-96833.0/40960.0+(144001.0/16384.0+(-7436429.0/663552.0+37182145.0/7962624.0*p2)*p2)*p2)*p2)*p4
        u5p := (59535.0/262144.0+(-67608983.0/9175040.0+(250881631.0/5898240.0+(-108313205.0/1179648.0+(5391411025.0/63700992.0-5391411025.0/191102976.0*p2)*p2)*p2)*p2)*p2)*p4*p
        hornerresult := horner([u5p,u4p,u3p,u2p,u1p,u0p],vinv)
        EXP(v*eta)/(SQRT(2.0*Pi*v)*SQRT(opzsqroh))*hornerresult


---See also Olver, pp. 376-382
BesselKAsymptOrder (v,vz) ==
        z := vz/v
        vinv := 1.0/v
        opzsqroh := SQRT(1.0+z*z)
        eta := opzsqroh + LOG(z/(1.0+opzsqroh))
        p := 1.0/opzsqroh
        p2 := p**2
        p4 := p2**2
        u0p := 1.
        u1p := (1.0/8.0*p-5.0/24.0*p**3)*(-1.0)
        u2p := (9.0/128.0+(-77.0/192.0+385.0/1152.0*p2)*p2)*p2
        u3p := ((75.0/1024.0+(-4563.0/5120.0+(17017.0/9216.0-85085.0/82944.0*p2)_
                *p2)*p2)*p2*p)*(-1.0)
        u4p := (3675.0/32768.0+(-96833.0/40960.0+(144001.0/16384.0+(-7436429.0/663552.0+37182145.0/7962624.0*p2)*p2)*p2)*p2)*p4
        u5p := ((59535.0/262144.0+(-67608983.0/9175040.0+(250881631.0/5898240.0+(-108313205.0/1179648.0+(5391411025.0/63700992.0-5391411025.0/191102976.0*p2)*p2)*p2)*p2)*p2)*p4*p)*(-1.0)
        hornerresult := horner([u5p,u4p,u3p,u2p,u1p,u0p],vinv)
        SQRT(PI/(2.0*v))*EXP(-v*eta)/(SQRT(opzsqroh))*hornerresult