e_lgammal_r.c

Glibc 2.3.2源代码(解压后有100多M)

static const long double RNr9[NRNr9 + 1] =
{
  4.441379198241760069548832023257571176884E5L,
  1.273072988367176540909122090089580368732E6L,
  9.732422305818501557502584486510048387724E5L,
  -5.040539994443998275271644292272870348684E5L,
  -1.208719055525609446357448132109723786736E6L,
  -7.434275365370936547146540554419058907156E5L,
  -2.075642969983377738209203358199008185741E5L,
  -2.565534860781128618589288075109372218042E4L,
  -1.032901669542994124131223797515913955938E3L,
};
#define NRDr9 8
static const long double RDr9[NRDr9 + 1] =
{
  -7.694488331323118759486182246005193998007E5L,
  -3.301918855321234414232308938454112213751E6L,
  -5.856830900232338906742924836032279404702E6L,
  -5.540672519616151584486240871424021377540E6L,
  -3.006530901041386626148342989181721176919E6L,
  -9.350378280513062139466966374330795935163E5L,
  -1.566179100031063346901755685375732739511E5L,
  -1.205016539620260779274902967231510804992E4L,
  -2.724583156305709733221564484006088794284E2L
/* 1.0E0 */
};


/* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */

static long double
neval (long double x, const long double *p, int n)
{
  long double y;

  p += n;
  y = *p--;
  do
    {
      y = y * x + *p--;
    }
  while (--n > 0);
  return y;
}


/* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */

static long double
deval (long double x, const long double *p, int n)
{
  long double y;

  p += n;
  y = x + *p--;
  do
    {
      y = y * x + *p--;
    }
  while (--n > 0);
  return y;
}


#ifdef __STDC__
long double
__ieee754_lgammal_r (long double x, int *signgamp)
#else
long double
__ieee754_lgammal_r (x, signgamp)
     long double x;
     int *signgamp;
#endif
{
  long double p, q, w, z, nx;
  int i, nn;

  *signgamp = 1;

  if (! __finitel (x))
    return x * x;

  if (x < 0.0L)
    {
      q = -x;
      p = __floorl (q);
      if (p == q)
	return (one / (p - p));
      i = p;
      if ((i & 1) == 0)
	*signgamp = -1;
      else
	*signgamp = 1;
      z = q - p;
      if (z > 0.5L)
	{
	  p += 1.0L;
	  z = p - q;
	}
      z = q * __sinl (PIL * z);
      if (z == 0.0L)
	return (*signgamp * huge * huge);
      w = __ieee754_lgammal_r (q, &i);
      z = __logl (PIL / z) - w;
      return (z);
    }

  if (x < 13.5L)
    {
      p = 0.0L;
      nx = __floorl (x + 0.5L);
      nn = nx;
      switch (nn)
	{
	case 0:
	  /* log gamma (x + 1) = log(x) + log gamma(x) */
	  if (x <= 0.125)
	    {
	      p = x * neval (x, RN1, NRN1) / deval (x, RD1, NRD1);
	    }
	  else if (x <= 0.375)
	    {
	      z = x - 0.25L;
	      p = z * neval (z, RN1r25, NRN1r25) / deval (z, RD1r25, NRD1r25);
	      p += lgam1r25b;
	      p += lgam1r25a;
	    }
	  else if (x <= 0.625)
	    {
	      z = x + (1.0L - x0a);
	      z = z - x0b;
	      p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
	      p = p * z * z;
	      p = p + y0b;
	      p = p + y0a;
	    }
	  else if (x <= 0.875)
	    {
	      z = x - 0.75L;
	      p = z * neval (z, RN1r75, NRN1r75) / deval (z, RD1r75, NRD1r75);
	      p += lgam1r75b;
	      p += lgam1r75a;
	    }
	  else
	    {
	      z = x - 1.0L;
	      p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2);
	    }
	  p = p - __logl (x);
	  break;

	case 1:
	  if (x < 0.875L)
	    {
	      if (x <= 0.625)
		{
		  z = x + (1.0L - x0a);
		  z = z - x0b;
		  p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
		  p = p * z * z;
		  p = p + y0b;
		  p = p + y0a;
		}
	      else if (x <= 0.875)
		{
		  z = x - 0.75L;
		  p = z * neval (z, RN1r75, NRN1r75)
		        / deval (z, RD1r75, NRD1r75);
		  p += lgam1r75b;
		  p += lgam1r75a;
		}
	      else
		{
		  z = x - 1.0L;
		  p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2);
		}
	      p = p - __logl (x);
	    }
	  else if (x < 1.0L)
	    {
	      z = x - 1.0L;
	      p = z * neval (z, RNr9, NRNr9) / deval (z, RDr9, NRDr9);
	    }
	  else if (x == 1.0L)
	    p = 0.0L;
	  else if (x <= 1.125L)
	    {
	      z = x - 1.0L;
	      p = z * neval (z, RN1, NRN1) / deval (z, RD1, NRD1);
	    }
	  else if (x <= 1.375)
	    {
	      z = x - 1.25L;
	      p = z * neval (z, RN1r25, NRN1r25) / deval (z, RD1r25, NRD1r25);
	      p += lgam1r25b;
	      p += lgam1r25a;
	    }
	  else
	    {
	      /* 1.375 <= x+x0 <= 1.625 */
	      z = x - x0a;
	      z = z - x0b;
	      p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
	      p = p * z * z;
	      p = p + y0b;
	      p = p + y0a;
	    }
	  break;

	case 2:
	  if (x < 1.625L)
	    {
	      z = x - x0a;
	      z = z - x0b;
	      p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
	      p = p * z * z;
	      p = p + y0b;
	      p = p + y0a;
	    }
	  else if (x < 1.875L)
	    {
	      z = x - 1.75L;
	      p = z * neval (z, RN1r75, NRN1r75) / deval (z, RD1r75, NRD1r75);
	      p += lgam1r75b;
	      p += lgam1r75a;
	    }
	  else if (x == 2.0L)
	    p = 0.0L;
	  else if (x < 2.375L)
	    {
	      z = x - 2.0L;
	      p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2);
	    }
	  else
	    {
	      z = x - 2.5L;
	      p = z * neval (z, RN2r5, NRN2r5) / deval (z, RD2r5, NRD2r5);
	      p += lgam2r5b;
	      p += lgam2r5a;
	    }
	  break;

	case 3:
	  if (x < 2.75)
	    {
	      z = x - 2.5L;
	      p = z * neval (z, RN2r5, NRN2r5) / deval (z, RD2r5, NRD2r5);
	      p += lgam2r5b;
	      p += lgam2r5a;
	    }
	  else
	    {
	      z = x - 3.0L;
	      p = z * neval (z, RN3, NRN3) / deval (z, RD3, NRD3);
	      p += lgam3b;
	      p += lgam3a;
	    }
	  break;

	case 4:
	  z = x - 4.0L;
	  p = z * neval (z, RN4, NRN4) / deval (z, RD4, NRD4);
	  p += lgam4b;
	  p += lgam4a;
	  break;

	case 5:
	  z = x - 5.0L;
	  p = z * neval (z, RN5, NRN5) / deval (z, RD5, NRD5);
	  p += lgam5b;
	  p += lgam5a;
	  break;

	case 6:
	  z = x - 6.0L;
	  p = z * neval (z, RN6, NRN6) / deval (z, RD6, NRD6);
	  p += lgam6b;
	  p += lgam6a;
	  break;

	case 7:
	  z = x - 7.0L;
	  p = z * neval (z, RN7, NRN7) / deval (z, RD7, NRD7);
	  p += lgam7b;
	  p += lgam7a;
	  break;

	case 8:
	  z = x - 8.0L;
	  p = z * neval (z, RN8, NRN8) / deval (z, RD8, NRD8);
	  p += lgam8b;
	  p += lgam8a;
	  break;

	case 9:
	  z = x - 9.0L;
	  p = z * neval (z, RN9, NRN9) / deval (z, RD9, NRD9);
	  p += lgam9b;
	  p += lgam9a;
	  break;

	case 10:
	  z = x - 10.0L;
	  p = z * neval (z, RN10, NRN10) / deval (z, RD10, NRD10);
	  p += lgam10b;
	  p += lgam10a;
	  break;

	case 11:
	  z = x - 11.0L;
	  p = z * neval (z, RN11, NRN11) / deval (z, RD11, NRD11);
	  p += lgam11b;
	  p += lgam11a;
	  break;

	case 12:
	  z = x - 12.0L;
	  p = z * neval (z, RN12, NRN12) / deval (z, RD12, NRD12);
	  p += lgam12b;
	  p += lgam12a;
	  break;

	case 13:
	  z = x - 13.0L;
	  p = z * neval (z, RN13, NRN13) / deval (z, RD13, NRD13);
	  p += lgam13b;
	  p += lgam13a;
	  break;
	}
      return p;
    }

  if (x > MAXLGM)
    return (*signgamp * huge * huge);

  q = ls2pi - x;
  q = (x - 0.5L) * __logl (x) + q;
  if (x > 1.0e18L)
    return (q);

  p = 1.0L / (x * x);
  q += neval (p, RASY, NRASY) / x;
  return (q);
}