gamma.c
Go to the documentation of this file.
1 /*
2  * Gamma function
3  *
4  *
5  *
6  * SYNOPSIS:
7  *
8  * double x, y, Gamma();
9  *
10  * y = Gamma( x );
11  *
12  *
13  *
14  * DESCRIPTION:
15  *
16  * Returns Gamma function of the argument. The result is
17  * correctly signed.
18  *
19  * Arguments |x| <= 34 are reduced by recurrence and the function
20  * approximated by a rational function of degree 6/7 in the
21  * interval (2,3). Large arguments are handled by Stirling's
22  * formula. Large negative arguments are made positive using
23  * a reflection formula.
24  *
25  *
26  * ACCURACY:
27  *
28  * Relative error:
29  * arithmetic domain # trials peak rms
30  * IEEE -170,-33 20000 2.3e-15 3.3e-16
31  * IEEE -33, 33 20000 9.4e-16 2.2e-16
32  * IEEE 33, 171.6 20000 2.3e-15 3.2e-16
33  *
34  * Error for arguments outside the test range will be larger
35  * owing to error amplification by the exponential function.
36  *
37  */
38 
39 /* lgam()
40  *
41  * Natural logarithm of Gamma function
42  *
43  *
44  *
45  * SYNOPSIS:
46  *
47  * double x, y, lgam();
48  *
49  * y = lgam( x );
50  *
51  *
52  *
53  * DESCRIPTION:
54  *
55  * Returns the base e (2.718...) logarithm of the absolute
56  * value of the Gamma function of the argument.
57  *
58  * For arguments greater than 13, the logarithm of the Gamma
59  * function is approximated by the logarithmic version of
60  * Stirling's formula using a polynomial approximation of
61  * degree 4. Arguments between -33 and +33 are reduced by
62  * recurrence to the interval [2,3] of a rational approximation.
63  * The cosecant reflection formula is employed for arguments
64  * less than -33.
65  *
66  * Arguments greater than MAXLGM return INFINITY and an error
67  * message. MAXLGM = 2.556348e305 for IEEE arithmetic.
68  *
69  *
70  *
71  * ACCURACY:
72  *
73  *
74  * arithmetic domain # trials peak rms
75  * IEEE 0, 3 28000 5.4e-16 1.1e-16
76  * IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17
77  * The error criterion was relative when the function magnitude
78  * was greater than one but absolute when it was less than one.
79  *
80  * The following test used the relative error criterion, though
81  * at certain points the relative error could be much higher than
82  * indicated.
83  * IEEE -200, -4 10000 4.8e-16 1.3e-16
84  *
85  */
86 
87 /*
88  * Cephes Math Library Release 2.2: July, 1992
89  * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
90  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
91  */
92 
93 
94 #include "mconf.h"
95 
96 static double P[] = {
97  1.60119522476751861407E-4,
98  1.19135147006586384913E-3,
99  1.04213797561761569935E-2,
100  4.76367800457137231464E-2,
101  2.07448227648435975150E-1,
102  4.94214826801497100753E-1,
103  9.99999999999999996796E-1
104 };
105 
106 static double Q[] = {
107  -2.31581873324120129819E-5,
108  5.39605580493303397842E-4,
109  -4.45641913851797240494E-3,
110  1.18139785222060435552E-2,
111  3.58236398605498653373E-2,
112  -2.34591795718243348568E-1,
113  7.14304917030273074085E-2,
114  1.00000000000000000320E0
115 };
116 
117 #define MAXGAM 171.624376956302725
118 static double LOGPI = 1.14472988584940017414;
119 
120 /* Stirling's formula for the Gamma function */
121 static double STIR[5] = {
122  7.87311395793093628397E-4,
123  -2.29549961613378126380E-4,
124  -2.68132617805781232825E-3,
125  3.47222221605458667310E-3,
126  8.33333333333482257126E-2,
127 };
128 
129 #define MAXSTIR 143.01608
130 static double SQTPI = 2.50662827463100050242E0;
131 
132 extern double MAXLOG;
133 static double stirf(double);
134 
135 /* Gamma function computed by Stirling's formula.
136  * The polynomial STIR is valid for 33 <= x <= 172.
137  */
138 static double stirf(double x)
139 {
140  double y, w, v;
141 
142  if (x >= MAXGAM) {
143  return (INFINITY);
144  }
145  w = 1.0 / x;
146  w = 1.0 + w * polevl(w, STIR, 4);
147  y = exp(x);
148  if (x > MAXSTIR) { /* Avoid overflow in pow() */
149  v = pow(x, 0.5 * x - 0.25);
150  y = v * (v / y);
151  }
152  else {
153  y = pow(x, x - 0.5) / y;
154  }
155  y = SQTPI * y * w;
156  return (y);
157 }
158 
159 
160 double Gamma(double x)
161 {
162  double p, q, z;
163  int i;
164  int sgngam = 1;
165 
166  if (!cephes_isfinite(x)) {
167  return x;
168  }
169  q = fabs(x);
170 
171  if (q > 33.0) {
172  if (x < 0.0) {
173  p = floor(q);
174  if (p == q) {
175  gamnan:
176  sf_error("Gamma", SF_ERROR_OVERFLOW, NULL);
177  return (INFINITY);
178  }
179  i = p;
180  if ((i & 1) == 0)
181  sgngam = -1;
182  z = q - p;
183  if (z > 0.5) {
184  p += 1.0;
185  z = q - p;
186  }
187  z = q * sin(M_PI * z);
188  if (z == 0.0) {
189  return (sgngam * INFINITY);
190  }
191  z = fabs(z);
192  z = M_PI / (z * stirf(q));
193  }
194  else {
195  z = stirf(x);
196  }
197  return (sgngam * z);
198  }
199 
200  z = 1.0;
201  while (x >= 3.0) {
202  x -= 1.0;
203  z *= x;
204  }
205 
206  while (x < 0.0) {
207  if (x > -1.E-9)
208  goto small;
209  z /= x;
210  x += 1.0;
211  }
212 
213  while (x < 2.0) {
214  if (x < 1.e-9)
215  goto small;
216  z /= x;
217  x += 1.0;
218  }
219 
220  if (x == 2.0)
221  return (z);
222 
223  x -= 2.0;
224  p = polevl(x, P, 6);
225  q = polevl(x, Q, 7);
226  return (z * p / q);
227 
228  small:
229  if (x == 0.0) {
230  goto gamnan;
231  }
232  else
233  return (z / ((1.0 + 0.5772156649015329 * x) * x));
234 }
235 
236 
237 
238 /* A[]: Stirling's formula expansion of log Gamma
239  * B[], C[]: log Gamma function between 2 and 3
240  */
241 static double A[] = {
242  8.11614167470508450300E-4,
243  -5.95061904284301438324E-4,
244  7.93650340457716943945E-4,
245  -2.77777777730099687205E-3,
246  8.33333333333331927722E-2
247 };
248 
249 static double B[] = {
250  -1.37825152569120859100E3,
251  -3.88016315134637840924E4,
252  -3.31612992738871184744E5,
253  -1.16237097492762307383E6,
254  -1.72173700820839662146E6,
255  -8.53555664245765465627E5
256 };
257 
258 static double C[] = {
259  /* 1.00000000000000000000E0, */
260  -3.51815701436523470549E2,
261  -1.70642106651881159223E4,
262  -2.20528590553854454839E5,
263  -1.13933444367982507207E6,
264  -2.53252307177582951285E6,
265  -2.01889141433532773231E6
266 };
267 
268 /* log( sqrt( 2*pi ) ) */
269 static double LS2PI = 0.91893853320467274178;
270 
271 #define MAXLGM 2.556348e305
272 
273 
274 /* Logarithm of Gamma function */
275 double lgam(double x)
276 {
277  int sign;
278  return lgam_sgn(x, &sign);
279 }
280 
281 double lgam_sgn(double x, int *sign)
282 {
283  double p, q, u, w, z;
284  int i;
285 
286  *sign = 1;
287 
288  if (!cephes_isfinite(x))
289  return x;
290 
291  if (x < -34.0) {
292  q = -x;
293  w = lgam_sgn(q, sign);
294  p = floor(q);
295  if (p == q) {
296  lgsing:
297  sf_error("lgam", SF_ERROR_SINGULAR, NULL);
298  return (INFINITY);
299  }
300  i = p;
301  if ((i & 1) == 0)
302  *sign = -1;
303  else
304  *sign = 1;
305  z = q - p;
306  if (z > 0.5) {
307  p += 1.0;
308  z = p - q;
309  }
310  z = q * sin(M_PI * z);
311  if (z == 0.0)
312  goto lgsing;
313  /* z = log(M_PI) - log( z ) - w; */
314  z = LOGPI - log(z) - w;
315  return (z);
316  }
317 
318  if (x < 13.0) {
319  z = 1.0;
320  p = 0.0;
321  u = x;
322  while (u >= 3.0) {
323  p -= 1.0;
324  u = x + p;
325  z *= u;
326  }
327  while (u < 2.0) {
328  if (u == 0.0)
329  goto lgsing;
330  z /= u;
331  p += 1.0;
332  u = x + p;
333  }
334  if (z < 0.0) {
335  *sign = -1;
336  z = -z;
337  }
338  else
339  *sign = 1;
340  if (u == 2.0)
341  return (log(z));
342  p -= 2.0;
343  x = x + p;
344  p = x * polevl(x, B, 5) / p1evl(x, C, 6);
345  return (log(z) + p);
346  }
347 
348  if (x > MAXLGM) {
349  return (*sign * INFINITY);
350  }
351 
352  q = (x - 0.5) * log(x) - x + LS2PI;
353  if (x > 1.0e8)
354  return (q);
355 
356  p = 1.0 / (x * x);
357  if (x >= 1000.0)
358  q += ((7.9365079365079365079365e-4 * p
359  - 2.7777777777777777777778e-3) * p
360  + 0.0833333333333333333333) / x;
361  else
362  q += polevl(p, A, 4) / x;
363  return (q);
364 }
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
P
static double P[]
Definition: gamma.c:96
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
ceres::sin
Jet< T, N > sin(const Jet< T, N > &f)
Definition: jet.h:439
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
Gamma
double Gamma(double x)
Definition: gamma.c:160
sign
const EIGEN_DEVICE_FUNC SignReturnType sign() const
Definition: ArrayCwiseUnaryOps.h:219
B
Definition: test_numpy_dtypes.cpp:299
SF_ERROR_OVERFLOW
@ SF_ERROR_OVERFLOW
Definition: sf_error.h:12
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
lgam
double lgam(double x)
Definition: gamma.c:275
C
static double C[]
Definition: gamma.c:258
SQTPI
static double SQTPI
Definition: gamma.c:130
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
stirf
static double stirf(double)
Definition: gamma.c:138
p1evl
static double p1evl(double x, const double coef[], int N)
Definition: polevl.h:97
lgam_sgn
double lgam_sgn(double x, int *sign)
Definition: gamma.c:281
A
Definition: test_numpy_dtypes.cpp:298
MAXLGM
#define MAXLGM
Definition: gamma.c:271
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:1984
MAXSTIR
#define MAXSTIR
Definition: gamma.c:129
polevl
static double polevl(double x, const double coef[], int N)
Definition: polevl.h:65
MAXLOG
double MAXLOG
Definition: const.c:57
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
y
Scalar * y
Definition: level1_cplx_impl.h:124
E
DiscreteKey E(5, 2)
Eigen::ArrayBase::pow
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_pow_op< typename Derived::Scalar, typename ExponentDerived::Scalar >, const Derived, const ExponentDerived > pow(const Eigen::ArrayBase< Derived > &x, const Eigen::ArrayBase< ExponentDerived > &exponents)
Definition: GlobalFunctions.h:143
Eigen::Quaternion
The quaternion class used to represent 3D orientations and rotations.
Definition: ForwardDeclarations.h:293
mconf.h
cephes_isfinite
#define cephes_isfinite(x)
Definition: mconf.h:101
p
float * p
Definition: Tutorial_Map_using.cpp:9
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
STIR
static double STIR[5]
Definition: gamma.c:121
SF_ERROR_SINGULAR
@ SF_ERROR_SINGULAR
Definition: sf_error.h:10
NULL
#define NULL
Definition: ccolamd.c:609
M_PI
#define M_PI
Definition: mconf.h:117
LS2PI
static double LS2PI
Definition: gamma.c:269
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
floor
const EIGEN_DEVICE_FUNC FloorReturnType floor() const
Definition: ArrayCwiseUnaryOps.h:481
MAXGAM
#define MAXGAM
Definition: gamma.c:117
LOGPI
static double LOGPI
Definition: gamma.c:118


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:34