math: cbrt fix (#14395)

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

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])