-
Notifications
You must be signed in to change notification settings - Fork 0
/
index.js
67 lines (60 loc) · 1.63 KB
/
index.js
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
// transliterated from the python snippet here:
// http://en.wikipedia.org/wiki/Lanczos_approximation
var g = 7;
var p = [
0.99999999999980993,
676.5203681218851,
-1259.1392167224028,
771.32342877765313,
-176.61502916214059,
12.507343278686905,
-0.13857109526572012,
9.9843695780195716e-6,
1.5056327351493116e-7
];
var g_ln = 607/128;
var p_ln = [
0.99999999999999709182,
57.156235665862923517,
-59.597960355475491248,
14.136097974741747174,
-0.49191381609762019978,
0.33994649984811888699e-4,
0.46523628927048575665e-4,
-0.98374475304879564677e-4,
0.15808870322491248884e-3,
-0.21026444172410488319e-3,
0.21743961811521264320e-3,
-0.16431810653676389022e-3,
0.84418223983852743293e-4,
-0.26190838401581408670e-4,
0.36899182659531622704e-5
];
// Spouge approximation (suitable for large arguments)
function lngamma(z) {
if(z < 0) return Number('0/0');
var x = p_ln[0];
for(var i = p_ln.length - 1; i > 0; --i) x += p_ln[i] / (z + i);
var t = z + g_ln + 0.5;
return .5*Math.log(2*Math.PI)+(z+.5)*Math.log(t)-t+Math.log(x)-Math.log(z);
}
module.exports = function gamma (z) {
if (z < 0.5) {
return Math.PI / (Math.sin(Math.PI * z) * gamma(1 - z));
}
else if(z > 100) return Math.exp(lngamma(z));
else {
z -= 1;
var x = p[0];
for (var i = 1; i < g + 2; i++) {
x += p[i] / (z + i);
}
var t = z + g + 0.5;
return Math.sqrt(2 * Math.PI)
* Math.pow(t, z + 0.5)
* Math.exp(-t)
* x
;
}
};
module.exports.log = lngamma;