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

[Peter.Bex%xs4all.nl@localhost]

>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
>Organization:
N/A
>Environment:
        
        
System: NetBSD frohike.homeunix.org 5.0.1_PATCH NetBSD 5.0.1_PATCH (FROHIKE) 
#0: Tue Aug 11 22:31:56 CEST 2009 
sjamaan%frohike.homeunix.org@localhost:/usr/obj/sys/arch/amd64/compile/FROHIKE 
amd64
Architecture: x86_64
Machine: amd64
>Description:
        The pow() function does not behave like the specification at
        http://pubs.opengroup.org/onlinepubs/9699919799/functions/pow.html
        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.

>How-To-Repeat:
        $ 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
        wrong:

        "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)....

>Fix:
        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).

>Unformatted: