R> round(100000/3, -2) - 33300
[1] -7.275958e-12
I would have hoped for 0. The problem seems to be specifically for negative
"digits". This is in 1.3.1 on Solaris 2.6 (full description at
bottom).
[Apologies for making everyone read this 3 times; my bug.report() burped.]
Peter Dalgaard <p.dalgaard@biostat.ku.dk> says the problem is in fround.c,
which reads in part:
pow10 = R_pow_di(10., dig)
...
return sgn * (intx + R_rint(x * pow10) / pow10);
and (says Peter),> the problem is most likely that pow10 (== 1/100) has no finite binary
> representation. Patches are welcome...
I'm not up to patching C code, but I will say that a workaround is:
R> round(round(100000/3, -2)) - 33300
[1] 0
and a smarter C-programmer could probably implement this concept:
if digits < 0, force the final result to an integer.
-- David Brahm (brahm@alum.mit.edu)
Version:
platform = sparc-sun-solaris2.6
arch = sparc
os = solaris2.6
system = sparc, solaris2.6
status =
major = 1
minor = 3.0
year = 2001
month = 06
day = 22
language = R
Search Path:
.GlobalEnv, package:misc, package:io, package:arrays, package:ls1, package:db,
package:ts, package:ctest, Autoloads, package:base
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
I reported this bug:> R> round(100000/3, -2) - 33300 > [1] -7.275958e-12Peter 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._