1 // Copyright 2009 The Go Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
4 5 package math
6 7 // The go code is a modified version of the original C code from
8 // http://www.netlib.org/fdlibm/s_cbrt.c and came with this notice.
9 //
10 // ====================================================
11 // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
12 //
13 // Developed at SunSoft, a Sun Microsystems, Inc. business.
14 // Permission to use, copy, modify, and distribute this
15 // software is freely granted, provided that this notice
16 // is preserved.
17 // ====================================================
18 19 // Cbrt returns the cube root of x.
20 //
21 // Special cases are:
22 //
23 // Cbrt(±0) = ±0
24 // Cbrt(±Inf) = ±Inf
25 // Cbrt(NaN) = NaN
26 func Cbrt(x float64) float64 {
27 if haveArchCbrt {
28 return archCbrt(x)
29 }
30 return cbrt(x)
31 }
32 33 func cbrt(x float64) float64 {
34 const (
35 B1 = 715094163 // (682-0.03306235651)*2**20
36 B2 = 696219795 // (664-0.03306235651)*2**20
37 C = 5.42857142857142815906e-01 // 19/35 = 0x3FE15F15F15F15F1
38 D = -7.05306122448979611050e-01 // -864/1225 = 0xBFE691DE2532C834
39 E = 1.41428571428571436819e+00 // 99/70 = 0x3FF6A0EA0EA0EA0F
40 F = 1.60714285714285720630e+00 // 45/28 = 0x3FF9B6DB6DB6DB6E
41 G = 3.57142857142857150787e-01 // 5/14 = 0x3FD6DB6DB6DB6DB7
42 SmallestNormal = 2.22507385850720138309e-308 // 2**-1022 = 0x0010000000000000
43 )
44 // special cases
45 switch {
46 case x == 0 || IsNaN(x) || IsInf(x, 0):
47 return x
48 }
49 50 sign := false
51 if x < 0 {
52 x = -x
53 sign = true
54 }
55 56 // rough cbrt to 5 bits
57 t := Float64frombits(Float64bits(x)/3 + B1<<32)
58 if x < SmallestNormal {
59 // subnormal number
60 t = float64(1 << 54) // set t= 2**54
61 t *= x
62 t = Float64frombits(Float64bits(t)/3 + B2<<32)
63 }
64 65 // new cbrt to 23 bits
66 r := t * t / x
67 s := C + r*t
68 t *= G + F/(s+E+D/s)
69 70 // chop to 22 bits, make larger than cbrt(x)
71 t = Float64frombits(Float64bits(t)&(0xFFFFFFFFC<<28) + 1<<30)
72 73 // one step newton iteration to 53 bits with error less than 0.667ulps
74 s = t * t // t*t is exact
75 r = x / s
76 w := t + t
77 r = (r - t) / (w + r) // r-s is exact
78 t = t + t*r
79 80 // restore the sign bit
81 if sign {
82 t = -t
83 }
84 return t
85 }
86