\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}