math: implement pow in pure V (#12105)

vlib/math/pow.v

@@ -34,3 +34,119 @@ pub fn pow10(n int) f64 {
 	// n < -323
 	return 0.0
 }
+
+// pow returns base raised to the provided power.
+//
+// todo(playXE): make this function work on JS backend, probably problem of JS codegen that it does not work.
+pub fn pow(x f64, y f64) f64 {
+	if y == 0 || x == 1 {
+		return 1
+	} else if y == 1 {
+		return x
+	} else if is_nan(x) || is_nan(y) {
+		return nan()
+	} else if x == 0 {
+		if y < 0 {
+			if is_odd_int(y) {
+				return copysign(inf(1), x)
+			}
+			return inf(1)
+		} else if y > 0 {
+			if is_odd_int(y) {
+				return x
+			}
+			return 0
+		}
+	} else if is_inf(y, 0) {
+		if x == -1 {
+			return 1
+		} else if (abs(x) < 1) == is_inf(y, 1) {
+			return 0
+		} else {
+			return inf(1)
+		}
+	} else if is_inf(x, 0) {
+		if is_inf(x, -1) {
+			return pow(1 / x, -y)
+		}
+
+		if y < 0 {
+			return 0
+		} else if y > 0 {
+			return inf(1)
+		}
+	} else if y == 0.5 {
+		return sqrt(x)
+	} else if y == -0.5 {
+		return 1 / sqrt(x)
+	}
+	mut yi, mut yf := modf(abs(y))
+
+	if yf != 0 && x < 0 {
+		return nan()
+	}
+	if yi >= (u64(1) << 63) {
+		// yi is a large even int that will lead to overflow (or underflow to 0)
+		// for all x except -1 (x == 1 was handled earlier)
+
+		if x == -1 {
+			return 1
+		} else if (abs(x) < 1) == (y > 0) {
+			return 0
+		} else {
+			return inf(1)
+		}
+	}
+
+	// ans = a1 * 2**ae (= 1 for now).
+	mut a1 := 1.0
+	mut ae := 0
+
+	// ans *= x**yf
+	if yf != 0 {
+		if yf > 0.5 {
+			yf--
+			yi++
+		}
+
+		a1 = exp(yf * log(x))
+	}
+
+	// ans *= x**yi
+	// by multiplying in successive squarings
+	// of x according to bits of yi.
+	// accumulate powers of two into exp.
+	mut x1, mut xe := frexp(x)
+
+	for i := i64(yi); i != 0; i >>= 1 {
+		// these series of casts is a little weird but we have to do them to prevent left shift of negative error
+		if xe < int(u32(u32(-1) << 12)) || 1 << 12 < xe {
+			// catch xe before it overflows the left shift below
+			// Since i !=0 it has at least one bit still set, so ae will accumulate xe
+			// on at least one more iteration, ae += xe is a lower bound on ae
+			// the lower bound on ae exceeds the size of a float64 exp
+			// so the final call to Ldexp will produce under/overflow (0/Inf)
+			ae += xe
+			break
+		}
+		if i & 1 == 1 {
+			a1 *= x1
+			ae += xe
+		}
+		x1 *= x1
+		xe <<= 1
+		if x1 < .5 {
+			x1 += x1
+			xe--
+		}
+	}
+
+	// ans = a1*2**ae
+	// if y < 0 { ans = 1 / ans }
+	// but in the opposite order
+	if y < 0 {
+		a1 = 1 / a1
+		ae = -ae
+	}
+	return ldexp(a1, ae)
+}