#2306 powl() fails to return Inf in some cases

[Dan Collins]

The powl() shipped with MinGW is defective. 9**(9**9) should return Inf, as this value is orders of magnitude above LDBL_MAX. Instead, this returns LDBL_MAX. Example:

lddemo.c:

#include <tgmath.h>
#include <stdio.h>

int main() {
  long double a, b, c;
  a = 9;
  b = powl(a, a);
  c = powl(a, b);
  __mingw_printf("powl(9, powl(9, 9)) = %Le\n", c);
}

Console session:

W:\buildbot>set PATH=C:\MinGW\bin;%PATH%

W:\buildbot>gcc -v
Using built-in specs.
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=c:/mingw/bin/../libexec/gcc/mingw32/4.9.3/lto-wrapper.exe
Target: mingw32
Configured with: ../src/gcc-4.9.3/configure --build=x86_64-pc-linux-gnu --host=mingw32 --prefix=/mingw --disable-win32-registry --target=mingw32 --with-arch=i586 --enable-languages=c,c++,objc,obj-c++,fortran,ada --enable-static --enable-shared --enable-threads --with-dwarf2 --disable-sjlj-exceptions --enable-version-specific-runtime-libs --enable-libstdcxx-debug --with-tune=generic --enable-nls
Thread model: win32
gcc version 4.9.3 (GCC)

W:\buildbot>gcc lddemo.c -Wall
lddemo.c: In function 'main':
lddemo.c:10:1: warning: control reaches end of non-void function [-Wreturn-type]
 }
 ^

W:\buildbot>a.exe
powl(9, powl(9, 9)) = 1.189731e+4932

_mingw.h:

#define __MINGW32_VERSION           3022000L
#define __MINGW32_MAJOR_VERSION           3
#define __MINGW32_MINOR_VERSION          22
#define __MINGW32_PATCHLEVEL              0

Expected output:

W:\buildbot>set PATH=C:\mingw-w64\x86_64-4.9.3\mingw64\bin;%PATH%

W:\buildbot>gcc lddemo.c -Wall
lddemo.c: In function 'main':
lddemo.c:10:1: warning: control reaches end of non-void function [-Wreturn-type]
 }
 ^

W:\buildbot>a.exe
powl(9, powl(9, 9)) = inf

Regards,
Dan Collins

[Keith Marshall]

Thanks for bringing this to our attention. FWIW, the powl() implementation in question is a quick and dirty adaptation of Stephen Moshier's Cephes maths library implementation, (which must come close to being the most hideously ugly code I've ever seen); it looks like this bug has been inherited from the upstream source, where I see:

/* Test the power of 2 for overflow */
if( w > MEXP )
    {
/*  printf( "w = %.4Le ", w ); */
    mtherr( fname, OVERFLOW );
    return( MAXNUML );
    }

Here, MAXNUML is equivalent to ISO-C's LDBL_MAX, which is the (incorrect) return value you are seeing; the correct return value should be HUGE_VALL, (which corresponds to infinity).

[Keith Marshall]

FWIW, if I modify this fragment of the upstream source:

/* Test the power of 2 for overflow */
if( w > MEXP )
    {
/*  printf( "w = %.4Le ", w ); */
    mtherr( fname, OVERFLOW );
    return( MAXNUML );
    }

such that it becomes:

/* Test the power of 2 for overflow */
if( w > MEXP )
    {
/*  printf( "w = %.4Le ", w ); */
    mtherr( fname, OVERFLOW );
    return( INFINITYL );
    }

then your test case returns the expected inf. I may push this into mingwrt-3.22.1, as a Q&D resolution for this issue, but I do feel rather uncomfortable publishing such Q&D modifications to such hideously ugly code. Ideally, this requires a complete overhaul, and furthermore, I am inclined to wonder if the complexity of table lookups for base two logarithms, and their anti-logarithms, really offers any advantage over a direct computation, (perhaps even within the FPU directly, via inline assembly code), of the IEEE 754 80-bit extended precision analytical expression:

x ** y = 2 ** (log2(x) * y)

(with ** representing the exponentiation operator, as it does in FORTRAN).

[Keith Marshall]

This particular aspect of this issue is corrected by commit [9dfbd7]; the modification is included within mingwrt-3.22.1. However, I am not closing the issue at this juncture, because there are other related boundary cases, for which POSIX prescribes a particular outcome for any IEC-60559 conforming implementation, (which any 8087 derivative FPU should deliver), and with which this implementation does not comply; there may be others, but in particular:

• For any value of y (including NaN), if x is +1, 1.0 shall be returned; implementation incorrectly returns NaN.
• If x is -1, and y is ±Inf, 1.0 shall be returned; implementation incorrectly returns ±Inf.

[Keith Marshall]

As I suspected, the particular failure initially reported here is only one of many non-compliances, with the boundary condition return values prescribed by POSIX. To demonstrate the issues, I've compiled the attached set of unit tests, which may be integrated into the autotest framework I committed as [f8d4c2]; when hooked into testsuite.at.in immediately after headers.at, and run with:

$ make -C mingwrt check TESTSUITEFLAGS='-k powl'

it exhibits 21 failures, out of 55 unit tests run; (that's a 38% failure rate, which I don't think should be acceptable).

[Keith Marshall]

I've now extended the previously proposed testsuite module, to 60 unit tests, and committed it to mingwrt as [5b105e] ; even with the adjustment of [9dfbd7], this still exhibits the 21 POSIX compliance failures, as indicated in the attached log file.

[Keith Marshall]

For comparison, when subjected to the same 60 unit tests, Microsoft's pow() implementation exhibits 11 departures from POSIXly correct behaviour, (and FWIW, the equivalent wine implementation stretches this to 21).

[Keith Marshall]

I've now committed [ae8972], which completely replaces the old powf() and powl() implementations from Stephen Moshier's Cephes Math Library, with a pure x87 FPU implementation. This new implementation is:--
a) Significantly faster than the older implementations, and
b) Passes all 60 of the specified unit tests, in respect of both functions.

Note that the Microsoft pow() implementation, from MSVCRT.DLL has not been replaced, and still fails 11 of the 60 unit tests, (21 on wine). However, the replacement implementation does include pow(), not currently added to libmingwex.a, but could be included in some other supplementary extension, (or maybe libm.a); this potential replacement does pass all 60 unit tests.