blob: d2b7e910b84868aaa807fc2dc5e0916e195f8447 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
|
// 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 its argument.
//
// 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
)
// TODO(rsc): Remove manual inlining of IsNaN, IsInf
// when compiler does it for us
// special cases
switch {
case x == 0 || x != x || x < -MaxFloat64 || x > MaxFloat64: // x == 0 || IsNaN(x) || IsInf(x, 0):
return x
}
sign := false
if x < 0 {
x = -x
sign = true
}
// Reduce argument
f, e := Frexp(x)
m := e % 3
if m > 0 {
m -= 3
e -= m // e is multiple of 3
}
f = Ldexp(f, m) // 0.125 <= f < 1.0
// Estimate cube root
switch m {
case 0: // 0.5 <= f < 1.0
f = A1*f + A2 - A3/(A4+f)
case -1: // 0.25 <= f < 0.5
f = B1*f + B2 - B3/(B4+f)
default: // 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
}
|
