[Rd] round() doesn't (PR#1139)

[brahm@alum.mit.edu]

Hi, all,

In PR#1138 and PR#1139, I pointed out that round() with digits<0 would not
actually round to an integer.  Brian D. Ripley (hereafter BDR) fixed it in
R-1.4.0, but the fix introduced a new error, PR#1254 (caught by Ole
Christensen).  BDR (I'm guessing) fixed the fix; here's the relevant line from
fround.c:

  if(dig <= 0 && fabs(res) < 1e9) res = (int)floor(res + 0.5);

This solution casts the result as int if it's under a billion, but it
re-introduces the original problem for large numbers:

  R-patched 01/13/2002> round(1e15, -2) - 1e15
                        [1] -0.25

So I'd ask to re-open PR#1139, and suggest another look at the "fround" C
function that I posted (in a followup) on 10/24/2001 (copied below).  I
substituted it into src/nmath/fround.c, successfully re-built R, and it works.
It handles the 3 cases dig=0, dig>0, and dig<0 separately, avoids the roundoff
errors that were occurring for dig<0, and avoids a call to R_pow_di altogether
when dig=0.  It gets the correct answer "0" in the above example.

Here are a few explicit test cases; the last one fails in R-patched 01/13/2002,
but all work with my "fround".
  all.equal(round(1234.5678,  2),         1234.57)
  all.equal(round(1234.5678,  0),         1235)
  all.equal(round(1234.5678, -2),         1200)
  all.equal(round(100000/3,  -2) - 33300, 0)
  all.equal(round(1e15,      -2) - 1e15,  0)  # Fails in R-patched

Thanks to all for a great program, and congratulations on the release of 1.4.0!

---- begin C code for fround.c: ----

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 */
    double pow10, sgn, intx;
    int dig;

#ifdef IEEE_754
    if (ISNAN(x) || ISNAN(digits))
	return x + digits;
    if(!R_FINITE(x)) return x;
#endif

    if (digits > MAX_DIGITS)
	digits = MAX_DIGITS;
    dig = (int)floor(digits + 0.5);
    if(x < 0.) {
	sgn = -1.;
	x = -x;
    } else
	sgn = 1.;
    if (dig == 0) {
      return sgn * R_rint(x);
    } else if (dig > 0) {
        pow10 = R_pow_di(10., dig);
	intx = floor(x);
	return sgn * (intx + R_rint((x-intx) * pow10) / pow10);
    } else {
        pow10 = R_pow_di(10., -dig);
        return sgn * R_rint(x/pow10) * pow10;
    }
}
---- end C code ----
                              -- David Brahm (brahm@alum.mit.edu)

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._