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