View Bug Activity | Format For Printing
The following program #include <math.h> #include <stdio.h> #define N_MAX 34 #undef M_PI #define M_PI 3.14159265358979323846 int main() { float c[N_MAX+1]; float xs, ts; float diff; int i; c[0] = 1.; for (i=1; i<=N_MAX; i++) c[i] = (2*i-1)*c[i-1]*0.5; for (i=1; i<=N_MAX; i++) { xs = i + 0.5f; /* argument to the gamma function */ ts = c[i] * sqrtf(M_PI); diff = fabsf(tgammaf(xs)-ts)/ts; if (diff > 1e-6) printf("Arg: %g, Exact: %12.8g, tgamma: %12.8g, diff: %12.8g\n", xs, ts, tgammaf(xs), diff); } return 0; } returns, on my Debian system, $ /usr/bin/gcc --std=c99 gamma.c -lm && ./a.out Arg: 18.5, Exact: 1.4986121e+15, tgamma: 1.4986145e+15, diff: 1.6121044e-06 Arg: 19.5, Exact: 2.7724323e+16, tgamma: 2.772436e+16, diff: 1.3167939e-06 Arg: 23.5, Exact: 5.3613034e+21, tgamma: 5.3612978e+21, diff: 1.0500245e-06 Arg: 29.5, Exact: 1.6348124e+30, tgamma: 1.6348178e+30, diff: 3.3277006e-06 Arg: 30.5, Exact: 4.8226972e+31, tgamma: 4.8226831e+31, diff: 2.9078208e-06 Arg: 31.5, Exact: 1.4709225e+33, tgamma: 1.4709265e+33, diff: 2.7352257e-06 Arg: 32.5, Exact: 4.633406e+34, tgamma: 4.6333912e+34, diff: 3.2061253e-06 Arg: 33.5, Exact: 1.5058569e+36, tgamma: 1.5058626e+36, diff: 3.7881605e-06 Arg: 34.5, Exact: 5.0446211e+37, tgamma: 5.0446287e+37, diff: 1.5077254e-06 We ran across this when a user reported inaccuracies in gfortran's gamma function, which in turn calls tgammaf() from the system library. The libc6 version of Debian is 2.6.1-1+b1; I am fairly certain that the Debian maintainers left the gamma function unchanged :-)
Two possible implementations at http://gcc.gnu.org/ml/fortran/2007-10/msg00197.html http://gcc.gnu.org/ml/fortran/2007-10/msg00201.html
There now is a C implementaton of the [tl]gamma* functions in libgfortran, contributed by FX Coudert: http://gcc.gnu.org/viewcvs/trunk/libgfortran/intrinsics/c99_functions.c?r1=128648&r2=130245 This might be OK for inclusion in glibc.