path: root/src/interp/sfsfun.boot

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