round(.): fix PR#17668

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

doc/NEWS.Rd

@@ -377,6 +377,10 @@
       argument which should fix behaviour during daylight savings
       changeover days, fixing \PR{16764}, thanks to proposals and analysis
       by Johannes Ranke and Kirill Müller.
+
+      \item \code{round()} correctly rounds \dQuote{to even} in more cases,
+      thanks to increased internal accuracy via code simplification
+      proposed by Adam Wheeler in his report (\PR{17668}).
     }
   }
 }

src/nmath/fround.c

@@ -1,7 +1,7 @@
 /*
  *  Mathlib : A C Library of Special Functions
+ *  Copyright (C) 2000-2019 The R Core Team
  *  Copyright (C) 1998 Ross Ihaka
- *  Copyright (C) 2000-2018 The R Core Team
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -35,8 +35,6 @@ double fround(double x, double digits) {
 #define MAX_DIGITS DBL_MAX_10_EXP
     /* = 308 (IEEE); was till R 0.99: (DBL_DIG - 1) */
     /* Note that large digits make sense for very small numbers */
-    LDOUBLE pow10, sgn, intx;
-    int dig;
 
     if (ISNAN(x) || ISNAN(digits))
 	return x + digits;
@@ -46,20 +44,26 @@ double fround(double x, double digits) {
     else if(digits == ML_NEGINF) return 0.0;
 
     if (digits > MAX_DIGITS) digits = MAX_DIGITS;
-    dig = (int)floor(digits + 0.5);
+
+    int dig = (int)floor(digits + 0.5);
+    LDOUBLE sgn = +1.;
     if(x < 0.) {
 	sgn = -1.;
 	x = -x;
-    } else
-	sgn = 1.;
+    }
     if (dig == 0) {
 	return (double)(sgn * nearbyint(x));
     } else if (dig > 0) {
-        pow10 = R_pow_di(10., dig);
-	intx = floor(x);
+	LDOUBLE pow10 = R_pow_di(10., dig)
+#ifdef R_3_AND_OLDER
+	    , intx = floor(x);
 	return (double)(sgn * (intx + nearbyint((double)((x-intx) * pow10)) / pow10));
+#else
+	;
+	return (double)(sgn * (nearbyint((double)(x * pow10)) / pow10));
+#endif
     } else {
-        pow10 = R_pow_di(10., -dig);
+	LDOUBLE pow10 = R_pow_di(10., -dig);
         return (double)(sgn * nearbyint((double)(x/pow10)) * pow10);
     }
 }

tests/reg-tests-1d.R

@@ -3463,6 +3463,32 @@ stopifnot(identical(w0[sel], w1[sel]), identical(w0[sel], wII[sel]))
 ## Inf-Inf  etc broken in paired case in R <= 3.6.x
 
 
+## round(x, n) "to even" fails in some cases -- PR#17668
+dd <- 0:12
+x55 <- 55 + as.numeric(vapply(dd+1, function(k) paste0(".", strrep("5",k)), ""))
+rnd.x <- vapply(dd+1L, function(k) round(x55[k], dd[k]), 1.1)
+noquote(formatC(cbind(x55, dd, rnd.x), w=1, digits=15))
+stopifnot(exprs = {
+      print (   rnd.x - x55) > 0
+      all.equal(rnd.x - x55, 5 * 10^-(dd+1), tol = 1e-11) # see diff. of 6.8e-13
+})
+## more than half of the above were rounded *down* in R <= 3.6.x
+## This must work, too, the range of 'e' depending on 'd'
+EE <- c(-307, -300, -250, -200,-100,-50, -20, -10, -2:2,
+        10, 20, 50, 100, 200, 250, 290:307)
+for(d in 0:16) { cat("digits 'd' = ", d, ": ")
+    for(e in EE[EE+d <= 308]) {
+        f <- 10^e
+        cat(".")
+        stopifnot(all.equal(tolerance = if(d < 14) 1e-15
+                                        else if(d == 14) 1e-14 else 1e-13,
+                            round(pi/f, e + d) * f,
+                            round(pi, d)))
+    };cat("\n")
+}
+## (2nd part: continued working)
+
+
 
 ## keep at end
 rbind(last =  proc.time() - .pt,