overby.me
this repo has no description
1/*
2 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
3 *
4 * @APPLE_LICENSE_HEADER_START@
5 *
6 * The contents of this file constitute Original Code as defined in and
7 * are subject to the Apple Public Source License Version 1.1 (the
8 * "License"). You may not use this file except in compliance with the
9 * License. Please obtain a copy of the License at
10 * http://www.apple.com/publicsource and read it before using this file.
11 *
12 * This Original Code and all software distributed under the License are
13 * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
14 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
15 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the
17 * License for the specific language governing rights and limitations
18 * under the License.
19 *
20 * @APPLE_LICENSE_HEADER_END@
21 */
22/*******************************************************************************
23* *
24* File gamma.c, *
25* Functions gamma(x), *
26* Implementation of the gamma function for the PowerPC. *
27* *
28* Copyright c 1993 Apple Computer, Inc. All rights reserved. *
29* *
30* Written by Ali Sazegari, started on January 1993, *
31* *
32* based on FORTRAN routines in SpecFun by J. W. Cody and L. Stoltz of *
33* Applied Mathematics Division of Argonne National Laboratory, *
34* Argonne, IL 60439. *
35* *
36* W A R N I N G: This routine expects a 64 bit double model. *
37* *
38* January 29 1993: First C implementation for PowerPC. *
39* April 26 1993: fixed a few bugs in the interval [xbig,inf]. *
40* July 14 1993: added #pragma fenv_access. This function is now *
41* using the the string oriented �nan�. replaced *
42* feholdexcept by _feprocentry to guard rounding. *
43* July 29 1993: corrected the string nan. *
44* October 07 1993: removed the raising of the overflow flag for arg= �. *
45* September19 1994: changed all environemental funtions to __setflm, *
46* changed HUGE_VAL to Huge.d for perfrormance. *
47* January 11 1995: a humilating error corrected. in the interval *
48* [12,178], the nonexistant array element c[7] is *
49* addressed. it should be c[6]. *
50* a short error analysis reveals that in double *
51* precision referencing c[7] instead of c[6] has no *
52* impact on the accuracy of the result, provided that *
53* the compiler assigns 0.0 to c[7], which the ibm xlc *
54* does. this explains why the double precision *
55* gamma function passed accuracy tests. the relative *
56* error in extended is at most 5.56E-17 and occurs *
57* for x=12.0. the calculation is no longer affected *
58* for arguments x�19. *
59* *
60********************************************************************************
61* *
62* G A M M A ( X ) *
63* *
64* The gamma function is an increasing function over [2,inf]. For large *
65* enough x, it behaves like [ e^( x ln ( x/ e ) ] / sqrt(x/pi). It may *
66* be more appropriate to work with the logorithm of Gamma. *
67* *
68* This routine calculates the gamma function for a real argument x. *
69* Computation is based on an algorithm outlined in reference 1 below. *
70* The program uses rational functions that approximate the gamma *
71* function to at least 20 significant decimal digits. Coefficients *
72* for the approximation over the interval (1,2) are unpublished. *
73* Those for the approximation for x >= 12 are from reference 2 below. *
74* *
75********************************************************************************
76* *
77* Explanation of machine-dependent constants: *
78* *
79* xbig - the largest argument for which gamma(x) is representable *
80* in the machine, i.e., the solution to the equation *
81* gamma ( xbig ) = 2 ** MaxExp, where MaxExp is the maximum *
82* power of 2 before infinity; *
83* xinf - the largest machine representable floating-point number *
84* before infinity, approximately 2 ** MaxExp; *
85* eps - the smallest positive floating-point number such that *
86* 1.0 + eps > 1.0; *
87* MinimumX - the smallest positive floating-point number such that *
88* 1/MinimumX is machine representable. *
89* *
90* Approximate values for the macintosh and the powerpc are: *
91* *
92* base MaxExp xbig *
93* *
94* Macintosh extended 2 16382 1755.36 *
95* PowerPC double 2 1023 171.624 *
96* *
97* xinf eps MinimumX *
98* *
99* Macintosh extended 1.19x+4932 5.42x-20 8.40x-4933 *
100* PowerPC double 1.79d+308 2.22d-16 2.23d-308 *
101* *
102********************************************************************************
103* *
104* The program returns a quiet NaN for singularities and infinity when *
105* overflow occurs. The computation is believed to be free of undeserved *
106* underflow and overflow. *
107* *
108* References: *
109* *
110* [1] "An Overview of Software Development for Special Functions", *
111* W. J. Cody, Lecture Notes in Mathematics, 506, Numerical Analysis *
112* Dundee, 1975, G. A. Watson (ed.), Springer Verlag, Berlin, 1976. *
113* *
114* [2] Computer Approximations, Hart, Et. Al., Wiley and sons, New York *
115* 1968. *
116* *
117*******************************************************************************/
118
119#include "math.h"
120#include "fenv.h"
121#include "fp_private.h"
122#include "fenv_private.h"
123
124/*******************************************************************************
125* Functions needed for the computation. *
126*******************************************************************************/
127
128/* the following fp.h functions are used: */
129/* exp, log, sin, __fpclassifyd and nan; */
130/* the following environmental functions are used: */
131/* __setflm. */
132
133
134/*******************************************************************************
135* Coefficients for P in gamma approximation over [1,2] in decreasing order.*
136*******************************************************************************/
137
138static const double p[8] = { -1.71618513886549492533811e+0,
139 2.47656508055759199108314e+1,
140 -3.79804256470945635097577e+2,
141 6.29331155312818442661052e+2,
142 8.66966202790413211295064e+2,
143 -3.14512729688483675254357e+4,
144 -3.61444134186911729807069e+4,
145 6.64561438202405440627855e+4 };
146
147/*******************************************************************************
148* Coefficients for Q in gamma approximation over [1,2] in decreasing order.*
149*******************************************************************************/
150
151static const double q[8] = { -3.08402300119738975254353e+1,
152 3.15350626979604161529144e+2,
153 -1.01515636749021914166146e+3,
154 -3.10777167157231109440444e+3,
155 2.25381184209801510330112e+4,
156 4.75584627752788110767815e+3,
157 -1.34659959864969306392456e+5,
158 -1.15132259675553483497211e+5 };
159
160/*******************************************************************************
161* Coefficients for minimax approximation over [12, INF]. *
162*******************************************************************************/
163
164static const double c[7] = { -1.910444077728e-03,
165 8.4171387781295e-04,
166 -5.952379913043012e-04,
167 7.93650793500350248e-04,
168 -2.777777777777681622553e-03,
169 8.333333333333333331554247e-02,
170 5.7083835261e-03 };
171
172static const double LogSqrt2pi = 0.9189385332046727417803297e+0;
173static const double pi = 3.1415926535897932384626434e+0;
174static const double xbig = 171.624e+0;
175static const double MinimumX = 2.23e-308;
176static const double eps = 2.22e-16;
177static const hexdouble Huge = HEXDOUBLE(0x7FF00000, 0x00000000);
178static const hexdouble MinusHuge = HEXDOUBLE(0xFFF00000, 0x00000000);
179
180#define GAMMA_NAN "42"
181#define SET_INVALID 0x01000000
182
183#pragma STDC FENV_ACCESS ON
184
185double tgamma ( double x )
186 {
187 register int n, parity, i;
188 register double y, y1, result, fact, IsItAnInt, z, numerator,
189 denominator, ysquared, sum;
190 hexdouble OldEnvironment;
191
192 FEGETENVD( OldEnvironment.d ); // save environment, set default
193 FESETENVD( 0.0 );
194
195/*******************************************************************************
196* The next switch will decipher what sort of argument we have. If argument *
197* is SNaN then a QNaN has to be returned and the invalid flag signaled. *
198*******************************************************************************/
199
200 switch ( __fpclassifyd ( x ) )
201 {
202 case FP_NAN:
203 x *= 2.0; /* quiets NaN */
204 FESETENVD( OldEnvironment.d ); // restore caller's environment
205 return x;
206
207 case FP_ZERO:
208 OldEnvironment.i.lo |= FE_DIVBYZERO;
209 FESETENVD( OldEnvironment.d );
210 return copysign( Huge.d, x);
211
212 case FP_INFINITE:
213 if ( x > 0.0 )
214 x = Huge.d;
215 else
216 {
217 x = nan ( GAMMA_NAN );
218 OldEnvironment.i.lo |= SET_INVALID;
219 }
220
221 FESETENVD( OldEnvironment.d );
222 return x;
223
224 default: /* NORMALNUM and DENORMALNUM */
225 break;
226 }
227
228 parity = 0;
229 fact = 1.0;
230 n = 0;
231 y = x;
232
233/*******************************************************************************
234* The argument is negative. *
235*******************************************************************************/
236
237 if ( y <= 0.0 )
238 {
239 y = - x;
240 if ( y < MinimumX )
241 {
242 OldEnvironment.i.lo |= FE_OVERFLOW;
243 FESETENVD( OldEnvironment.d );
244 return MinusHuge.d;
245 }
246 y1 = trunc ( y );
247 IsItAnInt = y - y1;
248 if ( IsItAnInt != 0.0 ) /* is it an integer? */
249 { /* is it odd or even? */
250 if ( y1 != trunc ( y1 * 0.5 ) * 2.0 ) parity = 1;
251 fact = - pi / sin ( pi * IsItAnInt );
252 y += 1.0;
253 }
254 else
255 {
256 OldEnvironment.i.lo |= SET_INVALID;
257 FESETENVD( OldEnvironment.d );
258 return nan ( GAMMA_NAN );
259 }
260 }
261
262/*******************************************************************************
263* The argument is positive. *
264*******************************************************************************/
265
266 if ( y < eps ) /* argument is less than epsilon. */
267 {
268 if ( y >= MinimumX ) /* x is in [MinimumX,eps]. */
269 result = 1.0 / y;
270 else /* othewise, x is in [0,MinimumX). */
271 {
272 OldEnvironment.i.lo |= FE_OVERFLOW;
273 FESETENVD( OldEnvironment.d );
274 return Huge.d;
275 }
276 }
277 else if ( y < 12.0 ) /* argument x is eps < x < 12.0. */
278 {
279 y1 = y;
280 if ( y < 1.0 ) /* x is in (eps, 1.0). */
281 {
282 z = y;
283 y += 1.0;
284 }
285 else /* x is in [1.0,12.0]. */
286 {
287 n = ( int ) y - 1;
288 y -= ( double ) n;
289 z = y - 1.0;
290 }
291 numerator = 0.0;
292 denominator = 1.0;
293 for ( i = 0; i < 8; i++ )
294 {
295 numerator = ( numerator + p[i] ) * z;
296 denominator = denominator * z + q[i];
297 }
298 result = numerator / denominator + 1.0;
299 if ( y1 < y )
300 result /= y1;
301 else if ( y1 > y )
302 {
303 for ( i = 0; i < n; i++ )
304 {
305 result *= y;
306 y += 1.0;
307 }
308 }
309 }
310 else
311 {
312 if ( x <= xbig )
313 {
314 ysquared = y * y;
315 sum = c[6];
316 for ( i = 0; i < 6; i++ )
317 sum = sum / ysquared + c[i];
318 sum = sum / y - y + LogSqrt2pi;
319 sum += ( y - 0.5 ) * log ( y );
320 result = exp ( sum );
321 }
322 else
323 {
324 OldEnvironment.i.lo |= FE_OVERFLOW;
325 FESETENVD( OldEnvironment.d ); // restore caller's environment
326 return Huge.d;
327 }
328 }
329 if ( parity ) result = - result;
330 if ( fact != 1.0 ) result = fact / result;
331 FESETENVD( OldEnvironment.d ); // restore caller's environment
332 return result;
333 }