diff options
Diffstat (limited to 'src/interp/sfsfun.boot.pamphlet')
| -rw-r--r-- | src/interp/sfsfun.boot.pamphlet | 1031 |
1 files changed, 1031 insertions, 0 deletions
diff --git a/src/interp/sfsfun.boot.pamphlet b/src/interp/sfsfun.boot.pamphlet new file mode 100644 index 00000000..50bd7b5b --- /dev/null +++ b/src/interp/sfsfun.boot.pamphlet @@ -0,0 +1,1031 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/interp sfsfun.boot} +\author{The Axiom Team} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +\begin{verbatim} +NOTEfrom TTT: at least BesselJAsymptOrder needs work + +1. This file contains the contents of BWC's original files: + floaterrors.boot + floatutils.boot + rgamma.boot + cgamma.boot + rpsi.boot + cpsi.boot + f01.boot + chebf01cmake.boot + chebevalsf.boot + besselIJ.boot + +2. All declarations have been commented out with "--@@" + since the boot translator is generating bad lisp code from them. + +3. The functions PsiAsymptotic, PsiEps and PsiAsymptoticOrder + had inconpatible definitions in rpsi.boot and cpsi.boot -- + the local variables were declared float in one file and COMPLEX in + the other. The type declarations have been commented out and the + duplicate definitions have been deleted. + +4. BesselIJ was not compiling. I have modified the code from that + file to make it compile. It should be checked for correctness. + +SMW June 25, 1991 + +"Fixes" to BesselJ, B. Char June 14, 1992. Needs extensive testing and + further fixes to BesselI and BesselJ. +More fixes to BesselJ, T. Tsikas 24 Feb, 1995. + +\end{verbatim} +\section{License} +<<license>>= +-- Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. +-- All rights reserved. +-- +-- Redistribution and use in source and binary forms, with or without +-- modification, are permitted provided that the following conditions are +-- met: +-- +-- - Redistributions of source code must retain the above copyright +-- notice, this list of conditions and the following disclaimer. +-- +-- - Redistributions in binary form must reproduce the above copyright +-- notice, this list of conditions and the following disclaimer in +-- the documentation and/or other materials provided with the +-- distribution. +-- +-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the +-- names of its contributors may be used to endorse or promote products +-- derived from this software without specific prior written permission. +-- +-- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +-- IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +-- TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +-- PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +-- OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +-- EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +-- PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +-- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +-- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +-- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +-- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +@ +<<*>>= +<<license>> + +-- Used to be SPECFNSF +)package "BOOT" + +FloatError(formatstring,arg) == +-- ERROR(formatstring,arg) + ERROR FORMAT([],formatstring,arg) + +nangenericcomplex () == + 1.0/COMPLEX(0.0) + + + +fracpart(x) == + CADR(MULTIPLE_-VALUE_-LIST(FLOOR(x))) + +intpart(x) == + CAR(MULTIPLE_-VALUE_-LIST(FLOOR(x))) + +negintp(x) == + if ZEROP IMAGPART(x) and x<0.0 and ZEROP fracpart(x) + then + true + else + false + +--- Small float implementation of Gamma function. Valid for +--- real arguments. Maximum error (relative) seems to be +--- 2-4 ulps for real x 2<x<9, and up to ten ulps for larger x +--- up to overflow. See Hart & Cheney. +--- Bruce Char, April, 1990. + +horner(l,x) == + result := 0 + for el in l repeat + result := result *x + el + return result + +rgamma (x) == + if COMPLEXP(x) then FloatError('"Gamma not implemented for complex value ~D",x) + ZEROP (x-1.0) => 1.0 + if x>20 then gammaStirling(x) else gammaRatapprox(x) + +lnrgamma (x) == + if x>20 then lnrgammaRatapprox(x) else LOG(gammaRatapprox(x)) + +cbeta(z,w) == + cgamma(z)*cgamma(w)/(cgamma(z+w)) + +gammaStirling(x) == + EXP(lnrgamma(x)) + +lnrgammaRatapprox(x) == + (x-.5)*LOG(x) - x + LOG(SQRT(2.0*PI)) + phiRatapprox(x) + +phiRatapprox(x) == + arg := 1/(x**2) + p := horner([.0666629070402007526,_ + .6450730291289920251389,_ + .670827838343321349614617,_ + .12398282342474941538685913],arg); + q := horner([1.0,7.996691123663644194772,_ + 8.09952718948975574728214,_ + 1.48779388109699298468156],arg); + result := p/(x*q) + result + +gammaRatapprox (x) == + if (x>=2 and x<=3) + then + result := gammaRatkernel(x) + else + if x>3 + then + n := FLOOR(x)-2 + a := x-n-2 + reducedarg := 2+a + prod := */[reducedarg+i for i in 0..n-1] + result := prod* gammaRatapprox(reducedarg) + else + if (2>x and x>0) + then + n := 2-FLOOR(x) + a := x-FLOOR(x) + reducedarg := 2+a + prod := */[x+i for i in 0..n-1] + result := gammaRatapprox(reducedarg)/prod + else + Pi := PI + lx := MULTIPLE_-VALUE_-LIST(FLOOR(x)) + intpartx := CAR(lx)+1 + restx := CADR(lx) + if ZEROP restx -- case of negative non-integer value + then + FloatError ('"Gamma undefined for non-positive integers: ~D",x) + result := nangenericcomplex () + else + result := Pi/(gammaRatapprox(1.0-x)*(-1.0)**(intpartx+1)*SIN(restx*Pi)) + result + +gammaRatkernel(x) == + p := horner(REVERSE([3786.01050348257245475108,_ + 2077.45979389418732098416,_ + 893.58180452374981423868,_ + 222.1123961680117948396,_ + 48.95434622790993805232,_ + 6.12606745033608429879,_ + .778079585613300575867]),x-2) + q := horner(REVERSE([3786.01050348257187258861,_ + 476.79386050368791516095,_ + -867.23098753110299445707,_ + 83.55005866791976957459,_ + 50.78847532889540973716,_ + -13.40041478578134826274,_ + 1]),x-2.0) + p/q + +-- cgamma(z) Gamma function for complex arguments. +-- Bruce Char April-May, 1990. +-- +-- Our text for complex gamma is H. Kuki's paper Complex Gamma +-- Function with Error Control", CACM vol. 15, no. 4, ppp. 262-267. +-- (April, 1972.) It uses the reflection formula and the basic +-- z+1 recurrence to transform the argument into something that +-- Stirling's asymptotic formula can handle. +-- +-- However along the way it does a few tricky things to reduce +-- problems due to roundoff/cancellation error for particular values. + +-- cgammat is auxiliary "t" function (see p. 263 Kuki) +cgammat(x) == + MAX(0.1, MIN(10.0, 10.0*SQRT(2.0) - ABS(x))) + +cgamma (z) == + z2 := IMAGPART(z) + z1 := REALPART(z) --- call real valued gamma if z is real + if ZEROP z2 + then result := rgamma(z1) + else + result := clngamma(z1,z2,z) + result := EXP(result) + result + +lncgamma(z) == + clngamma(REALPART z, IMAGPART z, z) + +clngamma(z1,z2,z) == + --- conjugate of gamma is gamma of conjugate. map 2nd and 4th quads + --- to first and third quadrants + if z1<0.0 + then if z2 > 0.0 + then result := CONJUGATE(clngammacase1(z1,-z2,COMPLEX(z1,-z2))) + else result := clngammacase1(z1,z2,z) + else if z2 < 0.0 + then result := CONJUGATE(clngammacase23(z1,-z2,_ + COMPLEX(z1,-z2))) + else result := clngammacase23(z1,z2,z) + result + +clngammacase1(z1,z2,z) == + result1 := PiMinusLogSinPi(z1,z2,z) + result2 := clngamma(1.0-z1,-z2,1.0-z) + result1-result2 + +PiMinusLogSinPi(z1,z2,z) == + cgammaG(z1,z2) - logH(z1,z2,z) + +cgammaG(z1,z2) == + LOG(2*PI) + PI*z2 - COMPLEX(0.0,1.0)*PI*(z1-.5) + +logH(z1,z2,z) == + z1bar := CADR(MULTIPLE_-VALUE_-LIST(FLOOR(z1))) ---frac part of z1 + piz1bar := PI*z1bar + piz2 := PI*z2 + twopiz2 := 2.0*piz2 + i := COMPLEX(0.0,1.0) + part2 := EXP(twopiz2)*(2.0*SIN(piz1bar)**2 + SIN(2.0*piz1bar)*i) + part1 := -TANH(piz2)*(1.0+EXP(twopiz2)) +--- part1 is another way of saying 1 - exp(2*Pi*z1bar) + LOG(part1+part2) + +clngammacase23(z1,z2,z) == + tz2 := cgammat(z2) + if (z1 < tz2) + then result:= clngammacase2(z1,z2,tz2,z) + else result:= clngammacase3(z) + result + +clngammacase2(z1,z2,tz2,z) == + n := float(CEILING(tz2-z1)) + zpn := z+n + (z-.5)*LOG(zpn) - (zpn) + cgammaBernsum(zpn) - cgammaAdjust(logS(z1,z2,z,n,zpn)) + +logS(z1,z2,z,n,zpn) == + sum := 0.0 + for k in 0..(n-1) repeat + if z1+k < 5.0 - 0.6*z2 + then sum := sum + LOG((z+k)/zpn) + else sum := sum + LOG(1.0 - (n-k)/zpn) + sum + +--- on p. 265, Kuki, logS result should have its imaginary part +--- adjusted by 2 Pi if it is negative. +cgammaAdjust(z) == + if IMAGPART(z)<0.0 + then result := z + COMPLEX(0.0, 2.0*PI) + else result := z + result + +clngammacase3(z) == + (z- .5)*LOG(z) - z + cgammaBernsum(z) + +cgammaBernsum (z) == + sum := LOG(2.0*PI)/2.0 + zterm := z + zsquaredinv := 1.0/(z*z) + l:= [.083333333333333333333, -.0027777777777777777778,_ + .00079365079365079365079, -.00059523809523809523810,_ + .00084175084175084175084, -.0019175269175269175269,_ + .0064102564102564102564] + for m in 1..7 for el in l repeat + zterm := zterm*zsquaredinv + sum := sum + el*zterm + sum + + + + +--- nth derivatives of ln gamma for real x, n = 0,1,.... +--- requires files floatutils, rgamma +$PsiAsymptoticBern := VECTOR(0.0, 0.1666666666666667, -0.03333333333333333, 0.02380952380952381,_ + -0.03333333333333333, 0.07575757575757576, -0.2531135531135531, 1.166666666666667,_ + -7.092156862745098, 54.97117794486216, -529.1242424242424, 6192.123188405797,_ + -86580.25311355311, 1425517.166666667, -27298231.06781609, 601580873.9006424,_ + -15116315767.09216, 429614643061.1667, -13711655205088.33, 488332318973593.2,_ + -19296579341940070.0, 841693047573682600.0, -40338071854059460000.0) + + +PsiAsymptotic(n,x) == + xn := x**n + xnp1 := xn*x + xsq := x*x + xterm := xsq + factterm := rgamma(n+2)/2.0/rgamma(float(n+1)) + --- initialize to 1/n! + sum := AREF($PsiAsymptoticBern,1)*factterm/xterm + for k in 2..22 repeat + xterm := xterm * xsq + if n=0 then factterm := 1.0/float(2*k) + else if n=1 then factterm := 1 + else factterm := factterm * float(2*k+n-1)*float(2*k+n-2)/(float(2*k)*float(2*k-1)) + sum := sum + AREF($PsiAsymptoticBern,k)*factterm/xterm + PsiEps(n,x) + 1.0/(2.0*xnp1) + 1.0/xn * sum + + +PsiEps(n,x) == + if n = 0 + then + result := -LOG(x) + else + result := 1.0/(float(n)*(x**n)) + result + +PsiAsymptoticOrder(n,x,nterms) == + sum := 0 + xterm := 1.0 + np1 := n+1 + for k in 0..nterms repeat + xterm := (x+float(k))**np1 + sum := sum + 1.0/xterm + sum + + +rPsi(n,x) == + if x<=0.0 + then + if ZEROP fracpart(x) + then FloatError('"singularity encountered at ~D",x) + else + m := MOD(n,2) + sign := (-1)**m + if fracpart(x)=.5 + then + skipit := 1 + else + skipit := 0 + sign*((PI**(n+1))*cotdiffeval(n,PI*(-x),skipit) + rPsi(n,1.0-x)) + else if n=0 + then + - rPsiW(n,x) + else + rgamma(float(n+1))*rPsiW(n,x)*(-1)**MOD(n+1,2) + +---Amos' w function, with w(0,x) picked to be -psi(x) for x>0 +rPsiW(n,x) == + if (x <=0 or n < 0) + then + HardError('"rPsiW not implemented for negative n or non-positive x") + nd := 6 -- magic number for number of digits in a word? + alpha := 3.5 + .40*nd + beta := 0.21 + (.008677e-3)*(nd-3) + (.0006038e-4)*(nd-3)**2 + xmin := float(FLOOR(alpha + beta*n) + 1) + if n>0 + then + a := MIN(0,1.0/float(n)*LOG(DOUBLE_-FLOAT_-EPSILON/MIN(1.0,x))) + c := EXP(a) + if ABS(a) >= 0.001 + then + fln := x/c*(1.0-c) + else + fln := -x*a/c + bign := FLOOR(fln) + 1 +--- Amos says to use alternative series for large order if ordinary +--- backwards recurrence too expensive + if (bign < 15) and (xmin > 7.0+x) + then + return PsiAsymptoticOrder(n,x,bign) + if x>= xmin + then + return PsiAsymptotic(n,x) +---ordinary case -- use backwards recursion + PsiBack(n,x,xmin) + +PsiBack(n,x,xmin) == + xintpart := PsiIntpart(x) + x0 := x-xintpart ---frac part of x + result := PsiAsymptotic(n,x0+xmin+1.0) + for k in xmin..xintpart by -1 repeat +--- Why not decrement from x? See Amos p. 498 + result := result + 1.0/((x0 + float(k))**(n+1)) + result + + +PsiIntpart(x) == + if x<0 + then + result := -PsiInpart(-x) + else + result := FLOOR(x) + return result + + +---Code for computation of derivatives of cot(z), necessary for +--- polygamma reflection formula. If you want to compute n-th derivatives of +---cot(Pi*x), you have to multiply the result of cotdiffeval by Pi**n. + +-- MCD: This is defined at the Lisp Level. +-- COT(z) == +-- 1.0/TAN(z) + +cotdiffeval(n,z,skipit) == +---skip=1 if arg z is known to be an exact multiple of Pi/2 + a := MAKE_-ARRAY(n+2) + SETF(AREF(a,0),0.0) + SETF(AREF(a,1),1.0) + for i in 2..n repeat + SETF(AREF(a,i),0.0) + for l in 1..n repeat + m := MOD(l+1,2) + for k in m..l+1 by 2 repeat + if k<1 + then + t1 := 0 + else + t1 := -AREF(a,k-1)*(k-1) + if k>l + then + t2 := 0 + else + t2 := -AREF(a,k+1)*(k+1) + SETF(AREF(a,k), t1+t2) + --- evaluate d^N/dX^N cot(z) via Horner-like rule + v := COT(z) + sq := v**2 + s := AREF(a,n+1) + for i in (n-1)..0 by -2 repeat + s := s*sq + AREF(a,i) + m := MOD(n,2) + if m=0 + then + s := s*v + if skipit=1 + then + if m=0 + then + return 0 + else + return AREF(a,0) + else + return s +--- nth derivatives of ln gamma for complex z, n=0,1,... +--- requires files rpsi (and dependents), floaterrors +--- currently defined only in right half plane until reflection formula +--- works + +--- B. Char, June, 1990. + +cPsi(n,z) == + x := REALPART(z) + y := IMAGPART(z) + if ZEROP y + then --- call real function if real + return rPsi(n,x) + if y<0.0 + then -- if imagpart(z) negative, take conjugate of conjugate + conjresult := cPsi(n,COMPLEX(x,-y)) + return COMPLEX(REALPART(conjresult),-IMAGPART(conjresult)) + nterms := 22 + bound := 10.0 + if x<0.0 --- and ABS(z)>bound and ABS(y)<bound + then + FloatError('"Psi implementation can't compute at ~S ",[n,z]) +--- return cpsireflect(n,x,y,z) + else if (x>0.0 and ABS(z)>bound ) --- or (x<0.0 and ABS(y)>bound) + then + return PsiXotic(n,PsiAsymptotic(n,z)) + else --- use recursion formula + m := CEILING(SQRT(bound*bound - y*y) - x + 1.0) + result := COMPLEX(0.0,0.0) + for k in 0..(m-1) repeat + result := result + 1.0/((z + float(k))**(n+1)) + return PsiXotic(n,result+PsiAsymptotic(n,z+m)) + +PsiXotic(n,result) == + rgamma(float(n+1))*(-1)**MOD(n+1,2)*result + +cpsireflect(n,z) == + m := MOD(n,2) + sign := (-1)**m + sign*PI**(n+1)*cotdiffeval(n,PI*z,0) + cPsi(n,1.0-z) + +--- c parameter to 0F1, possibly complex +--- z argument to 0F1 +--- Depends on files floaterror, floatutils + +--- Program transcribed from Fortran,, p. 80 Luke 1977 + +chebf01 (c,z) == +--- w scale factor so that 0<z/w<1 +--- n n+2 coefficients will be produced stored in an array +--- indexed from 0 to n+1. +--- See Luke's books for further explanation + n := 75 --- ad hoc decision +--- if ABS(z)/ABS(c) > 200.0 and ABS(z)>10000.0 +--- then +--- FloatError('"cheb0F1 not implemented for ~S < 1",[c,z]) + w := 2.0*z +--- arr will be used to store the Cheb. series coefficients + four:= 4.0 + start := EXPT(10.0, -200) + n1 := n+1 + n2 := n+2 + a3 := 0.0 + a2 := 0.0 + a1 := start -- arbitrary starting value + z1 := four/w + ncount := n1 + arr := MAKE_-ARRAY(n2) + SETF(AREF(arr,ncount) , start) -- start off + x1 := n2 + c1 := 1.0 - c + for ncount in n..0 by -1 repeat + divfac := 1.0/x1 + x1 := x1 -1.0 + SETF(AREF(arr,ncount) ,_ + x1*((divfac+z1*(x1-c1))*a1 +_ + (1.0/x1 + z1*(x1+c1+1.0))*a2-divfac*a3)) + a3 := a2 + a2 := a1 + a1 := AREF(arr,ncount) + SETF(AREF(arr,0),AREF(arr,0)/2.0) +-- compute scale factor + rho := AREF(arr,0) + sum := rho + p := 1.0 + for i in 1..n1 repeat + rho := rho - p*AREF(arr,i) + sum := sum+AREF(arr,i) + p := -p + for l in 0..n1 repeat + SETF(AREF(arr,l), AREF(arr,l)/rho) + sum := sum/rho +--- Now evaluate array at argument + b := 0.0 + temp := 0.0 + for i in (n+1)..0 by -1 repeat + cc := b + b := temp + temp := -cc + AREF(arr,i) + temp + + +brutef01(c,z) == +-- Use series definition. Won't work well if cancellation occurs + term := 1.0 + sum := term + n := 0.0 + oldsum := 0.0 + maxnterms := 10000 + for i in 0..maxnterms until oldsum=sum repeat + oldsum := sum + term := term*z/(c+n)/(n+1.0) + sum := sum + term + n := n+1.0 + sum + +f01(c,z)== + if negintp(c) + then + FloatError('"0F1 not defined for negative integer parameter value ~S",c) +-- conditions when we'll permit the computation + else if ABS(c)<1000.0 and ABS(z)<1000.0 + then + brutef01(c,z) + else if ZEROP IMAGPART(z) and ZEROP IMAGPART(c) and z>=0.0 and c>=0.0 + then + brutef01(c,z) +--- else +--- t := SQRT(-z) +--- c1 := c-1.0 +--- p := PHASE(c) +--- if ABS(c)>10.0*ABS(t) and p>=0.0 and PHASE(c)<.90*PI +--- then BesselJAsymptOrder(c1,2*t)*cgamma(c/(t**(c1))) +--- else if ABS(t)>10.0*ABS(c) and ABS(t)>50.0 +--- then BesselJAsympt(c1,2*t)*cgamma(c/(t**(c1))) +--- else +--- FloatError('"0F1 not implemented for ~S",[c,z]) + else if (10.0*ABS(c)>ABS(z)) and ABS(c)<1.0E4 and ABS(z)<1.0E4 + then + brutef01(c,z) + else + FloatError('"0F1 not implemented for ~S",[c,z]) + +--- c parameter to 0F1 +--- w scale factor so that 0<z/w<1 +--- n n+2 coefficients will be produced stored in an array +--- indexed from 0 to n+1. +--- See Luke's books for further explanation + +--- Program transcribed from Fortran,, p. 80 Luke 1977 +chebf01coefmake (c,w,n) == +--- arr will be used to store the Cheb. series coefficients + four:= 4.0 + start := EXPT(10.0, -200) + n1 := n+1 + n2 := n+2 + a3 := 0.0 + a2 := 0.0 + a1 := start -- arbitrary starting value + z1 := four/w + ncount := n1 + arr := MAKE_-ARRAY(n2) + SETF(AREF(arr,ncount) , start) -- start off + x1 := n2 + c1 := 1.0 - c + for ncount in n..0 by -1 repeat + divfac := 1.0/x1 + x1 := x1 -1.0 + SETF(AREF(arr,ncount) ,_ + x1*((divfac+z1*(x1-c1))*a1 +_ + (1.0/x1 + z1*(x1+c1+1.0))*a2-divfac*a3)) + a3 := a2 + a2 := a1 + a1 := AREF(arr,ncount) + SETF(AREF(arr,0),AREF(arr,0)/2.0) +-- compute scale factor + rho := AREF(arr,0) + sum := rho + p := 1.0 + for i in 1..n1 repeat + rho := rho - p*AREF(arr,i) + sum := sum+AREF(arr,i) + p := -p + for l in 0..n1 repeat + SETF(AREF(arr,l), AREF(arr,l)/rho) + sum := sum/rho + return([sum,arr]) + + + + +---evaluation of Chebychev series of degree n at x, where the series's +---coefficients are given by the list in descending order (coef. of highest +---power first) + +---May be numerically unstable for certain lists of coefficients; +--- could possibly reverse sequence of coefficients + +--- Cheney and Hart p. 15. + +--- B. Char, March 1990 + +chebeval (coeflist,x) == + b := 0; + temp := 0; + y := 2*x; + + for el in coeflist repeat + c := b; + b := temp + temp := y*b -c + el + (temp -c)/2 + + +chebevalarr(coefarr,x,n) == + b := 0; + temp := 0; + y := 2*x; + + for i in 1..n repeat + c := b; + b := temp + temp := y*b -c + coefarr.i + (temp -c)/2 + +--- If plist is a list of coefficients for the Chebychev approximation +--- of a function f(x), then chebderiveval computes the Chebychev approximation +--- of f'(x). See Luke, "Special Functions and their approximations, vol. 1 +--- Academic Press 1969., p. 329 (from Clenshaw and Cooper) + +--- < definition to be supplied> + +--- chebstareval(plist,n) computes a Chebychev approximation from a +--- coefficient list, using shifted Chebychev polynomials of the first kind +--- The defining relation is that T*(n,x) = T(n,2*x-1). Thus the interval +--- [0,1] of T*n is the interval [-1,1] of Tn. + +chebstareval(coeflist,x) == -- evaluation of T*(n,x) + b := 0 + temp := 0 + y := 2*(2*x-1) + + for el in coeflist repeat + c := b; + b := temp + temp := y*b -c + el + temp - y*b/2 + + +chebstarevalarr(coefarr,x,n) == -- evaluation of sum(C(n)*T*(n,x)) + + b := 0 + temp := 0 + y := 2*(2*x-1) + + for i in (n+1)..0 by -1 repeat + c := b + b := temp + temp := y*b -c + AREF(coefarr,i) + temp - y*b/2 + +--Float definitions for Bessel functions I and J. +--External references: cgamma, rgamma, chebf01coefmake, chebevalstarsf +-- floatutils + +---BesselJ works for complex and real values of v and z +BesselJ(v,z) == +---Ad hoc boundaries for approximation + B1:= 10 + B2:= 10 + n := 50 --- number of terms in Chebychev series. + --- tests for negative integer order + (FLOATP(v) and ZEROP fracpart(v) and (v<0)) or (COMPLEXP(v) and ZEROP IMAGPART(v) and ZEROP fracpart(REALPART(v)) and REALPART(v)<0.0) => + --- odd or even according to v (9.1.5 A&S) + --- $J_{-n}(z)=(-1)^{n} J_{n}(z)$ + BesselJ(-v,z)*EXPT(-1.0,v) + (FLOATP(z) and (z<0)) or (COMPLEXP(z) and REALPART(z)<0.0) => + --- negative argument (9.1.35 A&S) + --- $J_{\nu}(z e^{m \pi i}) = e^{m \nu \pi i} J_{\nu}(z)$ + BesselJ(v,-z)*EXPT(-1.0,v) + ZEROP z and ((FLOATP(v) and (v>=0.0)) or (COMPLEXP(v) and + ZEROP IMAGPART(v) and REALPART(v)>=0.0)) => --- zero arg, pos. real order + ZEROP v => 1.0 --- J(0,0)=1 + 0.0 --- J(v,0)=0 for real v>0 + rv := ABS(v) + rz := ABS(z) + (rz>B1) and (rz > B2*rv) => --- asymptotic argument + BesselJAsympt(v,z) + (rv>B1) and (rv > B2*rz) => --- asymptotic order + BesselJAsymptOrder(v,z) + (rz< B1) and (rv<B1) => --- small order and argument + arg := -(z*z)/4.0 + w := 2.0*arg + vp1 := v+1.0 + [sum,arr] := chebf01coefmake(vp1,w,n) + ---if we get NaNs then half n + while not _=(sum,sum) repeat + n:=FLOOR(n/2) + [sum,arr] := chebf01coefmake(vp1,w,n) + ---now n is safe, can we increase it (we know that 2*n is bad)? + chebstarevalarr(arr,arg/w,n)/cgamma(vp1)*EXPT(z/2.0,v) + true => BesselJRecur(v,z) + FloatError('"BesselJ not implemented for ~S", [v,z]) + +BesselJRecur(v,z) == + -- boost order + --Numerical.Recipes. suggest so:=v+sqrt(n.s.f.^2*v) + so:=15.0*z + -- reduce order until non-zero + while ZEROP ABS(BesselJAsymptOrder(so,z)) repeat so:=so/2.0 + if ABS(so)<ABS(z) then so:=v+18.*SQRT(v) + m:= FLOOR(ABS(so-v))+1 + w:=MAKE_-ARRAY(m) + SETF(AREF(w,m-1),BesselJAsymptOrder(v+m-1,z)) + SETF(AREF(w,m-2),BesselJAsymptOrder(v+m-2,z)) + for i in m-3 .. 0 by -1 repeat + SETF(AREF(w,i), 2.0 * (v+i+1.0) * AREF(w,i+1) /z -AREF(w,i+2)) + AREF(w,0) + +BesselI(v,z) == + B1 := 15.0 + B2 := 10.0 + ZEROP(z) and FLOATP(v) and (v>=0.0) => --- zero arg, pos. real order + ZEROP(v) => 1.0 --- I(0,0)=1 + 0.0 --- I(v,0)=0 for real v>0 +--- Transformations for negative integer orders + FLOATP(v) and ZEROP(fracpart(v)) and (v<0) => BesselI(-v,z) +--- Halfplane transformations for Re(z)<0 + REALPART(z)<0.0 => BesselI(v,-z)*EXPT(-1.0,v) +--- Conjugation for complex order and real argument + REALPART(v)<0.0 and not ZEROP IMAGPART(v) and FLOATP(z) => + CONJUGATE(BesselI(CONJUGATE(v),z)) +---We now know that Re(z)>= 0.0 + ABS(z) > B1 => --- asymptotic argument case + FloatError('"BesselI not implemented for ~S",[v,z]) + ABS(v) > B1 => + FloatError('"BesselI not implemented for ~S",[v,z]) +--- case of small argument and order + REALPART(v)>= 0.0 => besselIback(v,z) + REALPART(v)< 0.0 => + chebterms := 50 + besselIcheb(z,v,chebterms) + FloatError('"BesselI not implemented for ~S",[v,z]) + +--- Compute n terms of the chebychev series for f01 +besselIcheb(z,v,n) == + arg := (z*z)/4.0 + w := 2.0*arg; + vp1 := v+1.0; + [sum,arr] := chebf01coefmake(vp1,w,n) + result := chebstarevalarr(arr,arg/w,n)/cgamma(vp1)*EXPT(z/2.0,v) + +besselIback(v,z) == + ipv := IMAGPART(v) + rpv := REALPART(v) + lm := MULTIPLE_-VALUE_-LIST(FLOOR(rpv)) + m := CAR(lm) --- floor of real part of v + n := 2*MAX(20,m+10) --- how large the back recurrence should be + tv := CADR(lm)+(v-rpv) --- fractional part of real part of v + --- plus imaginary part of v + vp1 := tv+1.0; + result := BesselIBackRecur(v,m,tv,z,'"I",n) + result := result/cgamma(vp1)*EXPT(z/2.0,tv) + +--- Backward recurrence for Bessel functions. Luke (1975), p. 247. +--- works for -Pi< arg z <= Pi and -Pi < arg v <= Pi +BesselIBackRecur(largev,argm,v,z,type,n) == +--- v + m = largev + one := 1.0 + two := 2.0 + zero := 0.0 + start := EXPT(10.0,-40) + z2 := two/z + m2 := n+3 + w:=MAKE_-ARRAY(m2+1) + SETF(AREF(w,m2), zero) --- start off + if type = '"I" + then + val := one + else + val := -one + m1 := n+2 + SETF(AREF(w,m1), start) + m := n+1 + xm := float(m) + ct1 := z2*(xm+v) + --- initialize + for m in (n+1)..1 by -1 repeat + SETF(AREF(w,m), AREF(w,m+1)*ct1 + val*AREF(w,m+2)) + ct1 := ct1 - z2 + m := 1 + FLOOR(n/2) + m2 := m + m -1 + if (v=0) + then + pn := AREF(w, m2 + 2) + for m2 in (2*m-1)..3 by -2 repeat + pn := AREF(w, m2) - val *pn + pn := AREF(w,1) - val*(pn+pn) + else + v1 := v-one + xm := float(m) + ct1 := v + xm + xm + pn := ct1*AREF(w, m2 + 2) + for m2 in (m+m -1)..3 by -2 repeat + ct1 := ct1 - two + pn := ct1*AREF(w,m2) - val*pn/xm*(v1+xm) + xm := xm - one + pn := AREF(w,1) - val * pn + m1 := n+2 + for m in 1..m1 repeat + SETF(AREF(w,m), AREF(w,m)/pn) + AREF(w,argm+1) + + + + +---Asymptotic functions for large values of z. See p. 204 Luke 1969 vol. 1. + +--- mu is 4*v**2 +--- zsqr is z**2 +--- zfth is z**4 + +BesselasymptA(mu,zsqr,zfth) == + (mu -1)/(16.0*zsqr) * (1 + (mu - 13.0)/(8.0*zsqr) + _ + (mu**2 - 53.0*mu + 412.0)/(48.0*zfth)) + +BesselasymptB(mu,z,zsqr,zfth) == + musqr := mu*mu + z + (mu-1.0)/(8.0*z) *(1.0 + (mu - 25.0)/(48.0*zsqr) + _ + (musqr - 114.0*mu + 1073.0)/(640.0*zfth) +_ + (5.0*mu*musqr - 1535.0*musqr + 54703.0*mu - 375733.0)/(128.0*zsqr*zfth)) + +--- Asymptotic series only works when |v| < |z|. + +BesselJAsympt (v,z) == + pi := PI + mu := 4.0*v*v + zsqr := z*z + zfth := zsqr*zsqr + SQRT(2.0/(pi*z))*EXP(BesselasymptA(mu,zsqr,zfth))*_ + COS(BesselasymptB(mu,z,zsqr,zfth) - pi*v/2.0 - pi/4.0) + + +---Asymptotic series for I. See Whittaker, p. 373. +--- valid for -3/2 Pi < arg z < 1/2 Pi + +BesselIAsympt(v,z,n) == + i := COMPLEX(0.0, 1.0) + if (REALPART(z) = 0.0) + then return EXPT(i,v)*BesselJ(v,-IMAGPART(z)) + sum1 := 0.0 + sum2 := 0.0 + fourvsq := 4.0*v**2 + two := 2.0 + eight := 8.0 + term1 := 1.0 +--- sum1, sum2, fourvsq,two,i,eight,term1]) + for r in 1..n repeat + term1 := -term1 *(fourvsq-(two*float(r)-1.0)**2)/_ + (float(r)*eight*z) + sum1 := sum1 + term1 + sum2 := sum2 + ABS(term1) + sqrttwopiz := SQRT(two*PI*z) + EXP(z)/sqrttwopiz*(1.0 + sum1 ) +_ + EXP(-(float(n)+.5)*PI*i)*EXP(-z)/sqrttwopiz*(1.0+ sum2) + + +---Asymptotic formula for BesselJ when order is large comes from +---Debye (1909). See Olver, Asymptotics and Special Functions, p. 134. +---Expansion good for 0<=phase(v)<Pi +---A&S recommend "uniform expansion" with complicated coefficients and Airy function. +---Debye's Formula is in 9.3.7,9.3.9,9.3.10 of A&S +---AXIOM recurrence for u_{k} +---f(0)==1::EXPR INT +---f(n)== (t^2)*(1-t^2)*D(f(n-1),t)/2 + (1/8)*integrate( (1-5*t^2)*f(n-1),t) +BesselJAsymptOrder(v,z) == + sechalpha := z/v + alpha := ACOSH(1.0/sechalpha) + tanhalpha := SQRT(1.0-(sechalpha*sechalpha)) + -- cothalpha := 1.0/tanhalpha + ca := 1.0/tanhalpha + + Pi := PI + ca2:=ca*ca + ca4:=ca2*ca2 + ca8:=ca4*ca4 + EXP(-v*(alpha-tanhalpha))/SQRT(2.0*Pi*v*tanhalpha)*_ + (1.0+_ + horner([ -5.0, 3.0],_ + ca2)*ca/(v*24.0)+_ + horner([ 385.0, -462.0, 81.0],_ + ca2)*ca2/(1152.0*v*v)+_ + horner([ -425425.0, 765765.0, -369603.0, 30375.0],_ + ca2)*ca2*ca/(414720.0*v*v*v)+_ + horner([ 185910725.0, -446185740.0, 349922430.0, -94121676.0, 4465125.0],_ + ca2)*ca4/(39813120.0*v*v*v*v)+_ + horner([ -188699385875.0, 566098157625.0, -614135872350.0, 284499769554.0, -49286948607.0, 1519035525.0],_ + ca2)*ca4*ca/(6688604160.0*v*v*v*v*v)+_ + horner([1023694168371875.0,-3685299006138750.0,5104696716244125.0,-3369032068261860.0,1050760774457901.0,-127577298354750.0,2757049477875.0],_ + ca2)*ca4*ca2/(4815794995200.0*v*v*v*v*v*v)) + + +--- See Olver, p. 376-382. +BesselIAsymptOrder(v,vz) == + z := vz/v + Pi := PI +--- Use reflection formula (Atlas, p. 492) if v not in right half plane; Is this always accurate? + if REALPART(v)<0.0 + then return BesselIAsymptOrder(-v,vz) + 2.0/Pi*SIN(-v*Pi)*BesselKAsymptOrder(-v,vz) +--- Use the reflection formula (Atlas, p. 496) if z not in right half plane; + if REALPART(vz) < 0.0 + then return EXPT(-1.0,v)*BesselIAsymptOrder(v,-vz) + vinv := 1.0/v + opzsqroh := SQRT(1.0+z*z) + eta := opzsqroh + LOG(z/(1.0+opzsqroh)) + p := 1.0/opzsqroh + p2 := p*p + p4 := p2*p2 + u0p := 1. + u1p := 1.0/8.0*p-5.0/24.0*p*p2 + u2p := (9.0/128.0+(-77.0/192.0+385.0/1152.0*p2)*p2)*p2 + u3p := (75.0/1024.0+(-4563.0/5120.0+(17017.0/9216.0-85085.0/82944.0*p2)_ + *p2)*p2)*p2*p + u4p := (3675.0/32768.0+(-96833.0/40960.0+(144001.0/16384.0+(-7436429.0/663552.0+37182145.0/7962624.0*p2)*p2)*p2)*p2)*p4 + u5p := (59535.0/262144.0+(-67608983.0/9175040.0+(250881631.0/5898240.0+(-108313205.0/1179648.0+(5391411025.0/63700992.0-5391411025.0/191102976.0*p2)*p2)*p2)*p2)*p2)*p4*p + hornerresult := horner([u5p,u4p,u3p,u2p,u1p,u0p],vinv) + EXP(v*eta)/(SQRT(2.0*Pi*v)*SQRT(opzsqroh))*hornerresult + + +---See also Olver, pp. 376-382 +BesselKAsymptOrder (v,vz) == + z := vz/v + vinv := 1.0/v + opzsqroh := SQRT(1.0+z*z) + eta := opzsqroh + LOG(z/(1.0+opzsqroh)) + p := 1.0/opzsqroh + p2 := p**2 + p4 := p2**2 + u0p := 1. + u1p := (1.0/8.0*p-5.0/24.0*p**3)*(-1.0) + u2p := (9.0/128.0+(-77.0/192.0+385.0/1152.0*p2)*p2)*p2 + u3p := ((75.0/1024.0+(-4563.0/5120.0+(17017.0/9216.0-85085.0/82944.0*p2)_ + *p2)*p2)*p2*p)*(-1.0) + u4p := (3675.0/32768.0+(-96833.0/40960.0+(144001.0/16384.0+(-7436429.0/663552.0+37182145.0/7962624.0*p2)*p2)*p2)*p2)*p4 + u5p := ((59535.0/262144.0+(-67608983.0/9175040.0+(250881631.0/5898240.0+(-108313205.0/1179648.0+(5391411025.0/63700992.0-5391411025.0/191102976.0*p2)*p2)*p2)*p2)*p2)*p4*p)*(-1.0) + hornerresult := horner([u5p,u4p,u3p,u2p,u1p,u0p],vinv) + SQRT(PI/(2.0*v))*EXP(-v*eta)/(SQRT(opzsqroh))*hornerresult + + + + + + + +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document} |