Skip to content

Commit ed79cc2

Browse files
youknowoneyouknowone
committed
tgamma.c
1 parent 7fb326e commit ed79cc2

File tree

1 file changed

+316
-0
lines changed

1 file changed

+316
-0
lines changed

‎tgamma.c

+316
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,316 @@
1+
/*
2+
* Implementation of the gamma function based on CPython's implementation.
3+
* This is a standalone version without dependencies on CPython.
4+
*/
5+
6+
#include <math.h>
7+
#include <errno.h>
8+
#include <float.h>
9+
#include <assert.h>
10+
#include <stdbool.h>
11+
12+
/* Constants and definitions needed for gamma calculation */
13+
static const double pi = 3.141592653589793238462643383279502884197;
14+
static const double logpi = 1.144729885849400174143427351353058711647;
15+
16+
/* Check if a value is finite */
17+
#ifndef Py_IS_FINITE
18+
#define Py_IS_FINITE(x) isfinite(x)
19+
#endif
20+
21+
/* Check if a value is infinity */
22+
#ifndef Py_IS_INFINITY
23+
#define Py_IS_INFINITY(x) isinf(x)
24+
#endif
25+
26+
/* Check if a value is NaN */
27+
#ifndef Py_IS_NAN
28+
#define Py_IS_NAN(x) isnan(x)
29+
#endif
30+
31+
/* Define HUGE_VAL if not available */
32+
#ifndef Py_HUGE_VAL
33+
#define Py_HUGE_VAL HUGE_VAL
34+
#endif
35+
36+
/* Define mathematical constants required for Lanczos approximation */
37+
#define LANCZOS_N 13
38+
static const double lanczos_g = 6.024680040776729583740234375;
39+
static const double lanczos_g_minus_half = 5.524680040776729583740234375;
40+
static const double lanczos_num_coeffs[LANCZOS_N] = {
41+
23531376880.410759688572007674451636754734846804940,
42+
42919803642.649098768957899047001988850926355848959,
43+
35711959237.355668049440185451547166705960488635843,
44+
17921034426.037209699919755754458931112671403265390,
45+
6039542586.3520280050642916443072979210699388420708,
46+
1439720407.3117216736632230727949123939715485786772,
47+
248874557.86205415651146038641322942321632125127801,
48+
31426415.585400194380614231628318205362874684987640,
49+
2876370.6289353724412254090516208496135991145378768,
50+
186056.26539522349504029498971604569928220784236328,
51+
8071.6720023658162106380029022722506138218516325024,
52+
210.82427775157934587250973392071336271166969580291,
53+
2.5066282746310002701649081771338373386264310793408
54+
};
55+
56+
/* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */
57+
static const double lanczos_den_coeffs[LANCZOS_N] = {
58+
0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0,
59+
13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0};
60+
61+
/* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */
62+
#define NGAMMA_INTEGRAL 23
63+
static const double gamma_integral[NGAMMA_INTEGRAL] = {
64+
1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,
65+
3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,
66+
1307674368000.0, 20922789888000.0, 355687428096000.0,
67+
6402373705728000.0, 121645100408832000.0, 2432902008176640000.0,
68+
51090942171709440000.0, 1124000727777607680000.0,
69+
};
70+
71+
/* Helper function for sin(pi*x) used in gamma calculation */
72+
static double m_sinpi(double x)
73+
{
74+
double y, r;
75+
int n;
76+
/* this function should only ever be called for finite arguments */
77+
assert(Py_IS_FINITE(x));
78+
y = fmod(fabs(x), 2.0);
79+
n = (int)round(2.0*y);
80+
assert(0 <= n && n <= 4);
81+
switch (n) {
82+
case 0:
83+
r = sin(pi*y);
84+
break;
85+
case 1:
86+
r = cos(pi*(y-0.5));
87+
break;
88+
case 2:
89+
/* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give
90+
-0.0 instead of 0.0 when y == 1.0. */
91+
r = sin(pi*(1.0-y));
92+
break;
93+
case 3:
94+
r = -cos(pi*(y-1.5));
95+
break;
96+
case 4:
97+
r = sin(pi*(y-2.0));
98+
break;
99+
default:
100+
/* Should never happen due to the assert above */
101+
r = 0.0;
102+
}
103+
return copysign(1.0, x)*r;
104+
}
105+
106+
/* Implementation of Lanczos sum for gamma function */
107+
static double lanczos_sum(double x)
108+
{
109+
double num = 0.0, den = 0.0;
110+
int i;
111+
assert(x > 0.0);
112+
/* evaluate the rational function lanczos_sum(x). For large
113+
x, the obvious algorithm risks overflow, so we instead
114+
rescale the denominator and numerator of the rational
115+
function by x**(1-LANCZOS_N) and treat this as a
116+
rational function in 1/x. This also reduces the error for
117+
larger x values. The choice of cutoff point (5.0 below) is
118+
somewhat arbitrary; in tests, smaller cutoff values than
119+
this resulted in lower accuracy. */
120+
if (x < 5.0) {
121+
for (i = LANCZOS_N; --i >= 0; ) {
122+
num = num * x + lanczos_num_coeffs[i];
123+
den = den * x + lanczos_den_coeffs[i];
124+
}
125+
}
126+
else {
127+
for (i = 0; i < LANCZOS_N; i++) {
128+
num = num / x + lanczos_num_coeffs[i];
129+
den = den / x + lanczos_den_coeffs[i];
130+
}
131+
}
132+
return num/den;
133+
}
134+
135+
/**
136+
* Computes the gamma function value for x.
137+
*
138+
* The implementation is based on the Lanczos approximation with parameters
139+
* N=13 and g=6.024680040776729583740234375, which provides excellent accuracy
140+
* across the domain of the function.
141+
*
142+
* @param x The input value
143+
* @return The gamma function value
144+
*
145+
* Special cases:
146+
* - If x is NaN, returns NaN
147+
* - If x is +Inf, returns +Inf
148+
* - If x is -Inf, sets errno to EDOM and returns NaN
149+
* - If x is 0, sets errno to EDOM and returns +/-Inf (with the sign of x)
150+
* - If x is a negative integer, sets errno to EDOM and returns NaN
151+
*/
152+
double tgamma(double x) {
153+
double absx, r, y, z, sqrtpow;
154+
155+
/* special cases */
156+
if (!Py_IS_FINITE(x)) {
157+
if (Py_IS_NAN(x) || x > 0.0)
158+
return x; /* tgamma(nan) = nan, tgamma(inf) = inf */
159+
else {
160+
errno = EDOM;
161+
return NAN; /* tgamma(-inf) = nan, invalid */
162+
}
163+
}
164+
if (x == 0.0) {
165+
errno = EDOM;
166+
/* tgamma(+-0.0) = +-inf, divide-by-zero */
167+
return copysign(INFINITY, x);
168+
}
169+
170+
/* integer arguments */
171+
if (x == floor(x)) {
172+
if (x < 0.0) {
173+
errno = EDOM; /* tgamma(n) = nan, invalid for */
174+
return NAN; /* negative integers n */
175+
}
176+
if (x <= NGAMMA_INTEGRAL)
177+
return gamma_integral[(int)x - 1];
178+
}
179+
absx = fabs(x);
180+
181+
/* tiny arguments: tgamma(x) ~ 1/x for x near 0 */
182+
if (absx < 1e-20) {
183+
r = 1.0/x;
184+
if (Py_IS_INFINITY(r))
185+
errno = ERANGE;
186+
return r;
187+
}
188+
189+
/* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for
190+
x > 200, and underflows to +-0.0 for x < -200, not a negative
191+
integer. */
192+
if (absx > 200.0) {
193+
if (x < 0.0) {
194+
return 0.0/m_sinpi(x);
195+
}
196+
else {
197+
errno = ERANGE;
198+
return Py_HUGE_VAL;
199+
}
200+
}
201+
202+
y = absx + lanczos_g_minus_half;
203+
/* compute error in sum */
204+
if (absx > lanczos_g_minus_half) {
205+
/* note: the correction can be foiled by an optimizing
206+
compiler that (incorrectly) thinks that an expression like
207+
a + b - a - b can be optimized to 0.0. This shouldn't
208+
happen in a standards-conforming compiler. */
209+
double q = y - absx;
210+
z = q - lanczos_g_minus_half;
211+
}
212+
else {
213+
double q = y - lanczos_g_minus_half;
214+
z = q - absx;
215+
}
216+
z = z * lanczos_g / y;
217+
if (x < 0.0) {
218+
r = -pi / m_sinpi(absx) / absx * exp(y) / lanczos_sum(absx);
219+
r -= z * r;
220+
if (absx < 140.0) {
221+
r /= pow(y, absx - 0.5);
222+
}
223+
else {
224+
sqrtpow = pow(y, absx / 2.0 - 0.25);
225+
r /= sqrtpow;
226+
r /= sqrtpow;
227+
}
228+
}
229+
else {
230+
r = lanczos_sum(absx) / exp(y);
231+
r += z * r;
232+
if (absx < 140.0) {
233+
r *= pow(y, absx - 0.5);
234+
}
235+
else {
236+
sqrtpow = pow(y, absx / 2.0 - 0.25);
237+
r *= sqrtpow;
238+
r *= sqrtpow;
239+
}
240+
}
241+
if (Py_IS_INFINITY(r))
242+
errno = ERANGE;
243+
return r;
244+
}
245+
246+
/**
247+
* Computes the natural logarithm of the absolute value of the gamma function.
248+
*
249+
* @param x The input value
250+
* @return The log of the absolute gamma function value
251+
*
252+
* Special cases:
253+
* - If x is NaN, returns NaN
254+
* - If x is +/-Inf, returns +Inf
255+
* - If x is a non-positive integer, sets errno to EDOM and returns +Inf
256+
* - If x is 1 or 2, returns 0.0
257+
*/
258+
double lgamma(double x) {
259+
double r;
260+
double absx;
261+
262+
/* special cases */
263+
if (!Py_IS_FINITE(x)) {
264+
if (Py_IS_NAN(x))
265+
return x; /* lgamma(nan) = nan */
266+
else
267+
return HUGE_VAL; /* lgamma(+-inf) = +inf */
268+
}
269+
270+
/* integer arguments */
271+
if (x == floor(x) && x <= 2.0) {
272+
if (x <= 0.0) {
273+
errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */
274+
return HUGE_VAL; /* integers n <= 0 */
275+
}
276+
else {
277+
return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */
278+
}
279+
}
280+
281+
absx = fabs(x);
282+
/* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */
283+
if (absx < 1e-20)
284+
return -log(absx);
285+
286+
/* Lanczos' formula. We could save a fraction of a ulp in accuracy by
287+
having a second set of numerator coefficients for lanczos_sum that
288+
absorbed the exp(-lanczos_g) term, and throwing out the lanczos_g
289+
subtraction below; it's probably not worth it. */
290+
r = log(lanczos_sum(absx)) - lanczos_g;
291+
r += (absx - 0.5) * (log(absx + lanczos_g - 0.5) - 1);
292+
if (x < 0.0)
293+
/* Use reflection formula to get value for negative x. */
294+
r = logpi - log(fabs(m_sinpi(absx))) - log(absx) - r;
295+
if (Py_IS_INFINITY(r))
296+
errno = ERANGE;
297+
return r;
298+
}
299+
300+
#include <stdio.h>
301+
#include <stdint.h>
302+
303+
union Result {
304+
double d;
305+
uint64_t u;
306+
};
307+
308+
int main() {
309+
union Result result;
310+
result.d = tgamma(-3.8510064710745118);
311+
printf("The result of tgamma(-3.8510064710745118) is: %f\n", result.d);
312+
313+
// Print the binary representation of a double
314+
printf("Bit representation of result: %llu\n", result.u);
315+
return 0;
316+
}

0 commit comments

Comments
 (0)