forked from OSchip/llvm-project
77 lines
1.6 KiB
Go
77 lines
1.6 KiB
Go
// 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 algorithm is based in part on "Optimal Partitioning of
|
|
Newton's Method for Calculating Roots", by Gunter Meinardus
|
|
and G. D. Taylor, Mathematics of Computation © 1980 American
|
|
Mathematical Society.
|
|
(http://www.jstor.org/stable/2006387?seq=9, accessed 11-Feb-2010)
|
|
*/
|
|
|
|
// Cbrt returns the cube root of x.
|
|
//
|
|
// Special cases are:
|
|
// Cbrt(±0) = ±0
|
|
// Cbrt(±Inf) = ±Inf
|
|
// Cbrt(NaN) = NaN
|
|
func Cbrt(x float64) float64 {
|
|
const (
|
|
A1 = 1.662848358e-01
|
|
A2 = 1.096040958e+00
|
|
A3 = 4.105032829e-01
|
|
A4 = 5.649335816e-01
|
|
B1 = 2.639607233e-01
|
|
B2 = 8.699282849e-01
|
|
B3 = 1.629083358e-01
|
|
B4 = 2.824667908e-01
|
|
C1 = 4.190115298e-01
|
|
C2 = 6.904625373e-01
|
|
C3 = 6.46502159e-02
|
|
C4 = 1.412333954e-01
|
|
)
|
|
// special cases
|
|
switch {
|
|
case x == 0 || IsNaN(x) || IsInf(x, 0):
|
|
return x
|
|
}
|
|
sign := false
|
|
if x < 0 {
|
|
x = -x
|
|
sign = true
|
|
}
|
|
// Reduce argument and estimate cube root
|
|
f, e := Frexp(x) // 0.5 <= f < 1.0
|
|
m := e % 3
|
|
if m > 0 {
|
|
m -= 3
|
|
e -= m // e is multiple of 3
|
|
}
|
|
switch m {
|
|
case 0: // 0.5 <= f < 1.0
|
|
f = A1*f + A2 - A3/(A4+f)
|
|
case -1:
|
|
f *= 0.5 // 0.25 <= f < 0.5
|
|
f = B1*f + B2 - B3/(B4+f)
|
|
default: // m == -2
|
|
f *= 0.25 // 0.125 <= f < 0.25
|
|
f = C1*f + C2 - C3/(C4+f)
|
|
}
|
|
y := Ldexp(f, e/3) // e/3 = exponent of cube root
|
|
|
|
// Iterate
|
|
s := y * y * y
|
|
t := s + x
|
|
y *= (t + x) / (s + t)
|
|
// Reiterate
|
|
s = (y*y*y - x) / x
|
|
y -= y * (((14.0/81.0)*s-(2.0/9.0))*s + (1.0 / 3.0)) * s
|
|
if sign {
|
|
y = -y
|
|
}
|
|
return y
|
|
}
|