What's the difference between R's R_pow() and libc's pow()?

Question

The Writing R Extensions manual, Sec. 6.7.3., states that the R API function declared as double R_pow (double x, double y) computes x^y:

[...] using R_FINITE checks and returning the proper result (the same as R) for the cases where x, y or i are 0 or missing or infinite or NaN.

However, I cannot find such x and y for which the pow() function from the C library gives improper results. I tried various cases, like x, y beingInf,NA/NaN, integers, and so on, but I found no input data that generated the results different than those returned by ordinarypow()`.

Rcpp::evalCpp("::pow(1.124e-15, 2)", includes = "#include <cmath>") == Rcpp::evalCpp("R_pow(1.124e-15, 2)")
## [1] TRUE

Maybe you guys will provide me with some improper example.

BTW, I'm using gcc 4.8.2 with glibc 2.18 (Fedora 20, x86_64). For R_pow()'s source code search for R_pow here.

Does pow() give the same handling of special cases on all platforms? — David Heffernan, Jul 02 '14 at 20:10
@DavidHeffernan: I don't think so. Here's the source code for [glibc](https://sourceware.org/git/?p=glibc.git;a=tree;f=math). Basically, it's `exp(x*log(y))` or `__ieee754_powf` or any of them + some checks like the ones [here](https://sourceware.org/git/?p=glibc.git;a=blob;f=math/w_powf.c;h=fa9e0071c82939f8e3db87251328ab00f840a672;hb=HEAD). — gagolews, Jul 02 '14 at 20:14
`man 3 pow` also gives an exhaustive explanation of all potential singularity conditions. It would be hard to find a set of `x` or `y` that is not protected against. — David C. Rankin, Jul 02 '14 at 20:15
This [table](http://stackoverflow.com/a/19956041/1708801) may help, it should be the same as for C. — Shafik Yaghmour, Jul 02 '14 at 20:17
@DavidHeffernan: No, it isn't necessarily linked against glibc (man 3 pow on my old Solaris tells nothing about the *boundary* conditions). All right, thank you guys, I guess I have my answer now. — gagolews, Jul 02 '14 at 20:22

score 4 · Accepted Answer · answered Jul 05 '14 at 15:43

It may just be easiest to look at the source code from which I copied the actual function:

double R_pow(double x, double y) /* = x ^ y */
{
    /* squaring is the most common of the specially handled cases so
       check for it first. */
    if(y == 2.0)
        return x * x;
    if(x == 1. || y == 0.)
        return(1.);
    if(x == 0.) {
        if(y > 0.) return(0.);
        else if(y < 0) return(R_PosInf);
        else return(y); /* NA or NaN, we assert */
    }
    if (R_FINITE(x) && R_FINITE(y)) {
        /* There was a special case for y == 0.5 here, but
           gcc 4.3.0 -g -O2 mis-compiled it.  Showed up with
           100^0.5 as 3.162278, example(pbirthday) failed. */
        return pow(x, y);
    }
    if (ISNAN(x) || ISNAN(y))
        return(x + y);
    if(!R_FINITE(x)) {
        if(x > 0)               /* Inf ^ y */
            return (y < 0.)? 0. : R_PosInf;
        else {                  /* (-Inf) ^ y */
            if(R_FINITE(y) && y == floor(y)) /* (-Inf) ^ n */
                return (y < 0.) ? 0. : (myfmod(y, 2.) ? x  : -x);
        }
    }
    if(!R_FINITE(y)) {
        if(x >= 0) {
            if(y > 0)           /* y == +Inf */
                return (x >= 1) ? R_PosInf : 0.;
            else                /* y == -Inf */
                return (x < 1) ? R_PosInf : 0.;
        }
    }
    return R_NaN; // all other cases: (-Inf)^{+-Inf, non-int}; (neg)^{+-Inf}
}

which shows in which cases this collapses to pow(x, y).

What's the difference between R's R_pow() and libc's pow()?

1 Answers1