NetBSD-Bugs archive

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index][Old Index]

standards/46388: pow() behaves inconsistently and not POSIXly correct

>Number:         46388
>Category:       standards
>Synopsis:       pow(1, inf) results in NaN or 1.0 in different situations
>Confidential:   no
>Severity:       non-critical
>Priority:       low
>Responsible:    standards-manager
>State:          open
>Class:          sw-bug
>Submitter-Id:   net
>Arrival-Date:   Sun Apr 29 14:55:00 +0000 2012
>Originator:     Peter Bex
>Release:        NetBSD 5.0.1_PATCH
System: NetBSD 5.0.1_PATCH NetBSD 5.0.1_PATCH (FROHIKE) 
#0: Tue Aug 11 22:31:56 CEST 2009 
Architecture: x86_64
Machine: amd64
        The pow() function does not behave like the specification at
        says it should for infinity (and NaN?) values.

        Worse, gcc seems to optimize calls to pow with constant values away
        and replace it with POSIXly correct values, which means the result
        may differ depending on how the code is structured and possibly also
        on which optimization level is chosen.

        See below for the offending code.

        $ cat test.c
        #include <stdio.h>
        #include <math.h>
        int main(void) {
                double m1 = 1.0, m2 = -1.0/0.0;
                printf("directly: pow(%f, %f) = %lf\n", (double)1.0, 
(double)-1.0/0.0, pow(1.0, -1.0/0.0));
                printf("indirectly: pow(%f, %f) = %lf\n", m1, m2, pow(m1, m2));

                return 0;
        $ gcc -lm test.c
        $ ./a.out
        directly: pow(1.000000, -inf) = 1.000000
        indirectly: pow(1.000000, -inf) = nan

        This is obviously inconsistent, and SuS says the second line is

        "For any value of y (including NaN), if x is +1, 1.0 shall be returned."

        Looking at the generated assembly, there's only one call to
        the pow function (and when commenting out the "indirect" line and
        compiling/linking without libm I get no error), so it looks like
        GCC decides to optimize away this call, even though I specified
        no optimization level.

        By the way, the spec also says:
        "If x is -1, and y is ±Inf, 1.0 shall be returned."
        This goes wrong as well:

        directly: pow(-1.000000, -inf) = nan
        indirectly: pow(-1.000000, -inf) = nan

        Looking at the generated assembly, this contains two calls
        to pow().  I don't know why GCC decides not to optimize away
        pow(-1.0, -1.0/0.0) but it does for pow(1.0, -1.0/0.0)....

        No patch since this code is a little above my head, but
        it looks like the problem occurs in src/lib/libm/src/e_pow.c:145

        That line reads
                 return  y - y; /* inf**+-1 is NaN */
        but just above, it says
                if (iy==0x7ff00000) {   /* y is +-inf */
        this means that the comment on line 145 is wrong, and I think the
        returned value is wrong as well because this seems to be the
        branch that's taken for pow(1, -inf).


Home | Main Index | Thread Index | Old Index