lgamma function
Implementation
Lgamma lgamma(double x) {
const double Ymin = 1.461632144968362245;
final double Two52 = 4.5036e+15; // 0x4330000000000000 ~4.5036e+15
// final double Two53 = 9.0072e+15; // 0x4340000000000000 ~9.0072e+15
final double Two58 = 2.8823e+17; // 0x4390000000000000 ~2.8823e+17
const double Tiny = 8.47033e-22; // 0x3b90000000000000 ~8.47033e-22
const double Tc = 1.46163214496836224576e+00; // 0x3FF762D86356BE3F
const double Tf = -1.21486290535849611461e-01; // 0xBFBF19B9BCC38A42
// Tt = -(tail of Tf)
const double Tt = -3.63867699703950536541e-18; // 0xBC50C7CAA48A971F
// special cases
int sign = 1;
if (x.isNaN) {
return Lgamma(x, sign);
}
if (x.isInfinite) {
return Lgamma(x, sign);
}
if (x == 0.0) {
return Lgamma(double.infinity, sign);
}
bool neg = false;
if (x < 0) {
x = -x;
neg = true;
}
if (x < Tiny) {
// if |x| < 2**-70, return -log(|x|)
if (neg) {
sign = -1;
}
return Lgamma(-math.log(x), sign);
}
double nadj = 0;
if (neg) {
if (x >= Two52) {
// |x| >= 2**52, must be -integer
return Lgamma(double.infinity, sign);
}
final double t = sinPi(x);
if (t == 0) {
return Lgamma(double.infinity, sign);
}
nadj = math.log(math.pi / (t * x).abs());
if (t < 0) {
sign = -1;
}
}
// purge off 1 and 2
if (x == 1.0 || x == 2.0) {
return Lgamma(0.0, sign);
}
double lgamma;
// use lgamma(x) = lgamma(x+1) - log(x)
if (x < 2) {
double y;
int i;
if (x <= 0.9) {
lgamma = -math.log(x);
if (x >= (Ymin - 1 + 0.27)) {
// 0.7316 <= x <= 0.9
y = 1 - x;
i = 0;
} else if (x >= (Ymin - 1 - 0.27)) {
// 0.2316 <= x < 0.7316
y = x - (Tc - 1);
i = 1;
} else {
// 0 < x < 0.2316
y = x;
i = 2;
}
} else {
lgamma = 0.0;
if (x >= (Ymin + 0.27)) {
// 1.7316 <= x < 2
y = 2 - x;
i = 0;
} else if (x >= (Ymin - 0.27)) {
// 1.2316 <= x < 1.7316
y = x - Tc;
i = 1;
} else {
// 0.9 < x < 1.2316
y = x - 1;
i = 2;
}
}
switch (i) {
case 0:
final double z = y * y;
final double p1 = _lgamA[0] +
z *
(_lgamA[2] +
z *
(_lgamA[4] +
z *
(_lgamA[6] +
z * (_lgamA[8] + z * _lgamA[10]))));
final double p2 = z *
(_lgamA[1] +
z *
(_lgamA[3] +
z *
(_lgamA[5] +
z *
(_lgamA[7] +
z * (_lgamA[9] + z * _lgamA[11])))));
final double p = y * p1 + p2;
lgamma += (p - 0.5 * y);
break;
case 1:
final double z = y * y;
final double w = z * y;
final double p1 = _lgamT[0] +
w *
(_lgamT[3] +
w *
(_lgamT[6] +
w * (_lgamT[9] + w * _lgamT[12]))); // parallel comp
final double p2 = _lgamT[1] +
w *
(_lgamT[4] +
w * (_lgamT[7] + w * (_lgamT[10] + w * _lgamT[13])));
final double p3 = _lgamT[2] +
w *
(_lgamT[5] +
w * (_lgamT[8] + w * (_lgamT[11] + w * _lgamT[14])));
final double p = z * p1 - (Tt - w * (p2 + y * p3));
lgamma += (Tf + p);
break;
case 2:
final double p1 = y *
(_lgamU[0] +
y *
(_lgamU[1] +
y *
(_lgamU[2] +
y *
(_lgamU[3] +
y * (_lgamU[4] + y * _lgamU[5])))));
final double p2 = 1 +
y *
(_lgamV[1] +
y *
(_lgamV[2] +
y * (_lgamV[3] + y * (_lgamV[4] + y * _lgamV[5]))));
lgamma += (-0.5 * y + p1 / p2);
break;
}
} else if (x < 8) {
// 2 <= x < 8
int i = x.toInt();
final double y = x - i.toDouble();
final double p = y *
(_lgamS[0] +
y *
(_lgamS[1] +
y *
(_lgamS[2] +
y *
(_lgamS[3] +
y *
(_lgamS[4] +
y *
(_lgamS[5] +
y * _lgamS[6]))))));
final double q = 1 +
y *
(_lgamR[1] +
y *
(_lgamR[2] +
y *
(_lgamR[3] +
y *
(_lgamR[4] +
y * (_lgamR[5] + y * _lgamR[6])))));
lgamma = 0.5 * y + p / q;
double z = 1.0; // Lgamma(1+s) = Log(s) + Lgamma(s)
for (; i >= 3; i--) {
switch (i) {
case 7:
z *= (y + 6);
break;
case 6:
z *= (y + 5);
break;
case 5:
z *= (y + 4);
break;
case 4:
z *= (y + 3);
break;
case 3:
z *= (y + 2);
lgamma += math.log(z);
break;
}
}
} else if (x < Two58) {
// 8 <= x < 2**58
final double t = math.log(x);
final double z = 1 / x;
final double y = z * z;
final double w = _lgamW[0] +
z *
(_lgamW[1] +
y *
(_lgamW[2] +
y *
(_lgamW[3] +
y *
(_lgamW[4] +
y * (_lgamW[5] + y * _lgamW[6])))));
lgamma = (x - 0.5) * (t - 1) + w;
} else {
// 2**58 <= x <= Inf
lgamma = x * (math.log(x) - 1);
}
if (neg) {
lgamma = nadj - lgamma;
}
return Lgamma(lgamma, sign);
}
