1 module rmathd.toms708; 2 3 public import rmathd.common; 4 public import rmathd.gamma; 5 6 static auto gamln1(T)(T a) 7 { 8 /* ----------------------------------------------------------------------- */ 9 /* EVALUATION OF LN(GAMMA(1 + A)) FOR -0.2 <= A <= 1.25 */ 10 /* ----------------------------------------------------------------------- */ 11 12 T w; 13 if (a < 0.6) { 14 static T p0 = .577215664901533; 15 static T p1 = .844203922187225; 16 static T p2 = -.168860593646662; 17 static T p3 = -.780427615533591; 18 static T p4 = -.402055799310489; 19 static T p5 = -.0673562214325671; 20 static T p6 = -.00271935708322958; 21 static T q1 = 2.88743195473681; 22 static T q2 = 3.12755088914843; 23 static T q3 = 1.56875193295039; 24 static T q4 = .361951990101499; 25 static T q5 = .0325038868253937; 26 static T q6 = 6.67465618796164e-4; 27 w = ((((((p6 * a + p5)* a + p4)* a + p3)* a + p2)* a + p1)* a + p0) / 28 ((((((q6 * a + q5)* a + q4)* a + q3)* a + q2)* a + q1)* a + 1.); 29 return -(a) * w; 30 } 31 else { /* 0.6 <= a <= 1.25 */ 32 static T r0 = .422784335098467; 33 static T r1 = .848044614534529; 34 static T r2 = .565221050691933; 35 static T r3 = .156513060486551; 36 static T r4 = .017050248402265; 37 static T r5 = 4.97958207639485e-4; 38 static T s1 = 1.24313399877507; 39 static T s2 = .548042109832463; 40 static T s3 = .10155218743983; 41 static T s4 = .00713309612391; 42 static T s5 = 1.16165475989616e-4; 43 T x = a - 0.5 - 0.5; 44 w = (((((r5 * x + r4) * x + r3) * x + r2) * x + r1) * x + r0) / 45 (((((s5 * x + s4) * x + s3) * x + s2) * x + s1) * x + 1.); 46 return x * w; 47 } 48 } /* gamln1 */ 49 50 51 static auto gamln(T)(T a) 52 { 53 /* ----------------------------------------------------------------------- 54 * Evaluation of ln(gamma(a)) for positive a 55 * ----------------------------------------------------------------------- */ 56 /* Written by Alfred H. Morris */ 57 /* Naval Surface Warfare Center */ 58 /* Dahlgren, Virginia */ 59 /* ----------------------------------------------------------------------- */ 60 61 static T d = .418938533204673;/* d == 0.5*(LN(2*PI) - 1) */ 62 63 static T c0 = .0833333333333333; 64 static T c1 = -.00277777777760991; 65 static T c2 = 7.9365066682539e-4; 66 static T c3 = -5.9520293135187e-4; 67 static T c4 = 8.37308034031215e-4; 68 static T c5 = -.00165322962780713; 69 70 if (a <= 0.8) 71 return gamln1!T(a) - log(a); /* ln(G(a+1)) - ln(a) == ln(G(a+1)/a) = ln(G(a)) */ 72 else if (a <= 2.25) 73 return gamln1!T(a - 0.5 - 0.5); 74 75 else if (a < 10.) { 76 int i, n = cast(int)(a - 1.25); 77 T t = a; 78 T w = 1.; 79 for (i = 1; i <= n; ++i) { 80 t += -1.; 81 w *= t; 82 } 83 return gamln1!T(t - 1.) + log(w); 84 } 85 else { /* a >= 10 */ 86 T t = 1. / (a * a); 87 T w = (((((c5 * t + c4) * t + c3) * t + c2) * t + c1) * t + c0) / a; 88 return d + w + (a - 0.5) * (log(a) - 1.); 89 } 90 } /* gamln */ 91 92 93 static auto bcorr(T)(T a0, T b0) 94 { 95 /* ----------------------------------------------------------------------- */ 96 97 /* EVALUATION OF DEL(A0) + DEL(B0) - DEL(A0 + B0) WHERE */ 98 /* LN(GAMMA(A)) = (A - 0.5)*LN(A) - A + 0.5*LN(2*PI) + DEL(A). */ 99 /* IT IS ASSUMED THAT A0 >= 8 AND B0 >= 8. */ 100 101 /* ----------------------------------------------------------------------- */ 102 /* Initialized data */ 103 104 static T c0 = .0833333333333333; 105 static T c1 = -.00277777777760991; 106 static T c2 = 7.9365066682539e-4; 107 static T c3 = -5.9520293135187e-4; 108 static T c4 = 8.37308034031215e-4; 109 static T c5 = -.00165322962780713; 110 111 /* System generated locals */ 112 T ret_val, r1; 113 114 /* Local variables */ 115 T a, b, c, h, t, w, x, s3, s5, x2, s7, s9, s11; 116 /* ------------------------ */ 117 a = min!(T)(a0, b0); 118 b = max!(T)(a0, b0); 119 120 h = a / b; 121 c = h / (h + 1.); 122 x = 1. / (h + 1.); 123 x2 = x * x; 124 125 /* SET SN = (1 - X^N)/(1 - X) */ 126 127 s3 = x + x2 + 1.; 128 s5 = x + x2 * s3 + 1.; 129 s7 = x + x2 * s5 + 1.; 130 s9 = x + x2 * s7 + 1.; 131 s11 = x + x2 * s9 + 1.; 132 133 /* SET W = DEL(B) - DEL(A + B) */ 134 135 /* Computing 2nd power */ 136 r1 = 1. / b; 137 t = r1 * r1; 138 w = ((((c5 * s11 * t + c4 * s9) * t + c3 * s7) * t + c2 * s5) * t + c1 * 139 s3) * t + c0; 140 w *= c / b; 141 142 /* COMPUTE DEL(A) + W */ 143 144 /* Computing 2nd power */ 145 r1 = 1. / a; 146 t = r1 * r1; 147 ret_val = (((((c5 * t + c4) * t + c3) * t + c2) * t + c1) * t + c0) / a + 148 w; 149 return ret_val; 150 } /* bcorr */ 151 152 153 static auto psi(T)(T x) 154 { 155 /* --------------------------------------------------------------------- 156 157 * Evaluation of the Digamma function psi(x) 158 159 * ----------- 160 161 * Psi(xx) is assigned the value 0 when the digamma function cannot 162 * be computed. 163 164 * The main computation involves evaluation of rational Chebyshev 165 * approximations published in Math. Comp. 27, 123-127(1973) by 166 * Cody, Strecok and Thacher. */ 167 168 /* --------------------------------------------------------------------- */ 169 /* Psi was written at Argonne National Laboratory for the FUNPACK */ 170 /* package of special function subroutines. Psi was modified by */ 171 /* A.H. Morris (NSWC). */ 172 /* --------------------------------------------------------------------- */ 173 174 static T piov4 = .785398163397448; /* == pi / 4 */ 175 /* dx0 = zero of psi() to extended precision : */ 176 static T dx0 = 1.461632144968362341262659542325721325; 177 178 /* --------------------------------------------------------------------- */ 179 /* COEFFICIENTS FOR RATIONAL APPROXIMATION OF */ 180 /* PSI(X) / (X - X0), 0.5 <= X <= 3. */ 181 static T[7] p1 = [.0089538502298197,4.77762828042627, 182 142.441585084029,1186.45200713425,3633.51846806499, 183 4138.10161269013,1305.60269827897]; 184 static T[6] q1 = [44.8452573429826,520.752771467162, 185 2210.0079924783,3641.27349079381,1908.310765963, 186 6.91091682714533e-6]; 187 /* --------------------------------------------------------------------- */ 188 189 190 /* --------------------------------------------------------------------- */ 191 /* COEFFICIENTS FOR RATIONAL APPROXIMATION OF */ 192 /* PSI(X) - LN(X) + 1 / (2*X), X > 3. */ 193 194 static T[4] p2 = [-2.12940445131011,-7.01677227766759, 195 -4.48616543918019,-.648157123766197]; 196 static T[4] q2 = [32.2703493791143,89.2920700481861, 197 54.6117738103215,7.77788548522962]; 198 /* --------------------------------------------------------------------- */ 199 200 int i, m, n, nq; 201 T d2; 202 T w, z; 203 T den, aug, sgn, xmx0, xmax1, upper, xsmall; 204 205 /* --------------------------------------------------------------------- */ 206 207 208 /* MACHINE DEPENDENT CONSTANTS ... */ 209 210 /* --------------------------------------------------------------------- */ 211 /* XMAX1 = THE SMALLEST POSITIVE FLOATING POINT CONSTANT 212 WITH ENTIRELY INT REPRESENTATION. ALSO USED 213 AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE 214 ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH 215 PSI MAY BE REPRESENTED AS LOG(X). 216 * Originally: xmax1 = amin1(ipmpar(3), 1./spmpar(1)) */ 217 xmax1 = cast(T) INT_MAX; 218 d2 = 0.5/d1mach!T(3); /*= 0.5 / (0.5 * DBL_EPS) = 1/DBL_EPSILON = 2^52 */ 219 if(xmax1 > d2) xmax1 = d2; 220 221 /* --------------------------------------------------------------------- */ 222 /* XSMALL = ABSOLUTE ARGUMENT BELOW WHICH PI*COTAN(PI*X) */ 223 /* MAY BE REPRESENTED BY 1/X. */ 224 xsmall = 1e-9; 225 /* --------------------------------------------------------------------- */ 226 aug = 0.; 227 if (x < 0.5) { 228 /* --------------------------------------------------------------------- */ 229 /* X < 0.5, USE REFLECTION FORMULA */ 230 /* PSI(1-X) = PSI(X) + PI * COTAN(PI*X) */ 231 /* --------------------------------------------------------------------- */ 232 if (fabs(x) <= xsmall) { 233 234 if (x == 0.) { 235 goto L_err; 236 } 237 /* --------------------------------------------------------------------- */ 238 /* 0 < |X| <= XSMALL. USE 1/X AS A SUBSTITUTE */ 239 /* FOR PI*COTAN(PI*X) */ 240 /* --------------------------------------------------------------------- */ 241 aug = -1. / x; 242 } else { /* |x| > xsmall */ 243 /* --------------------------------------------------------------------- */ 244 /* REDUCTION OF ARGUMENT FOR COTAN */ 245 /* --------------------------------------------------------------------- */ 246 /* L100: */ 247 w = -x; 248 sgn = piov4; 249 if (w <= 0.) { 250 w = -w; 251 sgn = -sgn; 252 } 253 /* --------------------------------------------------------------------- */ 254 /* MAKE AN ERROR EXIT IF |X| >= XMAX1 */ 255 /* --------------------------------------------------------------------- */ 256 if (w >= xmax1) { 257 goto L_err; 258 } 259 nq = cast(int) w; 260 w -= cast(T) nq; 261 nq = cast(int) (w * 4.); 262 w = (w - cast(T) nq * 0.25) * 4.; 263 /* --------------------------------------------------------------------- */ 264 /* W IS NOW RELATED TO THE FRACTIONAL PART OF 4. * X. */ 265 /* ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST */ 266 /* QUADRANT AND DETERMINE SIGN */ 267 /* --------------------------------------------------------------------- */ 268 n = nq / 2; 269 if (n + n != nq) { 270 w = 1. - w; 271 } 272 z = piov4 * w; 273 m = n / 2; 274 if (m + m != n) { 275 sgn = -sgn; 276 } 277 /* --------------------------------------------------------------------- */ 278 /* DETERMINE FINAL VALUE FOR -PI*COTAN(PI*X) */ 279 /* --------------------------------------------------------------------- */ 280 n = (nq + 1) / 2; 281 m = n / 2; 282 m += m; 283 if (m == n) { 284 /* --------------------------------------------------------------------- */ 285 /* CHECK FOR SINGULARITY */ 286 /* --------------------------------------------------------------------- */ 287 if (z == 0.) { 288 goto L_err; 289 } 290 /* --------------------------------------------------------------------- */ 291 /* USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND */ 292 /* SIN/COS AS A SUBSTITUTE FOR TAN */ 293 /* --------------------------------------------------------------------- */ 294 aug = sgn*(cos(z) / sin(z) * 4.); 295 296 } else { /* L140: */ 297 aug = sgn*(sin(z) / cos(z) * 4.); 298 } 299 } 300 301 x = 1. - x; 302 303 } 304 /* L200: */ 305 if (x <= 3.) { 306 /* --------------------------------------------------------------------- */ 307 /* 0.5 <= X <= 3. */ 308 /* --------------------------------------------------------------------- */ 309 den = x; 310 upper = p1[0] * x; 311 312 for (i = 1; i <= 5; ++i) { 313 den = (den + q1[i - 1]) * x; 314 upper = (upper + p1[i]) * x; 315 } 316 317 den = (upper + p1[6]) / (den + q1[5]); 318 xmx0 = x - dx0; 319 return den * xmx0 + aug; 320 } 321 322 /* --------------------------------------------------------------------- */ 323 /* IF X >= XMAX1, PSI = LN(X) */ 324 /* --------------------------------------------------------------------- */ 325 if (x < xmax1) { 326 /* --------------------------------------------------------------------- */ 327 /* 3. < X < XMAX1 */ 328 /* --------------------------------------------------------------------- */ 329 w = 1. / (x * x); 330 den = w; 331 upper = p2[0] * w; 332 333 for (i = 1; i <= 3; ++i) { 334 den = (den + q2[i - 1]) * w; 335 upper = (upper + p2[i]) * w; 336 } 337 338 aug = upper / (den + q2[3]) - 0.5 / x + aug; 339 } 340 return aug + log(x); 341 342 /* --------------------------------------------------------------------- */ 343 /* ERROR RETURN */ 344 /* --------------------------------------------------------------------- */ 345 L_err: 346 return 0.; 347 } /* psi */ 348 349 350 static auto gam1(T)(T a) 351 { 352 /* ------------------------------------------------------------------ */ 353 /* COMPUTATION OF 1/GAMMA(A+1) - 1 FOR -0.5 <= A <= 1.5 */ 354 /* ------------------------------------------------------------------ */ 355 356 T d, t, w, bot, top; 357 358 t = a; 359 d = a - 0.5; 360 // t := if(a > 1/2) a-1 else a 361 if (d > 0.) 362 t = d - 0.5; 363 if (t < 0.) { /* L30: */ 364 static T[9] r = [-.422784335098468,-.771330383816272, 365 -.244757765222226,.118378989872749,9.30357293360349e-4, 366 -.0118290993445146,.00223047661158249,2.66505979058923e-4, 367 -1.32674909766242e-4]; 368 static T s1 = .273076135303957, s2 = .0559398236957378; 369 370 top = (((((((r[8] * t + r[7]) * t + r[6]) * t + r[5]) * t + r[4] 371 ) * t + r[3]) * t + r[2]) * t + r[1]) * t + r[0]; 372 bot = (s2 * t + s1) * t + 1.; 373 w = top/bot; 374 //R_ifDEBUG_printf(" gam1(a = %.15g): t < 0: w=%.15g\n", a, w); 375 if (d > 0.) 376 return t * w / a; 377 else 378 return a * (w + 0.5 + 0.5); 379 380 } else if (t == 0) { // L10: a in {0, 1} 381 return 0.; 382 383 } else { /* t > 0; L20: */ 384 static T[7] p = [.577215664901533,-.409078193005776, 385 -.230975380857675,.0597275330452234,.0076696818164949, 386 -.00514889771323592,5.89597428611429e-4]; 387 static T[5] q = [1.,.427569613095214,.158451672430138, .0261132021441447,.00423244297896961]; 388 389 top = (((((p[6] * t + p[5]) * t + p[4]) * t + p[3]) * t + p[2] 390 ) * t + p[1]) * t + p[0]; 391 bot = (((q[4] * t + q[3]) * t + q[2]) * t + q[1]) * t + 1.; 392 w = top/bot; 393 //R_ifDEBUG_printf(" gam1(a = %.15g): t > 0: (is a < 1.5 ?) w=%.15g\n", a, w); 394 if (d > 0.) /* L21: */ 395 return t / a * (w - 0.5 - 0.5); 396 else 397 return a * w; 398 } 399 } /* gam1 */ 400 401 402 static auto exparg(T)(int l) 403 { 404 /* -------------------------------------------------------------------- 405 * If l = 0 then exparg(l) = The largest positive W for which 406 * exp(W) can be computed. With 0.99999 fuzz ==> exparg(0) = 709.7756 nowadays 407 408 * if l = 1 (nonzero) then exparg(l) = the largest negative W for 409 * which the computed value of exp(W) is nonzero. 410 * With 0.99999 fuzz ==> exparg(1) = -709.0825 nowadays 411 412 * Note... only an approximate value for exparg(L) is needed. 413 * -------------------------------------------------------------------- */ 414 415 static const(T) lnb = .69314718055995; 416 int m = (l == 0) ? i1mach(16) : i1mach(15) - 1; 417 418 return m * lnb * .99999; 419 } /* exparg */ 420 421 422 static auto esum(T)(int mu, T x, int give_log) 423 { 424 /* ----------------------------------------------------------------------- */ 425 /* EVALUATION OF EXP(MU + X) */ 426 /* ----------------------------------------------------------------------- */ 427 428 if(give_log) 429 return x + cast(T) mu; 430 431 // else : 432 T w; 433 if(x > 0.) { /* L10: */ 434 if(mu > 0) 435 return exp(cast(T)mu)*exp(x); 436 w = mu + x; 437 if(w < 0.) 438 return exp(cast(T)mu)*exp(x); 439 } 440 else { /* x <= 0 */ 441 if(mu < 0) 442 return exp(cast(T)mu)*exp(x); 443 w = mu + x; 444 if(w > 0.) 445 return exp(cast(T)mu)*exp(x); 446 } 447 return exp(w); 448 449 } /* esum */ 450 451 452 auto rexpm1(T)(T x) 453 { 454 /* ----------------------------------------------------------------------- */ 455 /* EVALUATION OF THE FUNCTION EXP(X) - 1 */ 456 /* ----------------------------------------------------------------------- */ 457 458 static T p1 = 9.14041914819518e-10; 459 static T p2 = .0238082361044469; 460 static T q1 = -.499999999085958; 461 static T q2 = .107141568980644; 462 static T q3 = -.0119041179760821; 463 static T q4 = 5.95130811860248e-4; 464 465 if (fabs(x) <= 0.15) { 466 return x * (((p2 * x + p1) * x + 1.) / 467 ((((q4 * x + q3) * x + q2) * x + q1) * x + 1.)); 468 } 469 else { /* |x| > 0.15 : */ 470 T w = exp(x); 471 if (x > 0.) 472 return w * (0.5 - 1. / w + 0.5); 473 else 474 return w - 0.5 - 0.5; 475 } 476 477 } /* rexpm1 */ 478 479 static auto alnrel(T)(T a) 480 { 481 /* ----------------------------------------------------------------------- 482 * Evaluation of the function ln(1 + a) 483 * ----------------------------------------------------------------------- */ 484 485 if (fabs(a) > 0.375) 486 return log(1. + a); 487 // else : |a| <= 0.375 488 static T p1 = -1.29418923021993, 489 p2 = .405303492862024, 490 p3 = -.0178874546012214, 491 q1 = -1.62752256355323, 492 q2 = .747811014037616, 493 q3 = -.0845104217945565; 494 T t = a / (a + 2.), 495 t2 = t * t, 496 w = (((p3 * t2 + p2) * t2 + p1) * t2 + 1.) / 497 (((q3 * t2 + q2) * t2 + q1) * t2 + 1.); 498 return t * 2. * w; 499 500 } /* alnrel */ 501 502 503 static auto rlog1(T)(T x) 504 { 505 /* ----------------------------------------------------------------------- 506 * Evaluation of the function x - ln(1 + x) 507 * ----------------------------------------------------------------------- */ 508 509 static T a = .0566749439387324; 510 static T b = .0456512608815524; 511 static T p0 = .333333333333333; 512 static T p1 = -.224696413112536; 513 static T p2 = .00620886815375787; 514 static T q1 = -1.27408923933623; 515 static T q2 = .354508718369557; 516 517 T h, r, t, w, w1; 518 if (x < -0.39 || x > 0.57) { /* direct evaluation */ 519 w = x + 0.5 + 0.5; 520 return x - log(w); 521 } 522 /* else */ 523 if (x < -0.18) { /* L10: */ 524 h = x + .3; 525 h /= .7; 526 w1 = a - h * .3; 527 } else if (x > 0.18) { /* L20: */ 528 h = x * .75 - .25; 529 w1 = b + h / 3.; 530 } else { /* Argument Reduction */ 531 h = x; 532 w1 = 0.; 533 } 534 535 /* L30: Series Expansion */ 536 537 r = h/(h + 2.); 538 t = r*r; 539 w = ((p2 * t + p1) * t + p0) / ((q2 * t + q1) * t + 1.); 540 return t * 2. * (1. / (1. - r) - r * w) + w1; 541 542 } /* rlog1 */ 543 544 545 static auto erf__(T)(T x) 546 { 547 /* ----------------------------------------------------------------------- 548 * EVALUATION OF THE REAL ERROR FUNCTION 549 * ----------------------------------------------------------------------- */ 550 551 /* Initialized data */ 552 553 static T c = .564189583547756; 554 static T[5] a = [7.7105849500132e-5,-.00133733772997339, 555 .0323076579225834,.0479137145607681,.128379167095513]; 556 static T[3] b = [.00301048631703895,.0538971687740286, 557 .375795757275549]; 558 static T[8] p = [-1.36864857382717e-7,.564195517478974, 559 7.21175825088309,43.1622272220567,152.98928504694, 560 339.320816734344,451.918953711873,300.459261020162]; 561 static T[8] q = [1.,12.7827273196294,77.0001529352295, 562 277.585444743988,638.980264465631,931.35409485061, 563 790.950925327898,300.459260956983]; 564 static T[5] r = [2.10144126479064,26.2370141675169, 565 21.3688200555087,4.6580782871847,.282094791773523]; 566 static T[4] s = [94.153775055546,187.11481179959, 567 99.0191814623914,18.0124575948747]; 568 569 /* Local variables */ 570 T t, x2, ax, bot, top; 571 572 ax = fabs(x); 573 if (ax <= 0.5) { 574 t = x * x; 575 top = (((a[0] * t + a[1]) * t + a[2]) * t + a[3]) * t + a[4] + 1.; 576 bot = ((b[0] * t + b[1]) * t + b[2]) * t + 1.; 577 578 return x * (top / bot); 579 } 580 581 // else: |x| > 0.5 582 583 if (ax <= 4.) { // |x| in (0.5, 4] 584 top = ((((((p[0] * ax + p[1]) * ax + p[2]) * ax + p[3]) * ax + p[4]) * ax 585 + p[5]) * ax + p[6]) * ax + p[7]; 586 bot = ((((((q[0] * ax + q[1]) * ax + q[2]) * ax + q[3]) * ax + q[4]) * ax 587 + q[5]) * ax + q[6]) * ax + q[7]; 588 T R = 0.5 - exp(-x * x) * top / bot + 0.5; 589 return (x < 0) ? -R : R; 590 } 591 592 // else: |x| > 4 593 594 if (ax >= 5.8) { 595 return x > 0 ? 1 : -1; 596 } 597 598 // else: 4 < |x| < 5.8 599 x2 = x * x; 600 t = 1. / x2; 601 top = (((r[0] * t + r[1]) * t + r[2]) * t + r[3]) * t + r[4]; 602 bot = (((s[0] * t + s[1]) * t + s[2]) * t + s[3]) * t + 1.; 603 t = (c - top / (x2 * bot)) / ax; 604 T R = 0.5 - exp(-x2) * t + 0.5; 605 return (x < 0) ? -R : R; 606 } /* erf */ 607 608 609 static auto erfc1(T)(int ind, T x) 610 { 611 /* ----------------------------------------------------------------------- */ 612 /* EVALUATION OF THE COMPLEMENTARY ERROR FUNCTION */ 613 614 /* ERFC1(IND,X) = ERFC(X) IF IND = 0 */ 615 /* ERFC1(IND,X) = EXP(X*X)*ERFC(X) OTHERWISE */ 616 /* ----------------------------------------------------------------------- */ 617 618 /* Initialized data */ 619 620 static T c = .564189583547756; 621 static T[5] a = [7.7105849500132e-5,-.00133733772997339, 622 .0323076579225834,.0479137145607681,.128379167095513]; 623 static T[3] b = [.00301048631703895,.0538971687740286, .375795757275549]; 624 static T[8] p = [-1.36864857382717e-7,.564195517478974, 625 7.21175825088309,43.1622272220567,152.98928504694, 626 339.320816734344,451.918953711873,300.459261020162]; 627 static T[8] q = [1.,12.7827273196294,77.0001529352295, 628 277.585444743988,638.980264465631,931.35409485061, 629 790.950925327898,300.459260956983]; 630 static T[5] r = [2.10144126479064,26.2370141675169, 631 21.3688200555087,4.6580782871847,.282094791773523]; 632 static T[4] s = [94.153775055546,187.11481179959, 633 99.0191814623914,18.0124575948747]; 634 635 T ret_val; 636 T e, t, w, bot, top; 637 638 T ax = fabs(x); 639 // |X| <= 0.5 */ 640 if (ax <= 0.5) { 641 t = x * x, 642 top = (((a[0] * t + a[1]) * t + a[2]) * t + a[3]) * t + a[4] + 1., 643 bot = ((b[0] * t + b[1]) * t + b[2]) * t + 1.; 644 ret_val = 0.5 - x * (top / bot) + 0.5; 645 if (ind != 0) { 646 ret_val = exp(t) * ret_val; 647 } 648 return ret_val; 649 } 650 // else (L10:): 0.5 < |X| <= 4 651 if (ax <= 4.) { 652 top = ((((((p[0] * ax + p[1]) * ax + p[2]) * ax + p[3]) * ax + p[4]) * ax 653 + p[5]) * ax + p[6]) * ax + p[7]; 654 bot = ((((((q[0] * ax + q[1]) * ax + q[2]) * ax + q[3]) * ax + q[4]) * ax 655 + q[5]) * ax + q[6]) * ax + q[7]; 656 ret_val = top / bot; 657 } else { // |X| > 4 658 // L20: 659 if (x <= -5.6) { 660 // L50: LIMIT VALUE FOR "LARGE" NEGATIVE X 661 ret_val = 2.; 662 if (ind != 0) { 663 ret_val = exp(x * x) * 2.; 664 } 665 return ret_val; 666 } 667 if (ind == 0 && (x > 100. || x * x > -exparg!T(1))) { 668 // LIMIT VALUE FOR LARGE POSITIVE X WHEN IND = 0 669 // L60: 670 return 0.; 671 } 672 673 // L30: 674 t = 1. / (x * x); 675 top = (((r[0] * t + r[1]) * t + r[2]) * t + r[3]) * t + r[4]; 676 bot = (((s[0] * t + s[1]) * t + s[2]) * t + s[3]) * t + 1.; 677 ret_val = (c - t * top / bot) / ax; 678 } 679 680 // L40: FINAL ASSEMBLY 681 if (ind != 0) { 682 if (x < 0.) 683 ret_val = exp(x * x) * 2. - ret_val; 684 } else { 685 // L41: ind == 0 : 686 w = x * x; 687 t = w; 688 e = w - t; 689 ret_val = (0.5 - e + 0.5) * exp(-t) * ret_val; 690 if (x < 0.) 691 ret_val = 2. - ret_val; 692 } 693 return ret_val; 694 695 } /* erfc1 */ 696 697 698 static auto basym(T)(T a, T b, T lambda, T eps, int log_p) 699 { 700 /* ----------------------------------------------------------------------- */ 701 /* ASYMPTOTIC EXPANSION FOR I_x(A,B) FOR LARGE A AND B. */ 702 /* LAMBDA = (A + B)*Y - B AND EPS IS THE TOLERANCE USED. */ 703 /* IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT */ 704 /* A AND B ARE GREATER THAN OR EQUAL TO 15. */ 705 /* ----------------------------------------------------------------------- */ 706 707 708 /* ------------------------ */ 709 /* ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP */ 710 /* ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN. */ 711 enum num_IT = 20; 712 /* THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1. */ 713 714 static const(T) e0 = 1.12837916709551;/* e0 == 2/sqrt(pi) */ 715 static const(T) e1 = .353553390593274;/* e1 == 2^(-3/2) */ 716 static const(T) ln_e0 = 0.120782237635245; /* == ln(e0) */ 717 718 T[num_IT + 1] a0, b0, c, d; 719 720 T f = a * rlog1!T(-lambda/a) + b * rlog1!T(lambda/b), t; 721 if(log_p) 722 t = -f; 723 else { 724 t = exp(-f); 725 if (t == 0.) { 726 return 0; /* once underflow, always underflow .. */ 727 } 728 } 729 T z0 = sqrt(f), 730 z = z0 / e1 * 0.5, 731 z2 = f + f, 732 h, r0, r1, w0; 733 734 if (a < b) { 735 h = a / b; 736 r0 = 1. / (h + 1.); 737 r1 = (b - a) / b; 738 w0 = 1. / sqrt(a * (h + 1.)); 739 } else { 740 h = b / a; 741 r0 = 1. / (h + 1.); 742 r1 = (b - a) / a; 743 w0 = 1. / sqrt(b * (h + 1.)); 744 } 745 746 a0[0] = r1 * .66666666666666663; 747 c[0] = a0[0] * -0.5; 748 d[0] = -c[0]; 749 T j0 = 0.5 / e0 * erfc1!T(1, z0), 750 j1 = e1, 751 sum = j0 + d[0] * w0 * j1; 752 753 T s = 1., 754 h2 = h * h, 755 hn = 1., 756 w = w0, 757 znm1 = z, 758 zn = z2; 759 for (int n = 2; n <= num_IT; n += 2) { 760 hn *= h2; 761 a0[n - 1] = r0 * 2. * (h * hn + 1.) / (n + 2.); 762 int np1 = n + 1; 763 s += hn; 764 a0[np1 - 1] = r1 * 2. * s / (n + 3.); 765 766 for (int i = n; i <= np1; ++i) { 767 T r = (i + 1.) * -0.5; 768 b0[0] = r * a0[0]; 769 for (int m = 2; m <= i; ++m) { 770 T bsum = 0.; 771 for (int j = 1; j <= m-1; ++j) { 772 int mmj = m - j; 773 bsum += (j * r - mmj) * a0[j - 1] * b0[mmj - 1]; 774 } 775 b0[m - 1] = r * a0[m - 1] + bsum / m; 776 } 777 c[i - 1] = b0[i - 1] / (i + 1.); 778 779 T dsum = 0.; 780 for (int j = 1; j <= i-1; ++j) { 781 dsum += d[i - j - 1] * c[j - 1]; 782 } 783 d[i - 1] = -(dsum + c[i - 1]); 784 } 785 786 j0 = e1 * znm1 + (n - 1.) * j0; 787 j1 = e1 * zn + n * j1; 788 znm1 = z2 * znm1; 789 zn = z2 * zn; 790 w *= w0; 791 T t0 = d[n - 1] * w * j0; 792 w *= w0; 793 T t1 = d[np1 - 1] * w * j1; 794 sum += t0 + t1; 795 if (fabs(t0) + fabs(t1) <= eps * sum) { 796 break; 797 } 798 } 799 800 if(log_p) 801 return ln_e0 + t - bcorr!T(a, b) + log(sum); 802 else { 803 T u = exp(-bcorr!T(a, b)); 804 return e0 * t * u * sum; 805 } 806 807 } /* basym_ */ 808 809 810 // called only from bgrat() , as q_r = grat_r(b, z, log_r, eps) : 811 static auto grat_r(T)(T a, T x, T log_r, T eps) 812 { 813 /* ----------------------------------------------------------------------- 814 * Scaled complement of incomplete gamma ratio function 815 * grat_r(a,x,r) := Q(a,x) / r 816 * where 817 * Q(a,x) = pgamma(x,a, lower.tail=FALSE) 818 * and r = e^(-x)* x^a / Gamma(a) == exp(log_r) 819 * 820 * It is assumed that a <= 1. eps is the tolerance to be used. 821 * ----------------------------------------------------------------------- */ 822 823 if (a * x == 0.) { /* L130: */ 824 if (x <= a) { 825 /* L100: */ return exp(-log_r); 826 } else { 827 /* L110:*/ return 0.; 828 } 829 } else if (a == 0.5) { // e.g. when called from pt() 830 /* L120: */ 831 if (x < 0.25) { 832 T p = erf__!T(sqrt(x)); 833 //R_ifDEBUG_printf(" grat_r(a=%g, x=%g ..)): a=1/2 --> p=erf__(.)= %g\n", a, x, p); 834 return (0.5 - p + 0.5)*exp(-log_r); 835 836 } else { // 2013-02-27: improvement for "large" x: direct computation of q/r: 837 T sx = sqrt(x), 838 q_r = erfc1!T(1, sx)/sx * M_SQRT_PI; 839 //R_ifDEBUG_printf(" grat_r(a=%g, x=%g ..)): a=1/2 --> q_r=erfc1(..)/r= %g\n", a,x, q_r); 840 return q_r; 841 } 842 843 } else if (x < 1.1) { /* L10: Taylor series for P(a,x)/x^a */ 844 845 T an = 3., 846 c = x, 847 sum = x / (a + 3.), 848 tol = eps * 0.1 / (a + 1.), t; 849 do { 850 an += 1.; 851 c *= -(x / an); 852 t = c / (a + an); 853 sum += t; 854 } while (fabs(t) > tol); 855 856 //R_ifDEBUG_printf(" grat_r(a=%g, x=%g, log_r=%g): sum=%g; Taylor w/ %.0f terms", a,x,log_r, sum, an-3.); 857 T j = a * x * ((sum/6. - 0.5/(a + 2.)) * x + 1./(a + 1.)), 858 z = a * log(x), 859 h = gam1!T(a), 860 g = h + 1.; 861 862 if ((x >= 0.25 && (a < x / 2.59)) || (z > -0.13394)) { 863 // L40: 864 T l = rexpm1!T(z), 865 q = ((l + 0.5 + 0.5) * j - l) * g - h; 866 if (q <= 0.) { 867 //R_ifDEBUG_printf(" => q_r= 0.\n"); 868 /* L110:*/ return 0.; 869 } else { 870 //R_ifDEBUG_printf(" => q_r=%.15g\n", q * exp(-log_r)); 871 return q * exp(-log_r); 872 } 873 874 } else { 875 T p = exp(z) * g * (0.5 - j + 0.5); 876 //R_ifDEBUG_printf(" => q_r=%.15g\n", (0.5 - p + 0.5) * exp(-log_r)); 877 return /* q/r = */ (0.5 - p + 0.5) * exp(-log_r); 878 } 879 880 } else { 881 /* L50: ---- (x >= 1.1) ---- Continued Fraction Expansion */ 882 883 T a2n_1 = 1., 884 a2n = 1., 885 b2n_1 = x, 886 b2n = x + (1. - a), 887 c = 1., am0, an0; 888 889 do { 890 a2n_1 = x * a2n + c * a2n_1; 891 b2n_1 = x * b2n + c * b2n_1; 892 am0 = a2n_1 / b2n_1; 893 c += 1.; 894 T c_a = c - a; 895 a2n = a2n_1 + c_a * a2n; 896 b2n = b2n_1 + c_a * b2n; 897 an0 = a2n/b2n; 898 } while (fabs(an0 - am0) >= eps * an0); 899 900 //R_ifDEBUG_printf(" grat_r(a=%g, x=%g, log_r=%g): Cont.frac. %.0f terms => q_r=%.15g\n", a,x, log_r, c-1., an0); 901 return /* q/r = (r * an0)/r = */ an0; 902 } 903 } /* grat_r */ 904 905 906 static auto algdiv(T)(T a, T b) 907 { 908 /* ----------------------------------------------------------------------- */ 909 910 /* COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B >= 8 */ 911 912 /* -------- */ 913 914 /* IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY */ 915 /* LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X). */ 916 917 /* ----------------------------------------------------------------------- */ 918 919 /* Initialized data */ 920 921 static T c0 = .0833333333333333; 922 static T c1 = -.00277777777760991; 923 static T c2 = 7.9365066682539e-4; 924 static T c3 = -5.9520293135187e-4; 925 static T c4 = 8.37308034031215e-4; 926 static T c5 = -.00165322962780713; 927 928 T c, d, h, t, u, v, w, x, s3, s5, x2, s7, s9, s11; 929 930 /* ------------------------ */ 931 if (a > b) { 932 h = b / a; 933 c = 1. / (h + 1.); 934 x = h / (h + 1.); 935 d = a + (b - 0.5); 936 } else { 937 h = a / b; 938 c = h / (h + 1.); 939 x = 1. / (h + 1.); 940 d = b + (a - 0.5); 941 } 942 943 /* Set s<n> = (1 - x^n)/(1 - x) : */ 944 945 x2 = x * x; 946 s3 = x + x2 + 1.; 947 s5 = x + x2 * s3 + 1.; 948 s7 = x + x2 * s5 + 1.; 949 s9 = x + x2 * s7 + 1.; 950 s11 = x + x2 * s9 + 1.; 951 952 /* w := Del(b) - Del(a + b) */ 953 954 t = 1./ (b * b); 955 w = ((((c5 * s11 * t + c4 * s9) * t + c3 * s7) * t + c2 * s5) * t + c1 * 956 s3) * t + c0; 957 w *= c / b; 958 959 /* COMBINE THE RESULTS */ 960 961 u = d * alnrel!T(a/b); 962 v = a * (log(b) - 1.); 963 if (u > v) 964 return w - v - u; 965 else 966 return w - u - v; 967 } /* algdiv */ 968 969 970 static void bgrat(T)(T a, T b, T x, T y, T* w, T eps, int* ierr, int log_w) 971 { 972 /* ----------------------------------------------------------------------- 973 * Asymptotic Expansion for I_x(a,b) when a is larger than b. 974 * Compute w := w + I_x(a,b) 975 * It is assumed a >= 15 and b <= 1. 976 * eps is the tolerance used. 977 * ierr is a variable that reports the status of the results. 978 * 979 * if(log_w), *w itself must be in log-space; 980 * compute w := w + I_x(a,b) but return *w = log(w): 981 * *w := log(exp(*w) + I_x(a,b)) = logspace_add(*w, log( I_x(a,b) )) 982 * ----------------------------------------------------------------------- */ 983 984 enum n_terms_bgrat = 30; 985 T[n_terms_bgrat] c, d; 986 T bm1 = b - 0.5 - 0.5, nu = a + bm1 * 0.5, /* nu = a + (b-1)/2 =: T, in (9.1) of * Didonato & Morris(1992), p.362 */ 987 lnx = (y > 0.375) ? log(x) : alnrel!T(-y), z = -nu * lnx; // z =: u in (9.1) of D.&M.(1992) 988 989 if (b*z == 0.) { 990 // should not happen, but does, e.g., 991 // for pbeta(1e-320, 1e-5, 0.5) i.e., _subnormal_ x, 992 // Warning ... bgrat(a=20.5, b=1e-05, x=1, y=9.99989e-321): .. 993 // MATHLIB_WARNING5( 994 // "bgrat(a=%g, b=%g, x=%g, y=%g): z=%g, b*z == 0 underflow, hence inaccurate pbeta()", a, b, x, y, z); 995 /* L_Error: THE EXPANSION CANNOT BE COMPUTED */ 996 *ierr = 1; return; 997 } 998 999 /* COMPUTATION OF THE EXPANSION */ 1000 T 1001 /* r1 = b * (gam1(b) + 1.) * exp(b * log(z)),// = b/gamma(b+1) z^b = z^b / gamma(b) 1002 * set r := exp(-z) * z^b / gamma(b) ; 1003 * gam1(b) = 1/gamma(b+1) - 1 , b in [-1/2, 3/2] */ 1004 // exp(a*lnx) underflows for large (a * lnx); e.g. large a ==> using log_r := log(r): 1005 // r = r1 * exp(a * lnx) * exp(bm1 * 0.5 * lnx); 1006 // log(r)=log(b) + log1p(gam1(b)) + b * log(z) + (a * lnx) + (bm1 * 0.5 * lnx), 1007 log_r = log(b) + log1p(gam1!T(b)) + b*log(z) + nu*lnx, 1008 // FIXME work with log_u = log(u) also when log_p=FALSE (??) 1009 // u is 'factored out' from the expansion {and multiplied back, at the end}: 1010 log_u = log_r - (algdiv!T(b, a) + b*log(nu)),// algdiv(b,a) = log(gamma(a)/gamma(a+b)) 1011 /* u = (log_p) ? log_r - u : exp(log_r-u); // =: M in (9.2) of {reference above} */ 1012 /* u = algdiv(b, a) + b * log(nu);// algdiv(b,a) = log(gamma(a)/gamma(a+b)) */ 1013 // u = (log_p) ? log_u : exp(log_u); // =: M in (9.2) of {reference above} 1014 u = exp(log_u); 1015 1016 if (log_u == -T.infinity) { 1017 //R_ifDEBUG_printf(" bgrat(*): underflow log_u = -Inf = log_r -u', log_r = %g ", log_r); 1018 /* L_Error: THE EXPANSION CANNOT BE COMPUTED */ *ierr = 2; return; 1019 } 1020 1021 // Rboolean 1022 int u_0 = (u == 0.); // underflow --> do work with log(u) == log_u ! 1023 T l = // := *w/u .. but with care: such that it also works when u underflows to 0: 1024 log_w ? ((*w == -T.infinity) ? 0. : exp(*w - log_u)) : ((*w == 0.) ? 0. : exp(log(*w) - log_u)); 1025 1026 //R_ifDEBUG_printf(" bgrat(a=%g, b=%g, x=%g, *)\n -> u=%g, l='w/u'=%g, ", a,b,x, u, l); 1027 T q_r = grat_r(b, z, log_r, eps), // = q/r of former grat1(b,z, r, &p, &q) 1028 v = 0.25 / (nu * nu), 1029 t2 = lnx * 0.25 * lnx, 1030 j = q_r, 1031 sum = j, 1032 t = 1., cn = 1., n2 = 0.; 1033 for (int n = 1; n <= n_terms_bgrat; ++n) { 1034 T bp2n = b + n2; 1035 j = (bp2n * (bp2n + 1.) * j + (z + bp2n + 1.) * t) * v; 1036 n2 += 2.; 1037 t *= t2; 1038 cn /= n2 * (n2 + 1.); 1039 int nm1 = n - 1; 1040 c[nm1] = cn; 1041 T s = 0.; 1042 if (n > 1) { 1043 T coef = b - n; 1044 for (int i = 1; i <= nm1; ++i) { 1045 s += coef * c[i - 1] * d[nm1 - i]; 1046 coef += b; 1047 } 1048 } 1049 d[nm1] = bm1 * cn + s / n; 1050 T dj = d[nm1] * j; 1051 sum += dj; 1052 if (sum <= 0.) { 1053 //R_ifDEBUG_printf(" bgrat(*): sum_n(..) <= 0; should not happen (n=%d)\n", n); 1054 /* L_Error: THE EXPANSION CANNOT BE COMPUTED */ *ierr = 3; return; 1055 } 1056 if (fabs(dj) <= eps * (sum + l)) { 1057 *ierr = 0; 1058 break; 1059 } else if(n == n_terms_bgrat) { // never? ; please notify R-core if seen: 1060 *ierr = 4; 1061 //MATHLIB_WARNING5( 1062 //"bgrat(a=%g, b=%g, x=%g) *no* convergence: NOTIFY R-core!\n dj=%g, rel.err=%g\n", a,b,x, dj, fabs(dj) /(sum + l)); 1063 } 1064 } // for(n .. n_terms..) 1065 1066 /* ADD THE RESULTS TO W */ 1067 1068 if(log_w) // *w is in log space already: 1069 *w = logspace_add!T(*w, log_u + log(sum)); 1070 else 1071 *w += (u_0 ? exp(log_u + log(sum)) : u * sum); 1072 return; 1073 } /* bgrat */ 1074 1075 1076 static auto gsumln(T)(T a, T b) 1077 { 1078 /* ----------------------------------------------------------------------- */ 1079 /* EVALUATION OF THE FUNCTION LN(GAMMA(A + B)) */ 1080 /* FOR 1 <= A <= 2 AND 1 <= B <= 2 */ 1081 /* ----------------------------------------------------------------------- */ 1082 1083 T x = a + b - 2.;/* in [0, 2] */ 1084 1085 if (x <= 0.25) 1086 return gamln1!T(x + 1.); 1087 1088 /* else */ 1089 if (x <= 1.25) 1090 return gamln1!T(x) + alnrel!T(x); 1091 /* else x > 1.25 : */ 1092 return gamln1!T(x - 1.) + log(x * (x + 1.)); 1093 1094 } /* gsumln */ 1095 1096 1097 static auto betaln(T)(T a0, T b0) 1098 { 1099 /* ----------------------------------------------------------------------- 1100 * Evaluation of the logarithm of the beta function ln(beta(a0,b0)) 1101 * ----------------------------------------------------------------------- */ 1102 1103 static T e = .918938533204673;/* e == 0.5*LN(2*PI) */ 1104 1105 T a = min!T(a0 ,b0), b = max!T(a0, b0); 1106 int n; 1107 1108 if (a < 8.) { 1109 if (a < 1.) { 1110 /* ----------------------------------------------------------------------- */ 1111 // A < 1 1112 /* ----------------------------------------------------------------------- */ 1113 if (b < 8.) 1114 return gamln!T(a) + (gamln!T(b) - gamln!T(a+b)); 1115 else 1116 return gamln!T(a) + algdiv!T(a, b); 1117 } 1118 /* else */ 1119 /* ----------------------------------------------------------------------- */ 1120 // 1 <= A < 8 1121 /* ----------------------------------------------------------------------- */ 1122 T w; 1123 if (a < 2.) { 1124 if (b <= 2.) { 1125 return gamln!T(a) + gamln!T(b) - gsumln!T(a, b); 1126 } 1127 /* else */ 1128 1129 if (b < 8.) { 1130 w = 0.; 1131 goto L40; 1132 } 1133 return gamln!T(a) + algdiv!T(a, b); 1134 } 1135 // else L30: REDUCTION OF A WHEN B <= 1000 1136 1137 if (b <= 1e3) { 1138 n = cast(int)(a - 1.); 1139 w = 1.; 1140 for (int i = 1; i <= n; ++i) { 1141 a += -1.; 1142 T h = a / b; 1143 w *= h / (h + 1.); 1144 } 1145 w = log(w); 1146 1147 if (b >= 8.) 1148 return w + gamln!T(a) + algdiv!T(a, b); 1149 1150 // else 1151 L40: 1152 // 1 < A <= B < 8 : reduction of B 1153 n = cast(int)(b - 1.); 1154 T z = 1.; 1155 for (int i = 1; i <= n; ++i) { 1156 b += -1.; 1157 z *= b / (a + b); 1158 } 1159 return w + log(z) + (gamln!T(a) + (gamln!T(b) - gsumln!T(a, b))); 1160 } else { // L50: reduction of A when B > 1000 1161 n = cast(int)(a - 1.); 1162 w = 1.; 1163 for (int i = 1; i <= n; ++i) { 1164 a += -1.; 1165 w *= a / (a / b + 1.); 1166 } 1167 return log(w) - n * log(b) + (gamln!T(a) + algdiv!T(a, b)); 1168 } 1169 1170 } else { 1171 /* ----------------------------------------------------------------------- */ 1172 // L60: A >= 8 1173 /* ----------------------------------------------------------------------- */ 1174 1175 T w = bcorr!T(a, b), h = a / b, 1176 u = -(a - 0.5) * log(h / (h + 1.)), 1177 v = b * alnrel!T(h); 1178 if (u > v) 1179 return log(b) * -0.5 + e + w - v - u; 1180 else 1181 return log(b) * -0.5 + e + w - u - v; 1182 } 1183 1184 } /* betaln */ 1185 1186 1187 1188 // called only once from bup(), as r = brcmp1(mu, a, b, x, y, FALSE) / a; 1189 static auto brcmp1(T)(int mu, T a, T b, T x, T y, int give_log) 1190 { 1191 /* ----------------------------------------------------------------------- 1192 * Evaluation of exp(mu) * x^a * y^b / beta(a,b) 1193 * ----------------------------------------------------------------------- */ 1194 1195 static T const__ = .398942280401433; /* == 1/sqrt(2*pi); */ 1196 /* R has M_1_SQRT_2PI */ 1197 1198 /* Local variables */ 1199 T c, t, u, v, z, a0, b0, apb; 1200 1201 a0 = min!T(a,b); 1202 if (a0 < 8.) { 1203 T lnx, lny; 1204 if (x <= .375) { 1205 lnx = log(x); 1206 lny = alnrel!T(-x); 1207 } else if (y > .375) { 1208 // L11: 1209 lnx = log(x); 1210 lny = log(y); 1211 } else { 1212 lnx = alnrel!T(-y); 1213 lny = log(y); 1214 } 1215 1216 // L20: 1217 z = a * lnx + b * lny; 1218 if (a0 >= 1.) { 1219 z -= betaln!T(a, b); 1220 return esum!T(mu, z, give_log); 1221 } 1222 // else : 1223 /* ----------------------------------------------------------------------- */ 1224 /* PROCEDURE FOR A < 1 OR B < 1 */ 1225 /* ----------------------------------------------------------------------- */ 1226 // L30: 1227 b0 = max!T(a,b); 1228 if (b0 >= 8.) { 1229 /* L80: ALGORITHM FOR b0 >= 8 */ 1230 u = gamln1!T(a0) + algdiv!T(a0, b0); 1231 //R_ifDEBUG_printf(" brcmp1(mu,a,b,*): a0 < 1, b0 >= 8; z=%.15g\n", z); 1232 return give_log ? log(a0) + esum!T(mu, z - u, 1) : a0 * esum!T(mu, z - u, 0); 1233 1234 } else if (b0 <= 1.) { 1235 // a0 < 1, b0 <= 1 1236 T ans = esum!T(mu, z, give_log); 1237 if (ans == (give_log ? T.infinity : 0.)) 1238 return ans; 1239 1240 apb = a + b; 1241 if (apb > 1.) { 1242 // L40: 1243 u = a + b - 1.; 1244 z = (gam1!T(u) + 1.) / apb; 1245 } else { 1246 z = gam1!T(apb) + 1.; 1247 } 1248 // L50: 1249 c = give_log ? log1p(gam1!T(a)) + log1p(gam1!T(b)) - log(z) : (gam1!T(a) + 1.) * (gam1!T(b) + 1.) / z; 1250 //R_ifDEBUG_printf(" brcmp1(mu,a,b,*): a0 < 1, b0 <= 1; c=%.15g\n", c); 1251 return give_log ? ans + log(a0) + c - log1p(a0 / b0) : ans * (a0 * c) / (a0 / b0 + 1.); 1252 } 1253 // else: algorithm for a0 < 1 < b0 < 8 1254 // L60: 1255 u = gamln1!T(a0); 1256 int n = cast(int)(b0 - 1.); 1257 if (n >= 1) { 1258 c = 1.; 1259 for (int i = 1; i <= n; ++i) { 1260 b0 += -1.; 1261 c *= b0 / (a0 + b0); 1262 /* L61: */ 1263 } 1264 u += log(c); // TODO?: log(c) = log( prod(...) ) = sum( log(...) ) 1265 } 1266 // L70: 1267 z -= u; 1268 b0 += -1.; 1269 apb = a0 + b0; 1270 if (apb > 1.) { 1271 // L71: 1272 t = (gam1!T(apb - 1.) + 1.) / apb; 1273 } else { 1274 t = gam1!T(apb) + 1.; 1275 } 1276 //R_ifDEBUG_printf(" brcmp1(mu,a,b,*): a0 < 1 < b0 < 8; t=%.15g\n", t); 1277 // L72: 1278 return give_log 1279 ? log(a0)+ esum!T(mu, z, 1) + log1p(gam1!T(b0)) - log(t) // TODO? log(t) = log1p(..) 1280 : a0 * esum!T(mu, z, 0) * (gam1!T(b0) + 1.) / t; 1281 1282 } else { 1283 1284 /* ----------------------------------------------------------------------- */ 1285 /* PROCEDURE FOR A >= 8 AND B >= 8 */ 1286 /* ----------------------------------------------------------------------- */ 1287 // L100: 1288 T h, x0, y0, lambda; 1289 if (a > b) { 1290 // L101: 1291 h = b / a; 1292 x0 = 1. / (h + 1.);// => lx0 := log(x0) = 0 - log1p(h) 1293 y0 = h / (h + 1.); 1294 lambda = (a + b) * y - b; 1295 } else { 1296 h = a / b; 1297 x0 = h / (h + 1.); // => lx0 := log(x0) = - log1p(1/h) 1298 y0 = 1. / (h + 1.); 1299 lambda = a - (a + b) * x; 1300 } 1301 T lx0 = -log1p(b/a); // in both cases 1302 1303 //R_ifDEBUG_printf(" brcmp1(mu,a,b,*): a,b >= 8; x0=%.15g, lx0=log(x0)=%.15g\n", x0, lx0); 1304 // L110: 1305 T e = -lambda / a; 1306 if (fabs(e) > 0.6) { 1307 // L111: 1308 u = e - log(x / x0); 1309 } else { 1310 u = rlog1!T(e); 1311 } 1312 1313 // L120: 1314 e = lambda / b; 1315 if (fabs(e) > 0.6) { 1316 // L121: 1317 v = e - log(y / y0); 1318 } else { 1319 v = rlog1!T(e); 1320 } 1321 1322 // L130: 1323 z = esum!T(mu, -(a*u + b*v), give_log); 1324 return give_log ? log(const__) + (log(b) + lx0)/2. + z - bcorr!T(a, b) 1325 : const__ * sqrt(b * x0)*z*exp(-bcorr!T(a, b)); 1326 } 1327 1328 } /* brcmp1 */ 1329 1330 1331 static auto brcomp(T)(T a, T b, T x, T y, int log_p) 1332 { 1333 /* ----------------------------------------------------------------------- 1334 * Evaluation of x^a * y^b / Beta(a,b) 1335 * ----------------------------------------------------------------------- */ 1336 1337 static T const__ = .398942280401433; /* == 1/sqrt(2*pi); */ 1338 /* R has M_1_SQRT_2PI , and M_LN_SQRT_2PI = ln(sqrt(2*pi)) = 0.918938.. */ 1339 int i, n; 1340 T c, e, u, v, z, a0, b0, apb; 1341 1342 mixin R_D__0!log_p; 1343 if (x == 0. || y == 0.) { 1344 return R_D__0; 1345 } 1346 a0 = min!T(a, b); 1347 if (a0 < 8.) { 1348 T lnx, lny; 1349 if (x <= .375) { 1350 lnx = log(x); 1351 lny = alnrel!T(-x); 1352 } 1353 else { 1354 if (y > .375) { 1355 lnx = log(x); 1356 lny = log(y); 1357 } else { 1358 lnx = alnrel!T(-y); 1359 lny = log(y); 1360 } 1361 } 1362 1363 z = a * lnx + b * lny; 1364 if (a0 >= 1.) { 1365 z -= betaln!T(a, b); 1366 return R_D_exp!T(z, log_p); 1367 } 1368 1369 /* ----------------------------------------------------------------------- */ 1370 /* PROCEDURE FOR a < 1 OR b < 1 */ 1371 /* ----------------------------------------------------------------------- */ 1372 1373 b0 = max(a, b); 1374 if (b0 >= 8.) { /* L80: */ 1375 u = gamln1!T(a0) + algdiv!T(a0, b0); 1376 1377 return (log_p ? log(a0) + (z - u) : a0 * exp(z - u)); 1378 } 1379 /* else : */ 1380 1381 if (b0 <= 1.) { /* algorithm for max(a,b) = b0 <= 1 */ 1382 1383 T e_z = R_D_exp!T(z, log_p); 1384 1385 if (!log_p && e_z == 0.) /* exp() underflow */ 1386 return 0.; 1387 1388 apb = a + b; 1389 if (apb > 1.) { 1390 u = a + b - 1.; 1391 z = (gam1!T(u) + 1.) / apb; 1392 } else { 1393 z = gam1!T(apb) + 1.; 1394 } 1395 1396 c = (gam1!T(a) + 1.) * (gam1!T(b) + 1.) / z; 1397 /* FIXME? log(a0*c)= log(a0)+ log(c) and that is improvable */ 1398 return (log_p ? e_z + log(a0 * c) - log1p(a0/b0) 1399 : e_z * (a0 * c) / (a0 / b0 + 1.)); 1400 } 1401 1402 /* else : ALGORITHM FOR 1 < b0 < 8 */ 1403 1404 u = gamln1!T(a0); 1405 n = cast(int)(b0 - 1.); 1406 if (n >= 1) { 1407 c = 1.; 1408 for (i = 1; i <= n; ++i) { 1409 b0 += -1.; 1410 c *= b0 / (a0 + b0); 1411 } 1412 u = log(c) + u; 1413 } 1414 z -= u; 1415 b0 += -1.; 1416 apb = a0 + b0; 1417 T t; 1418 if (apb > 1.) { 1419 u = a0 + b0 - 1.; 1420 t = (gam1!T(u) + 1.) / apb; 1421 } else { 1422 t = gam1!T(apb) + 1.; 1423 } 1424 1425 return (log_p ? log(a0) + z + log1p(gam1!T(b0)) - log(t) : a0 * exp(z) * (gam1!T(b0) + 1.) / t); 1426 1427 } else { 1428 /* ----------------------------------------------------------------------- */ 1429 /* PROCEDURE FOR A >= 8 AND B >= 8 */ 1430 /* ----------------------------------------------------------------------- */ 1431 T h, x0, y0, lambda; 1432 if (a <= b) { 1433 h = a / b; 1434 x0 = h / (h + 1.); 1435 y0 = 1. / (h + 1.); 1436 lambda = a - (a + b) * x; 1437 } else { 1438 h = b / a; 1439 x0 = 1. / (h + 1.); 1440 y0 = h / (h + 1.); 1441 lambda = (a + b) * y - b; 1442 } 1443 1444 e = -lambda / a; 1445 if (fabs(e) > .6) 1446 u = e - log(x / x0); 1447 else 1448 u = rlog1!T(e); 1449 1450 e = lambda / b; 1451 if (fabs(e) <= .6) 1452 v = rlog1!T(e); 1453 else 1454 v = e - log(y / y0); 1455 1456 z = log_p ? -(a * u + b * v) : exp(-(a * u + b * v)); 1457 1458 return(log_p ? -M_LN_SQRT_2PI + .5*log(b * x0) + z - bcorr!T(a,b) 1459 : const__ * sqrt(b * x0) * z * exp(-bcorr!T(a, b))); 1460 } 1461 } /* brcomp */ 1462 1463 1464 1465 static auto bfrac(T)(T a, T b, T x, T y, T lambda, T eps, int log_p) 1466 { 1467 /* ----------------------------------------------------------------------- 1468 Continued fraction expansion for I_x(a,b) when a, b > 1. 1469 It is assumed that lambda = (a + b)*y - b. 1470 -----------------------------------------------------------------------*/ 1471 1472 T c, e, n, p, r, s, t, w, c0, c1, r0, an, bn, yp1, anp1, bnp1, beta, alpha, brc; 1473 1474 if(!isFinite(lambda)) 1475 return T.nan;// TODO: can return 0 or 1 (?) 1476 //R_ifDEBUG_printf(" bfrac(a=%g, b=%g, x=%g, y=%g, lambda=%g, eps=%g, log_p=%d):", a, b, x, y, lambda, eps, log_p); 1477 brc = brcomp!T(a, b, x, y, log_p); 1478 if(isNaN(brc)) { // e.g. from L <- 1e308; pnbinom(L, L, mu = 5) 1479 //R_ifDEBUG_printf(" --> brcomp(a,b,x,y) = NaN\n"); 1480 return T.nan; // TODO: could we know better? 1481 } 1482 if (!log_p && brc == 0.) { 1483 //R_ifDEBUG_printf(" --> brcomp(a,b,x,y) underflowed to 0.\n"); 1484 return 0.; 1485 } 1486 //#ifdef DEBUG_bratio 1487 // else 1488 // REprintf("\n"); 1489 //#endif 1490 1491 c = lambda + 1.; 1492 c0 = b / a; 1493 c1 = 1. / a + 1.; 1494 yp1 = y + 1.; 1495 1496 n = 0.; 1497 p = 1.; 1498 s = a + 1.; 1499 an = 0.; 1500 bn = 1.; 1501 anp1 = 1.; 1502 bnp1 = c / c1; 1503 r = c1 / c; 1504 1505 /* CONTINUED FRACTION CALCULATION */ 1506 1507 do { 1508 n += 1.; 1509 t = n/a; 1510 w = n*(b - n)*x; 1511 e = a / s; 1512 alpha = p * (p + c0)*e*e*(w*x); 1513 e = (t + 1.) / (c1 + t + t); 1514 beta = n + w / s + e*(c + n*yp1); 1515 p = t + 1.; 1516 s += 2.; 1517 1518 /* update an, bn, anp1, and bnp1 */ 1519 1520 t = alpha*an + beta*anp1; an = anp1; anp1 = t; 1521 t = alpha*bn + beta*bnp1; bn = bnp1; bnp1 = t; 1522 1523 r0 = r; 1524 r = anp1 / bnp1; 1525 //#ifdef _not_normally_DEBUG_bfrac 1526 // R_ifDEBUG_printf(" n=%5.0f, a_{n,n+1}= (%12g,%12g), b_{n,n+1} = (%12g,%12g) => r0,r = (%14g,%14g)\n", 1527 // n, an,anp1, bn,bnp1, r0, r); 1528 //#endif 1529 if (fabs(r - r0) <= eps * r) 1530 break; 1531 1532 /* rescale an, bn, anp1, and bnp1 */ 1533 1534 an /= bnp1; 1535 bn /= bnp1; 1536 anp1 = r; 1537 bnp1 = 1.; 1538 } while (n < 10000);// arbitrary; had '1' --> infinite loop for lambda = Inf 1539 //R_ifDEBUG_printf(" in bfrac(): n=%.0f terms cont.frac.; brc=%g, r=%g\n", n, brc, r); 1540 //if(n >= 10000 && fabs(r - r0) > eps * r) 1541 // MATHLIB_WARNING5(" bfrac(a=%g, b=%g, x=%g, y=%g, lambda=%g) did *not* converge (in 10000 steps)\n", a, b, x, y, lambda); 1542 return (log_p ? brc + log(r) : brc*r); 1543 } /* bfrac */ 1544 1545 1546 1547 static auto bup(T)(T a, T b, T x, T y, int n, T eps, int give_log) 1548 { 1549 /* ----------------------------------------------------------------------- */ 1550 /* EVALUATION OF I_x(A,B) - I_x(A+N,B) WHERE N IS A POSITIVE INT. */ 1551 /* EPS IS THE TOLERANCE USED. */ 1552 /* ----------------------------------------------------------------------- */ 1553 1554 T ret_val; 1555 int i, k, mu; 1556 T d, l; 1557 1558 // Obtain the scaling factor exp(-mu) and exp(mu)*(x^a * y^b / beta(a,b))/a 1559 1560 T apb = a + b, 1561 ap1 = a + 1.; 1562 if (n > 1 && a >= 1. && apb >= ap1 * 1.1) { 1563 mu = cast(int)fabs(exparg!T(1)); 1564 k = cast(int)exparg!T(0); 1565 if (mu > k) 1566 mu = k; 1567 d = exp(-cast(T) mu); 1568 } else { 1569 mu = 0; 1570 d = 1.; 1571 } 1572 1573 /* L10: */ 1574 ret_val = give_log ? brcmp1!T(mu, a, b, x, y, 1) - log(a) : brcmp1!T(mu, a, b, x, y, 0) / a; 1575 if(n == 1 || (give_log && ret_val == -T.infinity) || (!give_log && ret_val == 0.)) 1576 return ret_val; 1577 1578 int nm1 = n - 1; 1579 T w = d; 1580 1581 /* LET K BE THE INDEX OF THE MAXIMUM TERM */ 1582 1583 k = 0; 1584 if (b > 1.) { 1585 if (y > 1e-4) { 1586 T r = (b - 1.) * x / y - a; 1587 if (r >= 1.) 1588 k = (r < nm1) ? cast(int)r : nm1; 1589 } else 1590 k = nm1; 1591 1592 // ADD THE INCREASING TERMS OF THE SERIES - if k > 0 1593 /* L30: */ 1594 for (i = 0; i < k; ++i) { 1595 l = cast(T) i; 1596 d *= (apb + l) / (ap1 + l) * x; 1597 w += d; 1598 } 1599 } 1600 1601 // L40: ADD THE REMAINING TERMS OF THE SERIES 1602 1603 for (i = k; i < nm1; ++i) { 1604 l = cast(T) i; 1605 d *= (apb + l) / (ap1 + l) * x; 1606 w += d; 1607 if (d <= eps * w) /* relativ convergence (eps) */ 1608 break; 1609 } 1610 1611 // L50: TERMINATE THE PROCEDURE 1612 if(give_log) { 1613 ret_val += log(w); 1614 } else 1615 ret_val *= w; 1616 1617 return ret_val; 1618 } /* bup */ 1619 1620 1621 static auto bpser(T)(T a, T b, T x, T eps, int log_p) 1622 { 1623 /* ----------------------------------------------------------------------- 1624 * Power SERies expansion for evaluating I_x(a,b) when 1625 * b <= 1 or b*x <= 0.7. eps is the tolerance used. 1626 * NB: if log_p is TRUE, also use it if (b < 40 & lambda > 650) 1627 * ----------------------------------------------------------------------- */ 1628 1629 int i, m; 1630 T ans, c, t, u, z, a0, b0, apb; 1631 mixin R_D__0!log_p; 1632 1633 if (x == 0.) { 1634 return R_D__0; 1635 } 1636 /* ----------------------------------------------------------------------- */ 1637 /* compute the factor x^a/(a*Beta(a,b)) */ 1638 /* ----------------------------------------------------------------------- */ 1639 a0 = min!T(a,b); 1640 if (a0 >= 1.) { /* ------ 1 <= a0 <= b0 ------ */ 1641 z = a * log(x) - betaln!T(a, b); 1642 ans = log_p ? z - log(a) : exp(z) / a; 1643 } 1644 else { 1645 b0 = max!T(a,b); 1646 1647 if (b0 < 8.) { 1648 1649 if (b0 <= 1.) { /* ------ a0 < 1 and b0 <= 1 ------ */ 1650 1651 if(log_p) { 1652 ans = a * log(x); 1653 } else { 1654 ans = pow(x, a); 1655 if (ans == 0.) /* once underflow, always underflow .. */ 1656 return ans; 1657 } 1658 apb = a + b; 1659 if (apb > 1.) { 1660 u = a + b - 1.; 1661 z = (gam1!T(u) + 1.) / apb; 1662 } else { 1663 z = gam1!T(apb) + 1.; 1664 } 1665 c = (gam1!T(a) + 1.) * (gam1!T(b) + 1.) / z; 1666 1667 if(log_p) /* FIXME ? -- improve quite a bit for c ~= 1 */ 1668 ans += log(c * (b / apb)); 1669 else 1670 ans *= c * (b / apb); 1671 1672 } else { /* ------ a0 < 1 < b0 < 8 ------ */ 1673 1674 u = gamln1!T(a0); 1675 m = cast(int)(b0 - 1.); 1676 if (m >= 1) { 1677 c = 1.; 1678 for (i = 1; i <= m; ++i) { 1679 b0 += -1.; 1680 c *= b0 / (a0 + b0); 1681 } 1682 u += log(c); 1683 } 1684 1685 z = a * log(x) - u; 1686 b0 += -1.; // => b0 in (0, 7) 1687 apb = a0 + b0; 1688 if (apb > 1.) { 1689 u = a0 + b0 - 1.; 1690 t = (gam1!T(u) + 1.) / apb; 1691 } else { 1692 t = gam1!T(apb) + 1.; 1693 } 1694 1695 if(log_p) /* FIXME? potential for improving log(t) */ 1696 ans = z + log(a0 / a) + log1p(gam1!T(b0)) - log(t); 1697 else 1698 ans = exp(z) * (a0 / a) * (gam1!T(b0) + 1.) / t; 1699 } 1700 1701 } else { /* ------ a0 < 1 < 8 <= b0 ------ */ 1702 1703 u = gamln1!T(a0) + algdiv!T(a0, b0); 1704 z = a * log(x) - u; 1705 1706 if(log_p) 1707 ans = z + log(a0 / a); 1708 else 1709 ans = a0 / a * exp(z); 1710 } 1711 } 1712 //R_ifDEBUG_printf(" bpser(a=%g, b=%g, x=%g, log=%d): prelim.ans = %.14g;\n", a,b,x, log_p, ans); 1713 if (ans == R_D__0 || (!log_p && a <= eps * 0.1)) { 1714 return ans; 1715 } 1716 1717 /* ----------------------------------------------------------------------- */ 1718 /* COMPUTE THE SERIES */ 1719 /* ----------------------------------------------------------------------- */ 1720 T tol = eps / a, 1721 n = 0., 1722 sum = 0., w; 1723 c = 1.; 1724 do { // sum is alternating as long as n < b (<==> 1 - b/n < 0) 1725 n += 1.; 1726 c *= (0.5 - b / n + 0.5) * x; 1727 w = c / (a + n); 1728 sum += w; 1729 } while (n < 1e7 && fabs(w) > tol); 1730 if(fabs(w) > tol) { // the series did not converge (in time) 1731 // warn only when the result seems to matter: 1732 //if(( log_p && !(a*sum > -1. && fabs(log1p(a * sum)) < eps*fabs(ans))) || (!log_p && fabs(a*sum + 1.) != 1.)) 1733 //MATHLIB_WARNING5(" bpser(a=%g, b=%g, x=%g,...) did not converge (n=1e7, |w|/tol=%g > 1; A=%g)", a,b,x, fabs(w)/tol, ans); 1734 } 1735 //R_ifDEBUG_printf(" -> n=%.0f iterations, |w|=%g %s %g=tol:=eps/a ==> a*sum=%g\n", n, fabs(w), (fabs(w) > tol) ? ">!!>" : "<=", tol, a*sum); 1736 if(log_p) { 1737 if (a*sum > -1.) 1738 ans += log1p(a*sum); 1739 else { 1740 //if(ans > ML_NEGINF) 1741 //MATHLIB_WARNING3("pbeta(*, log.p=TRUE) -> bpser(a=%g, b=%g, x=%g,...) underflow to -Inf", a, b, x); 1742 ans = -T.infinity; 1743 } 1744 } else if (a*sum > -1.) 1745 ans *= (a*sum + 1.); 1746 else // underflow to 1747 ans = 0.; 1748 return ans; 1749 } /* bpser */ 1750 1751 1752 1753 static auto apser(T)(T a, T b, T x, T eps) 1754 { 1755 /* ----------------------------------------------------------------------- 1756 * apser() yields the incomplete beta ratio I_{1-x}(b,a) for 1757 * a <= min(eps,eps*b), b*x <= 1, and x <= 0.5, i.e., a is very small. 1758 * Use only if above inequalities are satisfied. 1759 * ----------------------------------------------------------------------- */ 1760 1761 static const(T) g = .577215664901533; 1762 1763 T tol, c, j, s, t, aj; 1764 T bx = b * x; 1765 1766 t = x - bx; 1767 if (b * eps <= 0.02) 1768 c = log(x) + psi!T(b) + g + t; 1769 else // b > 2e13 : psi(b) ~= log(b) 1770 c = log(bx) + g + t; 1771 1772 tol = eps * 5. * fabs(c); 1773 j = 1.; 1774 s = 0.; 1775 do { 1776 j += 1.; 1777 t *= x - bx / j; 1778 aj = t / j; 1779 s += aj; 1780 } while (fabs(aj) > tol); 1781 1782 return -a*(c + s); 1783 } /* apser */ 1784 1785 1786 auto fpser(T)(T a, T b, T x, T eps, int log_p) 1787 { 1788 /* ----------------------------------------------------------------------- * 1789 1790 * EVALUATION OF I (A,B) 1791 * X 1792 1793 * FOR B < MIN(EPS, EPS*A) AND X <= 0.5 1794 1795 * ----------------------------------------------------------------------- */ 1796 1797 T ans, c, s, t, an, tol; 1798 1799 /* SET ans := x^a : */ 1800 if (log_p) { 1801 ans = a * log(x); 1802 } else if (a > eps * 0.001) { 1803 t = a * log(x); 1804 if (t < exparg!T(1)) { /* exp(t) would underflow */ 1805 return 0.; 1806 } 1807 ans = exp(t); 1808 } else 1809 ans = 1.; 1810 1811 /* NOTE THAT 1/B(A,B) = B */ 1812 1813 if (log_p) 1814 ans += log(b) - log(a); 1815 else 1816 ans *= b / a; 1817 1818 tol = eps / a; 1819 an = a + 1.; 1820 t = x; 1821 s = t / an; 1822 do { 1823 an += 1.; 1824 t = x * t; 1825 c = t / an; 1826 s += c; 1827 } while (fabs(c) > tol); 1828 1829 if (log_p) 1830 ans += log1p(a * s); 1831 else 1832 ans *= a * s + 1.; 1833 return ans; 1834 } /* fpser */ 1835 1836 immutable(string) SET_0_noswap(){ 1837 return `a0 = a; x0 = x; b0 = b; y0 = y;`; 1838 } 1839 1840 immutable(string) SET_0_swap(){ 1841 return `a0 = b; x0 = y, b0 = a; y0 = x;`; 1842 } 1843 1844 void bratio(T)(T a, T b, T x, T y, T* w, T* w1, int *ierr, int log_p) 1845 { 1846 /* ----------------------------------------------------------------------- 1847 1848 * Evaluation of the Incomplete Beta function I_x(a,b) 1849 1850 * -------------------- 1851 1852 * It is assumed that a and b are nonnegative, and that x <= 1 1853 * and y = 1 - x. Bratio assigns w and w1 the values 1854 1855 * w = I_x(a,b) 1856 * w1 = 1 - I_x(a,b) 1857 1858 * ierr is a variable that reports the status of the results. 1859 * If no input errors are detected then ierr is set to 0 and 1860 * w and w1 are computed. otherwise, if an error is detected, 1861 * then w and w1 are assigned the value 0 and ierr is set to 1862 * one of the following values ... 1863 1864 * ierr = 1 if a or b is negative 1865 * ierr = 2 if a = b = 0 1866 * ierr = 3 if x < 0 or x > 1 1867 * ierr = 4 if y < 0 or y > 1 1868 * ierr = 5 if x + y != 1 1869 * ierr = 6 if x = a = 0 1870 * ierr = 7 if y = b = 0 1871 * ierr = 8 (not used currently) 1872 * ierr = 9 NaN in a, b, x, or y 1873 * ierr = 10 (not used currently) 1874 * ierr = 11 bgrat() error code 1 [+ warning in bgrat()] 1875 * ierr = 12 bgrat() error code 2 (no warning here) 1876 * ierr = 13 bgrat() error code 3 (no warning here) 1877 * ierr = 14 bgrat() error code 4 [+ WARNING in bgrat()] 1878 1879 1880 * -------------------- 1881 * Written by Alfred H. Morris, Jr. 1882 * Naval Surface Warfare Center 1883 * Dahlgren, Virginia 1884 * Revised ... Nov 1991 1885 * ----------------------------------------------------------------------- */ 1886 1887 // Rboolean 1888 int do_swap, a_lt_b; 1889 int n, ierr1 = 0; 1890 T z, a0, b0, x0, y0, lambda; 1891 1892 /* eps is a machine dependent constant: the smallest 1893 * floating point number for which 1. + eps > 1. 1894 * NOTE: for almost all purposes it is replaced by 1e-15 (~= 4.5 times larger) below */ 1895 T eps = 2.*d1mach!T(3); /* == DBL_EPSILON (in R, Rmath) */ 1896 1897 /* ----------------------------------------------------------------------- */ 1898 mixin R_D__0!log_p; 1899 mixin R_D__1!log_p; 1900 *w = R_D__0; 1901 *w1 = R_D__0; 1902 1903 //#ifdef IEEE_754 1904 // safeguard, preventing infinite loops further down 1905 if (isNaN(x) || isNaN(y) || isNaN(a) || isNaN(b)){ 1906 *ierr = 9; 1907 return; 1908 } 1909 //#endif 1910 if(a < 0. || b < 0.){ 1911 *ierr = 1; 1912 return; 1913 } 1914 if(a == 0. && b == 0.){ 1915 *ierr = 2; return; 1916 } 1917 if(x < 0. || x > 1.){ 1918 *ierr = 3; return; 1919 } 1920 if(y < 0. || y > 1.){ 1921 *ierr = 4; return; 1922 } 1923 1924 /* check that 'y == 1 - x' : */ 1925 z = x + y - 0.5 - 0.5; 1926 1927 if (fabs(z) > eps * 3.){ 1928 *ierr = 5; return; 1929 } 1930 1931 //R_ifDEBUG_printf("bratio(a=%g, b=%g, x=%9g, y=%9g, .., log_p=%d): ", 1932 // a,b,x,y, log_p); 1933 *ierr = 0; 1934 if(x == 0.) 1935 goto L200; 1936 if(y == 0.) 1937 goto L210; 1938 1939 if (a == 0.) goto L211; 1940 if (b == 0.) goto L201; 1941 1942 eps = max!T(eps, 1e-15); 1943 //Rboolean 1944 a_lt_b = (a < b); 1945 if (/* max(a,b) */ (a_lt_b ? b : a) < eps * .001){ /* procedure for a and b < 0.001 * eps */ 1946 // L230: -- result *independent* of x (!) 1947 // *w = a/(a+b) and w1 = b/(a+b) : 1948 if(log_p) { 1949 if(a_lt_b) { 1950 *w = log1p(-a/(a+b)); // notably if a << b 1951 *w1 = log(a/(a+b)); 1952 } else { // b <= a 1953 *w = log(b/(a+b)); 1954 *w1 = log1p(-b/(a+b)); 1955 } 1956 } else { 1957 *w = b/(a + b); 1958 *w1 = a/(a + b); 1959 } 1960 1961 //R_ifDEBUG_printf("a & b very small -> simple ratios (%g,%g)\n", *w,*w1); 1962 return; 1963 } 1964 1965 //#define SET_0_noswap \ 1966 // a0 = a; x0 = x; \ 1967 // b0 = b; y0 = y; 1968 1969 //#define SET_0_swap \ 1970 // a0 = b; x0 = y; \ 1971 // b0 = a; y0 = x; 1972 1973 if (min!T(a,b) <= 1.) { /*------------------------ a <= 1 or b <= 1 ---- */ 1974 1975 do_swap = (x > 0.5); 1976 if (do_swap) { 1977 mixin(SET_0_swap()); 1978 } else { 1979 mixin(SET_0_noswap()); 1980 } 1981 /* now have x0 <= 1/2 <= y0 (still x0+y0 == 1) */ 1982 1983 //R_ifDEBUG_printf(" min(a,b) <= 1, do_swap=%d;", do_swap); 1984 1985 if (b0 < min!T(eps, eps * a0)) { /* L80: */ 1986 *w = fpser!T(a0, b0, x0, eps, log_p); 1987 *w1 = log_p ? R_Log1_Exp!T(*w) : 0.5 - *w + 0.5; 1988 //R_ifDEBUG_printf(" b0 small -> w := fpser(*) = %.15g\n", *w); 1989 goto L_end; 1990 } 1991 1992 if (a0 < min!T(eps, eps * b0) && b0 * x0 <= 1.) { /* L90: */ 1993 *w1 = apser!T(a0, b0, x0, eps); 1994 //R_ifDEBUG_printf(" a0 small -> w1 := apser(*) = %.15g\n", *w1); 1995 goto L_end_from_w1; 1996 } 1997 1998 //Rboolean 1999 int did_bup = 0; 2000 if (max!T(a0,b0) > 1.) { /* L20: min(a,b) <= 1 < max(a,b) */ 2001 //R_ifDEBUG_printf("\n L20: min(a,b) <= 1 < max(a,b); "); 2002 if (b0 <= 1.) goto L_w_bpser; 2003 2004 if (x0 >= 0.29) /* was 0.3, PR#13786 */ goto L_w1_bpser; 2005 2006 if (x0 < 0.1 && pow(x0*b0, a0) <= 0.7) goto L_w_bpser; 2007 2008 if (b0 > 15.) { 2009 *w1 = 0.; 2010 goto L131; 2011 } 2012 } else { /* a, b <= 1 */ 2013 //R_ifDEBUG_printf("\n both a,b <= 1; "); 2014 if (a0 >= min!T(0.2, b0)) goto L_w_bpser; 2015 2016 if (pow(x0, a0) <= 0.9) goto L_w_bpser; 2017 2018 if (x0 >= 0.3) goto L_w1_bpser; 2019 } 2020 n = 20; /* goto L130; */ 2021 *w1 = bup!T(b0, a0, y0, x0, n, eps, 0); did_bup = 1; 2022 //R_ifDEBUG_printf(" ... n=20 and *w1 := bup(*) = %.15g; ", *w1); 2023 b0 += n; 2024 L131: 2025 //R_ifDEBUG_printf(" L131: bgrat(*, w1=%.15g) ", *w1); 2026 bgrat!T(b0, a0, y0, x0, w1, 15*eps, &ierr1, 0); 2027 //#ifdef DEBUG_bratio 2028 // REprintf(" ==> new w1=%.15g", *w1); 2029 // if(ierr1) REprintf(" ERROR(code=%d)\n", ierr1) ; else REprintf("\n"); 2030 //#endif 2031 if(*w1 == 0 || (0 < *w1 && *w1 < DBL_MIN)) { // w1=0 or very close: 2032 // "almost surely" from underflow, try more: [2013-03-04] 2033 // FIXME: it is even better to do this in bgrat *directly* at least for the case 2034 // !did_bup, i.e., where *w1 = (0 or -Inf) on entry 2035 //R_ifDEBUG_printf(" denormalized or underflow (?) -> retrying: "); 2036 if(did_bup) { // re-do that part on log scale: 2037 *w1 = bup!T(b0-n, a0, y0, x0, n, eps, 1); 2038 } 2039 else *w1 = -T.infinity; // = 0 on log-scale 2040 bgrat!T(b0, a0, y0, x0, w1, 15*eps, &ierr1, 1); 2041 if(ierr1) 2042 *ierr = 10 + ierr1; 2043 //#ifdef DEBUG_bratio 2044 // REprintf(" ==> new log(w1)=%.15g", *w1); 2045 // if(ierr1) REprintf(" Error(code=%d)\n", ierr1) ; else REprintf("\n"); 2046 //#endif 2047 goto L_end_from_w1_log; 2048 } 2049 // else 2050 if(ierr1) 2051 *ierr = 10 + ierr1; 2052 //if(*w1 < 0) 2053 // MATHLIB_WARNING4("bratio(a=%g, b=%g, x=%g): bgrat() -> w1 = %g", a, b, x, *w1); 2054 goto L_end_from_w1; 2055 } else { /* L30: -------------------- both a, b > 1 {a0 > 1 & b0 > 1} ---*/ 2056 2057 /* lambda := a y - b x = (a + b)y = a - (a+b)x {using x + y == 1}, 2058 * ------ using the numerically best version : */ 2059 lambda = isFinite(a + b) ? ((a > b) ? (a + b) * y - b : a - (a + b) * x) : a*y - b*x; 2060 do_swap = (lambda < 0.); 2061 if (do_swap) { 2062 lambda = -lambda; 2063 mixin (SET_0_swap()); 2064 } else { 2065 mixin (SET_0_noswap()); 2066 } 2067 2068 //R_ifDEBUG_printf(" L30: both a, b > 1; |lambda| = %#g, do_swap = %d\n", lambda, do_swap); 2069 2070 if (b0 < 40.) { 2071 //R_ifDEBUG_printf(" b0 < 40;"); 2072 if (b0 * x0 <= 0.7 || (log_p && lambda > 650.)) // << added 2010-03; svn r51327 2073 goto L_w_bpser; 2074 else 2075 goto L140; 2076 } 2077 else if (a0 > b0) { /* ---- a0 > b0 >= 40 ---- */ 2078 //R_ifDEBUG_printf(" a0 > b0 >= 40;"); 2079 if (b0 <= 100. || lambda > b0 * 0.03) 2080 goto L_bfrac; 2081 2082 } else if (a0 <= 100.) { 2083 //R_ifDEBUG_printf(" a0 <= 100; a0 <= b0 >= 40;"); 2084 goto L_bfrac; 2085 } 2086 else if (lambda > a0 * 0.03) { 2087 //R_ifDEBUG_printf(" b0 >= a0 > 100; lambda > a0 * 0.03 "); 2088 goto L_bfrac; 2089 } 2090 2091 /* else if none of the above L180: */ 2092 *w = basym!T(a0, b0, lambda, eps * 100., log_p); 2093 *w1 = log_p ? R_Log1_Exp!T(*w) : 0.5 - *w + 0.5; 2094 //R_ifDEBUG_printf(" b0 >= a0 > 100; lambda <= a0 * 0.03: *w:= basym(*) =%.15g\n", *w); 2095 goto L_end; 2096 2097 } /* else: a, b > 1 */ 2098 2099 /* EVALUATION OF THE APPROPRIATE ALGORITHM */ 2100 2101 L_w_bpser: // was L100 2102 *w = bpser!T(a0, b0, x0, eps, log_p); 2103 *w1 = log_p ? R_Log1_Exp!T(*w) : 0.5 - *w + 0.5; 2104 //R_ifDEBUG_printf(" L_w_bpser: *w := bpser(*) = %.15g\n", *w); 2105 goto L_end; 2106 2107 L_w1_bpser: // was L110 2108 *w1 = bpser!T(b0, a0, y0, eps, log_p); 2109 *w = log_p ? R_Log1_Exp!T(*w1) : 0.5 - *w1 + 0.5; 2110 //R_ifDEBUG_printf(" L_w1_bpser: *w1 := bpser(*) = %.15g\n", *w1); 2111 goto L_end; 2112 2113 L_bfrac: 2114 *w = bfrac!T(a0, b0, x0, y0, lambda, eps * 15., log_p); 2115 *w1 = log_p ? R_Log1_Exp!T(*w) : 0.5 - *w + 0.5; 2116 //R_ifDEBUG_printf(" L_bfrac: *w := bfrac(*) = %g\n", *w); 2117 goto L_end; 2118 2119 L140: 2120 /* b0 := fractional_part( b0 ) in (0, 1] */ 2121 n = cast(int) b0; 2122 b0 -= n; 2123 if (b0 == 0.) { 2124 --n; b0 = 1.; 2125 } 2126 2127 *w = bup!T(b0, a0, y0, x0, n, eps, 0); 2128 2129 if(*w < DBL_MIN && log_p) { /* do not believe it; try bpser() : */ 2130 //R_ifDEBUG_printf(" L140: bup(b0=%g,..)=%.15g < DBL_MIN - not used; ", b0, *w); 2131 /*revert: */ b0 += n; 2132 /* which is only valid if b0 <= 1 || b0*x0 <= 0.7 */ 2133 goto L_w_bpser; 2134 } 2135 //R_ifDEBUG_printf(" L140: *w := bup(b0=%g,..) = %.15g; ", b0, *w); 2136 if (x0 <= 0.7) { 2137 /* log_p : TODO: w = bup(.) + bpser(.) -- not so easy to use log-scale */ 2138 *w += bpser!T(a0, b0, x0, eps, /* log_p = */ 0); 2139 //R_ifDEBUG_printf(" x0 <= 0.7: *w := *w + bpser(*) = %.15g\n", *w); 2140 goto L_end_from_w; 2141 } 2142 /* L150: */ 2143 if (a0 <= 15.) { 2144 n = 20; 2145 *w += bup!T(a0, b0, x0, y0, n, eps, 0); 2146 //R_ifDEBUG_printf("\n a0 <= 15: *w := *w + bup(*) = %.15g;", *w); 2147 a0 += n; 2148 } 2149 //R_ifDEBUG_printf(" bgrat(*, w=%.15g) ", *w); 2150 bgrat!T(a0, b0, x0, y0, w, 15*eps, &ierr1, 0); 2151 if(ierr1) *ierr = 10 + ierr1; 2152 //#ifdef DEBUG_bratio 2153 // REprintf("==> new w=%.15g", *w); 2154 // if(ierr1) REprintf(" Error(code=%d)\n", ierr1) ; else REprintf("\n"); 2155 //#endif 2156 goto L_end_from_w; 2157 2158 2159 /* TERMINATION OF THE PROCEDURE */ 2160 2161 L200: 2162 if (a == 0.) { *ierr = 6; return; } 2163 // else: 2164 L201: *w = R_D__0; *w1 = R_D__1; return; 2165 2166 L210: 2167 if (b == 0.) { *ierr = 7; return; } 2168 // else: 2169 L211: *w = R_D__1; *w1 = R_D__0; return; 2170 2171 L_end_from_w: 2172 if(log_p) { 2173 *w1 = log1p(-*w); 2174 *w = log(*w); 2175 } else { 2176 *w1 = 0.5 - *w + 0.5; 2177 } 2178 goto L_end; 2179 2180 L_end_from_w1: 2181 if(log_p) { 2182 *w = log1p(-*w1); 2183 *w1 = log(*w1); 2184 } else { 2185 *w = 0.5 - *w1 + 0.5; 2186 } 2187 goto L_end; 2188 2189 L_end_from_w1_log: 2190 // *w1 = log(w1) already; w = 1 - w1 ==> log(w) = log(1 - w1) = log(1 - exp(*w1)) 2191 if(log_p) { 2192 *w = R_Log1_Exp!T(*w1); 2193 } else { 2194 *w = /* 1 - exp(*w1) */ -expm1(*w1); 2195 *w1 = exp(*w1); 2196 } 2197 goto L_end; 2198 2199 2200 L_end: 2201 if (do_swap) { /* swap */ 2202 T t = *w; *w = *w1; *w1 = t; 2203 } 2204 return; 2205 2206 } /* bratio */ 2207