From: Stephen Milborrow <milbo_at_sonic.net>

Date: Sun, 26 Oct 2008 14:35:50 +0000

}

R-devel_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Sun 26 Oct 2008 - 14:37:34 GMT

Date: Sun, 26 Oct 2008 14:35:50 +0000

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" <jmc_at_r-project.org> To: "Stephen Milborrow" <milbo_at_sonic.net> Cc: <r-devel_at_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" <jmc_at_r-project.org> To: "Stephen Milborrow" <milbo_at_sonic.net> Cc: <r-devel_at_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_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Sun 26 Oct 2008 - 14:37:34 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Sun 26 Oct 2008 - 15:30:27 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel.
Please read the posting
guide before posting to the list.
*