forked from chewxy/math32
-
Notifications
You must be signed in to change notification settings - Fork 0
/
sqrt.go
53 lines (49 loc) · 1.03 KB
/
sqrt.go
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
package math32
func Sqrt(x float32) float32
// TODO: add assembly for !build noasm
func sqrt(x float32) float32 {
// special cases
switch {
case x == 0 || IsNaN(x) || IsInf(x, 1):
return x
case x < 0:
return NaN()
}
ix := Float32bits(x)
// normalize x
exp := int((ix >> shift) & mask)
if exp == 0 { // subnormal x
for ix&(1<<shift) == 0 {
ix <<= 1
exp--
}
exp++
}
exp -= bias // unbias exponent
ix &^= mask << shift
ix |= 1 << shift
if exp&1 == 1 { // odd exp, double x to make it even
ix <<= 1
}
exp >>= 1 // exp = exp/2, exponent of square root
// generate sqrt(x) bit by bit
ix <<= 1
var q, s uint32 // q = sqrt(x)
r := uint32(1 << (shift + 1)) // r = moving bit from MSB to LSB
for r != 0 {
t := s + r
if t <= ix {
s = t + r
ix -= t
q += r
}
ix <<= 1
r >>= 1
}
// final rounding
if ix != 0 { // remainder, result not exact
q += q & 1 // round according to extra bit
}
ix = q>>1 + uint32(exp-1+bias)<<shift // significand + biased exponent
return Float32frombits(ix)
}