blob: 45c8ecb3a8cf25974abf58a9a31084030e59388c (
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
80
81
82
83
84
|
// 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
}
|
