-- 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 sys_-macros namespace BOOT FloatError(formatstring,arg) == -- ERROR(formatstring,arg) ERROR formatToString(formatstring,arg) nangenericcomplex () == COMPLEX NaNQ() fracpart(x) == second integerAndFractionalParts x intpart(x) == first integerAndFractionalParts 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 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 := integerAndFractionalParts 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 integerAndFractionalParts 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 := newVector(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)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 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 := newVector(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 --- 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 --- 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)=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 := integerAndFractionalParts 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 := newVector(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)