1 files changed, 214 insertions, 0 deletions
diff --git a/v_windows/v/vlib/math/exp.v b/v_windows/v/vlib/math/exp.v
new file mode 100644
index 0000000..7a11eb6
--- /dev/null
+++ b/v_windows/v/vlib/math/exp.v
@@ -0,0 +1,214 @@
+module math
+
+import math.internal
+
+const (
+	f64_max_exp = f64(1024)
+	f64_min_exp = f64(-1021)
+	threshold   = 7.09782712893383973096e+02 // 0x40862E42FEFA39EF
+	ln2_x56     = 3.88162421113569373274e+01 // 0x4043687a9f1af2b1
+	ln2_halfx3  = 1.03972077083991796413e+00 // 0x3ff0a2b23f3bab73
+	ln2_half    = 3.46573590279972654709e-01 // 0x3fd62e42fefa39ef
+	ln2hi       = 6.93147180369123816490e-01 // 0x3fe62e42fee00000
+	ln2lo       = 1.90821492927058770002e-10 // 0x3dea39ef35793c76
+	inv_ln2     = 1.44269504088896338700e+00 // 0x3ff71547652b82fe
+	// scaled coefficients related to expm1
+	expm1_q1    = -3.33333333333331316428e-02 // 0xBFA11111111110F4
+	expm1_q2    = 1.58730158725481460165e-03 // 0x3F5A01A019FE5585
+	expm1_q3    = -7.93650757867487942473e-05 // 0xBF14CE199EAADBB7
+	expm1_q4    = 4.00821782732936239552e-06 // 0x3ED0CFCA86E65239
+	expm1_q5    = -2.01099218183624371326e-07 // 0xBE8AFDB76E09C32D
+)
+
+// exp returns e**x, the base-e exponential of x.
+//
+// special cases are:
+// exp(+inf) = +inf
+// exp(nan) = nan
+// Very large values overflow to 0 or +inf.
+// Very small values underflow to 1.
+pub fn exp(x f64) f64 {
+	log2e := 1.44269504088896338700e+00
+	overflow := 7.09782712893383973096e+02
+	underflow := -7.45133219101941108420e+02
+	near_zero := 1.0 / (1 << 28) // 2**-28
+	// special cases
+	if is_nan(x) || is_inf(x, 1) {
+		return x
+	}
+	if is_inf(x, -1) {
+		return 0.0
+	}
+	if x > overflow {
+		return inf(1)
+	}
+	if x < underflow {
+		return 0.0
+	}
+	if -near_zero < x && x < near_zero {
+		return 1.0 + x
+	}
+	// reduce; computed as r = hi - lo for extra precision.
+	mut k := 0
+	if x < 0 {
+		k = int(log2e * x - 0.5)
+	}
+	if x > 0 {
+		k = int(log2e * x + 0.5)
+	}
+	hi := x - f64(k) * math.ln2hi
+	lo := f64(k) * math.ln2lo
+	// compute
+	return expmulti(hi, lo, k)
+}
+
+// exp2 returns 2**x, the base-2 exponential of x.
+//
+// special cases are the same as exp.
+pub fn exp2(x f64) f64 {
+	overflow := 1.0239999999999999e+03
+	underflow := -1.0740e+03
+	if is_nan(x) || is_inf(x, 1) {
+		return x
+	}
+	if is_inf(x, -1) {
+		return 0
+	}
+	if x > overflow {
+		return inf(1)
+	}
+	if x < underflow {
+		return 0
+	}
+	// argument reduction; x = r×lg(e) + k with |r| ≤ ln(2)/2.
+	// computed as r = hi - lo for extra precision.
+	mut k := 0
+	if x > 0 {
+		k = int(x + 0.5)
+	}
+	if x < 0 {
+		k = int(x - 0.5)
+	}
+	mut t := x - f64(k)
+	hi := t * math.ln2hi
+	lo := -t * math.ln2lo
+	// compute
+	return expmulti(hi, lo, k)
+}
+
+pub fn ldexp(x f64, e int) f64 {
+	if x == 0.0 {
+		return x
+	} else {
+		mut y, ex := frexp(x)
+		mut e2 := f64(e + ex)
+		if e2 >= math.f64_max_exp {
+			y *= pow(2.0, e2 - math.f64_max_exp + 1.0)
+			e2 = math.f64_max_exp - 1.0
+		} else if e2 <= math.f64_min_exp {
+			y *= pow(2.0, e2 - math.f64_min_exp - 1.0)
+			e2 = math.f64_min_exp + 1.0
+		}
+		return y * pow(2.0, e2)
+	}
+}
+
+// frexp breaks f into a normalized fraction
+// and an integral power of two.
+// It returns frac and exp satisfying f == frac × 2**exp,
+// with the absolute value of frac in the interval [½, 1).
+//
+// special cases are:
+// frexp(±0) = ±0, 0
+// frexp(±inf) = ±inf, 0
+// frexp(nan) = nan, 0
+// pub fn frexp(f f64) (f64, int) {
+// // special cases
+// if f == 0.0 {
+// return f, 0 // correctly return -0
+// }
+// if is_inf(f, 0) || is_nan(f) {
+// return f, 0
+// }
+// f_norm, mut exp := normalize(f)
+// mut x := f64_bits(f_norm)
+// exp += int((x>>shift)&mask) - bias + 1
+// x &= ~(mask << shift)
+// x |= (-1 + bias) << shift
+// return f64_from_bits(x), exp
+pub fn frexp(x f64) (f64, int) {
+	if x == 0.0 {
+		return 0.0, 0
+	} else if !is_finite(x) {
+		return x, 0
+	} else if abs(x) >= 0.5 && abs(x) < 1 { // Handle the common case
+		return x, 0
+	} else {
+		ex := ceil(log(abs(x)) / ln2)
+		mut ei := int(ex) // Prevent underflow and overflow of 2**(-ei)
+		if ei < int(math.f64_min_exp) {
+			ei = int(math.f64_min_exp)
+		}
+		if ei > -int(math.f64_min_exp) {
+			ei = -int(math.f64_min_exp)
+		}
+		mut f := x * pow(2.0, -ei)
+		if !is_finite(f) { // This should not happen
+			return f, 0
+		}
+		for abs(f) >= 1.0 {
+			ei++
+			f /= 2.0
+		}
+		for abs(f) > 0 && abs(f) < 0.5 {
+			ei--
+			f *= 2.0
+		}
+		return f, ei
+	}
+}
+
+// special cases are:
+// expm1(+inf) = +inf
+// expm1(-inf) = -1
+// expm1(nan) = nan
+pub fn expm1(x f64) f64 {
+	if is_inf(x, 1) || is_nan(x) {
+		return x
+	}
+	if is_inf(x, -1) {
+		return f64(-1)
+	}
+	// FIXME: this should be improved
+	if abs(x) < ln2 { // Compute the taylor series S = x + (1/2!) x^2 + (1/3!) x^3 + ...
+		mut i := 1.0
+		mut sum := x
+		mut term := x / 1.0
+		i++
+		term *= x / f64(i)
+		sum += term
+		for abs(term) > abs(sum) * internal.f64_epsilon {
+			i++
+			term *= x / f64(i)
+			sum += term
+		}
+		return sum
+	} else {
+		return exp(x) - 1
+	}
+}
+
+// exp1 returns e**r × 2**k where r = hi - lo and |r| ≤ ln(2)/2.
+fn expmulti(hi f64, lo f64, k int) f64 {
+	exp_p1 := 1.66666666666666657415e-01 // 0x3FC55555; 0x55555555
+	exp_p2 := -2.77777777770155933842e-03 // 0xBF66C16C; 0x16BEBD93
+	exp_p3 := 6.61375632143793436117e-05 // 0x3F11566A; 0xAF25DE2C
+	exp_p4 := -1.65339022054652515390e-06 // 0xBEBBBD41; 0xC5D26BF1
+	exp_p5 := 4.13813679705723846039e-08 // 0x3E663769; 0x72BEA4D0
+	r := hi - lo
+	t := r * r
+	c := r - t * (exp_p1 + t * (exp_p2 + t * (exp_p3 + t * (exp_p4 + t * exp_p5))))
+	y := 1 - ((lo - (r * c) / (2 - c)) - hi)
+	// TODO(rsc): make sure ldexp can handle boundary k
+	return ldexp(y, k)
+}