Re: [Rd] 0 ^ NaN == Inf, why?

From: Stephen Milborrow <milbo_at_sonic.net>
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

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

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.

list of date sections of archive