The R Project SVN – R – Diff – /trunk/src/nmath/toms708.c

    to make log1p and expm1 declared */
 /**----------- DEBUGGING -------------
+ *
  *	make CFLAGS='-DDEBUG_bratio  ...'
+ *MM:
- *MM: (cd `R-devel 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)
+ (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)
 */
 #ifdef DEBUG_bratio
 /* for REprintf */
 # include <R_ext/Print.h>
 #endif
 #undef R_Log1_Exp
 #define R_Log1_Exp(x)   ((x) > -M_LN2 ? log(-rexpm1(x)) : log1p(-exp(x)))
 static double bfrac(double, double, double, double, double, double, int log_p);
-static void bgrat(double, double, double, double, double *, double, int *);
+static void bgrat(double, double, double, double, double *, double, int *, Rboolean log_w);
-static void grat1(double, double, double, double *, double *, double);
+static double grat_r(double a, double x, double r, double eps);
 static double apser(double, double, double, double);
 static double bpser(double, double, double, double, int log_p);
 static double basym(double, double, double, double, int log_p);
 static double fpser(double, double, double, double, int log_p);
-static double bup(double, double, double, double, int, double);
+static double bup(double, double, double, double, int, double, int give_log);
 static double exparg(int);
 static double psi(double);
 static double gam1(double);
 static double gamln1(double);
 static double betaln(double, double);
 	if (b0 < min(eps, eps * a0)) { /* L80: */
 	    *w = fpser(a0, b0, x0, eps, log_p);
 	    *w1 = log_p ? R_Log1_Exp(*w) : 0.5 - *w + 0.5;
 #ifdef DEBUG_bratio
-	    REprintf("  b0 small -> w := fpser(*) = %15g\n", *w);
+	    REprintf("  b0 small -> w := fpser(*) = %.15g\n", *w);
 #endif
-	    goto L_end_after_log;
+	    goto L_end;
+	}
 	if (a0 < min(eps, eps * b0) && b0 * x0 <= 1.0) { /* L90: */
 	    *w1 = apser(a0, b0, x0, eps);
 #ifdef DEBUG_bratio
-	    REprintf("  a0 small -> w1 := apser(*) = %15g\n", *w1);
+	    REprintf("  a0 small -> w1 := apser(*) = %.15g\n", *w1);
 #endif
 	    goto L_end_from_w1;
+	}
+	Rboolean did_bup = FALSE;
 	if (max(a0,b0) > 1.0) { /* L20:  min(a,b) <= 1 < max(a,b)  */
 #ifdef DEBUG_bratio
-	    REprintf(" L20:  min(a,b) <= 1 < max(a,b); ");
+	    REprintf("\n L20:  min(a,b) <= 1 < max(a,b); ");
 #endif
 	    if (b0 <= 1.0) goto L_w_bpser;
 	    if (x0 >= 0.29) /* was 0.3, PR#13786 */	goto L_w1_bpser;
 		*w1 = 0.;
 		goto L131;
+	    }
 	} else { /*  a, b <= 1 */
 #ifdef DEBUG_bratio
-	    REprintf("  both a,b <= 1; ");
+	    REprintf("\n      both a,b <= 1; ");
 #endif
 	    if (a0 >= min(0.2, b0))	goto L_w_bpser;
 	    if (pow(x0, a0) <= 0.9) 	goto L_w_bpser;
 	    if (x0 >= 0.3)		goto L_w1_bpser;
+	}
 	n = 20; /* goto L130; */
-	*w1 = bup(b0, a0, y0, x0, n, eps);
+	*w1 = bup(b0, a0, y0, x0, n, eps, FALSE); did_bup = TRUE;
 #ifdef DEBUG_bratio
-	REprintf("  ... n=20 and *w1 := bup(*) = %19.15g; ", *w1);
+	REprintf("  ... n=20 and *w1 := bup(*) = %.15g; ", *w1);
 #endif
 	b0 += n;
     L131:
+#ifdef DEBUG_bratio
+	REprintf(" L131: bgrat(*, w1=%.15g) ", *w1);
+#endif
-	bgrat(b0, a0, y0, x0, w1, 15*eps, &ierr1);//TODO: <-- log_p=TRUE  option !!
+	bgrat(b0, a0, y0, x0, w1, 15*eps, &ierr1, FALSE);
 	if(ierr1) *ierr = 8;
 #ifdef DEBUG_bratio
-	REprintf(" L131: bgrat(*, w1) ==> %s w1=%19.15g\n",
+	REprintf(" ==> new w1=%.15g", *w1);
-		 (ierr1)? "ERROR" : "", *w1);
+	if(ierr1) REprintf(" ERROR(code=%d)\n", ierr1) ; else REprintf("\n");
 #endif
+	if(*w1 == 0.) { // "almost surely" from underflow, try more: [2013-03-04]
+// FIXME: it is even better to do this in bgrat *directly* at least for the case
+//  !did_bup, i.e., where *w1 = (0 or -Inf) on entry
+#ifdef DEBUG_bratio
+	    REprintf(" underflow? -> retrying: ");
+#endif
+	    if(did_bup) { // re-do that part on log scale:
+		*w1 = bup(b0-n, a0, y0, x0, n, eps, TRUE);
+	    }
+	    else *w1 = ML_NEGINF; // = 0 on log-scale
+	    bgrat(b0, a0, y0, x0, w1, 15*eps, &ierr1, TRUE);
+	    if(ierr1) *ierr = 8;
+#ifdef DEBUG_bratio
+	    REprintf(" ==> new log(w1)=%.15g", *w1);
+	    if(ierr1) REprintf(" Error(code=%d)\n", ierr1) ; else REprintf("\n");
+#endif
+	    goto L_end_from_w1_log;
+	}
+	// else
 	if(*w1 < 0)
 	    MATHLIB_WARNING4("bratio(a=%g, b=%g, x=%g): bgrat() -> w1 = %g",
 			     a,b,x, *w1);
 	goto L_end_from_w1;
+    }
 	/* else if none of the above    L180: */
 	*w = basym(a0, b0, lambda, eps * 100.0, log_p);
 	*w1 = log_p ? R_Log1_Exp(*w) : 0.5 - *w + 0.5;
 #ifdef DEBUG_bratio
-	REprintf("  b0 >= a0 > 100; lambda <= a0 * 0.03: *w:= basym(*) =%15g\n",
+	REprintf("  b0 >= a0 > 100; lambda <= a0 * 0.03: *w:= basym(*) =%.15g\n",
 		 *w);
 #endif
-	goto L_end_after_log;
+	goto L_end;
     } /* else: a, b > 1 */
 /*            EVALUATION OF THE APPROPRIATE ALGORITHM */
 L_w_bpser: // was L100
     *w = bpser(a0, b0, x0, eps, log_p);
     *w1 = log_p ? R_Log1_Exp(*w) : 0.5 - *w + 0.5;
 #ifdef DEBUG_bratio
-    REprintf(" L_w_bpser: *w := bpser(*) = %15g\n", *w);
+    REprintf(" L_w_bpser: *w := bpser(*) = %.15g\n", *w);
 #endif
-    goto L_end_after_log;
+    goto L_end;
 L_w1_bpser:  // was L110
     *w1 = bpser(b0, a0, y0, eps, log_p);
     *w  = log_p ? R_Log1_Exp(*w1) : 0.5 - *w1 + 0.5;
 #ifdef DEBUG_bratio
-    REprintf(" L_w1_bpser: *w1 := bpser(*) = %15g\n", *w1);
+    REprintf(" L_w1_bpser: *w1 := bpser(*) = %.15g\n", *w1);
 #endif
-    goto L_end_after_log;
+    goto L_end;
 L_bfrac:
     *w = bfrac(a0, b0, x0, y0, lambda, eps * 15.0, log_p);
     *w1 = log_p ? R_Log1_Exp(*w) : 0.5 - *w + 0.5;
 #ifdef DEBUG_bratio
     REprintf(" L_bfrac: *w := bfrac(*) = %g\n", *w);
 #endif
-    goto L_end_after_log;
+    goto L_end;
 L140:
     /* b0 := fractional_part( b0 )  in (0, 1]  */
     n = (int) b0;
     b0 -= n;
     if (b0 == 0.) {
 	--n; b0 = 1.;
+    }
-    *w = bup(b0, a0, y0, x0, n, eps);
+    *w = bup(b0, a0, y0, x0, n, eps, FALSE);
 #ifdef DEBUG_bratio
-    REprintf(" L140: *w := bup(b0=%g,..) = %19.15g; ", b0, *w);
+    REprintf(" L140: *w := bup(b0=%g,..) = %.15g; ", b0, *w);
 #endif
     if(*w < DBL_MIN && log_p) { /* do not believe it; try bpser() : */
 	/*revert: */ b0 += n;
 	/* which is only valid if b0 <= 1 || b0*x0 <= 0.7 */
 	goto L_w_bpser;
+    }
     if (x0 <= 0.7) {
 	/* log_p :  TODO:  w = bup(.) + bpser(.)  -- not so easy to use log-scale */
 	*w += bpser(a0, b0, x0, eps, /* log_p = */ FALSE);
 #ifdef DEBUG_bratio
-	REprintf(" x0 <= 0.7: *w := *w + bpser(*) = %19.15g\n", *w);
+	REprintf(" x0 <= 0.7: *w := *w + bpser(*) = %.15g\n", *w);
 #endif
 	goto L_end_from_w;
+    }
     /* L150: */
     if (a0 <= 15.0) {
 	n = 20;
-	*w += bup(a0, b0, x0, y0, n, eps);
+	*w += bup(a0, b0, x0, y0, n, eps, FALSE);
 #ifdef DEBUG_bratio
-	REprintf("\n a0 <= 15: *w := *w + bup(*) = %19.15g;", *w);
+	REprintf("\n a0 <= 15: *w := *w + bup(*) = %.15g;", *w);
 #endif
 	a0 += n;
+    }
+#ifdef DEBUG_bratio
+    REprintf(" bgrat(*, w=%.15g) ", *w);
+#endif
-    bgrat(a0, b0, x0, y0, w, 15*eps, &ierr1);
+    bgrat(a0, b0, x0, y0, w, 15*eps, &ierr1, FALSE);
     if(ierr1) *ierr = 8;
 #ifdef DEBUG_bratio
-    REprintf(" bgrat(*) ==> %s w=%g\n",
+    REprintf("==> new w=%.15g", *w);
-		 (ierr1)? "error" : "", *w);
+    if(ierr1) REprintf(" Error(code=%d)\n", ierr1) ; else REprintf("\n");
 #endif
     goto L_end_from_w;
 /* TERMINATION OF THE PROCEDURE */
 	*w1 = log1p(-*w);
 	*w  = log(*w);
     } else {
 	*w1 = 0.5 - *w + 0.5;
+    }
-    goto L_end_after_log;
+    goto L_end;
 L_end_from_w1:
     if(log_p) {
 	*w  = log1p(-*w1);
 	*w1 = log(*w1);
     } else {
 	*w = 0.5 - *w1 + 0.5;
+    }
+    goto L_end;
+L_end_from_w1_log:
+    // *w1 = log(w1) already; w = 1 - w1  ==> log(w) = log(1 - w1) = log(1 - exp(*w1))
+    if(log_p) {
+	*w = R_Log1_Exp(*w1);
+    } else {
+	*w  = /* 1 - exp(*w1) */ -expm1(*w1);
+	*w1 = exp(*w1);
+    }
+    goto L_end;
-L_end_after_log:
+L_end:
     if (do_swap) { /* swap */
 	double t = *w; *w = *w1; *w1 = t;
+    }
     return;
     } else
 	ans *= a * sum + 1.0;
     return ans;
 } /* bpser */
-static double bup(double a, double b, double x, double y, int n, double eps)
+static double bup(double a, double b, double x, double y, int n, double eps,
+		  int give_log)
+{
 /* ----------------------------------------------------------------------- */
 /*     EVALUATION OF I_x(A,B) - I_x(A+N,B) WHERE N IS A POSITIVE INT. */
 /*     EPS IS THE TOLERANCE USED. */
 /* ----------------------------------------------------------------------- */
 	mu = 0;
 	d = 1.0;
+    }
     /* L10: */
+    ret_val = log_p
+	? brcmp1(mu, a, b, x, y, TRUE) - log(a)
-    ret_val = brcmp1(mu, a, b, x, y, FALSE) / a;
+	: brcmp1(mu, a, b, x, y, FALSE)  / a;
-    if (n == 1 || ret_val == 0.0)
+    if (n == 1 || (log_p && ret_val == ML_NEGINF) || (!log_p && ret_val == 0.))
 	return ret_val;
     int nm1 = n - 1;
     double w = d;
 		break;
+	}
+    }
     // L50: TERMINATE THE PROCEDURE
+    if(give_log) {
+	ret_val += log(w);
+    } else
-    ret_val *= w;
+	ret_val *= w;
     return ret_val;
 } /* bup */
 static double bfrac(double a, double b, double x, double y, double lambda,
 		    double eps, int log_p)
 	       ? -M_LN_SQRT_2PI + .5*log(b * x0) + z - bcorr(a,b)
 	       : const__ * sqrt(b * x0) * z * exp(-bcorr(a, b)));
+    }
 } /* brcomp */
+// called only once from  bup(),  as   r = brcmp1(mu, a, b, x, y, FALSE) / a;
+//                        -----
 static double brcmp1(int mu, double a, double b, double x, double y, int give_log)
+{
 /* -----------------------------------------------------------------------
  *          Evaluation of    exp(mu) * x^a * y^b / beta(a,b)
  * ----------------------------------------------------------------------- */
 	b0 = max(a,b);
 	if (b0 >= 8.0) {
 	/* L80:                  ALGORITHM FOR b0 >= 8 */
 	    u = gamln1(a0) + algdiv(a0, b0);
 #ifdef DEBUG_bratio
-	    REprintf(" brcmp1(mu,a,b,*): a0 < 1, b0 >= 8;  z=%19.15g\n", z);
+	    REprintf(" brcmp1(mu,a,b,*): a0 < 1, b0 >= 8;  z=%.15g\n", z);
 #endif
 	    return give_log
 		? log(a0) + esum(mu, z - u, TRUE)
 		:     a0  * esum(mu, z - u, FALSE);
 	    // L50:
 	    c = give_log
 		? log1p(gam1(a)) + log1p(gam1(b)) - log(z)
 		: (gam1(a) + 1.0) * (gam1(b) + 1.0) / z;
 #ifdef DEBUG_bratio
-	    REprintf(" brcmp1(mu,a,b,*): a0 < 1, b0 <= 1;  c=%19.15g\n", c);
+	    REprintf(" brcmp1(mu,a,b,*): a0 < 1, b0 <= 1;  c=%.15g\n", c);
 #endif
 	    return give_log
 		? ans + log(a0) + c - log1p(a0 / b0)
 		: ans * (a0 * c) / (a0 / b0 + 1.0);
+	}
 	    t = (gam1(apb - 1.) + 1.0) / apb;
 	} else {
 	    t = gam1(apb) + 1.0;
+	}
 #ifdef DEBUG_bratio
-	REprintf(" brcmp1(mu,a,b,*): a0 < 1 < b0 < 8;  t=%19.15g\n", t);
+	REprintf(" brcmp1(mu,a,b,*): a0 < 1 < b0 < 8;  t=%.15g\n", t);
 #endif
 	// L72:
 	return give_log
 	    ? log(a0)+ esum(mu, z, TRUE) + log1p(gam1(b0)) - log(t) // TODO? log(t) = log1p(..)
 	    :     a0 * esum(mu, z, FALSE) * (gam1(b0) + 1.0) / t;
 	    lambda = a - (a + b) * x;
+	}
 	double lx0 = -log1p(b/a); // in both cases
 #ifdef DEBUG_bratio
-	REprintf(" brcmp1(mu,a,b,*): a,b >= 8;  x0=%19.15g, lx0=log(x0)=%19.15g\n", x0, lx0);
+	REprintf(" brcmp1(mu,a,b,*): a,b >= 8;	x0=%.15g, lx0=log(x0)=%.15g\n",
+		 x0, lx0);
 #endif
 	// L110:
 	double e = -lambda / a;
 	if (fabs(e) > 0.6) {
 	    // L111:
+    }
 } /* brcmp1 */
 static void bgrat(double a, double b, double x, double y, double *w,
-		  double eps, int *ierr)
+		  double eps, int *ierr, Rboolean log_w)
+{
 /* -----------------------------------------------------------------------
 *     Asymptotic Expansion for I_x(a,b)  when a is larger than b.
 *     Compute   w := w + I_x(a,b)
 *     It is assumed a >= 15 and b <= 1.
 *     eps is the tolerance used.
 *     ierr is a variable that reports the status of the results.
+*
+* TODO:  if(log_w) ,  *w  itself must be in log-space;
+*     compute   w := w + I_x(a,b)  but return *w = log(w):
+*          *w := log(exp(*w) + I_x(a,b)) = logspace_add(*w, log( I_x(a,b) )
 * ----------------------------------------------------------------------- */
 #define n_terms_bgrat 30
     double c[n_terms_bgrat], d[n_terms_bgrat];
     double bm1 = b - 0.5 - 0.5,
 	    a,b,x,y);
 	/* L_Error:    THE EXPANSION CANNOT BE COMPUTED */ *ierr = 1; return;
+    }
 /*                 COMPUTATION OF THE EXPANSION */
+    double
+	/* r1 = b * (gam1(b) + 1.0) * exp(b * log(z)),// = b/gamma(b+1) z^b = z^b / gamma(b)
+	 * set r := exp(-z) * z^b / gamma(b) ;
+	 *          gam1(b) = 1/gamma(b+1) - 1 , b in [-1/2, 3/2] */
+	// exp(a*lnx) underflows for large (a * lnx); e.g. large a ==> using log_r := log(r):
+	// r = r1 * exp(a * lnx) * exp(bm1 * 0.5 * lnx);
+	// log(r)=log(b) + log1p(gam1(b)) + b * log(z) + (a * lnx) + (bm1 * 0.5 * lnx),
+	log_r = log(b) + log1p(gam1(b)) + b * log(z) + nu * lnx,
+	// FIXME work with  log_u = log(u)  also when log_p=FALSE  (??)
+	// u is 'factored out' from the expansion {and multiplied back, at the end}:
+	log_u = log_r - (algdiv(b, a) + b * log(nu)),// algdiv(b,a) = log(gamma(a)/gamma(a+b))
+	/* u = (log_p) ? log_r - u : exp(log_r-u); // =: M  in (9.2) of {reference above} */
+	/* u = algdiv(b, a) + b * log(nu);// algdiv(b,a) = log(gamma(a)/gamma(a+b)) */
+	// u = (log_p) ? log_u : exp(log_u); // =: M  in (9.2) of {reference above}
+	u = exp(log_u);
-    double r, u, p, q;
-    // set r := exp(-z) * z^b / Gamma(b) : gam1(b) = 1/gamma(b+1) - 1 , b in [-1/2, 3/2]
-    r = b * (gam1(b) + 1.0) * exp(b * log(z));// = b/gamma(b+1) z^b =  z^b / gamma(b)
-    r = r * exp(a * lnx) * exp(bm1 * 0.5 * lnx);
-    // u is 'factored out' from the expansion {and multiplied back, at the end}:
-    u = algdiv(b, a) + b * log(nu);// algdiv(b,a) = log(gamma(a)/gamma(a+b))
-    u = r * exp(-u); // =: M  in (9.2) of {reference above}
-    if (u == 0.0) {
+    if (log_u == ML_NEGINF) {
 #ifdef DEBUG_bratio
-	REprintf(" bgrat(*): underflow u = 0 = r*exp(-u'), r = %g ", r);
+	REprintf(" bgrat(*): underflow log_u = -Inf  = log_r -u', log_r = %g ",
+		 log_r);
 #endif
 	/* L_Error:    THE EXPANSION CANNOT BE COMPUTED */ *ierr = 2; return;
+    }
+    Rboolean u_0 = (u == 0.); // underflow --> do work with log(u) == log_u !
+    double l = // := *w/u .. but with care: such that it also works when u underflows to 0:
+	log_w
+	? ((*w == ML_NEGINF) ? 0. : exp(  *w    - log_u))
-    grat1(b, z, r, &p, &q, eps);
+	: ((*w == 0.)        ? 0. : exp(log(*w) - log_u));
+#ifdef DEBUG_bratio
+    REprintf(" bgrat(a=%g, b=%g, x=%g, *)\n -> u=%g, l='w/u'=%g, ", a,b,x, u, l);
+#endif
+    double
-//  -----     |--> (p , q)  {p + q = 1}  and we only use  q/r :
+	q_r = grat_r(b, z, log_r, eps), // = q/r of former grat1(b,z, r, &p, &q)
-    double v = 0.25 / (nu * nu),
+	v = 0.25 / (nu * nu),
 	t2 = lnx * 0.25 * lnx,
-	l = *w / u,
-	j = q / r,
+	j = q_r,
 	sum = j,
 	t = 1.0, cn = 1.0, n2 = 0.;
     for (int n = 1; n <= n_terms_bgrat; ++n) {
 	double bp2n = b + n2;
 	j = (bp2n * (bp2n + 1.0) * j + (z + bp2n + 1.0) * t) * v;
 		a,b,x, fabs(dj) /(sum + l));
+    }
 /*                    ADD THE RESULTS TO W */
     *ierr = 0;
+    if(log_w) // *w is in log space already:
+	*w = logspace_add(*w, log_u + log(sum));
-    *w += u * sum;
+    else
+	*w += (u_0 ? exp(log_u + log(sum)) : u * sum);
     return;
 } /* bgrat */
+// called only from bgrat() , as   q_r = grat_r(b, z, log_r, eps)  :
-static void grat1(double a, double x, double r, double *p, double *q,
+static double grat_r(double a, double x, double log_r, double eps)
-		  double eps)
+{
 /* -----------------------------------------------------------------------
- *        Evaluation of the incomplete gamma ratio functions
+ *        Scaled complement of incomplete gamma ratio function
- *                      P(a,x) and Q(a,x)
+ *                   grat_r(a,x,r) :=  Q(a,x) / r
+ * where
- *     P(a,x) = pgamma(x,a)  ---  Q(a,x) = pgamma(x,a, lower.tail=FALSE)
+ *               Q(a,x) = pgamma(x,a, lower.tail=FALSE)
+ *     and            r = e^(-x)* x^a / Gamma(a) ==  exp(log_r)
+ *
  *     It is assumed that a <= 1.  eps is the tolerance to be used.
- *     the input argument r has the value  r = e^(-x)* x^a / Gamma(a).
  * ----------------------------------------------------------------------- */
-// called only from bgrat() , as
-//  grat1(b, z, r, &p, &q, eps);  .. |-->  q  and we only use  q/r :
     if (a * x == 0.0) { /* L130: */
 	if (x <= a) {
-	    /* L100: */ *p = 0.0; *q = 1.0;
+	    /* L100: */ return exp(-log_r);
 	} else {
-	    /* L110:*/ *p = 1.0; *q = 0.0;
+	    /* L110:*/  return 0.;
+	}
+    }
-    else if (a == 0.5) { //--- FIXME not always good if (log_p): underflow!
+    else if (a == 0.5) { // e.g. when called from pt()
 	/* L120: */
 	if (x < 0.25) {
-	    *p = erf__(sqrt(x));
+	    double p = erf__(sqrt(x));
-	    *q = 0.5 - *p + 0.5;
 #ifdef DEBUG_bratio
-	    REprintf(" grat1(a=%g, x=%g ..)): a=1/2 --> p=erf__(.)= %g\n", a,x, *p);
+	    REprintf(" grat_r(a=%g, x=%g ..)): a=1/2 --> p=erf__(.)= %g\n", a,x, p);
 #endif
-	} else {
+	    return (0.5 - p + 0.5)*exp(-log_r);
+        } else { // 2013-02-27: improvement for "large" x: direct computation of q/r:
-	    *q = erfc1(0, sqrt(x));
+	    double sx = sqrt(x),
-	    *p = 0.5 - *q + 0.5;
+		q_r = erfc1(1, sx)/sx * M_SQRT_PI;
 #ifdef DEBUG_bratio
-	    REprintf(" grat1(a=%g, x=%g ..)): a=1/2 --> q=erfc1(..)= %g\n", a,x, *q);
+	    REprintf(" grat_r(a=%g, x=%g ..)): a=1/2 --> q_r=erfc1(..)/r= %g\n",
+		     a,x, q_r);
 #endif
+	    return q_r;
+	}
     } else if (x < 1.1) { /* L10:  Taylor series for  P(a,x)/x^a */
 	double an = 3.,
 	    t = c / (a + an);
 	    sum += t;
 	} while (fabs(t) > tol);
 #ifdef DEBUG_bratio
-	REprintf(" grat1(a=%g, x=%g, r=%g): sum=%g; Taylor used %.0f terms\n",
+	REprintf(" grat_r(a=%g, x=%g, log_r=%g): sum=%g; Taylor w/ %.0f terms",
-		 a,x,r, sum, an-3.);
+		 a,x,log_r, sum, an-3.);
 #endif
 	double j = a * x * ((sum/6. - 0.5/(a + 2.)) * x + 1./(a + 1.)),
 	    z = a * log(x),
 	    h = gam1(a),
 	    g = h + 1.0;
 	if ((x >= 0.25 && (a < x / 2.59)) || (z > -0.13394)) {
 	    // L40:
-	    double l = rexpm1(z);
+	    double l = rexpm1(z),
-	    *q = ((l + 0.5 + 0.5) * j - l) * g - h;
+		q = ((l + 0.5 + 0.5) * j - l) * g - h;
-	    if (*q < 0.0) {
+	    if (q <= 0.0) {
+#ifdef DEBUG_bratio
+		REprintf(" => q_r= 0.\n");
+#endif
-		/* L110:*/ *p = 1.0; *q = 0.0;
+		/* L110:*/ return 0.;
 	    } else {
+#ifdef DEBUG_bratio
+		REprintf(" => q_r=%.15g\n", q * exp(-log_r));
+#endif
-		*p = 0.5 - *q + 0.5;
+		return q * exp(-log_r);
+	    }
 	} else {
-	    *p = exp(z) * g * (0.5 - j + 0.5);
+	    double p = exp(z) * g * (0.5 - j + 0.5);
+#ifdef DEBUG_bratio
+	    REprintf(" => q_r=%.15g\n", (0.5 - p + 0.5) * exp(-log_r));
+#endif
-	    *q = 0.5 - *p + 0.5;
+	    return /* q/r = */ (0.5 - p + 0.5) * exp(-log_r);
+	}
     } else {
 	/* L50: ----  (x >= 1.1)  ---- Continued Fraction Expansion */
 	    b2n = b2n_1 + c_a * b2n;
 	    an0 = a2n / b2n;
 	} while (fabs(an0 - am0) >= eps * an0);
 #ifdef DEBUG_bratio
-	REprintf(" grat1(a=%g, x=%g, r=%g): Cont.Fraction used %.0f terms\n",
+	REprintf(" grat_r(a=%g, x=%g, log_r=%g): Cont.frac. %.0f terms => q_r=%.15g\n",
-		 a,x,r, c-1.);
+		 a,x, log_r, c-1., an0);
 #endif
-	*q = r * an0;
+	return /* q/r = (r * an0)/r = */ an0;
-	*p = 0.5 - *q + 0.5;
+    }
+} /* grat_r */
-    return;
-} /* grat1 */
 static double basym(double a, double b, double lambda, double eps, int log_p)
+{
 /* ----------------------------------------------------------------------- */
 /*     ASYMPTOTIC EXPANSION FOR I_x(A,B) FOR LARGE A AND B. */

The R Project SVN R

(root)/trunk/src/nmath/toms708.c – Rev 62109 → 62123