1 files changed, 85 insertions, 0 deletions
diff --git a/src/math/cbrt.go b/src/math/cbrt.go
new file mode 100644
index 0000000..e5e9548
--- /dev/null
+++ b/src/math/cbrt.go
@@ -0,0 +1,85 @@
+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+// The go 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 x.
+//
+// Special cases are:
+//
+//	Cbrt(±0) = ±0
+//	Cbrt(±Inf) = ±Inf
+//	Cbrt(NaN) = NaN
+func Cbrt(x float64) float64 {
+	if haveArchCbrt {
+		return archCbrt(x)
+	}
+	return cbrt(x)
+}
+
+func cbrt(x float64) float64 {
+	const (
+		B1             = 715094163                   // (682-0.03306235651)*2**20
+		B2             = 696219795                   // (664-0.03306235651)*2**20
+		C              = 5.42857142857142815906e-01  // 19/35     = 0x3FE15F15F15F15F1
+		D              = -7.05306122448979611050e-01 // -864/1225 = 0xBFE691DE2532C834
+		E              = 1.41428571428571436819e+00  // 99/70     = 0x3FF6A0EA0EA0EA0F
+		F              = 1.60714285714285720630e+00  // 45/28     = 0x3FF9B6DB6DB6DB6E
+		G              = 3.57142857142857150787e-01  // 5/14      = 0x3FD6DB6DB6DB6DB7
+		SmallestNormal = 2.22507385850720138309e-308 // 2**-1022  = 0x0010000000000000
+	)
+	// special cases
+	switch {
+	case x == 0 || IsNaN(x) || IsInf(x, 0):
+		return x
+	}
+
+	sign := false
+	if x < 0 {
+		x = -x
+		sign = true
+	}
+
+	// rough cbrt to 5 bits
+	t := Float64frombits(Float64bits(x)/3 + B1<<32)
+	if x < SmallestNormal {
+		// subnormal number
+		t = float64(1 << 54) // set t= 2**54
+		t *= x
+		t = Float64frombits(Float64bits(t)/3 + B2<<32)
+	}
+
+	// new cbrt to 23 bits
+	r := t * t / x
+	s := C + r*t
+	t *= G + F/(s+E+D/s)
+
+	// chop to 22 bits, make larger than cbrt(x)
+	t = Float64frombits(Float64bits(t)&(0xFFFFFFFFC<<28) + 1<<30)
+
+	// one step newton 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 the sign bit
+	if sign {
+		t = -t
+	}
+	return t
+}