Adds most of the toolsHEADmaster

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
+}