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