1 /* Copyright (C) 2007-2024 Free Software Foundation, Inc.
3 This file is part of GCC.
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
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
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.
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/>. */
25 #include "bid_div_macros.h"
26 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
29 #define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT
32 extern UINT32 convert_table
[5][128][2];
33 extern SINT8 factors
[][2];
34 extern UINT8 packed_10000_zeros
[];
36 BID128_FUNCTION_ARG2 (bid128_div
, x
, y
)
38 UINT256 CA4
, CA4r
, P256
;
39 UINT128 CX
, CY
, T128
, CQ
, CR
, CA
, TP128
, Qh
, res
;
40 UINT64 sign_x
, sign_y
, T
, carry64
, D
, Q_high
, Q_low
, QX
, PD
,
42 int_float fx
, fy
, f64
;
43 UINT32 QX32
, tdigit
[3], digit
, digit_h
, digit_low
;
44 int exponent_x
, exponent_y
, bin_index
, bin_expon
, diff_expon
, ed2
,
46 int nzeros
, i
, j
, k
, d5
;
48 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
49 fexcept_t binaryflags
= 0;
52 valid_y
= unpack_BID128_value (&sign_y
, &exponent_y
, &CY
, y
);
54 // unpack arguments, check for NaN or Infinity
55 if (!unpack_BID128_value (&sign_x
, &exponent_x
, &CX
, x
)) {
57 if ((x
.w
[1] & 0x7c00000000000000ull
) == 0x7c00000000000000ull
) {
58 #ifdef SET_STATUS_FLAGS
59 if ((x
.w
[1] & 0x7e00000000000000ull
) == 0x7e00000000000000ull
|| // sNaN
60 (y
.w
[1] & 0x7e00000000000000ull
) == 0x7e00000000000000ull
)
61 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
63 res
.w
[1] = (CX
.w
[1]) & QUIET_MASK64
;
68 if ((x
.w
[1] & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
70 if (((y
.w
[1] & 0x7c00000000000000ull
) == 0x7800000000000000ull
))
73 #ifdef SET_STATUS_FLAGS
74 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
76 res
.w
[1] = 0x7c00000000000000ull
;
81 if (((y
.w
[1] & 0x7c00000000000000ull
) != 0x7c00000000000000ull
))
85 res
.w
[1] = ((x
.w
[1] ^ y
.w
[1]) & 0x8000000000000000ull
) |
86 0x7800000000000000ull
;
92 if ((y
.w
[1] & 0x7800000000000000ull
) < 0x7800000000000000ull
) {
93 if ((!CY
.w
[0]) && !(CY
.w
[1] & 0x0001ffffffffffffull
)) {
94 #ifdef SET_STATUS_FLAGS
95 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
98 res
.w
[1] = 0x7c00000000000000ull
;
103 res
.w
[1] = (x
.w
[1] ^ y
.w
[1]) & 0x8000000000000000ull
;
104 exponent_x
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS_128
;
105 if (exponent_x
> DECIMAL_MAX_EXPON_128
)
106 exponent_x
= DECIMAL_MAX_EXPON_128
;
107 else if (exponent_x
< 0)
109 res
.w
[1] |= (((UINT64
) exponent_x
) << 49);
118 if ((y
.w
[1] & 0x7c00000000000000ull
) == 0x7c00000000000000ull
) {
119 #ifdef SET_STATUS_FLAGS
120 if ((y
.w
[1] & 0x7e00000000000000ull
) == 0x7e00000000000000ull
) // sNaN
121 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
123 res
.w
[1] = CY
.w
[1] & QUIET_MASK64
;
128 if ((y
.w
[1] & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
130 res
.w
[1] = sign_x
^ sign_y
;
134 // y is 0, return +/-Inf
135 #ifdef SET_STATUS_FLAGS
136 __set_status_flags (pfpsf
, ZERO_DIVIDE_EXCEPTION
);
139 ((x
.w
[1] ^ y
.w
[1]) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
143 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
144 (void) fegetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
146 diff_expon
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS_128
;
148 if (__unsigned_compare_gt_128 (CY
, CX
)) {
155 fx
.d
= (float) CX
.w
[1] * f64
.d
+ (float) CX
.w
[0];
156 fy
.d
= (float) CY
.w
[1] * f64
.d
+ (float) CY
.w
[0];
157 // expon_cy - expon_cx
158 bin_index
= (fy
.i
- fx
.i
) >> 23;
161 T
= power10_index_binexp_128
[bin_index
].w
[0];
162 __mul_64x128_short (CA
, T
, CX
);
164 T128
= power10_index_binexp_128
[bin_index
];
165 __mul_64x128_short (CA
, CX
.w
[0], T128
);
169 if (__unsigned_compare_gt_128 (CY
, CA
))
172 T128
= power10_table_128
[ed2
];
173 __mul_128x128_to_256 (CA4
, CA
, T128
);
175 ed2
+= estimate_decimal_digits
[bin_index
];
176 CQ
.w
[0] = CQ
.w
[1] = 0;
177 diff_expon
= diff_expon
- ed2
;
181 __div_128_by_128 (&CQ
, &CR
, CX
, CY
);
183 if (!CR
.w
[1] && !CR
.w
[0]) {
184 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
,
186 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
187 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
191 // get number of decimal digits in CQ
194 fx
.d
= (float) CQ
.w
[1] * f64
.d
+ (float) CQ
.w
[0];
195 // binary expon. of CQ
196 bin_expon
= (fx
.i
- 0x3f800000) >> 23;
198 digits_q
= estimate_decimal_digits
[bin_expon
];
199 TP128
.w
[0] = power10_index_binexp_128
[bin_expon
].w
[0];
200 TP128
.w
[1] = power10_index_binexp_128
[bin_expon
].w
[1];
201 if (__unsigned_compare_ge_128 (CQ
, TP128
))
205 T128
.w
[0] = power10_table_128
[ed2
].w
[0];
206 T128
.w
[1] = power10_table_128
[ed2
].w
[1];
207 __mul_128x128_to_256 (CA4
, CR
, T128
);
208 diff_expon
= diff_expon
- ed2
;
209 __mul_128x128_low (CQ
, CQ
, T128
);
213 __div_256_by_128 (&CQ
, &CA4
, CY
);
215 #ifdef SET_STATUS_FLAGS
216 if (CA4
.w
[0] || CA4
.w
[1]) {
218 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
220 #ifndef LEAVE_TRAILING_ZEROS
224 #ifndef LEAVE_TRAILING_ZEROS
225 if (!CA4
.w
[0] && !CA4
.w
[1])
228 #ifndef LEAVE_TRAILING_ZEROS
229 // check whether result is exact
231 // check whether CX, CY are short
232 if (!CX
.w
[1] && !CY
.w
[1] && (CX
.w
[0] <= 1024) && (CY
.w
[0] <= 1024)) {
233 i
= (int) CY
.w
[0] - 1;
234 j
= (int) CX
.w
[0] - 1;
235 // difference in powers of 2 factors for Y and X
236 nzeros
= ed2
- factors
[i
][0] + factors
[j
][0];
237 // difference in powers of 5 factors
238 d5
= ed2
- factors
[i
][1] + factors
[j
][1];
241 // get P*(2^M[extra_digits])/10^extra_digits
242 __mul_128x128_high (Qh
, CQ
, reciprocals10_128
[nzeros
]);
244 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
245 amount
= recip_scale
[nzeros
];
246 __shr_128_long (CQ
, Qh
, amount
);
248 diff_expon
+= nzeros
;
250 // decompose Q as Qh*10^17 + Ql
251 //T128 = reciprocals10_128[17];
252 T128
.w
[0] = 0x44909befeb9fad49ull
;
253 T128
.w
[1] = 0x000b877aa3236a4bull
;
254 __mul_128x128_to_256 (P256
, CQ
, T128
);
255 //amount = recip_scale[17];
256 Q_high
= (P256
.w
[2] >> 44) | (P256
.w
[3] << (64 - 44));
257 Q_low
= CQ
.w
[0] - Q_high
* 100000000000000000ull;
262 tdigit
[0] = Q_high
& 0x3ffffff;
268 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
270 tdigit
[0] += convert_table
[j
][k
][0];
271 tdigit
[1] += convert_table
[j
][k
][1];
272 if (tdigit
[0] >= 100000000) {
273 tdigit
[0] -= 100000000;
278 if (tdigit
[1] >= 100000000) {
279 tdigit
[1] -= 100000000;
280 if (tdigit
[1] >= 100000000)
281 tdigit
[1] -= 100000000;
285 if (!digit
&& !tdigit
[1])
293 PD
= (UINT64
) digit
*0x068DB8BBull
;
294 digit_h
= (UINT32
) (PD
>> 40);
295 digit_low
= digit
- digit_h
* 10000;
304 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
309 __mul_64x64_to_128 (CQ
, Q_high
, reciprocals10_64
[nzeros
]);
311 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
312 amount
= short_recip_scale
[nzeros
];
313 CQ
.w
[0] = CQ
.w
[1] >> amount
;
318 diff_expon
+= nzeros
;
320 tdigit
[0] = Q_low
& 0x3ffffff;
326 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
328 tdigit
[0] += convert_table
[j
][k
][0];
329 tdigit
[1] += convert_table
[j
][k
][1];
330 if (tdigit
[0] >= 100000000) {
331 tdigit
[0] -= 100000000;
336 if (tdigit
[1] >= 100000000) {
337 tdigit
[1] -= 100000000;
338 if (tdigit
[1] >= 100000000)
339 tdigit
[1] -= 100000000;
343 if (!digit
&& !tdigit
[1])
351 PD
= (UINT64
) digit
*0x068DB8BBull
;
352 digit_h
= (UINT32
) (PD
>> 40);
353 digit_low
= digit
- digit_h
* 10000;
362 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
367 // get P*(2^M[extra_digits])/10^extra_digits
368 __mul_128x128_high (Qh
, CQ
, reciprocals10_128
[nzeros
]);
370 //now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
371 amount
= recip_scale
[nzeros
];
372 __shr_128 (CQ
, Qh
, amount
);
374 diff_expon
+= nzeros
;
378 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
, pfpsf
);
379 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
380 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
386 if (diff_expon
>= 0) {
387 #ifdef IEEE_ROUND_NEAREST
390 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
391 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
392 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
393 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
395 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
396 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
399 if (CQ
.w
[0] < carry64
)
402 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
405 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
406 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
407 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
408 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
410 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
411 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
414 if (CQ
.w
[0] < carry64
)
418 if (sign_x
^ sign_y
&& (unsigned) (rmode
- 1) < 2)
421 case ROUNDING_TO_NEAREST
: // round to nearest code
424 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
425 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
426 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
427 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
428 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
429 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
431 if (CQ
.w
[0] < carry64
)
434 case ROUNDING_TIES_AWAY
:
437 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
438 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
439 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
440 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
441 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
442 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
444 if (CQ
.w
[0] < carry64
)
448 case ROUNDING_TO_ZERO
:
450 default: // rounding up
460 #ifdef SET_STATUS_FLAGS
461 if (CA4
.w
[0] || CA4
.w
[1]) {
463 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
467 handle_UF_128_rem (&res
, sign_x
^ sign_y
, diff_expon
, CQ
,
468 CA4
.w
[1] | CA4
.w
[0], &rnd_mode
, pfpsf
);
469 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
470 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
476 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
, pfpsf
);
477 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
478 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
484 //#define LEAVE_TRAILING_ZEROS
486 TYPE0_FUNCTION_ARGTYPE1_ARGTYPE2 (UINT128
, bid128dd_div
, UINT64
, x
,
489 UINT256 CA4
, CA4r
, P256
;
490 UINT128 CX
, CY
, T128
, CQ
, CR
, CA
, TP128
, Qh
, res
;
491 UINT64 sign_x
, sign_y
, T
, carry64
, D
, Q_high
, Q_low
, QX
, PD
,
493 int_float fx
, fy
, f64
;
494 UINT32 QX32
, tdigit
[3], digit
, digit_h
, digit_low
;
495 int exponent_x
, exponent_y
, bin_index
, bin_expon
, diff_expon
, ed2
,
497 int nzeros
, i
, j
, k
, d5
;
499 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
500 fexcept_t binaryflags
= 0;
503 valid_y
= unpack_BID64 (&sign_y
, &exponent_y
, &CY
.w
[0], y
);
505 // unpack arguments, check for NaN or Infinity
507 if (!unpack_BID64 (&sign_x
, &exponent_x
, &CX
.w
[0], (x
))) {
508 #ifdef SET_STATUS_FLAGS
509 if ((y
& SNAN_MASK64
) == SNAN_MASK64
) // y is sNaN
510 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
514 if ((x
& NAN_MASK64
) == NAN_MASK64
) {
515 #ifdef SET_STATUS_FLAGS
516 if ((x
& SNAN_MASK64
) == SNAN_MASK64
) // sNaN
517 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
519 res
.w
[0] = (CX
.w
[0] & 0x0003ffffffffffffull
);
520 __mul_64x64_to_128 (res
, res
.w
[0], power10_table_128
[18].w
[0]);
521 res
.w
[1] |= ((CX
.w
[0]) & 0xfc00000000000000ull
);
525 if (((x
) & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
526 // check if y is Inf.
527 if ((((y
) & 0x7c00000000000000ull
) == 0x7800000000000000ull
))
530 #ifdef SET_STATUS_FLAGS
531 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
533 res
.w
[1] = 0x7c00000000000000ull
;
537 if ((((y
) & 0x7c00000000000000ull
) != 0x7c00000000000000ull
)) {
538 // otherwise return +/-Inf
540 (((x
) ^ (y
)) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
546 if ((((y
) & 0x7800000000000000ull
) != 0x7800000000000000ull
)) {
548 #ifdef SET_STATUS_FLAGS
549 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
552 res
.w
[1] = 0x7c00000000000000ull
;
557 res
.w
[1] = ((x
) ^ (y
)) & 0x8000000000000000ull
;
558 if (((y
) & 0x6000000000000000ull
) == 0x6000000000000000ull
)
559 exponent_y
= ((UINT32
) ((y
) >> 51)) & 0x3ff;
561 exponent_y
= ((UINT32
) ((y
) >> 53)) & 0x3ff;
562 exponent_x
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS_128
;
563 if (exponent_x
> DECIMAL_MAX_EXPON_128
)
564 exponent_x
= DECIMAL_MAX_EXPON_128
;
565 else if (exponent_x
< 0)
567 res
.w
[1] |= (((UINT64
) exponent_x
) << 49);
578 if ((y
& NAN_MASK64
) == NAN_MASK64
) {
579 #ifdef SET_STATUS_FLAGS
580 if ((y
& SNAN_MASK64
) == SNAN_MASK64
) // sNaN
581 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
583 res
.w
[0] = (CY
.w
[0] & 0x0003ffffffffffffull
);
584 __mul_64x64_to_128 (res
, res
.w
[0], power10_table_128
[18].w
[0]);
585 res
.w
[1] |= ((CY
.w
[0]) & 0xfc00000000000000ull
);
589 if (((y
) & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
591 res
.w
[1] = sign_x
^ sign_y
;
595 // y is 0, return +/-Inf
597 (((x
) ^ (y
)) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
599 #ifdef SET_STATUS_FLAGS
600 __set_status_flags (pfpsf
, ZERO_DIVIDE_EXCEPTION
);
604 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
605 (void) fegetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
607 diff_expon
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS_128
;
609 if (__unsigned_compare_gt_128 (CY
, CX
)) {
616 fx
.d
= (float) CX
.w
[1] * f64
.d
+ (float) CX
.w
[0];
617 fy
.d
= (float) CY
.w
[1] * f64
.d
+ (float) CY
.w
[0];
618 // expon_cy - expon_cx
619 bin_index
= (fy
.i
- fx
.i
) >> 23;
622 T
= power10_index_binexp_128
[bin_index
].w
[0];
623 __mul_64x128_short (CA
, T
, CX
);
625 T128
= power10_index_binexp_128
[bin_index
];
626 __mul_64x128_short (CA
, CX
.w
[0], T128
);
630 if (__unsigned_compare_gt_128 (CY
, CA
))
633 T128
= power10_table_128
[ed2
];
634 __mul_128x128_to_256 (CA4
, CA
, T128
);
636 ed2
+= estimate_decimal_digits
[bin_index
];
637 CQ
.w
[0] = CQ
.w
[1] = 0;
638 diff_expon
= diff_expon
- ed2
;
642 __div_128_by_128 (&CQ
, &CR
, CX
, CY
);
644 if (!CR
.w
[1] && !CR
.w
[0]) {
645 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
,
647 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
648 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
652 // get number of decimal digits in CQ
655 fx
.d
= (float) CQ
.w
[1] * f64
.d
+ (float) CQ
.w
[0];
656 // binary expon. of CQ
657 bin_expon
= (fx
.i
- 0x3f800000) >> 23;
659 digits_q
= estimate_decimal_digits
[bin_expon
];
660 TP128
.w
[0] = power10_index_binexp_128
[bin_expon
].w
[0];
661 TP128
.w
[1] = power10_index_binexp_128
[bin_expon
].w
[1];
662 if (__unsigned_compare_ge_128 (CQ
, TP128
))
666 T128
.w
[0] = power10_table_128
[ed2
].w
[0];
667 T128
.w
[1] = power10_table_128
[ed2
].w
[1];
668 __mul_128x128_to_256 (CA4
, CR
, T128
);
669 diff_expon
= diff_expon
- ed2
;
670 __mul_128x128_low (CQ
, CQ
, T128
);
674 __div_256_by_128 (&CQ
, &CA4
, CY
);
677 #ifdef SET_STATUS_FLAGS
678 if (CA4
.w
[0] || CA4
.w
[1]) {
680 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
682 #ifndef LEAVE_TRAILING_ZEROS
686 #ifndef LEAVE_TRAILING_ZEROS
687 if (!CA4
.w
[0] && !CA4
.w
[1])
690 #ifndef LEAVE_TRAILING_ZEROS
691 // check whether result is exact
693 // check whether CX, CY are short
694 if (!CX
.w
[1] && !CY
.w
[1] && (CX
.w
[0] <= 1024) && (CY
.w
[0] <= 1024)) {
695 i
= (int) CY
.w
[0] - 1;
696 j
= (int) CX
.w
[0] - 1;
697 // difference in powers of 2 factors for Y and X
698 nzeros
= ed2
- factors
[i
][0] + factors
[j
][0];
699 // difference in powers of 5 factors
700 d5
= ed2
- factors
[i
][1] + factors
[j
][1];
703 // get P*(2^M[extra_digits])/10^extra_digits
704 __mul_128x128_high (Qh
, CQ
, reciprocals10_128
[nzeros
]);
705 //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
707 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
708 amount
= recip_scale
[nzeros
];
709 __shr_128_long (CQ
, Qh
, amount
);
711 diff_expon
+= nzeros
;
713 // decompose Q as Qh*10^17 + Ql
714 //T128 = reciprocals10_128[17];
715 T128
.w
[0] = 0x44909befeb9fad49ull
;
716 T128
.w
[1] = 0x000b877aa3236a4bull
;
717 __mul_128x128_to_256 (P256
, CQ
, T128
);
718 //amount = recip_scale[17];
719 Q_high
= (P256
.w
[2] >> 44) | (P256
.w
[3] << (64 - 44));
720 Q_low
= CQ
.w
[0] - Q_high
* 100000000000000000ull;
725 tdigit
[0] = Q_high
& 0x3ffffff;
731 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
733 tdigit
[0] += convert_table
[j
][k
][0];
734 tdigit
[1] += convert_table
[j
][k
][1];
735 if (tdigit
[0] >= 100000000) {
736 tdigit
[0] -= 100000000;
742 if (tdigit
[1] >= 100000000) {
743 tdigit
[1] -= 100000000;
744 if (tdigit
[1] >= 100000000)
745 tdigit
[1] -= 100000000;
749 if (!digit
&& !tdigit
[1])
757 PD
= (UINT64
) digit
*0x068DB8BBull
;
758 digit_h
= (UINT32
) (PD
>> 40);
759 digit_low
= digit
- digit_h
* 10000;
768 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
773 __mul_64x64_to_128 (CQ
, Q_high
, reciprocals10_64
[nzeros
]);
775 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
776 amount
= short_recip_scale
[nzeros
];
777 CQ
.w
[0] = CQ
.w
[1] >> amount
;
782 diff_expon
+= nzeros
;
784 tdigit
[0] = Q_low
& 0x3ffffff;
790 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
792 tdigit
[0] += convert_table
[j
][k
][0];
793 tdigit
[1] += convert_table
[j
][k
][1];
794 if (tdigit
[0] >= 100000000) {
795 tdigit
[0] -= 100000000;
800 if (tdigit
[1] >= 100000000) {
801 tdigit
[1] -= 100000000;
802 if (tdigit
[1] >= 100000000)
803 tdigit
[1] -= 100000000;
807 if (!digit
&& !tdigit
[1])
815 PD
= (UINT64
) digit
*0x068DB8BBull
;
816 digit_h
= (UINT32
) (PD
>> 40);
817 digit_low
= digit
- digit_h
* 10000;
826 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
831 // get P*(2^M[extra_digits])/10^extra_digits
832 __mul_128x128_high (Qh
, CQ
, reciprocals10_128
[nzeros
]);
834 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
835 amount
= recip_scale
[nzeros
];
836 __shr_128 (CQ
, Qh
, amount
);
838 diff_expon
+= nzeros
;
842 get_BID128(&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
,pfpsf
);
843 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
844 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
850 if (diff_expon
>= 0) {
851 #ifdef IEEE_ROUND_NEAREST
854 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
855 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
856 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
857 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
859 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
860 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
863 if (CQ
.w
[0] < carry64
)
866 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
869 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
870 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
871 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
872 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
874 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
875 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
878 if (CQ
.w
[0] < carry64
)
882 if (sign_x
^ sign_y
&& (unsigned) (rmode
- 1) < 2)
885 case ROUNDING_TO_NEAREST
: // round to nearest code
888 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
889 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
890 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
891 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
892 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
893 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
895 if (CQ
.w
[0] < carry64
)
898 case ROUNDING_TIES_AWAY
:
901 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
902 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
903 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
904 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
905 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
906 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
908 if (CQ
.w
[0] < carry64
)
912 case ROUNDING_TO_ZERO
:
914 default: // rounding up
924 #ifdef SET_STATUS_FLAGS
925 if (CA4
.w
[0] || CA4
.w
[1]) {
927 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
930 handle_UF_128_rem (&res
, sign_x
^ sign_y
, diff_expon
, CQ
,
931 CA4
.w
[1] | CA4
.w
[0], &rnd_mode
, pfpsf
);
932 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
933 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
939 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
, pfpsf
);
940 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
941 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
947 BID128_FUNCTION_ARGTYPE1_ARG128 (bid128dq_div
, UINT64
, x
, y
)
948 UINT256 CA4
, CA4r
, P256
;
949 UINT128 CX
, CY
, T128
, CQ
, CR
, CA
, TP128
, Qh
, res
;
950 UINT64 sign_x
, sign_y
, T
, carry64
, D
, Q_high
, Q_low
, QX
, valid_y
,
952 int_float fx
, fy
, f64
;
953 UINT32 QX32
, tdigit
[3], digit
, digit_h
, digit_low
;
954 int exponent_x
, exponent_y
, bin_index
, bin_expon
, diff_expon
, ed2
,
956 int nzeros
, i
, j
, k
, d5
;
958 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
959 fexcept_t binaryflags
= 0;
962 valid_y
= unpack_BID128_value (&sign_y
, &exponent_y
, &CY
, y
);
964 // unpack arguments, check for NaN or Infinity
966 if (!unpack_BID64 (&sign_x
, &exponent_x
, &CX
.w
[0], x
)) {
967 #ifdef SET_STATUS_FLAGS
968 if ((y
.w
[1] & SNAN_MASK64
) == SNAN_MASK64
) // y is sNaN
969 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
973 if ((x
& NAN_MASK64
) == NAN_MASK64
) {
974 #ifdef SET_STATUS_FLAGS
975 if ((x
& SNAN_MASK64
) == SNAN_MASK64
) // sNaN
976 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
978 res
.w
[0] = (CX
.w
[0] & 0x0003ffffffffffffull
);
979 __mul_64x64_to_128 (res
, res
.w
[0], power10_table_128
[18].w
[0]);
980 res
.w
[1] |= ((CX
.w
[0]) & 0xfc00000000000000ull
);
984 if ((x
& 0x7800000000000000ull
) == 0x7800000000000000ull
) {
985 // check if y is Inf.
986 if (((y
.w
[1] & 0x7c00000000000000ull
) == 0x7800000000000000ull
))
989 #ifdef SET_STATUS_FLAGS
990 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
992 res
.w
[1] = 0x7c00000000000000ull
;
996 if (((y
.w
[1] & 0x7c00000000000000ull
) != 0x7c00000000000000ull
)) {
997 // otherwise return +/-Inf
999 ((x
^ y
.w
[1]) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
1005 if ((y
.w
[1] & INFINITY_MASK64
) != INFINITY_MASK64
) {
1006 if ((!CY
.w
[0]) && !(CY
.w
[1] & 0x0001ffffffffffffull
)) {
1007 #ifdef SET_STATUS_FLAGS
1008 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
1010 // x=y=0, return NaN
1011 res
.w
[1] = 0x7c00000000000000ull
;
1016 res
.w
[1] = (x
^ y
.w
[1]) & 0x8000000000000000ull
;
1017 exponent_x
= exponent_x
- exponent_y
+ (DECIMAL_EXPONENT_BIAS_128
<<1) - DECIMAL_EXPONENT_BIAS
;
1018 if (exponent_x
> DECIMAL_MAX_EXPON_128
)
1019 exponent_x
= DECIMAL_MAX_EXPON_128
;
1020 else if (exponent_x
< 0)
1022 res
.w
[1] |= (((UINT64
) exponent_x
) << 49);
1027 exponent_x
+= (DECIMAL_EXPONENT_BIAS_128
- DECIMAL_EXPONENT_BIAS
);
1033 if ((y
.w
[1] & 0x7c00000000000000ull
) == 0x7c00000000000000ull
) {
1034 #ifdef SET_STATUS_FLAGS
1035 if ((y
.w
[1] & 0x7e00000000000000ull
) == 0x7e00000000000000ull
) // sNaN
1036 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
1038 res
.w
[1] = CY
.w
[1] & QUIET_MASK64
;
1043 if ((y
.w
[1] & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
1045 res
.w
[1] = sign_x
^ sign_y
;
1049 // y is 0, return +/-Inf
1051 ((x
^ y
.w
[1]) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
1053 #ifdef SET_STATUS_FLAGS
1054 __set_status_flags (pfpsf
, ZERO_DIVIDE_EXCEPTION
);
1058 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1059 (void) fegetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1061 diff_expon
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS_128
;
1063 if (__unsigned_compare_gt_128 (CY
, CX
)) {
1070 fx
.d
= (float) CX
.w
[1] * f64
.d
+ (float) CX
.w
[0];
1071 fy
.d
= (float) CY
.w
[1] * f64
.d
+ (float) CY
.w
[0];
1072 // expon_cy - expon_cx
1073 bin_index
= (fy
.i
- fx
.i
) >> 23;
1076 T
= power10_index_binexp_128
[bin_index
].w
[0];
1077 __mul_64x128_short (CA
, T
, CX
);
1079 T128
= power10_index_binexp_128
[bin_index
];
1080 __mul_64x128_short (CA
, CX
.w
[0], T128
);
1084 if (__unsigned_compare_gt_128 (CY
, CA
))
1087 T128
= power10_table_128
[ed2
];
1088 __mul_128x128_to_256 (CA4
, CA
, T128
);
1090 ed2
+= estimate_decimal_digits
[bin_index
];
1091 CQ
.w
[0] = CQ
.w
[1] = 0;
1092 diff_expon
= diff_expon
- ed2
;
1096 __div_128_by_128 (&CQ
, &CR
, CX
, CY
);
1098 if (!CR
.w
[1] && !CR
.w
[0]) {
1099 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
,
1101 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1102 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1106 // get number of decimal digits in CQ
1109 fx
.d
= (float) CQ
.w
[1] * f64
.d
+ (float) CQ
.w
[0];
1110 // binary expon. of CQ
1111 bin_expon
= (fx
.i
- 0x3f800000) >> 23;
1113 digits_q
= estimate_decimal_digits
[bin_expon
];
1114 TP128
.w
[0] = power10_index_binexp_128
[bin_expon
].w
[0];
1115 TP128
.w
[1] = power10_index_binexp_128
[bin_expon
].w
[1];
1116 if (__unsigned_compare_ge_128 (CQ
, TP128
))
1119 ed2
= 34 - digits_q
;
1120 T128
.w
[0] = power10_table_128
[ed2
].w
[0];
1121 T128
.w
[1] = power10_table_128
[ed2
].w
[1];
1122 __mul_128x128_to_256 (CA4
, CR
, T128
);
1123 diff_expon
= diff_expon
- ed2
;
1124 __mul_128x128_low (CQ
, CQ
, T128
);
1128 __div_256_by_128 (&CQ
, &CA4
, CY
);
1130 #ifdef SET_STATUS_FLAGS
1131 if (CA4
.w
[0] || CA4
.w
[1]) {
1133 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
1135 #ifndef LEAVE_TRAILING_ZEROS
1139 #ifndef LEAVE_TRAILING_ZEROS
1140 if (!CA4
.w
[0] && !CA4
.w
[1])
1143 #ifndef LEAVE_TRAILING_ZEROS
1144 // check whether result is exact
1146 //printf("ed2=%d,nz=%d,a=%d,CQ="LX16","LX16", RH="LX16", RL="LX16"\n",ed2,nzeros,amount,CQ.w[1],CQ.w[0],reciprocals10_128[nzeros].w[1],reciprocals10_128[nzeros].w[0]);fflush(stdout);
1147 // check whether CX, CY are short
1148 if (!CX
.w
[1] && !CY
.w
[1] && (CX
.w
[0] <= 1024) && (CY
.w
[0] <= 1024)) {
1149 i
= (int) CY
.w
[0] - 1;
1150 j
= (int) CX
.w
[0] - 1;
1151 // difference in powers of 2 factors for Y and X
1152 nzeros
= ed2
- factors
[i
][0] + factors
[j
][0];
1153 // difference in powers of 5 factors
1154 d5
= ed2
- factors
[i
][1] + factors
[j
][1];
1157 // get P*(2^M[extra_digits])/10^extra_digits
1158 __mul_128x128_high (Qh
, CQ
, reciprocals10_128
[nzeros
]);
1159 //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
1161 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1162 amount
= recip_scale
[nzeros
];
1163 __shr_128_long (CQ
, Qh
, amount
);
1165 diff_expon
+= nzeros
;
1167 // decompose Q as Qh*10^17 + Ql
1168 //T128 = reciprocals10_128[17];
1169 T128
.w
[0] = 0x44909befeb9fad49ull
;
1170 T128
.w
[1] = 0x000b877aa3236a4bull
;
1171 __mul_128x128_to_256 (P256
, CQ
, T128
);
1172 //amount = recip_scale[17];
1173 Q_high
= (P256
.w
[2] >> 44) | (P256
.w
[3] << (64 - 44));
1174 Q_low
= CQ
.w
[0] - Q_high
* 100000000000000000ull;
1179 tdigit
[0] = Q_high
& 0x3ffffff;
1185 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
1187 tdigit
[0] += convert_table
[j
][k
][0];
1188 tdigit
[1] += convert_table
[j
][k
][1];
1189 if (tdigit
[0] >= 100000000) {
1190 tdigit
[0] -= 100000000;
1196 if (tdigit
[1] >= 100000000) {
1197 tdigit
[1] -= 100000000;
1198 if (tdigit
[1] >= 100000000)
1199 tdigit
[1] -= 100000000;
1203 if (!digit
&& !tdigit
[1])
1211 PD
= (UINT64
) digit
*0x068DB8BBull
;
1212 digit_h
= (UINT32
) (PD
>> 40);
1213 //printf("i=%d, nz=%d, digit=%d (%d, %016I64x %016I64x)\n",i,nzeros,digit_h,digit,PD,digit_h);fflush(stdout);
1214 digit_low
= digit
- digit_h
* 10000;
1219 digit_h
= digit_low
;
1223 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
1228 __mul_64x64_to_128 (CQ
, Q_high
, reciprocals10_64
[nzeros
]);
1230 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
1231 amount
= short_recip_scale
[nzeros
];
1232 CQ
.w
[0] = CQ
.w
[1] >> amount
;
1237 diff_expon
+= nzeros
;
1239 tdigit
[0] = Q_low
& 0x3ffffff;
1245 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
1247 tdigit
[0] += convert_table
[j
][k
][0];
1248 tdigit
[1] += convert_table
[j
][k
][1];
1249 if (tdigit
[0] >= 100000000) {
1250 tdigit
[0] -= 100000000;
1255 if (tdigit
[1] >= 100000000) {
1256 tdigit
[1] -= 100000000;
1257 if (tdigit
[1] >= 100000000)
1258 tdigit
[1] -= 100000000;
1262 if (!digit
&& !tdigit
[1])
1270 PD
= (UINT64
) digit
*0x068DB8BBull
;
1271 digit_h
= (UINT32
) (PD
>> 40);
1272 //printf("i=%d, nz=%d, digit=%d (%d, %016I64x %016I64x)\n",i,nzeros,digit_h,digit,PD,digit_h);fflush(stdout);
1273 digit_low
= digit
- digit_h
* 10000;
1278 digit_h
= digit_low
;
1282 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
1287 // get P*(2^M[extra_digits])/10^extra_digits
1288 __mul_128x128_high (Qh
, CQ
, reciprocals10_128
[nzeros
]);
1290 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1291 amount
= recip_scale
[nzeros
];
1292 __shr_128 (CQ
, Qh
, amount
);
1294 diff_expon
+= nzeros
;
1298 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
,
1300 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1301 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1307 if (diff_expon
>= 0) {
1308 #ifdef IEEE_ROUND_NEAREST
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
;
1316 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
1317 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
1320 if (CQ
.w
[0] < carry64
)
1323 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
1326 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
1327 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
1328 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
1329 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
1331 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
1332 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
1335 if (CQ
.w
[0] < carry64
)
1339 if (sign_x
^ sign_y
&& (unsigned) (rmode
- 1) < 2)
1342 case ROUNDING_TO_NEAREST
: // round to nearest code
1345 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
1346 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
1347 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
1348 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
1349 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
1350 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
1352 if (CQ
.w
[0] < carry64
)
1355 case ROUNDING_TIES_AWAY
:
1358 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
1359 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
1360 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
1361 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
1362 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
1363 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
1365 if (CQ
.w
[0] < carry64
)
1369 case ROUNDING_TO_ZERO
:
1371 default: // rounding up
1381 #ifdef SET_STATUS_FLAGS
1382 if (CA4
.w
[0] || CA4
.w
[1]) {
1384 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
1387 handle_UF_128_rem (&res
, sign_x
^ sign_y
, diff_expon
, CQ
,
1388 CA4
.w
[1] | CA4
.w
[0], &rnd_mode
, pfpsf
);
1389 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1390 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1395 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
, pfpsf
);
1396 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1397 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1404 BID128_FUNCTION_ARG128_ARGTYPE2 (bid128qd_div
, x
, UINT64
, y
)
1405 UINT256 CA4
, CA4r
, P256
;
1406 UINT128 CX
, CY
, T128
, CQ
, CR
, CA
, TP128
, Qh
, res
;
1407 UINT64 sign_x
, sign_y
, T
, carry64
, D
, Q_high
, Q_low
, QX
, PD
,
1409 int_float fx
, fy
, f64
;
1410 UINT32 QX32
, tdigit
[3], digit
, digit_h
, digit_low
;
1411 int exponent_x
, exponent_y
, bin_index
, bin_expon
, diff_expon
, ed2
,
1413 int nzeros
, i
, j
, k
, d5
, rmode
;
1414 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1415 fexcept_t binaryflags
= 0;
1419 valid_y
= unpack_BID64 (&sign_y
, &exponent_y
, &CY
.w
[0], y
);
1420 // unpack arguments, check for NaN or Infinity
1421 if (!unpack_BID128_value (&sign_x
, &exponent_x
, &CX
, x
)) {
1423 if ((x
.w
[1] & 0x7c00000000000000ull
) == 0x7c00000000000000ull
) {
1424 #ifdef SET_STATUS_FLAGS
1425 if ((x
.w
[1] & 0x7e00000000000000ull
) == 0x7e00000000000000ull
|| // sNaN
1426 (y
& 0x7e00000000000000ull
) == 0x7e00000000000000ull
)
1427 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
1429 res
.w
[1] = (CX
.w
[1]) & QUIET_MASK64
;
1434 if ((x
.w
[1] & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
1435 // check if y is Inf.
1436 if (((y
& 0x7c00000000000000ull
) == 0x7800000000000000ull
))
1439 #ifdef SET_STATUS_FLAGS
1440 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
1442 res
.w
[1] = 0x7c00000000000000ull
;
1447 if (((y
& 0x7c00000000000000ull
) != 0x7c00000000000000ull
))
1451 res
.w
[1] = ((x
.w
[1] ^ y
) & 0x8000000000000000ull
) |
1452 0x7800000000000000ull
;
1458 if ((y
& 0x7800000000000000ull
) < 0x7800000000000000ull
) {
1460 #ifdef SET_STATUS_FLAGS
1461 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
1463 // x=y=0, return NaN
1464 res
.w
[1] = 0x7c00000000000000ull
;
1469 res
.w
[1] = (x
.w
[1] ^ y
) & 0x8000000000000000ull
;
1470 exponent_x
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS
;
1471 if (exponent_x
> DECIMAL_MAX_EXPON_128
)
1472 exponent_x
= DECIMAL_MAX_EXPON_128
;
1473 else if (exponent_x
< 0)
1475 res
.w
[1] |= (((UINT64
) exponent_x
) << 49);
1485 if ((y
& NAN_MASK64
) == NAN_MASK64
) {
1486 #ifdef SET_STATUS_FLAGS
1487 if ((y
& SNAN_MASK64
) == SNAN_MASK64
) // sNaN
1488 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
1490 res
.w
[0] = (CY
.w
[0] & 0x0003ffffffffffffull
);
1491 __mul_64x64_to_128 (res
, res
.w
[0], power10_table_128
[18].w
[0]);
1492 res
.w
[1] |= ((CY
.w
[0]) & 0xfc00000000000000ull
);
1496 if ((y
& INFINITY_MASK64
) == INFINITY_MASK64
) {
1498 res
.w
[1] = ((x
.w
[1] ^ y
) & 0x8000000000000000ull
);
1503 #ifdef SET_STATUS_FLAGS
1504 __set_status_flags (pfpsf
, ZERO_DIVIDE_EXCEPTION
);
1506 res
.w
[1] = (sign_x
^ sign_y
) | INFINITY_MASK64
;
1510 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1511 (void) fegetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1513 diff_expon
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS
;
1515 if (__unsigned_compare_gt_128 (CY
, CX
)) {
1522 fx
.d
= (float) CX
.w
[1] * f64
.d
+ (float) CX
.w
[0];
1523 fy
.d
= (float) CY
.w
[1] * f64
.d
+ (float) CY
.w
[0];
1524 // expon_cy - expon_cx
1525 bin_index
= (fy
.i
- fx
.i
) >> 23;
1528 T
= power10_index_binexp_128
[bin_index
].w
[0];
1529 __mul_64x128_short (CA
, T
, CX
);
1531 T128
= power10_index_binexp_128
[bin_index
];
1532 __mul_64x128_short (CA
, CX
.w
[0], T128
);
1536 if (__unsigned_compare_gt_128 (CY
, CA
))
1539 T128
= power10_table_128
[ed2
];
1540 __mul_128x128_to_256 (CA4
, CA
, T128
);
1542 ed2
+= estimate_decimal_digits
[bin_index
];
1543 CQ
.w
[0] = CQ
.w
[1] = 0;
1544 diff_expon
= diff_expon
- ed2
;
1548 __div_128_by_128 (&CQ
, &CR
, CX
, CY
);
1550 if (!CR
.w
[1] && !CR
.w
[0]) {
1551 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
,
1553 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1554 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1558 // get number of decimal digits in CQ
1561 fx
.d
= (float) CQ
.w
[1] * f64
.d
+ (float) CQ
.w
[0];
1562 // binary expon. of CQ
1563 bin_expon
= (fx
.i
- 0x3f800000) >> 23;
1565 digits_q
= estimate_decimal_digits
[bin_expon
];
1566 TP128
.w
[0] = power10_index_binexp_128
[bin_expon
].w
[0];
1567 TP128
.w
[1] = power10_index_binexp_128
[bin_expon
].w
[1];
1568 if (__unsigned_compare_ge_128 (CQ
, TP128
))
1571 ed2
= 34 - digits_q
;
1572 T128
.w
[0] = power10_table_128
[ed2
].w
[0];
1573 T128
.w
[1] = power10_table_128
[ed2
].w
[1];
1574 __mul_128x128_to_256 (CA4
, CR
, T128
);
1575 diff_expon
= diff_expon
- ed2
;
1576 __mul_128x128_low (CQ
, CQ
, T128
);
1580 __div_256_by_128 (&CQ
, &CA4
, CY
);
1583 #ifdef SET_STATUS_FLAGS
1584 if (CA4
.w
[0] || CA4
.w
[1]) {
1586 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
1588 #ifndef LEAVE_TRAILING_ZEROS
1592 #ifndef LEAVE_TRAILING_ZEROS
1593 if (!CA4
.w
[0] && !CA4
.w
[1])
1596 #ifndef LEAVE_TRAILING_ZEROS
1597 // check whether result is exact
1599 // check whether CX, CY are short
1600 if (!CX
.w
[1] && !CY
.w
[1] && (CX
.w
[0] <= 1024) && (CY
.w
[0] <= 1024)) {
1601 i
= (int) CY
.w
[0] - 1;
1602 j
= (int) CX
.w
[0] - 1;
1603 // difference in powers of 2 factors for Y and X
1604 nzeros
= ed2
- factors
[i
][0] + factors
[j
][0];
1605 // difference in powers of 5 factors
1606 d5
= ed2
- factors
[i
][1] + factors
[j
][1];
1609 // get P*(2^M[extra_digits])/10^extra_digits
1610 __mul_128x128_high (Qh
, CQ
, reciprocals10_128
[nzeros
]);
1611 //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
1613 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1614 amount
= recip_scale
[nzeros
];
1615 __shr_128_long (CQ
, Qh
, amount
);
1617 diff_expon
+= nzeros
;
1619 // decompose Q as Qh*10^17 + Ql
1620 //T128 = reciprocals10_128[17];
1621 T128
.w
[0] = 0x44909befeb9fad49ull
;
1622 T128
.w
[1] = 0x000b877aa3236a4bull
;
1623 __mul_128x128_to_256 (P256
, CQ
, T128
);
1624 //amount = recip_scale[17];
1625 Q_high
= (P256
.w
[2] >> 44) | (P256
.w
[3] << (64 - 44));
1626 Q_low
= CQ
.w
[0] - Q_high
* 100000000000000000ull;
1631 tdigit
[0] = Q_high
& 0x3ffffff;
1637 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
1639 tdigit
[0] += convert_table
[j
][k
][0];
1640 tdigit
[1] += convert_table
[j
][k
][1];
1641 if (tdigit
[0] >= 100000000) {
1642 tdigit
[0] -= 100000000;
1648 if (tdigit
[1] >= 100000000) {
1649 tdigit
[1] -= 100000000;
1650 if (tdigit
[1] >= 100000000)
1651 tdigit
[1] -= 100000000;
1655 if (!digit
&& !tdigit
[1])
1663 PD
= (UINT64
) digit
*0x068DB8BBull
;
1664 digit_h
= (UINT32
) (PD
>> 40);
1665 digit_low
= digit
- digit_h
* 10000;
1670 digit_h
= digit_low
;
1674 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
1679 __mul_64x64_to_128 (CQ
, Q_high
, reciprocals10_64
[nzeros
]);
1681 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
1682 amount
= short_recip_scale
[nzeros
];
1683 CQ
.w
[0] = CQ
.w
[1] >> amount
;
1688 diff_expon
+= nzeros
;
1690 tdigit
[0] = Q_low
& 0x3ffffff;
1696 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
1698 tdigit
[0] += convert_table
[j
][k
][0];
1699 tdigit
[1] += convert_table
[j
][k
][1];
1700 if (tdigit
[0] >= 100000000) {
1701 tdigit
[0] -= 100000000;
1706 if (tdigit
[1] >= 100000000) {
1707 tdigit
[1] -= 100000000;
1708 if (tdigit
[1] >= 100000000)
1709 tdigit
[1] -= 100000000;
1713 if (!digit
&& !tdigit
[1])
1721 PD
= (UINT64
) digit
*0x068DB8BBull
;
1722 digit_h
= (UINT32
) (PD
>> 40);
1723 digit_low
= digit
- digit_h
* 10000;
1728 digit_h
= digit_low
;
1732 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
1737 // get P*(2^M[extra_digits])/10^extra_digits
1738 __mul_128x128_high (Qh
, CQ
, reciprocals10_128
[nzeros
]);
1740 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1741 amount
= recip_scale
[nzeros
];
1742 __shr_128 (CQ
, Qh
, amount
);
1744 diff_expon
+= nzeros
;
1748 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
,pfpsf
);
1749 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1750 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1756 if (diff_expon
>= 0) {
1757 #ifdef IEEE_ROUND_NEAREST
1760 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
1761 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
1762 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
1763 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
1765 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
1766 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
1769 if (CQ
.w
[0] < carry64
)
1772 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
1775 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
1776 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
1777 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
1778 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
1780 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
1781 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
1784 if (CQ
.w
[0] < carry64
)
1788 if (sign_x
^ sign_y
&& (unsigned) (rmode
- 1) < 2)
1791 case ROUNDING_TO_NEAREST
: // round to nearest code
1794 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
1795 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
1796 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
1797 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
1798 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
1799 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
1801 if (CQ
.w
[0] < carry64
)
1804 case ROUNDING_TIES_AWAY
:
1807 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
1808 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
1809 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
1810 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
1811 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
1812 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
1814 if (CQ
.w
[0] < carry64
)
1818 case ROUNDING_TO_ZERO
:
1820 default: // rounding up
1830 #ifdef SET_STATUS_FLAGS
1831 if (CA4
.w
[0] || CA4
.w
[1]) {
1833 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
1836 handle_UF_128_rem (&res
, sign_x
^ sign_y
, diff_expon
, CQ
,
1837 CA4
.w
[1] | CA4
.w
[0], &rnd_mode
, pfpsf
);
1838 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1839 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1845 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
, pfpsf
);
1846 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1847 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);