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.
You didn't mention what architecture this was on. Anyway, confirmed with current sources on x86, and with more and greater errors on x86_64 (I suppose excess precision is helping a bit on x86). Presumably the same cause as other tgamma bugs - all five versions (flt-32, dbl-64, ldbl-96, ldbl-128, ldbl-128ibm) follow the same approach of exp (lgamma) which is problematic for results with large (positive or negative) exponent and needs better implementations for all five cases.
Fixed for 2.18 by: commit d8cd06db62d92f86cc8cc3c0d6a489bd207bb834 Author: Joseph Myers <joseph@codesourcery.com> Date: Wed May 8 11:58:18 2013 +0000 Improve tgamma accuracy (bugs 2546, 2560, 5159, 15426).