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