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