From a3fa575cf61c4aef694d80bf16f7f2b526b504c6 Mon Sep 17 00:00:00 2001
From: Leo Developer <LeoDeveloper@protonmail.com>
Date: Mon, 4 Sep 2023 17:37:29 +0200
Subject: [PATCH] math: speedup the pure V math.pow implementation for
 non-fractional powers (#19270)

---
 vlib/math/ROADMAP.md  |  4 ++--
 vlib/math/math_test.v |  8 +++++++-
 vlib/math/pow.v       | 15 +++++++++++++++
 3 files changed, 24 insertions(+), 3 deletions(-)

diff --git a/vlib/math/ROADMAP.md b/vlib/math/ROADMAP.md
index 8564244179d21a..4668bb7106402f 100644
--- a/vlib/math/ROADMAP.md
+++ b/vlib/math/ROADMAP.md
@@ -1,4 +1,4 @@
 - [x] Move `vsl/vmath` to `vlib/math` as default backend
-- [ ] Implement `log` in pure V
-- [ ] Implement `pow` in pure V
+- [x] Implement `log` in pure V
+- [x] Implement `pow` in pure V
 - [ ] Define functions for initial release of hardware implementations
diff --git a/vlib/math/math_test.v b/vlib/math/math_test.v
index 9e681f7ed301ae..84c47a402af9e8 100644
--- a/vlib/math/math_test.v
+++ b/vlib/math/math_test.v
@@ -776,7 +776,8 @@ fn test_pow() {
 		[pi, nan()], [inf(1), -pi], [inf(1), -0.0], [inf(1), 0],
 		[inf(1), 1], [inf(1), pi], [inf(1), nan()], [nan(), -pi],
 		[nan(), -0.0], [nan(), 0], [nan(), 1], [nan(), pi], [nan(),
-			nan()]]
+			nan()],
+		[5.0, 2.0], [5.0, 3.0], [5.0, 10.0], [5.0, -2.0], [-5.0, -2.0]]
 	pow_sc_ := [f64(0), // pow(-inf, -pi)
 	 	-0.0, // pow(-inf, -3)
 	 	1, // pow(-inf, -0)
@@ -839,6 +840,11 @@ fn test_pow() {
 	 	nan(), // pow(nan, 1)
 	 	nan(), // pow(nan, pi)
 	 	nan(), // pow(nan, nan)
+	 	25, // pow(5, 2) => 5 * 5
+	 	125, // pow(5, 3) => 5 * 5 * 5
+	 	9765625, // pow(5, 10)
+	 	0.04, // pow(5, -2)
+	 	-0.04, // pow(-5, -2)
 	]
 	for i := 0; i < vfpow_sc_.len; i++ {
 		f := pow(vfpow_sc_[i][0], vfpow_sc_[i][1])
diff --git a/vlib/math/pow.v b/vlib/math/pow.v
index 007ebdd4a91bcc..f045362f750dc9 100644
--- a/vlib/math/pow.v
+++ b/vlib/math/pow.v
@@ -80,6 +80,10 @@ pub fn pow(x f64, y f64) f64 {
 		return x
 	} else if is_nan(x) || is_nan(y) {
 		return nan()
+	} else if y == 2 {
+		return x * x
+	} else if y == 3 {
+		return x * x * x
 	} else if x == 0 {
 		if y < 0 {
 			if is_odd_int(y) {
@@ -133,6 +137,17 @@ pub fn pow(x f64, y f64) f64 {
 		}
 	}
 
+	if yf == 0.0 {
+		mut result := x
+		for _ in 1 .. i64(yi) {
+			result *= x
+		}
+		if y > 0 {
+			return result
+		}
+		return copysign(1, x) / abs(result)
+	}
+
 	// ans = a1 * 2**ae (= 1 for now).
 	mut a1 := 1.0
 	mut ae := 0