Ruby  3.1.4p223 (2023-03-30 revision HEAD)
tgamma.c
1 /* tgamma.c - public domain implementation of function tgamma(3m)
2 
3 reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten
4  (New Algorithm handbook in C language) (Gijyutsu hyouron
5  sha, Tokyo, 1991) [in Japanese]
6  http://oku.edu.mie-u.ac.jp/~okumura/algo/
7 */
8 
9 /***********************************************************
10  gamma.c -- Gamma function
11 ***********************************************************/
12 #include "ruby/internal/config.h"
13 #include "ruby/missing.h"
14 #include <math.h>
15 #include <errno.h>
16 
17 #ifndef HAVE_LGAMMA_R
18 
19 #include <errno.h>
20 #define PI 3.14159265358979324 /* $\pi$ */
21 #define LOG_2PI 1.83787706640934548 /* $\log 2\pi$ */
22 #define N 8
23 
24 #define B0 1 /* Bernoulli numbers */
25 #define B1 (-1.0 / 2.0)
26 #define B2 ( 1.0 / 6.0)
27 #define B4 (-1.0 / 30.0)
28 #define B6 ( 1.0 / 42.0)
29 #define B8 (-1.0 / 30.0)
30 #define B10 ( 5.0 / 66.0)
31 #define B12 (-691.0 / 2730.0)
32 #define B14 ( 7.0 / 6.0)
33 #define B16 (-3617.0 / 510.0)
34 
35 static double
36 loggamma(double x) /* the natural logarithm of the Gamma function. */
37 {
38  double v, w;
39 
40  v = 1;
41  while (x < N) { v *= x; x++; }
42  w = 1 / (x * x);
43  return ((((((((B16 / (16 * 15)) * w + (B14 / (14 * 13))) * w
44  + (B12 / (12 * 11))) * w + (B10 / (10 * 9))) * w
45  + (B8 / ( 8 * 7))) * w + (B6 / ( 6 * 5))) * w
46  + (B4 / ( 4 * 3))) * w + (B2 / ( 2 * 1))) / x
47  + 0.5 * LOG_2PI - log(v) - x + (x - 0.5) * log(x);
48 }
49 #endif
50 
51 double tgamma(double x) /* Gamma function */
52 {
53  int sign;
54  if (x == 0.0) { /* Pole Error */
55  errno = ERANGE;
56  return 1/x < 0 ? -HUGE_VAL : HUGE_VAL;
57  }
58  if (isinf(x)) {
59  if (x < 0) goto domain_error;
60  return x;
61  }
62  if (x < 0) {
63  static double zero = 0.0;
64  double i, f;
65  f = modf(-x, &i);
66  if (f == 0.0) { /* Domain Error */
67  domain_error:
68  errno = EDOM;
69  return zero/zero;
70  }
71 #ifndef HAVE_LGAMMA_R
72  sign = (fmod(i, 2.0) != 0.0) ? 1 : -1;
73  return sign * PI / (sin(PI * f) * exp(loggamma(1 - x)));
74 #endif
75  }
76 #ifndef HAVE_LGAMMA_R
77  return exp(loggamma(x));
78 #else
79  x = lgamma_r(x, &sign);
80  return sign * exp(x);
81 #endif
82 }