gamma.c   [plain text]


/*
 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
 *
 * @APPLE_LICENSE_HEADER_START@
 * 
 * Copyright (c) 1999-2003 Apple Computer, Inc.  All Rights Reserved.
 * 
 * This file contains Original Code and/or Modifications of Original Code
 * as defined in and that are subject to the Apple Public Source License
 * Version 2.0 (the 'License'). You may not use this file except in
 * compliance with the License. Please obtain a copy of the License at
 * http://www.opensource.apple.com/apsl/ and read it before using this
 * file.
 * 
 * The Original Code and all software distributed under the License are
 * distributed on an 'AS IS' basis, WITHOUT WARRANTY OF ANY KIND, EITHER
 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE, QUIET ENJOYMENT OR NON-INFRINGEMENT.
 * Please see the License for the specific language governing rights and
 * limitations under the License.
 * 
 * @APPLE_LICENSE_HEADER_END@
 */
/*******************************************************************************
*                                                                              *
*      File gamma.c,                                                           *
*      Functions gamma(x),                                                     *
*      Implementation of the gamma function for the PowerPC.                   *
*                                                                              *
*      Copyright c 1993 Apple Computer, Inc.  All rights reserved.             *
*                                                                              *
*      Written by Ali Sazegari, started on January 1993,                       *
*                                                                              *
*      based on FORTRAN routines in SpecFun by J. W. Cody and L. Stoltz of     *
*      Applied Mathematics Division of Argonne National Laboratory,            *
*      Argonne, IL 60439.                                                      *
*                                                                              *
*      W A R N I N G:  This routine expects a 64 bit double model.             *
*                                                                              *
*      January  29  1993: First C implementation for PowerPC.                  *
*      April    26  1993: fixed a few bugs in the interval [xbig,inf].         *
*      July     14  1993: added #pragma fenv_access. This function is now      *
*                         using the the string oriented nan.  replaced       *
*                         feholdexcept by _feprocentry to guard rounding.      *
*      July     29  1993: corrected the string nan.                            *
*      October  07  1993: removed the raising of the overflow flag for arg= . *
*      September19  1994: changed all environemental funtions to __setflm,     *
*                         changed HUGE_VAL to Huge.d for perfrormance.       *
*      January  11  1995: a humilating error corrected.  in the interval       *
*                         [12,178], the nonexistant array element c[7] is      *
*                         addressed.  it should be c[6].                       *
*                         a short error analysis reveals that in double        *
*                         precision referencing c[7] instead of c[6] has no    *
*                         impact on the accuracy of the result, provided that  *
*                         the compiler assigns 0.0 to c[7], which the ibm xlc  *
*                         does.  this explains why the double precision        *
*                         gamma function passed accuracy tests.  the relative  *
*                         error in extended is at most 5.56E-17 and occurs     *
*                         for x=12.0.  the calculation is no longer affected   *
*                         for arguments x19.                                  *
*                                                                              *
********************************************************************************
*                                                                              *
*                              G  A  M  M  A  (  X  )                          *
*                                                                              *
*      The gamma function is an increasing function over [2,inf].  For large   *
*      enough x, it behaves like [ e^( x ln ( x/ e ) ] / sqrt(x/pi).  It may   *
*      be more appropriate to work with the logorithm of Gamma.                *
*                                                                              *
*      This routine calculates the gamma function for a real argument x.       *
*      Computation is based on an algorithm outlined in reference 1 below.     *
*      The program uses rational functions that approximate the gamma          *
*      function to at least 20 significant decimal digits.  Coefficients       *
*      for the approximation over the interval (1,2) are unpublished.          *
*      Those for the approximation for x >= 12 are from reference 2 below.     *
*                                                                              *
********************************************************************************
*                                                                              *
*      Explanation of machine-dependent constants:                             *
*                                                                              *
*      xbig     - the largest argument for which gamma(x) is representable     *
*                 in the machine, i.e., the solution to the equation           *
*                 gamma ( xbig ) = 2 ** MaxExp, where MaxExp is the maximum    *
*                 power of 2 before infinity;                                  *
*      xinf     - the largest machine representable floating-point number      *
*                 before infinity, approximately 2 ** MaxExp;                  *
*      eps      - the smallest positive floating-point number such that        *
*                 1.0 + eps > 1.0;                                             *
*      MinimumX - the smallest positive floating-point number such that        *
*                 1/MinimumX is machine representable.                         *
*                                                                              *
*      Approximate values for the macintosh and the powerpc are:               *
*                                                                              *
*                                base       MaxExp        xbig                 *
*                                                                              *
*      Macintosh     extended     2         16382        1755.36               *
*      PowerPC        double      2          1023        171.624               *
*                                                                              *
*                                 xinf         eps        MinimumX             *
*                                                                              *
*      Macintosh     extended   1.19x+4932   5.42x-20    8.40x-4933            *
*      PowerPC        double    1.79d+308    2.22d-16    2.23d-308             *
*                                                                              *
********************************************************************************
*                                                                              *
*      The program returns a quiet NaN for singularities and infinity when     *
*      overflow occurs.  The computation is believed to be free of undeserved  *
*      underflow and overflow.                                                 *
*                                                                              *
*      References:                                                             *
*                                                                              *
*      [1] "An Overview of Software Development for Special Functions",        *
*          W. J. Cody, Lecture Notes in Mathematics, 506, Numerical Analysis   *
*          Dundee, 1975, G. A. Watson (ed.), Springer Verlag, Berlin, 1976.    *
*                                                                              *
*      [2] Computer Approximations, Hart, Et. Al., Wiley and sons, New York    *
*          1968.                                                               *
*                                                                              *
*******************************************************************************/
#ifdef      __APPLE_CC__
#if         __APPLE_CC__ > 930

#include      "math.h"
#include      "fenv.h"
#include      "fp_private.h"
#include      "fenv_private.h"

/*******************************************************************************
*            Functions needed for the computation.                             *
*******************************************************************************/

/*     the following fp.h functions are used:                                 */
/*     exp, log, sin, __fpclassifyd and nan;                                  */
/*     the following environmental functions are used:                        */
/*     __setflm.                                                              */


/*******************************************************************************
*     Coefficients for P in gamma approximation over [1,2] in decreasing order.*
*******************************************************************************/

static double p[8] = { -1.71618513886549492533811e+0,
                        2.47656508055759199108314e+1,
                       -3.79804256470945635097577e+2,
                        6.29331155312818442661052e+2,
                        8.66966202790413211295064e+2,
                       -3.14512729688483675254357e+4,
                       -3.61444134186911729807069e+4,
                        6.64561438202405440627855e+4 };

/*******************************************************************************
*     Coefficients for Q in gamma approximation over [1,2] in decreasing order.*
*******************************************************************************/

static double q[8] = { -3.08402300119738975254353e+1,
                        3.15350626979604161529144e+2,
                       -1.01515636749021914166146e+3,
                       -3.10777167157231109440444e+3,
                        2.25381184209801510330112e+4,
                        4.75584627752788110767815e+3,
                       -1.34659959864969306392456e+5,
                       -1.15132259675553483497211e+5 };
                         
/*******************************************************************************
*     Coefficients for minimax approximation over [12, INF].                   *
*******************************************************************************/

static double c[7] = { -1.910444077728e-03,
                        8.4171387781295e-04,
                       -5.952379913043012e-04,
                        7.93650793500350248e-04,
                       -2.777777777777681622553e-03,
                        8.333333333333333331554247e-02,
                        5.7083835261e-03 };

static double LogSqrt2pi = 0.9189385332046727417803297e+0;
static double pi         = 3.1415926535897932384626434e+0;
static double xbig       = 171.624e+0;
static double MinimumX   = 2.23e-308;
static double eps        = 2.22e-16;
static hexdouble Huge    = HEXDOUBLE(0x7FF00000, 0x00000000);

#define      GAMMA_NAN      "42"
#define      SET_INVALID    0x01000000

#pragma fenv_access on

double tgamma ( double x )
      {
      register int n, parity, i;
      register double y, y1, result, fact, IsItAnInt, z, numerator, 
	                denominator, ysquared, sum; 
      hexdouble OldEnvironment;
    
      fegetenvd( OldEnvironment.d );               // save environment, set default
      fesetenvd( 0.0 );
	
/*******************************************************************************
*     The next switch will decipher what sort of argument we have. If argument *
*     is SNaN then a QNaN has to be returned and the invalid flag signaled.    * 
*******************************************************************************/

      switch ( __fpclassifyd ( x ) )
            {
            case FP_NAN:
                  x *= 2.0;                  /* quiets NaN */
                  fesetenvd( OldEnvironment.d ); //   restore caller's environment
                  return x;
                  
            case FP_ZERO:
                  OldEnvironment.i.lo |= FE_DIVBYZERO;
                  fesetenvd( OldEnvironment.d );
                  return copysign( Huge.d, x);

             case FP_INFINITE:
                  if ( x > 0.0 )
                        x = Huge.d;
                  else
                        {
                        x = nan ( GAMMA_NAN );
                        OldEnvironment.i.lo |= SET_INVALID;
                        }

                  fesetenvd( OldEnvironment.d );
                  return x;
                  
            default:                  /*      NORMALNUM and DENORMALNUM      */
                  break;
            }
      
      parity = 0;
      fact = 1.0;
      n = 0;
      y = x;

/*******************************************************************************
*     The argument is negative.                                                *
*******************************************************************************/

      if ( y <= 0.0 )
            {
            y = - x;
            y1 = trunc ( y );
            IsItAnInt = y - y1;
            if ( IsItAnInt != 0.0 )                   /* is it an integer?   */
                  {                                   /* is it odd or even?  */
                  if ( y1 != trunc ( y1 * 0.5 ) * 2.0 ) parity = 1;
                  fact = - pi / sin ( pi * IsItAnInt );
                  y += 1.0;
                  }
            else
                  {
                  OldEnvironment.i.lo |= SET_INVALID;
                  fesetenvd( OldEnvironment.d );
                  return nan ( GAMMA_NAN );
                  }
            }

/*******************************************************************************
*     The argument is positive.                                                *
*******************************************************************************/

     if ( y < eps )                         /* argument is less than epsilon. */
            {
            if ( y >= MinimumX )          /* x is in [MinimumX,eps].          */
                  result = 1.0 / y;
            else                          /* othewise, x is in [0,MinimumX).  */
                  {
                  OldEnvironment.i.lo |= FE_OVERFLOW;
                  fesetenvd( OldEnvironment.d );
                  return Huge.d;
                  }
            }
      else if ( y < 12.0 )                 /* argument x is eps < x < 12.0.  */
            {
            y1 = y;
            if ( y < 1.0 )                 /* x is in (eps, 1.0).            */
                  {
                  z = y;
                  y += 1.0;
                  }
            else                           /* x is in [1.0,12.0].            */
                  {
                  n = ( int ) y - 1;
                  y -= ( double ) n;
                  z = y - 1.0;
                  }
            numerator = 0.0;
            denominator = 1.0;
            for ( i = 0; i < 8; i++ )
                  {
                  numerator = ( numerator + p[i] ) * z;
                  denominator = denominator * z + q[i];
                  }
            result = numerator / denominator + 1.0;
            if ( y1 < y )
                  result /= y1;
            else if ( y1 > y )
                  {
                  for ( i = 0; i < n; i++ )
                        {
                        result *= y;
                        y += 1.0;
                        }
                  }
            }
      else
            {
            if ( x <= xbig )
                  {
                  ysquared = y * y;
                  sum = c[6];
                  for ( i = 0; i < 6; i++ )
                        sum = sum / ysquared + c[i];
                  sum = sum / y - y + LogSqrt2pi;
                  sum += ( y - 0.5 ) * log ( y );
                  result = exp ( sum );
                  }
            else
                  {
                  OldEnvironment.i.lo |= FE_OVERFLOW;
                  fesetenvd( OldEnvironment.d );         //   restore caller's environment
                  return Huge.d;
                  }
            }
      if ( parity ) result = - result;
      if ( fact != 1.0 ) result = fact / result;
      fesetenvd( OldEnvironment.d );         //   restore caller's environment
      return result;
      }
      
#ifdef notdef
float tgammaf ( float x )
{
    return (float)tgamma ( x );
}
#endif

#else       /* __APPLE_CC__ version */
#error Version gcc-932 or higher required.  Compilation terminated.
#endif      /* __APPLE_CC__ version */
#endif      /* __APPLE_CC__ */