Logo Search packages:      
Sourcecode: mingw-w64 version File versions  Download package

tgammaf.c
/**
 * This file has no copyright assigned and is placed in the Public Domain.
 * This file is part of the w64 mingw-runtime package.
 * No warranty is given; refer to the file DISCLAIMER.PD within this package.
 */
#include "cephes_mconf.h"

/* define MAXGAM 34.84425627277176174 */

/* Stirling's formula for the gamma function
 * gamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) ( 1 + 1/x P(1/x) )
 * .028 < 1/x < .1
 * relative error < 1.9e-11
 */
static const float STIR[] = {
  -2.705194986674176E-003,
   3.473255786154910E-003,
   8.333331788340907E-002,
};
static const float MAXSTIR = 26.77;
static const float SQTPIF = 2.50662827463100050242; /* sqrt( 2 pi ) */

static float stirf(float);

/* Gamma function computed by Stirling's formula,
 * sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
 * The polynomial STIR is valid for 33 <= x <= 172.
 */
static float stirf( float x )
{
      float  y, w, v;

      w = 1.0/x;
      w = 1.0 + w * polevlf(w, STIR, 2);
      y = expf(-x);
      if (x > MAXSTIR)
      { /* Avoid overflow in pow() */
            v = powf(x, 0.5 * x - 0.25);
            y *= v;
            y *= v;
      }
      else
      {
            y = powf(x, x - 0.5) * y;
      }
      y = SQTPIF * y * w;
      return (y);
}


/* gamma(x+2), 0 < x < 1 */
static const float P[] = {
  1.536830450601906E-003,
  5.397581592950993E-003,
  4.130370201859976E-003,
  7.232307985516519E-002,
  8.203960091619193E-002,
  4.117857447645796E-001,
  4.227867745131584E-001,
  9.999999822945073E-001,
};

float __tgammaf_r( float x, int* sgngamf);

float __tgammaf_r( float x, int* sgngamf)
{
      float p, q, z, nz;
      int i, direction, negative;

#ifdef NANS
      if (isnan(x))
            return (x);
#endif
#ifdef INFINITIES
#ifdef NANS
      if (x == INFINITYF)
            return (x);
      if (x == -INFINITYF)
            return (NANF);
#else
      if (!isfinite(x))
            return (x);
#endif
#endif

      *sgngamf = 1;
      negative = 0;
      nz = 0.0;
      if (x < 0.0)
      {
            negative = 1;
            q = -x;
            p = floorf(q);
            if (p == q)
            {
gsing:
                  _SET_ERRNO(EDOM);
                  mtherr("tgammaf", SING);
#ifdef INFINITIES
                  return (INFINITYF);
#else
                  return (MAXNUMF);
#endif
            }
            i = p;
            if ((i & 1) == 0)
                  *sgngamf = -1;
            nz = q - p;
            if (nz > 0.5)
            {
                  p += 1.0;
                  nz = q - p;
            }
            nz = q * sinf(PIF * nz);
            if (nz == 0.0)
            {
                  _SET_ERRNO(ERANGE);
                  mtherr("tgamma", OVERFLOW);
#ifdef INFINITIES
                  return(*sgngamf * INFINITYF);
#else
                  return(*sgngamf * MAXNUMF);
#endif
            }
            if (nz < 0)
                  nz = -nz;
            x = q;
      }
      if (x >= 10.0)
      {
            z = stirf(x);
      }
      if (x < 2.0)
            direction = 1;
      else
            direction = 0;
      z = 1.0;
      while (x >= 3.0)
      {
            x -= 1.0;
            z *= x;
      }
      /*
      while (x < 0.0)
      {
            if (x > -1.E-4)
                  goto Small;
            z *=x;
            x += 1.0;
      }
      */
      while (x < 2.0)
      {
            if (x < 1.e-4)
                  goto Small;
            z *=x;
            x += 1.0;
      }

      if (direction)
            z = 1.0/z;

      if (x == 2.0)
            return (z);

      x -= 2.0;
      p = z * polevlf(x, P, 7);

gdone:
      if (negative)
      {
            p = *sgngamf * PIF/(nz * p );
      }
      return (p);

Small:
      if (x == 0.0)
      {
            goto gsing;
      }
      else
      {
            p = z / ((1.0 + 0.5772156649015329 * x) * x);
            goto gdone;
      }
}

/* This is the C99 version */
float tgammaf(float x)
{
      int local_sgngamf = 0;
      return (__tgammaf_r(x, &local_sgngamf));
}


Generated by  Doxygen 1.6.0   Back to index