diff --git a/changelog b/changelog index de5301a..2a03627 100644 --- a/changelog +++ b/changelog @@ -1,3 +1,7 @@ +20090828 tpd src/axiom-website/patches.html 20090828.04.tpd.patch +20090828 tpd src/interp/Makefile move sfsfun.boot to sfsfun.lisp +20090828 tpd src/interp/sfsfun.lisp added, rewritten from sfsfun.boot +20090828 tpd src/interp/sfsfun.boot removed, rewritten to sfsfun.lisp 20090828 tpd src/axiom-website/patches.html 20090828.03.tpd.patch 20090828 tpd src/interp/Makefile move compiler.boot to compiler.lisp 20090828 tpd src/interp/compiler.lisp added, rewritten from compiler.boot diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html index 280f14b..e367381 100644 --- a/src/axiom-website/patches.html +++ b/src/axiom-website/patches.html @@ -1934,5 +1934,7 @@ package.lisp rewrite from boot to lisp
htcheck.lisp rewrite from boot to lisp
20090828.03.tpd.patch compiler.lisp rewrite from boot to lisp
+20090828.04.tpd.patch +sfsfun.lisp rewrite from boot to lisp
diff --git a/src/interp/Makefile.pamphlet b/src/interp/Makefile.pamphlet index 6b57e43..ce7e6bc 100644 --- a/src/interp/Makefile.pamphlet +++ b/src/interp/Makefile.pamphlet @@ -4669,44 +4669,26 @@ ${MID}/pf2sex.lisp: ${IN}/pf2sex.lisp.pamphlet @ -\subsection{sfsfun.boot} +\subsection{sfsfun.lisp} <>= -${OUT}/sfsfun.${O}: ${MID}/sfsfun.clisp - @ echo 576 making ${OUT}/sfsfun.${O} from ${MID}/sfsfun.clisp - @ if [ -z "${NOISE}" ] ; then \ - echo '(progn (compile-file "${MID}/sfsfun.clisp"' \ +${OUT}/sfsfun.${O}: ${MID}/sfsfun.lisp + @ echo 136 making ${OUT}/sfsfun.${O} from ${MID}/sfsfun.lisp + @ ( cd ${MID} ; \ + if [ -z "${NOISE}" ] ; then \ + echo '(progn (compile-file "${MID}/sfsfun.lisp"' \ ':output-file "${OUT}/sfsfun.${O}") (${BYE}))' | ${DEPSYS} ; \ else \ - echo '(progn (compile-file "${MID}/sfsfun.clisp"' \ + echo '(progn (compile-file "${MID}/sfsfun.lisp"' \ ':output-file "${OUT}/sfsfun.${O}") (${BYE}))' | ${DEPSYS} \ >${TMP}/trace ; \ - fi + fi ) @ -<>= -${MID}/sfsfun.clisp: ${IN}/sfsfun.boot.pamphlet - @ echo 577 making ${MID}/sfsfun.clisp from ${IN}/sfsfun.boot.pamphlet +<>= +${MID}/sfsfun.lisp: ${IN}/sfsfun.lisp.pamphlet + @ echo 137 making ${MID}/sfsfun.lisp from ${IN}/sfsfun.lisp.pamphlet @ (cd ${MID} ; \ - ${TANGLE} ${IN}/sfsfun.boot.pamphlet >sfsfun.boot ; \ - if [ -z "${NOISE}" ] ; then \ - echo '(progn (boottran::boottocl "${MID}/sfsfun.boot") (${BYE}))' \ - | ${BOOTSYS} ; \ - else \ - echo '(progn (boottran::boottocl "${MID}/sfsfun.boot") (${BYE}))' \ - | ${BOOTSYS} >${TMP}/trace ; \ - fi ; \ - rm sfsfun.boot ) - -@ -<>= -${DOC}/sfsfun.boot.dvi: ${IN}/sfsfun.boot.pamphlet - @echo 578 making ${DOC}/sfsfun.boot.dvi from ${IN}/sfsfun.boot.pamphlet - @(cd ${DOC} ; \ - cp ${IN}/sfsfun.boot.pamphlet ${DOC} ; \ - ${DOCUMENT} ${NOISE} sfsfun.boot ; \ - rm -f ${DOC}/sfsfun.boot.pamphlet ; \ - rm -f ${DOC}/sfsfun.boot.tex ; \ - rm -f ${DOC}/sfsfun.boot ) + ${TANGLE} ${IN}/sfsfun.lisp.pamphlet >sfsfun.lisp ) @ @@ -5609,8 +5591,7 @@ clean: <> <> -<> -<> +<> <> <> diff --git a/src/interp/sfsfun.boot.pamphlet b/src/interp/sfsfun.boot.pamphlet deleted file mode 100644 index 6536e2e..0000000 --- a/src/interp/sfsfun.boot.pamphlet +++ /dev/null @@ -1,1036 +0,0 @@ -\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} -<>= --- 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. - -@ -<<*>>= -<> - --- 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 - --- Lisp PI is a long float and causes type errors, here we give --- enough digits to have double accuracy even after conversion --- to binary -DEFCONSTANT(dfPi,3.14159265358979323846264338328) - ---- 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*dfPi)) + 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 := dfPi - 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*dfPi) + dfPi*z2 - COMPLEX(0.0,1.0)*dfPi*(z1-.5) - -logH(z1,z2,z) == - z1bar := CADR(MULTIPLE_-VALUE_-LIST(FLOOR(z1))) ---frac part of z1 - piz1bar := dfPi*z1bar - piz2 := dfPi*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*dfPi) - else result := z - result - -clngammacase3(z) == - (z- .5)*LOG(z) - z + cgammaBernsum(z) - -cgammaBernsum (z) == - sum := LOG(2.0*dfPi)/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*((dfPi**(n+1))*cotdiffeval(n,dfPi*(-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)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*dfPi**(n+1)*cotdiffeval(n,dfPi*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 := 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*dfPi ---- 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 - (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 --- 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 - 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 := dfPi - 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*dfPi*z) - EXP(z)/sqrttwopiz*(1.0 + sum1 ) +_ - EXP(-(float(n)+.5)*dfPi*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)>= + +(IN-PACKAGE "BOOT") + +;-- Used to be SPECFNSF + +;FloatError(formatstring,arg) == +;-- ERROR(formatstring,arg) +; ERROR FORMAT([],formatstring,arg) + +(DEFUN |FloatError| (|formatstring| |arg|) + (PROG () (RETURN (ERROR (FORMAT NIL |formatstring| |arg|))))) + +;nangenericcomplex () == +; 1.0/COMPLEX(0.0) + +(DEFUN |nangenericcomplex| () + (PROG () (RETURN (/ 1.0 (COMPLEX 0.0))))) + +;fracpart(x) == +; CADR(MULTIPLE_-VALUE_-LIST(FLOOR(x))) + +(DEFUN |fracpart| (|x|) + (PROG () (RETURN (CADR (MULTIPLE-VALUE-LIST (FLOOR |x|)))))) + +;intpart(x) == +; CAR(MULTIPLE_-VALUE_-LIST(FLOOR(x))) + +(DEFUN |intpart| (|x|) + (PROG () (RETURN (CAR (MULTIPLE-VALUE-LIST (FLOOR |x|)))))) + +;negintp(x) == +; if ZEROP IMAGPART(x) and x<0.0 and ZEROP fracpart(x) +; then +; true +; else +; false + +(DEFUN |negintp| (|x|) + (PROG () + (RETURN + (COND + ((AND (ZEROP (IMAGPART |x|)) (< |x| 0.0) + (ZEROP (|fracpart| |x|))) + T) + ('T NIL))))) + +;-- Lisp PI is a long float and causes type errors, here we give +;-- enough digits to have double accuracy even after conversion +;-- to binary +;DEFCONSTANT(dfPi,3.14159265358979323846264338328) + +(EVAL-WHEN (EVAL LOAD) + (PROG () (RETURN (DEFCONSTANT |dfPi| 3.1415926535897936)))) + +;--- 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) + +(DEFUN |rgamma| (|x|) + (PROG () + (RETURN + (PROGN + (COND + ((COMPLEXP |x|) + (|FloatError| "Gamma not implemented for complex value ~D" + |x|))) + (COND + ((ZEROP (- |x| 1.0)) 1.0) + ('T + (COND + ((< 20 |x|) (|gammaStirling| |x|)) + ('T (|gammaRatapprox| |x|))))))))) + +;lnrgamma (x) == +; if x>20 then lnrgammaRatapprox(x) else LOG(gammaRatapprox(x)) + +(DEFUN |lnrgamma| (|x|) + (PROG () + (RETURN + (COND + ((< 20 |x|) (|lnrgammaRatapprox| |x|)) + ('T (LOG (|gammaRatapprox| |x|))))))) + +;cbeta(z,w) == +; cgamma(z)*cgamma(w)/(cgamma(z+w)) + +(DEFUN |cbeta| (|z| |w|) + (PROG () + (RETURN + (/ (* (|cgamma| |z|) (|cgamma| |w|)) (|cgamma| (+ |z| |w|)))))) + +;gammaStirling(x) == +; EXP(lnrgamma(x)) + +(DEFUN |gammaStirling| (|x|) + (PROG () (RETURN (EXP (|lnrgamma| |x|))))) + +;lnrgammaRatapprox(x) == +; (x-.5)*LOG(x) - x + LOG(SQRT(2.0*dfPi)) + phiRatapprox(x) + +(DEFUN |lnrgammaRatapprox| (|x|) + (PROG () + (RETURN + (+ (+ (- (* (- |x| 0.5) (LOG |x|)) |x|) + (LOG (SQRT (* 2.0 |dfPi|)))) + (|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 + +(DEFUN |phiRatapprox| (|x|) + (PROG (|result| |q| |p| |arg|) + (RETURN + (PROGN + (SETQ |arg| (/ 1 (EXPT |x| 2))) + (SETQ |p| + (|horner| + (LIST 0.066662907040200753 0.64507302912899211 + 0.67082783834332138 0.12398282342474939) + |arg|)) + (SETQ |q| + (|horner| + (LIST 1.0 7.9966911236636431 8.0995271894897574 + 1.4877938810969931) + |arg|)) + (SETQ |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 := dfPi +; 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 + +(DEFUN |gammaRatapprox| (|x|) + (PROG (|restx| |intpartx| |lx| |Pi| |prod| |reducedarg| |a| |n| + |result|) + (RETURN + (PROGN + (COND + ((AND (NOT (< |x| 2)) (NOT (< 3 |x|))) + (SETQ |result| (|gammaRatkernel| |x|))) + ((< 3 |x|) (SETQ |n| (- (FLOOR |x|) 2)) + (SETQ |a| (- (- |x| |n|) 2)) (SETQ |reducedarg| (+ 2 |a|)) + (SETQ |prod| + ((LAMBDA (|bfVar#3| |bfVar#2| |i|) + (LOOP + (COND + ((> |i| |bfVar#2|) (RETURN |bfVar#3|)) + ('T + (SETQ |bfVar#3| + (* |bfVar#3| (+ |reducedarg| |i|))))) + (SETQ |i| (+ |i| 1)))) + 1 (- |n| 1) 0)) + (SETQ |result| (* |prod| (|gammaRatapprox| |reducedarg|)))) + ((AND (< |x| 2) (< 0 |x|)) (SETQ |n| (- 2 (FLOOR |x|))) + (SETQ |a| (- |x| (FLOOR |x|))) (SETQ |reducedarg| (+ 2 |a|)) + (SETQ |prod| + ((LAMBDA (|bfVar#5| |bfVar#4| |i|) + (LOOP + (COND + ((> |i| |bfVar#4|) (RETURN |bfVar#5|)) + ('T (SETQ |bfVar#5| (* |bfVar#5| (+ |x| |i|))))) + (SETQ |i| (+ |i| 1)))) + 1 (- |n| 1) 0)) + (SETQ |result| (/ (|gammaRatapprox| |reducedarg|) |prod|))) + ('T (SETQ |Pi| |dfPi|) + (SETQ |lx| (MULTIPLE-VALUE-LIST (FLOOR |x|))) + (SETQ |intpartx| (+ (CAR |lx|) 1)) + (SETQ |restx| (CADR |lx|)) + (COND + ((ZEROP |restx|) + (|FloatError| + "Gamma undefined for non-positive integers: ~D" |x|) + (SETQ |result| (|nangenericcomplex|))) + ('T + (SETQ |result| + (/ |Pi| + (* (* (|gammaRatapprox| (- 1.0 |x|)) + (EXPT (- 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 + +(DEFUN |gammaRatkernel| (|x|) + (PROG (|q| |p|) + (RETURN + (PROGN + (SETQ |p| + (|horner| + (REVERSE (LIST 3786.0105034825724 2077.4597938941874 + 893.58180452374972 222.11239616801177 + 48.95434622790993 6.1260674503360839 + 0.77807958561330059)) + (- |x| 2))) + (SETQ |q| + (|horner| + (REVERSE (LIST 3786.0105034825719 476.79386050368794 + (- 867.23098753110298) + 83.550058667919771 50.788475328895402 + (- 13.400414785781347) 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))) + +(DEFUN |cgammat| (|x|) + (PROG () + (RETURN + (MAX 0.10000000000000001 + (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 + +(DEFUN |cgamma| (|z|) + (PROG (|result| |z1| |z2|) + (RETURN + (PROGN + (SETQ |z2| (IMAGPART |z|)) + (SETQ |z1| (REALPART |z|)) + (COND + ((ZEROP |z2|) (SETQ |result| (|rgamma| |z1|))) + ('T (SETQ |result| (|clngamma| |z1| |z2| |z|)) + (SETQ |result| (EXP |result|)))) + |result|)))) + +;lncgamma(z) == +; clngamma(REALPART z, IMAGPART z, z) + +(DEFUN |lncgamma| (|z|) + (PROG () (RETURN (|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 + +(DEFUN |clngamma| (|z1| |z2| |z|) + (PROG (|result|) + (RETURN + (PROGN + (COND + ((< |z1| 0.0) + (COND + ((< 0.0 |z2|) + (SETQ |result| + (CONJUGATE + (|clngammacase1| |z1| (- |z2|) + (COMPLEX |z1| (- |z2|)))))) + ('T (SETQ |result| (|clngammacase1| |z1| |z2| |z|))))) + ((< |z2| 0.0) + (SETQ |result| + (CONJUGATE + (|clngammacase23| |z1| (- |z2|) + (COMPLEX |z1| (- |z2|)))))) + ('T (SETQ |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 + +(DEFUN |clngammacase1| (|z1| |z2| |z|) + (PROG (|result2| |result1|) + (RETURN + (PROGN + (SETQ |result1| (|PiMinusLogSinPi| |z1| |z2| |z|)) + (SETQ |result2| (|clngamma| (- 1.0 |z1|) (- |z2|) (- 1.0 |z|))) + (- |result1| |result2|))))) + +;PiMinusLogSinPi(z1,z2,z) == +; cgammaG(z1,z2) - logH(z1,z2,z) + +(DEFUN |PiMinusLogSinPi| (|z1| |z2| |z|) + (PROG () (RETURN (- (|cgammaG| |z1| |z2|) (|logH| |z1| |z2| |z|))))) + +;cgammaG(z1,z2) == +; LOG(2*dfPi) + dfPi*z2 - COMPLEX(0.0,1.0)*dfPi*(z1-.5) + +(DEFUN |cgammaG| (|z1| |z2|) + (PROG () + (RETURN + (- (+ (LOG (* 2 |dfPi|)) (* |dfPi| |z2|)) + (* (* (COMPLEX 0.0 1.0) |dfPi|) (- |z1| 0.5)))))) + +;logH(z1,z2,z) == +; z1bar := CADR(MULTIPLE_-VALUE_-LIST(FLOOR(z1))) ---frac part of z1 +; piz1bar := dfPi*z1bar +; piz2 := dfPi*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) + +(DEFUN |logH| (|z1| |z2| |z|) + (PROG (|part1| |part2| |i| |twopiz2| |piz2| |piz1bar| |z1bar|) + (RETURN + (PROGN + (SETQ |z1bar| (CADR (MULTIPLE-VALUE-LIST (FLOOR |z1|)))) + (SETQ |piz1bar| (* |dfPi| |z1bar|)) + (SETQ |piz2| (* |dfPi| |z2|)) + (SETQ |twopiz2| (* 2.0 |piz2|)) + (SETQ |i| (COMPLEX 0.0 1.0)) + (SETQ |part2| + (* (EXP |twopiz2|) + (+ (* 2.0 (EXPT (SIN |piz1bar|) 2)) + (* (SIN (* 2.0 |piz1bar|)) |i|)))) + (SETQ |part1| (- (* (TANH |piz2|) (+ 1.0 (EXP |twopiz2|))))) + (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 + +(DEFUN |clngammacase23| (|z1| |z2| |z|) + (PROG (|result| |tz2|) + (RETURN + (PROGN + (SETQ |tz2| (|cgammat| |z2|)) + (COND + ((< |z1| |tz2|) + (SETQ |result| (|clngammacase2| |z1| |z2| |tz2| |z|))) + ('T (SETQ |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)) + +(DEFUN |clngammacase2| (|z1| |z2| |tz2| |z|) + (PROG (|zpn| |n|) + (RETURN + (PROGN + (SETQ |n| (|float| (CEILING (- |tz2| |z1|)))) + (SETQ |zpn| (+ |z| |n|)) + (- (+ (- (* (- |z| 0.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 + +(DEFUN |logS| (|z1| |z2| |z| |n| |zpn|) + (PROG (|sum|) + (RETURN + (PROGN + (SETQ |sum| 0.0) + ((LAMBDA (|bfVar#6| |k|) + (LOOP + (COND + ((> |k| |bfVar#6|) (RETURN NIL)) + ('T + (COND + ((< (+ |z1| |k|) + (- 5.0 (* 0.60000000000000009 |z2|))) + (SETQ |sum| (+ |sum| (LOG (/ (+ |z| |k|) |zpn|))))) + ('T + (SETQ |sum| + (+ |sum| (LOG (- 1.0 (/ (- |n| |k|) |zpn|))))))))) + (SETQ |k| (+ |k| 1)))) + (- |n| 1) 0) + |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*dfPi) +; else result := z +; result + +(DEFUN |cgammaAdjust| (|z|) + (PROG (|result|) + (RETURN + (PROGN + (COND + ((< (IMAGPART |z|) 0.0) + (SETQ |result| (+ |z| (COMPLEX 0.0 (* 2.0 |dfPi|))))) + ('T (SETQ |result| |z|))) + |result|)))) + +;clngammacase3(z) == +; (z- .5)*LOG(z) - z + cgammaBernsum(z) + +(DEFUN |clngammacase3| (|z|) + (PROG () + (RETURN + (+ (- (* (- |z| 0.5) (LOG |z|)) |z|) (|cgammaBernsum| |z|))))) + +;cgammaBernsum (z) == +; sum := LOG(2.0*dfPi)/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 + +(DEFUN |cgammaBernsum| (|z|) + (PROG (|l| |zsquaredinv| |zterm| |sum|) + (RETURN + (PROGN + (SETQ |sum| (/ (LOG (* 2.0 |dfPi|)) 2.0)) + (SETQ |zterm| |z|) + (SETQ |zsquaredinv| (/ 1.0 (* |z| |z|))) + (SETQ |l| + (LIST 0.083333333333333315 (- 0.0027777777777777779) + 7.9365079365079376E-4 (- 5.9523809523809529E-4) + 8.4175084175084182E-4 (- 0.0019175269175269176) + 0.0064102564102564109)) + ((LAMBDA (|m| |bfVar#7| |el|) + (LOOP + (COND + ((OR (> |m| 7) (ATOM |bfVar#7|) + (PROGN (SETQ |el| (CAR |bfVar#7|)) NIL)) + (RETURN NIL)) + ('T + (PROGN + (SETQ |zterm| (* |zterm| |zsquaredinv|)) + (SETQ |sum| (+ |sum| (* |el| |zterm|)))))) + (SETQ |m| (+ |m| 1)) + (SETQ |bfVar#7| (CDR |bfVar#7|)))) + 1 |l| NIL) + |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) +; + +(EVAL-WHEN (EVAL LOAD) + (SETQ |$PsiAsymptoticBern| + (VECTOR 0.0 0.16666666666666669 (- 0.033333333333333333) + 0.023809523809523812 (- 0.033333333333333333) + 0.07575757575757576 (- 0.2531135531135531) + 1.1666666666666672 (- 7.0921568627450986) + 54.971177944862163 (- 529.12424242424242) + 6192.123188405797 (- 86580.253113553103) + 1425517.166666667 (- 2.729823106781609E7) + 6.015808739006424E8 (- 1.5116315767092161E10) + 4.2961464306116675E11 (- 1.371165520508833E13) + 4.8833231897359325E14 (- 1.9296579341940072E16) + 8.4169304757368269E17 (- 4.0338071854059463E19)))) + +;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 + +(DEFUN |PsiAsymptotic| (|n| |x|) + (PROG (|sum| |factterm| |xterm| |xsq| |xnp1| |xn|) + (DECLARE (SPECIAL |$PsiAsymptoticBern|)) + (RETURN + (PROGN + (SETQ |xn| (EXPT |x| |n|)) + (SETQ |xnp1| (* |xn| |x|)) + (SETQ |xsq| (* |x| |x|)) + (SETQ |xterm| |xsq|) + (SETQ |factterm| + (/ (/ (|rgamma| (+ |n| 2)) 2.0) + (|rgamma| (|float| (+ |n| 1))))) + (SETQ |sum| + (/ (* (AREF |$PsiAsymptoticBern| 1) |factterm|) |xterm|)) + ((LAMBDA (|k|) + (LOOP + (COND + ((> |k| 22) (RETURN NIL)) + ('T + (PROGN + (SETQ |xterm| (* |xterm| |xsq|)) + (COND + ((EQL |n| 0) + (SETQ |factterm| (/ 1.0 (|float| (* 2 |k|))))) + ((EQL |n| 1) (SETQ |factterm| 1)) + ('T + (SETQ |factterm| + (/ (* (* |factterm| + (|float| (- (+ (* 2 |k|) |n|) 1))) + (|float| (- (+ (* 2 |k|) |n|) 2))) + (* (|float| (* 2 |k|)) + (|float| (- (* 2 |k|) 1))))))) + (SETQ |sum| + (+ |sum| + (/ (* (AREF |$PsiAsymptoticBern| |k|) + |factterm|) + |xterm|)))))) + (SETQ |k| (+ |k| 1)))) + 2) + (+ (+ (|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 + +(DEFUN |PsiEps| (|n| |x|) + (PROG (|result|) + (RETURN + (PROGN + (COND + ((EQL |n| 0) (SETQ |result| (- (LOG |x|)))) + ('T (SETQ |result| (/ 1.0 (* (|float| |n|) (EXPT |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 + +(DEFUN |PsiAsymptoticOrder| (|n| |x| |nterms|) + (PROG (|np1| |xterm| |sum|) + (RETURN + (PROGN + (SETQ |sum| 0) + (SETQ |xterm| 1.0) + (SETQ |np1| (+ |n| 1)) + ((LAMBDA (|k|) + (LOOP + (COND + ((> |k| |nterms|) (RETURN NIL)) + ('T + (PROGN + (SETQ |xterm| (EXPT (+ |x| (|float| |k|)) |np1|)) + (SETQ |sum| (+ |sum| (/ 1.0 |xterm|)))))) + (SETQ |k| (+ |k| 1)))) + 0) + |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*((dfPi**(n+1))*cotdiffeval(n,dfPi*(-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) + +(DEFUN |rPsi| (|n| |x|) + (PROG (|skipit| |sign| |m|) + (RETURN + (COND + ((NOT (< 0.0 |x|)) + (COND + ((ZEROP (|fracpart| |x|)) + (|FloatError| "singularity encountered at ~D" |x|)) + ('T (SETQ |m| (MOD |n| 2)) (SETQ |sign| (EXPT (- 1) |m|)) + (COND + ((EQUAL (|fracpart| |x|) 0.5) (SETQ |skipit| 1)) + ('T (SETQ |skipit| 0))) + (* |sign| + (+ (* (EXPT |dfPi| (+ |n| 1)) + (|cotdiffeval| |n| (* |dfPi| (- |x|)) |skipit|)) + (|rPsi| |n| (- 1.0 |x|))))))) + ((EQL |n| 0) (- (|rPsiW| |n| |x|))) + ('T + (* (* (|rgamma| (|float| (+ |n| 1))) (|rPsiW| |n| |x|)) + (EXPT (- 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) + +(DEFUN |rPsiW| (|n| |x|) + (PROG (|bign| |fln| |c| |a| |xmin| |beta| |alpha| |nd|) + (RETURN + (PROGN + (COND + ((OR (NOT (< 0 |x|)) (MINUSP |n|)) + (|HardError| + "rPsiW not implemented for negative n or non-positive x"))) + (SETQ |nd| 6) + (SETQ |alpha| (+ 3.5 (* 0.40000000000000002 |nd|))) + (SETQ |beta| + (+ (+ 0.20999999999999999 + (* 8.6770000000000001E-6 (- |nd| 3))) + (* 6.0380000000000002E-8 (EXPT (- |nd| 3) 2)))) + (SETQ |xmin| + (|float| (+ (FLOOR (+ |alpha| (* |beta| |n|))) 1))) + (COND + ((< 0 |n|) + (SETQ |a| + (MIN 0 + (* (/ 1.0 (|float| |n|)) + (LOG (/ DOUBLE-FLOAT-EPSILON (MIN 1.0 |x|)))))) + (SETQ |c| (EXP |a|)) + (COND + ((NOT (< (ABS |a|) 0.001)) + (SETQ |fln| (* (/ |x| |c|) (- 1.0 |c|)))) + ('T (SETQ |fln| (- (/ (* |x| |a|) |c|))))) + (SETQ |bign| (+ (FLOOR |fln|) 1)) + (COND + ((AND (< |bign| 15) (< (+ 7.0 |x|) |xmin|)) + (RETURN (|PsiAsymptoticOrder| |n| |x| |bign|)))))) + (COND + ((NOT (< |x| |xmin|)) (RETURN (|PsiAsymptotic| |n| |x|)))) + (|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 + +(DEFUN |PsiBack| (|n| |x| |xmin|) + (PROG (|result| |x0| |xintpart|) + (RETURN + (PROGN + (SETQ |xintpart| (|PsiIntpart| |x|)) + (SETQ |x0| (- |x| |xintpart|)) + (SETQ |result| (|PsiAsymptotic| |n| (+ (+ |x0| |xmin|) 1.0))) + ((LAMBDA (|bfVar#8| |k|) + (LOOP + (COND + ((COND + ((MINUSP |bfVar#8|) (< |k| |xintpart|)) + (T (> |k| |xintpart|))) + (RETURN NIL)) + ('T + (SETQ |result| + (+ |result| + (/ 1.0 + (EXPT (+ |x0| (|float| |k|)) (+ |n| 1))))))) + (SETQ |k| (+ |k| |bfVar#8|)))) + (- 1) |xmin|) + |result|)))) + +;PsiIntpart(x) == +; if x<0 +; then +; result := -PsiInpart(-x) +; else +; result := FLOOR(x) +; return result + +(DEFUN |PsiIntpart| (|x|) + (PROG (|result|) + (RETURN + (PROGN + (COND + ((MINUSP |x|) (SETQ |result| (- (|PsiInpart| (- |x|))))) + ('T (SETQ |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 + +(DEFUN |cotdiffeval| (|n| |z| |skipit|) + (PROG (|s| |sq| |v| |t2| |t1| |m| |a|) + (RETURN + (PROGN + (SETQ |a| (MAKE-ARRAY (+ |n| 2))) + (SETF (AREF |a| 0) 0.0) + (SETF (AREF |a| 1) 1.0) + ((LAMBDA (|i|) + (LOOP + (COND + ((> |i| |n|) (RETURN NIL)) + ('T (SETF (AREF |a| |i|) 0.0))) + (SETQ |i| (+ |i| 1)))) + 2) + ((LAMBDA (|l|) + (LOOP + (COND + ((> |l| |n|) (RETURN NIL)) + ('T + (PROGN + (SETQ |m| (MOD (+ |l| 1) 2)) + ((LAMBDA (|bfVar#9| |k|) + (LOOP + (COND + ((> |k| |bfVar#9|) (RETURN NIL)) + ('T + (PROGN + (COND + ((< |k| 1) (SETQ |t1| 0)) + ('T + (SETQ |t1| + (- + (* (AREF |a| (- |k| 1)) + (- |k| 1)))))) + (COND + ((< |l| |k|) (SETQ |t2| 0)) + ('T + (SETQ |t2| + (- + (* (AREF |a| (+ |k| 1)) + (+ |k| 1)))))) + (SETF (AREF |a| |k|) (+ |t1| |t2|))))) + (SETQ |k| (+ |k| 2)))) + (+ |l| 1) |m|)))) + (SETQ |l| (+ |l| 1)))) + 1) + (SETQ |v| (COT |z|)) + (SETQ |sq| (EXPT |v| 2)) + (SETQ |s| (AREF |a| (+ |n| 1))) + ((LAMBDA (|bfVar#10| |i|) + (LOOP + (COND + ((COND ((MINUSP |bfVar#10|) (< |i| 0)) (T (> |i| 0))) + (RETURN NIL)) + ('T (SETQ |s| (+ (* |s| |sq|) (AREF |a| |i|))))) + (SETQ |i| (+ |i| |bfVar#10|)))) + (- 2) (- |n| 1)) + (SETQ |m| (MOD |n| 2)) + (COND ((EQL |m| 0) (SETQ |s| (* |s| |v|)))) + (COND + ((EQL |skipit| 1) + (COND ((EQL |m| 0) (RETURN 0)) ('T (RETURN (AREF |a| 0))))) + ('T (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)) + +(DEFUN |cPsi| (|n| |z|) + (PROG (|result| |m| |bound| |nterms| |conjresult| |y| |x|) + (RETURN + (PROGN + (SETQ |x| (REALPART |z|)) + (SETQ |y| (IMAGPART |z|)) + (COND ((ZEROP |y|) (RETURN (|rPsi| |n| |x|)))) + (COND + ((< |y| 0.0) + (SETQ |conjresult| (|cPsi| |n| (COMPLEX |x| (- |y|)))) + (RETURN + (COMPLEX (REALPART |conjresult|) + (- (IMAGPART |conjresult|)))))) + (SETQ |nterms| 22) + (SETQ |bound| 10.0) + (COND + ((< |x| 0.0) + (|FloatError| "Psi implementation can't compute at ~S " + (LIST |n| |z|))) + ((AND (< 0.0 |x|) (< |bound| (ABS |z|))) + (RETURN (|PsiXotic| |n| (|PsiAsymptotic| |n| |z|)))) + ('T + (SETQ |m| + (CEILING (+ (- (SQRT (- (* |bound| |bound|) + (* |y| |y|))) + |x|) + 1.0))) + (SETQ |result| (COMPLEX 0.0 0.0)) + ((LAMBDA (|bfVar#11| |k|) + (LOOP + (COND + ((> |k| |bfVar#11|) (RETURN NIL)) + ('T + (SETQ |result| + (+ |result| + (/ 1.0 + (EXPT (+ |z| (|float| |k|)) (+ |n| 1))))))) + (SETQ |k| (+ |k| 1)))) + (- |m| 1) 0) + (RETURN + (|PsiXotic| |n| + (+ |result| (|PsiAsymptotic| |n| (+ |z| |m|))))))))))) + +;PsiXotic(n,result) == +; rgamma(float(n+1))*(-1)**MOD(n+1,2)*result + +(DEFUN |PsiXotic| (|n| |result|) + (PROG () + (RETURN + (* (* (|rgamma| (|float| (+ |n| 1))) + (EXPT (- 1) (MOD (+ |n| 1) 2))) + |result|)))) + +;cpsireflect(n,z) == +; m := MOD(n,2) +; sign := (-1)**m +; sign*dfPi**(n+1)*cotdiffeval(n,dfPi*z,0) + cPsi(n,1.0-z) + +(DEFUN |cpsireflect| (|n| |z|) + (PROG (|sign| |m|) + (RETURN + (PROGN + (SETQ |m| (MOD |n| 2)) + (SETQ |sign| (EXPT (- 1) |m|)) + (+ (* (* |sign| (EXPT |dfPi| (+ |n| 1))) + (|cotdiffeval| |n| (* |dfPi| |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 := 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 + +(DEFUN |chebf01| (|c| |z|) + (PROG (|cc| |temp| |b| |p| |sum| |rho| |divfac| |c1| |x1| |arr| + |ncount| |z1| |a1| |a2| |a3| |n2| |n1| |start| |four| |w| + |n|) + (RETURN + (PROGN + (SETQ |n| 75) + (SETQ |w| (* 2.0 |z|)) + (SETQ |four| 4.0) + (SETQ |start| (EXPT 10.0 (- 200))) + (SETQ |n1| (+ |n| 1)) + (SETQ |n2| (+ |n| 2)) + (SETQ |a3| 0.0) + (SETQ |a2| 0.0) + (SETQ |a1| |start|) + (SETQ |z1| (/ |four| |w|)) + (SETQ |ncount| |n1|) + (SETQ |arr| (MAKE-ARRAY |n2|)) + (SETF (AREF |arr| |ncount|) |start|) + (SETQ |x1| |n2|) + (SETQ |c1| (- 1.0 |c|)) + ((LAMBDA (|bfVar#12| |ncount|) + (LOOP + (COND + ((COND + ((MINUSP |bfVar#12|) (< |ncount| 0)) + (T (> |ncount| 0))) + (RETURN NIL)) + ('T + (PROGN + (SETQ |divfac| (/ 1.0 |x1|)) + (SETQ |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|)))) + (SETQ |a3| |a2|) + (SETQ |a2| |a1|) + (SETQ |a1| (AREF |arr| |ncount|))))) + (SETQ |ncount| (+ |ncount| |bfVar#12|)))) + (- 1) |n|) + (SETF (AREF |arr| 0) (/ (AREF |arr| 0) 2.0)) + (SETQ |rho| (AREF |arr| 0)) + (SETQ |sum| |rho|) + (SETQ |p| 1.0) + ((LAMBDA (|i|) + (LOOP + (COND + ((> |i| |n1|) (RETURN NIL)) + ('T + (PROGN + (SETQ |rho| (- |rho| (* |p| (AREF |arr| |i|)))) + (SETQ |sum| (+ |sum| (AREF |arr| |i|))) + (SETQ |p| (- |p|))))) + (SETQ |i| (+ |i| 1)))) + 1) + ((LAMBDA (|l|) + (LOOP + (COND + ((> |l| |n1|) (RETURN NIL)) + ('T (SETF (AREF |arr| |l|) (/ (AREF |arr| |l|) |rho|)))) + (SETQ |l| (+ |l| 1)))) + 0) + (SETQ |sum| (/ |sum| |rho|)) + (SETQ |b| 0.0) + (SETQ |temp| 0.0) + ((LAMBDA (|bfVar#13| |i|) + (LOOP + (COND + ((COND ((MINUSP |bfVar#13|) (< |i| 0)) (T (> |i| 0))) + (RETURN NIL)) + ('T + (PROGN + (SETQ |cc| |b|) + (SETQ |b| |temp|) + (SETQ |temp| (+ (- |cc|) (AREF |arr| |i|)))))) + (SETQ |i| (+ |i| |bfVar#13|)))) + (- 1) (+ |n| 1)) + |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 + +(DEFUN |brutef01| (|c| |z|) + (PROG (|maxnterms| |oldsum| |n| |sum| |term|) + (RETURN + (PROGN + (SETQ |term| 1.0) + (SETQ |sum| |term|) + (SETQ |n| 0.0) + (SETQ |oldsum| 0.0) + (SETQ |maxnterms| 10000) + ((LAMBDA (|i| |bfVar#14|) + (LOOP + (COND + ((OR (> |i| |maxnterms|) |bfVar#14|) (RETURN NIL)) + ('T + (PROGN + (SETQ |oldsum| |sum|) + (SETQ |term| + (/ (/ (* |term| |z|) (+ |c| |n|)) (+ |n| 1.0))) + (SETQ |sum| (+ |sum| |term|)) + (SETQ |n| (+ |n| 1.0))))) + (SETQ |i| (+ |i| 1)) + (SETQ |bfVar#14| (EQUAL |oldsum| |sum|)))) + 0 NIL) + |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*dfPi +;--- 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]) + +(DEFUN |f01| (|c| |z|) + (PROG () + (RETURN + (COND + ((|negintp| |c|) + (|FloatError| + "0F1 not defined for negative integer parameter value ~S" + |c|)) + ((AND (< (ABS |c|) 1000.0) (< (ABS |z|) 1000.0)) + (|brutef01| |c| |z|)) + ((AND (ZEROP (IMAGPART |z|)) (ZEROP (IMAGPART |c|)) + (NOT (< |z| 0.0)) (NOT (< |c| 0.0))) + (|brutef01| |c| |z|)) + ((AND (< (ABS |z|) (* 10.0 (ABS |c|))) (< (ABS |c|) 10000.0) + (< (ABS |z|) 10000.0)) + (|brutef01| |c| |z|)) + ('T (|FloatError| "0F1 not implemented for ~S" (LIST |c| |z|))))))) + +;--- c parameter to 0F1 +;--- w scale factor so that 0 |ncount| 0))) + (RETURN NIL)) + ('T + (PROGN + (SETQ |divfac| (/ 1.0 |x1|)) + (SETQ |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|)))) + (SETQ |a3| |a2|) + (SETQ |a2| |a1|) + (SETQ |a1| (AREF |arr| |ncount|))))) + (SETQ |ncount| (+ |ncount| |bfVar#15|)))) + (- 1) |n|) + (SETF (AREF |arr| 0) (/ (AREF |arr| 0) 2.0)) + (SETQ |rho| (AREF |arr| 0)) + (SETQ |sum| |rho|) + (SETQ |p| 1.0) + ((LAMBDA (|i|) + (LOOP + (COND + ((> |i| |n1|) (RETURN NIL)) + ('T + (PROGN + (SETQ |rho| (- |rho| (* |p| (AREF |arr| |i|)))) + (SETQ |sum| (+ |sum| (AREF |arr| |i|))) + (SETQ |p| (- |p|))))) + (SETQ |i| (+ |i| 1)))) + 1) + ((LAMBDA (|l|) + (LOOP + (COND + ((> |l| |n1|) (RETURN NIL)) + ('T (SETF (AREF |arr| |l|) (/ (AREF |arr| |l|) |rho|)))) + (SETQ |l| (+ |l| 1)))) + 0) + (SETQ |sum| (/ |sum| |rho|)) + (RETURN (LIST |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 + +(DEFUN |chebeval| (|coeflist| |x|) + (PROG (|c| |y| |temp| |b|) + (RETURN + (PROGN + (SETQ |b| 0) + (SETQ |temp| 0) + (SETQ |y| (* 2 |x|)) + ((LAMBDA (|bfVar#16| |el|) + (LOOP + (COND + ((OR (ATOM |bfVar#16|) + (PROGN (SETQ |el| (CAR |bfVar#16|)) NIL)) + (RETURN NIL)) + ('T + (PROGN + (SETQ |c| |b|) + (SETQ |b| |temp|) + (SETQ |temp| (+ (- (* |y| |b|) |c|) |el|))))) + (SETQ |bfVar#16| (CDR |bfVar#16|)))) + |coeflist| NIL) + (/ (- |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 + +(DEFUN |chebevalarr| (|coefarr| |x| |n|) + (PROG (|c| |y| |temp| |b|) + (RETURN + (PROGN + (SETQ |b| 0) + (SETQ |temp| 0) + (SETQ |y| (* 2 |x|)) + ((LAMBDA (|i|) + (LOOP + (COND + ((> |i| |n|) (RETURN NIL)) + ('T + (PROGN + (SETQ |c| |b|) + (SETQ |b| |temp|) + (SETQ |temp| + (+ (- (* |y| |b|) |c|) (ELT |coefarr| |i|)))))) + (SETQ |i| (+ |i| 1)))) + 1) + (/ (- |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 + +(DEFUN |chebstareval| (|coeflist| |x|) + (PROG (|c| |y| |temp| |b|) + (RETURN + (PROGN + (SETQ |b| 0) + (SETQ |temp| 0) + (SETQ |y| (* 2 (- (* 2 |x|) 1))) + ((LAMBDA (|bfVar#17| |el|) + (LOOP + (COND + ((OR (ATOM |bfVar#17|) + (PROGN (SETQ |el| (CAR |bfVar#17|)) NIL)) + (RETURN NIL)) + ('T + (PROGN + (SETQ |c| |b|) + (SETQ |b| |temp|) + (SETQ |temp| (+ (- (* |y| |b|) |c|) |el|))))) + (SETQ |bfVar#17| (CDR |bfVar#17|)))) + |coeflist| NIL) + (- |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 + +(DEFUN |chebstarevalarr| (|coefarr| |x| |n|) + (PROG (|c| |y| |temp| |b|) + (RETURN + (PROGN + (SETQ |b| 0) + (SETQ |temp| 0) + (SETQ |y| (* 2 (- (* 2 |x|) 1))) + ((LAMBDA (|bfVar#18| |i|) + (LOOP + (COND + ((COND ((MINUSP |bfVar#18|) (< |i| 0)) (T (> |i| 0))) + (RETURN NIL)) + ('T + (PROGN + (SETQ |c| |b|) + (SETQ |b| |temp|) + (SETQ |temp| + (+ (- (* |y| |b|) |c|) (AREF |coefarr| |i|)))))) + (SETQ |i| (+ |i| |bfVar#18|)))) + (- 1) (+ |n| 1)) + (- |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 --- 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]) + +(DEFUN |BesselJ| (|v| |z|) + (PROG (|arr| |sum| |LETTMP#1| |vp1| |w| |arg| |rz| |rv| |n| B2 B1) + (RETURN + (PROGN + (SETQ B1 10) + (SETQ B2 10) + (SETQ |n| 50) + (COND + ((OR (AND (FLOATP |v|) (ZEROP (|fracpart| |v|)) (MINUSP |v|)) + (AND (COMPLEXP |v|) (ZEROP (IMAGPART |v|)) + (ZEROP (|fracpart| (REALPART |v|))) + (< (REALPART |v|) 0.0))) + (* (|BesselJ| (- |v|) |z|) (EXPT (- 1.0) |v|))) + ((OR (AND (FLOATP |z|) (MINUSP |z|)) + (AND (COMPLEXP |z|) (< (REALPART |z|) 0.0))) + (* (|BesselJ| |v| (- |z|)) (EXPT (- 1.0) |v|))) + ((AND (ZEROP |z|) + (OR (AND (FLOATP |v|) (NOT (< |v| 0.0))) + (AND (COMPLEXP |v|) (ZEROP (IMAGPART |v|)) + (NOT (< (REALPART |v|) 0.0))))) + (COND ((ZEROP |v|) 1.0) ('T 0.0))) + ('T + (PROGN + (SETQ |rv| (ABS |v|)) + (SETQ |rz| (ABS |z|)) + (COND + ((AND (< B1 |rz|) (< (* B2 |rv|) |rz|)) + (|BesselJAsympt| |v| |z|)) + ((AND (< B1 |rv|) (< (* B2 |rz|) |rv|)) + (|BesselJAsymptOrder| |v| |z|)) + ((AND (< |rz| B1) (< |rv| B1)) + (PROGN + (SETQ |arg| (- (/ (* |z| |z|) 4.0))) + (SETQ |w| (* 2.0 |arg|)) + (SETQ |vp1| (+ |v| 1.0)) + (SETQ |LETTMP#1| (|chebf01coefmake| |vp1| |w| |n|)) + (SETQ |sum| (CAR |LETTMP#1|)) + (SETQ |arr| (CADR |LETTMP#1|)) + ((LAMBDA () + (LOOP + (COND + ((= |sum| |sum|) (RETURN NIL)) + ('T + (PROGN + (SETQ |n| (FLOOR (/ |n| 2))) + (SETQ |LETTMP#1| + (|chebf01coefmake| |vp1| |w| |n|)) + (SETQ |sum| (CAR |LETTMP#1|)) + (SETQ |arr| (CADR |LETTMP#1|)) + |LETTMP#1|)))))) + (* (/ (|chebstarevalarr| |arr| (/ |arg| |w|) |n|) + (|cgamma| |vp1|)) + (EXPT (/ |z| 2.0) |v|)))) + (T (|BesselJRecur| |v| |z|)) + ('T + (|FloatError| "BesselJ not implemented for ~S" + (LIST |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) |i| 0))) + (RETURN NIL)) + ('T + (SETF (AREF |w| |i|) + (- (/ (* (* 2.0 (+ (+ |v| |i|) 1.0)) + (AREF |w| (+ |i| 1))) + |z|) + (AREF |w| (+ |i| 2)))))) + (SETQ |i| (+ |i| |bfVar#19|)))) + (- 1) (- |m| 3)) + (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]) + +(DEFUN |BesselI| (|v| |z|) + (PROG (|chebterms| B2 B1) + (RETURN + (PROGN + (SETQ B1 15.0) + (SETQ B2 10.0) + (COND + ((AND (ZEROP |z|) (FLOATP |v|) (NOT (< |v| 0.0))) + (COND ((ZEROP |v|) 1.0) ('T 0.0))) + ((AND (FLOATP |v|) (ZEROP (|fracpart| |v|)) (MINUSP |v|)) + (|BesselI| (- |v|) |z|)) + ((< (REALPART |z|) 0.0) + (* (|BesselI| |v| (- |z|)) (EXPT (- 1.0) |v|))) + ((AND (< (REALPART |v|) 0.0) (NULL (ZEROP (IMAGPART |v|))) + (FLOATP |z|)) + (CONJUGATE (|BesselI| (CONJUGATE |v|) |z|))) + ((< B1 (ABS |z|)) + (|FloatError| "BesselI not implemented for ~S" + (LIST |v| |z|))) + ((< B1 (ABS |v|)) + (|FloatError| "BesselI not implemented for ~S" + (LIST |v| |z|))) + ((NOT (< (REALPART |v|) 0.0)) (|besselIback| |v| |z|)) + ((< (REALPART |v|) 0.0) + (PROGN + (SETQ |chebterms| 50) + (|besselIcheb| |z| |v| |chebterms|))) + ('T + (|FloatError| "BesselI not implemented for ~S" + (LIST |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) + +(DEFUN |besselIcheb| (|z| |v| |n|) + (PROG (|result| |arr| |sum| |LETTMP#1| |vp1| |w| |arg|) + (RETURN + (PROGN + (SETQ |arg| (/ (* |z| |z|) 4.0)) + (SETQ |w| (* 2.0 |arg|)) + (SETQ |vp1| (+ |v| 1.0)) + (SETQ |LETTMP#1| (|chebf01coefmake| |vp1| |w| |n|)) + (SETQ |sum| (CAR |LETTMP#1|)) + (SETQ |arr| (CADR |LETTMP#1|)) + (SETQ |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) + +(DEFUN |besselIback| (|v| |z|) + (PROG (|result| |vp1| |tv| |n| |m| |lm| |rpv| |ipv|) + (RETURN + (PROGN + (SETQ |ipv| (IMAGPART |v|)) + (SETQ |rpv| (REALPART |v|)) + (SETQ |lm| (MULTIPLE-VALUE-LIST (FLOOR |rpv|))) + (SETQ |m| (CAR |lm|)) + (SETQ |n| (* 2 (MAX 20 (+ |m| 10)))) + (SETQ |tv| (+ (CADR |lm|) (- |v| |rpv|))) + (SETQ |vp1| (+ |tv| 1.0)) + (SETQ |result| (|BesselIBackRecur| |v| |m| |tv| |z| "I" |n|)) + (SETQ |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) + +(DEFUN |BesselIBackRecur| (|largev| |argm| |v| |z| |type| |n|) + (PROG (|v1| |pn| |ct1| |xm| |m| |m1| |val| |w| |m2| |z2| |start| + |zero| |two| |one|) + (RETURN + (PROGN + (SETQ |one| 1.0) + (SETQ |two| 2.0) + (SETQ |zero| 0.0) + (SETQ |start| (EXPT 10.0 (- 40))) + (SETQ |z2| (/ |two| |z|)) + (SETQ |m2| (+ |n| 3)) + (SETQ |w| (MAKE-ARRAY (+ |m2| 1))) + (SETF (AREF |w| |m2|) |zero|) + (COND + ((EQUAL |type| "I") (SETQ |val| |one|)) + ('T (SETQ |val| (- |one|)))) + (SETQ |m1| (+ |n| 2)) + (SETF (AREF |w| |m1|) |start|) + (SETQ |m| (+ |n| 1)) + (SETQ |xm| (|float| |m|)) + (SETQ |ct1| (* |z2| (+ |xm| |v|))) + ((LAMBDA (|bfVar#20| |m|) + (LOOP + (COND + ((COND ((MINUSP |bfVar#20|) (< |m| 1)) (T (> |m| 1))) + (RETURN NIL)) + ('T + (PROGN + (SETF (AREF |w| |m|) + (+ (* (AREF |w| (+ |m| 1)) |ct1|) + (* |val| (AREF |w| (+ |m| 2))))) + (SETQ |ct1| (- |ct1| |z2|))))) + (SETQ |m| (+ |m| |bfVar#20|)))) + (- 1) (+ |n| 1)) + (SETQ |m| (+ 1 (FLOOR (/ |n| 2)))) + (SETQ |m2| (- (+ |m| |m|) 1)) + (COND + ((EQL |v| 0) (SETQ |pn| (AREF |w| (+ |m2| 2))) + ((LAMBDA (|bfVar#21| |m2|) + (LOOP + (COND + ((COND + ((MINUSP |bfVar#21|) (< |m2| 3)) + (T (> |m2| 3))) + (RETURN NIL)) + ('T (SETQ |pn| (- (AREF |w| |m2|) (* |val| |pn|))))) + (SETQ |m2| (+ |m2| |bfVar#21|)))) + (- 2) (- (* 2 |m|) 1)) + (SETQ |pn| (- (AREF |w| 1) (* |val| (+ |pn| |pn|))))) + ('T (SETQ |v1| (- |v| |one|)) (SETQ |xm| (|float| |m|)) + (SETQ |ct1| (+ (+ |v| |xm|) |xm|)) + (SETQ |pn| (* |ct1| (AREF |w| (+ |m2| 2)))) + ((LAMBDA (|bfVar#22| |m2|) + (LOOP + (COND + ((COND + ((MINUSP |bfVar#22|) (< |m2| 3)) + (T (> |m2| 3))) + (RETURN NIL)) + ('T + (PROGN + (SETQ |ct1| (- |ct1| |two|)) + (SETQ |pn| + (- (* |ct1| (AREF |w| |m2|)) + (* (/ (* |val| |pn|) |xm|) (+ |v1| |xm|)))) + (SETQ |xm| (- |xm| |one|))))) + (SETQ |m2| (+ |m2| |bfVar#22|)))) + (- 2) (- (+ |m| |m|) 1)) + (SETQ |pn| (- (AREF |w| 1) (* |val| |pn|))))) + (SETQ |m1| (+ |n| 2)) + ((LAMBDA (|m|) + (LOOP + (COND + ((> |m| |m1|) (RETURN NIL)) + ('T (SETF (AREF |w| |m|) (/ (AREF |w| |m|) |pn|)))) + (SETQ |m| (+ |m| 1)))) + 1) + (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)) + +(DEFUN |BesselasymptA| (|mu| |zsqr| |zfth|) + (PROG () + (RETURN + (* (/ (- |mu| 1) (* 16.0 |zsqr|)) + (+ (+ 1 (/ (- |mu| 13.0) (* 8.0 |zsqr|))) + (/ (+ (- (EXPT |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)) + +(DEFUN |BesselasymptB| (|mu| |z| |zsqr| |zfth|) + (PROG (|musqr|) + (RETURN + (PROGN + (SETQ |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 := dfPi +; 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) + +(DEFUN |BesselJAsympt| (|v| |z|) + (PROG (|zfth| |zsqr| |mu| |pi|) + (RETURN + (PROGN + (SETQ |pi| |dfPi|) + (SETQ |mu| (* (* 4.0 |v|) |v|)) + (SETQ |zsqr| (* |z| |z|)) + (SETQ |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*dfPi*z) +; EXP(z)/sqrttwopiz*(1.0 + sum1 ) +_ +; EXP(-(float(n)+.5)*dfPi*i)*EXP(-z)/sqrttwopiz*(1.0+ sum2) + +(DEFUN |BesselIAsympt| (|v| |z| |n|) + (PROG (|sqrttwopiz| |term1| |eight| |two| |fourvsq| |sum2| |sum1| + |i|) + (RETURN + (PROGN + (SETQ |i| (COMPLEX 0.0 1.0)) + (COND + ((EQUAL (REALPART |z|) 0.0) + (RETURN + (* (EXPT |i| |v|) (|BesselJ| |v| (- (IMAGPART |z|))))))) + (SETQ |sum1| 0.0) + (SETQ |sum2| 0.0) + (SETQ |fourvsq| (* 4.0 (EXPT |v| 2))) + (SETQ |two| 2.0) + (SETQ |eight| 8.0) + (SETQ |term1| 1.0) + ((LAMBDA (|r|) + (LOOP + (COND + ((> |r| |n|) (RETURN NIL)) + ('T + (PROGN + (SETQ |term1| + (- (/ (* |term1| + (- |fourvsq| + (EXPT + (- (* |two| (|float| |r|)) 1.0) 2))) + (* (* (|float| |r|) |eight|) |z|)))) + (SETQ |sum1| (+ |sum1| |term1|)) + (SETQ |sum2| (+ |sum2| (ABS |term1|)))))) + (SETQ |r| (+ |r| 1)))) + 1) + (SETQ |sqrttwopiz| (SQRT (* (* |two| |dfPi|) |z|))) + (+ (* (/ (EXP |z|) |sqrttwopiz|) (+ 1.0 |sum1|)) + (* (/ (* (EXP (- (* (* (+ (|float| |n|) 0.5) |dfPi|) |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)