|
0
|
1 /* Copyright (C) 2007, 2009 Free Software Foundation, Inc.
|
|
|
2
|
|
|
3 This file is part of GCC.
|
|
|
4
|
|
|
5 GCC is free software; you can redistribute it and/or modify it under
|
|
|
6 the terms of the GNU General Public License as published by the Free
|
|
|
7 Software Foundation; either version 3, or (at your option) any later
|
|
|
8 version.
|
|
|
9
|
|
|
10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
|
11 WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
|
12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
|
13 for more details.
|
|
|
14
|
|
|
15 Under Section 7 of GPL version 3, you are granted additional
|
|
|
16 permissions described in the GCC Runtime Library Exception, version
|
|
|
17 3.1, as published by the Free Software Foundation.
|
|
|
18
|
|
|
19 You should have received a copy of the GNU General Public License and
|
|
|
20 a copy of the GCC Runtime Library Exception along with this program;
|
|
|
21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
|
|
|
22 <http://www.gnu.org/licenses/>. */
|
|
|
23
|
|
|
24 #ifndef _SQRT_MACROS_H_
|
|
|
25 #define _SQRT_MACROS_H_
|
|
|
26
|
|
|
27 #define FENCE __fence
|
|
|
28
|
|
|
29 #if DOUBLE_EXTENDED_ON
|
|
|
30
|
|
|
31 extern BINARY80 SQRT80 (BINARY80);
|
|
|
32
|
|
|
33
|
|
|
34 __BID_INLINE__ UINT64
|
|
|
35 short_sqrt128 (UINT128 A10) {
|
|
|
36 BINARY80 lx, ly, l64;
|
|
|
37 int_float f64;
|
|
|
38
|
|
|
39 // 2^64
|
|
|
40 f64.i = 0x5f800000;
|
|
|
41 l64 = (BINARY80) f64.d;
|
|
|
42 lx = (BINARY80) A10.w[1] * l64 + (BINARY80) A10.w[0];
|
|
|
43 ly = SQRT80 (lx);
|
|
|
44 return (UINT64) ly;
|
|
|
45 }
|
|
|
46
|
|
|
47
|
|
|
48 __BID_INLINE__ void
|
|
|
49 long_sqrt128 (UINT128 * pCS, UINT256 C256) {
|
|
|
50 UINT256 C4;
|
|
|
51 UINT128 CS;
|
|
|
52 UINT64 X;
|
|
|
53 SINT64 SE;
|
|
|
54 BINARY80 l64, lm64, l128, lxL, lx, ly, lS, lSH, lSL, lE, l3, l2,
|
|
|
55 l1, l0, lp, lCl;
|
|
|
56 int_float fx, f64, fm64;
|
|
|
57 int *ple = (int *) &lx;
|
|
|
58
|
|
|
59 // 2^64
|
|
|
60 f64.i = 0x5f800000;
|
|
|
61 l64 = (BINARY80) f64.d;
|
|
|
62
|
|
|
63 l128 = l64 * l64;
|
|
|
64 lx = l3 = (BINARY80) C256.w[3] * l64 * l128;
|
|
|
65 l2 = (BINARY80) C256.w[2] * l128;
|
|
|
66 lx = FENCE (lx + l2);
|
|
|
67 l1 = (BINARY80) C256.w[1] * l64;
|
|
|
68 lx = FENCE (lx + l1);
|
|
|
69 l0 = (BINARY80) C256.w[0];
|
|
|
70 lx = FENCE (lx + l0);
|
|
|
71 // sqrt(C256)
|
|
|
72 lS = SQRT80 (lx);
|
|
|
73
|
|
|
74 // get coefficient
|
|
|
75 // 2^(-64)
|
|
|
76 fm64.i = 0x1f800000;
|
|
|
77 lm64 = (BINARY80) fm64.d;
|
|
|
78 CS.w[1] = (UINT64) (lS * lm64);
|
|
|
79 CS.w[0] = (UINT64) (lS - (BINARY80) CS.w[1] * l64);
|
|
|
80
|
|
|
81 ///////////////////////////////////////
|
|
|
82 // CAUTION!
|
|
|
83 // little endian code only
|
|
|
84 // add solution for big endian
|
|
|
85 //////////////////////////////////////
|
|
|
86 lSH = lS;
|
|
|
87 *((UINT64 *) & lSH) &= 0xffffffff00000000ull;
|
|
|
88
|
|
|
89 // correction for C256 rounding
|
|
|
90 lCl = FENCE (l3 - lx);
|
|
|
91 lCl = FENCE (lCl + l2);
|
|
|
92 lCl = FENCE (lCl + l1);
|
|
|
93 lCl = FENCE (lCl + l0);
|
|
|
94
|
|
|
95 lSL = lS - lSH;
|
|
|
96
|
|
|
97 //////////////////////////////////////////
|
|
|
98 // Watch for compiler re-ordering
|
|
|
99 //
|
|
|
100 /////////////////////////////////////////
|
|
|
101 // C256-S^2
|
|
|
102 lxL = FENCE (lx - lSH * lSH);
|
|
|
103 lp = lSH * lSL;
|
|
|
104 lp += lp;
|
|
|
105 lxL = FENCE (lxL - lp);
|
|
|
106 lSL *= lSL;
|
|
|
107 lxL = FENCE (lxL - lSL);
|
|
|
108 lCl += lxL;
|
|
|
109
|
|
|
110 // correction term
|
|
|
111 lE = lCl / (lS + lS);
|
|
|
112
|
|
|
113 // get low part of coefficient
|
|
|
114 X = CS.w[0];
|
|
|
115 if (lCl >= 0) {
|
|
|
116 SE = (SINT64) (lE);
|
|
|
117 CS.w[0] += SE;
|
|
|
118 if (CS.w[0] < X)
|
|
|
119 CS.w[1]++;
|
|
|
120 } else {
|
|
|
121 SE = (SINT64) (-lE);
|
|
|
122 CS.w[0] -= SE;
|
|
|
123 if (CS.w[0] > X)
|
|
|
124 CS.w[1]--;
|
|
|
125 }
|
|
|
126
|
|
|
127 pCS->w[0] = CS.w[0];
|
|
|
128 pCS->w[1] = CS.w[1];
|
|
|
129 }
|
|
|
130
|
|
|
131 #else
|
|
|
132
|
|
|
133 extern double sqrt (double);
|
|
|
134
|
|
|
135 __BID_INLINE__ UINT64
|
|
|
136 short_sqrt128 (UINT128 A10) {
|
|
|
137 UINT256 ARS, ARS0, AE0, AE, S;
|
|
|
138
|
|
|
139 UINT64 MY, ES, CY;
|
|
|
140 double lx, l64;
|
|
|
141 int_double f64, ly;
|
|
|
142 int ey, k;
|
|
|
143
|
|
|
144 // 2^64
|
|
|
145 f64.i = 0x43f0000000000000ull;
|
|
|
146 l64 = f64.d;
|
|
|
147 lx = (double) A10.w[1] * l64 + (double) A10.w[0];
|
|
|
148 ly.d = 1.0 / sqrt (lx);
|
|
|
149
|
|
|
150 MY = (ly.i & 0x000fffffffffffffull) | 0x0010000000000000ull;
|
|
|
151 ey = 0x3ff - (ly.i >> 52);
|
|
|
152
|
|
|
153 // A10*RS^2
|
|
|
154 __mul_64x128_to_192 (ARS0, MY, A10);
|
|
|
155 __mul_64x192_to_256 (ARS, MY, ARS0);
|
|
|
156
|
|
|
157 // shr by 2*ey+40, to get a 64-bit value
|
|
|
158 k = (ey << 1) + 104 - 64;
|
|
|
159 if (k >= 128) {
|
|
|
160 if (k > 128)
|
|
|
161 ES = (ARS.w[2] >> (k - 128)) | (ARS.w[3] << (192 - k));
|
|
|
162 else
|
|
|
163 ES = ARS.w[2];
|
|
|
164 } else {
|
|
|
165 if (k >= 64) {
|
|
|
166 ARS.w[0] = ARS.w[1];
|
|
|
167 ARS.w[1] = ARS.w[2];
|
|
|
168 k -= 64;
|
|
|
169 }
|
|
|
170 if (k) {
|
|
|
171 __shr_128 (ARS, ARS, k);
|
|
|
172 }
|
|
|
173 ES = ARS.w[0];
|
|
|
174 }
|
|
|
175
|
|
|
176 ES = ((SINT64) ES) >> 1;
|
|
|
177
|
|
|
178 if (((SINT64) ES) < 0) {
|
|
|
179 ES = -ES;
|
|
|
180
|
|
|
181 // A*RS*eps (scaled by 2^64)
|
|
|
182 __mul_64x192_to_256 (AE0, ES, ARS0);
|
|
|
183
|
|
|
184 AE.w[0] = AE0.w[1];
|
|
|
185 AE.w[1] = AE0.w[2];
|
|
|
186 AE.w[2] = AE0.w[3];
|
|
|
187
|
|
|
188 __add_carry_out (S.w[0], CY, ARS0.w[0], AE.w[0]);
|
|
|
189 __add_carry_in_out (S.w[1], CY, ARS0.w[1], AE.w[1], CY);
|
|
|
190 S.w[2] = ARS0.w[2] + AE.w[2] + CY;
|
|
|
191 } else {
|
|
|
192 // A*RS*eps (scaled by 2^64)
|
|
|
193 __mul_64x192_to_256 (AE0, ES, ARS0);
|
|
|
194
|
|
|
195 AE.w[0] = AE0.w[1];
|
|
|
196 AE.w[1] = AE0.w[2];
|
|
|
197 AE.w[2] = AE0.w[3];
|
|
|
198
|
|
|
199 __sub_borrow_out (S.w[0], CY, ARS0.w[0], AE.w[0]);
|
|
|
200 __sub_borrow_in_out (S.w[1], CY, ARS0.w[1], AE.w[1], CY);
|
|
|
201 S.w[2] = ARS0.w[2] - AE.w[2] - CY;
|
|
|
202 }
|
|
|
203
|
|
|
204 k = ey + 51;
|
|
|
205
|
|
|
206 if (k >= 64) {
|
|
|
207 if (k >= 128) {
|
|
|
208 S.w[0] = S.w[2];
|
|
|
209 S.w[1] = 0;
|
|
|
210 k -= 128;
|
|
|
211 } else {
|
|
|
212 S.w[0] = S.w[1];
|
|
|
213 S.w[1] = S.w[2];
|
|
|
214 }
|
|
|
215 k -= 64;
|
|
|
216 }
|
|
|
217 if (k) {
|
|
|
218 __shr_128 (S, S, k);
|
|
|
219 }
|
|
|
220
|
|
|
221
|
|
|
222 return (UINT64) ((S.w[0] + 1) >> 1);
|
|
|
223
|
|
|
224 }
|
|
|
225
|
|
|
226
|
|
|
227
|
|
|
228 __BID_INLINE__ void
|
|
|
229 long_sqrt128 (UINT128 * pCS, UINT256 C256) {
|
|
|
230 UINT512 ARS0, ARS;
|
|
|
231 UINT256 ARS00, AE, AE2, S;
|
|
|
232 UINT128 ES, ES2, ARS1;
|
|
|
233 UINT64 ES32, CY, MY;
|
|
|
234 double l64, l128, lx, l2, l1, l0;
|
|
|
235 int_double f64, ly;
|
|
|
236 int ey, k, k2;
|
|
|
237
|
|
|
238 // 2^64
|
|
|
239 f64.i = 0x43f0000000000000ull;
|
|
|
240 l64 = f64.d;
|
|
|
241
|
|
|
242 l128 = l64 * l64;
|
|
|
243 lx = (double) C256.w[3] * l64 * l128;
|
|
|
244 l2 = (double) C256.w[2] * l128;
|
|
|
245 lx = FENCE (lx + l2);
|
|
|
246 l1 = (double) C256.w[1] * l64;
|
|
|
247 lx = FENCE (lx + l1);
|
|
|
248 l0 = (double) C256.w[0];
|
|
|
249 lx = FENCE (lx + l0);
|
|
|
250 // sqrt(C256)
|
|
|
251 ly.d = 1.0 / sqrt (lx);
|
|
|
252
|
|
|
253 MY = (ly.i & 0x000fffffffffffffull) | 0x0010000000000000ull;
|
|
|
254 ey = 0x3ff - (ly.i >> 52);
|
|
|
255
|
|
|
256 // A10*RS^2, scaled by 2^(2*ey+104)
|
|
|
257 __mul_64x256_to_320 (ARS0, MY, C256);
|
|
|
258 __mul_64x320_to_384 (ARS, MY, ARS0);
|
|
|
259
|
|
|
260 // shr by k=(2*ey+104)-128
|
|
|
261 // expect k is in the range (192, 256) if result in [10^33, 10^34)
|
|
|
262 // apply an additional signed shift by 1 at the same time (to get eps=eps0/2)
|
|
|
263 k = (ey << 1) + 104 - 128 - 192;
|
|
|
264 k2 = 64 - k;
|
|
|
265 ES.w[0] = (ARS.w[3] >> (k + 1)) | (ARS.w[4] << (k2 - 1));
|
|
|
266 ES.w[1] = (ARS.w[4] >> k) | (ARS.w[5] << k2);
|
|
|
267 ES.w[1] = ((SINT64) ES.w[1]) >> 1;
|
|
|
268
|
|
|
269 // A*RS >> 192 (for error term computation)
|
|
|
270 ARS1.w[0] = ARS0.w[3];
|
|
|
271 ARS1.w[1] = ARS0.w[4];
|
|
|
272
|
|
|
273 // A*RS>>64
|
|
|
274 ARS00.w[0] = ARS0.w[1];
|
|
|
275 ARS00.w[1] = ARS0.w[2];
|
|
|
276 ARS00.w[2] = ARS0.w[3];
|
|
|
277 ARS00.w[3] = ARS0.w[4];
|
|
|
278
|
|
|
279 if (((SINT64) ES.w[1]) < 0) {
|
|
|
280 ES.w[0] = -ES.w[0];
|
|
|
281 ES.w[1] = -ES.w[1];
|
|
|
282 if (ES.w[0])
|
|
|
283 ES.w[1]--;
|
|
|
284
|
|
|
285 // A*RS*eps
|
|
|
286 __mul_128x128_to_256 (AE, ES, ARS1);
|
|
|
287
|
|
|
288 __add_carry_out (S.w[0], CY, ARS00.w[0], AE.w[0]);
|
|
|
289 __add_carry_in_out (S.w[1], CY, ARS00.w[1], AE.w[1], CY);
|
|
|
290 __add_carry_in_out (S.w[2], CY, ARS00.w[2], AE.w[2], CY);
|
|
|
291 S.w[3] = ARS00.w[3] + AE.w[3] + CY;
|
|
|
292 } else {
|
|
|
293 // A*RS*eps
|
|
|
294 __mul_128x128_to_256 (AE, ES, ARS1);
|
|
|
295
|
|
|
296 __sub_borrow_out (S.w[0], CY, ARS00.w[0], AE.w[0]);
|
|
|
297 __sub_borrow_in_out (S.w[1], CY, ARS00.w[1], AE.w[1], CY);
|
|
|
298 __sub_borrow_in_out (S.w[2], CY, ARS00.w[2], AE.w[2], CY);
|
|
|
299 S.w[3] = ARS00.w[3] - AE.w[3] - CY;
|
|
|
300 }
|
|
|
301
|
|
|
302 // 3/2*eps^2, scaled by 2^128
|
|
|
303 ES32 = ES.w[1] + (ES.w[1] >> 1);
|
|
|
304 __mul_64x64_to_128 (ES2, ES32, ES.w[1]);
|
|
|
305 // A*RS*3/2*eps^2
|
|
|
306 __mul_128x128_to_256 (AE2, ES2, ARS1);
|
|
|
307
|
|
|
308 // result, scaled by 2^(ey+52-64)
|
|
|
309 __add_carry_out (S.w[0], CY, S.w[0], AE2.w[0]);
|
|
|
310 __add_carry_in_out (S.w[1], CY, S.w[1], AE2.w[1], CY);
|
|
|
311 __add_carry_in_out (S.w[2], CY, S.w[2], AE2.w[2], CY);
|
|
|
312 S.w[3] = S.w[3] + AE2.w[3] + CY;
|
|
|
313
|
|
|
314 // k in (0, 64)
|
|
|
315 k = ey + 51 - 128;
|
|
|
316 k2 = 64 - k;
|
|
|
317 S.w[0] = (S.w[1] >> k) | (S.w[2] << k2);
|
|
|
318 S.w[1] = (S.w[2] >> k) | (S.w[3] << k2);
|
|
|
319
|
|
|
320 // round to nearest
|
|
|
321 S.w[0]++;
|
|
|
322 if (!S.w[0])
|
|
|
323 S.w[1]++;
|
|
|
324
|
|
|
325 pCS->w[0] = (S.w[1] << 63) | (S.w[0] >> 1);
|
|
|
326 pCS->w[1] = S.w[1] >> 1;
|
|
|
327
|
|
|
328 }
|
|
|
329
|
|
|
330 #endif
|
|
|
331 #endif
|
