aboutsummaryrefslogtreecommitdiff
path: root/v_windows/v/vlib/math/gamma.v
diff options
context:
space:
mode:
Diffstat (limited to 'v_windows/v/vlib/math/gamma.v')
-rw-r--r--v_windows/v/vlib/math/gamma.v335
1 files changed, 335 insertions, 0 deletions
diff --git a/v_windows/v/vlib/math/gamma.v b/v_windows/v/vlib/math/gamma.v
new file mode 100644
index 0000000..e0061db
--- /dev/null
+++ b/v_windows/v/vlib/math/gamma.v
@@ -0,0 +1,335 @@
+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
+}