From brahm@alum.mit.edu  Wed Oct 24 19:52:45 2001
From: brahm@alum.mit.edu (brahm@alum.mit.edu)
Date: Wed, 24 Oct 2001 20:52:45 +0200 (MET DST)
Subject: [Rd] round() doesn't (PR#1139)
Message-ID: <200110241852.UAA10218@pubhealth.ku.dk>

I reported this bug:
> R> round(100000/3, -2) - 33300
>    [1] -7.275958e-12

Peter Dalgaard <p.dalgaard@biostat.ku.dk> pointed out the problem is in
nmath/fround.c and said
> Patches are welcome...
so here goes!  I've only tested it stand-alone, not compiled into the full R
system.  I handle the 3 cases dig=0, dig>0, and dig<0 separately, and note I
avoid a call to R_pow_di altogether in the common case of dig=0.  Hope this is
useful.
			 -- David Brahm (brahm@alum.mit.edu)


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;
    }
}

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._


[Rd] round() doesn't (PR#1139)
brahm@alum.mit.edu