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 original C code and the long comment below are
8 // from FreeBSD's /usr/src/lib/msun/src/e_sqrt.c and
9 // came with this notice. The go code is a simplified
10 // version of the original C.
11 //
12 // ====================================================
13 // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
14 //
15 // Developed at SunPro, a Sun Microsystems, Inc. business.
16 // Permission to use, copy, modify, and distribute this
17 // software is freely granted, provided that this notice
18 // is preserved.
19 // ====================================================
20 //
21 // __ieee754_sqrt(x)
22 // Return correctly rounded sqrt.
23 // -----------------------------------------
24 // | Use the hardware sqrt if you have one |
25 // -----------------------------------------
26 // Method:
27 // Bit by bit method using integer arithmetic. (Slow, but portable)
28 // 1. Normalization
29 // Scale x to y in [1,4) with even powers of 2:
30 // find an integer k such that 1 <= (y=x*2**(2k)) < 4, then
31 // sqrt(x) = 2**k * sqrt(y)
32 // 2. Bit by bit computation
33 // Let q = sqrt(y) truncated to i bit after binary point (q = 1),
34 // i 0
35 // i+1 2
36 // s = 2*q , and y = 2 * ( y - q ). (1)
37 // i i i i
38 //
39 // To compute q from q , one checks whether
40 // i+1 i
41 //
42 // -(i+1) 2
43 // (q + 2 ) <= y. (2)
44 // i
45 // -(i+1)
46 // If (2) is false, then q = q ; otherwise q = q + 2 .
47 // i+1 i i+1 i
48 //
49 // With some algebraic manipulation, it is not difficult to see
50 // that (2) is equivalent to
51 // -(i+1)
52 // s + 2 <= y (3)
53 // i i
54 //
55 // The advantage of (3) is that s and y can be computed by
56 // i i
57 // the following recurrence formula:
58 // if (3) is false
59 //
60 // s = s , y = y ; (4)
61 // i+1 i i+1 i
62 //
63 // otherwise,
64 // -i -(i+1)
65 // s = s + 2 , y = y - s - 2 (5)
66 // i+1 i i+1 i i
67 //
68 // One may easily use induction to prove (4) and (5).
69 // Note. Since the left hand side of (3) contain only i+2 bits,
70 // it is not necessary to do a full (53-bit) comparison
71 // in (3).
72 // 3. Final rounding
73 // After generating the 53 bits result, we compute one more bit.
74 // Together with the remainder, we can decide whether the
75 // result is exact, bigger than 1/2ulp, or less than 1/2ulp
76 // (it will never equal to 1/2ulp).
77 // The rounding mode can be detected by checking whether
78 // huge + tiny is equal to huge, and whether huge - tiny is
79 // equal to huge for some floating point number "huge" and "tiny".
80 //
81 //
82 // Notes: Rounding mode detection omitted. The constants "mask", "shift",
83 // and "bias" are found in src/math/bits.go
84 85 // Sqrt returns the square root of x.
86 //
87 // Special cases are:
88 //
89 // Sqrt(+Inf) = +Inf
90 // Sqrt(±0) = ±0
91 // Sqrt(x < 0) = NaN
92 // Sqrt(NaN) = NaN
93 func Sqrt(x float64) float64 {
94 return sqrt(x)
95 }
96 97 // Note: On systems where Sqrt is a single instruction, the compiler
98 // may turn a direct call into a direct use of that instruction instead.
99 100 func sqrt(x float64) float64 {
101 // special cases
102 switch {
103 case x == 0 || IsNaN(x) || IsInf(x, 1):
104 return x
105 case x < 0:
106 return NaN()
107 }
108 ix := Float64bits(x)
109 // normalize x
110 exp := int((ix >> shift) & mask)
111 if exp == 0 { // subnormal x
112 for ix&(1<<shift) == 0 {
113 ix <<= 1
114 exp--
115 }
116 exp++
117 }
118 exp -= bias // unbias exponent
119 ix &^= mask << shift
120 ix |= 1 << shift
121 if exp&1 == 1 { // odd exp, double x to make it even
122 ix <<= 1
123 }
124 exp >>= 1 // exp = exp/2, exponent of square root
125 // generate sqrt(x) bit by bit
126 ix <<= 1
127 var q, s uint64 // q = sqrt(x)
128 r := uint64(1 << (shift + 1)) // r = moving bit from MSB to LSB
129 for r != 0 {
130 t := s + r
131 if t <= ix {
132 s = t + r
133 ix -= t
134 q += r
135 }
136 ix <<= 1
137 r >>= 1
138 }
139 // final rounding
140 if ix != 0 { // remainder, result not exact
141 q += q & 1 // round according to extra bit
142 }
143 ix = q>>1 + uint64(exp-1+bias)<<shift // significand + biased exponent
144 return Float64frombits(ix)
145 }
146