floor(NaN) problem fixed in massdist.c (PR#291)

[maechler@stat.math.ethz.ch]

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