[Bug libfortran/68652] gamma function hangs on some arguments, returns NaN on other ones

["arbautjc at gmail dot com" <gcc-bugzilla at gcc dot gnu dot org>]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=68652

--- Comment #5 from Jean-Claude Arbaut <arbautjc at gmail dot com> ---
I posted the bug report here, I hope it's the right place now:
http://sourceforge.net/p/mingw-w64/bugs/517/

Nothing yet.

However, I had a look at the assembly output, and it's interesting: the Fortran
gamma function translates to a call to tgamma.

tgamma is a C99 function, which can be found in several places: Apple libm, gcc
source, glibc, netbsd... and the mingw-w64 source.

Actually, the assembly output also tells me it's a function in tgamma.c that
merely calls __tgamma_r. This is in the mingw-w64 source
(mingw-w64-v4.0.4\mingw-w64-crt\math\tgamma.c and tgammaf.c).

I tried a C program, and it shows the same bug, with some variations, with
mingw-w64 4.9.3 and 5.2.0 (the recent rev1), and also with Rtools (mingw-w64
4.6.3), the toolchain for the statistical package R. Similar program for
tgammaf:

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

int main() {
    double x, y;
    for(;;) {
        scanf("%lf", &x);
        y = tgamma(x);
        printf("%.4e %.4e\n", x, y);
    }
    return 0;
}

All in all:

* 32 bits GCC

tgamma returns NaN for integer x >= 710 (the program prints -1.#INDe+000), but
it's even worse than that: for 710.5 for instance, I get zero (also for 2000.5,
2100.5, and probably many other values of x).

tgammaf does *not* fail for large values, and it returns Inf. However, it gets
slower and slower for large values.

* 64 bits gcc

same problem with tgamma, though I don't get zeros.

tgammaf hangs for integer x >= 1677720 (either I was off by 1 in my preceding
post, or it's a slightly different version of mingw-w64, I'm not sure since it
wasn't on the same machine).

Just a note about the tgamma implementation of mingw-w64: it seems to use some
code from Cephes (a library by Stephen Moshier, see www.moshier.net), but
apparently only functions to evaluate polynomials.

I'll have to investigate further what really goes wrong in tgamma.c and
tgammaf.c. There is also a tgammal.c that is probably buggy as well.