v/vlib/math/gamma.v

336 lines
7.7 KiB
V
Raw Normal View History

module math
// gamma function computed by Stirling's formula.
// The pair of results must be multiplied together to get the actual answer.
// The multiplication is left to the caller so that, if careful, the caller can avoid
// infinity for 172 <= x <= 180.
// The polynomial is valid for 33 <= x <= 172 larger values are only used
// in reciprocal and produce denormalized floats. The lower precision there
// masks any imprecision in the polynomial.
fn stirling(x f64) (f64, f64) {
if x > 200 {
return inf(1), 1.0
}
sqrt_two_pi := 2.506628274631000502417
max_stirling := 143.01608
mut w := 1.0 / x
w = 1.0 + w * ((((gamma_s[0] * w + gamma_s[1]) * w + gamma_s[2]) * w + gamma_s[3]) * w +
gamma_s[4])
mut y1 := exp(x)
mut y2 := 1.0
if x > max_stirling { // avoid Pow() overflow
v := pow(x, 0.5 * x - 0.25)
y1_ := y1
y1 = v
y2 = v / y1_
} else {
y1 = pow(x, x - 0.5) / y1
}
return y1, f64(sqrt_two_pi) * w * y2
}
// gamma returns the gamma function of x.
//
// special ifs are:
// gamma(+inf) = +inf
// gamma(+0) = +inf
// gamma(-0) = -inf
// gamma(x) = nan for integer x < 0
// gamma(-inf) = nan
// gamma(nan) = nan
pub fn gamma(a f64) f64 {
mut x := a
euler := 0.57721566490153286060651209008240243104215933593992 // A001620
if is_neg_int(x) || is_inf(x, -1) || is_nan(x) {
return nan()
}
if is_inf(x, 1) {
return inf(1)
}
if x == 0.0 {
return copysign(inf(1), x)
}
mut q := abs(x)
mut p := floor(q)
if q > 33 {
if x >= 0 {
y1, y2 := stirling(x)
return y1 * y2
}
// Note: x is negative but (checked above) not a negative integer,
// so x must be small enough to be in range for conversion to i64.
// If |x| were >= 2⁶³ it would have to be an integer.
mut signgam := 1
ip := i64(p)
if (ip & 1) == 0 {
signgam = -1
}
mut z := q - p
if z > 0.5 {
p = p + 1
z = q - p
}
z = q * sin(pi * z)
if z == 0 {
return inf(signgam)
}
sq1, sq2 := stirling(q)
absz := abs(z)
d := absz * sq1 * sq2
if is_inf(d, 0) {
z = pi / absz / sq1 / sq2
} else {
z = pi / d
}
return f64(signgam) * z
}
// Reduce argument
mut z := 1.0
for x >= 3 {
x = x - 1
z = z * x
}
for x < 0 {
if x > -1e-09 {
unsafe {
goto small
}
}
z = z / x
x = x + 1
}
for x < 2 {
if x < 1e-09 {
unsafe {
goto small
}
}
z = z / x
x = x + 1
}
if x == 2 {
return z
}
x = x - 2
p = (((((x * gamma_p[0] + gamma_p[1]) * x + gamma_p[2]) * x + gamma_p[3]) * x +
gamma_p[4]) * x + gamma_p[5]) * x + gamma_p[6]
q = ((((((x * gamma_q[0] + gamma_q[1]) * x + gamma_q[2]) * x + gamma_q[3]) * x +
gamma_q[4]) * x + gamma_q[5]) * x + gamma_q[6]) * x + gamma_q[7]
if true {
return z * p / q
}
small:
if x == 0 {
return inf(1)
}
return z / ((1.0 + euler * x) * x)
}
// log_gamma returns the natural logarithm and sign (-1 or +1) of Gamma(x).
//
// special ifs are:
// log_gamma(+inf) = +inf
// log_gamma(0) = +inf
// log_gamma(-integer) = +inf
// log_gamma(-inf) = -inf
// log_gamma(nan) = nan
pub fn log_gamma(x f64) f64 {
y, _ := log_gamma_sign(x)
return y
}
pub fn log_gamma_sign(a f64) (f64, int) {
mut x := a
ymin := 1.461632144968362245
tiny := exp2(-70)
two52 := exp2(52) // 0x4330000000000000 ~4.5036e+15
two58 := exp2(58) // 0x4390000000000000 ~2.8823e+17
tc := 1.46163214496836224576e+00 // 0x3FF762D86356BE3F
tf := -1.21486290535849611461e-01 // 0xBFBF19B9BCC38A42
// tt := -(tail of tf)
tt := -3.63867699703950536541e-18 // 0xBC50C7CAA48A971F
mut sign := 1
if is_nan(x) {
return x, sign
}
if is_inf(x, 1) {
return x, sign
}
if x == 0.0 {
return inf(1), sign
}
mut neg := false
if x < 0 {
x = -x
neg = true
}
if x < tiny { // if |x| < 2**-70, return -log(|x|)
if neg {
sign = -1
}
return -log(x), sign
}
mut nadj := 0.0
if neg {
if x >= two52 {
// x| >= 2**52, must be -integer
return inf(1), sign
}
t := sin_pi(x)
if t == 0 {
return inf(1), sign
}
nadj = log(pi / abs(t * x))
if t < 0 {
sign = -1
}
}
mut lgamma := 0.0
if x == 1 || x == 2 { // purge off 1 and 2
return 0.0, sign
} else if x < 2 { // use lgamma(x) = lgamma(x+1) - log(x)
mut y := 0.0
mut i := 0
if x <= 0.9 {
lgamma = -log(x)
if x >= (ymin - 1 + 0.27) { // 0.7316 <= x <= 0.9
y = 1.0 - 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
if x >= (ymin + 0.27) { // 1.7316 <= x < 2
y = f64(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
}
}
if i == 0 {
z := y * y
gamma_p1 := lgamma_a[0] + z * (lgamma_a[2] + z * (lgamma_a[4] + z * (lgamma_a[6] +
z * (lgamma_a[8] + z * lgamma_a[10]))))
gamma_p2 := z * (lgamma_a[1] + z * (lgamma_a[3] + z * (lgamma_a[5] + z * (lgamma_a[7] +
z * (lgamma_a[9] + z * lgamma_a[11])))))
p := y * gamma_p1 + gamma_p2
lgamma += (p - 0.5 * y)
} else if i == 1 {
z := y * y
w := z * y
gamma_p1 := lgamma_t[0] + w * (lgamma_t[3] + w * (lgamma_t[6] + w * (lgamma_t[9] +
w * lgamma_t[12]))) // parallel comp
gamma_p2 := lgamma_t[1] + w * (lgamma_t[4] + w * (lgamma_t[7] + w * (lgamma_t[10] +
w * lgamma_t[13])))
gamma_p3 := lgamma_t[2] + w * (lgamma_t[5] + w * (lgamma_t[8] + w * (lgamma_t[11] +
w * lgamma_t[14])))
p := z * gamma_p1 - (tt - w * (gamma_p2 + y * gamma_p3))
lgamma += (tf + p)
} else if i == 2 {
gamma_p1 := y * (lgamma_u[0] + y * (lgamma_u[1] + y * (lgamma_u[2] + y * (lgamma_u[3] +
y * (lgamma_u[4] + y * lgamma_u[5])))))
gamma_p2 := 1.0 + y * (lgamma_v[1] + y * (lgamma_v[2] + y * (lgamma_v[3] +
y * (lgamma_v[4] + y * lgamma_v[5]))))
lgamma += (-0.5 * y + gamma_p1 / gamma_p2)
}
} else if x < 8 { // 2 <= x < 8
i := int(x)
y := x - f64(i)
p := y * (lgamma_s[0] + y * (lgamma_s[1] + y * (lgamma_s[2] + y * (lgamma_s[3] +
y * (lgamma_s[4] + y * (lgamma_s[5] + y * lgamma_s[6]))))))
q := 1.0 + y * (lgamma_r[1] + y * (lgamma_r[2] + y * (lgamma_r[3] + y * (lgamma_r[4] +
y * (lgamma_r[5] + y * lgamma_r[6])))))
lgamma = 0.5 * y + p / q
mut z := 1.0 // lgamma(1+s) = log(s) + lgamma(s)
if i == 7 {
z *= (y + 6)
z *= (y + 5)
z *= (y + 4)
z *= (y + 3)
z *= (y + 2)
lgamma += log(z)
} else if i == 6 {
z *= (y + 5)
z *= (y + 4)
z *= (y + 3)
z *= (y + 2)
lgamma += log(z)
} else if i == 5 {
z *= (y + 4)
z *= (y + 3)
z *= (y + 2)
lgamma += log(z)
} else if i == 4 {
z *= (y + 3)
z *= (y + 2)
lgamma += log(z)
} else if i == 3 {
z *= (y + 2)
lgamma += log(z)
}
} else if x < two58 { // 8 <= x < 2**58
t := log(x)
z := 1.0 / x
y := z * z
w := lgamma_w[0] + z * (lgamma_w[1] + y * (lgamma_w[2] + y * (lgamma_w[3] +
y * (lgamma_w[4] + y * (lgamma_w[5] + y * lgamma_w[6])))))
lgamma = (x - 0.5) * (t - 1.0) + w
} else { // 2**58 <= x <= Inf
lgamma = x * (log(x) - 1.0)
}
if neg {
lgamma = nadj - lgamma
}
return lgamma, sign
}
// sin_pi(x) is a helper function for negative x
fn sin_pi(x_ f64) f64 {
mut x := x_
two52 := exp2(52) // 0x4330000000000000 ~4.5036e+15
two53 := exp2(53) // 0x4340000000000000 ~9.0072e+15
if x < 0.25 {
return -sin(pi * x)
}
// argument reduction
mut z := floor(x)
mut n := 0
if z != x { // inexact
x = mod(x, 2)
n = int(x * 4)
} else {
if x >= two53 { // x must be even
x = 0
n = 0
} else {
if x < two52 {
z = x + two52 // exact
}
n = 1 & int(f64_bits(z))
x = f64(n)
n <<= 2
}
}
if n == 0 {
x = sin(pi * x)
} else if n == 1 || n == 2 {
x = cos(pi * (0.5 - x))
} else if n == 3 || n == 4 {
x = sin(pi * (1.0 - x))
} else if n == 5 || n == 6 {
x = -cos(pi * (x - 1.5))
} else {
x = sin(pi * (x - 2))
}
return -x
}