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