From 8d141878ce06e39a1d8a5ef988c7c3502ce84511 Mon Sep 17 00:00:00 2001
From: David 'Epper' Marshall <epper.marshall@gmail.com>
Date: Sat, 14 May 2022 14:06:38 -0400
Subject: [PATCH] math: cbrt fix (#14395)

---
 vlib/math/cbrt.v      | 16 ++++++++++++++--
 vlib/math/math_test.v |  7 +++++++
 2 files changed, 21 insertions(+), 2 deletions(-)

diff --git a/vlib/math/cbrt.v b/vlib/math/cbrt.v
index 2c34ef2c3..a9d883ee2 100644
--- a/vlib/math/cbrt.v
+++ b/vlib/math/cbrt.v
@@ -1,5 +1,17 @@
 module math
 
+// The vlang code is a modified version of the original C code from
+// http://www.netlib.org/fdlibm/s_cbrt.c and came with this notice.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunSoft, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+
 // cbrt returns the cube root of a.
 //
 // special cases are:
@@ -25,12 +37,12 @@ pub fn cbrt(a f64) f64 {
 		sign = true
 	}
 	// rough cbrt to 5 bits
-	mut t := f64_from_bits(f64_bits(x) / u64(3 + (u64(b1) << 32)))
+	mut t := f64_from_bits(f64_bits(x) / u64(3) + (u64(b1) << 32))
 	if x < smallest_normal {
 		// subnormal number
 		t = f64(u64(1) << 54) // set t= 2**54
 		t *= x
-		t = f64_from_bits(f64_bits(t) / u64(3 + (u64(b2) << 32)))
+		t = f64_from_bits(f64_bits(t) / u64(3) + (u64(b2) << 32))
 	}
 	// new cbrt to 23 bits
 	mut r := t * t / x
diff --git a/vlib/math/math_test.v b/vlib/math/math_test.v
index 652cd4c7f..aeaabf0a1 100644
--- a/vlib/math/math_test.v
+++ b/vlib/math/math_test.v
@@ -506,6 +506,13 @@ fn test_mod() {
 	assert (0.6447968302508578) == f
 }
 
+fn test_cbrt() {
+	cbrts := [2.0, 10, 56]
+	for idx, i in [8.0, 1000, 175_616] {
+		assert cbrt(i) == cbrts[idx]
+	}
+}
+
 fn test_exp() {
 	for i := 0; i < math.vf_.len; i++ {
 		f := exp(math.vf_[i])
-- 
2.30.2