blob: a28a71b0e03fb54963f3e674fff9090d89699144 (
plain) (
blame)
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
|
#include "ymath.h"
double Exp2(double x) {
return pow(2.0, x);
}
float Exp2f(float x) {
return powf(2.0f, x);
}
#ifdef _MSC_VER
double Erf(double x) {
static constexpr double _M_2_SQRTPI = 1.12837916709551257390;
static constexpr double eps = 1.0e-7;
if (fabs(x) >= 3.75) {
return x > 0 ? 1.0 : -1.0;
}
double r = _M_2_SQRTPI * x;
double f = r;
for (int i = 1;; ++i) {
r *= -x * x / i;
f += r / (2 * i + 1);
if (fabs(r) < eps * (2 * i + 1)) {
break;
}
}
return f;
}
#endif // _MSC_VER
double LogGammaImpl(double x) {
static constexpr double lnSqrt2Pi = 0.91893853320467274178; // log(sqrt(2.0 * PI))
static constexpr double coeff9 = 1.0 / 1188.0;
static constexpr double coeff7 = -1.0 / 1680.0;
static constexpr double coeff5 = 1.0 / 1260.0;
static constexpr double coeff3 = -1.0 / 360.0;
static constexpr double coeff1 = 1.0 / 12.0;
if ((x == 1.0) || (x == 2.0)) {
return 0.0; // 0! = 1
}
double bonus = 0.0;
while (x < 3.0) {
bonus -= log(x);
x += 1.0;
}
double lnX = log(x);
double sqrXInv = 1.0 / (x * x);
double res = coeff9 * sqrXInv + coeff7;
res = res * sqrXInv + coeff5;
res = res * sqrXInv + coeff3;
res = res * sqrXInv + coeff1;
res /= x;
res += x * lnX - x + lnSqrt2Pi - 0.5 * lnX;
return res + bonus;
}
|
