aboutsummaryrefslogtreecommitdiffstats
path: root/util/generic/ymath.cpp
blob: d6a0038b5d03ea8a8b3f649611acbb6fe2406417 (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
#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; 
}