Re: Problems with GAMMA functions?

[Angelo Graziosi <Angelo dot Graziosi at roma1 dot infn dot it>]

On Sat, 13 Oct 2007, Thomas Koenig wrote:

>
> I have submitted this inaccuracy as a bug to the glibc people
> as http://sources.redhat.com/bugzilla/show_bug.cgi?id=5159 .
> Gcc, g++ and gfortran on glibc-based systems will all
> benefit if a solution is found there.
>
> If this is important for you, I would suggest that you contribute
> on this bug report.  In the meantime, you are probably better
> off using the external gamma function from CERNLIB.
>

Thomas,

I am not very expert of these things, but if I have seen correctly, the
core of GLIBC (and NEWLIB?) implementation of 'tgammaf' is the function
'__ieee754_gamma_r ' in e_gamma_r.c file. It computes the Gamma function
as exp(ln(...)): now it is clear why 'tgammaf' is not very accurate!

Computing something as exp(ln) is simple but not efficient, and some
comment in e_gamma_r.c confirms this.

The CERNLIB implementation is totally different and based on
"Approximation by truncated Chebyshev series and functional relations"
(http://wwwasdoc.web.cern.ch/wwwasdoc/cernlib.html,
http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/c303/top.html).

In C it can be summaryzed in this way:

==================================================================
double dgammf(double x)
{
   static const double PI = 3.14159265358979324;

   static const double C[] =
   {
      3.65738772508338244E0,
      1.95754345666126827E0,
      0.33829711382616039E0,
      0.04208951276557549E0,
      0.00428765048212909E0,
      0.00036521216929462E0,
      0.00002740064222642E0,
      0.00000181240233365E0,
      0.00000010965775866E0,
      0.00000000598718405E0,
      0.00000000030769081E0,
      0.00000000001431793E0,
      0.00000000000065109E0,
      0.00000000000002596E0,
      0.00000000000000111E0,
      0.00000000000000004E0
   };

   static int i;
   static double u,f,h,alfa,b0,b1,b2;

   u = x;

   if (u <= 0)
   {
      /* For non positive integer there are poles, return 0 is a choice */
      if (u == (int)x)
      {

	 h = 0;

	 goto LABEL_9;
      }
      else
      {
	 u = 1-x;
      }
   }

   f = 1;

   if (u < 3)
   {
      i = (int)(4-u);
      while (i--)
      {
	 f /= u;
	 u += 1;
      }
   }
   else
   {
      i = (int)(u-3);
      while (i--)
      {
	 u -= 1;
	 f *= u;
      }
   }

   h = u+u-7;
   alfa = h+h;
   b1 = 0;
   b2 = 0;

   for (i = 15; i >= 0; i--)
   {
      b0 = C[i]+alfa*b1-b2;
      b2 = b1;
      b1 = b0;
   }

   h = f*(b0-h*b2);

   if (x < 0) h = PI/(sin(PI*x)*h);

   LABEL_9:
      return h;
}

------------------------------------------------------------------

#define SROUND(D) ((D)+((D)-((float)(D))))

float tgammaf(float x)
{
   double dgammf(double);

   return (SROUND(dgammf((double)x)));
}

==================================================================

Using this implementation, your test case in
http://sources.redhat.com/bugzilla/show_bug.cgi?id=5159 does not print any
differences!


   Angelo.