ntakebay@bio.indiana.edu
1999-Oct-08 02:57 UTC
floor(NaN) problem fixed in massdist.c (PR#291)
Full_Name: Naoki Takebayashi
Version: 0.65.0+R-release.diff (Oct 6, 1999)
OS: Linux/Alpha
Submission from: (NULL) (129.79.224.171)
This will fix the "problem 2 (crash in fft)" in Bug ID #277
On Linux/Alpha, make check failed because R could not handle the following
example in base-Ex.R
##___ Examples ___:
# The Old Faithful geyser data
data(faithful)
:
:
## Missing values:
x <- xx <- faithful$eruptions
x[i.out <- sample(length(x), 10)] <- NA
doR <- density(x, bw=0.15, na.rm = TRUE)
doN <- density(x, bw=0.15, na.rm = FALSE)
### it fails at this statement ###
lines(doR, col="blue")
lines(doN, col="red")
points(xx[i.out], rep(.01,10))
The following statement in density() was not getting the right output:
y <- .C("massdist", x = as.double(x), nx = N, xlo =
as.double(lo),
xhi = as.double(up), y = double(2 * n), ny = as.integer(n),
NAOK = has.na, PACKAGE = "base")$y
y[1] and y[2] should be 0 and 0 but both of them became NaN after this
statement.
So I looked at src/appl/masdist.c with a debugger. The bug is caused
by floor(). On Alpha/Linux, floot(NaN) returns 0. x[i] in the
following code is the list passed from the above statement in
density(). The array x contains a couple of NaN. When x[i] is NaN,
xpos become NaN. Then ix = floor(xpos) returns an unexpected result.
On linux/i386, ix become -2147483648 when xpos = NaN(0x80000000007a2).
However ix become 0 on Linux/Alpha when xpos=NaN. When xpos = NaN, it
shouldn't go into the following "if" statement. But it goes
inside on
Alpha/Linux (because ixmin is 0), and set y[0] and y[1] to NaN.
Following is extracted from massdist.c
for(i=0 ; i<*nx ; i++) {
xpos = (x[i] - *xlow) / xdelta;
ix = floor(xpos);
fx = xpos - ix;
if(ixmin <= ix && ix <= ixmax) {
y[ix] += (1.0 - fx);
y[ix + 1] += fx;
}
else if(ix == -1) {
y[0] += fx;
}
else if(ix == ixmax + 1) {
y[ixmax + 1] += (1.0 - fx);
}
}
So I think the statement
if(ixmin <= ix && ix <= ixmax) {
should be changed to
if((! isnan(ix)) && ixmin <= ix && ix <= ixmax) {
or
everything inside of this for loop should be encolsed with
for(i=0 ; i<*nx ; i++) {
if ( ! isnan( x[i] ) ) {
:
:
}
}
A patch to fix this bug is included.
--- R-0.65.0/src/appl/massdist.c.orig Thu Oct 7 21:26:53 1999
+++ R-0.65.0/src/appl/massdist.c Thu Oct 7 21:47:55 1999
@@ -38,6 +38,7 @@
y[i] = 0;
for(i=0 ; i<*nx ; i++) {
+ if ( ! isnan( x[i] ) ) {
xpos = (x[i] - *xlow) / xdelta;
ix = floor(xpos);
fx = xpos - ix;
@@ -51,6 +52,7 @@
else if(ix == ixmax + 1) {
y[ixmax + 1] += (1.0 - fx);
}
+ }
}
for(i=0 ; i<*ny ; i++)
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
maechler@stat.math.ethz.ch
1999-Oct-08 16:54 UTC
floor(NaN) problem fixed in massdist.c (PR#291)
Dear Naoki, Thank you very much, for your detailed bug report and your proposed patch.>> This will fix the "problem 2 (crash in fft)" in Bug ID #277 >> >> On Linux/Alpha, make check failed because R could not handle the following >> example in base-Ex.R >> >> ##___ Examples ___: >> >> # The Old Faithful geyser data >> data(faithful) >> : >> : >> ## Missing values: >> x <- xx <- faithful$eruptions >> x[i.out <- sample(length(x), 10)] <- NA >> doR <- density(x, bw=0.15, na.rm = TRUE) >> doN <- density(x, bw=0.15, na.rm = FALSE) >> ### it fails at this statement ### >> lines(doR, col="blue") >> ...>> The following statement in density() was not getting the right output: >> >> y <- .C("massdist", x = as.double(x), nx = N, xlo = as.double(lo), >> xhi = as.double(up), y = double(2 * n), ny = as.integer(n), >> NAOK = has.na, PACKAGE = "base")$y >> >> y[1] and y[2] should be 0 and 0 but both of them became NaN after this >> statement. >> >> So I looked at src/appl/masdist.c with a debugger. The bug is caused >> by floor(). On Alpha/Linux, floor(NaN) returns 0.Thank you so much for finding this! However, this is clearly a bug of (your version of) gcc, respectively libm on Alpha/Linux. In general, we do rely on IEEE Arithmetic working correctly.>> x[i] in the >> following code is the list passed from the above statement in >> density(). The array x contains a couple of NaN. When x[i] is NaN, >> xpos become NaN. Then ix = floor(xpos) returns an unexpected result. >> On linux/i386, ix become -2147483648 when xpos = NaN(0x80000000007a2). >> However ix become 0 on Linux/Alpha when xpos=NaN. When xpos = NaN, it >> shouldn't go into the following "if" statement. But it goes inside on >> Alpha/Linux (because ixmin is 0), and set y[0] and y[1] to NaN. >> >> Following is extracted from massdist.c >> >> for(i=0 ; i<*nx ; i++) { >> xpos = (x[i] - *xlow) / xdelta; >> ix = floor(xpos); >> fx = xpos - ix; >> if(ixmin <= ix && ix <= ixmax) { >> y[ix] += (1.0 - fx); >> y[ix + 1] += fx; >> } >> else if(ix == -1) { >> y[0] += fx; >> } >> else if(ix == ixmax + 1) { >> y[ixmax + 1] += (1.0 - fx); >> } >> } >> >> >> So I think the statement >> >> if(ixmin <= ix && ix <= ixmax) { >> >> should be changed to >> >> if((! isnan(ix)) && ixmin <= ix && ix <= ixmax) { >> >> or >> ........Instead of isnan () , it should rather be ISNAN() [an R defined macro] which is true for both NaN and NA. (note that we also try to have a working ISNAN() when there's no isnan() available from C (note that ANSI or POSIX C unfortunately do NOT require IEEE arithmetic!) I'll have a look at the code and your proposed patch. Thanks again! -- Martin Maechler <maechler@stat.math.ethz.ch> http://stat.ethz.ch/~maechler/ Seminar fuer Statistik, ETH-Zentrum LEO D10 Leonhardstr. 27 ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND phone: x-41-1-632-3408 fax: ...-1228 <>< -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._