tweak bd0(<large>,<small>) to prevent overflow in x/np, also fixing dbinom(*, prob=<subnormal_small>, log=TRUE)

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

Author: maechler

src/nmath/bd0.c

@@ -59,7 +59,7 @@ attribute_hidden double bd0(double x, double np)
 	    v = (x_ - n_)/(x_ + n_);
 	}
 	double s = ldexp(d, -1) * v; // was d * v
-	if(fabs(s) < DBL_MIN*0.5) return ldexp(s, 1);// CARE: assumes subnormal numbers
+	if(fabs(ldexp(s, 1)) < DBL_MIN) return ldexp(s, 1);
 	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
@@ -79,8 +79,11 @@ attribute_hidden double bd0(double x, double np)
 			 x, np, s, ej/((1000<<1)+1));
     }
     /* else:  | x - np |  is not too small */
-    return (x > np) ? x*(log(x/np) -1.) + np
-	            : x* log(x/np) + np -x;
+    /* NB: x/np |--> inf (overflow)  doesn't happen when called from rpois_raw() */
+#define lg_x_n (R_FINITE(x/np) ? log(x/np) : (log(x) - log(np)))
+    return (x > np) ? x*(lg_x_n -1.) + np
+	            : x* lg_x_n + np -x;
+#undef lg_x_n
 }
 
 

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

@@ -665,11 +665,19 @@ dbArg <- cbind(
 (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
+    all.equal(db, tolerance = 1e-5,
+	      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)))
+## all == -Inf   in R <= 4.5.0
+
+x <- c(12:20, 100*c(1,3,10,20,50)) # dbinom for very small prob (did overflow in bd0())
+(db <- dbinom(x, size=x+1, prob = 2^-1024.1, log = TRUE))
+stopifnot(
+    all.equal(c(-8515.66, -9225.44, -9935.22, -10645, -11354.8, -12064.6, -12774.4,
+                -13484.2, -14194, -70980.6, -212950, -709845, -1419700, -3549250),
+              db, tolerance = 1e-5))
+## all but the first two where -Inf in R <= 4.5.0
 
 
 ## qpois(p, *) for invalid 'p' should give NaN -- PR#16972