YA tweak to bd0(<large>,<large>): prevent overflow in 2*x*v

git-svn-id: https://svn.r-project.org/R/trunk@88193 00db46b3-68df-0310-9c12-caf00c1e9a41

Author: maechler

src/nmath/bd0.c

@@ -58,20 +58,20 @@ attribute_hidden double bd0(double x, double np)
 		n_ = ldexp(np,-2);
 	    v = (x_ - n_)/(x_ + n_);
 	}
-	double s = d * v;
-	if(fabs(s) < DBL_MIN) return s;
-	double ej = x * ldexp(v, 1); // = 2*x*v;
+	double s = ldexp(d, -1) * v; // was d * v
+	if(fabs(s) < DBL_MIN*0.5) return ldexp(s, 1);// CARE: assumes subnormal numbers
+	double ej = x * v; // as 2*x*v could overflow:  v > 1/2  <==> ej = 2xv > x
 	v *= v; // "v = v^2"
 	for (int j = 1; j < 1000; j++) { /* Taylor series; 1000: no infinite loop
-					as |v| < .1,  v^2000 is "zero" */
-	    ej *= v;// = 2 x v^(2j+1)
+					    as |v| < .1,  v^2000 is "zero" */
+	    ej *= v;// = x v^(2j+1)
 	    double s_ = s;
 	    s += ej/((j<<1)+1);
 	    if (s == s_) { /* last term was effectively 0 */
 #ifdef DEBUG_bd0
-		REprintf("bd0(%g, %g): T.series w/ %d terms -> bd0=%g\n", x, np, j, s);
+		REprintf("bd0(%g, %g): T.series w/ %d terms -> bd0=%g\n", x, np, j, ldexp(s, 1));
 #endif
-		return s;
+		return ldexp(s, 1); // 2*s ; as we dropped '2 *' above
 	    }
 	}
 	/* ---- the following should _never_ happen ------------ */

tests/d-p-q-r-tst-2.R

@@ -655,6 +655,23 @@ stopifnot(all.equal(tolerance = 1e-11,
                 134.598198071, 55.4909088367, 10.1000515546, 1.6625043907, 36.710476479,
                 127.476528317, 297.82558994, 600.93919587, 1195.32578926, 3368.52998705), db))
 
+## db0() tweak (against overflow of ej = 2*x*v):
+dbArg <- cbind(
+    x = c(1.465e+308, 1.4715e+308, 1.4762e+308, 1.4822e+308, 1.4869e+308, 1.4949e+308, 1.5034e+308,
+          1.5137e+308, 1.523e+308, 1.5305e+308, 1.5416e+308, 1.5486e+308, 1.5653e+308, 1.5853e+308, 1.639e+308),
+    size=c(1.6574e+308, 1.6514e+308, 1.7035e+308, 1.679e+308, 1.6531e+308, 1.6285e+308, 1.6993e+308,
+           1.6661e+308, 1.6801e+308, 1.6873e+308, 1.7506e+308, 1.7052e+308, 1.6752e+308, 1.6885e+308, 1.731e+308),
+    prob=c(27, 27, 26, 27, 28, 29, 28, 29, 29, 30, 29, 30, 32, 33, 35)/128)
+(db <- apply(dbArg, 1, \(v3) do.call(dbinom, c(v3, list(log=TRUE))))) ## was all  -Inf !
+## dput(signif(db, 6))
+stopifnot(
+    all.equal(db,
+              c(-1.73031e+308, -1.76399e+308, -1.73534e+308, -1.74653e+308, -1.7615e+308,
+                -1.79181e+308, -1.7259e+308,  -1.77688e+308, -1.77981e+308, -1.74055e+308,
+                -1.70236e+308, -1.76548e+308, -1.79598e+308, -1.79126e+308, -1.79514e+308), tolerance = 1e-5))
+## all == -Inf previously
+
+
 ## qpois(p, *) for invalid 'p' should give NaN -- PR#16972
 stopifnot(is.nan(suppressWarnings(c(qpois(c(-2,3, NaN), 3), qpois(1, 3, log.p=TRUE),
                                     qpois(.5, 0, log.p=TRUE), qpois(c(-1,pi), 0)))))