signif(1***e308, dig) no longer truncates wrongly

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

Author: maechler

doc/NEWS.Rd

@@ -118,6 +118,9 @@
       \code{dbinom()}, \code{dnbinom()}, and via \code{dbinom_raw()}
       potentially \code{dgeom()}, \code{dhyper()}, \code{dbeta()}, and
       \code{df()}.
+
+      \item \code{signif(1.**e308, digits)} no longer truncates unnecessarily
+      (but still to prevent overflow to \code{Inf}), fixing \PR{18889}.
     }
   }
 }

src/nmath/fprec.c

@@ -1,6 +1,6 @@
 /*
  *  Mathlib : A C Library of Special Functions
- *  Copyright (C) 2000-2019 The R Core Team
+ *  Copyright (C) 2000-2025 The R Core Team
  *  Copyright (C) 1998 Ross Ihaka
  *
  *  This program is free software; you can redistribute it and/or modify
@@ -57,38 +57,37 @@
 // R's  signif(x, digits)   via   Math2(args, fprec) in  ../main/arithmetic.c :
 double fprec(double x, double digits)
 {
-    double l10, pow10, sgn, p10, P10;
-    int e10, e2, do_round, dig;
-    // Max.expon. of 10 (w/o denormalizing or overflow; = R's  trunc( log10(.Machine$double.xmax) )
-    const static int max10e = (int) DBL_MAX_10_EXP; // == 308 ("IEEE")
-
     if (ISNAN(x) || ISNAN(digits))
 	return x + digits;
     if (!R_FINITE(x)) return x;
     if (!R_FINITE(digits)) {
-	if(digits > 0.0) return x;
-	else digits = 1.0;
+	if(digits > 0.) return x;
+	else digits = 1.;
     }
     if(x == 0) return x;
-    dig = (int)round(digits);
+
+    int dig = (int)round(digits);
     if (dig > MAX_DIGITS) {
 	return x;
     } else if (dig < 1)
 	dig = 1;
 
-    sgn = 1.0;
+    double sgn = 1.;
     if(x < 0.0) {
 	sgn = -sgn;
 	x = -x;
     }
-    l10 = log10(x);
-    e10 = (int)(dig-1-floor(l10));
+    double l10 = log10(x);
+    int e10 = dig-1 - (int)floor(l10);
+    // Max.expon. of 10 (w/o denormalizing or overflow; = R's trunc( log10(.Machine$double.xmax) ):
+    const static int max10e = (int) DBL_MAX_10_EXP; // == 308 ("IEEE")
     if(fabs(l10) < max10e - 2) {
-	p10 = 1.0;
+	double p10 = 1.;
 	if(e10 > max10e) { /* numbers less than 10^(dig-1) * 1e-308 */
 	    p10 =  R_pow_di(10., e10-max10e);
 	    e10 = max10e;
 	}
+	double pow10;
 	if(e10 > 0) { /* Try always to have pow >= 1
 			 and so exactly representable */
 	    pow10 = R_pow_di(10., e10);
@@ -98,10 +97,13 @@ double fprec(double x, double digits)
 	    return(sgn*(nearbyint((x/pow10))*pow10));
 	}
     } else { /* -- LARGE or small -- */
-	do_round = max10e - l10	 >= R_pow_di(10., -dig);
-	e2 = dig + ((e10>0)? 1 : -1) * MAX_DIGITS;
-	p10 = R_pow_di(10., e2);	x *= p10;
-	P10 = R_pow_di(10., e10-e2);	x *= P10;
+	bool do_round = log10(DBL_MAX) - l10 >= R_pow_di(10., -dig); /* e.g. signif(1.09e308, 2) */
+	int e2 = dig + ((e10 > 0)? 1 : -1) * MAX_DIGITS;
+	double
+	    p10 = R_pow_di(10., e2),
+	    P10 = R_pow_di(10., e10-e2);
+	x *= p10;
+	x *= P10;
 	/*-- p10 * P10 = 10 ^ e10 */
 	if(do_round) x += 0.5;
 	x = floor(x) / p10;

tests/reg-tests-1e.R

@@ -1949,6 +1949,11 @@ r <- return; assertWarnV(environment(r) <- baseenv())
 ## then an error for about one day; now is deprecated (and no longer mutating).
 
 
+## signif(<LRG>, dig) -- PR#18889
+stopifnot(all.equal(1.1e308,  signif(1.06e308,   2)),
+          all.equal(1.01e308, signif(1.0055e308, 3)))
+
+
 
 ## keep at end
 rbind(last =  proc.time() - .pt,