1-- Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. 2-- All rights reserved. 3-- 4-- Redistribution and use in source and binary forms, with or without 5-- modification, are permitted provided that the following conditions are 6-- met: 7-- 8-- - Redistributions of source code must retain the above copyright 9-- notice, this list of conditions and the following disclaimer. 10-- 11-- - Redistributions in binary form must reproduce the above copyright 12-- notice, this list of conditions and the following disclaimer in 13-- the documentation and/or other materials provided with the 14-- distribution. 15-- 16-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the 17-- names of its contributors may be used to endorse or promote products 18-- derived from this software without specific prior written permission. 19-- 20-- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS 21-- IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED 22-- TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 23-- PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER 24-- OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 25-- EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 26-- PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 27-- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 28-- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 29-- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 30-- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 31 32 33-- Used to be SPECFNSF 34)package "BOOT" 35 36FloatError(formatstring,arg) == 37-- ERROR(formatstring,arg) 38 ERROR FORMAT([],formatstring,arg) 39 40float(x) == FLOAT(x, 0.0) 41 42fracpart(x) == 43 CADR(MULTIPLE_-VALUE_-LIST(FLOOR(x))) 44 45intpart(x) == 46 first(MULTIPLE_-VALUE_-LIST(FLOOR(x))) 47 48negintp(x) == 49 if ZEROP IMAGPART(x) and x<0.0 and ZEROP fracpart(x) 50 then 51 true 52 else 53 false 54 55-- Lisp PI is a long float and causes type errors, here we give 56-- enough digits to have double accuracy even after conversion 57-- to binary 58DEFCONSTANT(dfPi, 3.14159265358979323846264338328) 59 60--- Small float implementation of Gamma function. Valid for 61--- real arguments. Maximum error (relative) seems to be 62--- 2-4 ulps for real x 2<x<9, and up to ten ulps for larger x 63--- up to overflow. See Hart & Cheney. 64--- Bruce Char, April, 1990. 65 66horner(l,x) == 67 result := 0 68 for el in l repeat 69 result := result *x + el 70 return result 71 72r_gamma (x) == 73 if COMPLEXP(x) then FloatError('"Gamma not implemented for complex value ~D",x) 74 ZEROP (x-1.0) => 1.0 75 if x>20 then gammaStirling(x) else gammaRatapprox(x) 76 77r_lngamma (x) == 78 if x>20 then lnrgammaRatapprox(x) else LOG(gammaRatapprox(x)) 79 80cbeta(z,w) == 81 cgamma(z)*cgamma(w)/(cgamma(z+w)) 82 83gammaStirling(x) == 84 EXP(r_lngamma(x)) 85 86lnrgammaRatapprox(x) == 87 (x-.5)*LOG(x) - x + LOG(SQRT(2.0*dfPi)) + phiRatapprox(x) 88 89phiRatapprox(x) == 90 arg := 1/(x^2) 91 p := horner([.0666629070402007526,_ 92 .6450730291289920251389,_ 93 .670827838343321349614617,_ 94 .12398282342474941538685913],arg); 95 q := horner([1.0,7.996691123663644194772,_ 96 8.09952718948975574728214,_ 97 1.48779388109699298468156],arg); 98 result := p/(x*q) 99 result 100 101gammaRatapprox (x) == 102 if (x>=2 and x<=3) 103 then 104 result := gammaRatkernel(x) 105 else 106 if x>3 107 then 108 n := FLOOR(x)-2 109 a := x-n-2 110 reducedarg := 2+a 111 prod := */[reducedarg+i for i in 0..n-1] 112 result := prod* gammaRatapprox(reducedarg) 113 else 114 if (2>x and x>0) 115 then 116 n := 2-FLOOR(x) 117 a := x-FLOOR(x) 118 reducedarg := 2+a 119 prod := */[x+i for i in 0..n-1] 120 result := gammaRatapprox(reducedarg)/prod 121 else 122 Pi := dfPi 123 lx := MULTIPLE_-VALUE_-LIST(FLOOR(x)) 124 intpartx := first(lx)+1 125 restx := CADR(lx) 126 if ZEROP restx -- case of negative non-integer value 127 then 128 FloatError ('"Gamma undefined for non-positive integers: ~D",x) 129 else 130 result := Pi/(gammaRatapprox(1.0-x)*(-1.0)^(intpartx+1)*SIN(restx*Pi)) 131 result 132 133gammaRatkernel(x) == 134 p := horner(REVERSE([3786.01050348257245475108,_ 135 2077.45979389418732098416,_ 136 893.58180452374981423868,_ 137 222.1123961680117948396,_ 138 48.95434622790993805232,_ 139 6.12606745033608429879,_ 140 .778079585613300575867]),x-2) 141 q := horner(REVERSE([3786.01050348257187258861,_ 142 476.79386050368791516095,_ 143 -867.23098753110299445707,_ 144 83.55005866791976957459,_ 145 50.78847532889540973716,_ 146 -13.40041478578134826274,_ 147 1]),x-2.0) 148 p/q 149 150-- cgamma(z) Gamma function for complex arguments. 151-- Bruce Char April-May, 1990. 152-- 153-- Our text for complex gamma is H. Kuki's paper Complex Gamma 154-- Function with Error Control", CACM vol. 15, no. 4, ppp. 262-267. 155-- (April, 1972.) It uses the reflection formula and the basic 156-- z+1 recurrence to transform the argument into something that 157-- Stirling's asymptotic formula can handle. 158-- 159-- However along the way it does a few tricky things to reduce 160-- problems due to roundoff/cancellation error for particular values. 161 162-- cgammat is auxiliary "t" function (see p. 263 Kuki) 163cgammat(x) == 164 MAX(0.1, MIN(10.0, 10.0*SQRT(2.0) - ABS(x))) 165 166cgamma (z) == 167 z2 := IMAGPART(z) 168 z1 := REALPART(z) --- call real valued gamma if z is real 169 if ZEROP z2 170 then result := r_gamma(z1) 171 else 172 result := clngamma(z1,z2,z) 173 result := EXP(result) 174 result 175 176lncgamma(z) == 177 clngamma(REALPART z, IMAGPART z, z) 178 179clngamma(z1,z2,z) == 180 --- conjugate of gamma is gamma of conjugate. map 2nd and 4th quads 181 --- to first and third quadrants 182 if z1<0.0 183 then if z2 > 0.0 184 then result := CONJUGATE(clngammacase1(z1,-z2,COMPLEX(z1,-z2))) 185 else result := clngammacase1(z1,z2,z) 186 else if z2 < 0.0 187 then result := CONJUGATE(clngammacase23(z1,-z2,_ 188 COMPLEX(z1,-z2))) 189 else result := clngammacase23(z1,z2,z) 190 result 191 192clngammacase1(z1,z2,z) == 193 result1 := PiMinusLogSinPi(z1,z2,z) 194 result2 := clngamma(1.0-z1,-z2,1.0-z) 195 result1-result2 196 197PiMinusLogSinPi(z1,z2,z) == 198 cgammaG(z1,z2) - logH(z1,z2,z) 199 200cgammaG(z1,z2) == 201 LOG(2*dfPi) + dfPi*z2 - COMPLEX(0.0,1.0)*dfPi*(z1-.5) 202 203logH(z1,z2,z) == 204 z1bar := CADR(MULTIPLE_-VALUE_-LIST(FLOOR(z1))) ---frac part of z1 205 piz1bar := dfPi*z1bar 206 piz2 := dfPi*z2 207 twopiz2 := 2.0*piz2 208 i := COMPLEX(0.0,1.0) 209 part2 := EXP(twopiz2)*(2.0*SIN(piz1bar)^2 + SIN(2.0*piz1bar)*i) 210 part1 := -TANH(piz2)*(1.0+EXP(twopiz2)) 211--- part1 is another way of saying 1 - exp(2*Pi*z1bar) 212 LOG(part1+part2) 213 214clngammacase23(z1,z2,z) == 215 tz2 := cgammat(z2) 216 if (z1 < tz2) 217 then result:= clngammacase2(z1,z2,tz2,z) 218 else result:= clngammacase3(z) 219 result 220 221clngammacase2(z1,z2,tz2,z) == 222 n := float(CEILING(tz2-z1)) 223 zpn := z+n 224 (z-.5)*LOG(zpn) - (zpn) + cgammaBernsum(zpn) - cgammaAdjust(logS(z1,z2,z,n,zpn)) 225 226logS(z1,z2,z,n,zpn) == 227 sum := 0.0 228 for k in 0..(n-1) repeat 229 if z1+k < 5.0 - 0.6*z2 230 then sum := sum + LOG((z+k)/zpn) 231 else sum := sum + LOG(1.0 - (n-k)/zpn) 232 sum 233 234--- on p. 265, Kuki, logS result should have its imaginary part 235--- adjusted by 2 Pi if it is negative. 236cgammaAdjust(z) == 237 if IMAGPART(z)<0.0 238 then result := z + COMPLEX(0.0, 2.0*dfPi) 239 else result := z 240 result 241 242clngammacase3(z) == 243 (z- .5)*LOG(z) - z + cgammaBernsum(z) 244 245cgammaBernsum (z) == 246 sum := LOG(2.0*dfPi)/2.0 247 zterm := z 248 zsquaredinv := 1.0/(z*z) 249 l:= [.083333333333333333333, -.0027777777777777777778,_ 250 .00079365079365079365079, -.00059523809523809523810,_ 251 .00084175084175084175084, -.0019175269175269175269,_ 252 .0064102564102564102564] 253 for m in 1..7 for el in l repeat 254 zterm := zterm*zsquaredinv 255 sum := sum + el*zterm 256 sum 257 258 259 260 261--- nth derivatives of ln gamma for real x, n = 0,1,.... 262--- requires files floatutils, rgamma 263$PsiAsymptoticBern := VECTOR(0.0, 0.1666666666666667, -0.03333333333333333, 0.02380952380952381,_ 264 -0.03333333333333333, 0.07575757575757576, -0.2531135531135531, 1.166666666666667,_ 265 -7.092156862745098, 54.97117794486216, -529.1242424242424, 6192.123188405797,_ 266 -86580.25311355311, 1425517.166666667, -27298231.06781609, 601580873.9006424,_ 267 -15116315767.09216, 429614643061.1667, -13711655205088.33, 488332318973593.2,_ 268 -19296579341940070.0, 841693047573682600.0, -40338071854059460000.0) 269 270 271PsiAsymptotic(n,x) == 272 xn := x^n 273 xnp1 := xn*x 274 xsq := x*x 275 xterm := xsq 276 factterm := r_gamma(n+2)/2.0/r_gamma(float(n+1)) 277 --- initialize to 1/n! 278 sum := AREF($PsiAsymptoticBern,1)*factterm/xterm 279 for k in 2..22 repeat 280 xterm := xterm * xsq 281 if n=0 then factterm := 1.0/float(2*k) 282 else if n=1 then factterm := 1 283 else factterm := factterm * float(2*k+n-1)*float(2*k+n-2)/(float(2*k)*float(2*k-1)) 284 sum := sum + AREF($PsiAsymptoticBern,k)*factterm/xterm 285 PsiEps(n,x) + 1.0/(2.0*xnp1) + 1.0/xn * sum 286 287 288PsiEps(n,x) == 289 if n = 0 290 then 291 result := -LOG(x) 292 else 293 result := 1.0/(float(n)*(x^n)) 294 result 295 296PsiAsymptoticOrder(n,x,nterms) == 297 sum := 0 298 xterm := 1.0 299 np1 := n+1 300 for k in 0..nterms repeat 301 xterm := (x+float(k))^np1 302 sum := sum + 1.0/xterm 303 sum 304 305 306r_psi(n,x) == 307 if x<=0.0 308 then 309 if ZEROP fracpart(x) 310 then FloatError('"singularity encountered at ~D",x) 311 else 312 m := MOD(n,2) 313 sign := (-1)^m 314 if fracpart(x)=.5 315 then 316 skipit := 1 317 else 318 skipit := 0 319 sign*((dfPi^(n + 1))*cotdiffeval(n, dfPi*(-x), skipit) 320 + r_psi(n, 1.0 - x)) 321 else if n=0 322 then 323 - rPsiW(n,x) 324 else 325 r_gamma(float(n+1))*rPsiW(n,x)*(-1)^MOD(n+1,2) 326 327---Amos' w function, with w(0,x) picked to be -psi(x) for x>0 328rPsiW(n,x) == 329 if (x <=0 or n < 0) 330 then 331 HardError('"rPsiW not implemented for negative n or non-positive x") 332 nd := 6 -- magic number for number of digits in a word? 333 alpha := 3.5 + .40*nd 334 beta := 0.21 + (.008677e-3)*(nd-3) + (.0006038e-4)*(nd-3)^2 335 xmin := float(FLOOR(alpha + beta*n) + 1) 336 if n>0 337 then 338 a := MIN(0,1.0/float(n)*LOG(DOUBLE_-FLOAT_-EPSILON/MIN(1.0,x))) 339 c := EXP(a) 340 if ABS(a) >= 0.001 341 then 342 fln := x/c*(1.0-c) 343 else 344 fln := -x*a/c 345 bign := FLOOR(fln) + 1 346--- Amos says to use alternative series for large order if ordinary 347--- backwards recurrence too expensive 348 if (bign < 15) and (xmin > 7.0+x) 349 then 350 return PsiAsymptoticOrder(n,x,bign) 351 if x>= xmin 352 then 353 return PsiAsymptotic(n,x) 354---ordinary case -- use backwards recursion 355 PsiBack(n,x,xmin) 356 357PsiBack(n,x,xmin) == 358 xintpart := PsiIntpart(x) 359 x0 := x-xintpart ---frac part of x 360 result := PsiAsymptotic(n,x0+xmin+1.0) 361 for k in xmin..xintpart by -1 repeat 362--- Why not decrement from x? See Amos p. 498 363 result := result + 1.0/((x0 + float(k))^(n+1)) 364 result 365 366 367PsiIntpart(x) == 368 if x<0 369 then 370 result := -PsiInpart(-x) 371 else 372 result := FLOOR(x) 373 return result 374 375 376---Code for computation of derivatives of cot(z), necessary for 377--- polygamma reflection formula. If you want to compute n-th derivatives of 378---cot(Pi*x), you have to multiply the result of cotdiffeval by Pi^n. 379 380-- MCD: This is defined at the Lisp Level. 381-- COT(z) == 382-- 1.0/TAN(z) 383 384cotdiffeval(n,z,skipit) == 385---skip=1 if arg z is known to be an exact multiple of Pi/2 386 a := MAKE_-ARRAY(n+2) 387 SETF(AREF(a,0),0.0) 388 SETF(AREF(a,1),1.0) 389 for i in 2..n repeat 390 SETF(AREF(a,i),0.0) 391 for l in 1..n repeat 392 m := MOD(l+1,2) 393 for k in m..l+1 by 2 repeat 394 if k<1 395 then 396 t1 := 0 397 else 398 t1 := -AREF(a,k-1)*(k-1) 399 if k>l 400 then 401 t2 := 0 402 else 403 t2 := -AREF(a,k+1)*(k+1) 404 SETF(AREF(a,k), t1+t2) 405 --- evaluate d^N/dX^N cot(z) via Horner-like rule 406 v := COT(z) 407 sq := v^2 408 s := AREF(a,n+1) 409 for i in (n-1)..0 by -2 repeat 410 s := s*sq + AREF(a,i) 411 m := MOD(n,2) 412 if m=0 413 then 414 s := s*v 415 if skipit=1 416 then 417 if m=0 418 then 419 return 0 420 else 421 return AREF(a,0) 422 else 423 return s 424--- nth derivatives of ln gamma for complex z, n=0,1,... 425--- requires files rpsi (and dependents), floaterrors 426--- currently defined only in right half plane until reflection formula 427--- works 428 429--- B. Char, June, 1990. 430 431cPsi(n,z) == 432 x := REALPART(z) 433 y := IMAGPART(z) 434 if ZEROP y 435 then --- call real function if real 436 return r_psi(n, x) 437 if y<0.0 438 then -- if imagpart(z) negative, take conjugate of conjugate 439 conjresult := cPsi(n,COMPLEX(x,-y)) 440 return COMPLEX(REALPART(conjresult),-IMAGPART(conjresult)) 441 nterms := 22 442 bound := 10.0 443 if x<0.0 --- and ABS(z)>bound and ABS(y)<bound 444 then 445 FloatError('"Psi implementation can't compute at ~S ",[n,z]) 446--- return cpsireflect(n,x,y,z) 447 else if (x>0.0 and ABS(z)>bound ) --- or (x<0.0 and ABS(y)>bound) 448 then 449 return PsiXotic(n,PsiAsymptotic(n,z)) 450 else --- use recursion formula 451 m := CEILING(SQRT(bound*bound - y*y) - x + 1.0) 452 result := COMPLEX(0.0,0.0) 453 for k in 0..(m-1) repeat 454 result := result + 1.0/((z + float(k))^(n+1)) 455 return PsiXotic(n,result+PsiAsymptotic(n,z+m)) 456 457PsiXotic(n,result) == 458 r_gamma(float(n+1))*(-1)^MOD(n+1,2)*result 459 460--- c parameter to 0F1, possibly complex 461--- z argument to 0F1 462--- Depends on files floaterror, floatutils 463 464--- Program transcribed from Fortran,, p. 80 Luke 1977 465 466chebf01 (c,z) == 467--- w scale factor so that 0<z/w<1 468--- n n+2 coefficients will be produced stored in an array 469--- indexed from 0 to n+1. 470--- See Luke's books for further explanation 471 n := 75 --- ad hoc decision 472--- if ABS(z)/ABS(c) > 200.0 and ABS(z)>10000.0 473--- then 474--- FloatError('"cheb0F1 not implemented for ~S < 1",[c,z]) 475 w := 2.0*z 476--- arr will be used to store the Cheb. series coefficients 477 four:= 4.0 478 start := EXPT(10.0, -200) 479 n1 := n+1 480 n2 := n+2 481 a3 := 0.0 482 a2 := 0.0 483 a1 := start -- arbitrary starting value 484 z1 := four/w 485 ncount := n1 486 arr := MAKE_-ARRAY(n2) 487 SETF(AREF(arr,ncount) , start) -- start off 488 x1 := n2 489 c1 := 1.0 - c 490 for ncount in n..0 by -1 repeat 491 divfac := 1.0/x1 492 x1 := x1 -1.0 493 SETF(AREF(arr,ncount) ,_ 494 x1*((divfac+z1*(x1-c1))*a1 +_ 495 (1.0/x1 + z1*(x1+c1+1.0))*a2-divfac*a3)) 496 a3 := a2 497 a2 := a1 498 a1 := AREF(arr,ncount) 499 SETF(AREF(arr,0),AREF(arr,0)/2.0) 500-- compute scale factor 501 rho := AREF(arr,0) 502 sum := rho 503 p := 1.0 504 for i in 1..n1 repeat 505 rho := rho - p*AREF(arr,i) 506 sum := sum+AREF(arr,i) 507 p := -p 508 for l in 0..n1 repeat 509 SETF(AREF(arr,l), AREF(arr,l)/rho) 510 sum := sum/rho 511--- Now evaluate array at argument 512 b := 0.0 513 temp := 0.0 514 for i in (n+1)..0 by -1 repeat 515 cc := b 516 b := temp 517 temp := -cc + AREF(arr,i) 518 temp 519 520 521--- c parameter to 0F1 522--- w scale factor so that 0<z/w<1 523--- n n+2 coefficients will be produced stored in an array 524--- indexed from 0 to n+1. 525--- See Luke's books for further explanation 526 527--- Program transcribed from Fortran,, p. 80 Luke 1977 528chebf01coefmake (c,w,n) == 529--- arr will be used to store the Cheb. series coefficients 530 four:= 4.0 531 start := EXPT(10.0, -200) 532 n1 := n+1 533 n2 := n+2 534 a3 := 0.0 535 a2 := 0.0 536 a1 := start -- arbitrary starting value 537 z1 := four/w 538 ncount := n1 539 arr := MAKE_-ARRAY(n2) 540 SETF(AREF(arr,ncount) , start) -- start off 541 x1 := n2 542 c1 := 1.0 - c 543 for ncount in n..0 by -1 repeat 544 divfac := 1.0/x1 545 x1 := x1 -1.0 546 SETF(AREF(arr,ncount) ,_ 547 x1*((divfac+z1*(x1-c1))*a1 +_ 548 (1.0/x1 + z1*(x1+c1+1.0))*a2-divfac*a3)) 549 a3 := a2 550 a2 := a1 551 a1 := AREF(arr,ncount) 552 SETF(AREF(arr,0),AREF(arr,0)/2.0) 553-- compute scale factor 554 rho := AREF(arr,0) 555 sum := rho 556 p := 1.0 557 for i in 1..n1 repeat 558 rho := rho - p*AREF(arr,i) 559 sum := sum+AREF(arr,i) 560 p := -p 561 for l in 0..n1 repeat 562 SETF(AREF(arr,l), AREF(arr,l)/rho) 563 sum := sum/rho 564 return([sum,arr]) 565 566 567 568 569---evaluation of Chebychev series of degree n at x, where the series's 570---coefficients are given by the list in descending order (coef. of highest 571---power first) 572 573---May be numerically unstable for certain lists of coefficients; 574--- could possibly reverse sequence of coefficients 575 576--- Cheney and Hart p. 15. 577 578--- B. Char, March 1990 579 580--- If plist is a list of coefficients for the Chebychev approximation 581--- of a function f(x), then chebderiveval computes the Chebychev approximation 582--- of f'(x). See Luke, "Special Functions and their approximations, vol. 1 583--- Academic Press 1969., p. 329 (from Clenshaw and Cooper) 584 585--- < definition to be supplied> 586 587--- chebstareval(plist,n) computes a Chebychev approximation from a 588--- coefficient list, using shifted Chebychev polynomials of the first kind 589--- The defining relation is that T*(n,x) = T(n,2*x-1). Thus the interval 590--- [0,1] of T*n is the interval [-1,1] of Tn. 591 592chebstarevalarr(coefarr,x,n) == -- evaluation of sum(C(n)*T*(n,x)) 593 594 b := 0 595 temp := 0 596 y := 2*(2*x-1) 597 598 for i in (n+1)..0 by -1 repeat 599 c := b 600 b := temp 601 temp := y*b -c + AREF(coefarr,i) 602 temp - y*b/2 603 604--Float definitions for Bessel functions I and J. 605--External references: cgamma, rgamma, chebf01coefmake, chebevalstarsf 606-- floatutils 607 608---BesselJ works for complex and real values of v and z 609BesselJ(v,z) == 610---Ad hoc boundaries for approximation 611 B1:= 10 612 B2:= 10 613 n := 50 --- number of terms in Chebychev series. 614 --- tests for negative integer order 615 (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) => 616 --- odd or even according to v (9.1.5 A&S) 617 --- $J_{-n}(z)=(-1)^{n} J_{n}(z)$ 618 BesselJ(-v,z)*EXPT(-1.0,v) 619 (FLOATP(z) and (z<0)) or (COMPLEXP(z) and REALPART(z)<0.0) => 620 --- negative argument (9.1.35 A&S) 621 --- $J_{\nu}(z e^{m \pi i}) = e^{m \nu \pi i} J_{\nu}(z)$ 622 BesselJ(v,-z)*EXPT(-1.0,v) 623 ZEROP z and ((FLOATP(v) and (v>=0.0)) or (COMPLEXP(v) and 624 ZEROP IMAGPART(v) and REALPART(v)>=0.0)) => --- zero arg, pos. real order 625 ZEROP v => 1.0 --- J(0,0)=1 626 0.0 --- J(v,0)=0 for real v>0 627 rv := ABS(v) 628 rz := ABS(z) 629 (rz>B1) and (rz > B2*rv) => --- asymptotic argument 630 BesselJAsympt(v,z) 631 (rv>B1) and (rv > B2*rz) => --- asymptotic order 632 BesselJAsymptOrder(v,z) 633 (rz< B1) and (rv<B1) => --- small order and argument 634 arg := -(z*z)/4.0 635 w := 2.0*arg 636 vp1 := v+1.0 637 [sum,arr] := chebf01coefmake(vp1,w,n) 638 ---if we get NaNs then half n 639 while not _=(sum,sum) repeat 640 n:=FLOOR(n/2) 641 [sum,arr] := chebf01coefmake(vp1,w,n) 642 ---now n is safe, can we increase it (we know that 2*n is bad)? 643 chebstarevalarr(arr,arg/w,n)/cgamma(vp1)*EXPT(z/2.0,v) 644 true => BesselJRecur(v,z) 645 FloatError('"BesselJ not implemented for ~S", [v,z]) 646 647BesselJRecur(v,z) == 648 -- boost order 649 --Numerical.Recipes. suggest so:=v+sqrt(n.s.f.^2*v) 650 so:=15.0*z 651 -- reduce order until non-zero 652 while ZEROP ABS(BesselJAsymptOrder(so,z)) repeat so:=so/2.0 653 if ABS(so)<ABS(z) then so:=v+18.*SQRT(v) 654 m:= FLOOR(ABS(so-v))+1 655 w:=MAKE_-ARRAY(m) 656 SETF(AREF(w,m-1),BesselJAsymptOrder(v+m-1,z)) 657 SETF(AREF(w,m-2),BesselJAsymptOrder(v+m-2,z)) 658 for i in m-3 .. 0 by -1 repeat 659 SETF(AREF(w,i), 2.0 * (v+i+1.0) * AREF(w,i+1) /z -AREF(w,i+2)) 660 AREF(w,0) 661 662BesselI(v,z) == 663 B1 := 15.0 664 B2 := 10.0 665 ZEROP(z) and FLOATP(v) and (v>=0.0) => --- zero arg, pos. real order 666 ZEROP(v) => 1.0 --- I(0,0)=1 667 0.0 --- I(v,0)=0 for real v>0 668--- Transformations for negative integer orders 669 FLOATP(v) and ZEROP(fracpart(v)) and (v<0) => BesselI(-v,z) 670--- Halfplane transformations for Re(z)<0 671 REALPART(z)<0.0 => BesselI(v,-z)*EXPT(-1.0,v) 672--- Conjugation for complex order and real argument 673 REALPART(v)<0.0 and not ZEROP IMAGPART(v) and FLOATP(z) => 674 CONJUGATE(BesselI(CONJUGATE(v),z)) 675---We now know that Re(z)>= 0.0 676 ABS(z) > B1 => --- asymptotic argument case 677 FloatError('"BesselI not implemented for ~S",[v,z]) 678 ABS(v) > B1 => 679 FloatError('"BesselI not implemented for ~S",[v,z]) 680--- case of small argument and order 681 REALPART(v)>= 0.0 => besselIback(v,z) 682 REALPART(v)< 0.0 => 683 chebterms := 50 684 besselIcheb(z,v,chebterms) 685 FloatError('"BesselI not implemented for ~S",[v,z]) 686 687--- Compute n terms of the chebychev series for f01 688besselIcheb(z,v,n) == 689 arg := (z*z)/4.0 690 w := 2.0*arg; 691 vp1 := v+1.0; 692 [sum,arr] := chebf01coefmake(vp1,w,n) 693 result := chebstarevalarr(arr,arg/w,n)/cgamma(vp1)*EXPT(z/2.0,v) 694 695besselIback(v,z) == 696 ipv := IMAGPART(v) 697 rpv := REALPART(v) 698 lm := MULTIPLE_-VALUE_-LIST(FLOOR(rpv)) 699 m := first(lm) --- floor of real part of v 700 n := 2*MAX(20,m+10) --- how large the back recurrence should be 701 tv := CADR(lm)+(v-rpv) --- fractional part of real part of v 702 --- plus imaginary part of v 703 vp1 := tv+1.0; 704 result := BesselIBackRecur(v,m,tv,z,'"I",n) 705 result := result/cgamma(vp1)*EXPT(z/2.0,tv) 706 707--- Backward recurrence for Bessel functions. Luke (1975), p. 247. 708--- works for -Pi< arg z <= Pi and -Pi < arg v <= Pi 709BesselIBackRecur(largev,argm,v,z,type,n) == 710--- v + m = largev 711 one := 1.0 712 two := 2.0 713 zero := 0.0 714 start := EXPT(10.0,-40) 715 z2 := two/z 716 m2 := n+3 717 w:=MAKE_-ARRAY(m2+1) 718 SETF(AREF(w,m2), zero) --- start off 719 if type = '"I" 720 then 721 val := one 722 else 723 val := -one 724 m1 := n+2 725 SETF(AREF(w,m1), start) 726 m := n+1 727 xm := float(m) 728 ct1 := z2*(xm+v) 729 --- initialize 730 for m in (n+1)..1 by -1 repeat 731 SETF(AREF(w,m), AREF(w,m+1)*ct1 + val*AREF(w,m+2)) 732 ct1 := ct1 - z2 733 m := 1 + FLOOR(n/2) 734 m2 := m + m -1 735 if (v=0) 736 then 737 pn := AREF(w, m2 + 2) 738 for m2 in (2*m-1)..3 by -2 repeat 739 pn := AREF(w, m2) - val *pn 740 pn := AREF(w,1) - val*(pn+pn) 741 else 742 v1 := v-one 743 xm := float(m) 744 ct1 := v + xm + xm 745 pn := ct1*AREF(w, m2 + 2) 746 for m2 in (m+m -1)..3 by -2 repeat 747 ct1 := ct1 - two 748 pn := ct1*AREF(w,m2) - val*pn/xm*(v1+xm) 749 xm := xm - one 750 pn := AREF(w,1) - val * pn 751 m1 := n+2 752 for m in 1..m1 repeat 753 SETF(AREF(w,m), AREF(w,m)/pn) 754 AREF(w,argm+1) 755 756 757 758 759---Asymptotic functions for large values of z. See p. 204 Luke 1969 vol. 1. 760 761--- mu is 4*v^2 762--- zsqr is z^2 763--- zfth is z^4 764 765BesselasymptA(mu,zsqr,zfth) == 766 (mu -1)/(16.0*zsqr) * (1 + (mu - 13.0)/(8.0*zsqr) + _ 767 (mu^2 - 53.0*mu + 412.0)/(48.0*zfth)) 768 769BesselasymptB(mu,z,zsqr,zfth) == 770 musqr := mu*mu 771 z + (mu-1.0)/(8.0*z) *(1.0 + (mu - 25.0)/(48.0*zsqr) + _ 772 (musqr - 114.0*mu + 1073.0)/(640.0*zfth) +_ 773 (5.0*mu*musqr - 1535.0*musqr + 54703.0*mu - 375733.0)/(128.0*zsqr*zfth)) 774 775--- Asymptotic series only works when |v| < |z|. 776 777BesselJAsympt (v,z) == 778 pi := dfPi 779 mu := 4.0*v*v 780 zsqr := z*z 781 zfth := zsqr*zsqr 782 SQRT(2.0/(pi*z))*EXP(BesselasymptA(mu,zsqr,zfth))*_ 783 COS(BesselasymptB(mu,z,zsqr,zfth) - pi*v/2.0 - pi/4.0) 784 785 786---Asymptotic series for I. See Whittaker, p. 373. 787--- valid for -3/2 Pi < arg z < 1/2 Pi 788 789BesselIAsympt(v,z,n) == 790 i := COMPLEX(0.0, 1.0) 791 if (REALPART(z) = 0.0) 792 then return EXPT(i,v)*BesselJ(v,-IMAGPART(z)) 793 sum1 := 0.0 794 sum2 := 0.0 795 fourvsq := 4.0*v^2 796 two := 2.0 797 eight := 8.0 798 term1 := 1.0 799--- sum1, sum2, fourvsq,two,i,eight,term1]) 800 for r in 1..n repeat 801 term1 := -term1 *(fourvsq-(two*float(r)-1.0)^2)/_ 802 (float(r)*eight*z) 803 sum1 := sum1 + term1 804 sum2 := sum2 + ABS(term1) 805 sqrttwopiz := SQRT(two*dfPi*z) 806 EXP(z)/sqrttwopiz*(1.0 + sum1 ) +_ 807 EXP(-(float(n)+.5)*dfPi*i)*EXP(-z)/sqrttwopiz*(1.0+ sum2) 808 809 810---Asymptotic formula for BesselJ when order is large comes from 811---Debye (1909). See Olver, Asymptotics and Special Functions, p. 134. 812---Expansion good for 0<=phase(v)<Pi 813---A&S recommend "uniform expansion" with complicated coefficients and Airy function. 814---Debye's Formula is in 9.3.7,9.3.9,9.3.10 of A&S 815---FriCAS recurrence for u_{k} 816---f(0)==1::EXPR INT 817---f(n)== (t^2)*(1-t^2)*D(f(n-1),t)/2 + (1/8)*integrate( (1-5*t^2)*f(n-1),t) 818BesselJAsymptOrder(v,z) == 819 sechalpha := z/v 820 alpha := ACOSH(1.0/sechalpha) 821 tanhalpha := SQRT(1.0-(sechalpha*sechalpha)) 822 -- cothalpha := 1.0/tanhalpha 823 ca := 1.0/tanhalpha 824 825 Pi := dfPi 826 ca2:=ca*ca 827 ca4:=ca2*ca2 828 ca8:=ca4*ca4 829 EXP(-v*(alpha-tanhalpha))/SQRT(2.0*Pi*v*tanhalpha)*_ 830 (1.0+_ 831 horner([ -5.0, 3.0],_ 832 ca2)*ca/(v*24.0)+_ 833 horner([ 385.0, -462.0, 81.0],_ 834 ca2)*ca2/(1152.0*v*v)+_ 835 horner([ -425425.0, 765765.0, -369603.0, 30375.0],_ 836 ca2)*ca2*ca/(414720.0*v*v*v)+_ 837 horner([ 185910725.0, -446185740.0, 349922430.0, -94121676.0, 4465125.0],_ 838 ca2)*ca4/(39813120.0*v*v*v*v)+_ 839 horner([ -188699385875.0, 566098157625.0, -614135872350.0, 284499769554.0, -49286948607.0, 1519035525.0],_ 840 ca2)*ca4*ca/(6688604160.0*v*v*v*v*v)+_ 841 horner([1023694168371875.0,-3685299006138750.0,5104696716244125.0,-3369032068261860.0,1050760774457901.0,-127577298354750.0,2757049477875.0],_ 842 ca2)*ca4*ca2/(4815794995200.0*v*v*v*v*v*v)) 843 844 845--- See Olver, p. 376-382. 846BesselIAsymptOrder(v,vz) == 847 z := vz/v 848 Pi := dfPi 849--- Use reflection formula (Atlas, p. 492) if v not in right half plane; Is this always accurate? 850 if REALPART(v)<0.0 851 then return BesselIAsymptOrder(-v,vz) + 2.0/Pi*SIN(-v*Pi)*BesselKAsymptOrder(-v,vz) 852--- Use the reflection formula (Atlas, p. 496) if z not in right half plane; 853 if REALPART(vz) < 0.0 854 then return EXPT(-1.0,v)*BesselIAsymptOrder(v,-vz) 855 vinv := 1.0/v 856 opzsqroh := SQRT(1.0+z*z) 857 eta := opzsqroh + LOG(z/(1.0+opzsqroh)) 858 p := 1.0/opzsqroh 859 p2 := p*p 860 p4 := p2*p2 861 u0p := 1. 862 u1p := 1.0/8.0*p-5.0/24.0*p*p2 863 u2p := (9.0/128.0+(-77.0/192.0+385.0/1152.0*p2)*p2)*p2 864 u3p := (75.0/1024.0+(-4563.0/5120.0+(17017.0/9216.0-85085.0/82944.0*p2)_ 865 *p2)*p2)*p2*p 866 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 867 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 868 hornerresult := horner([u5p,u4p,u3p,u2p,u1p,u0p],vinv) 869 EXP(v*eta)/(SQRT(2.0*Pi*v)*SQRT(opzsqroh))*hornerresult 870 871 872---See also Olver, pp. 376-382 873BesselKAsymptOrder (v,vz) == 874 z := vz/v 875 vinv := 1.0/v 876 opzsqroh := SQRT(1.0+z*z) 877 eta := opzsqroh + LOG(z/(1.0+opzsqroh)) 878 p := 1.0/opzsqroh 879 p2 := p^2 880 p4 := p2^2 881 u0p := 1. 882 u1p := (1.0/8.0*p-5.0/24.0*p^3)*(-1.0) 883 u2p := (9.0/128.0+(-77.0/192.0+385.0/1152.0*p2)*p2)*p2 884 u3p := ((75.0/1024.0+(-4563.0/5120.0+(17017.0/9216.0-85085.0/82944.0*p2)_ 885 *p2)*p2)*p2*p)*(-1.0) 886 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 887 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) 888 hornerresult := horner([u5p,u4p,u3p,u2p,u1p,u0p],vinv) 889 SQRT(dfPi/(2.0*v))*EXP(-v*eta)/(SQRT(opzsqroh))*hornerresult 890 891 892-- Conversion between spad and lisp complex representations 893s_to_c(c) == COMPLEX(first c, CDR c) 894c_to_s(c) == CONS(REALPART c, IMAGPART c) 895c_to_r(c) == 896 r := REALPART c 897 i := IMAGPART c 898 if ZEROP(i) or (ABS(i) < 1.0E-10*(ABS r)) then 899 r 900 else 901 error "Result is not real." 902 903c_to_rf(c) == COERCE(c_to_r(c), 'DOUBLE_-FLOAT) 904 905-- Wrappers for functions in the special function package 906c_lngamma(z) == c_to_s(lncgamma(s_to_c z)) 907 908c_gamma(z) == c_to_s(cgamma (s_to_c z)) 909 910c_psi(n, z) == c_to_s(cPsi(n, s_to_c(z))) 911 912r_besselj(n, x) == c_to_r(BesselJ(n, x)) 913c_besselj(v, z) == c_to_s(BesselJ(s_to_c(v), s_to_c(z))) 914 915r_besseli(n, x) == c_to_r(BesselI(n, x)) 916c_besseli(v, z) == c_to_s(BesselI(s_to_c(v), s_to_c(z))) 917 918c_hyper0f1(a, z) == c_to_s(chebf01(s_to_c(a), s_to_c(z))) 919