Initial population.

1 files changed, 1031 insertions, 0 deletions

diff --git a/src/interp/sfsfun.boot.pamphlet b/src/interp/sfsfun.boot.pamphlet
new file mode 100644
index 00000000..50bd7b5b
--- /dev/null
+++ b/src/interp/sfsfun.boot.pamphlet
@@ -0,0 +1,1031 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/interp sfsfun.boot}
+\author{The Axiom Team}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\begin{verbatim}
+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.
+
+\end{verbatim}
+\section{License}
+<<license>>=
+-- Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
+-- 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.
+
+@
+<<*>>=
+<<license>>
+
+-- Used to be SPECFNSF
+)package "BOOT" 
+
+FloatError(formatstring,arg) ==
+--        ERROR(formatstring,arg)
+        ERROR FORMAT([],formatstring,arg)
+
+nangenericcomplex () ==
+        1.0/COMPLEX(0.0)
+
+
+
+fracpart(x) ==
+        CADR(MULTIPLE_-VALUE_-LIST(FLOOR(x)))
+
+intpart(x) ==
+        CAR(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 := CAR(lx)+1
+                        restx := CADR(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 := CADR(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(DOUBLE_-FLOAT_-EPSILON/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)
+        SETF(AREF(a,0),0.0)
+        SETF(AREF(a,1),1.0)
+        for i in 2..n repeat
+                SETF(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)
+                        SETF(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)
+        SETF(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
+                SETF(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)
+        SETF(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
+                SETF(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)
+        SETF(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
+                SETF(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)
+        SETF(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
+                SETF(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
+        (FLOATP(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)
+        (FLOATP(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 ((FLOATP(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:=MAKE_-ARRAY(m)
+	SETF(AREF(w,m-1),BesselJAsymptOrder(v+m-1,z))
+        SETF(AREF(w,m-2),BesselJAsymptOrder(v+m-2,z))
+        for i in m-3 .. 0 by -1 repeat
+          SETF(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 FLOATP(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
+        FLOATP(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 FLOATP(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 := CAR(lm)    --- floor of real part of v
+        n := 2*MAX(20,m+10)  --- how large the back recurrence should be
+        tv := CADR(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)
+        SETF(AREF(w,m2), zero) --- start off
+        if type = '"I"
+        then
+                val := one
+        else
+                val := -one
+        m1 := n+2
+        SETF(AREF(w,m1), start)
+        m := n+1
+        xm := float(m)
+        ct1 := z2*(xm+v)
+        --- initialize
+        for m in (n+1)..1 by -1 repeat
+                SETF(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
+                SETF(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
+
+
+
+
+
+
+
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}