1 /* Copyright (C) 2007-2023 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 #define BID_128RES
25 #include "bid_internal.h"
26 #include "bid_sqrt_macros.h"
27 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
28 #include <fenv.h>
29
30 #define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT
31 #endif
32
33 BID128_FUNCTION_ARG1 (bid128_sqrt, x)
34
35 UINT256 M256, C256, C4, C8;
36 UINT128 CX, CX1, CX2, A10, S2, T128, TP128, CS, CSM, res;
37 UINT64 sign_x, Carry;
38 SINT64 D;
39 int_float fx, f64;
40 int exponent_x, bin_expon_cx;
41 int digits, scale, exponent_q;
42 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
43 fexcept_t binaryflags = 0;
44 #endif
45
46 // unpack arguments, check for NaN or Infinity
47 if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
48 res.w[1] = CX.w[1];
49 res.w[0] = CX.w[0];
50 // NaN ?
51 if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
52 #ifdef SET_STATUS_FLAGS
53 if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN
54 __set_status_flags (pfpsf, INVALID_EXCEPTION);
55 #endif
56 res.w[1] = CX.w[1] & QUIET_MASK64;
57 BID_RETURN (res);
58 }
59 // x is Infinity?
60 if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
61 res.w[1] = CX.w[1];
62 if (sign_x) {
63 // -Inf, return NaN
64 res.w[1] = 0x7c00000000000000ull;
65 #ifdef SET_STATUS_FLAGS
66 __set_status_flags (pfpsf, INVALID_EXCEPTION);
67 #endif
68 }
69 BID_RETURN (res);
70 }
71 // x is 0 otherwise
72
73 res.w[1] =
74 sign_x |
75 ((((UINT64) (exponent_x + DECIMAL_EXPONENT_BIAS_128)) >> 1) << 49);
76 res.w[0] = 0;
77 BID_RETURN (res);
78 }
79 if (sign_x) {
80 res.w[1] = 0x7c00000000000000ull;
81 res.w[0] = 0;
82 #ifdef SET_STATUS_FLAGS
83 __set_status_flags (pfpsf, INVALID_EXCEPTION);
84 #endif
85 BID_RETURN (res);
86 }
87 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
88 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
89 #endif
90 // 2^64
91 f64.i = 0x5f800000;
92
93 // fx ~ CX
94 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
95 bin_expon_cx = ((fx.i >> 23) & 0xff) - 0x7f;
96 digits = estimate_decimal_digits[bin_expon_cx];
97
98 A10 = CX;
99 if (exponent_x & 1) {
100 A10.w[1] = (CX.w[1] << 3) | (CX.w[0] >> 61);
101 A10.w[0] = CX.w[0] << 3;
102 CX2.w[1] = (CX.w[1] << 1) | (CX.w[0] >> 63);
103 CX2.w[0] = CX.w[0] << 1;
104 __add_128_128 (A10, A10, CX2);
105 }
106
107 CS.w[0] = short_sqrt128 (A10);
108 CS.w[1] = 0;
109 // check for exact result
110 if (CS.w[0] * CS.w[0] == A10.w[0]) {
111 __mul_64x64_to_128_fast (S2, CS.w[0], CS.w[0]);
112 if (S2.w[1] == A10.w[1]) // && S2.w[0]==A10.w[0])
113 {
114 get_BID128_very_fast (&res, 0,
115 (exponent_x +
116 DECIMAL_EXPONENT_BIAS_128) >> 1, CS);
117 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
118 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
119 #endif
120 BID_RETURN (res);
121 }
122 }
123 // get number of digits in CX
124 D = CX.w[1] - power10_index_binexp_128[bin_expon_cx].w[1];
125 if (D > 0
126 || (!D && CX.w[0] >= power10_index_binexp_128[bin_expon_cx].w[0]))
127 digits++;
128
129 // if exponent is odd, scale coefficient by 10
130 scale = 67 - digits;
131 exponent_q = exponent_x - scale;
132 scale += (exponent_q & 1); // exp. bias is even
133
134 if (scale > 38) {
135 T128 = power10_table_128[scale - 37];
136 __mul_128x128_low (CX1, CX, T128);
137
138 TP128 = power10_table_128[37];
139 __mul_128x128_to_256 (C256, CX1, TP128);
140 } else {
141 T128 = power10_table_128[scale];
142 __mul_128x128_to_256 (C256, CX, T128);
143 }
144
145
146 // 4*C256
147 C4.w[3] = (C256.w[3] << 2) | (C256.w[2] >> 62);
148 C4.w[2] = (C256.w[2] << 2) | (C256.w[1] >> 62);
149 C4.w[1] = (C256.w[1] << 2) | (C256.w[0] >> 62);
150 C4.w[0] = C256.w[0] << 2;
151
152 long_sqrt128 (&CS, C256);
153
154 #ifndef IEEE_ROUND_NEAREST
155 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
156 if (!((rnd_mode) & 3)) {
157 #endif
158 #endif
159 // compare to midpoints
160 CSM.w[1] = (CS.w[1] << 1) | (CS.w[0] >> 63);
161 CSM.w[0] = (CS.w[0] + CS.w[0]) | 1;
162 // CSM^2
163 //__mul_128x128_to_256(M256, CSM, CSM);
164 __sqr128_to_256 (M256, CSM);
165
166 if (C4.w[3] > M256.w[3]
167 || (C4.w[3] == M256.w[3]
168 && (C4.w[2] > M256.w[2]
169 || (C4.w[2] == M256.w[2]
170 && (C4.w[1] > M256.w[1]
171 || (C4.w[1] == M256.w[1]
172 && C4.w[0] > M256.w[0])))))) {
173 // round up
174 CS.w[0]++;
175 if (!CS.w[0])
176 CS.w[1]++;
177 } else {
178 C8.w[1] = (CS.w[1] << 3) | (CS.w[0] >> 61);
179 C8.w[0] = CS.w[0] << 3;
180 // M256 - 8*CSM
181 __sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]);
182 __sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);
183 __sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);
184 M256.w[3] = M256.w[3] - Carry;
185
186 // if CSM' > C256, round up
187 if (M256.w[3] > C4.w[3]
188 || (M256.w[3] == C4.w[3]
189 && (M256.w[2] > C4.w[2]
190 || (M256.w[2] == C4.w[2]
191 && (M256.w[1] > C4.w[1]
192 || (M256.w[1] == C4.w[1]
193 && M256.w[0] > C4.w[0])))))) {
194 // round down
195 if (!CS.w[0])
196 CS.w[1]--;
197 CS.w[0]--;
198 }
199 }
200 #ifndef IEEE_ROUND_NEAREST
201 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
202 } else {
203 __sqr128_to_256 (M256, CS);
204 C8.w[1] = (CS.w[1] << 1) | (CS.w[0] >> 63);
205 C8.w[0] = CS.w[0] << 1;
206 if (M256.w[3] > C256.w[3]
207 || (M256.w[3] == C256.w[3]
208 && (M256.w[2] > C256.w[2]
209 || (M256.w[2] == C256.w[2]
210 && (M256.w[1] > C256.w[1]
211 || (M256.w[1] == C256.w[1]
212 && M256.w[0] > C256.w[0])))))) {
213 __sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]);
214 __sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);
215 __sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);
216 M256.w[3] = M256.w[3] - Carry;
217 M256.w[0]++;
218 if (!M256.w[0]) {
219 M256.w[1]++;
220 if (!M256.w[1]) {
221 M256.w[2]++;
222 if (!M256.w[2])
223 M256.w[3]++;
224 }
225 }
226
227 if (!CS.w[0])
228 CS.w[1]--;
229 CS.w[0]--;
230
231 if (M256.w[3] > C256.w[3]
232 || (M256.w[3] == C256.w[3]
233 && (M256.w[2] > C256.w[2]
234 || (M256.w[2] == C256.w[2]
235 && (M256.w[1] > C256.w[1]
236 || (M256.w[1] == C256.w[1]
237 && M256.w[0] > C256.w[0])))))) {
238
239 if (!CS.w[0])
240 CS.w[1]--;
241 CS.w[0]--;
242 }
243 }
244
245 else {
246 __add_carry_out (M256.w[0], Carry, M256.w[0], C8.w[0]);
247 __add_carry_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);
248 __add_carry_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);
249 M256.w[3] = M256.w[3] + Carry;
250 M256.w[0]++;
251 if (!M256.w[0]) {
252 M256.w[1]++;
253 if (!M256.w[1]) {
254 M256.w[2]++;
255 if (!M256.w[2])
256 M256.w[3]++;
257 }
258 }
259 if (M256.w[3] < C256.w[3]
260 || (M256.w[3] == C256.w[3]
261 && (M256.w[2] < C256.w[2]
262 || (M256.w[2] == C256.w[2]
263 && (M256.w[1] < C256.w[1]
264 || (M256.w[1] == C256.w[1]
265 && M256.w[0] <= C256.w[0])))))) {
266
267 CS.w[0]++;
268 if (!CS.w[0])
269 CS.w[1]++;
270 }
271 }
272 // RU?
273 if ((rnd_mode) == ROUNDING_UP) {
274 CS.w[0]++;
275 if (!CS.w[0])
276 CS.w[1]++;
277 }
278
279 }
280 #endif
281 #endif
282
283 #ifdef SET_STATUS_FLAGS
284 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
285 #endif
286 get_BID128_fast (&res, 0,
287 (exponent_q + DECIMAL_EXPONENT_BIAS_128) >> 1, CS);
288 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
289 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
290 #endif
291 BID_RETURN (res);
292 }
293
294
295
296 BID128_FUNCTION_ARGTYPE1 (bid128d_sqrt, UINT64, x)
297
298 UINT256 M256, C256, C4, C8;
299 UINT128 CX, CX1, CX2, A10, S2, T128, TP128, CS, CSM, res;
300 UINT64 sign_x, Carry;
301 SINT64 D;
302 int_float fx, f64;
303 int exponent_x, bin_expon_cx;
304 int digits, scale, exponent_q;
305 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
306 fexcept_t binaryflags = 0;
307 #endif
308
309 // unpack arguments, check for NaN or Infinity
310 // unpack arguments, check for NaN or Infinity
311 CX.w[1] = 0;
312 if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], x)) {
313 res.w[1] = CX.w[0];
314 res.w[0] = 0;
315 // NaN ?
316 if ((x & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
317 #ifdef SET_STATUS_FLAGS
318 if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN
319 __set_status_flags (pfpsf, INVALID_EXCEPTION);
320 #endif
321 res.w[0] = (CX.w[0] & 0x0003ffffffffffffull);
322 __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
323 res.w[1] |= ((CX.w[0]) & 0xfc00000000000000ull);
324 BID_RETURN (res);
325 }
326 // x is Infinity?
327 if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {
328 if (sign_x) {
329 // -Inf, return NaN
330 res.w[1] = 0x7c00000000000000ull;
331 #ifdef SET_STATUS_FLAGS
332 __set_status_flags (pfpsf, INVALID_EXCEPTION);
333 #endif
334 }
335 BID_RETURN (res);
336 }
337 // x is 0 otherwise
338
339 exponent_x =
340 exponent_x - DECIMAL_EXPONENT_BIAS + DECIMAL_EXPONENT_BIAS_128;
341 res.w[1] =
342 sign_x | ((((UINT64) (exponent_x + DECIMAL_EXPONENT_BIAS_128)) >> 1)
343 << 49);
344 res.w[0] = 0;
345 BID_RETURN (res);
346 }
347 if (sign_x) {
348 res.w[1] = 0x7c00000000000000ull;
349 res.w[0] = 0;
350 #ifdef SET_STATUS_FLAGS
351 __set_status_flags (pfpsf, INVALID_EXCEPTION);
352 #endif
353 BID_RETURN (res);
354 }
355 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
356 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
357 #endif
358 exponent_x =
359 exponent_x - DECIMAL_EXPONENT_BIAS + DECIMAL_EXPONENT_BIAS_128;
360
361 // 2^64
362 f64.i = 0x5f800000;
363
364 // fx ~ CX
365 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
366 bin_expon_cx = ((fx.i >> 23) & 0xff) - 0x7f;
367 digits = estimate_decimal_digits[bin_expon_cx];
368
369 A10 = CX;
370 if (exponent_x & 1) {
371 A10.w[1] = (CX.w[1] << 3) | (CX.w[0] >> 61);
372 A10.w[0] = CX.w[0] << 3;
373 CX2.w[1] = (CX.w[1] << 1) | (CX.w[0] >> 63);
374 CX2.w[0] = CX.w[0] << 1;
375 __add_128_128 (A10, A10, CX2);
376 }
377
378 CS.w[0] = short_sqrt128 (A10);
379 CS.w[1] = 0;
380 // check for exact result
381 if (CS.w[0] * CS.w[0] == A10.w[0]) {
382 __mul_64x64_to_128_fast (S2, CS.w[0], CS.w[0]);
383 if (S2.w[1] == A10.w[1]) {
384 get_BID128_very_fast (&res, 0,
385 (exponent_x + DECIMAL_EXPONENT_BIAS_128) >> 1,
386 CS);
387 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
388 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
389 #endif
390 BID_RETURN (res);
391 }
392 }
393 // get number of digits in CX
394 D = CX.w[1] - power10_index_binexp_128[bin_expon_cx].w[1];
395 if (D > 0
396 || (!D && CX.w[0] >= power10_index_binexp_128[bin_expon_cx].w[0]))
397 digits++;
398
399 // if exponent is odd, scale coefficient by 10
400 scale = 67 - digits;
401 exponent_q = exponent_x - scale;
402 scale += (exponent_q & 1); // exp. bias is even
403
404 if (scale > 38) {
405 T128 = power10_table_128[scale - 37];
406 __mul_128x128_low (CX1, CX, T128);
407
408 TP128 = power10_table_128[37];
409 __mul_128x128_to_256 (C256, CX1, TP128);
410 } else {
411 T128 = power10_table_128[scale];
412 __mul_128x128_to_256 (C256, CX, T128);
413 }
414
415
416 // 4*C256
417 C4.w[3] = (C256.w[3] << 2) | (C256.w[2] >> 62);
418 C4.w[2] = (C256.w[2] << 2) | (C256.w[1] >> 62);
419 C4.w[1] = (C256.w[1] << 2) | (C256.w[0] >> 62);
420 C4.w[0] = C256.w[0] << 2;
421
422 long_sqrt128 (&CS, C256);
423
424 #ifndef IEEE_ROUND_NEAREST
425 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
426 if (!((rnd_mode) & 3)) {
427 #endif
428 #endif
429 // compare to midpoints
430 CSM.w[1] = (CS.w[1] << 1) | (CS.w[0] >> 63);
431 CSM.w[0] = (CS.w[0] + CS.w[0]) | 1;
432 // CSM^2
433 //__mul_128x128_to_256(M256, CSM, CSM);
434 __sqr128_to_256 (M256, CSM);
435
436 if (C4.w[3] > M256.w[3]
437 || (C4.w[3] == M256.w[3]
438 && (C4.w[2] > M256.w[2]
439 || (C4.w[2] == M256.w[2]
440 && (C4.w[1] > M256.w[1]
441 || (C4.w[1] == M256.w[1]
442 && C4.w[0] > M256.w[0])))))) {
443 // round up
444 CS.w[0]++;
445 if (!CS.w[0])
446 CS.w[1]++;
447 } else {
448 C8.w[1] = (CS.w[1] << 3) | (CS.w[0] >> 61);
449 C8.w[0] = CS.w[0] << 3;
450 // M256 - 8*CSM
451 __sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]);
452 __sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);
453 __sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);
454 M256.w[3] = M256.w[3] - Carry;
455
456 // if CSM' > C256, round up
457 if (M256.w[3] > C4.w[3]
458 || (M256.w[3] == C4.w[3]
459 && (M256.w[2] > C4.w[2]
460 || (M256.w[2] == C4.w[2]
461 && (M256.w[1] > C4.w[1]
462 || (M256.w[1] == C4.w[1]
463 && M256.w[0] > C4.w[0])))))) {
464 // round down
465 if (!CS.w[0])
466 CS.w[1]--;
467 CS.w[0]--;
468 }
469 }
470 #ifndef IEEE_ROUND_NEAREST
471 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
472 } else {
473 __sqr128_to_256 (M256, CS);
474 C8.w[1] = (CS.w[1] << 1) | (CS.w[0] >> 63);
475 C8.w[0] = CS.w[0] << 1;
476 if (M256.w[3] > C256.w[3]
477 || (M256.w[3] == C256.w[3]
478 && (M256.w[2] > C256.w[2]
479 || (M256.w[2] == C256.w[2]
480 && (M256.w[1] > C256.w[1]
481 || (M256.w[1] == C256.w[1]
482 && M256.w[0] > C256.w[0])))))) {
483 __sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]);
484 __sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);
485 __sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);
486 M256.w[3] = M256.w[3] - Carry;
487 M256.w[0]++;
488 if (!M256.w[0]) {
489 M256.w[1]++;
490 if (!M256.w[1]) {
491 M256.w[2]++;
492 if (!M256.w[2])
493 M256.w[3]++;
494 }
495 }
496
497 if (!CS.w[0])
498 CS.w[1]--;
499 CS.w[0]--;
500
501 if (M256.w[3] > C256.w[3]
502 || (M256.w[3] == C256.w[3]
503 && (M256.w[2] > C256.w[2]
504 || (M256.w[2] == C256.w[2]
505 && (M256.w[1] > C256.w[1]
506 || (M256.w[1] == C256.w[1]
507 && M256.w[0] > C256.w[0])))))) {
508
509 if (!CS.w[0])
510 CS.w[1]--;
511 CS.w[0]--;
512 }
513 }
514
515 else {
516 __add_carry_out (M256.w[0], Carry, M256.w[0], C8.w[0]);
517 __add_carry_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);
518 __add_carry_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);
519 M256.w[3] = M256.w[3] + Carry;
520 M256.w[0]++;
521 if (!M256.w[0]) {
522 M256.w[1]++;
523 if (!M256.w[1]) {
524 M256.w[2]++;
525 if (!M256.w[2])
526 M256.w[3]++;
527 }
528 }
529 if (M256.w[3] < C256.w[3]
530 || (M256.w[3] == C256.w[3]
531 && (M256.w[2] < C256.w[2]
532 || (M256.w[2] == C256.w[2]
533 && (M256.w[1] < C256.w[1]
534 || (M256.w[1] == C256.w[1]
535 && M256.w[0] <= C256.w[0])))))) {
536
537 CS.w[0]++;
538 if (!CS.w[0])
539 CS.w[1]++;
540 }
541 }
542 // RU?
543 if ((rnd_mode) == ROUNDING_UP) {
544 CS.w[0]++;
545 if (!CS.w[0])
546 CS.w[1]++;
547 }
548
549 }
550 #endif
551 #endif
552
553 #ifdef SET_STATUS_FLAGS
554 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
555 #endif
556 get_BID128_fast (&res, 0, (exponent_q + DECIMAL_EXPONENT_BIAS_128) >> 1,
557 CS);
558 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
559 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
560 #endif
561 BID_RETURN (res);
562
563
564 }