Add math.cbrt

Based on the https://www.netlib.org/fdlibm/s_cbrt.c

core/math/cbrt.odin

@@ -0,0 +1,103 @@
+package math
+
+import "base:intrinsics"
+
+@(require_results) cbrt_f16le :: proc "contextless" (x: f16le) -> f16le { return #force_inline f16le(cbrt_f16(f16(x))) }
+@(require_results) cbrt_f16be :: proc "contextless" (x: f16be) -> f16be { return #force_inline f16be(cbrt_f16(f16(x))) }
+@(require_results) cbrt_f32le :: proc "contextless" (x: f32le) -> f32le { return #force_inline f32le(cbrt_f32(f32(x))) }
+@(require_results) cbrt_f32be :: proc "contextless" (x: f32be) -> f32be { return #force_inline f32be(cbrt_f32(f32(x))) }
+@(require_results) cbrt_f64le :: proc "contextless" (x: f64le) -> f64le { return #force_inline f64le(cbrt_f64(f64(x))) }
+@(require_results) cbrt_f64be :: proc "contextless" (x: f64be) -> f64be { return #force_inline f64be(cbrt_f64(f64(x))) }
+
+// cbrt returns the cube root of x.
+//
+// Special cases are:
+//
+//	cbrt(±0) = ±0
+//	cbrt(±Inf) = ±Inf
+//	cbrt(NaN) = NaN
+cbrt :: proc{
+	cbrt_f16, cbrt_f16le, cbrt_f16be,
+	cbrt_f32, cbrt_f32le, cbrt_f32be,
+	cbrt_f64, cbrt_f64le, cbrt_f64be,
+}
+
+
+
+@(require_results)
+cbrt_f16 :: proc "contextless" (x: f16) -> f16 { return #force_inline f16(cbrt_f64(f64(x))) }
+@(require_results)
+cbrt_f32 :: proc "contextless" (x: f32) -> f32 { return #force_inline f32(cbrt_f64(f64(x))) }
+
+// cbrt returns the cube root of x.
+//
+// Special cases are:
+//
+//	cbrt(±0) = ±0
+//	cbrt(±Inf) = ±Inf
+//	cbrt(NaN) = NaN
+@(require_results)
+cbrt_f64 :: proc "contextless" (x: f64) -> f64 {
+	// 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.
+	// ====================================================
+
+	B1              :: 715094163          // (682-0.03306235651)*2**20
+	B2              :: 696219795          // (664-0.03306235651)*2**20
+	SMALLEST_NORMAL :: 0h0010000000000000 // 2.22507385850720138309e-308 == 2**-1022
+
+	C               :: 0h3FE15F15F15F15F1 // 5.42857142857142815906e-01  == 19/35
+	D               :: 0hBFE691DE2532C834 // -7.05306122448979611050e-01 == -864/1225
+	E               :: 0h3FF6A0EA0EA0EA0F // 1.41428571428571436819e+00  == 99/70
+	F               :: 0h3FF9B6DB6DB6DB6E // 1.60714285714285720630e+00  == 45/28
+	G               :: 0h3FD6DB6DB6DB6DB7 // 3.57142857142857150787e-01  == 5/14
+
+	x := x
+
+	switch {
+	case x == 0 || is_nan(x) || is_inf(x, 0):
+		return x
+	}
+
+	sign := false
+	if x < 0 {
+		x = -x
+		sign = true
+	}
+
+	// Approximate cbrt (5-bits)
+	t := transmute(f64)((transmute(u64)x)/3 + B1<<32)
+	if x < SMALLEST_NORMAL {
+		t = f64(1 << 54)
+		t *= x
+		t = transmute(f64)((transmute(u64)x)/3 + B2<<32)
+	}
+
+	// Approximate cbrt (23-bits)
+	r := t * t / x
+	s := C + r*t
+	t *= G + F/(s+E + D/s)
+
+	// Truncate to 22 bits, make larger than cbrt(x)
+	t = transmute(f64)((transmute(u64)t)&(0xffffffffc<<28) + 1<<30)
+
+	// Single step of Newton-Raphson iteration to 53 bits with error less than 0.667ulps
+	s = t * t // t*t is exact
+	r = x / s
+	w := t + t
+	r = (r - t) / (w + r) // r-s is exact
+	t = t + t*r
+
+	// Restore sign
+	if sign {
+		t = -t
+	}
+	return t
+}