]> git.ipfire.org Git - thirdparty/gcc.git/blob - libgcc/config/libbid/bid64_div.c
projects / thirdparty / gcc.git / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
Update copyright years.
[thirdparty/gcc.git] / libgcc / config / libbid / bid64_div.c
1 /* Copyright (C) 2007-2020 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 divide
26 *****************************************************************************
27 *
28 * Algorithm description:
29 *
30 * if(coefficient_x<coefficient_y)
31 * p = number_digits(coefficient_y) - number_digits(coefficient_x)
32 * A = coefficient_x*10^p
33 * B = coefficient_y
34 * CA= A*10^(15+j), j=0 for A>=B, 1 otherwise
35 * Q = 0
36 * else
37 * get Q=(int)(coefficient_x/coefficient_y)
38 * (based on double precision divide)
39 * check for exact divide case
40 * Let R = coefficient_x - Q*coefficient_y
41 * Let m=16-number_digits(Q)
42 * CA=R*10^m, Q=Q*10^m
43 * B = coefficient_y
44 * endif
45 * if (CA<2^64)
46 * Q += CA/B (64-bit unsigned divide)
47 * else
48 * get final Q using double precision divide, followed by 3 integer
49 * iterations
50 * if exact result, eliminate trailing zeros
51 * check for underflow
52 * round coefficient to nearest
53 *
54 ****************************************************************************/
55
56 #include "bid_internal.h"
57 #include "bid_div_macros.h"
58 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
59 #include <fenv.h>
60
61 #define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT
62 #endif
63
64 extern UINT32 convert_table[5][128][2];
65 extern SINT8 factors[][2];
66 extern UINT8 packed_10000_zeros[];
67
68
69 #if DECIMAL_CALL_BY_REFERENCE
70
71 void
72 bid64_div (UINT64 * pres, UINT64 * px,
73 UINT64 *
74 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
75 _EXC_INFO_PARAM) {
76 UINT64 x, y;
77 #else
78
79 UINT64
80 bid64_div (UINT64 x,
81 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
82 _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
83 #endif
84 UINT128 CA, CT;
85 UINT64 sign_x, sign_y, coefficient_x, coefficient_y, A, B, QX, PD;
86 UINT64 A2, Q, Q2, B2, B4, B5, R, T, DU, res;
87 UINT64 valid_x, valid_y;
88 SINT64 D;
89 int_double t_scale, tempq, temp_b;
90 int_float tempx, tempy;
91 double da, db, dq, da_h, da_l;
92 int exponent_x, exponent_y, bin_expon_cx;
93 int diff_expon, ed1, ed2, bin_index;
94 int rmode, amount;
95 int nzeros, i, j, k, d5;
96 UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
97 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
98 fexcept_t binaryflags = 0;
99 #endif
100
101 #if DECIMAL_CALL_BY_REFERENCE
102 #if !DECIMAL_GLOBAL_ROUNDING
103 _IDEC_round rnd_mode = *prnd_mode;
104 #endif
105 x = *px;
106 y = *py;
107 #endif
108
109 valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
110 valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
111
112 // unpack arguments, check for NaN or Infinity
113 if (!valid_x) {
114 // x is Inf. or NaN
115 #ifdef SET_STATUS_FLAGS
116 if ((y & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
117 __set_status_flags (pfpsf, INVALID_EXCEPTION);
118 #endif
119
120 // test if x is NaN
121 if ((x & NAN_MASK64) == NAN_MASK64) {
122 #ifdef SET_STATUS_FLAGS
123 if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN
124 __set_status_flags (pfpsf, INVALID_EXCEPTION);
125 #endif
126 BID_RETURN (coefficient_x & QUIET_MASK64);
127 }
128 // x is Infinity?
129 if ((x & INFINITY_MASK64) == INFINITY_MASK64) {
130 // check if y is Inf or NaN
131 if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
132 // y==Inf, return NaN
133 if ((y & NAN_MASK64) == INFINITY_MASK64) { // Inf/Inf
134 #ifdef SET_STATUS_FLAGS
135 __set_status_flags (pfpsf, INVALID_EXCEPTION);
136 #endif
137 BID_RETURN (NAN_MASK64);
138 }
139 } else {
140 // otherwise return +/-Inf
141 BID_RETURN (((x ^ y) & 0x8000000000000000ull) |
142 INFINITY_MASK64);
143 }
144 }
145 // x==0
146 if (((y & INFINITY_MASK64) != INFINITY_MASK64)
147 && !(coefficient_y)) {
148 // y==0 , return NaN
149 #ifdef SET_STATUS_FLAGS
150 __set_status_flags (pfpsf, INVALID_EXCEPTION);
151 #endif
152 BID_RETURN (NAN_MASK64);
153 }
154 if (((y & INFINITY_MASK64) != INFINITY_MASK64)) {
155 if ((y & SPECIAL_ENCODING_MASK64) == SPECIAL_ENCODING_MASK64)
156 exponent_y = ((UINT32) (y >> 51)) & 0x3ff;
157 else
158 exponent_y = ((UINT32) (y >> 53)) & 0x3ff;
159 sign_y = y & 0x8000000000000000ull;
160
161 exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS;
162 if (exponent_x > DECIMAL_MAX_EXPON_64)
163 exponent_x = DECIMAL_MAX_EXPON_64;
164 else if (exponent_x < 0)
165 exponent_x = 0;
166 BID_RETURN ((sign_x ^ sign_y) | (((UINT64) exponent_x) << 53));
167 }
168
169 }
170 if (!valid_y) {
171 // y is Inf. or NaN
172
173 // test if y is NaN
174 if ((y & NAN_MASK64) == NAN_MASK64) {
175 #ifdef SET_STATUS_FLAGS
176 if ((y & SNAN_MASK64) == SNAN_MASK64) // sNaN
177 __set_status_flags (pfpsf, INVALID_EXCEPTION);
178 #endif
179 BID_RETURN (coefficient_y & QUIET_MASK64);
180 }
181 // y is Infinity?
182 if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
183 // return +/-0
184 BID_RETURN (((x ^ y) & 0x8000000000000000ull));
185 }
186 // y is 0
187 #ifdef SET_STATUS_FLAGS
188 __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
189 #endif
190 BID_RETURN ((sign_x ^ sign_y) | INFINITY_MASK64);
191 }
192 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
193 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
194 #endif
195 diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS;
196
197 if (coefficient_x < coefficient_y) {
198 // get number of decimal digits for c_x, c_y
199
200 //--- get number of bits in the coefficients of x and y ---
201 tempx.d = (float) coefficient_x;
202 tempy.d = (float) coefficient_y;
203 bin_index = (tempy.i - tempx.i) >> 23;
204
205 A = coefficient_x * power10_index_binexp[bin_index];
206 B = coefficient_y;
207
208 temp_b.d = (double) B;
209
210 // compare A, B
211 DU = (A - B) >> 63;
212 ed1 = 15 + (int) DU;
213 ed2 = estimate_decimal_digits[bin_index] + ed1;
214 T = power10_table_128[ed1].w[0];
215 __mul_64x64_to_128 (CA, A, T);
216
217 Q = 0;
218 diff_expon = diff_expon - ed2;
219
220 // adjust double precision db, to ensure that later A/B - (int)(da/db) > -1
221 if (coefficient_y < 0x0020000000000000ull) {
222 temp_b.i += 1;
223 db = temp_b.d;
224 } else
225 db = (double) (B + 2 + (B & 1));
226
227 } else {
228 // get c_x/c_y
229
230 // set last bit before conversion to DP
231 A2 = coefficient_x | 1;
232 da = (double) A2;
233
234 db = (double) coefficient_y;
235
236 tempq.d = da / db;
237 Q = (UINT64) tempq.d;
238
239 R = coefficient_x - coefficient_y * Q;
240
241 // will use to get number of dec. digits of Q
242 bin_expon_cx = (tempq.i >> 52) - 0x3ff;
243
244 // R<0 ?
245 D = ((SINT64) R) >> 63;
246 Q += D;
247 R += (coefficient_y & D);
248
249 // exact result ?
250 if (((SINT64) R) <= 0) {
251 // can have R==-1 for coeff_y==1
252 res =
253 get_BID64 (sign_x ^ sign_y, diff_expon, (Q + R), rnd_mode,
254 pfpsf);
255 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
256 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
257 #endif
258 BID_RETURN (res);
259 }
260 // get decimal digits of Q
261 DU = power10_index_binexp[bin_expon_cx] - Q - 1;
262 DU >>= 63;
263
264 ed2 = 16 - estimate_decimal_digits[bin_expon_cx] - (int) DU;
265
266 T = power10_table_128[ed2].w[0];
267 __mul_64x64_to_128 (CA, R, T);
268 B = coefficient_y;
269
270 Q *= power10_table_128[ed2].w[0];
271 diff_expon -= ed2;
272
273 }
274
275 if (!CA.w[1]) {
276 Q2 = CA.w[0] / B;
277 B2 = B + B;
278 B4 = B2 + B2;
279 R = CA.w[0] - Q2 * B;
280 Q += Q2;
281 } else {
282
283 // 2^64
284 t_scale.i = 0x43f0000000000000ull;
285 // convert CA to DP
286 da_h = CA.w[1];
287 da_l = CA.w[0];
288 da = da_h * t_scale.d + da_l;
289
290 // quotient
291 dq = da / db;
292 Q2 = (UINT64) dq;
293
294 // get w[0] remainder
295 R = CA.w[0] - Q2 * B;
296
297 // R<0 ?
298 D = ((SINT64) R) >> 63;
299 Q2 += D;
300 R += (B & D);
301
302 // now R<6*B
303
304 // quick divide
305
306 // 4*B
307 B2 = B + B;
308 B4 = B2 + B2;
309
310 R = R - B4;
311 // R<0 ?
312 D = ((SINT64) R) >> 63;
313 // restore R if negative
314 R += (B4 & D);
315 Q2 += ((~D) & 4);
316
317 R = R - B2;
318 // R<0 ?
319 D = ((SINT64) R) >> 63;
320 // restore R if negative
321 R += (B2 & D);
322 Q2 += ((~D) & 2);
323
324 R = R - B;
325 // R<0 ?
326 D = ((SINT64) R) >> 63;
327 // restore R if negative
328 R += (B & D);
329 Q2 += ((~D) & 1);
330
331 Q += Q2;
332 }
333
334 #ifdef SET_STATUS_FLAGS
335 if (R) {
336 // set status flags
337 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
338 }
339 #ifndef LEAVE_TRAILING_ZEROS
340 else
341 #endif
342 #else
343 #ifndef LEAVE_TRAILING_ZEROS
344 if (!R)
345 #endif
346 #endif
347 #ifndef LEAVE_TRAILING_ZEROS
348 {
349 // eliminate trailing zeros
350
351 // check whether CX, CY are short
352 if ((coefficient_x <= 1024) && (coefficient_y <= 1024)) {
353 i = (int) coefficient_y - 1;
354 j = (int) coefficient_x - 1;
355 // difference in powers of 2 factors for Y and X
356 nzeros = ed2 - factors[i][0] + factors[j][0];
357 // difference in powers of 5 factors
358 d5 = ed2 - factors[i][1] + factors[j][1];
359 if (d5 < nzeros)
360 nzeros = d5;
361
362 __mul_64x64_to_128 (CT, Q, reciprocals10_64[nzeros]);
363
364 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
365 amount = short_recip_scale[nzeros];
366 Q = CT.w[1] >> amount;
367
368 diff_expon += nzeros;
369 } else {
370 tdigit[0] = Q & 0x3ffffff;
371 tdigit[1] = 0;
372 QX = Q >> 26;
373 QX32 = QX;
374 nzeros = 0;
375
376 for (j = 0; QX32; j++, QX32 >>= 7) {
377 k = (QX32 & 127);
378 tdigit[0] += convert_table[j][k][0];
379 tdigit[1] += convert_table[j][k][1];
380 if (tdigit[0] >= 100000000) {
381 tdigit[0] -= 100000000;
382 tdigit[1]++;
383 }
384 }
385
386 digit = tdigit[0];
387 if (!digit && !tdigit[1])
388 nzeros += 16;
389 else {
390 if (!digit) {
391 nzeros += 8;
392 digit = tdigit[1];
393 }
394 // decompose digit
395 PD = (UINT64) digit *0x068DB8BBull;
396 digit_h = (UINT32) (PD >> 40);
397 digit_low = digit - digit_h * 10000;
398
399 if (!digit_low)
400 nzeros += 4;
401 else
402 digit_h = digit_low;
403
404 if (!(digit_h & 1))
405 nzeros +=
406 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
407 (digit_h & 7));
408 }
409
410 if (nzeros) {
411 __mul_64x64_to_128 (CT, Q, reciprocals10_64[nzeros]);
412
413 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
414 amount = short_recip_scale[nzeros];
415 Q = CT.w[1] >> amount;
416 }
417 diff_expon += nzeros;
418
419 }
420 if (diff_expon >= 0) {
421 res =
422 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, Q,
423 rnd_mode, pfpsf);
424 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
425 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
426 #endif
427 BID_RETURN (res);
428 }
429 }
430 #endif
431
432 if (diff_expon >= 0) {
433 #ifdef IEEE_ROUND_NEAREST
434 // round to nearest code
435 // R*10
436 R += R;
437 R = (R << 2) + R;
438 B5 = B4 + B;
439
440 // compare 10*R to 5*B
441 R = B5 - R;
442 // correction for (R==0 && (Q&1))
443 R -= (Q & 1);
444 // R<0 ?
445 D = ((UINT64) R) >> 63;
446 Q += D;
447 #else
448 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
449 // round to nearest code
450 // R*10
451 R += R;
452 R = (R << 2) + R;
453 B5 = B4 + B;
454
455 // compare 10*R to 5*B
456 R = B5 - R;
457 // correction for (R==0 && (Q&1))
458 R -= (Q & 1);
459 // R<0 ?
460 D = ((UINT64) R) >> 63;
461 Q += D;
462 #else
463 rmode = rnd_mode;
464 if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
465 rmode = 3 - rmode;
466 switch (rmode) {
467 case 0: // round to nearest code
468 case ROUNDING_TIES_AWAY:
469 // R*10
470 R += R;
471 R = (R << 2) + R;
472 B5 = B4 + B;
473 // compare 10*R to 5*B
474 R = B5 - R;
475 // correction for (R==0 && (Q&1))
476 R -= ((Q | (rmode >> 2)) & 1);
477 // R<0 ?
478 D = ((UINT64) R) >> 63;
479 Q += D;
480 break;
481 case ROUNDING_DOWN:
482 case ROUNDING_TO_ZERO:
483 break;
484 default: // rounding up
485 Q++;
486 break;
487 }
488 #endif
489 #endif
490
491 res =
492 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, Q, rnd_mode,
493 pfpsf);
494 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
495 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
496 #endif
497 BID_RETURN (res);
498 } else {
499 // UF occurs
500
501 #ifdef SET_STATUS_FLAGS
502 if ((diff_expon + 16 < 0)) {
503 // set status flags
504 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
505 }
506 #endif
507 rmode = rnd_mode;
508 res =
509 get_BID64_UF (sign_x ^ sign_y, diff_expon, Q, R, rmode, pfpsf);
510 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
511 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
512 #endif
513 BID_RETURN (res);
514
515 }
516 }
517
518
519
520 TYPE0_FUNCTION_ARGTYPE1_ARG128 (UINT64, bid64dq_div, UINT64, x, y)
521 UINT256 CA4 =
522 { {0x0ull, 0x0ull, 0x0ull, 0x0ull} }, CA4r, P256, QB256;
523 UINT128 CX, CY, T128, CQ, CQ2, CR, CA, TP128, Qh, Tmp;
524 UINT64 sign_x, sign_y, T, carry64, D, Q_low, QX, valid_y, PD, res;
525 int_float fx, fy, f64;
526 UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
527 int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
528 digits_q, amount;
529 int nzeros, i, j, k, d5, done = 0;
530 unsigned rmode;
531 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
532 fexcept_t binaryflags = 0;
533 #endif
534
535 valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
536
537 // unpack arguments, check for NaN or Infinity
538 CX.w[1] = 0;
539 if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], (x))) {
540 #ifdef SET_STATUS_FLAGS
541 if (((y.w[1] & SNAN_MASK64) == SNAN_MASK64) || // y is sNaN
542 ((x & SNAN_MASK64) == SNAN_MASK64))
543 __set_status_flags (pfpsf, INVALID_EXCEPTION);
544 #endif
545 // test if x is NaN
546 if (((x) & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
547 res = CX.w[0];
548 BID_RETURN (res & QUIET_MASK64);
549 }
550 // x is Infinity?
551 if (((x) & 0x7800000000000000ull) == 0x7800000000000000ull) {
552 // check if y is Inf.
553 if (((y.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull))
554 // return NaN
555 {
556 #ifdef SET_STATUS_FLAGS
557 __set_status_flags (pfpsf, INVALID_EXCEPTION);
558 #endif
559 res = 0x7c00000000000000ull;
560 BID_RETURN (res);
561 }
562 if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
563 // otherwise return +/-Inf
564 res =
565 (((x) ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
566 BID_RETURN (res);
567 }
568 }
569 // x is 0
570 if ((y.w[1] & INFINITY_MASK64) != INFINITY_MASK64) {
571 if ((!CY.w[0]) && !(CY.w[1] & 0x0001ffffffffffffull)) {
572 #ifdef SET_STATUS_FLAGS
573 __set_status_flags (pfpsf, INVALID_EXCEPTION);
574 #endif
575 // x=y=0, return NaN
576 res = 0x7c00000000000000ull;
577 BID_RETURN (res);
578 }
579 // return 0
580 res = ((x) ^ y.w[1]) & 0x8000000000000000ull;
581 exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
582 if (exponent_x > DECIMAL_MAX_EXPON_64)
583 exponent_x = DECIMAL_MAX_EXPON_64;
584 else if (exponent_x < 0)
585 exponent_x = 0;
586 res |= (((UINT64) exponent_x) << 53);
587 BID_RETURN (res);
588 }
589 }
590 exponent_x += (DECIMAL_EXPONENT_BIAS_128 - DECIMAL_EXPONENT_BIAS);
591 if (!valid_y) {
592 // y is Inf. or NaN
593
594 // test if y is NaN
595 if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
596 #ifdef SET_STATUS_FLAGS
597 if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN
598 __set_status_flags (pfpsf, INVALID_EXCEPTION);
599 #endif
600 Tmp.w[1] = (CY.w[1] & 0x00003fffffffffffull);
601 Tmp.w[0] = CY.w[0];
602 TP128 = reciprocals10_128[18];
603 __mul_128x128_high (Qh, Tmp, TP128);
604 amount = recip_scale[18];
605 __shr_128 (Tmp, Qh, amount);
606 res = (CY.w[1] & 0xfc00000000000000ull) | Tmp.w[0];
607 BID_RETURN (res);
608 }
609 // y is Infinity?
610 if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
611 // return +/-0
612 res = sign_x ^ sign_y;
613 BID_RETURN (res);
614 }
615 // y is 0, return +/-Inf
616 res =
617 (((x) ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
618 #ifdef SET_STATUS_FLAGS
619 __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
620 #endif
621 BID_RETURN (res);
622 }
623 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
624 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
625 #endif
626 diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS;
627
628 if (__unsigned_compare_gt_128 (CY, CX)) {
629 // CX < CY
630
631 // 2^64
632 f64.i = 0x5f800000;
633
634 // fx ~ CX, fy ~ CY
635 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
636 fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
637 // expon_cy - expon_cx
638 bin_index = (fy.i - fx.i) >> 23;
639
640 if (CX.w[1]) {
641 T = power10_index_binexp_128[bin_index].w[0];
642 __mul_64x128_short (CA, T, CX);
643 } else {
644 T128 = power10_index_binexp_128[bin_index];
645 __mul_64x128_short (CA, CX.w[0], T128);
646 }
647
648 ed2 = 15;
649 if (__unsigned_compare_gt_128 (CY, CA))
650 ed2++;
651
652 T128 = power10_table_128[ed2];
653 __mul_128x128_to_256 (CA4, CA, T128);
654
655 ed2 += estimate_decimal_digits[bin_index];
656 CQ.w[0] = CQ.w[1] = 0;
657 diff_expon = diff_expon - ed2;
658
659 } else {
660 // get CQ = CX/CY
661 __div_128_by_128 (&CQ, &CR, CX, CY);
662
663 // get number of decimal digits in CQ
664 // 2^64
665 f64.i = 0x5f800000;
666 fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
667 // binary expon. of CQ
668 bin_expon = (fx.i - 0x3f800000) >> 23;
669
670 digits_q = estimate_decimal_digits[bin_expon];
671 TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
672 TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
673 if (__unsigned_compare_ge_128 (CQ, TP128))
674 digits_q++;
675
676 if (digits_q <= 16) {
677 if (!CR.w[1] && !CR.w[0]) {
678 res = get_BID64 (sign_x ^ sign_y, diff_expon,
679 CQ.w[0], rnd_mode, pfpsf);
680 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
681 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
682 #endif
683 BID_RETURN (res);
684 }
685
686 ed2 = 16 - digits_q;
687 T128.w[0] = power10_table_128[ed2].w[0];
688 __mul_64x128_to_192 (CA4, (T128.w[0]), CR);
689 diff_expon = diff_expon - ed2;
690 CQ.w[0] *= T128.w[0];
691 } else {
692 ed2 = digits_q - 16;
693 diff_expon += ed2;
694 T128 = reciprocals10_128[ed2];
695 __mul_128x128_to_256 (P256, CQ, T128);
696 amount = recip_scale[ed2];
697 CQ.w[0] = (P256.w[2] >> amount) | (P256.w[3] << (64 - amount));
698 CQ.w[1] = 0;
699
700 __mul_64x64_to_128 (CQ2, CQ.w[0], (power10_table_128[ed2].w[0]));
701
702 __mul_64x64_to_128 (QB256, CQ2.w[0], CY.w[0]);
703 QB256.w[1] += CQ2.w[0] * CY.w[1] + CQ2.w[1] * CY.w[0];
704
705 CA4.w[1] = CX.w[1] - QB256.w[1];
706 CA4.w[0] = CX.w[0] - QB256.w[0];
707 if (CX.w[0] < QB256.w[0])
708 CA4.w[1]--;
709 if (CR.w[0] || CR.w[1])
710 CA4.w[0] |= 1;
711 done = 1;
712
713 }
714
715 }
716 if (!done) {
717 __div_256_by_128 (&CQ, &CA4, CY);
718 }
719
720
721
722 #ifdef SET_STATUS_FLAGS
723 if (CA4.w[0] || CA4.w[1]) {
724 // set status flags
725 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
726 }
727 #ifndef LEAVE_TRAILING_ZEROS
728 else
729 #endif
730 #else
731 #ifndef LEAVE_TRAILING_ZEROS
732 if (!CA4.w[0] && !CA4.w[1])
733 #endif
734 #endif
735 #ifndef LEAVE_TRAILING_ZEROS
736 // check whether result is exact
737 {
738 // check whether CX, CY are short
739 if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
740 i = (int) CY.w[0] - 1;
741 j = (int) CX.w[0] - 1;
742 // difference in powers of 2 factors for Y and X
743 nzeros = ed2 - factors[i][0] + factors[j][0];
744 // difference in powers of 5 factors
745 d5 = ed2 - factors[i][1] + factors[j][1];
746 if (d5 < nzeros)
747 nzeros = d5;
748 // get P*(2^M[extra_digits])/10^extra_digits
749 __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]);
750
751 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
752 amount = recip_scale[nzeros];
753 __shr_128_long (CQ, Qh, amount);
754
755 diff_expon += nzeros;
756 } else {
757 // decompose Q as Qh*10^17 + Ql
758 Q_low = CQ.w[0];
759
760 {
761 tdigit[0] = Q_low & 0x3ffffff;
762 tdigit[1] = 0;
763 QX = Q_low >> 26;
764 QX32 = QX;
765 nzeros = 0;
766
767 for (j = 0; QX32; j++, QX32 >>= 7) {
768 k = (QX32 & 127);
769 tdigit[0] += convert_table[j][k][0];
770 tdigit[1] += convert_table[j][k][1];
771 if (tdigit[0] >= 100000000) {
772 tdigit[0] -= 100000000;
773 tdigit[1]++;
774 }
775 }
776
777 if (tdigit[1] >= 100000000) {
778 tdigit[1] -= 100000000;
779 if (tdigit[1] >= 100000000)
780 tdigit[1] -= 100000000;
781 }
782
783 digit = tdigit[0];
784 if (!digit && !tdigit[1])
785 nzeros += 16;
786 else {
787 if (!digit) {
788 nzeros += 8;
789 digit = tdigit[1];
790 }
791 // decompose digit
792 PD = (UINT64) digit *0x068DB8BBull;
793 digit_h = (UINT32) (PD >> 40);
794 digit_low = digit - digit_h * 10000;
795
796 if (!digit_low)
797 nzeros += 4;
798 else
799 digit_h = digit_low;
800
801 if (!(digit_h & 1))
802 nzeros +=
803 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
804 (digit_h & 7));
805 }
806
807 if (nzeros) {
808 // get P*(2^M[extra_digits])/10^extra_digits
809 __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]);
810
811 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
812 amount = recip_scale[nzeros];
813 __shr_128 (CQ, Qh, amount);
814 }
815 diff_expon += nzeros;
816
817 }
818 }
819 if(diff_expon>=0){
820 res =
821 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0],
822 rnd_mode, pfpsf);
823 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
824 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
825 #endif
826 BID_RETURN (res);
827 }
828 }
829 #endif
830
831 if (diff_expon >= 0) {
832 #ifdef IEEE_ROUND_NEAREST
833 // rounding
834 // 2*CA4 - CY
835 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
836 CA4r.w[0] = CA4.w[0] + CA4.w[0];
837 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
838 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
839
840 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
841 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
842
843 CQ.w[0] += carry64;
844 #else
845 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
846 // rounding
847 // 2*CA4 - CY
848 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
849 CA4r.w[0] = CA4.w[0] + CA4.w[0];
850 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
851 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
852
853 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
854 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
855
856 CQ.w[0] += carry64;
857 if (CQ.w[0] < carry64)
858 CQ.w[1]++;
859 #else
860 rmode = rnd_mode;
861 if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
862 rmode = 3 - rmode;
863 switch (rmode) {
864 case ROUNDING_TO_NEAREST: // round to nearest code
865 // rounding
866 // 2*CA4 - CY
867 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
868 CA4r.w[0] = CA4.w[0] + CA4.w[0];
869 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
870 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
871 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
872 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
873 CQ.w[0] += carry64;
874 if (CQ.w[0] < carry64)
875 CQ.w[1]++;
876 break;
877 case ROUNDING_TIES_AWAY:
878 // rounding
879 // 2*CA4 - CY
880 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
881 CA4r.w[0] = CA4.w[0] + CA4.w[0];
882 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
883 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
884 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
885 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
886 CQ.w[0] += carry64;
887 if (CQ.w[0] < carry64)
888 CQ.w[1]++;
889 break;
890 case ROUNDING_DOWN:
891 case ROUNDING_TO_ZERO:
892 break;
893 default: // rounding up
894 CQ.w[0]++;
895 if (!CQ.w[0])
896 CQ.w[1]++;
897 break;
898 }
899 #endif
900 #endif
901
902 res =
903 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0], rnd_mode,
904 pfpsf);
905 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
906 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
907 #endif
908 BID_RETURN (res);
909 } else {
910 // UF occurs
911
912 #ifdef SET_STATUS_FLAGS
913 if ((diff_expon + 16 < 0)) {
914 // set status flags
915 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
916 }
917 #endif
918 rmode = rnd_mode;
919 res =
920 get_BID64_UF (sign_x ^ sign_y, diff_expon, CQ.w[0], CA4.w[1] | CA4.w[0], rmode, pfpsf);
921 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
922 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
923 #endif
924 BID_RETURN (res);
925
926 }
927
928 }
929
930
931 //#define LEAVE_TRAILING_ZEROS
932
933 TYPE0_FUNCTION_ARG128_ARGTYPE2 (UINT64, bid64qd_div, x, UINT64, y)
934
935 UINT256 CA4 =
936 { {0x0ull, 0x0ull, 0x0ull, 0x0ull} }, CA4r, P256, QB256;
937 UINT128 CX, CY, T128, CQ, CQ2, CR, CA, TP128, Qh, Tmp;
938 UINT64 sign_x, sign_y, T, carry64, D, Q_low, QX, PD, res, valid_y;
939 int_float fx, fy, f64;
940 UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
941 int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
942 digits_q, amount;
943 int nzeros, i, j, k, d5, done = 0;
944 unsigned rmode;
945 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
946 fexcept_t binaryflags = 0;
947 #endif
948
949 valid_y = unpack_BID64 (&sign_y, &exponent_y, &CY.w[0], (y));
950
951 // unpack arguments, check for NaN or Infinity
952 if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
953 // test if x is NaN
954 if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
955 #ifdef SET_STATUS_FLAGS
956 if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull || // sNaN
957 (y & 0x7e00000000000000ull) == 0x7e00000000000000ull)
958 __set_status_flags (pfpsf, INVALID_EXCEPTION);
959 #endif
960 Tmp.w[1] = (CX.w[1] & 0x00003fffffffffffull);
961 Tmp.w[0] = CX.w[0];
962 TP128 = reciprocals10_128[18];
963 __mul_128x128_high (Qh, Tmp, TP128);
964 amount = recip_scale[18];
965 __shr_128 (Tmp, Qh, amount);
966 res = (CX.w[1] & 0xfc00000000000000ull) | Tmp.w[0];
967 BID_RETURN (res);
968 }
969 // x is Infinity?
970 if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
971 // check if y is Inf.
972 if (((y & 0x7c00000000000000ull) == 0x7800000000000000ull))
973 // return NaN
974 {
975 #ifdef SET_STATUS_FLAGS
976 __set_status_flags (pfpsf, INVALID_EXCEPTION);
977 #endif
978 res = 0x7c00000000000000ull;
979 BID_RETURN (res);
980 }
981 if (((y & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
982 // otherwise return +/-Inf
983 res =
984 ((x.w[1] ^ (y)) & 0x8000000000000000ull) | 0x7800000000000000ull;
985 BID_RETURN (res);
986 }
987 }
988 // x is 0
989 if (((y & INFINITY_MASK64) != INFINITY_MASK64) &&
990 !(CY.w[0])) {
991 #ifdef SET_STATUS_FLAGS
992 __set_status_flags (pfpsf, INVALID_EXCEPTION);
993 #endif
994 // x=y=0, return NaN
995 res = 0x7c00000000000000ull;
996 BID_RETURN (res);
997 }
998 // return 0
999 if (((y & 0x7800000000000000ull) != 0x7800000000000000ull)) {
1000 if (!CY.w[0]) {
1001 #ifdef SET_STATUS_FLAGS
1002 __set_status_flags (pfpsf, INVALID_EXCEPTION);
1003 #endif
1004 res = 0x7c00000000000000ull;
1005 BID_RETURN (res);
1006 }
1007 exponent_x =
1008 exponent_x - exponent_y - DECIMAL_EXPONENT_BIAS_128 +
1009 (DECIMAL_EXPONENT_BIAS << 1);
1010 if (exponent_x > DECIMAL_MAX_EXPON_64)
1011 exponent_x = DECIMAL_MAX_EXPON_64;
1012 else if (exponent_x < 0)
1013 exponent_x = 0;
1014 res = (sign_x ^ sign_y) | (((UINT64) exponent_x) << 53);
1015 BID_RETURN (res);
1016 }
1017 }
1018 CY.w[1] = 0;
1019 if (!valid_y) {
1020 // y is Inf. or NaN
1021
1022 // test if y is NaN
1023 if ((y & NAN_MASK64) == NAN_MASK64) {
1024 #ifdef SET_STATUS_FLAGS
1025 if ((y & SNAN_MASK64) == SNAN_MASK64) // sNaN
1026 __set_status_flags (pfpsf, INVALID_EXCEPTION);
1027 #endif
1028 BID_RETURN (CY.w[0] & QUIET_MASK64);
1029 }
1030 // y is Infinity?
1031 if (((y) & 0x7800000000000000ull) == 0x7800000000000000ull) {
1032 // return +/-0
1033 res = sign_x ^ sign_y;
1034 BID_RETURN (res);
1035 }
1036 // y is 0, return +/-Inf
1037 res =
1038 ((x.w[1] ^ (y)) & 0x8000000000000000ull) | 0x7800000000000000ull;
1039 #ifdef SET_STATUS_FLAGS
1040 __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
1041 #endif
1042 BID_RETURN (res);
1043 }
1044 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1045 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
1046 #endif
1047 diff_expon =
1048 exponent_x - exponent_y - DECIMAL_EXPONENT_BIAS_128 +
1049 (DECIMAL_EXPONENT_BIAS << 1);
1050
1051 if (__unsigned_compare_gt_128 (CY, CX)) {
1052 // CX < CY
1053
1054 // 2^64
1055 f64.i = 0x5f800000;
1056
1057 // fx ~ CX, fy ~ CY
1058 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
1059 fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
1060 // expon_cy - expon_cx
1061 bin_index = (fy.i - fx.i) >> 23;
1062
1063 if (CX.w[1]) {
1064 T = power10_index_binexp_128[bin_index].w[0];
1065 __mul_64x128_short (CA, T, CX);
1066 } else {
1067 T128 = power10_index_binexp_128[bin_index];
1068 __mul_64x128_short (CA, CX.w[0], T128);
1069 }
1070
1071 ed2 = 15;
1072 if (__unsigned_compare_gt_128 (CY, CA))
1073 ed2++;
1074
1075 T128 = power10_table_128[ed2];
1076 __mul_128x128_to_256 (CA4, CA, T128);
1077
1078 ed2 += estimate_decimal_digits[bin_index];
1079 CQ.w[0] = CQ.w[1] = 0;
1080 diff_expon = diff_expon - ed2;
1081
1082 } else {
1083 // get CQ = CX/CY
1084 __div_128_by_128 (&CQ, &CR, CX, CY);
1085
1086 // get number of decimal digits in CQ
1087 // 2^64
1088 f64.i = 0x5f800000;
1089 fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
1090 // binary expon. of CQ
1091 bin_expon = (fx.i - 0x3f800000) >> 23;
1092
1093 digits_q = estimate_decimal_digits[bin_expon];
1094 TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
1095 TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
1096 if (__unsigned_compare_ge_128 (CQ, TP128))
1097 digits_q++;
1098
1099 if (digits_q <= 16) {
1100 if (!CR.w[1] && !CR.w[0]) {
1101 res = get_BID64 (sign_x ^ sign_y, diff_expon,
1102 CQ.w[0], rnd_mode, pfpsf);
1103 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1104 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1105 #endif
1106 BID_RETURN (res);
1107 }
1108
1109 ed2 = 16 - digits_q;
1110 T128.w[0] = power10_table_128[ed2].w[0];
1111 __mul_64x128_to_192 (CA4, (T128.w[0]), CR);
1112 diff_expon = diff_expon - ed2;
1113 CQ.w[0] *= T128.w[0];
1114 } else {
1115 ed2 = digits_q - 16;
1116 diff_expon += ed2;
1117 T128 = reciprocals10_128[ed2];
1118 __mul_128x128_to_256 (P256, CQ, T128);
1119 amount = recip_scale[ed2];
1120 CQ.w[0] = (P256.w[2] >> amount) | (P256.w[3] << (64 - amount));
1121 CQ.w[1] = 0;
1122
1123 __mul_64x64_to_128 (CQ2, CQ.w[0], (power10_table_128[ed2].w[0]));
1124
1125 __mul_64x64_to_128 (QB256, CQ2.w[0], CY.w[0]);
1126 QB256.w[1] += CQ2.w[0] * CY.w[1] + CQ2.w[1] * CY.w[0];
1127
1128 CA4.w[1] = CX.w[1] - QB256.w[1];
1129 CA4.w[0] = CX.w[0] - QB256.w[0];
1130 if (CX.w[0] < QB256.w[0])
1131 CA4.w[1]--;
1132 if (CR.w[0] || CR.w[1])
1133 CA4.w[0] |= 1;
1134 done = 1;
1135 if(CA4.w[1]|CA4.w[0]) {
1136 __mul_64x128_low(CY, (power10_table_128[ed2].w[0]),CY);
1137 }
1138
1139 }
1140
1141 }
1142
1143 if (!done) {
1144 __div_256_by_128 (&CQ, &CA4, CY);
1145 }
1146
1147 #ifdef SET_STATUS_FLAGS
1148 if (CA4.w[0] || CA4.w[1]) {
1149 // set status flags
1150 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
1151 }
1152 #ifndef LEAVE_TRAILING_ZEROS
1153 else
1154 #endif
1155 #else
1156 #ifndef LEAVE_TRAILING_ZEROS
1157 if (!CA4.w[0] && !CA4.w[1])
1158 #endif
1159 #endif
1160 #ifndef LEAVE_TRAILING_ZEROS
1161 // check whether result is exact
1162 {
1163 if(!done) {
1164 // check whether CX, CY are short
1165 if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
1166 i = (int) CY.w[0] - 1;
1167 j = (int) CX.w[0] - 1;
1168 // difference in powers of 2 factors for Y and X
1169 nzeros = ed2 - factors[i][0] + factors[j][0];
1170 // difference in powers of 5 factors
1171 d5 = ed2 - factors[i][1] + factors[j][1];
1172 if (d5 < nzeros)
1173 nzeros = d5;
1174 // get P*(2^M[extra_digits])/10^extra_digits
1175 __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]);
1176 //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
1177
1178 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1179 amount = recip_scale[nzeros];
1180 __shr_128_long (CQ, Qh, amount);
1181
1182 diff_expon += nzeros;
1183 } else {
1184 // decompose Q as Qh*10^17 + Ql
1185 //T128 = reciprocals10_128[17];
1186 Q_low = CQ.w[0];
1187
1188 {
1189 tdigit[0] = Q_low & 0x3ffffff;
1190 tdigit[1] = 0;
1191 QX = Q_low >> 26;
1192 QX32 = QX;
1193 nzeros = 0;
1194
1195 for (j = 0; QX32; j++, QX32 >>= 7) {
1196 k = (QX32 & 127);
1197 tdigit[0] += convert_table[j][k][0];
1198 tdigit[1] += convert_table[j][k][1];
1199 if (tdigit[0] >= 100000000) {
1200 tdigit[0] -= 100000000;
1201 tdigit[1]++;
1202 }
1203 }
1204
1205 if (tdigit[1] >= 100000000) {
1206 tdigit[1] -= 100000000;
1207 if (tdigit[1] >= 100000000)
1208 tdigit[1] -= 100000000;
1209 }
1210
1211 digit = tdigit[0];
1212 if (!digit && !tdigit[1])
1213 nzeros += 16;
1214 else {
1215 if (!digit) {
1216 nzeros += 8;
1217 digit = tdigit[1];
1218 }
1219 // decompose digit
1220 PD = (UINT64) digit *0x068DB8BBull;
1221 digit_h = (UINT32) (PD >> 40);
1222 digit_low = digit - digit_h * 10000;
1223
1224 if (!digit_low)
1225 nzeros += 4;
1226 else
1227 digit_h = digit_low;
1228
1229 if (!(digit_h & 1))
1230 nzeros +=
1231 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
1232 (digit_h & 7));
1233 }
1234
1235 if (nzeros) {
1236 // get P*(2^M[extra_digits])/10^extra_digits
1237 __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]);
1238
1239 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1240 amount = recip_scale[nzeros];
1241 __shr_128 (CQ, Qh, amount);
1242 }
1243 diff_expon += nzeros;
1244
1245 }
1246 }
1247 }
1248 if(diff_expon>=0){
1249 res =
1250 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0],
1251 rnd_mode, pfpsf);
1252 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1253 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1254 #endif
1255 BID_RETURN (res);
1256 }
1257 }
1258 #endif
1259
1260 if (diff_expon >= 0) {
1261 #ifdef IEEE_ROUND_NEAREST
1262 // rounding
1263 // 2*CA4 - CY
1264 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1265 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1266 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1267 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1268
1269 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
1270 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
1271
1272 CQ.w[0] += carry64;
1273 //if(CQ.w[0]<carry64)
1274 //CQ.w[1] ++;
1275 #else
1276 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
1277 // rounding
1278 // 2*CA4 - CY
1279 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1280 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1281 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1282 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1283
1284 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
1285 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
1286
1287 CQ.w[0] += carry64;
1288 if (CQ.w[0] < carry64)
1289 CQ.w[1]++;
1290 #else
1291 rmode = rnd_mode;
1292 if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
1293 rmode = 3 - rmode;
1294 switch (rmode) {
1295 case ROUNDING_TO_NEAREST: // round to nearest code
1296 // rounding
1297 // 2*CA4 - CY
1298 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1299 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1300 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1301 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1302 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
1303 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
1304 CQ.w[0] += carry64;
1305 if (CQ.w[0] < carry64)
1306 CQ.w[1]++;
1307 break;
1308 case ROUNDING_TIES_AWAY:
1309 // rounding
1310 // 2*CA4 - CY
1311 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1312 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1313 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1314 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1315 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
1316 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
1317 CQ.w[0] += carry64;
1318 if (CQ.w[0] < carry64)
1319 CQ.w[1]++;
1320 break;
1321 case ROUNDING_DOWN:
1322 case ROUNDING_TO_ZERO:
1323 break;
1324 default: // rounding up
1325 CQ.w[0]++;
1326 if (!CQ.w[0])
1327 CQ.w[1]++;
1328 break;
1329 }
1330 #endif
1331 #endif
1332
1333
1334 res =
1335 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0], rnd_mode,
1336 pfpsf);
1337 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1338 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1339 #endif
1340 BID_RETURN (res);
1341 } else {
1342 // UF occurs
1343
1344 #ifdef SET_STATUS_FLAGS
1345 if ((diff_expon + 16 < 0)) {
1346 // set status flags
1347 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
1348 }
1349 #endif
1350 rmode = rnd_mode;
1351 res =
1352 get_BID64_UF (sign_x ^ sign_y, diff_expon, CQ.w[0], CA4.w[1] | CA4.w[0], rmode, pfpsf);
1353 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1354 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1355 #endif
1356 BID_RETURN (res);
1357
1358 }
1359
1360 }
1361
1362 //#define LEAVE_TRAILING_ZEROS
1363
1364 extern UINT32 convert_table[5][128][2];
1365 extern SINT8 factors[][2];
1366 extern UINT8 packed_10000_zeros[];
1367
1368
1369 //UINT64* bid64_div128x128(UINT64 res, UINT128 *px, UINT128 *py, unsigned rnd_mode, unsigned *pfpsf)
1370
1371 TYPE0_FUNCTION_ARG128_ARG128 (UINT64, bid64qq_div, x, y)
1372 UINT256 CA4 =
1373 { {0x0ull, 0x0ull, 0x0ull, 0x0ull} }, CA4r, P256, QB256;
1374 UINT128 CX, CY, T128, CQ, CQ2, CR, CA, TP128, Qh, Tmp;
1375 UINT64 sign_x, sign_y, T, carry64, D, Q_low, QX, valid_y, PD, res;
1376 int_float fx, fy, f64;
1377 UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
1378 int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2,
1379 digits_q, amount;
1380 int nzeros, i, j, k, d5, done = 0;
1381 unsigned rmode;
1382 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1383 fexcept_t binaryflags = 0;
1384 #endif
1385
1386 valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
1387
1388 // unpack arguments, check for NaN or Infinity
1389 if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
1390 // test if x is NaN
1391 if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
1392 #ifdef SET_STATUS_FLAGS
1393 if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull || // sNaN
1394 (y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull)
1395 __set_status_flags (pfpsf, INVALID_EXCEPTION);
1396 #endif
1397 Tmp.w[1] = (CX.w[1] & 0x00003fffffffffffull);
1398 Tmp.w[0] = CX.w[0];
1399 TP128 = reciprocals10_128[18];
1400 __mul_128x128_high (Qh, Tmp, TP128);
1401 amount = recip_scale[18];
1402 __shr_128 (Tmp, Qh, amount);
1403 res = (CX.w[1] & 0xfc00000000000000ull) | Tmp.w[0];
1404 BID_RETURN (res);
1405 }
1406 // x is Infinity?
1407 if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
1408 // check if y is Inf.
1409 if (((y.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull))
1410 // return NaN
1411 {
1412 #ifdef SET_STATUS_FLAGS
1413 __set_status_flags (pfpsf, INVALID_EXCEPTION);
1414 #endif
1415 res = 0x7c00000000000000ull;
1416 BID_RETURN (res);
1417 }
1418 if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull)) {
1419 // otherwise return +/-Inf
1420 res =
1421 ((x.w[1] ^ y.
1422 w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
1423 BID_RETURN (res);
1424 }
1425 }
1426 // x is 0
1427 if (((y.w[1] & 0x7800000000000000ull) != 0x7800000000000000ull)) {
1428 if ((!CY.w[0]) && !(CY.w[1] & 0x0001ffffffffffffull)) {
1429 #ifdef SET_STATUS_FLAGS
1430 __set_status_flags (pfpsf, INVALID_EXCEPTION);
1431 #endif
1432 // x=y=0, return NaN
1433 res = 0x7c00000000000000ull;
1434 BID_RETURN (res);
1435 }
1436 // return 0
1437 res = (x.w[1] ^ y.w[1]) & 0x8000000000000000ull;
1438 exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS;
1439 if (exponent_x > DECIMAL_MAX_EXPON_64)
1440 exponent_x = DECIMAL_MAX_EXPON_64;
1441 else if (exponent_x < 0)
1442 exponent_x = 0;
1443 res |= (((UINT64) exponent_x) << 53);
1444 BID_RETURN (res);
1445 }
1446 }
1447 if (!valid_y) {
1448 // y is Inf. or NaN
1449
1450 // test if y is NaN
1451 if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
1452 #ifdef SET_STATUS_FLAGS
1453 if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN
1454 __set_status_flags (pfpsf, INVALID_EXCEPTION);
1455 #endif
1456 Tmp.w[1] = (CY.w[1] & 0x00003fffffffffffull);
1457 Tmp.w[0] = CY.w[0];
1458 TP128 = reciprocals10_128[18];
1459 __mul_128x128_high (Qh, Tmp, TP128);
1460 amount = recip_scale[18];
1461 __shr_128 (Tmp, Qh, amount);
1462 res = (CY.w[1] & 0xfc00000000000000ull) | Tmp.w[0];
1463 BID_RETURN (res);
1464 }
1465 // y is Infinity?
1466 if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
1467 // return +/-0
1468 res = sign_x ^ sign_y;
1469 BID_RETURN (res);
1470 }
1471 // y is 0, return +/-Inf
1472 res =
1473 ((x.w[1] ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
1474 #ifdef SET_STATUS_FLAGS
1475 __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
1476 #endif
1477 BID_RETURN (res);
1478 }
1479 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1480 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
1481 #endif
1482 diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS;
1483
1484 if (__unsigned_compare_gt_128 (CY, CX)) {
1485 // CX < CY
1486
1487 // 2^64
1488 f64.i = 0x5f800000;
1489
1490 // fx ~ CX, fy ~ CY
1491 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
1492 fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
1493 // expon_cy - expon_cx
1494 bin_index = (fy.i - fx.i) >> 23;
1495
1496 if (CX.w[1]) {
1497 T = power10_index_binexp_128[bin_index].w[0];
1498 __mul_64x128_short (CA, T, CX);
1499 } else {
1500 T128 = power10_index_binexp_128[bin_index];
1501 __mul_64x128_short (CA, CX.w[0], T128);
1502 }
1503
1504 ed2 = 15;
1505 if (__unsigned_compare_gt_128 (CY, CA))
1506 ed2++;
1507
1508 T128 = power10_table_128[ed2];
1509 __mul_128x128_to_256 (CA4, CA, T128);
1510
1511 ed2 += estimate_decimal_digits[bin_index];
1512 CQ.w[0] = CQ.w[1] = 0;
1513 diff_expon = diff_expon - ed2;
1514
1515 } else {
1516 // get CQ = CX/CY
1517 __div_128_by_128 (&CQ, &CR, CX, CY);
1518
1519 // get number of decimal digits in CQ
1520 // 2^64
1521 f64.i = 0x5f800000;
1522 fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
1523 // binary expon. of CQ
1524 bin_expon = (fx.i - 0x3f800000) >> 23;
1525
1526 digits_q = estimate_decimal_digits[bin_expon];
1527 TP128.w[0] = power10_index_binexp_128[bin_expon].w[0];
1528 TP128.w[1] = power10_index_binexp_128[bin_expon].w[1];
1529 if (__unsigned_compare_ge_128 (CQ, TP128))
1530 digits_q++;
1531
1532 if (digits_q <= 16) {
1533 if (!CR.w[1] && !CR.w[0]) {
1534 res = get_BID64 (sign_x ^ sign_y, diff_expon,
1535 CQ.w[0], rnd_mode, pfpsf);
1536 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1537 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1538 #endif
1539 BID_RETURN (res);
1540 }
1541
1542 ed2 = 16 - digits_q;
1543 T128.w[0] = power10_table_128[ed2].w[0];
1544 __mul_64x128_to_192 (CA4, (T128.w[0]), CR);
1545 diff_expon = diff_expon - ed2;
1546 CQ.w[0] *= T128.w[0];
1547 } else {
1548 ed2 = digits_q - 16;
1549 diff_expon += ed2;
1550 T128 = reciprocals10_128[ed2];
1551 __mul_128x128_to_256 (P256, CQ, T128);
1552 amount = recip_scale[ed2];
1553 CQ.w[0] = (P256.w[2] >> amount) | (P256.w[3] << (64 - amount));
1554 CQ.w[1] = 0;
1555
1556 __mul_64x64_to_128 (CQ2, CQ.w[0], (power10_table_128[ed2].w[0]));
1557
1558 __mul_64x64_to_128 (QB256, CQ2.w[0], CY.w[0]);
1559 QB256.w[1] += CQ2.w[0] * CY.w[1] + CQ2.w[1] * CY.w[0];
1560
1561 CA4.w[1] = CX.w[1] - QB256.w[1];
1562 CA4.w[0] = CX.w[0] - QB256.w[0];
1563 if (CX.w[0] < QB256.w[0])
1564 CA4.w[1]--;
1565 if (CR.w[0] || CR.w[1])
1566 CA4.w[0] |= 1;
1567 done = 1;
1568 if(CA4.w[1]|CA4.w[0]) {
1569 __mul_64x128_low(CY, (power10_table_128[ed2].w[0]),CY);
1570 }
1571 }
1572
1573 }
1574
1575 if (!done) {
1576 __div_256_by_128 (&CQ, &CA4, CY);
1577 }
1578
1579
1580
1581 #ifdef SET_STATUS_FLAGS
1582 if (CA4.w[0] || CA4.w[1]) {
1583 // set status flags
1584 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
1585 }
1586 #ifndef LEAVE_TRAILING_ZEROS
1587 else
1588 #endif
1589 #else
1590 #ifndef LEAVE_TRAILING_ZEROS
1591 if (!CA4.w[0] && !CA4.w[1])
1592 #endif
1593 #endif
1594 #ifndef LEAVE_TRAILING_ZEROS
1595 // check whether result is exact
1596 {
1597 if(!done) {
1598 // check whether CX, CY are short
1599 if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
1600 i = (int) CY.w[0] - 1;
1601 j = (int) CX.w[0] - 1;
1602 // difference in powers of 2 factors for Y and X
1603 nzeros = ed2 - factors[i][0] + factors[j][0];
1604 // difference in powers of 5 factors
1605 d5 = ed2 - factors[i][1] + factors[j][1];
1606 if (d5 < nzeros)
1607 nzeros = d5;
1608 // get P*(2^M[extra_digits])/10^extra_digits
1609 __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]);
1610 //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
1611
1612 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1613 amount = recip_scale[nzeros];
1614 __shr_128_long (CQ, Qh, amount);
1615
1616 diff_expon += nzeros;
1617 } else {
1618 // decompose Q as Qh*10^17 + Ql
1619 //T128 = reciprocals10_128[17];
1620 Q_low = CQ.w[0];
1621
1622 {
1623 tdigit[0] = Q_low & 0x3ffffff;
1624 tdigit[1] = 0;
1625 QX = Q_low >> 26;
1626 QX32 = QX;
1627 nzeros = 0;
1628
1629 for (j = 0; QX32; j++, QX32 >>= 7) {
1630 k = (QX32 & 127);
1631 tdigit[0] += convert_table[j][k][0];
1632 tdigit[1] += convert_table[j][k][1];
1633 if (tdigit[0] >= 100000000) {
1634 tdigit[0] -= 100000000;
1635 tdigit[1]++;
1636 }
1637 }
1638
1639 if (tdigit[1] >= 100000000) {
1640 tdigit[1] -= 100000000;
1641 if (tdigit[1] >= 100000000)
1642 tdigit[1] -= 100000000;
1643 }
1644
1645 digit = tdigit[0];
1646 if (!digit && !tdigit[1])
1647 nzeros += 16;
1648 else {
1649 if (!digit) {
1650 nzeros += 8;
1651 digit = tdigit[1];
1652 }
1653 // decompose digit
1654 PD = (UINT64) digit *0x068DB8BBull;
1655 digit_h = (UINT32) (PD >> 40);
1656 digit_low = digit - digit_h * 10000;
1657
1658 if (!digit_low)
1659 nzeros += 4;
1660 else
1661 digit_h = digit_low;
1662
1663 if (!(digit_h & 1))
1664 nzeros +=
1665 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >>
1666 (digit_h & 7));
1667 }
1668
1669 if (nzeros) {
1670 // get P*(2^M[extra_digits])/10^extra_digits
1671 __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]);
1672
1673 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1674 amount = recip_scale[nzeros];
1675 __shr_128 (CQ, Qh, amount);
1676 }
1677 diff_expon += nzeros;
1678
1679 }
1680 }
1681 }
1682 if(diff_expon>=0){
1683 res =
1684 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0],
1685 rnd_mode, pfpsf);
1686 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1687 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1688 #endif
1689 BID_RETURN (res);
1690 }
1691 }
1692 #endif
1693
1694 if(diff_expon>=0) {
1695
1696 #ifdef IEEE_ROUND_NEAREST
1697 // rounding
1698 // 2*CA4 - CY
1699 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1700 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1701 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1702 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1703
1704 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
1705 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
1706
1707 CQ.w[0] += carry64;
1708 //if(CQ.w[0]<carry64)
1709 //CQ.w[1] ++;
1710 #else
1711 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
1712 // rounding
1713 // 2*CA4 - CY
1714 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1715 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1716 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1717 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1718
1719 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
1720 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
1721
1722 CQ.w[0] += carry64;
1723 if (CQ.w[0] < carry64)
1724 CQ.w[1]++;
1725 #else
1726 rmode = rnd_mode;
1727 if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
1728 rmode = 3 - rmode;
1729 switch (rmode) {
1730 case ROUNDING_TO_NEAREST: // round to nearest code
1731 // rounding
1732 // 2*CA4 - CY
1733 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1734 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1735 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1736 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1737 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
1738 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
1739 CQ.w[0] += carry64;
1740 if (CQ.w[0] < carry64)
1741 CQ.w[1]++;
1742 break;
1743 case ROUNDING_TIES_AWAY:
1744 // rounding
1745 // 2*CA4 - CY
1746 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
1747 CA4r.w[0] = CA4.w[0] + CA4.w[0];
1748 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
1749 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
1750 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
1751 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
1752 CQ.w[0] += carry64;
1753 if (CQ.w[0] < carry64)
1754 CQ.w[1]++;
1755 break;
1756 case ROUNDING_DOWN:
1757 case ROUNDING_TO_ZERO:
1758 break;
1759 default: // rounding up
1760 CQ.w[0]++;
1761 if (!CQ.w[0])
1762 CQ.w[1]++;
1763 break;
1764 }
1765 #endif
1766 #endif
1767
1768
1769 res =
1770 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0], rnd_mode,
1771 pfpsf);
1772 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1773 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1774 #endif
1775 BID_RETURN (res);
1776 } else {
1777 // UF occurs
1778
1779 #ifdef SET_STATUS_FLAGS
1780 if ((diff_expon + 16 < 0)) {
1781 // set status flags
1782 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
1783 }
1784 #endif
1785 rmode = rnd_mode;
1786 res =
1787 get_BID64_UF (sign_x ^ sign_y, diff_expon, CQ.w[0], CA4.w[1] | CA4.w[0], rmode, pfpsf);
1788 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1789 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
1790 #endif
1791 BID_RETURN (res);
1792
1793 }
1794
1795 }
Mirror of https://gcc.gnu.org/git/gcc.git