1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
|
extend namespace Math {
public real gamma(real n)
/*
* Stieljes continuing fraction method of computing
* gamma
*/
{
real Stieltjes(real n, int ord, int bits) {
/*
* Compute the continuing fraction coefficients
*/
real[*] StieltjesCF(int len, int bits) {
/* Compute a set of Bernouli numbers
* using the Akiyama-Tanigawa algorithm
*/
real[*] AkiyamaTanigawa(int l, int bits) {
int n = 2 * l + 1;
real[n] t = { imprecise(1, bits) };
real[l] a;
int k = 1;
for (int m = 2; m <= n; m++) {
t[m-1] = 1/imprecise(m, bits);
for (int j = m-1; j >= 1; j--)
t[j-1] = j * (t[j-1] - t[j]);
if ((m & 1) != 0) {
real rk = imprecise(k, bits);
real v = t[0]/((2*rk-1)*(2*rk));
if ((k & 1) == 0)
a[k-1] = -v;
else
a[k-1] = v;
k++;
}
}
return a;
}
real[*] s = AkiyamaTanigawa(len, bits);
real[len,len] m;
for (int n = 0; n < len; n++)
m[n,0] = 0;
for (int n = 0; n < len-1; n++)
m[n,1] = s[n+1]/s[n];
for (int k = 3; k <= len; k++) {
for (int n = 1; n <= len - k + 1; n++) {
real a = m[n,k-3];
real b = m[n,k-2];
real c = m[n-1,k-2];
if ((k & 1) != 0)
m[n-1,k-1] = a + b - c;
else
m[n-1,k-1] = a * b / c;
}
}
m[0,0] = s[0];
return (real[len]) { [k] = m[0,k] };
}
real N = imprecise(n + 1, bits);
real q = N;
real[*] c = StieltjesCF(ord, bits);
real one = imprecise(1, bits);
for (int i = ord; i >= 2; i--)
q = N + c[i-1] / q;
return sqrt(2 * pi_value(bits)/N) * (N/exp(one)) ** N * exp(one/(12*q));
}
/*
* For positive integers, just use factorial
*/
if (is_int(n)) {
if (n <= 0)
raise invalid_argument("gamma of non-positive integer", 0, n);
return (floor(n)-1)!;
}
n = imprecise(n);
int bits = precision(n);
/*
* Use Euler's reflection formula to solve for arguments < -½
*
* Γ(z) = π/(sin(πz)Γ(1-z))
*/
if (n < -0.5) {
int ibits = bits + 15;
real n_ibits = imprecise(n, ibits);
real pi_ibits = pi_value(ibits);
real sin_n_pi = sin(pi_ibits * n_ibits);
real inv_gamma = gamma(1 - n_ibits);
/* convert divide-by-zero into invalid_argument */
real div = sin_n_pi * inv_gamma;
if (div == 0)
raise invalid_argument("gamma of non-positive integer", 0, n);
return imprecise(pi_ibits / div, bits);
}
/*
* Smaller numbers take a lot more coefficients to
* get the desired number of bits. Make the value
* larger, and increase the desired precision to match,
* to make the result converge faster.
*/
if (n < 10000) {
int new_bits = bits + 20 - exponent(n);
real n_new = imprecise(n, new_bits) + 10000;
real g = gamma(n_new);
for (int i = 0; i < 10000; i++) {
n_new -= 1;
g = g / n_new;
}
return imprecise(g, bits);
}
/* This is a rough approximation of the
* number of coefficients in the fraction
* needed to produce the desired precision, it's
* good for any value larger than 10000. Larger numbers
* could use a smaller number of coefficients, but we
* don't know how much smaller
*/
int ord = ceil(bits / 20);
return imprecise(Stieltjes(n-1, ord, bits + 20), bits);
}
public real(real n) Γ = gamma;
}
|
