Ruby
1.9.3p448(2013-06-27revision41675)
Main Page
Modules
Data Structures
Files
File List
Globals
missing
lgamma_r.c
Go to the documentation of this file.
1
/* lgamma_r.c - public domain implementation of 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
#include "
ruby/missing.h
"
12
/***********************************************************
13
gamma.c -- Gamma function
14
***********************************************************/
15
#include <math.h>
16
#include <errno.h>
17
#define PI 3.14159265358979324
/* $\pi$ */
18
#define LOG_2PI 1.83787706640934548
/* $\log 2\pi$ */
19
#define LOG_PI 1.14472988584940017
/* $\log_e \pi$ */
20
#define N 8
21
22
#define B0 1
/* Bernoulli numbers */
23
#define B1 (-1.0 / 2.0)
24
#define B2 ( 1.0 / 6.0)
25
#define B4 (-1.0 / 30.0)
26
#define B6 ( 1.0 / 42.0)
27
#define B8 (-1.0 / 30.0)
28
#define B10 ( 5.0 / 66.0)
29
#define B12 (-691.0 / 2730.0)
30
#define B14 ( 7.0 / 6.0)
31
#define B16 (-3617.0 / 510.0)
32
33
static
double
34
loggamma
(
double
x)
/* the natural logarithm of the Gamma function. */
35
{
36
double
v
, w;
37
38
if
(x == 1.0 || x == 2.0)
return
0.0;
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
50
51
#ifdef __MINGW_ATTRIB_PURE
52
/* get rid of bugs in math.h of mingw */
53
#define modf(_X, _Y) __extension__ ({\
54
double intpart_modf_bug = intpart_modf_bug;\
55
double result_modf_bug = modf((_X), &intpart_modf_bug);\
56
*(_Y) = intpart_modf_bug;\
57
result_modf_bug;\
58
})
59
#endif
60
61
/* the natural logarithm of the absolute value of the Gamma function */
62
double
63
lgamma_r
(
double
x,
int
*signp)
64
{
65
if
(x <= 0) {
66
double
i
, f, s;
67
f = modf(-x, &i);
68
if
(f == 0.0) {
/* pole error */
69
*signp = 1;
70
errno
= ERANGE;
71
return
HUGE_VAL;
72
}
73
*signp = (fmod(i, 2.0) != 0.0) ? 1 : -1;
74
s = sin(
PI
* f);
75
if
(s < 0) s = -s;
76
return
LOG_PI
- log(s) -
loggamma
(1 - x);
77
}
78
*signp = 1;
79
return
loggamma
(x);
80
}
81
Generated on Fri Jun 28 2013 02:34:41 for Ruby by
1.8.3