#517 gamma function hangs on some arguments, returns NaN on other ones

[Jean-Claude Arbaut]

First posted at https://gcc.gnu.org/bugzilla/show_bug.cgi?id=68652. Probably better here.

This happens with gfortran in MinGW-W64 (precisely x86_64-4.9.3-release-win32-seh-rt_v4-rev1.7z at http://sourceforge.net/projects/mingw-w64/files/Toolchains%20targetting%20Win64/Personal%20Builds/mingw-builds/4.9.3/threads-win32/seh/). So it may be a problem with components external to gfortran. The OS is Windows 7 Pro, 64 bits.

program bug
implicit none
real(4) :: x, y
do
read , x
y = gamma(x)
print , y
end do
end program

With real(4), and an argument that is an integer value (in a real variable, of course):

x=1 to 35: y has a finite real value, no problem
x=36 to 16777218: y is Infinity, as expected
x=16777219 or above: program hangs.
Notice that 2^24=16777216, by the way.

With real(8)
x=1 to 171: finite real
x=172 to 709: Infinity
x=710 or above: NaN, which is not expected (it should be Infinity).

[Jean-Claude Arbaut]

The gamma function in gfortran is simply a call to the C99 function tgamma (or tgammaf in real(4) precision). A similar program in C shows the same bug. It seems that MinGW-w64 is using a gamma implementation from the CEPHES library by Stephen Moshier, hence the bug may lie here.
I have this bug with gcc 4.9.3 and 5.2.0 with the sourceforge binaries, and with Rtools, which uses MinGW-w64/gcc 4.6.3.

I have another concern: it's not clear that the CEPHES license is compatible with GPL. See for example here: https://lists.debian.org/debian-legal/2004/12/msg00295.html It was in 2004, but it does not look like this problem has been solved since. And in the latest source zip files on Moshier's site (2014-10-04), there is absolutely no mention that the code is licensed under a GPL compatible license, so I assume this is not the case.

[Jean-Claude Arbaut]

Same program with appropriate quotes so that the code is still readable.

program bug
    implicit none
    real(4) :: x, y
    do
        read *, x
        y = gamma(x)
        print *, y
    end do
end program

Notice the same bug happens with gcc.

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

int main() {
    float x, y;
    for(;;) {
        scanf("%f", &x);
        y = tgammaf(x);
        printf("%f\n", y);
    }
}

With the C version, the program hangs for 16777220. The "double" version also returns NaN.

This has been tested with the compilers cited above, and also with the dongsheng build of MinGW-W64 (4.9.4, 5.3.1, 6.0.0).

The first link in the original message does not work (because of a dot at the end of the link), here is a working link: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=68652.

[WinPORTS]

Using GCC 6.1.0, we get the following:

16777220
1.#INF00