It may be too late to make such changes now, but would it not be simpler if
R ^ and R pow (and related functions like log) when the return type is
double simply call the standard C library routine pow? That way we would be
compatible with IEEE 754, with the fair assumption that the C library's pow
is compatible with IEEE 754 (this I believe is a fair assumption in 2008,
may not have been historically). Having special cases in the code for x==0
etc. is unnecessarily complicated, and returning NA in some cases and NaN in
others introduces further complication. But I may be missing something.
I think this is basically what you are saying with "Other things being
equal, `^` should follow the C pow() function for numeric arguments."
When I was building my just-in-time compiler for R, I looked into some
discrepancies between R and IEEE 754. My tests were not exhaustive because
the jit compiler in many cases uses the same C math routines as standard R,
so I didn't bother to check those for compatibility. But if you are
interested, the results are summarized in tables in the jit man page
www.milbo.users.sonic.net/ra/jit.html.
I would be prepared to sign up for doing a full test of incompatiblities
between R math and IEEE 754, if people think time spent doing that is worth
it.
Steve
www.milbo.users.sonic.net
----- Original Message -----
From: "John Chambers" <[EMAIL PROTECTED]>
To: "Stephen Milborrow" <[EMAIL PROTECTED]>
Cc: <r-devel@r-project.org>
Sent: Saturday, October 25, 2008 7:55 PM
Subject: Re: [Rd] 0 ^ NaN == Inf, why?
A small PS:
John Chambers wrote:
>
> Along the line, notice that both R_pow and pow give 0^0 as 1. (Just at
> a guess, C might give 0^-0 as Inf, but I don't know how to test that in
> R.)
>
I tried a little harder, and apparently the guess is wrong. It seems
that pow(0, -0) is 1 in C. Would seem better to either have pow(0,0)
and pow(0,-0) both be NaN or else 1 and Inf, but ...
> John
----- Original Message -----
From: "John Chambers" <[EMAIL PROTECTED]>
To: "Stephen Milborrow" <[EMAIL PROTECTED]>
Cc: <r-devel@r-project.org>
Subject: Re: [Rd] 0 ^ NaN == Inf, why?
Stephen Milborrow wrote:
In R, 0 ^ NaN yields Inf. I would have expected NaN or perhaps 0. Is this
behaviour intended?
Well, it certainly follows from the implementation. In the R_pow C routine,
double R_pow(double x, double y) /* = x ^ y */
{
if(x == 1. || y == 0.)
return(1.);
if(x == 0.) {
if(y > 0.) return(0.);
/* y < 0 */ return(R_PosInf);
}
It does seem, however, that NaN is the logical result here, which I think
results from changing the implementation to:
if(x == 0.) {
if(y > 0.) return(0.);
else if(y < 0) return(R_PosInf);
else return(y); /* NA or NaN, we assert */
}
Other things being equal, `^` should follow the C pow() function for numeric
arguments. After writing pow() as an R function that calls this C function:
> pow(0,NaN)
[1] NaN
> pow(0,NA)
[1] NA
> pow(0,0)
[1] 1
The second example presumably falls out because the C function returns its
second argument if that is a NaN (which a numeric NA is, in C but not in R).
The modified code in R_pow mimics that behavior.
Along the line, notice that both R_pow and pow give 0^0 as 1. (Just at a
guess, C might give 0^-0 as Inf, but I don't know how to test that in R.)
Other opinions?
John
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
