В чем разница между R_pow() в R и pow() в libc?

Руководство по Writing R Extensions, разд. 6.7.3., утверждает, что функция R API, объявленная как double R_pow (double x, double y), вычисляет x^y:

[...] используя проверки R_FINITE и возвращая правильный результат (такой же, как R) для случаев, когда x, y или i равны 0 или отсутствуют, или бесконечны, или NaN.

Однако я не могу найти такие x и y, для которых функция pow() из библиотеки C дает неверные результаты. Я пробовал разные случаи, например 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

Может быть, вы, ребята, предоставите мне какой-нибудь неправильный пример.

Кстати, я использую gcc 4.8.2 с glibc 2.18 (Fedora 20, x86_64). Исходный код R_pow() найдите R_pow здесь.


c r rcpp glibc pow
person gagolews    schedule 02.07.2014    source источник
comment
Обеспечивает ли pow() одинаковую обработку особых случаев на всех платформах?   -  person David Heffernan    schedule 03.07.2014
comment
@DavidHeffernan: я так не думаю. Вот исходный код для glibc. По сути, это exp(x*log(y)) или __ieee754_powf или любой из них + некоторые проверки, такие как здесь.   -  person gagolews    schedule 03.07.2014
comment
man 3 pow также дает исчерпывающее объяснение всех потенциальных условий сингулярности. Было бы трудно найти набор x или y, который не был бы защищен от них.   -  person David C. Rankin    schedule 03.07.2014
comment
R не всегда связан с glibc, не так ли?   -  person David Heffernan    schedule 03.07.2014
comment
Эта таблица может помочь, она должна быть такой же, как и для C.   -  person Shafik Yaghmour    schedule 03.07.2014
comment
@DavidHeffernan: Нет, это не обязательно связано с glibc (man 3 pow на моем старом Solaris ничего не говорит о граничных условиях). Хорошо, спасибо, ребята, я думаю, теперь у меня есть ответ.   -  person gagolews    schedule 03.07.2014

Ответы (1)


arrow_upward
4
arrow_downward

Возможно, проще всего посмотреть исходный код из которого я скопировал настоящую функцию:

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}
}

который показывает, в каких случаях это сворачивается до pow(x, y).

person Dirk Eddelbuettel    schedule 05.07.2014