Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Newer
Older
100644 69 lines (60 sloc) 2.103 kb
9c1d230 committing experimental branch content
Laurent Sansonetti authored
1 /* lgamma_r.c - public domain implementation of error function lgamma_r(3m)
2
3 lgamma_r() is based on gamma(). modified by Tanaka Akira.
4
5 reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten
6 (New Algorithm handbook in C language) (Gijyutsu hyouron
7 sha, Tokyo, 1991) [in Japanese]
8 http://oku.edu.mie-u.ac.jp/~okumura/algo/
9 */
10
11 /***********************************************************
12 gamma.c -- Gamma function
13 ***********************************************************/
14 #include <math.h>
15 #include <errno.h>
16 #define PI 3.14159265358979324 /* $\pi$ */
17 #define LOG_2PI 1.83787706640934548 /* $\log 2\pi$ */
18 #define LOG_PI 1.14472988584940017 /* $\log_e \pi$ */
19 #define N 8
20
21 #define B0 1 /* Bernoulli numbers */
22 #define B1 (-1.0 / 2.0)
23 #define B2 ( 1.0 / 6.0)
24 #define B4 (-1.0 / 30.0)
25 #define B6 ( 1.0 / 42.0)
26 #define B8 (-1.0 / 30.0)
27 #define B10 ( 5.0 / 66.0)
28 #define B12 (-691.0 / 2730.0)
29 #define B14 ( 7.0 / 6.0)
30 #define B16 (-3617.0 / 510.0)
31
32 static double
33 loggamma(double x) /* the natural logarithm of the Gamma function. */
34 {
35 double v, w;
36
37 if (x == 1.0 || x == 2.0) return 0.0;
38
39 v = 1;
40 while (x < N) { v *= x; x++; }
41 w = 1 / (x * x);
42 return ((((((((B16 / (16 * 15)) * w + (B14 / (14 * 13))) * w
43 + (B12 / (12 * 11))) * w + (B10 / (10 * 9))) * w
44 + (B8 / ( 8 * 7))) * w + (B6 / ( 6 * 5))) * w
45 + (B4 / ( 4 * 3))) * w + (B2 / ( 2 * 1))) / x
46 + 0.5 * LOG_2PI - log(v) - x + (x - 0.5) * log(x);
47 }
48
49 /* the natural logarithm of the absolute value of the Gamma function */
50 double
51 lgamma_r(double x, int *signp)
52 {
53 if (x <= 0) {
54 double i, f, s;
55 f = modf(-x, &i);
56 if (f == 0.0) { /* pole error */
57 *signp = 1;
58 errno = ERANGE;
59 return HUGE_VAL;
60 }
61 *signp = (fmod(i, 2.0) != 0.0) ? 1 : -1;
62 s = sin(PI * f);
63 if (s < 0) s = -s;
64 return LOG_PI - log(s) - loggamma(1 - x);
65 }
66 *signp = 1;
67 return loggamma(x);
68 }
Something went wrong with that request. Please try again.