1 /** 2 * Mathlib : A C Library of Special Functions 3 * Copyright (C) 1998-2014 Ross Ihaka and the R Core team. 4 * Copyright (C) 2002-3 The R Foundation 5 * Copyright (C) 2014-3 Ilya Yaroshenko (Ported to D) 6 * 7 * This program is free software; you can redistribute it and/or modify 8 * it under the terms of the GNU General Public License as published b 9 * the Free Software Foundation; either version 2 of the License, or 10 * (at your option) any later version. 11 * 12 * This program is distributed in the hope that it will be useful, 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of 14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 * GNU General Public License for more details. 16 * 17 * You should have received a copy of the GNU General Public License 18 * along with this program; if not, a copy is available at 19 * http://www.r-project.org/Licenses/ 20 */ 21 module bessel; 22 23 static import std.math; 24 import core.stdc.stdlib; 25 import core.stdc.tgmath; 26 import std.algorithm; 27 import std.exception; 28 import std.string; 29 import std.traits; 30 import std.typecons; 31 32 public import std.typecons : Flag; 33 34 /* Documentation that apply to several of the bessel[IJKY] functions */ 35 36 /* ******************************************************************* 37 38 Explanation of machine-dependent constants 39 40 beta = Radix for the floating-point system 41 minexp = Smallest representable power of beta 42 maxexp = Smallest power of beta that overflows 43 it = p = Number of bits (base-beta digits) in the mantissa 44 (significand) of a working precision (floating-point) variable 45 NSIG = Decimal significance desired. Should be set to 46 INT(LOG10(2)*it+1). Setting NSIG lower will result 47 in decreased accuracy while setting NSIG higher will 48 increase CPU time without increasing accuracy. The 49 truncation error is limited to a relative error of 50 T=.5*10^(-NSIG). 51 ENTEN = 10 ^ K, where K is the largest int such that 52 ENTEN is machine-representable in working precision 53 ENSIG = 10 ^ NSIG 54 RTNSIG = 10 ^ (-K) for the smallest int K such that K >= NSIG/4 55 ENMTEN = Smallest ABS(X) such that X/4 does not underflow 56 XINF = Largest positive machine number; approximately beta ^ maxexp 57 == double.max (defined in #include <float.h>) 58 SQXMIN = Square root of beta ^ minexp = sqrt(double.min_normal) 59 60 EPS = The smallest positive floating-point number such that 1.0+EPS > 1.0 61 = beta ^ (-p) == double.epsilon 62 63 64 For I : 65 66 EXPARG = Largest working precision argument that the liaary 67 EXP routine can handle and upper limit on the 68 magnitude of X when IZE=1; approximately LOG(beta ^ maxexp) 69 70 For I and J : 71 72 xlrg_IJ = xlrg_BESS_IJ (was = XLARGE). Upper limit on the magnitude of X 73 (when IZE=2 for I()). Bear in mind that if floor(abs(x)) =: N, then 74 at least N iterations of the backward recursion will be executed. 75 The value of 10 ^ 4 was used till Feb.2009, when it was increased 76 to 10 ^ 5 (= 1e5). 77 78 For j : 79 XMIN_J = Smallest acceptable argument for RBESY; approximately 80 max(2*beta ^ minexp, 2/XINF), rounded up 81 82 For Y : 83 84 xlrg_Y = (was = XLARGE). Upper bound on X; 85 approximately 1/DEL, because the sine and cosine functions 86 have lost about half of their precision at that point. 87 88 EPS_SINC = Machine number below which sin(x)/x = 1; approximately SQRT(EPS). 89 THRESH = Lower bound for use of the asymptotic form; 90 approximately AINT(-LOG10(EPS/2.0))+1.0 91 92 93 For K : 94 95 xmax_k = (was = XMAX). Upper limit on the magnitude of X when ize = 1; 96 i.e. maximal x for UNscaled answer. 97 98 Solution to equation: 99 W(X) * (1 -1/8 X + 9/128 X^2) = beta ^ minexp 100 where W(X) = EXP(-X)*SQRT(PI/2X) 101 102 -------------------------------------------------------------------- 103 104 Approximate values for some important machines are: 105 106 beta minexp maxexp it NSIG ENTEN ENSIG RTNSIG ENMTEN EXPARG 107 IEEE (IBM/XT, 108 SUN, etc.) (S.P.) 2 -126 128 24 8 1e38 1e8 1e-2 4.70e-38 88 109 IEEE (...) (D.P.) 2 -1022 1024 53 16 1e308 1e16 1e-4 8.90e-308 709 110 CRAY-1 (S.P.) 2 -8193 8191 48 15 1e2465 1e15 1e-4 1.84e-2466 5677 111 Cyber 180/855 112 under NOS (S.P.) 2 -975 1070 48 15 1e322 1e15 1e-4 1.25e-293 741 113 IBM 3033 (D.P.) 16 -65 63 14 5 1e75 1e5 1e-2 2.16e-78 174 114 VAX (S.P.) 2 -128 127 24 8 1e38 1e8 1e-2 1.17e-38 88 115 VAX D-Format (D.P.) 2 -128 127 56 17 1e38 1e17 1e-5 1.17e-38 88 116 VAX G-Format (D.P.) 2 -1024 1023 53 16 1e307 1e16 1e-4 2.22e-308 709 117 118 119 And routine specific : 120 121 xlrg_IJ xlrg_Y xmax_k EPS_SINC XMIN_J XINF THRESH 122 IEEE (IBM/XT, 123 SUN, etc.) (S.P.) 1e4 1e4 85.337 1e-4 2.36e-38 3.40e38 8. 124 IEEE (...) (D.P.) 1e4 1e8 705.342 1e-8 4.46e-308 1.79e308 16. 125 CRAY-1 (S.P.) 1e4 2e7 5674.858 5e-8 3.67e-2466 5.45e2465 15. 126 Cyber 180/855 127 under NOS (S.P.) 1e4 2e7 672.788 5e-8 6.28e-294 1.26e322 15. 128 IBM 3033 (D.P.) 1e4 1e8 177.852 1e-8 2.77e-76 7.23e75 17. 129 VAX (S.P.) 1e4 1e4 86.715 1e-4 1.18e-38 1.70e38 8. 130 VAX e-Format (D.P.) 1e4 1e9 86.715 1e-9 1.18e-38 1.70e38 17. 131 VAX G-Format (D.P.) 1e4 1e8 706.728 1e-8 2.23e-308 8.98e307 16. 132 133 */ 134 private enum double nsig_BESS = 16; 135 private enum double ensig_BESS = 1e16; 136 private enum double rtnsig_BESS = 1e-4; 137 private enum double enmten_BESS = 8.9e-308; 138 private enum double enten_BESS = 1e308; 139 private enum double exparg_BESS = 709.; 140 private enum double xlrg_BESS_IJ = 1e5; 141 private enum double xlrg_BESS_Y = 1e8; 142 private enum double thresh_BESS_Y = 16.; 143 private enum double xmax_BESS_K = 705.342 /* maximal x for UNscaled answer */ ; 144 /* sqrt(double.min_normal) = 1.491668e-154 */ 145 private enum double sqxmin_BESS_K = 1.49e-154; 146 /* x < eps_sinc <==> sin(x)/x == 1 (particularly "==>"); 147 Linux (around 2001-02) gives 2.14946906753213e-08 148 Solaris 2.5.1 gives 2.14911933289084e-08 149 */ 150 private enum double M_eps_sinc = 2.149e-8; 151 private enum double M_SQRT_2dPI = 0.797884560802865355879892119869; 152 /// 153 double besselK(double x, double alpha, 154 Flag!"ExponentiallyScaled" expFlag = Flag!"ExponentiallyScaled".no) 155 { 156 ptrdiff_t nb; 157 ptrdiff_t ncalc; 158 double* b; 159 /* NaNs propagated correctly */ 160 if (std.math.isNaN(x) || std.math.isNaN(alpha)) 161 return x + alpha; 162 if (x < 0) 163 return double.nan; 164 if (alpha < 0) 165 alpha = -alpha; 166 nb = 1 + cast(int) floor(alpha); /* b.length-1 <= |alpha| < b.length */ 167 alpha -= nb - 1; 168 b = cast(double* ) calloc(nb, double.sizeof).enforce("besselK allocation error"); 169 scope(exit) b.free; 170 besselK(x, alpha, expFlag, b[0 .. nb], ncalc); 171 if (ncalc != nb) 172 { 173 /* error input */ 174 if (ncalc < 0) 175 throw new Exception(format("besselK(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n", 176 x, ncalc, nb, alpha)); 177 else throw new Exception(format( 178 "besselK(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1))); 179 } 180 x = b[nb - 1]; 181 return x; 182 } 183 184 /// 185 unittest 186 { 187 assert(std.math.approxEqual(besselK(2, 1), 0.139865881816522427284)); 188 assert(std.math.approxEqual(besselK(20, 1), 5.88305796955703817765028e-10)); 189 assert(std.math.approxEqual(besselK(1, 1.2), 0.701080)); 190 } 191 192 /* modified version of besselK that accepts a work array instead of 193 allocating one. */ 194 double besselK(double x, double alpha, Flag!"ExponentiallyScaled" expFlag, double[] b) 195 { 196 ptrdiff_t nb; 197 ptrdiff_t ncalc; 198 /* NaNs propagated correctly */ 199 if (std.math.isNaN(x) || std.math.isNaN(alpha)) 200 return x + alpha; 201 if (x < 0) 202 return double.nan; 203 if (alpha < 0) 204 alpha = -alpha; 205 nb = 1 + cast(int) floor(alpha); /* b.length-1 <= |alpha| < b.length */ 206 alpha -= nb - 1; 207 besselK(x, alpha, expFlag, b[0 .. nb], ncalc); 208 if (ncalc != nb) 209 { 210 /* error input */ 211 if (ncalc < 0) 212 throw new Exception(format("besselK(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n", 213 x, ncalc, nb, alpha)); 214 else throw new Exception(format( 215 "besselK(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1))); 216 } 217 x = b[nb - 1]; 218 return x; 219 } 220 221 void besselK(double x, double alpha, Flag!"ExponentiallyScaled" expFlag, double[] b, 222 ref ptrdiff_t ncalc) 223 { 224 /*------------------------------------------------------------------- 225 226 This routine calculates modified Bessel functions 227 of the third kind, K_(N+ALPHA) (X), for non-negative 228 argument X, and non-negative order N+ALPHA, with or without 229 exponential scaling. 230 231 Explanation of variables in the calling sequence 232 233 X - Non-negative argument for which 234 K's or exponentially scaled K's (K*EXP(X)) 235 are to be calculated. If K's are to be calculated, 236 X must not be greater than XMAX_BESS_K. 237 ALPHA - Fractional part of order for which 238 K's or exponentially scaled K's (K*EXP(X)) are 239 to be calculated. 0 <= ALPHA < 1.0. 240 NB - Number of functions to be calculated, NB > 0. 241 The first function calculated is of order ALPHA, and the 242 last is of order (NB - 1 + ALPHA). 243 expFlag - Type. expFlag = No if unscaled K's are to be calculated, 244 = Yes if exponentially scaled K's are to be calculated. 245 BK - Output vector of length NB. If the 246 routine terminates normally (NCALC=NB), the vector BK 247 contains the functions K(ALPHA,X), ... , K(NB-1+ALPHA,X), 248 or the corresponding exponentially scaled functions. 249 If (0 < NCALC < NB), BK(I) contains correct function 250 values for I <= NCALC, and contains the ratios 251 K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array. 252 NCALC - Output variable indicating possible errors. 253 Before using the vector BK, the user should check that 254 NCALC=NB, i.e., all orders have been calculated to 255 the desired accuracy. See error returns below. 256 257 258 ******************************************************************* 259 260 Error returns 261 262 In case of an error, NCALC != NB, and not all K's are 263 calculated to the desired accuracy. 264 265 NCALC < -1: An argument is out of range. For example, 266 NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >= XMAX_BESS_K. 267 In this case, the B-vector is not calculated, 268 and NCALC is set to MIN0(NB,0)-2 so that NCALC != NB. 269 NCALC = -1: Either K(ALPHA,X) >= XINF or 270 K(ALPHA+NB-1,X)/K(ALPHA+NB-2,X) >= XINF. In this case, 271 the B-vector is not calculated. Note that again 272 NCALC != NB. 273 274 0 < NCALC < NB: Not all requested function values could 275 be calculated accurately. BK(I) contains correct function 276 values for I <= NCALC, and contains the ratios 277 K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array. 278 279 280 Intrinsic functions required are: 281 282 ABS, AINT, EXP, INT, LOG, MAX, MIN, SINH, SQRT 283 284 285 Acknowledgement 286 287 This program is based on a program written b J. B. Campbell 288 (2) that computes values of the Bessel functions K of float 289 argument and float order. Modifications include the addition 290 of non-scaled functions, parameterization of machine 291 dependencies, and the use of more accurate approximations 292 for SINH and SIN. 293 294 References: "On Temme's Algorithm for the Modified Bessel 295 Functions of the Third Kind," Campbell, J. B., 296 TOMS 6(4), Dec. 1980, pp. 581-586. 297 298 "A FORTRAN IV Suaoutine for the Modified Bessel 299 Functions of the Third Kind of Real Order and Real 300 Argument," Campbell, J. B., Report NRC/ERB-925, 301 National Research Council, Canada. 302 303 Latest modification: May 30, 1989 304 305 Modified b: W. J. Cody and L. Stoltz 306 Applied Mathematics Division 307 Argonne National Laboratory 308 Argonne, IL 60439 309 310 ------------------------------------------------------------------- 311 */ 312 313 /*--------------------------------------------------------------------- 314 * Mathematical constants 315 * A = LOG(2) - Euler's constant 316 * D = SQRT(2/PI) 317 ---------------------------------------------------------------------*/ 318 static immutable double a = .11593151565841244881; 319 /*--------------------------------------------------------------------- 320 P, Q - Approximation for LOG(GAMMA(1+ALPHA))/ALPHA + Euler's constant 321 Coefficients converted from hex to decimal and modified 322 b W. J. Cody, 2/26/82 */ 323 static immutable double[8] p = [.805629875690432845, 20.4045500205365151, 324 157.705605106676174, 536.671116469207504, 900.382759291288778, 325 730.923886650660393, 229.299301509425145, .822467033424113231]; 326 static immutable double[7] q = [29.4601986247850434, 277.577868510221208, 327 1206.70325591027438, 2762.91444159791519, 3443.74050506564618, 328 2210.63190113378647, 572.267338359892221]; 329 /* R, S - Approximation for (1-ALPHA*PI/SIN(ALPHA*PI))/(2.D0*ALPHA) */ 330 static immutable double[5] r = [-.48672575865218401848, 331 13.079485869097804016, 332 -101.96490580880537526, 347.65409106507813131, 3.495898124521934782e-4]; 333 static immutable double[4] s = [-25.579105509976461286, 334 212.57260432226544008, -610.69018684944109624, 422.69668805777760407]; 335 /* T - Approximation for SINH(Y)/Y */ 336 static immutable double[6] t = [1.6125990452916363814e-10, 337 2.5051878502858255354e-8, 2.7557319615147964774e-6, 338 1.9841269840928373686e-4, .0083333333333334751799, 339 .16666666666666666446]; 340 /*---------------------------------------------------------------------*/ 341 static immutable double[6] estm = [52.0583, 5.7607, 2.7782, 14.4303, 185.3004, 342 9.3715]; 343 static immutable double[7] estf = [41.8341, 7.1075, 6.4306, 42.511, 1.35633, 344 84.5096, 20.]; 345 /* Local variables */ 346 ptrdiff_t iend, i, j, k, m, ii, mplus1; 347 double x2b4, twox, c, blpha, ratio, wminf; 348 double d1, d2, d3, f0, f1, f2, p0, q0, t1, t2, twonu; 349 double dm, ex, b1, b2, nu; 350 ii = 0; /* -Wall */ 351 ex = x; 352 nu = alpha; 353 ncalc = -2; 354 if (b.length > 0 && (0. <= nu && nu < 1.)) 355 { 356 if (ex <= 0 || (expFlag == false && ex > xmax_BESS_K)) 357 { 358 if (ex <= 0) 359 { 360 enforce(ex == 0, "besselK Range Error"); 361 for (i = 0; i < b.length; i++) 362 b[i] = double.infinity; 363 } 364 else /* would only have underflow */ 365 for (i = 0; i < b.length; i++) 366 b[i] = 0; 367 ncalc = b.length; 368 return; 369 } 370 k = 0; 371 if (nu < sqxmin_BESS_K) 372 { 373 nu = 0; 374 } 375 else if (nu > .5) 376 { 377 k = 1; 378 nu -= 1.; 379 } 380 twonu = nu + nu; 381 iend = b.length + k - 1; 382 c = nu * nu; 383 d3 = -c; 384 if (ex <= 1.) 385 { 386 /* ------------------------------------------------------------ 387 Calculation of P0 = GAMMA(1+ALPHA) * (2/X)**ALPHA 388 Q0 = GAMMA(1-ALPHA) * (X/2)**ALPHA 389 ------------------------------------------------------------ */ 390 d1 = 0; 391 d2 = p[0]; 392 t1 = 1.; 393 t2 = q[0]; 394 for (i = 2; i <= 7; i += 2) 395 { 396 d1 = c * d1 + p[i - 1]; 397 d2 = c * d2 + p[i]; 398 t1 = c * t1 + q[i - 1]; 399 t2 = c * t2 + q[i]; 400 } 401 d1 = nu * d1; 402 t1 = nu * t1; 403 f1 = log(ex); 404 f0 = a + nu * (p[7] - nu * (d1 + d2) / (t1 + t2)) - f1; 405 q0 = exp(-nu * (a - nu * (p[7] + nu * (d1 - d2) / (t1 - t2)) - f1)); 406 f1 = nu * f0; 407 p0 = exp(f1); 408 /* ----------------------------------------------------------- 409 Calculation of F0 = 410 ----------------------------------------------------------- */ 411 d1 = r[4]; 412 t1 = 1.; 413 for (i = 0; i < 4; ++i) 414 { 415 d1 = c * d1 + r[i]; 416 t1 = c * t1 + s[i]; 417 } 418 /* d2 := sinh(f1)/ nu = sinh(f1)/(f1/f0) 419 * = f0 * sinh(f1)/f1 */ 420 if (fabs(f1) <= .5) 421 { 422 f1 *= f1; 423 d2 = 0; 424 for (i = 0; i < 6; ++i) 425 { 426 d2 = f1 * d2 + t[i]; 427 } 428 d2 = f0 + f0 * f1 * d2; 429 } 430 else 431 { 432 d2 = sinh(f1) / nu; 433 } 434 f0 = d2 - nu * d1 / (t1 * p0); 435 if (ex <= 1e-10) 436 { 437 /* --------------------------------------------------------- 438 X <= 1.0E-10 439 Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X) 440 --------------------------------------------------------- */ 441 b[0] = f0 + ex * f0; 442 if (expFlag == false) 443 { 444 b[0] -= ex * b[0]; 445 } 446 ratio = p0 / f0; 447 c = ex * double.max; 448 if (k != 0) 449 { 450 /* --------------------------------------------------- 451 Calculation of K(ALPHA,X) 452 and X*K(ALPHA+1,X)/K(ALPHA,X), ALPHA >= 1/2 453 --------------------------------------------------- */ 454 ncalc = -1; 455 if (b[0] >= c / ratio) 456 { 457 return; 458 } 459 b[0] = ratio * b[0] / ex; 460 twonu += 2.; 461 ratio = twonu; 462 } 463 ncalc = 1; 464 if (b.length == 1) 465 return; 466 /* ----------------------------------------------------- 467 Calculate K(ALPHA+L,X)/K(ALPHA+L-1,X), 468 L = 1, 2, ... , NB-1 469 ----------------------------------------------------- */ 470 ncalc = -1; 471 for (i = 1; i < b.length; ++i) 472 { 473 if (ratio >= c) 474 return; 475 b[i] = ratio / ex; 476 twonu += 2.; 477 ratio = twonu; 478 } 479 ncalc = 1; 480 goto L420; 481 } 482 else 483 { 484 /* ------------------------------------------------------ 485 10^-10 < X <= 1.0 486 ------------------------------------------------------ */ 487 c = 1.; 488 x2b4 = ex * ex / 4.; 489 p0 = .5 * p0; 490 q0 = .5 * q0; 491 d1 = -1.; 492 d2 = 0; 493 b1 = 0; 494 b2 = 0; 495 f1 = f0; 496 f2 = p0; 497 do 498 { 499 d1 += 2.; 500 d2 += 1.; 501 d3 = d1 + d3; 502 c = x2b4 * c / d2; 503 f0 = (d2 * f0 + p0 + q0) / d3; 504 p0 /= d2 - nu; 505 q0 /= d2 + nu; 506 t1 = c * f0; 507 t2 = c * (p0 - d2 * f0); 508 b1 += t1; 509 b2 += t2; 510 } 511 while (fabs(t1 / (f1 + b1)) > double.epsilon || fabs(t2 / ( 512 f2 + b2)) > double.epsilon); 513 b1 = f1 + b1; 514 b2 = 2. * (f2 + b2) / ex; 515 if (expFlag) 516 { 517 d1 = exp(ex); 518 b1 *= d1; 519 b2 *= d1; 520 } 521 wminf = estf[0] * ex + estf[1]; 522 } 523 } 524 else if (double.epsilon * ex > 1.) 525 { 526 /* ------------------------------------------------- 527 X > 1./EPS 528 ------------------------------------------------- */ 529 ncalc = b.length; 530 b1 = 1. / (M_SQRT_2dPI * sqrt(ex)); 531 for (i = 0; i < b.length; ++i) 532 b[i] = b1; 533 return; 534 } 535 else 536 { 537 /* ------------------------------------------------------- 538 X > 1.0 539 ------------------------------------------------------- */ 540 twox = ex + ex; 541 blpha = 0; 542 ratio = 0; 543 if (ex <= 4.) 544 { 545 /* ---------------------------------------------------------- 546 Calculation of K(ALPHA+1,X)/K(ALPHA,X), 1.0 <= X <= 4.0 547 ----------------------------------------------------------*/ 548 d2 = trunc(estm[0] / ex + estm[1]); 549 m = cast(int) d2; 550 d1 = d2 + d2; 551 d2 -= .5; 552 d2 *= d2; 553 for (i = 2; i <= m; ++i) 554 { 555 d1 -= 2.; 556 d2 -= d1; 557 ratio = (d3 + d2) / (twox + d1 - ratio); 558 } 559 /* ----------------------------------------------------------- 560 Calculation of I(|ALPHA|,X) and I(|ALPHA|+1,X) by backward 561 recurrence and K(ALPHA,X) from the wronskian 562 -----------------------------------------------------------*/ 563 d2 = trunc(estm[2] * ex + estm[3]); 564 m = cast(int) d2; 565 c = fabs(nu); 566 d3 = c + c; 567 d1 = d3 - 1.; 568 f1 = double.min_normal; 569 f0 = (2. * (c + d2) / ex + .5 * ex / (c + d2 + 1.)) * double.min_normal; 570 for (i = 3; i <= m; ++i) 571 { 572 d2 -= 1.; 573 f2 = (d3 + d2 + d2) * f0; 574 blpha = (1. + d1 / d2) * (f2 + blpha); 575 f2 = f2 / ex + f1; 576 f1 = f0; 577 f0 = f2; 578 } 579 f1 = (d3 + 2.) * f0 / ex + f1; 580 d1 = 0; 581 t1 = 1.; 582 for (i = 1; i <= 7; ++i) 583 { 584 d1 = c * d1 + p[i - 1]; 585 t1 = c * t1 + q[i - 1]; 586 } 587 p0 = exp(c * (a + c * (p[7] - c * d1 / t1) - log(ex))) / ex; 588 f2 = (c + .5 - ratio) * f1 / ex; 589 b1 = p0 + (d3 * f0 - f2 + f0 + blpha) / (f2 + f1 + f0) * p0; 590 if (expFlag == false) 591 { 592 b1 *= exp(-ex); 593 } 594 wminf = estf[2] * ex + estf[3]; 595 } 596 else 597 { 598 /* --------------------------------------------------------- 599 Calculation of K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X), by 600 backward recurrence, for X > 4.0 601 ----------------------------------------------------------*/ 602 dm = trunc(estm[4] / ex + estm[5]); 603 m = cast(int) dm; 604 d2 = dm - .5; 605 d2 *= d2; 606 d1 = dm + dm; 607 for (i = 2; i <= m; ++i) 608 { 609 dm -= 1.; 610 d1 -= 2.; 611 d2 -= d1; 612 ratio = (d3 + d2) / (twox + d1 - ratio); 613 blpha = (ratio + ratio * blpha) / dm; 614 } 615 b1 = 1. / ((M_SQRT_2dPI + M_SQRT_2dPI * blpha) * sqrt(ex)); 616 if (expFlag == false) 617 b1 *= exp(-ex); 618 wminf = estf[4] * (ex - fabs(ex - estf[6])) + estf[5]; 619 } 620 /* --------------------------------------------------------- 621 Calculation of K(ALPHA+1,X) 622 from K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X) 623 --------------------------------------------------------- */ 624 b2 = b1 + b1 * (nu + .5 - ratio) / ex; 625 } 626 /*-------------------------------------------------------------------- 627 Calculation of 'NCALC', K(ALPHA+I,X), I = 0, 1, ... , NCALC-1, 628 & K(ALPHA+I,X)/K(ALPHA+I-1,X), I = NCALC, NCALC+1, ... , NB-1 629 -------------------------------------------------------------------*/ 630 ncalc = b.length; 631 b[0] = b1; 632 if (iend == 0) 633 return; 634 j = 1 - k; 635 if (j >= 0) 636 b[j] = b2; 637 if (iend == 1) 638 return; 639 m = min(cast(int)(wminf - nu), iend); 640 for (i = 2; i <= m; ++i) 641 { 642 t1 = b1; 643 b1 = b2; 644 twonu += 2.; 645 if (ex < 1.) 646 { 647 if (b1 >= double.max / twonu * ex) 648 break; 649 } 650 else 651 { 652 if (b1 / ex >= double.max / twonu) 653 break; 654 } 655 b2 = twonu / ex * b1 + t1; 656 ii = i; 657 ++j; 658 if (j >= 0) 659 { 660 b[j] = b2; 661 } 662 } 663 m = ii; 664 if (m == iend) 665 { 666 return; 667 } 668 ratio = b2 / b1; 669 mplus1 = m + 1; 670 ncalc = -1; 671 for (i = mplus1; i <= iend; ++i) 672 { 673 twonu += 2.; 674 ratio = twonu / ex + 1. / ratio; 675 ++j; 676 if (j >= 1) 677 { 678 b[j] = ratio; 679 } 680 else 681 { 682 if (b2 >= double.max / ratio) 683 return; 684 b2 *= ratio; 685 } 686 } 687 ncalc = max(1, mplus1 - k); 688 if (ncalc == 1) 689 b[0] = b2; 690 if (b.length == 1) 691 return; 692 L420: for (i = ncalc; i < b.length; ++i) 693 { 694 /* i == ncalc */ 695 if (b[i - 1] >= double.max / b[i]) 696 return; 697 b[i] *= b[i - 1]; 698 ncalc++; 699 } 700 } 701 } 702 703 /// 704 double besselI(double x, double alpha, 705 Flag!"ExponentiallyScaled" expFlag = Flag!"ExponentiallyScaled".no) 706 { 707 ptrdiff_t nb; 708 ptrdiff_t ncalc; 709 double na; 710 double* b; 711 /* NaNs propagated correctly */ 712 if (std.math.isNaN(x) || std.math.isNaN(alpha)) 713 return x + alpha; 714 if (x < 0) 715 return double.nan; 716 na = floor(alpha); 717 if (alpha < 0) 718 { 719 /* Using Aaamowitz & Stegun 9.6.2 & 9.6.6 720 * this may not be quite optimal (CPU and accuracy wise) */ 721 return (besselI(x, -alpha, expFlag) + ((alpha == na) ? /* sin(pi * alpha) = 0 */ 0 : besselK(x, -alpha, 722 expFlag) * (expFlag == true ? 2. : 2. * exp(-2.0 * x)) / cast(double) std.math.PI * sinPi(-alpha))); 723 } 724 nb = 1 + cast(int) na; /* b.length-1 <= alpha < b.length */ 725 alpha -= nb - 1; 726 b = cast(double* ) calloc(nb, double.sizeof).enforce("besselI allocation error"); 727 scope(exit) b.free; 728 besselI(x, alpha, expFlag, b[0 .. nb], ncalc); 729 if (ncalc != nb) 730 { 731 /* error input */ 732 if (ncalc < 0) 733 throw new Exception(format("besselI(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n", 734 x, ncalc, nb, alpha)); 735 else throw new Exception(format( 736 "besselI(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1))); 737 } 738 x = b[nb - 1]; 739 return x; 740 } 741 742 /// 743 unittest 744 { 745 assert(std.math.approxEqual(besselI(2, 1), 746 1.590636854637329063382254424999666)); 747 assert(std.math.approxEqual(besselI(20, 1), 748 4.245497338512777018140990665e+7)); 749 assert(std.math.approxEqual(besselI(1, 1.2), 0.441739)); 750 } 751 752 /* modified version of besselI that accepts a work array instead of 753 allocating one. */ 754 double besselI(double x, double alpha, Flag!"ExponentiallyScaled" expFlag, double[] b) 755 { 756 ptrdiff_t nb; 757 ptrdiff_t ncalc; 758 double na; 759 /* NaNs propagated correctly */ 760 if (std.math.isNaN(x) || std.math.isNaN(alpha)) 761 return x + alpha; 762 if (x < 0) 763 return double.nan; 764 na = floor(alpha); 765 if (alpha < 0) 766 { 767 /* Using Aaamowitz & Stegun 9.6.2 & 9.6.6 768 * this may not be quite optimal (CPU and accuracy wise) */ 769 return (besselI(x, -alpha, expFlag, b) + ((alpha == na) ? 0 : besselK(x, 770 -alpha, expFlag, b) * (expFlag == true ? 2. : 2. * exp(-2.0 * x)) / cast( 771 double) std.math.PI * sinPi(-alpha))); 772 } 773 nb = 1 + cast(int) na; /* b.length-1 <= alpha < b.length */ 774 alpha -= nb - 1; 775 besselI(x, alpha, expFlag, b[0 .. nb], ncalc); 776 if (ncalc != nb) 777 { 778 /* error input */ 779 if (ncalc < 0) 780 throw new Exception(format("besselI(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n", 781 x, ncalc, nb, alpha)); 782 else throw new Exception(format( 783 "besselI(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1))); 784 } 785 x = b[nb - 1]; 786 return x; 787 } 788 789 void besselI(double x, double alpha, Flag!"ExponentiallyScaled" expFlag, double[] b, 790 ref ptrdiff_t ncalc) 791 { 792 /* ------------------------------------------------------------------- 793 794 This routine calculates Bessel functions I_(N+ALPHA) (X) 795 for non-negative argument X, and non-negative order N+ALPHA, 796 with or without exponential scaling. 797 798 799 Explanation of variables in the calling sequence 800 801 X - Non-negative argument for which 802 I's or exponentially scaled I's (I*EXP(-X)) 803 are to be calculated. If I's are to be calculated, 804 X must be less than exparg_BESS (IZE=1) or xlrg_BESS_IJ (IZE=2), 805 (see bessel.h). 806 ALPHA - Fractional part of order for which 807 I's or exponentially scaled I's (I*EXP(-X)) are 808 to be calculated. 0 <= ALPHA < 1.0. 809 NB - Number of functions to be calculated, NB > 0. 810 The first function calculated is of order ALPHA, and the 811 last is of order (NB - 1 + ALPHA). 812 IZE - Type. IZE = 1 if unscaled I's are to be calculated, 813 = 2 if exponentially scaled I's are to be calculated. 814 BI - Output vector of length NB. If the routine 815 terminates normally (NCALC=NB), the vector BI contains the 816 functions I(ALPHA,X) through I(NB-1+ALPHA,X), or the 817 corresponding exponentially scaled functions. 818 NCALC - Output variable indicating possible errors. 819 Before using the vector BI, the user should check that 820 NCALC=NB, i.e., all orders have been calculated to 821 the desired accuracy. See error returns below. 822 823 824 ******************************************************************* 825 ******************************************************************* 826 827 Error returns 828 829 In case of an error, NCALC != NB, and not all I's are 830 calculated to the desired accuracy. 831 832 NCALC < 0: An argument is out of range. For example, 833 NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >= EXPARG_BESS. 834 In this case, the BI-vector is not calculated, and NCALC is 835 set to MIN0(NB,0)-1 so that NCALC != NB. 836 837 NB > NCALC > 0: Not all requested function values could 838 be calculated accurately. This usually occurs because NB is 839 much larger than ABS(X). In this case, BI[N] is calculated 840 to the desired accuracy for N <= NCALC, but precision 841 is lost for NCALC < N <= NB. If BI[N] does not vanish 842 for N > NCALC (because it is too small to be represented), 843 and BI[N]/BI[NCALC] = 10**(-K), then only the first NSIG-K 844 significant figures of BI[N] can be trusted. 845 846 847 Intrinsic functions required are: 848 849 DBLE, EXP, gamma_cody, INT, MAX, MIN, REAL, SQRT 850 851 852 Acknowledgement 853 854 This program is based on a program written b David J. 855 Sookne (2) that computes values of the Bessel functions J or 856 I of float argument and long order. Modifications include 857 the restriction of the computation to the I Bessel function 858 of non-negative float argument, the extension of the computation 859 to arbtrary positive order, the inclusion of optional 860 exponential scaling, and the elimination of most underflow. 861 An earlier version was published in (3). 862 863 References: "A Note on Backward Recurrence Algorithms," Olver, 864 F. W. J., and Sookne, D. J., Math. Comp. 26, 1972, 865 pp 941-947. 866 867 "Bessel Functions of Real Argument and Integer Order," 868 Sookne, D. J., NBS Jour. of Res. B. 77B, 1973, pp 869 125-132. 870 871 "ALGORITHM 597, Sequence of Modified Bessel Functions 872 of the First Kind," Cody, W. J., Trans. Math. Soft., 873 1983, pp. 242-245. 874 875 Latest modification: May 30, 1989 876 877 Modified b: W. J. Cody and L. Stoltz 878 Applied Mathematics Division 879 Argonne National Laboratory 880 Argonne, IL 60439 881 */ 882 883 /*------------------------------------------------------------------- 884 Mathematical constants 885 -------------------------------------------------------------------*/ 886 enum double const__ = 1.585; 887 /* Local variables */ 888 ptrdiff_t nend, intx, nbmx, k, l, n, nstart; 889 double pold, test, p, em, en, empal, emp2al, halfx, aa, bb, cc, psave, plast, tover, 890 psavel, sum, nu, twonu; 891 /*Parameter adjustments */ 892 // --b; 893 nu = alpha; 894 twonu = nu + nu; 895 /*------------------------------------------------------------------- 896 Check for X, NB, OR IZE out of range. 897 ------------------------------------------------------------------- */ 898 if (b.length > 0 && x >= 0. && (0. <= nu && nu < 1.)) 899 { 900 ncalc = b.length; 901 if (expFlag == false && x > exparg_BESS) 902 { 903 b[] = double.infinity; /* the limit *is* = Inf */ 904 return; 905 } 906 if (expFlag && x > xlrg_BESS_IJ) 907 { 908 b[] = 0; /* The limit exp(-x) * I_nu(x) --> 0 : */ 909 return; 910 } 911 intx = cast(int)(x); /* fine, since x <= xlrg_BESS_IJ <<< LONG_MAX */ 912 if (x >= rtnsig_BESS) 913 { 914 /* "non-small" x ( >= 1e-4 ) */ 915 /* ------------------------------------------------------------------- 916 Initialize the forward sweep, the P-sequence of Olver 917 ------------------------------------------------------------------- */ 918 nbmx = b.length - intx; 919 n = intx + 1; 920 en = cast(double)(n + n) + twonu; 921 plast = 1.; 922 p = en / x; 923 /* ------------------------------------------------ 924 Calculate general significance test 925 ------------------------------------------------ */ 926 test = ensig_BESS + ensig_BESS; 927 if (intx << 1 > nsig_BESS * 5) 928 { 929 test = sqrt(test * p); 930 } 931 else 932 { 933 test /= const__ ^^ intx; 934 } 935 if (nbmx >= 3) 936 { 937 /* -------------------------------------------------- 938 Calculate P-sequence until N = NB-1 939 Check for possible overflow. 940 ------------------------------------------------ */ 941 tover = enten_BESS / ensig_BESS; 942 nstart = intx + 2; 943 nend = b.length - 1; 944 for (k = nstart; k <= nend; ++k) 945 { 946 n = k; 947 en += 2.; 948 pold = plast; 949 plast = p; 950 p = en * plast / x + pold; 951 if (p > tover) 952 { 953 /* ------------------------------------------------ 954 To avoid overflow, divide P-sequence b TOVER. 955 Calculate P-sequence until ABS(P) > 1. 956 ---------------------------------------------- */ 957 tover = enten_BESS; 958 p /= tover; 959 plast /= tover; 960 psave = p; 961 psavel = plast; 962 nstart = n + 1; 963 do 964 { 965 ++n; 966 en += 2.; 967 pold = plast; 968 plast = p; 969 p = en * plast / x + pold; 970 } 971 while (p <= 1.); 972 bb = en / x; 973 /* ------------------------------------------------ 974 Calculate backward test, and find NCALC, 975 the highest N such that the test is passed. 976 ------------------------------------------------ */ 977 test = pold * plast / ensig_BESS; 978 test *= .5 - .5 / (bb * bb); 979 p = plast * tover; 980 --n; 981 en -= 2.; 982 nend = min(b.length, n); 983 for (l = nstart; l <= nend; ++l) 984 { 985 ncalc = l; 986 pold = psavel; 987 psavel = psave; 988 psave = en * psavel / x + pold; 989 if (psave * psavel > test) 990 { 991 goto L90; 992 } 993 } 994 ncalc = nend + 1; 995 L90: ncalc--; 996 goto L120; 997 } 998 } 999 n = nend; 1000 en = cast(double)(n + n) + twonu; 1001 /*--------------------------------------------------- 1002 Calculate special significance test for NBMX > 2. 1003 --------------------------------------------------- */ 1004 test = max(test, sqrt(plast * ensig_BESS) * sqrt(p + p)); 1005 } 1006 /* -------------------------------------------------------- 1007 Calculate P-sequence until significance test passed. 1008 -------------------------------------------------------- */ 1009 do 1010 { 1011 ++n; 1012 en += 2.; 1013 pold = plast; 1014 plast = p; 1015 p = en * plast / x + pold; 1016 } 1017 while (p < test); 1018 L120: /* ------------------------------------------------------------------- 1019 Initialize the backward recursion and the normalization sum. 1020 ------------------------------------------------------------------- */ 1021 ++n; 1022 en += 2.; 1023 bb = 0; 1024 aa = 1. / p; 1025 em = cast(double) n - 1.; 1026 empal = em + nu; 1027 emp2al = em - 1. + twonu; 1028 sum = aa * empal * emp2al / em; 1029 nend = n - b.length; 1030 if (nend < 0) 1031 { 1032 /* ----------------------------------------------------- 1033 N < NB, so store BI[N] and set higher orders to 0.. 1034 ----------------------------------------------------- */ 1035 b[n - 1] = aa; 1036 nend = -nend; 1037 b[n .. n + nend] = 0; 1038 } 1039 else 1040 { 1041 if (nend > 0) 1042 { 1043 /* ----------------------------------------------------- 1044 Recur backward via difference equation, 1045 calculating (but not storing) BI[N], until N = NB. 1046 --------------------------------------------------- */ 1047 for (l = 1; l <= nend; ++l) 1048 { 1049 --n; 1050 en -= 2.; 1051 cc = bb; 1052 bb = aa; 1053 /* for x ~= 1500, sum would overflow to 'inf' here, 1054 * and the final b[] /= sum would give 0 wrongly; 1055 * RE-normalize (aa, sum) here -- no need to undo */ 1056 if (nend > 100 && aa > 1e200) 1057 { 1058 /* multiply b 2^-900 = 1.18e-271 */ 1059 cc = ldexp(cc, -900); 1060 bb = ldexp(bb, -900); 1061 sum = ldexp(sum, -900); 1062 } 1063 aa = en * bb / x + cc; 1064 em -= 1.; 1065 emp2al -= 1.; 1066 if (n == 1) 1067 { 1068 break; 1069 } 1070 if (n == 2) 1071 { 1072 emp2al = 1.; 1073 } 1074 empal -= 1.; 1075 sum = (sum + aa * empal) * emp2al / em; 1076 } 1077 } 1078 /* --------------------------------------------------- 1079 Store BI[NB] 1080 --------------------------------------------------- */ 1081 b[n - 1] = aa; 1082 if (b.length <= 1) 1083 { 1084 sum = sum + sum + aa; 1085 goto L230; 1086 } 1087 /* ------------------------------------------------- 1088 Calculate and Store BI[NB-1] 1089 ------------------------------------------------- */ 1090 --n; 1091 en -= 2.; 1092 b[n - 1] = en * aa / x + bb; 1093 if (n == 1) 1094 { 1095 goto L220; 1096 } 1097 em -= 1.; 1098 if (n == 2) 1099 emp2al = 1.; 1100 else emp2al -= 1.; 1101 empal -= 1.; 1102 sum = (sum + b[n - 1] * empal) * emp2al / em; 1103 } 1104 nend = n - 2; 1105 if (nend > 0) 1106 { 1107 /* -------------------------------------------- 1108 Calculate via difference equation 1109 and store BI[N], until N = 2. 1110 ------------------------------------------ */ 1111 for (l = 1; l <= nend; ++l) 1112 { 1113 --n; 1114 en -= 2.; 1115 b[n - 1] = en * b[n + 0] / x + b[n + 1]; 1116 em -= 1.; 1117 if (n == 2) 1118 emp2al = 1.; 1119 else emp2al -= 1.; 1120 empal -= 1.; 1121 sum = (sum + b[n - 1] * empal) * emp2al / em; 1122 } 1123 } 1124 /* ---------------------------------------------- 1125 Calculate BI[1] 1126 -------------------------------------------- */ 1127 b[0] = 2. * empal * b[1] / x + b[2]; 1128 L220: sum = sum + sum + b[0]; 1129 L230: /* --------------------------------------------------------- 1130 Normalize. Divide all BI[N] b sum. 1131 --------------------------------------------------------- */ 1132 if (nu != 0.) 1133 sum *= (gamma_cody(1. + nu) * pow(x * .5, -nu)); 1134 if (expFlag == false) 1135 sum *= exp(-(x)); 1136 aa = enmten_BESS; 1137 if (sum > 1.) 1138 aa *= sum; 1139 for (n = 1; n <= b.length; ++n) 1140 { 1141 if (b[n - 1] < aa) 1142 b[n - 1] = 0; 1143 else b[n - 1] /= sum; 1144 } 1145 return; 1146 } 1147 else 1148 { 1149 /* small x < 1e-4 */ 1150 /* ----------------------------------------------------------- 1151 Two-term ascending series for small X. 1152 -----------------------------------------------------------*/ 1153 aa = 1.; 1154 empal = 1. + nu; 1155 /* No need to check for underflow */ 1156 halfx = .5 * x; 1157 if (nu != 0.) 1158 aa = pow(halfx, nu) / gamma_cody(empal); 1159 if (expFlag) 1160 aa *= exp(-(x)); 1161 bb = halfx * halfx; 1162 b[0] = aa + aa * bb / empal; 1163 if (x != 0. && b[0] == 0.) 1164 ncalc = 0; 1165 if (b.length > 1) 1166 { 1167 if (x == 0.) 1168 { 1169 b[1 .. $] = 0; 1170 } 1171 else 1172 { 1173 /* ------------------------------------------------- 1174 Calculate higher-order functions. 1175 ------------------------------------------------- */ 1176 cc = halfx; 1177 tover = (enmten_BESS + enmten_BESS) / x; 1178 if (bb != 0.) 1179 tover = enmten_BESS / bb; 1180 for (n = 2; n <= b.length; ++n) 1181 { 1182 aa /= empal; 1183 empal += 1.; 1184 aa *= cc; 1185 if (aa <= tover * empal) 1186 b[n - 1] = aa = 0; 1187 else b[n - 1] = aa + aa * bb / empal; 1188 if (b[n - 1] == 0. && ncalc > n) 1189 ncalc = n - 1; 1190 } 1191 } 1192 } 1193 } 1194 } 1195 else 1196 { 1197 /* argument out of range */ 1198 ncalc = min(b.length, 0) - 1; 1199 } 1200 } 1201 1202 /// 1203 double besselJ(double x, double alpha) 1204 { 1205 ptrdiff_t nb, ncalc; 1206 double na; 1207 double* b; 1208 /* NaNs propagated correctly */ 1209 if (std.math.isNaN(x) || std.math.isNaN(alpha)) 1210 return x + alpha; 1211 if (x < 0) 1212 return double.nan; 1213 na = floor(alpha); 1214 if (alpha < 0) 1215 { 1216 /* Using Aaamowitz & Stegun 9.1.2 1217 * this may not be quite optimal (CPU and accuracy wise) */ 1218 return (besselJ(x, -alpha) * cosPi(alpha) + ((alpha == na) ? 0 : besselY(x, 1219 -alpha) * sinPi(alpha))); 1220 } 1221 nb = 1 + cast(int) na; /* b.length-1 <= alpha < b.length */ 1222 alpha -= nb - 1; 1223 b = cast(double* ) calloc(nb, double.sizeof).enforce("besselJ allocation error"); 1224 scope(exit) b.free; 1225 besselJ(x, alpha, b[0 .. nb], ncalc); 1226 if (ncalc != nb) 1227 { 1228 /* error input */ 1229 if (ncalc < 0) 1230 throw new Exception(format("besselJ(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n", 1231 x, ncalc, nb, alpha)); 1232 else throw new Exception(format( 1233 "besselJ(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1))); 1234 } 1235 x = b[nb - 1]; 1236 return x; 1237 } 1238 1239 /// 1240 unittest 1241 { 1242 assert(std.math.approxEqual(besselJ(2, 1), 1243 0.5767248077568733872024482422691)); 1244 assert(std.math.approxEqual(besselJ(20, 1), 1245 0.066833124175850045578992974193)); 1246 assert(std.math.approxEqual(besselJ(1, 1.2), 0.351884)); 1247 } 1248 1249 /* modified version of besselJ that accepts a work array instead of 1250 allocating one. */ 1251 double besselJ(double x, double alpha, double[] b) 1252 { 1253 ptrdiff_t nb, ncalc; 1254 double na; 1255 /* NaNs propagated correctly */ 1256 if (std.math.isNaN(x) || std.math.isNaN(alpha)) 1257 return x + alpha; 1258 if (x < 0) 1259 return double.nan; 1260 na = floor(alpha); 1261 if (alpha < 0) 1262 { 1263 /* Using Aaamowitz & Stegun 9.1.2 1264 * this may not be quite optimal (CPU and accuracy wise) */ 1265 return (besselJ(x, -alpha, b) * cosPi(alpha) + ((alpha == na) ? 0 : besselY(x, 1266 -alpha, b) * sinPi(alpha))); 1267 } 1268 nb = 1 + cast(int) na; /* b.length-1 <= alpha < b.length */ 1269 alpha -= nb - 1; 1270 besselJ(x, alpha, b[0 .. nb], ncalc); 1271 if (ncalc != nb) 1272 { 1273 /* error input */ 1274 if (ncalc < 0) 1275 throw new Exception(format("besselJ(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n", 1276 x, ncalc, nb, alpha)); 1277 else throw new Exception(format( 1278 "besselJ(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1))); 1279 } 1280 x = b[b.length - 1]; 1281 return x; 1282 } 1283 1284 void besselJ(double x, double alpha, double[] b, ref ptrdiff_t ncalc) 1285 { 1286 /* 1287 Calculates Bessel functions J_{n+alpha} (x) 1288 for non-negative argument x, and non-negative order n+alpha, n = 0,1,..,b.length-1. 1289 1290 Explanation of variables in the calling sequence. 1291 1292 X - Non-negative argument for which J's are to be calculated. 1293 ALPHA - Fractional part of order for which 1294 J's are to be calculated. 0 <= ALPHA < 1. 1295 NB - Number of functions to be calculated, NB >= 1. 1296 The first function calculated is of order ALPHA, and the 1297 last is of order (NB - 1 + ALPHA). 1298 B - Output vector of length NB. If RJBESL 1299 terminates normally (NCALC=NB), the vector B contains the 1300 functions J/ALPHA/(X) through J/NB-1+ALPHA/(X). 1301 NCALC - Output variable indicating possible errors. 1302 Before using the vector B, the user should check that 1303 NCALC=NB, i.e., all orders have been calculated to 1304 the desired accuracy. See the following 1305 1306 **************************************************************** 1307 1308 Error return codes 1309 1310 In case of an error, NCALC != NB, and not all J's are 1311 calculated to the desired accuracy. 1312 1313 NCALC < 0: An argument is out of range. For example, 1314 NBES <= 0, ALPHA < 0 or > 1, or X is too large. 1315 In this case, b[1] is set to zero, the remainder of the 1316 B-vector is not calculated, and NCALC is set to 1317 MIN(NB,0)-1 so that NCALC != NB. 1318 1319 NB > NCALC > 0: Not all requested function values could 1320 be calculated accurately. This usually occurs because NB is 1321 much larger than ABS(X). In this case, b[N] is calculated 1322 to the desired accuracy for N <= NCALC, but precision 1323 is lost for NCALC < N <= NB. If b[N] does not vanish 1324 for N > NCALC (because it is too small to be represented), 1325 and b[N]/b[NCALC] = 10^(-K), then only the first NSIG - K 1326 significant figures of b[N] can be trusted. 1327 1328 1329 Acknowledgement 1330 1331 This program is based on a program written b David J. Sookne 1332 (2) that computes values of the Bessel functions J or I of float 1333 argument and long order. Modifications include the restriction 1334 of the computation to the J Bessel function of non-negative float 1335 argument, the extension of the computation to arbtrary positive 1336 order, and the elimination of most underflow. 1337 1338 References: 1339 1340 Olver, F.W.J., and Sookne, D.J. (1972) 1341 "A Note on Backward Recurrence Algorithms"; 1342 Math. Comp. 26, 941-947. 1343 1344 Sookne, D.J. (1973) 1345 "Bessel Functions of Real Argument and Integer Order"; 1346 NBS Jour. of Res. B. 77B, 125-132. 1347 1348 Latest modification: March 19, 1990 1349 1350 Author: W. J. Cody 1351 Applied Mathematics Division 1352 Argonne National Laboratory 1353 Argonne, IL 60439 1354 ******************************************************************* 1355 */ 1356 1357 /* --------------------------------------------------------------------- 1358 Mathematical constants 1359 1360 PI2 = 2 / PI 1361 TWOPI1 = first few significant digits of 2 * PI 1362 TWOPI2 = (2*PI - TWOPI1) to working precision, i.e., 1363 TWOPI1 + TWOPI2 = 2 * PI to extra precision. 1364 --------------------------------------------------------------------- */ 1365 static immutable double pi2 = .636619772367581343075535; 1366 static immutable double twopi1 = 6.28125; 1367 static immutable double twopi2 = .001935307179586476925286767; 1368 /*--------------------------------------------------------------------- 1369 * Factorial(N) 1370 *--------------------------------------------------------------------- */ 1371 static immutable double[25] fact = [1., 1., 2., 6., 24., 120., 720., 5040., 40320., 1372 362880., 3628800., 39916800., 479001600., 6227020800., 87178291200., 1373 1.307674368e12, 2.0922789888e13, 3.55687428096e14, 6.402373705728e15, 1374 1.21645100408832e17, 2.43290200817664e18, 5.109094217170944e19, 1375 1.12400072777760768e21, 1376 2.585201673888497664e22, 6.2044840173323943936e23]; 1377 /* Local variables */ 1378 ptrdiff_t nend, intx, nbmx, i, j, k, l, m, n, nstart; 1379 double nu, twonu, capp, capq, pold, vcos, test, vsin; 1380 double p, s, t, z, alpem, halfx, aa, bb, cc, psave, plast; 1381 double tover, t1, alp2em, em, en, xc, xk, xm, psavel, gnu, xin, sum; 1382 /* Parameter adjustment */ 1383 // --b; 1384 nu = alpha; 1385 twonu = nu + nu; 1386 /*------------------------------------------------------------------- 1387 Check for out of range arguments. 1388 -------------------------------------------------------------------*/ 1389 if (b.length > 0 && x >= 0. && 0. <= nu && nu < 1.) 1390 { 1391 ncalc = b.length; 1392 if (x > xlrg_BESS_IJ) 1393 { 1394 /* indeed, the limit is 0, 1395 * but the cutoff happens too early */ 1396 b[] = double.infinity; 1397 return; 1398 } 1399 intx = cast(int)(x); 1400 /* Initialize result array to zero. */ 1401 b[] = 0; 1402 /*=================================================================== 1403 Branch into 3 cases : 1404 1) use 2-term ascending series for small X 1405 2) use asymptotic form for large X when NB is not too large 1406 3) use recursion otherwise 1407 ===================================================================*/ 1408 if (x < rtnsig_BESS) 1409 { 1410 /* --------------------------------------------------------------- 1411 Two-term ascending series for small X. 1412 --------------------------------------------------------------- */ 1413 alpem = 1. + nu; 1414 halfx = (x > enmten_BESS) ? .5 * x : 0; 1415 aa = (nu != 0.) ? pow(halfx, nu) / (nu * gamma_cody(nu)) : 1.; 1416 bb = (x + 1. > 1.) ? -halfx * halfx : 0; 1417 b[0] = aa + aa * bb / alpem; 1418 if (x != 0. && b[0] == 0.) 1419 ncalc = 0; 1420 if (b.length != 1) 1421 { 1422 if (x <= 0.) 1423 { 1424 for (n = 2; n <= b.length; ++n) 1425 b[n - 1] = 0; 1426 } 1427 else 1428 { 1429 /* ---------------------------------------------- 1430 Calculate higher order functions. 1431 ---------------------------------------------- */ 1432 if (bb == 0.) 1433 tover = (enmten_BESS + enmten_BESS) / x; 1434 else tover = enmten_BESS / bb; 1435 cc = halfx; 1436 for (n = 2; n <= b.length; ++n) 1437 { 1438 aa /= alpem; 1439 alpem += 1.; 1440 aa *= cc; 1441 if (aa <= tover * alpem) 1442 aa = 0; 1443 b[n - 1] = aa + aa * bb / alpem; 1444 if (b[n - 1] == 0. && ncalc > n) 1445 ncalc = n - 1; 1446 } 1447 } 1448 } 1449 } 1450 else if (x > 25. && b.length <= intx + 1) 1451 { 1452 /* ------------------------------------------------------------ 1453 Asymptotic series for X > 25 (and not too large b.length) 1454 ------------------------------------------------------------ */ 1455 xc = sqrt(pi2 / x); 1456 xin = 1 / (64 * x * x); 1457 if (x >= 130.) 1458 m = 4; 1459 else if (x >= 35.) 1460 m = 8; 1461 else m = 11; 1462 xm = 4. * cast(double) m; 1463 /* ------------------------------------------------ 1464 Argument reduction for SIN and COS routines. 1465 ------------------------------------------------ */ 1466 t = trunc(x / (twopi1 + twopi2) + .5); 1467 z = (x - t * twopi1) - t * twopi2 - (nu + .5) / pi2; 1468 vsin = sin(z); 1469 vcos = cos(z); 1470 gnu = twonu; 1471 for (i = 1; i <= 2; ++i) 1472 { 1473 s = (xm - 1. - gnu) * (xm - 1. + gnu) * xin * .5; 1474 t = (gnu - (xm - 3.)) * (gnu + (xm - 3.)); 1475 t1 = (gnu - (xm + 1.)) * (gnu + (xm + 1.)); 1476 k = m + m; 1477 capp = s * t / fact[k]; 1478 capq = s * t1 / fact[k + 1]; 1479 xk = xm; 1480 for (; k >= 4; k -= 2) 1481 { 1482 /* k + 2(j-2) == 2m */ 1483 xk -= 4.; 1484 s = (xk - 1. - gnu) * (xk - 1. + gnu); 1485 t1 = t; 1486 t = (gnu - (xk - 3.)) * (gnu + (xk - 3.)); 1487 capp = (capp + 1. / fact[k - 2]) * s * t * xin; 1488 capq = (capq + 1. / fact[k - 1]) * s * t1 * xin; 1489 } 1490 capp += 1.; 1491 capq = (capq + 1.) * (gnu * gnu - 1.) * (.125 / x); 1492 b[i - 1] = xc * (capp * vcos - capq * vsin); 1493 if (b.length == 1) 1494 return; 1495 /* vsin <--> vcos */ t = vsin; 1496 vsin = -vcos; 1497 vcos = t; 1498 gnu += 2.; 1499 } 1500 /* ----------------------------------------------- 1501 If NB > 2, compute J(X,ORDER+I) for I = 2, NB-1 1502 ----------------------------------------------- */ 1503 if (b.length > 2) 1504 for (gnu = twonu + 2., j = 2; j <= b.length; j++, gnu += 2.) 1505 b[j] = gnu * b[j - 1] / x - b[j - 2]; 1506 } 1507 else 1508 { 1509 /* rtnsig_BESS <= x && ( x <= 25 || intx+1 < b.length ) : 1510 -------------------------------------------------------- 1511 Use recurrence to generate results. 1512 First initialize the calculation of P*S. 1513 -------------------------------------------------------- */ 1514 nbmx = b.length - intx; 1515 n = intx + 1; 1516 en = cast(double)(n + n) + twonu; 1517 plast = 1.; 1518 p = en / x; 1519 /* --------------------------------------------------- 1520 Calculate general significance test. 1521 --------------------------------------------------- */ 1522 test = ensig_BESS + ensig_BESS; 1523 if (nbmx >= 3) 1524 { 1525 /* ------------------------------------------------------------ 1526 Calculate P*S until N = NB-1. Check for possible overflow. 1527 ---------------------------------------------------------- */ 1528 tover = enten_BESS / ensig_BESS; 1529 nstart = intx + 2; 1530 nend = b.length - 1; 1531 en = cast(double)(nstart + nstart) - 2. + twonu; 1532 for (k = nstart; k <= nend; ++k) 1533 { 1534 n = k; 1535 en += 2.; 1536 pold = plast; 1537 plast = p; 1538 p = en * plast / x - pold; 1539 if (p > tover) 1540 { 1541 /* ------------------------------------------- 1542 To avoid overflow, divide P*S b TOVER. 1543 Calculate P*S until ABS(P) > 1. 1544 -------------------------------------------*/ 1545 tover = enten_BESS; 1546 p /= tover; 1547 plast /= tover; 1548 psave = p; 1549 psavel = plast; 1550 nstart = n + 1; 1551 do 1552 { 1553 ++n; 1554 en += 2.; 1555 pold = plast; 1556 plast = p; 1557 p = en * plast / x - pold; 1558 } 1559 while (p <= 1.); 1560 bb = en / x; 1561 /* ----------------------------------------------- 1562 Calculate backward test and find NCALC, 1563 the highest N such that the test is passed. 1564 ----------------------------------------------- */ 1565 test = pold * plast * (.5 - .5 / (bb * bb)); 1566 test /= ensig_BESS; 1567 p = plast * tover; 1568 --n; 1569 en -= 2.; 1570 nend = min(b.length, n); 1571 for (l = nstart; l <= nend; ++l) 1572 { 1573 pold = psavel; 1574 psavel = psave; 1575 psave = en * psavel / x - pold; 1576 if (psave * psavel > test) 1577 { 1578 ncalc = l - 1; 1579 goto L190; 1580 } 1581 } 1582 ncalc = nend; 1583 goto L190; 1584 } 1585 } 1586 n = nend; 1587 en = cast(double)(n + n) + twonu; 1588 /* ----------------------------------------------------- 1589 Calculate special significance test for NBMX > 2. 1590 -----------------------------------------------------*/ 1591 test = max(test, sqrt(plast * ensig_BESS) * sqrt(p + p)); 1592 } 1593 /* ------------------------------------------------ 1594 Calculate P*S until significance test passes. */ 1595 do 1596 { 1597 ++n; 1598 en += 2.; 1599 pold = plast; 1600 plast = p; 1601 p = en * plast / x - pold; 1602 } 1603 while (p < test); 1604 L190: /*--------------------------------------------------------------- 1605 Initialize the backward recursion and the normalization sum. 1606 --------------------------------------------------------------- */ 1607 ++n; 1608 en += 2.; 1609 bb = 0; 1610 aa = 1. / p; 1611 m = n / 2; 1612 em = cast(double) m; 1613 m = (n << 1) - (m << 2); /* = 2 n - 4 (n/2) 1614 = 0 for even, 2 for odd n */ 1615 if (m == 0) 1616 sum = 0; 1617 else 1618 { 1619 alpem = em - 1. + nu; 1620 alp2em = em + em + nu; 1621 sum = aa * alpem * alp2em / em; 1622 } 1623 nend = n - b.length; 1624 /* if (nend > 0) */ 1625 /* -------------------------------------------------------- 1626 Recur backward via difference equation, calculating 1627 (but not storing) b[N], until N = NB. 1628 -------------------------------------------------------- */ 1629 for (l = 1; l <= nend; ++l) 1630 { 1631 --n; 1632 en -= 2.; 1633 cc = bb; 1634 bb = aa; 1635 aa = en * bb / x - cc; 1636 m = m ? 0 : 2; /* m = 2 - m failed on gcc4-20041019 */ 1637 if (m != 0) 1638 { 1639 em -= 1.; 1640 alp2em = em + em + nu; 1641 if (n == 1) 1642 break; 1643 alpem = em - 1. + nu; 1644 if (alpem == 0.) 1645 alpem = 1.; 1646 sum = (sum + aa * alp2em) * alpem / em; 1647 } 1648 } 1649 /*-------------------------------------------------- 1650 Store b[NB]. 1651 --------------------------------------------------*/ 1652 b[n - 1] = aa; 1653 if (nend >= 0) 1654 { 1655 if (b.length <= 1) 1656 { 1657 if (nu + 1. == 1.) 1658 alp2em = 1.; 1659 else alp2em = nu; 1660 sum += b[0] * alp2em; 1661 goto L250; 1662 } 1663 else 1664 { 1665 /*-- b.length >= 2 : --------------------------- 1666 Calculate and store b[NB-1]. 1667 ----------------------------------------*/ 1668 --n; 1669 en -= 2.; 1670 b[n - 1] = en * aa / x - bb; 1671 if (n == 1) 1672 goto L240; 1673 m = m ? 0 : 2; /* m = 2 - m failed on gcc4-20041019 */ 1674 if (m != 0) 1675 { 1676 em -= 1.; 1677 alp2em = em + em + nu; 1678 alpem = em - 1. + nu; 1679 if (alpem == 0.) 1680 alpem = 1.; 1681 sum = (sum + b[n - 1] * alp2em) * alpem / em; 1682 } 1683 } 1684 } 1685 /* if (n - 2 != 0) */ 1686 /* -------------------------------------------------------- 1687 Calculate via difference equation and store b[N], 1688 until N = 2. 1689 -------------------------------------------------------- */ 1690 for (n = n - 1; n >= 2; n--) 1691 { 1692 en -= 2.; 1693 b[n - 1] = en * b[n + 0] / x - b[n + 1]; 1694 m = m ? 0 : 2; /* m = 2 - m failed on gcc4-20041019 */ 1695 if (m != 0) 1696 { 1697 em -= 1.; 1698 alp2em = em + em + nu; 1699 alpem = em - 1. + nu; 1700 if (alpem == 0.) 1701 alpem = 1.; 1702 sum = (sum + b[n - 1] * alp2em) * alpem / em; 1703 } 1704 } 1705 /* --------------------------------------- 1706 Calculate b[1]. 1707 -----------------------------------------*/ 1708 b[0] = 2. * (nu + 1.) * b[1] / x - b[2]; 1709 L240: em -= 1.; 1710 alp2em = em + em + nu; 1711 if (alp2em == 0.) 1712 alp2em = 1.; 1713 sum += b[0] * alp2em; 1714 L250: /* --------------------------------------------------- 1715 Normalize. Divide all b[N] b sum. 1716 ---------------------------------------------------*/ 1717 1718 /* if (nu + 1. != 1.) poor test */ 1719 if (fabs(nu) > 1e-15) 1720 sum *= (gamma_cody(nu) * pow(.5 * x, -nu)); 1721 aa = enmten_BESS; 1722 if (sum > 1.) 1723 aa *= sum; 1724 for (n = 1; n <= b.length; ++n) 1725 { 1726 if (fabs(b[n - 1]) < aa) 1727 b[n - 1] = 0; 1728 else b[n - 1] /= sum; 1729 } 1730 } 1731 } 1732 else 1733 { 1734 /* Error return -- X, NB, or ALPHA is out of range : */ 1735 b[0] = 0; 1736 ncalc = min(b.length, 0) - 1; 1737 } 1738 } 1739 1740 /// 1741 double besselY(double x, double alpha) 1742 { 1743 ptrdiff_t nb, ncalc; 1744 double na; 1745 double* b; 1746 /* NaNs propagated correctly */ 1747 if (std.math.isNaN(x) || std.math.isNaN(alpha)) 1748 return x + alpha; 1749 if (x < 0) 1750 return double.nan; 1751 na = floor(alpha); 1752 if (alpha < 0) 1753 { 1754 /* Using Aaamowitz & Stegun 9.1.2 1755 * this may not be quite optimal (CPU and accuracy wise) */ 1756 return (besselY(x, -alpha) * cosPi(alpha) - ((alpha == na) ? 0 : besselJ(x, 1757 -alpha) * sinPi(alpha))); 1758 } 1759 nb = 1 + cast(int) na; /* b.length-1 <= alpha < b.length */ 1760 alpha -= nb - 1; 1761 b = cast(double* ) calloc(nb, double.sizeof).enforce("besselY allocation error"); 1762 scope(exit) b.free; 1763 besselY(x, alpha, b[0 .. nb], ncalc); 1764 if (ncalc != nb) 1765 { 1766 /* error input */ 1767 if (ncalc == -1) 1768 return double.infinity; 1769 else if (ncalc < -1) 1770 throw new Exception(format("besselY(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n", 1771 x, ncalc, nb, alpha)); 1772 else /* ncalc >= 0 */ 1773 throw new Exception(format( 1774 "besselY(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1))); 1775 } 1776 x = b[nb - 1]; 1777 return x; 1778 } 1779 1780 /// 1781 unittest 1782 { 1783 assert(std.math.approxEqual(besselY(2, 1), 1784 -0.1070324315409375468883707722774)); 1785 assert(std.math.approxEqual(besselY(20, 1), 1786 -0.165511614362521295863976023243)); 1787 assert(std.math.approxEqual(besselY(1, 1.2), -0.901215)); 1788 } 1789 1790 /* modified version of besselY that accepts a work array instead of 1791 allocating one. */ 1792 double besselY(double x, double alpha, double[] b) 1793 { 1794 ptrdiff_t nb, ncalc; 1795 double na; 1796 /* NaNs propagated correctly */ 1797 if (std.math.isNaN(x) || std.math.isNaN(alpha)) 1798 return x + alpha; 1799 if (x < 0) 1800 return double.nan; 1801 na = floor(alpha); 1802 if (alpha < 0) 1803 { 1804 /* Using Aaamowitz & Stegun 9.1.2 1805 * this may not be quite optimal (CPU and accuracy wise) */ 1806 return (besselY(x, -alpha, b) * cosPi(alpha) - ((alpha == na) ? 0 : besselJ(x, 1807 -alpha, b) * sinPi(alpha))); 1808 } 1809 nb = 1 + cast(int) na; /* b.length-1 <= alpha < b.length */ 1810 alpha -= nb - 1; 1811 besselY(x, alpha, b[0 .. nb], ncalc); 1812 if (ncalc != nb) 1813 { 1814 /* error input */ 1815 if (ncalc == -1) 1816 return double.infinity; 1817 else if (ncalc < -1) 1818 throw new Exception(format("besselY(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n", 1819 x, ncalc, nb, alpha)); 1820 else /* ncalc >= 0 */ 1821 throw new Exception(format( 1822 "besselY(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1))); 1823 } 1824 x = b[nb - 1]; 1825 return x; 1826 } 1827 1828 void besselY(double x, double alpha, double[] b, ref ptrdiff_t ncalc) 1829 { 1830 /* ---------------------------------------------------------------------- 1831 1832 This routine calculates Bessel functions Y_(N+ALPHA) (X) 1833 v for non-negative argument X, and non-negative order N+ALPHA. 1834 1835 1836 Explanation of variables in the calling sequence 1837 1838 X - Non-negative argument for which 1839 Y's are to be calculated. 1840 ALPHA - Fractional part of order for which 1841 Y's are to be calculated. 0 <= ALPHA < 1.0. 1842 NB - Number of functions to be calculated, NB > 0. 1843 The first function calculated is of order ALPHA, and the 1844 last is of order (NB - 1 + ALPHA). 1845 BY - Output vector of length NB. If the 1846 routine terminates normally (NCALC=NB), the vector BY 1847 contains the functions Y(ALPHA,X), ... , Y(NB-1+ALPHA,X), 1848 If (0 < NCALC < NB), BY(I) contains correct function 1849 values for I <= NCALC, and contains the ratios 1850 Y(ALPHA+I-1,X)/Y(ALPHA+I-2,X) for the rest of the array. 1851 NCALC - Output variable indicating possible errors. 1852 Before using the vector BY, the user should check that 1853 NCALC=NB, i.e., all orders have been calculated to 1854 the desired accuracy. See error returns below. 1855 1856 1857 ******************************************************************* 1858 1859 Error returns 1860 1861 In case of an error, NCALC != NB, and not all Y's are 1862 calculated to the desired accuracy. 1863 1864 NCALC < -1: An argument is out of range. For example, 1865 NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >= 1866 XMAX. In this case, BY[0] = 0.0, the remainder of the 1867 BY-vector is not calculated, and NCALC is set to 1868 MIN0(NB,0)-2 so that NCALC != NB. 1869 NCALC = -1: Y(ALPHA,X) >= XINF. The requested function 1870 values are set to 0.0. 1871 1 < NCALC < NB: Not all requested function values could 1872 be calculated accurately. BY(I) contains correct function 1873 values for I <= NCALC, and and the remaining NB-NCALC 1874 array elements contain 0.0. 1875 1876 1877 Intrinsic functions required are: 1878 1879 DBLE, EXP, INT, MAX, MIN, REAL, SQRT 1880 1881 1882 Acknowledgement 1883 1884 This program draws heavily on Temme's Algol program for Y(a,x) 1885 and Y(a+1,x) and on Campbell's programs for Y_nu(x). Temme's 1886 scheme is used for x < THRESH, and Campbell's scheme is used 1887 in the asymptotic region. Segments of code from both sources 1888 have been translated into Fortran 77, merged, and heavily modified. 1889 Modifications include parameterization of machine dependencies, 1890 use of a new approximation for ln(gamma(x)), and built-in 1891 protection against over/underflow. 1892 1893 References: "Bessel functions J_nu(x) and Y_nu(x) of float 1894 order and float argument," Campbell, J. B., 1895 Comp. Phy. Comm. 18, 1979, pp. 133-142. 1896 1897 "On the numerical evaluation of the ordinary 1898 Bessel function of the second kind," Temme, 1899 N. M., J. Comput. Phys. 21, 1976, pp. 343-350. 1900 1901 Latest modification: March 19, 1990 1902 1903 Modified b: W. J. Cody 1904 Applied Mathematics Division 1905 Argonne National Laboratory 1906 Argonne, IL 60439 1907 ----------------------------------------------------------------------*/ 1908 1909 /* ---------------------------------------------------------------------- 1910 Mathematical constants 1911 FIVPI = 5*PI 1912 PIM5 = 5*PI - 15 1913 ----------------------------------------------------------------------*/ 1914 static immutable double fivpi = 15.707963267948966192; 1915 static immutable double pim5 = .70796326794896619231; 1916 /*---------------------------------------------------------------------- 1917 Coefficients for Chebshev polynomial expansion of 1918 1/gamma(1-x), abs(x) <= .5 1919 ----------------------------------------------------------------------*/ 1920 static immutable double[21] ch = [-6.7735241822398840964e-24, 1921 -6.1455180116049879894e-23, 2.9017595056104745456e-21, 1922 1.3639417919073099464e-19, 2.3826220476859635824e-18, 1923 -9.0642907957550702534e-18, 1924 -1.4943667065169001769e-15, 1925 -3.3919078305362211264e-14, 1926 -1.7023776642512729175e-13, 9.1609750938768647911e-12, 1927 2.4230957900482704055e-10, 1.7451364971382984243e-9, 1928 -3.3126119768180852711e-8, -8.6592079961391259661e-7, 1929 -4.9717367041957398581e-6, 1930 7.6309597585908126618e-5, 1931 .0012719271366545622927, 1932 .0017063050710955562222, 1933 -.07685284084478667369, 1934 -.28387654227602353814, .92187029365045265648]; 1935 /* Local variables */ 1936 ptrdiff_t i, k, na; 1937 double alfa, div, ddiv, even, gamma, term, cosmu, sinmu, a, c, d, e, f, g, h, p, q, r, 1938 s, d1, d2, q0, pa, pa1, qa, qa1, en, en1, nu, ex, ya, ya1, twobx, den, odd, aye, 1939 dmu, x2, xna; 1940 en1 = ya = ya1 = 0; /* -Wall */ 1941 ex = x; 1942 nu = alpha; 1943 if (b.length > 0 && 0. <= nu && nu < 1.) 1944 { 1945 if (ex < double.min_normal || ex > xlrg_BESS_Y) 1946 { 1947 /* Warning is not really appropriate, give 1948 * proper limit: 1949 * ML_ERROR(ME_RANGE, "besselY"); */ 1950 ncalc = b.length; 1951 if (ex > xlrg_BESS_Y) 1952 b[0] = 0; /*was double.infinity */ 1953 else if (ex < double.min_normal) 1954 b[0] = -double.infinity; 1955 for (i = 0; i < b.length; i++) 1956 b[i] = b[0]; 1957 return; 1958 } 1959 xna = trunc(nu + .5); 1960 na = cast(int) xna; 1961 if (na == 1) 1962 { 1963 /* <==> .5 <= alpha < 1 <==> -5. <= nu < 0 */ 1964 nu -= xna; 1965 } 1966 if (nu == -.5) 1967 { 1968 p = M_SQRT_2dPI / sqrt(ex); 1969 ya = p * sin(ex); 1970 ya1 = -p * cos(ex); 1971 } 1972 else if (ex < 3.) 1973 { 1974 /* ------------------------------------------------------------- 1975 Use Temme's scheme for small X 1976 ------------------------------------------------------------- */ 1977 a = ex * .5; 1978 d = -log(a); 1979 f = nu * d; 1980 e = pow(a, -nu); 1981 if (fabs(nu) < M_eps_sinc) 1982 c = std.math.M_1_PI; 1983 else c = nu / sinPi(nu); 1984 /* ------------------------------------------------------------ 1985 Computation of sinh(f)/f 1986 ------------------------------------------------------------ */ 1987 if (fabs(f) < 1.) 1988 { 1989 x2 = f * f; 1990 en = 19.; 1991 s = 1.; 1992 for (i = 1; i <= 9; ++i) 1993 { 1994 s = s * x2 / en / (en - 1.) + 1.; 1995 en -= 2.; 1996 } 1997 } 1998 else 1999 { 2000 s = (e - 1. / e) * .5 / f; 2001 } 2002 /* -------------------------------------------------------- 2003 Computation of 1/gamma(1-a) using Chebshev polynomials */ 2004 x2 = nu * nu * 8.; 2005 aye = ch[0]; 2006 even = 0; 2007 alfa = ch[1]; 2008 odd = 0; 2009 for (i = 3; i <= 19; i += 2) 2010 { 2011 even = -(aye + aye + even); 2012 aye = -even * x2 - aye + ch[i - 1]; 2013 odd = -(alfa + alfa + odd); 2014 alfa = -odd * x2 - alfa + ch[i]; 2015 } 2016 even = (even * .5 + aye) * x2 - aye + ch[20]; 2017 odd = (odd + alfa) * 2.; 2018 gamma = odd * nu + even; 2019 /* End of computation of 1/gamma(1-a) 2020 ----------------------------------------------------------- */ 2021 g = e * gamma; 2022 e = (e + 1. / e) * .5; 2023 f = 2. * c * (odd * e + even * s * d); 2024 e = nu * nu; 2025 p = g * c; 2026 q = std.math.M_1_PI / g; 2027 c = nu * std.math.PI_2; 2028 if (fabs(c) < M_eps_sinc) 2029 r = 1.; 2030 else r = sinPi(nu / 2) / c; 2031 r = cast(double) std.math.PI * c * r * r; 2032 c = 1.; 2033 d = -(a * a); 2034 h = 0; 2035 ya = f + r * q; 2036 ya1 = p; 2037 en = 1.; 2038 while (fabs(g / (1. + fabs(ya))) + fabs(h / (1. + fabs(ya1))) > double.epsilon) 2039 { 2040 f = (f * en + p + q) / (en * en - e); 2041 c *= (d / en); 2042 p /= en - nu; 2043 q /= en + nu; 2044 g = c * (f + r * q); 2045 h = c * p - en * g; 2046 ya += g; 2047 ya1 += h; 2048 en += 1.; 2049 } 2050 ya = -ya; 2051 ya1 = -ya1 / a; 2052 } 2053 else if (ex < thresh_BESS_Y) 2054 { 2055 /* -------------------------------------------------------------- 2056 Use Temme's scheme for moderate X : 3 <= x < 16 2057 -------------------------------------------------------------- */ 2058 c = (.5 - nu) * (.5 + nu); 2059 a = ex + ex; 2060 e = ex * std.math.M_1_PI * cosPi(nu) / double.epsilon; 2061 e *= e; 2062 p = 1.; 2063 q = -ex; 2064 r = 1. + ex * ex; 2065 s = r; 2066 en = 2.; 2067 while (r * en * en < e) 2068 { 2069 en1 = en + 1.; 2070 d = (en - 1. + c / en) / s; 2071 p = (en + en - p * d) / en1; 2072 q = (-a + q * d) / en1; 2073 s = p * p + q * q; 2074 r *= s; 2075 en = en1; 2076 } 2077 f = p / s; 2078 p = f; 2079 g = -q / s; 2080 q = g; 2081 L220: en -= 1.; 2082 if (en > 0.) 2083 { 2084 r = en1 * (2. - p) - 2.; 2085 s = a + en1 * q; 2086 d = (en - 1. + c / en) / (r * r + s * s); 2087 p = d * r; 2088 q = d * s; 2089 e = f + 1.; 2090 f = p * e - g * q; 2091 g = q * e + p * g; 2092 en1 = en; 2093 goto L220; 2094 } 2095 f = 1. + f; 2096 d = f * f + g * g; 2097 pa = f / d; 2098 qa = -g / d; 2099 d = nu + .5 - p; 2100 q += ex; 2101 pa1 = (pa * q - qa * d) / ex; 2102 qa1 = (qa * q + pa * d) / ex; 2103 a = ex - std.math.PI_2 * (nu + .5); 2104 c = cos(a); 2105 s = sin(a); 2106 d = M_SQRT_2dPI / sqrt(ex); 2107 ya = d * (pa * s + qa * c); 2108 ya1 = d * (qa1 * s - pa1 * c); 2109 } 2110 else 2111 { 2112 /* x > thresh_BESS_Y */ 2113 /* ---------------------------------------------------------- 2114 Use Campbell's asymptotic scheme. 2115 ---------------------------------------------------------- */ 2116 na = 0; 2117 d1 = trunc(ex / fivpi); 2118 i = cast(int) d1; 2119 dmu = ex - 15. * d1 - d1 * pim5 - (alpha + .5) * std.math.PI_2; 2120 if (i - (i / 2 << 1) == 0) 2121 { 2122 cosmu = cos(dmu); 2123 sinmu = sin(dmu); 2124 } 2125 else 2126 { 2127 cosmu = -cos(dmu); 2128 sinmu = -sin(dmu); 2129 } 2130 ddiv = 8. * ex; 2131 dmu = alpha; 2132 den = sqrt(ex); 2133 for (k = 1; k <= 2; ++k) 2134 { 2135 p = cosmu; 2136 cosmu = sinmu; 2137 sinmu = -p; 2138 d1 = (2. * dmu - 1.) * (2. * dmu + 1.); 2139 d2 = 0; 2140 div = ddiv; 2141 p = 0; 2142 q = 0; 2143 q0 = d1 / div; 2144 term = q0; 2145 for (i = 2; i <= 20; ++i) 2146 { 2147 d2 += 8.; 2148 d1 -= d2; 2149 div += ddiv; 2150 term = -term * d1 / div; 2151 p += term; 2152 d2 += 8.; 2153 d1 -= d2; 2154 div += ddiv; 2155 term *= (d1 / div); 2156 q += term; 2157 if (fabs(term) <= double.epsilon) 2158 { 2159 break; 2160 } 2161 } 2162 p += 1.; 2163 q += q0; 2164 if (k == 1) 2165 ya = M_SQRT_2dPI * (p * cosmu - q * sinmu) / den; 2166 else ya1 = M_SQRT_2dPI * (p * cosmu - q * sinmu) / den; 2167 dmu += 1.; 2168 } 2169 } 2170 if (na == 1) 2171 { 2172 h = 2. * (nu + 1.) / ex; 2173 if (h > 1.) 2174 { 2175 if (fabs(ya1) > double.max / h) 2176 { 2177 h = 0; 2178 ya = 0; 2179 } 2180 } 2181 h = h * ya1 - ya; 2182 ya = ya1; 2183 ya1 = h; 2184 } 2185 /* --------------------------------------------------------------- 2186 Now have first one or two Y's 2187 --------------------------------------------------------------- */ 2188 b[0] = ya; 2189 ncalc = 1; 2190 if (b.length > 1) 2191 { 2192 b[1] = ya1; 2193 if (ya1 != 0.) 2194 { 2195 aye = 1. + alpha; 2196 twobx = 2. / ex; 2197 ncalc = 2; 2198 for (i = 2; i < b.length; ++i) 2199 { 2200 if (twobx < 1.) 2201 { 2202 if (fabs(b[i - 1]) * twobx >= double.max / aye) 2203 goto L450; 2204 } 2205 else 2206 { 2207 if (fabs(b[i - 1]) >= double.max / aye / twobx) 2208 goto L450; 2209 } 2210 b[i] = twobx * aye * b[i - 1] - b[i - 2]; 2211 aye += 1.; 2212 ncalc++; 2213 } 2214 } 2215 } 2216 L450: for (i = ncalc; i < b.length; ++i) 2217 b[i] = -double.infinity; /* was 0 */ 2218 } 2219 else 2220 { 2221 b[0] = 0; 2222 ncalc = min(b.length, 0) - 1; 2223 } 2224 } 2225 2226 private: 2227 double gamma_cody(double x) 2228 { 2229 /* ---------------------------------------------------------------------- 2230 2231 This routine calculates the GAMMA function for a float argument X. 2232 Computation is based on an algorithm outlined in reference [1]. 2233 The program uses rational functions that approximate the GAMMA 2234 function to at least 20 significant decimal digits. Coefficients 2235 for the approximation over the interval (1,2) are unpublished. 2236 Those for the approximation for X >= 12 are from reference [2]. 2237 The accuracy achieved depends on the arithmetic system, the 2238 compiler, the intrinsic functions, and proper selection of the 2239 machine-dependent constants. 2240 2241 ******************************************************************* 2242 2243 Error returns 2244 2245 The program returns the value XINF for singularities or 2246 when overflow would occur. The computation is believed 2247 to be free of underflow and overflow. 2248 2249 Intrinsic functions required are: 2250 2251 INT, DBLE, EXP, LOG, REAL, SIN 2252 2253 2254 References: 2255 [1] "An Overview of Software Development for Special Functions", 2256 W. J. Cody, Lecture Notes in Mathematics, 506, 2257 Numerical Analysis Dundee, 1975, G. A. Watson (ed.), 2258 Springer Verlag, Berlin, 1976. 2259 2260 [2] Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968. 2261 2262 Latest modification: October 12, 1989 2263 2264 Authors: W. J. Cody and L. Stoltz 2265 Applied Mathematics Division 2266 Argonne National Laboratory 2267 Argonne, IL 60439 2268 ----------------------------------------------------------------------*/ 2269 2270 /* ---------------------------------------------------------------------- 2271 Mathematical constants 2272 ----------------------------------------------------------------------*/ 2273 enum double sqrtpi = .9189385332046727417803297; /* == ??? */ 2274 2275 /* ******************************************************************* 2276 2277 Explanation of machine-dependent constants 2278 2279 beta - radix for the floating-point representation 2280 maxexp - the smallest positive power of beta that overflows 2281 XBIG - the largest argument for which GAMMA(X) is representable 2282 in the machine, i.e., the solution to the equation 2283 GAMMA(XBIG) = beta**maxexp 2284 XINF - the largest machine representable floating-point number; 2285 approximately beta**maxexp 2286 EPS - the smallest positive floating-point number such that 1.0+EPS > 1.0 2287 XMININ - the smallest positive floating-point number such that 2288 1/XMININ is machine representable 2289 2290 Approximate values for some important machines are: 2291 2292 beta maxexp XBIG 2293 2294 CRAY-1 (S.P.) 2 8191 966.961 2295 Cyber 180/855 2296 under NOS (S.P.) 2 1070 177.803 2297 IEEE (IBM/XT, 2298 SUN, etc.) (S.P.) 2 128 35.040 2299 IEEE (IBM/XT, 2300 SUN, etc.) (D.P.) 2 1024 171.624 2301 IBM 3033 (D.P.) 16 63 57.574 2302 VAX D-Format (D.P.) 2 127 34.844 2303 VAX G-Format (D.P.) 2 1023 171.489 2304 2305 XINF EPS XMININ 2306 2307 CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466 2308 Cyber 180/855 2309 under NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294 2310 IEEE (IBM/XT, 2311 SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.18E-38 2312 IEEE (IBM/XT, 2313 SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.23D-308 2314 IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76 2315 VAX D-Format (D.P.) 1.70D+38 1.39D-17 5.88D-39 2316 VAX G-Format (D.P.) 8.98D+307 1.11D-16 1.12D-308 2317 2318 ******************************************************************* 2319 2320 ---------------------------------------------------------------------- 2321 Machine dependent parameters 2322 ---------------------------------------------------------------------- 2323 */ 2324 enum double xbg = 171.624; 2325 /* ML_POSINF == const double xinf = 1.79e308;*/ 2326 /* DBL_EPSILON = const double eps = 2.22e-16;*/ 2327 /* DBL_MIN == const double xminin = 2.23e-308;*/ 2328 2329 /*---------------------------------------------------------------------- 2330 Numerator and denominator coefficients for rational minimax 2331 approximation over (1,2). 2332 ----------------------------------------------------------------------*/ 2333 static immutable double[8] p = [-1.71618513886549492533811, 2334 24.7656508055759199108314, -379.804256470945635097577, 2335 629.331155312818442661052, 866.966202790413211295064, 2336 -31451.2729688483675254357, 2337 -36144.4134186911729807069, 2338 66456.1438202405440627855]; 2339 static immutable double[8] q = [-30.8402300119738975254353, 2340 315.350626979604161529144, -1015.15636749021914166146, 2341 -3107.77167157231109440444, 22538.1184209801510330112, 2342 4755.84627752788110767815, 2343 -134659.959864969306392456, 2344 -115132.259675553483497211]; 2345 /*---------------------------------------------------------------------- 2346 Coefficients for minimax approximation over (12, INF). 2347 ----------------------------------------------------------------------*/ 2348 static immutable double[7] c = [-.001910444077728, 8.4171387781295e-4, 2349 -5.952379913043012e-4, 7.93650793500350248e-4, -.002777777777777681622553, 2350 .08333333333333333331554247, .0057083835261]; 2351 /* Local variables */ 2352 int i, n; 2353 int parity; /*logical*/ 2354 double fact, xden, xnum, y, z, yi, res, sum, ysq; 2355 parity = (0); 2356 fact = 1.; 2357 n = 0; 2358 y = x; 2359 if (y <= 0.) 2360 { 2361 /* ------------------------------------------------------------- 2362 Argument is negative 2363 ------------------------------------------------------------- */ 2364 y = -x; 2365 yi = trunc(y); 2366 res = y - yi; 2367 if (res != 0.) 2368 { 2369 if (yi != trunc(yi * .5) * 2.) 2370 parity = (1); 2371 fact = -std.math.PI / sinPi(res); 2372 y += 1.; 2373 } 2374 else 2375 { 2376 return (double.infinity); 2377 } 2378 } 2379 /* ----------------------------------------------------------------- 2380 Argument is positive 2381 -----------------------------------------------------------------*/ 2382 if (y < double.epsilon) 2383 { 2384 /* -------------------------------------------------------------- 2385 Argument < EPS 2386 -------------------------------------------------------------- */ 2387 if (y >= double.min_normal) 2388 { 2389 res = 1. / y; 2390 } 2391 else 2392 { 2393 return (double.infinity); 2394 } 2395 } 2396 else if (y < 12.) 2397 { 2398 yi = y; 2399 if (y < 1.) 2400 { 2401 /* --------------------------------------------------------- 2402 EPS < argument < 1 2403 --------------------------------------------------------- */ 2404 z = y; 2405 y += 1.; 2406 } 2407 else 2408 { 2409 /* ----------------------------------------------------------- 2410 1 <= argument < 12, reduce argument if necessary 2411 ----------------------------------------------------------- */ 2412 n = cast(int) y - 1; 2413 y -= cast(double) n; 2414 z = y - 1.; 2415 } 2416 /* --------------------------------------------------------- 2417 Evaluate approximation for 1. < argument < 2. 2418 ---------------------------------------------------------*/ 2419 xnum = 0; 2420 xden = 1.; 2421 for (i = 0; i < 8; ++i) 2422 { 2423 xnum = (xnum + p[i]) * z; 2424 xden = xden * z + q[i]; 2425 } 2426 res = xnum / xden + 1.; 2427 if (yi < y) 2428 { 2429 /* -------------------------------------------------------- 2430 Adjust result for case 0. < argument < 1. 2431 -------------------------------------------------------- */ 2432 res /= yi; 2433 } 2434 else if (yi > y) 2435 { 2436 /* ---------------------------------------------------------- 2437 Adjust result for case 2. < argument < 12. 2438 ---------------------------------------------------------- */ 2439 for (i = 0; i < n; ++i) 2440 { 2441 res *= y; 2442 y += 1.; 2443 } 2444 } 2445 } 2446 else 2447 { 2448 /* ------------------------------------------------------------- 2449 Evaluate for argument >= 12., 2450 ------------------------------------------------------------- */ 2451 if (y <= xbg) 2452 { 2453 ysq = y * y; 2454 sum = c[6]; 2455 for (i = 0; i < 6; ++i) 2456 { 2457 sum = sum / ysq + c[i]; 2458 } 2459 sum = sum / y - y + sqrtpi; 2460 sum += (y - .5) * log(y); 2461 res = exp(sum); 2462 } 2463 else 2464 { 2465 return double.infinity; 2466 } 2467 } 2468 /* ---------------------------------------------------------------------- 2469 Final adjustments and return 2470 ----------------------------------------------------------------------*/ 2471 if (parity) 2472 res = -res; 2473 if (fact != 1.) 2474 res = fact / res; 2475 return res; 2476 } 2477 2478 /// sin of pix: 2479 T sinPi(T)(T x) if (isFloatingPoint!T) 2480 { 2481 if (x < 0) 2482 return -sinPi(-x); 2483 // sin of pix: 2484 bool invert; 2485 if (x < cast(T)0.5) 2486 return sin(T(std.math.PI) * x); 2487 if (x < 1) 2488 { 2489 invert = true; 2490 x = -x; 2491 } 2492 else invert = false; 2493 T rem = floor(x); 2494 if (cast(int)rem & 1) 2495 invert = !invert; 2496 rem = x - rem; 2497 if (rem > cast(T)0.5) 2498 rem = 1 - rem; 2499 if (rem == cast(T)0.5) 2500 return invert ? -1 : 1; 2501 rem = sin(T(std.math.PI) * rem); 2502 return invert ? T(-rem) : rem; 2503 } 2504 2505 unittest 2506 { 2507 import std.math : approxEqual; 2508 2509 assert(sinPi(0.0) == 0); 2510 assert(sinPi(1.0) == 0); 2511 assert(sinPi(+1.0 / 6).approxEqual(+0.5)); 2512 assert(sinPi(-1.0 / 6).approxEqual(-0.5)); 2513 assert(sinPi(+5.0 / 6).approxEqual(+0.5)); 2514 assert(sinPi(-5.0 / 6).approxEqual(-0.5)); 2515 assert(sinPi(+2.0) == 0); 2516 assert(sinPi(-2.0) == 0); 2517 assert(sinPi(+0.5) == +1); 2518 assert(sinPi(-0.5) == -1); 2519 } 2520 2521 /// cos of pix: 2522 T cosPi(T)(T x) if (isFloatingPoint!T) 2523 { 2524 // cos of pix: 2525 bool invert = false; 2526 if (fabs(x) < 0.5) 2527 return cos(T(std.math.PI) * x); 2528 if (x < cast(T)1) 2529 { 2530 x = -x; 2531 } 2532 T rem = floor(x); 2533 if (cast(int)rem & 1) 2534 invert = !invert; 2535 rem = x - rem; 2536 if (rem > cast(T)0.5) 2537 { 2538 rem = 1 - rem; 2539 invert = !invert; 2540 } 2541 if (rem == cast(T)0.5) 2542 return 0; 2543 rem = cos(T(std.math.PI) * rem); 2544 return invert ? T(-rem) : rem; 2545 } 2546 2547 unittest 2548 { 2549 import std.math : approxEqual; 2550 2551 assert(cosPi(0.0) == 1); 2552 assert(cosPi(1.0) == -1); 2553 assert(cosPi(+1.0 / 3).approxEqual(+0.5)); 2554 assert(cosPi(-1.0 / 3).approxEqual(+0.5)); 2555 assert(cosPi(+2.0 / 3).approxEqual(-0.5)); 2556 assert(cosPi(-2.0 / 3).approxEqual(-0.5)); 2557 assert(cosPi(2.0) == 1); 2558 assert(cosPi(0.5) == 0); 2559 }