lib/41931: log(-INFINITY) or log(-1.0) gives uexpected results

[he%NetBSD.org@localhost]

>Number:         41931
>Category:       lib
>Synopsis:       log(-INFINITY) or log(-1.0) gives uexpected results
>Confidential:   no
>Severity:       non-critical
>Priority:       medium
>Responsible:    lib-bug-people
>State:          open
>Class:          sw-bug
>Submitter-Id:   net
>Arrival-Date:   Mon Aug 24 13:40:00 +0000 2009
>Originator:     Havard Eidnes
>Release:        NetBSD 5.0_RC3
>Organization:
        I Try...
>Environment:
System: NetBSD smistad.uninett.no 5.0_RC3 NetBSD 5.0_RC3 (MAANEN) #4: Mon Mar 
30 02:52:40 CEST 2009 
he%maanen.uninett.no@localhost:/usr/obj/sys/arch/i386/compile/MAANEN i386
Architecture: i386
Machine: i386
>Description:
        It appears that ANSI C specifies that log(-INFINITY) or
        log(negative value) is supposed to return a NaN.
        Instead, -INFINITY is returned.

>How-To-Repeat:
        
------------------------------
#include <math.h>
#include <stdio.h>
#include <errno.h>

/* Ref. a.o. http://www.ansic.net/showlog.php?id=31&res=log%20function */

int
main(int argc, char **argv)
{
        double a, b;

        a = INFINITY;
        b = log(a);
        printf("log(Inf)  = %5f errno=%2d (expected: inf / %d)\n",
                b, errno, 0);

        a = -INFINITY;
        b = log(a);
        printf("log(-Inf) = %5f errno=%2d (expected: NaN / %d)\n",
                b, errno, EDOM);

        a = -1;
        b = log(a);
        printf("log(-1)   = %5f errno=%2d (expected: NaN / %d)\n",
                b, errno, EDOM);

        return 0;
}
------------------------------

        When run on my host, this program produces:

log(Inf)  =   inf errno= 0 (expected: inf / 0)
log(-Inf) =  -inf errno=33 (expected: NaN / 33)
log(-1)   =  -inf errno=33 (expected: NaN / 33)

        The same is the result when run on a sparc64/5.0 host.

        I've looked at lib/libm/src/w_log.c, I think the __ieee754_log()
        function is chosen, and I've looked at
        lib/libm/src/e_log.c which implements that function, and
        replicated the "hx/lx" extraction in my own code, and hx
        definitely ends up as negative when -INFINITY is passed,
        and a NaN results from the division by zero.

        The comments in e_log.c even say:

 * Special cases:
 *      log(x) is NaN with signal if x < 0 (including -INF) ;
 *      log(+INF) is +INF; log(0) is -INF with signal;
 *      log(NaN) is that NaN with no signal.

        which conforms with the expected behaviour and with
        section F.9.3.7 in ISO/IEC 9899 committee draft at
        http://www.open-std.org/Jtc1/sc22/wg14/www/docs/n1124.pdf

        My small test program:

------------------------------
#include <sys/types.h>

#include <math.h>
#include <stdio.h>

#if (BYTE_ORDER == BIG_ENDIAN) || (defined(__arm__) && !defined(__VFP_FP__))

typedef union
{
  double value;
  struct
  {
    u_int32_t msw;
    u_int32_t lsw;
  } parts;
} ieee_double_shape_type;

#endif

#if (BYTE_ORDER == LITTLE_ENDIAN) && \
    !(defined(__arm__) && !defined(__VFP_FP__))

typedef union
{
  double value;
  struct
  {
    u_int32_t lsw;
    u_int32_t msw;
  } parts;
} ieee_double_shape_type;

#endif

/* Get two 32 bit ints from a double.  */

#define EXTRACT_WORDS(ix0,ix1,d)                                \
do {                                                            \
  ieee_double_shape_type ew_u;                                  \
  ew_u.value = (d);                                             \
  (ix0) = ew_u.parts.msw;                                       \
  (ix1) = ew_u.parts.lsw;                                       \
} while (/*CONSTCOND*/0)

static const double zero = 0.0;


int
main(int argc, char **argv)
{
        double a, b;

        int32_t hx;
        uint32_t lx;

        a = -INFINITY;

        printf("Decode of -INFINITY:\n");
        EXTRACT_WORDS(hx, lx, a);

        if (hx < 0) {
                printf("hx < 0 (as expected)\n");
        } else {
                printf("hx >= 0 (NOT expected)\n");
        }

        a = -INFINITY;
        b = (a-a)/zero;
        printf("(-Inf - -Inf)/0.0 = %f (expected NaN)\n", b);

        return 0;
}
------------------------------

        which when run on my host produces

Decode of -INFINITY:
hx < 0 (as expected)
(-Inf - -Inf)/0.0 = nan (expected NaN)

        (The root trigger of this problem was an attempted upgrade
        of pkgsrc/lang/parrot to version 1.5.0, and a self-test tripped
        over this unexpected result.)

>Fix:
        Sorry, I don't know.