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

[Stephen Milborrow]

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