Bug 5159 - Inaccuracies in tgammaf
Summary: Inaccuracies in tgammaf
Status: RESOLVED FIXED
Alias: None
Product: glibc
Classification: Unclassified
Component: math (show other bugs)
Version: unspecified
: P2 normal
Target Milestone: ---
Assignee: Not yet assigned to anyone
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2007-10-10 18:48 UTC by Thomas Koenig
Modified: 2014-07-04 07:21 UTC (History)
1 user (show)

See Also:
Host:
Target:
Build:
Last reconfirmed:
fweimer: security-


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description Thomas Koenig 2007-10-10 18:48:12 UTC
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 :-)
Comment 1 Thomas Koenig 2007-10-14 16:56:10 UTC
Two possible implementations at

http://gcc.gnu.org/ml/fortran/2007-10/msg00197.html
http://gcc.gnu.org/ml/fortran/2007-10/msg00201.html


Comment 2 Thomas Koenig 2007-11-17 17:07:35 UTC
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.
Comment 3 Joseph Myers 2012-02-29 19:57:17 UTC
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.
Comment 4 Joseph Myers 2013-05-08 12:00:12 UTC
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).