Something went wrong. Try again.
Cider Isn't Darwin Emulation, Really
Something went wrong. Try again.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334/* * Copyright (c) 2002 Apple Computer, Inc. All rights reserved. * * @APPLE_LICENSE_HEADER_START@ * * The contents of this file constitute Original Code as defined in and * are subject to the Apple Public Source License Version 1.1 (the * "License"). You may not use this file except in compliance with the * License. Please obtain a copy of the License at * http://www.apple.com/publicsource and read it before using this file. * * This 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 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 x�19. ** ********************************************************************************** ** 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. ** ********************************************************************************/
#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 const 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 const 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 const 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 const double LogSqrt2pi = 0.9189385332046727417803297e+0;static const double pi = 3.1415926535897932384626434e+0;static const double xbig = 171.624e+0;static const double MinimumX = 2.23e-308;static const double eps = 2.22e-16;static const hexdouble Huge = HEXDOUBLE(0x7FF00000, 0x00000000);static const hexdouble MinusHuge = HEXDOUBLE(0xFFF00000, 0x00000000);
#define GAMMA_NAN "42"#define SET_INVALID 0x01000000
#pragma STDC 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; if ( y < MinimumX ) { OldEnvironment.i.lo |= FE_OVERFLOW; FESETENVD( OldEnvironment.d ); return MinusHuge.d; } 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; }