1 | |
---|
2 | /* @(#)w_gamma.c 5.1 93/09/24 */ |
---|
3 | /* |
---|
4 | * ==================================================== |
---|
5 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
---|
6 | * |
---|
7 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
---|
8 | * Permission to use, copy, modify, and distribute this |
---|
9 | * software is freely granted, provided that this notice |
---|
10 | * is preserved. |
---|
11 | * ==================================================== |
---|
12 | * |
---|
13 | */ |
---|
14 | |
---|
15 | /* BUG: FIXME? |
---|
16 | According to Linux man pages for tgamma, lgamma, and gamma, the gamma |
---|
17 | function was originally defined in BSD as implemented here--the log of the gamma |
---|
18 | function. BSD 4.3 changed the name to lgamma, apparently removing gamma. BSD |
---|
19 | 4.4 re-introduced the gamma name with the more intuitive, without logarithm, |
---|
20 | plain gamma function. The C99 standard apparently wanted to avoid a problem |
---|
21 | with the poorly-named earlier gamma and used tgamma when adding a plain |
---|
22 | gamma function. |
---|
23 | So the current gamma is matching an old, bad definition, and not |
---|
24 | matching a newer, better definition. */ |
---|
25 | /* |
---|
26 | FUNCTION |
---|
27 | <<gamma>>, <<gammaf>>, <<lgamma>>, <<lgammaf>>, <<gamma_r>>, <<gammaf_r>>, <<lgamma_r>>, <<lgammaf_r>>, <<tgamma>>, and <<tgammaf>>---logarithmic and plain gamma functions |
---|
28 | |
---|
29 | INDEX |
---|
30 | gamma |
---|
31 | INDEX |
---|
32 | gammaf |
---|
33 | INDEX |
---|
34 | lgamma |
---|
35 | INDEX |
---|
36 | lgammaf |
---|
37 | INDEX |
---|
38 | gamma_r |
---|
39 | INDEX |
---|
40 | gammaf_r |
---|
41 | INDEX |
---|
42 | lgamma_r |
---|
43 | INDEX |
---|
44 | lgammaf_r |
---|
45 | INDEX |
---|
46 | tgamma |
---|
47 | INDEX |
---|
48 | tgammaf |
---|
49 | |
---|
50 | SYNOPSIS |
---|
51 | #include <math.h> |
---|
52 | double gamma(double <[x]>); |
---|
53 | float gammaf(float <[x]>); |
---|
54 | double lgamma(double <[x]>); |
---|
55 | float lgammaf(float <[x]>); |
---|
56 | double gamma_r(double <[x]>, int *<[signgamp]>); |
---|
57 | float gammaf_r(float <[x]>, int *<[signgamp]>); |
---|
58 | double lgamma_r(double <[x]>, int *<[signgamp]>); |
---|
59 | float lgammaf_r(float <[x]>, int *<[signgamp]>); |
---|
60 | double tgamma(double <[x]>); |
---|
61 | float tgammaf(float <[x]>); |
---|
62 | |
---|
63 | DESCRIPTION |
---|
64 | <<gamma>> calculates |
---|
65 | @tex |
---|
66 | $\mit ln\bigl(\Gamma(x)\bigr)$, |
---|
67 | @end tex |
---|
68 | the natural logarithm of the gamma function of <[x]>. The gamma function |
---|
69 | (<<exp(gamma(<[x]>))>>) is a generalization of factorial, and retains |
---|
70 | the property that |
---|
71 | @ifnottex |
---|
72 | <<exp(gamma(N))>> is equivalent to <<N*exp(gamma(N-1))>>. |
---|
73 | @end ifnottex |
---|
74 | @tex |
---|
75 | $\mit \Gamma(N)\equiv N\times\Gamma(N-1)$. |
---|
76 | @end tex |
---|
77 | Accordingly, the results of the gamma function itself grow very |
---|
78 | quickly. <<gamma>> is defined as |
---|
79 | @tex |
---|
80 | $\mit ln\bigl(\Gamma(x)\bigr)$ rather than simply $\mit \Gamma(x)$ |
---|
81 | @end tex |
---|
82 | @ifnottex |
---|
83 | the natural log of the gamma function, rather than the gamma function |
---|
84 | itself, |
---|
85 | @end ifnottex |
---|
86 | to extend the useful range of results representable. |
---|
87 | |
---|
88 | The sign of the result is returned in the global variable <<signgam>>, |
---|
89 | which is declared in math.h. |
---|
90 | |
---|
91 | <<gammaf>> performs the same calculation as <<gamma>>, but uses and |
---|
92 | returns <<float>> values. |
---|
93 | |
---|
94 | <<lgamma>> and <<lgammaf>> are alternate names for <<gamma>> and |
---|
95 | <<gammaf>>. The use of <<lgamma>> instead of <<gamma>> is a reminder |
---|
96 | that these functions compute the log of the gamma function, rather |
---|
97 | than the gamma function itself. |
---|
98 | |
---|
99 | The functions <<gamma_r>>, <<gammaf_r>>, <<lgamma_r>>, and |
---|
100 | <<lgammaf_r>> are just like <<gamma>>, <<gammaf>>, <<lgamma>>, and |
---|
101 | <<lgammaf>>, respectively, but take an additional argument. This |
---|
102 | additional argument is a pointer to an integer. This additional |
---|
103 | argument is used to return the sign of the result, and the global |
---|
104 | variable <<signgam>> is not used. These functions may be used for |
---|
105 | reentrant calls (but they will still set the global variable <<errno>> |
---|
106 | if an error occurs). |
---|
107 | |
---|
108 | <<tgamma>> and <<tgammaf>> are the "true gamma" functions, returning |
---|
109 | @tex |
---|
110 | $\mit \Gamma(x)$, |
---|
111 | @end tex |
---|
112 | the gamma function of <[x]>--without a logarithm. |
---|
113 | (They are apparently so named because of the prior existence of the old, |
---|
114 | poorly-named <<gamma>> functions which returned the log of gamma up |
---|
115 | through BSD 4.2.) |
---|
116 | |
---|
117 | RETURNS |
---|
118 | Normally, the computed result is returned. |
---|
119 | |
---|
120 | When <[x]> is a nonpositive integer, <<gamma>> returns <<HUGE_VAL>> |
---|
121 | and <<errno>> is set to <<EDOM>>. If the result overflows, <<gamma>> |
---|
122 | returns <<HUGE_VAL>> and <<errno>> is set to <<ERANGE>>. |
---|
123 | |
---|
124 | You can modify this error treatment using <<matherr>>. |
---|
125 | |
---|
126 | PORTABILITY |
---|
127 | Neither <<gamma>> nor <<gammaf>> is ANSI C. It is better not to use either |
---|
128 | of these; use <<lgamma>> or <<tgamma>> instead.@* |
---|
129 | <<lgamma>>, <<lgammaf>>, <<tgamma>>, and <<tgammaf>> are nominally C standard |
---|
130 | in terms of the base return values, although the <<matherr>> error-handling |
---|
131 | is not standard, nor is the <[signgam]> global for <<lgamma>>. |
---|
132 | */ |
---|
133 | |
---|
134 | /* double gamma(double x) |
---|
135 | * Return the logarithm of the Gamma function of x. |
---|
136 | * |
---|
137 | * Method: call gamma_r |
---|
138 | */ |
---|
139 | |
---|
140 | #include "fdlibm.h" |
---|
141 | #include <reent.h> |
---|
142 | #include <errno.h> |
---|
143 | |
---|
144 | #ifndef _DOUBLE_IS_32BITS |
---|
145 | |
---|
146 | #ifdef __STDC__ |
---|
147 | double gamma(double x) |
---|
148 | #else |
---|
149 | double gamma(x) |
---|
150 | double x; |
---|
151 | #endif |
---|
152 | { |
---|
153 | #ifdef _IEEE_LIBM |
---|
154 | return __ieee754_gamma_r(x,&(_REENT_SIGNGAM(_REENT))); |
---|
155 | #else |
---|
156 | double y; |
---|
157 | struct exception exc; |
---|
158 | y = __ieee754_gamma_r(x,&(_REENT_SIGNGAM(_REENT))); |
---|
159 | if(_LIB_VERSION == _IEEE_) return y; |
---|
160 | if(!finite(y)&&finite(x)) { |
---|
161 | #ifndef HUGE_VAL |
---|
162 | #define HUGE_VAL inf |
---|
163 | double inf = 0.0; |
---|
164 | |
---|
165 | SET_HIGH_WORD(inf,0x7ff00000); /* set inf to infinite */ |
---|
166 | #endif |
---|
167 | exc.name = "gamma"; |
---|
168 | exc.err = 0; |
---|
169 | exc.arg1 = exc.arg2 = x; |
---|
170 | if (_LIB_VERSION == _SVID_) |
---|
171 | exc.retval = HUGE; |
---|
172 | else |
---|
173 | exc.retval = HUGE_VAL; |
---|
174 | if(floor(x)==x&&x<=0.0) { |
---|
175 | /* gamma(-integer) or gamma(0) */ |
---|
176 | exc.type = SING; |
---|
177 | if (_LIB_VERSION == _POSIX_) |
---|
178 | errno = EDOM; |
---|
179 | else if (!matherr(&exc)) { |
---|
180 | errno = EDOM; |
---|
181 | } |
---|
182 | } else { |
---|
183 | /* gamma(finite) overflow */ |
---|
184 | exc.type = OVERFLOW; |
---|
185 | if (_LIB_VERSION == _POSIX_) |
---|
186 | errno = ERANGE; |
---|
187 | else if (!matherr(&exc)) { |
---|
188 | errno = ERANGE; |
---|
189 | } |
---|
190 | } |
---|
191 | if (exc.err != 0) |
---|
192 | errno = exc.err; |
---|
193 | return exc.retval; |
---|
194 | } else |
---|
195 | return y; |
---|
196 | #endif |
---|
197 | } |
---|
198 | |
---|
199 | #endif /* defined(_DOUBLE_IS_32BITS) */ |
---|