1 /* Copyright (C) 2007 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 2, or (at your option) any later
10 In addition to the permissions in the GNU General Public License, the
11 Free Software Foundation gives you unlimited permission to link the
12 compiled version of this file into combinations with other programs,
13 and to distribute those combinations without any restriction coming
14 from the use of this file. (The General Public License restrictions
15 do apply in other respects; for example, they cover modification of
16 the file, and distribution when not linked into a combine
19 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
20 WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
24 You should have received a copy of the GNU General Public License
25 along with GCC; see the file COPYING. If not, write to the Free
26 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
30 #include "bid_div_macros.h"
31 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
34 #define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT
37 extern UINT32 convert_table
[5][128][2];
38 extern SINT8 factors
[][2];
39 extern UINT8 packed_10000_zeros
[];
41 BID128_FUNCTION_ARG2 (bid128_div
, x
, y
)
43 UINT256 CA4
, CA4r
, P256
;
44 UINT128 CX
, CY
, T128
, CQ
, CR
, CA
, TP128
, Qh
, Ql
, res
;
45 UINT64 sign_x
, sign_y
, T
, carry64
, D
, Q_high
, Q_low
, QX
, PD
,
47 int_float fx
, fy
, f64
;
48 UINT32 QX32
, tdigit
[3], digit
, digit_h
, digit_low
;
49 int exponent_x
, exponent_y
, bin_index
, bin_expon
, diff_expon
, ed2
,
51 int nzeros
, i
, j
, k
, d5
;
53 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
54 fexcept_t binaryflags
= 0;
57 valid_y
= unpack_BID128_value (&sign_y
, &exponent_y
, &CY
, y
);
59 // unpack arguments, check for NaN or Infinity
60 if (!unpack_BID128_value (&sign_x
, &exponent_x
, &CX
, x
)) {
62 if ((x
.w
[1] & 0x7c00000000000000ull
) == 0x7c00000000000000ull
) {
63 #ifdef SET_STATUS_FLAGS
64 if ((x
.w
[1] & 0x7e00000000000000ull
) == 0x7e00000000000000ull
|| // sNaN
65 (y
.w
[1] & 0x7e00000000000000ull
) == 0x7e00000000000000ull
)
66 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
68 res
.w
[1] = (CX
.w
[1]) & QUIET_MASK64
;
73 if ((x
.w
[1] & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
75 if (((y
.w
[1] & 0x7c00000000000000ull
) == 0x7800000000000000ull
))
78 #ifdef SET_STATUS_FLAGS
79 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
81 res
.w
[1] = 0x7c00000000000000ull
;
86 if (((y
.w
[1] & 0x7c00000000000000ull
) != 0x7c00000000000000ull
))
90 res
.w
[1] = ((x
.w
[1] ^ y
.w
[1]) & 0x8000000000000000ull
) |
91 0x7800000000000000ull
;
97 if ((y
.w
[1] & 0x7800000000000000ull
) < 0x7800000000000000ull
) {
98 if ((!CY
.w
[0]) && !(CY
.w
[1] & 0x0001ffffffffffffull
)) {
99 #ifdef SET_STATUS_FLAGS
100 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
103 res
.w
[1] = 0x7c00000000000000ull
;
108 res
.w
[1] = (x
.w
[1] ^ y
.w
[1]) & 0x8000000000000000ull
;
109 exponent_x
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS_128
;
110 if (exponent_x
> DECIMAL_MAX_EXPON_128
)
111 exponent_x
= DECIMAL_MAX_EXPON_128
;
112 else if (exponent_x
< 0)
114 res
.w
[1] |= (((UINT64
) exponent_x
) << 49);
123 if ((y
.w
[1] & 0x7c00000000000000ull
) == 0x7c00000000000000ull
) {
124 #ifdef SET_STATUS_FLAGS
125 if ((y
.w
[1] & 0x7e00000000000000ull
) == 0x7e00000000000000ull
) // sNaN
126 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
128 res
.w
[1] = CY
.w
[1] & QUIET_MASK64
;
133 if ((y
.w
[1] & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
135 res
.w
[1] = sign_x
^ sign_y
;
139 // y is 0, return +/-Inf
140 #ifdef SET_STATUS_FLAGS
141 __set_status_flags (pfpsf
, ZERO_DIVIDE_EXCEPTION
);
144 ((x
.w
[1] ^ y
.w
[1]) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
148 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
149 (void) fegetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
151 diff_expon
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS_128
;
153 if (__unsigned_compare_gt_128 (CY
, CX
)) {
160 fx
.d
= (float) CX
.w
[1] * f64
.d
+ (float) CX
.w
[0];
161 fy
.d
= (float) CY
.w
[1] * f64
.d
+ (float) CY
.w
[0];
162 // expon_cy - expon_cx
163 bin_index
= (fy
.i
- fx
.i
) >> 23;
166 T
= power10_index_binexp_128
[bin_index
].w
[0];
167 __mul_64x128_short (CA
, T
, CX
);
169 T128
= power10_index_binexp_128
[bin_index
];
170 __mul_64x128_short (CA
, CX
.w
[0], T128
);
174 if (__unsigned_compare_gt_128 (CY
, CA
))
177 T128
= power10_table_128
[ed2
];
178 __mul_128x128_to_256 (CA4
, CA
, T128
);
180 ed2
+= estimate_decimal_digits
[bin_index
];
181 CQ
.w
[0] = CQ
.w
[1] = 0;
182 diff_expon
= diff_expon
- ed2
;
186 __div_128_by_128 (&CQ
, &CR
, CX
, CY
);
188 if (!CR
.w
[1] && !CR
.w
[0]) {
189 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
,
191 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
192 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
196 // get number of decimal digits in CQ
199 fx
.d
= (float) CQ
.w
[1] * f64
.d
+ (float) CQ
.w
[0];
200 // binary expon. of CQ
201 bin_expon
= (fx
.i
- 0x3f800000) >> 23;
203 digits_q
= estimate_decimal_digits
[bin_expon
];
204 TP128
.w
[0] = power10_index_binexp_128
[bin_expon
].w
[0];
205 TP128
.w
[1] = power10_index_binexp_128
[bin_expon
].w
[1];
206 if (__unsigned_compare_ge_128 (CQ
, TP128
))
210 T128
.w
[0] = power10_table_128
[ed2
].w
[0];
211 T128
.w
[1] = power10_table_128
[ed2
].w
[1];
212 __mul_128x128_to_256 (CA4
, CR
, T128
);
213 diff_expon
= diff_expon
- ed2
;
214 __mul_128x128_low (CQ
, CQ
, T128
);
218 __div_256_by_128 (&CQ
, &CA4
, CY
);
220 #ifdef SET_STATUS_FLAGS
221 if (CA4
.w
[0] || CA4
.w
[1]) {
223 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
225 #ifndef LEAVE_TRAILING_ZEROS
229 #ifndef LEAVE_TRAILING_ZEROS
230 if (!CA4
.w
[0] && !CA4
.w
[1])
233 #ifndef LEAVE_TRAILING_ZEROS
234 // check whether result is exact
236 // check whether CX, CY are short
237 if (!CX
.w
[1] && !CY
.w
[1] && (CX
.w
[0] <= 1024) && (CY
.w
[0] <= 1024)) {
238 i
= (int) CY
.w
[0] - 1;
239 j
= (int) CX
.w
[0] - 1;
240 // difference in powers of 2 factors for Y and X
241 nzeros
= ed2
- factors
[i
][0] + factors
[j
][0];
242 // difference in powers of 5 factors
243 d5
= ed2
- factors
[i
][1] + factors
[j
][1];
246 // get P*(2^M[extra_digits])/10^extra_digits
247 __mul_128x128_full (Qh
, Ql
, CQ
, reciprocals10_128
[nzeros
]);
249 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
250 amount
= recip_scale
[nzeros
];
251 __shr_128_long (CQ
, Qh
, amount
);
253 diff_expon
+= nzeros
;
255 // decompose Q as Qh*10^17 + Ql
256 //T128 = reciprocals10_128[17];
257 T128
.w
[0] = 0x44909befeb9fad49ull
;
258 T128
.w
[1] = 0x000b877aa3236a4bull
;
259 __mul_128x128_to_256 (P256
, CQ
, T128
);
260 //amount = recip_scale[17];
261 Q_high
= (P256
.w
[2] >> 44) | (P256
.w
[3] << (64 - 44));
262 Q_low
= CQ
.w
[0] - Q_high
* 100000000000000000ull;
267 tdigit
[0] = Q_high
& 0x3ffffff;
273 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
275 tdigit
[0] += convert_table
[j
][k
][0];
276 tdigit
[1] += convert_table
[j
][k
][1];
277 if (tdigit
[0] >= 100000000) {
278 tdigit
[0] -= 100000000;
283 if (tdigit
[1] >= 100000000) {
284 tdigit
[1] -= 100000000;
285 if (tdigit
[1] >= 100000000)
286 tdigit
[1] -= 100000000;
290 if (!digit
&& !tdigit
[1])
298 PD
= (UINT64
) digit
*0x068DB8BBull
;
299 digit_h
= (UINT32
) (PD
>> 40);
300 digit_low
= digit
- digit_h
* 10000;
309 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
314 __mul_64x64_to_128 (CQ
, Q_high
, reciprocals10_64
[nzeros
]);
316 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
317 amount
= short_recip_scale
[nzeros
];
318 CQ
.w
[0] = CQ
.w
[1] >> amount
;
323 diff_expon
+= nzeros
;
325 tdigit
[0] = Q_low
& 0x3ffffff;
331 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
333 tdigit
[0] += convert_table
[j
][k
][0];
334 tdigit
[1] += convert_table
[j
][k
][1];
335 if (tdigit
[0] >= 100000000) {
336 tdigit
[0] -= 100000000;
341 if (tdigit
[1] >= 100000000) {
342 tdigit
[1] -= 100000000;
343 if (tdigit
[1] >= 100000000)
344 tdigit
[1] -= 100000000;
348 if (!digit
&& !tdigit
[1])
356 PD
= (UINT64
) digit
*0x068DB8BBull
;
357 digit_h
= (UINT32
) (PD
>> 40);
358 digit_low
= digit
- digit_h
* 10000;
367 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
372 // get P*(2^M[extra_digits])/10^extra_digits
373 __mul_128x128_full (Qh
, Ql
, CQ
, reciprocals10_128
[nzeros
]);
375 //now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
376 amount
= recip_scale
[nzeros
];
377 __shr_128 (CQ
, Qh
, amount
);
379 diff_expon
+= nzeros
;
383 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
, pfpsf
);
384 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
385 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
391 if (diff_expon
>= 0) {
392 #ifdef IEEE_ROUND_NEAREST
395 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
396 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
397 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
398 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
400 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
401 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
404 if (CQ
.w
[0] < carry64
)
407 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
410 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
411 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
412 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
413 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
415 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
416 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
419 if (CQ
.w
[0] < carry64
)
423 if (sign_x
^ sign_y
&& (unsigned) (rmode
- 1) < 2)
426 case ROUNDING_TO_NEAREST
: // round to nearest code
429 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
430 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
431 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
432 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
433 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
434 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
436 if (CQ
.w
[0] < carry64
)
439 case ROUNDING_TIES_AWAY
:
442 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
443 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
444 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
445 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
446 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
447 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
449 if (CQ
.w
[0] < carry64
)
453 case ROUNDING_TO_ZERO
:
455 default: // rounding up
465 #ifdef SET_STATUS_FLAGS
466 if (CA4
.w
[0] || CA4
.w
[1]) {
468 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
472 handle_UF_128_rem (&res
, sign_x
^ sign_y
, diff_expon
, CQ
,
473 CA4
.w
[1] | CA4
.w
[0], &rnd_mode
, pfpsf
);
474 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
475 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
481 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
, pfpsf
);
482 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
483 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
489 //#define LEAVE_TRAILING_ZEROS
491 TYPE0_FUNCTION_ARGTYPE1_ARGTYPE2 (UINT128
, bid128dd_div
, UINT64
, x
,
494 UINT256 CA4
, CA4r
, P256
;
495 UINT128 CX
, CY
, T128
, CQ
, CR
, CA
, TP128
, Qh
, Ql
, res
;
496 UINT64 sign_x
, sign_y
, T
, carry64
, D
, Q_high
, Q_low
, QX
, PD
,
498 int_float fx
, fy
, f64
;
499 UINT32 QX32
, tdigit
[3], digit
, digit_h
, digit_low
;
500 int exponent_x
, exponent_y
, bin_index
, bin_expon
, diff_expon
, ed2
,
502 int nzeros
, i
, j
, k
, d5
;
504 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
505 fexcept_t binaryflags
= 0;
508 valid_y
= unpack_BID64 (&sign_y
, &exponent_y
, &CY
.w
[0], y
);
510 // unpack arguments, check for NaN or Infinity
512 if (!unpack_BID64 (&sign_x
, &exponent_x
, &CX
.w
[0], (x
))) {
513 #ifdef SET_STATUS_FLAGS
514 if ((y
& SNAN_MASK64
) == SNAN_MASK64
) // y is sNaN
515 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
519 if ((x
& NAN_MASK64
) == NAN_MASK64
) {
520 #ifdef SET_STATUS_FLAGS
521 if ((x
& SNAN_MASK64
) == SNAN_MASK64
) // sNaN
522 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
524 res
.w
[0] = (CX
.w
[0] & 0x0003ffffffffffffull
);
525 __mul_64x64_to_128 (res
, res
.w
[0], power10_table_128
[18].w
[0]);
526 res
.w
[1] |= ((CX
.w
[0]) & 0xfc00000000000000ull
);
530 if (((x
) & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
531 // check if y is Inf.
532 if ((((y
) & 0x7c00000000000000ull
) == 0x7800000000000000ull
))
535 #ifdef SET_STATUS_FLAGS
536 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
538 res
.w
[1] = 0x7c00000000000000ull
;
542 if ((((y
) & 0x7c00000000000000ull
) != 0x7c00000000000000ull
)) {
543 // otherwise return +/-Inf
545 (((x
) ^ (y
)) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
551 if ((((y
) & 0x7800000000000000ull
) != 0x7800000000000000ull
)) {
553 #ifdef SET_STATUS_FLAGS
554 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
557 res
.w
[1] = 0x7c00000000000000ull
;
562 res
.w
[1] = ((x
) ^ (y
)) & 0x8000000000000000ull
;
563 if (((y
) & 0x6000000000000000ull
) == 0x6000000000000000ull
)
564 exponent_y
= ((UINT32
) ((y
) >> 51)) & 0x3ff;
566 exponent_y
= ((UINT32
) ((y
) >> 53)) & 0x3ff;
567 exponent_x
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS_128
;
568 if (exponent_x
> DECIMAL_MAX_EXPON_128
)
569 exponent_x
= DECIMAL_MAX_EXPON_128
;
570 else if (exponent_x
< 0)
572 res
.w
[1] |= (((UINT64
) exponent_x
) << 49);
583 if ((y
& NAN_MASK64
) == NAN_MASK64
) {
584 #ifdef SET_STATUS_FLAGS
585 if ((y
& SNAN_MASK64
) == SNAN_MASK64
) // sNaN
586 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
588 res
.w
[0] = (CY
.w
[0] & 0x0003ffffffffffffull
);
589 __mul_64x64_to_128 (res
, res
.w
[0], power10_table_128
[18].w
[0]);
590 res
.w
[1] |= ((CY
.w
[0]) & 0xfc00000000000000ull
);
594 if (((y
) & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
596 res
.w
[1] = sign_x
^ sign_y
;
600 // y is 0, return +/-Inf
602 (((x
) ^ (y
)) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
604 #ifdef SET_STATUS_FLAGS
605 __set_status_flags (pfpsf
, ZERO_DIVIDE_EXCEPTION
);
609 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
610 (void) fegetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
612 diff_expon
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS_128
;
614 if (__unsigned_compare_gt_128 (CY
, CX
)) {
621 fx
.d
= (float) CX
.w
[1] * f64
.d
+ (float) CX
.w
[0];
622 fy
.d
= (float) CY
.w
[1] * f64
.d
+ (float) CY
.w
[0];
623 // expon_cy - expon_cx
624 bin_index
= (fy
.i
- fx
.i
) >> 23;
627 T
= power10_index_binexp_128
[bin_index
].w
[0];
628 __mul_64x128_short (CA
, T
, CX
);
630 T128
= power10_index_binexp_128
[bin_index
];
631 __mul_64x128_short (CA
, CX
.w
[0], T128
);
635 if (__unsigned_compare_gt_128 (CY
, CA
))
638 T128
= power10_table_128
[ed2
];
639 __mul_128x128_to_256 (CA4
, CA
, T128
);
641 ed2
+= estimate_decimal_digits
[bin_index
];
642 CQ
.w
[0] = CQ
.w
[1] = 0;
643 diff_expon
= diff_expon
- ed2
;
647 __div_128_by_128 (&CQ
, &CR
, CX
, CY
);
649 if (!CR
.w
[1] && !CR
.w
[0]) {
650 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
,
652 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
653 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
657 // get number of decimal digits in CQ
660 fx
.d
= (float) CQ
.w
[1] * f64
.d
+ (float) CQ
.w
[0];
661 // binary expon. of CQ
662 bin_expon
= (fx
.i
- 0x3f800000) >> 23;
664 digits_q
= estimate_decimal_digits
[bin_expon
];
665 TP128
.w
[0] = power10_index_binexp_128
[bin_expon
].w
[0];
666 TP128
.w
[1] = power10_index_binexp_128
[bin_expon
].w
[1];
667 if (__unsigned_compare_ge_128 (CQ
, TP128
))
671 T128
.w
[0] = power10_table_128
[ed2
].w
[0];
672 T128
.w
[1] = power10_table_128
[ed2
].w
[1];
673 __mul_128x128_to_256 (CA4
, CR
, T128
);
674 diff_expon
= diff_expon
- ed2
;
675 __mul_128x128_low (CQ
, CQ
, T128
);
679 __div_256_by_128 (&CQ
, &CA4
, CY
);
682 #ifdef SET_STATUS_FLAGS
683 if (CA4
.w
[0] || CA4
.w
[1]) {
685 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
687 #ifndef LEAVE_TRAILING_ZEROS
691 #ifndef LEAVE_TRAILING_ZEROS
692 if (!CA4
.w
[0] && !CA4
.w
[1])
695 #ifndef LEAVE_TRAILING_ZEROS
696 // check whether result is exact
698 // check whether CX, CY are short
699 if (!CX
.w
[1] && !CY
.w
[1] && (CX
.w
[0] <= 1024) && (CY
.w
[0] <= 1024)) {
700 i
= (int) CY
.w
[0] - 1;
701 j
= (int) CX
.w
[0] - 1;
702 // difference in powers of 2 factors for Y and X
703 nzeros
= ed2
- factors
[i
][0] + factors
[j
][0];
704 // difference in powers of 5 factors
705 d5
= ed2
- factors
[i
][1] + factors
[j
][1];
708 // get P*(2^M[extra_digits])/10^extra_digits
709 __mul_128x128_full (Qh
, Ql
, CQ
, reciprocals10_128
[nzeros
]);
710 //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
712 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
713 amount
= recip_scale
[nzeros
];
714 __shr_128_long (CQ
, Qh
, amount
);
716 diff_expon
+= nzeros
;
718 // decompose Q as Qh*10^17 + Ql
719 //T128 = reciprocals10_128[17];
720 T128
.w
[0] = 0x44909befeb9fad49ull
;
721 T128
.w
[1] = 0x000b877aa3236a4bull
;
722 __mul_128x128_to_256 (P256
, CQ
, T128
);
723 //amount = recip_scale[17];
724 Q_high
= (P256
.w
[2] >> 44) | (P256
.w
[3] << (64 - 44));
725 Q_low
= CQ
.w
[0] - Q_high
* 100000000000000000ull;
730 tdigit
[0] = Q_high
& 0x3ffffff;
736 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
738 tdigit
[0] += convert_table
[j
][k
][0];
739 tdigit
[1] += convert_table
[j
][k
][1];
740 if (tdigit
[0] >= 100000000) {
741 tdigit
[0] -= 100000000;
747 if (tdigit
[1] >= 100000000) {
748 tdigit
[1] -= 100000000;
749 if (tdigit
[1] >= 100000000)
750 tdigit
[1] -= 100000000;
754 if (!digit
&& !tdigit
[1])
762 PD
= (UINT64
) digit
*0x068DB8BBull
;
763 digit_h
= (UINT32
) (PD
>> 40);
764 digit_low
= digit
- digit_h
* 10000;
773 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
778 __mul_64x64_to_128 (CQ
, Q_high
, reciprocals10_64
[nzeros
]);
780 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
781 amount
= short_recip_scale
[nzeros
];
782 CQ
.w
[0] = CQ
.w
[1] >> amount
;
787 diff_expon
+= nzeros
;
789 tdigit
[0] = Q_low
& 0x3ffffff;
795 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
797 tdigit
[0] += convert_table
[j
][k
][0];
798 tdigit
[1] += convert_table
[j
][k
][1];
799 if (tdigit
[0] >= 100000000) {
800 tdigit
[0] -= 100000000;
805 if (tdigit
[1] >= 100000000) {
806 tdigit
[1] -= 100000000;
807 if (tdigit
[1] >= 100000000)
808 tdigit
[1] -= 100000000;
812 if (!digit
&& !tdigit
[1])
820 PD
= (UINT64
) digit
*0x068DB8BBull
;
821 digit_h
= (UINT32
) (PD
>> 40);
822 digit_low
= digit
- digit_h
* 10000;
831 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
836 // get P*(2^M[extra_digits])/10^extra_digits
837 __mul_128x128_full (Qh
, Ql
, CQ
, reciprocals10_128
[nzeros
]);
839 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
840 amount
= recip_scale
[nzeros
];
841 __shr_128 (CQ
, Qh
, amount
);
843 diff_expon
+= nzeros
;
847 get_BID128(&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
,pfpsf
);
848 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
849 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
855 if (diff_expon
>= 0) {
856 #ifdef IEEE_ROUND_NEAREST
859 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
860 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
861 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
862 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
864 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
865 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
868 if (CQ
.w
[0] < carry64
)
871 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
874 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
875 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
876 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
877 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
879 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
880 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
883 if (CQ
.w
[0] < carry64
)
887 if (sign_x
^ sign_y
&& (unsigned) (rmode
- 1) < 2)
890 case ROUNDING_TO_NEAREST
: // round to nearest code
893 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
894 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
895 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
896 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
897 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
898 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
900 if (CQ
.w
[0] < carry64
)
903 case ROUNDING_TIES_AWAY
:
906 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
907 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
908 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
909 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
910 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
911 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
913 if (CQ
.w
[0] < carry64
)
917 case ROUNDING_TO_ZERO
:
919 default: // rounding up
929 #ifdef SET_STATUS_FLAGS
930 if (CA4
.w
[0] || CA4
.w
[1]) {
932 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
935 handle_UF_128_rem (&res
, sign_x
^ sign_y
, diff_expon
, CQ
,
936 CA4
.w
[1] | CA4
.w
[0], &rnd_mode
, pfpsf
);
937 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
938 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
944 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
, pfpsf
);
945 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
946 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
952 BID128_FUNCTION_ARGTYPE1_ARG128 (bid128dq_div
, UINT64
, x
, y
)
953 UINT256 CA4
, CA4r
, P256
;
954 UINT128 CX
, CY
, T128
, CQ
, CR
, CA
, TP128
, Qh
, Ql
, res
;
955 UINT64 sign_x
, sign_y
, T
, carry64
, D
, Q_high
, Q_low
, QX
, valid_y
,
957 int_float fx
, fy
, f64
;
958 UINT32 QX32
, tdigit
[3], digit
, digit_h
, digit_low
;
959 int exponent_x
, exponent_y
, bin_index
, bin_expon
, diff_expon
, ed2
,
961 int nzeros
, i
, j
, k
, d5
;
963 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
964 fexcept_t binaryflags
= 0;
967 valid_y
= unpack_BID128_value (&sign_y
, &exponent_y
, &CY
, y
);
969 // unpack arguments, check for NaN or Infinity
971 if (!unpack_BID64 (&sign_x
, &exponent_x
, &CX
.w
[0], x
)) {
972 #ifdef SET_STATUS_FLAGS
973 if ((y
.w
[1] & SNAN_MASK64
) == SNAN_MASK64
) // y is sNaN
974 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
978 if ((x
& NAN_MASK64
) == NAN_MASK64
) {
979 #ifdef SET_STATUS_FLAGS
980 if ((x
& SNAN_MASK64
) == SNAN_MASK64
) // sNaN
981 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
983 res
.w
[0] = (CX
.w
[0] & 0x0003ffffffffffffull
);
984 __mul_64x64_to_128 (res
, res
.w
[0], power10_table_128
[18].w
[0]);
985 res
.w
[1] |= ((CX
.w
[0]) & 0xfc00000000000000ull
);
989 if ((x
& 0x7800000000000000ull
) == 0x7800000000000000ull
) {
990 // check if y is Inf.
991 if (((y
.w
[1] & 0x7c00000000000000ull
) == 0x7800000000000000ull
))
994 #ifdef SET_STATUS_FLAGS
995 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
997 res
.w
[1] = 0x7c00000000000000ull
;
1001 if (((y
.w
[1] & 0x7c00000000000000ull
) != 0x7c00000000000000ull
)) {
1002 // otherwise return +/-Inf
1004 ((x
^ y
.w
[1]) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
1010 if ((y
.w
[1] & INFINITY_MASK64
) != INFINITY_MASK64
) {
1011 if ((!CY
.w
[0]) && !(CY
.w
[1] & 0x0001ffffffffffffull
)) {
1012 #ifdef SET_STATUS_FLAGS
1013 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
1015 // x=y=0, return NaN
1016 res
.w
[1] = 0x7c00000000000000ull
;
1021 res
.w
[1] = (x
^ y
.w
[1]) & 0x8000000000000000ull
;
1022 exponent_x
= exponent_x
- exponent_y
+ (DECIMAL_EXPONENT_BIAS_128
<<1) - DECIMAL_EXPONENT_BIAS
;
1023 if (exponent_x
> DECIMAL_MAX_EXPON_128
)
1024 exponent_x
= DECIMAL_MAX_EXPON_128
;
1025 else if (exponent_x
< 0)
1027 res
.w
[1] |= (((UINT64
) exponent_x
) << 49);
1032 exponent_x
+= (DECIMAL_EXPONENT_BIAS_128
- DECIMAL_EXPONENT_BIAS
);
1038 if ((y
.w
[1] & 0x7c00000000000000ull
) == 0x7c00000000000000ull
) {
1039 #ifdef SET_STATUS_FLAGS
1040 if ((y
.w
[1] & 0x7e00000000000000ull
) == 0x7e00000000000000ull
) // sNaN
1041 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
1043 res
.w
[1] = CY
.w
[1] & QUIET_MASK64
;
1048 if ((y
.w
[1] & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
1050 res
.w
[1] = sign_x
^ sign_y
;
1054 // y is 0, return +/-Inf
1056 ((x
^ y
.w
[1]) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
1058 #ifdef SET_STATUS_FLAGS
1059 __set_status_flags (pfpsf
, ZERO_DIVIDE_EXCEPTION
);
1063 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1064 (void) fegetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1066 diff_expon
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS_128
;
1068 if (__unsigned_compare_gt_128 (CY
, CX
)) {
1075 fx
.d
= (float) CX
.w
[1] * f64
.d
+ (float) CX
.w
[0];
1076 fy
.d
= (float) CY
.w
[1] * f64
.d
+ (float) CY
.w
[0];
1077 // expon_cy - expon_cx
1078 bin_index
= (fy
.i
- fx
.i
) >> 23;
1081 T
= power10_index_binexp_128
[bin_index
].w
[0];
1082 __mul_64x128_short (CA
, T
, CX
);
1084 T128
= power10_index_binexp_128
[bin_index
];
1085 __mul_64x128_short (CA
, CX
.w
[0], T128
);
1089 if (__unsigned_compare_gt_128 (CY
, CA
))
1092 T128
= power10_table_128
[ed2
];
1093 __mul_128x128_to_256 (CA4
, CA
, T128
);
1095 ed2
+= estimate_decimal_digits
[bin_index
];
1096 CQ
.w
[0] = CQ
.w
[1] = 0;
1097 diff_expon
= diff_expon
- ed2
;
1101 __div_128_by_128 (&CQ
, &CR
, CX
, CY
);
1103 if (!CR
.w
[1] && !CR
.w
[0]) {
1104 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
,
1106 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1107 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1111 // get number of decimal digits in CQ
1114 fx
.d
= (float) CQ
.w
[1] * f64
.d
+ (float) CQ
.w
[0];
1115 // binary expon. of CQ
1116 bin_expon
= (fx
.i
- 0x3f800000) >> 23;
1118 digits_q
= estimate_decimal_digits
[bin_expon
];
1119 TP128
.w
[0] = power10_index_binexp_128
[bin_expon
].w
[0];
1120 TP128
.w
[1] = power10_index_binexp_128
[bin_expon
].w
[1];
1121 if (__unsigned_compare_ge_128 (CQ
, TP128
))
1124 ed2
= 34 - digits_q
;
1125 T128
.w
[0] = power10_table_128
[ed2
].w
[0];
1126 T128
.w
[1] = power10_table_128
[ed2
].w
[1];
1127 __mul_128x128_to_256 (CA4
, CR
, T128
);
1128 diff_expon
= diff_expon
- ed2
;
1129 __mul_128x128_low (CQ
, CQ
, T128
);
1133 __div_256_by_128 (&CQ
, &CA4
, CY
);
1135 #ifdef SET_STATUS_FLAGS
1136 if (CA4
.w
[0] || CA4
.w
[1]) {
1138 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
1140 #ifndef LEAVE_TRAILING_ZEROS
1144 #ifndef LEAVE_TRAILING_ZEROS
1145 if (!CA4
.w
[0] && !CA4
.w
[1])
1148 #ifndef LEAVE_TRAILING_ZEROS
1149 // check whether result is exact
1151 //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);
1152 // check whether CX, CY are short
1153 if (!CX
.w
[1] && !CY
.w
[1] && (CX
.w
[0] <= 1024) && (CY
.w
[0] <= 1024)) {
1154 i
= (int) CY
.w
[0] - 1;
1155 j
= (int) CX
.w
[0] - 1;
1156 // difference in powers of 2 factors for Y and X
1157 nzeros
= ed2
- factors
[i
][0] + factors
[j
][0];
1158 // difference in powers of 5 factors
1159 d5
= ed2
- factors
[i
][1] + factors
[j
][1];
1162 // get P*(2^M[extra_digits])/10^extra_digits
1163 __mul_128x128_full (Qh
, Ql
, CQ
, reciprocals10_128
[nzeros
]);
1164 //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
1166 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1167 amount
= recip_scale
[nzeros
];
1168 __shr_128_long (CQ
, Qh
, amount
);
1170 diff_expon
+= nzeros
;
1172 // decompose Q as Qh*10^17 + Ql
1173 //T128 = reciprocals10_128[17];
1174 T128
.w
[0] = 0x44909befeb9fad49ull
;
1175 T128
.w
[1] = 0x000b877aa3236a4bull
;
1176 __mul_128x128_to_256 (P256
, CQ
, T128
);
1177 //amount = recip_scale[17];
1178 Q_high
= (P256
.w
[2] >> 44) | (P256
.w
[3] << (64 - 44));
1179 Q_low
= CQ
.w
[0] - Q_high
* 100000000000000000ull;
1184 tdigit
[0] = Q_high
& 0x3ffffff;
1190 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
1192 tdigit
[0] += convert_table
[j
][k
][0];
1193 tdigit
[1] += convert_table
[j
][k
][1];
1194 if (tdigit
[0] >= 100000000) {
1195 tdigit
[0] -= 100000000;
1201 if (tdigit
[1] >= 100000000) {
1202 tdigit
[1] -= 100000000;
1203 if (tdigit
[1] >= 100000000)
1204 tdigit
[1] -= 100000000;
1208 if (!digit
&& !tdigit
[1])
1216 PD
= (UINT64
) digit
*0x068DB8BBull
;
1217 digit_h
= (UINT32
) (PD
>> 40);
1218 //printf("i=%d, nz=%d, digit=%d (%d, %016I64x %016I64x)\n",i,nzeros,digit_h,digit,PD,digit_h);fflush(stdout);
1219 digit_low
= digit
- digit_h
* 10000;
1224 digit_h
= digit_low
;
1228 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
1233 __mul_64x64_to_128 (CQ
, Q_high
, reciprocals10_64
[nzeros
]);
1235 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
1236 amount
= short_recip_scale
[nzeros
];
1237 CQ
.w
[0] = CQ
.w
[1] >> amount
;
1242 diff_expon
+= nzeros
;
1244 tdigit
[0] = Q_low
& 0x3ffffff;
1250 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
1252 tdigit
[0] += convert_table
[j
][k
][0];
1253 tdigit
[1] += convert_table
[j
][k
][1];
1254 if (tdigit
[0] >= 100000000) {
1255 tdigit
[0] -= 100000000;
1260 if (tdigit
[1] >= 100000000) {
1261 tdigit
[1] -= 100000000;
1262 if (tdigit
[1] >= 100000000)
1263 tdigit
[1] -= 100000000;
1267 if (!digit
&& !tdigit
[1])
1275 PD
= (UINT64
) digit
*0x068DB8BBull
;
1276 digit_h
= (UINT32
) (PD
>> 40);
1277 //printf("i=%d, nz=%d, digit=%d (%d, %016I64x %016I64x)\n",i,nzeros,digit_h,digit,PD,digit_h);fflush(stdout);
1278 digit_low
= digit
- digit_h
* 10000;
1283 digit_h
= digit_low
;
1287 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
1292 // get P*(2^M[extra_digits])/10^extra_digits
1293 __mul_128x128_full (Qh
, Ql
, CQ
, reciprocals10_128
[nzeros
]);
1295 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1296 amount
= recip_scale
[nzeros
];
1297 __shr_128 (CQ
, Qh
, amount
);
1299 diff_expon
+= nzeros
;
1303 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
,
1305 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1306 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1312 if (diff_expon
>= 0) {
1313 #ifdef IEEE_ROUND_NEAREST
1316 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
1317 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
1318 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
1319 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
1321 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
1322 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
1325 if (CQ
.w
[0] < carry64
)
1328 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
1331 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
1332 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
1333 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
1334 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
1336 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
1337 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
1340 if (CQ
.w
[0] < carry64
)
1344 if (sign_x
^ sign_y
&& (unsigned) (rmode
- 1) < 2)
1347 case ROUNDING_TO_NEAREST
: // round to nearest code
1350 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
1351 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
1352 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
1353 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
1354 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
1355 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
1357 if (CQ
.w
[0] < carry64
)
1360 case ROUNDING_TIES_AWAY
:
1363 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
1364 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
1365 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
1366 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
1367 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
1368 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
1370 if (CQ
.w
[0] < carry64
)
1374 case ROUNDING_TO_ZERO
:
1376 default: // rounding up
1386 #ifdef SET_STATUS_FLAGS
1387 if (CA4
.w
[0] || CA4
.w
[1]) {
1389 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
1392 handle_UF_128_rem (&res
, sign_x
^ sign_y
, diff_expon
, CQ
,
1393 CA4
.w
[1] | CA4
.w
[0], &rnd_mode
, pfpsf
);
1394 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1395 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1400 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
, pfpsf
);
1401 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1402 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1409 BID128_FUNCTION_ARG128_ARGTYPE2 (bid128qd_div
, x
, UINT64
, y
)
1410 UINT256 CA4
, CA4r
, P256
;
1411 UINT128 CX
, CY
, T128
, CQ
, CR
, CA
, TP128
, Qh
, Ql
, res
;
1412 UINT64 sign_x
, sign_y
, T
, carry64
, D
, Q_high
, Q_low
, QX
, PD
,
1414 int_float fx
, fy
, f64
;
1415 UINT32 QX32
, tdigit
[3], digit
, digit_h
, digit_low
;
1416 int exponent_x
, exponent_y
, bin_index
, bin_expon
, diff_expon
, ed2
,
1418 int nzeros
, i
, j
, k
, d5
, rmode
;
1419 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1420 fexcept_t binaryflags
= 0;
1424 valid_y
= unpack_BID64 (&sign_y
, &exponent_y
, &CY
.w
[0], y
);
1425 // unpack arguments, check for NaN or Infinity
1426 if (!unpack_BID128_value (&sign_x
, &exponent_x
, &CX
, x
)) {
1428 if ((x
.w
[1] & 0x7c00000000000000ull
) == 0x7c00000000000000ull
) {
1429 #ifdef SET_STATUS_FLAGS
1430 if ((x
.w
[1] & 0x7e00000000000000ull
) == 0x7e00000000000000ull
|| // sNaN
1431 (y
& 0x7e00000000000000ull
) == 0x7e00000000000000ull
)
1432 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
1434 res
.w
[1] = (CX
.w
[1]) & QUIET_MASK64
;
1439 if ((x
.w
[1] & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
1440 // check if y is Inf.
1441 if (((y
& 0x7c00000000000000ull
) == 0x7800000000000000ull
))
1444 #ifdef SET_STATUS_FLAGS
1445 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
1447 res
.w
[1] = 0x7c00000000000000ull
;
1452 if (((y
& 0x7c00000000000000ull
) != 0x7c00000000000000ull
))
1456 res
.w
[1] = ((x
.w
[1] ^ y
) & 0x8000000000000000ull
) |
1457 0x7800000000000000ull
;
1463 if ((y
& 0x7800000000000000ull
) < 0x7800000000000000ull
) {
1465 #ifdef SET_STATUS_FLAGS
1466 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
1468 // x=y=0, return NaN
1469 res
.w
[1] = 0x7c00000000000000ull
;
1474 res
.w
[1] = (x
.w
[1] ^ y
) & 0x8000000000000000ull
;
1475 exponent_x
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS
;
1476 if (exponent_x
> DECIMAL_MAX_EXPON_128
)
1477 exponent_x
= DECIMAL_MAX_EXPON_128
;
1478 else if (exponent_x
< 0)
1480 res
.w
[1] |= (((UINT64
) exponent_x
) << 49);
1490 if ((y
& NAN_MASK64
) == NAN_MASK64
) {
1491 #ifdef SET_STATUS_FLAGS
1492 if ((y
& SNAN_MASK64
) == SNAN_MASK64
) // sNaN
1493 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
1495 res
.w
[0] = (CY
.w
[0] & 0x0003ffffffffffffull
);
1496 __mul_64x64_to_128 (res
, res
.w
[0], power10_table_128
[18].w
[0]);
1497 res
.w
[1] |= ((CY
.w
[0]) & 0xfc00000000000000ull
);
1501 if ((y
& INFINITY_MASK64
) == INFINITY_MASK64
) {
1503 res
.w
[1] = ((x
.w
[1] ^ y
) & 0x8000000000000000ull
);
1508 #ifdef SET_STATUS_FLAGS
1509 __set_status_flags (pfpsf
, ZERO_DIVIDE_EXCEPTION
);
1511 res
.w
[1] = (sign_x
^ sign_y
) | INFINITY_MASK64
;
1515 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1516 (void) fegetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1518 diff_expon
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS
;
1520 if (__unsigned_compare_gt_128 (CY
, CX
)) {
1527 fx
.d
= (float) CX
.w
[1] * f64
.d
+ (float) CX
.w
[0];
1528 fy
.d
= (float) CY
.w
[1] * f64
.d
+ (float) CY
.w
[0];
1529 // expon_cy - expon_cx
1530 bin_index
= (fy
.i
- fx
.i
) >> 23;
1533 T
= power10_index_binexp_128
[bin_index
].w
[0];
1534 __mul_64x128_short (CA
, T
, CX
);
1536 T128
= power10_index_binexp_128
[bin_index
];
1537 __mul_64x128_short (CA
, CX
.w
[0], T128
);
1541 if (__unsigned_compare_gt_128 (CY
, CA
))
1544 T128
= power10_table_128
[ed2
];
1545 __mul_128x128_to_256 (CA4
, CA
, T128
);
1547 ed2
+= estimate_decimal_digits
[bin_index
];
1548 CQ
.w
[0] = CQ
.w
[1] = 0;
1549 diff_expon
= diff_expon
- ed2
;
1553 __div_128_by_128 (&CQ
, &CR
, CX
, CY
);
1555 if (!CR
.w
[1] && !CR
.w
[0]) {
1556 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
,
1558 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1559 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1563 // get number of decimal digits in CQ
1566 fx
.d
= (float) CQ
.w
[1] * f64
.d
+ (float) CQ
.w
[0];
1567 // binary expon. of CQ
1568 bin_expon
= (fx
.i
- 0x3f800000) >> 23;
1570 digits_q
= estimate_decimal_digits
[bin_expon
];
1571 TP128
.w
[0] = power10_index_binexp_128
[bin_expon
].w
[0];
1572 TP128
.w
[1] = power10_index_binexp_128
[bin_expon
].w
[1];
1573 if (__unsigned_compare_ge_128 (CQ
, TP128
))
1576 ed2
= 34 - digits_q
;
1577 T128
.w
[0] = power10_table_128
[ed2
].w
[0];
1578 T128
.w
[1] = power10_table_128
[ed2
].w
[1];
1579 __mul_128x128_to_256 (CA4
, CR
, T128
);
1580 diff_expon
= diff_expon
- ed2
;
1581 __mul_128x128_low (CQ
, CQ
, T128
);
1585 __div_256_by_128 (&CQ
, &CA4
, CY
);
1588 #ifdef SET_STATUS_FLAGS
1589 if (CA4
.w
[0] || CA4
.w
[1]) {
1591 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
1593 #ifndef LEAVE_TRAILING_ZEROS
1597 #ifndef LEAVE_TRAILING_ZEROS
1598 if (!CA4
.w
[0] && !CA4
.w
[1])
1601 #ifndef LEAVE_TRAILING_ZEROS
1602 // check whether result is exact
1604 // check whether CX, CY are short
1605 if (!CX
.w
[1] && !CY
.w
[1] && (CX
.w
[0] <= 1024) && (CY
.w
[0] <= 1024)) {
1606 i
= (int) CY
.w
[0] - 1;
1607 j
= (int) CX
.w
[0] - 1;
1608 // difference in powers of 2 factors for Y and X
1609 nzeros
= ed2
- factors
[i
][0] + factors
[j
][0];
1610 // difference in powers of 5 factors
1611 d5
= ed2
- factors
[i
][1] + factors
[j
][1];
1614 // get P*(2^M[extra_digits])/10^extra_digits
1615 __mul_128x128_full (Qh
, Ql
, CQ
, reciprocals10_128
[nzeros
]);
1616 //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2];
1618 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1619 amount
= recip_scale
[nzeros
];
1620 __shr_128_long (CQ
, Qh
, amount
);
1622 diff_expon
+= nzeros
;
1624 // decompose Q as Qh*10^17 + Ql
1625 //T128 = reciprocals10_128[17];
1626 T128
.w
[0] = 0x44909befeb9fad49ull
;
1627 T128
.w
[1] = 0x000b877aa3236a4bull
;
1628 __mul_128x128_to_256 (P256
, CQ
, T128
);
1629 //amount = recip_scale[17];
1630 Q_high
= (P256
.w
[2] >> 44) | (P256
.w
[3] << (64 - 44));
1631 Q_low
= CQ
.w
[0] - Q_high
* 100000000000000000ull;
1636 tdigit
[0] = Q_high
& 0x3ffffff;
1642 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
1644 tdigit
[0] += convert_table
[j
][k
][0];
1645 tdigit
[1] += convert_table
[j
][k
][1];
1646 if (tdigit
[0] >= 100000000) {
1647 tdigit
[0] -= 100000000;
1653 if (tdigit
[1] >= 100000000) {
1654 tdigit
[1] -= 100000000;
1655 if (tdigit
[1] >= 100000000)
1656 tdigit
[1] -= 100000000;
1660 if (!digit
&& !tdigit
[1])
1668 PD
= (UINT64
) digit
*0x068DB8BBull
;
1669 digit_h
= (UINT32
) (PD
>> 40);
1670 digit_low
= digit
- digit_h
* 10000;
1675 digit_h
= digit_low
;
1679 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
1684 __mul_64x64_to_128 (CQ
, Q_high
, reciprocals10_64
[nzeros
]);
1686 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
1687 amount
= short_recip_scale
[nzeros
];
1688 CQ
.w
[0] = CQ
.w
[1] >> amount
;
1693 diff_expon
+= nzeros
;
1695 tdigit
[0] = Q_low
& 0x3ffffff;
1701 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
1703 tdigit
[0] += convert_table
[j
][k
][0];
1704 tdigit
[1] += convert_table
[j
][k
][1];
1705 if (tdigit
[0] >= 100000000) {
1706 tdigit
[0] -= 100000000;
1711 if (tdigit
[1] >= 100000000) {
1712 tdigit
[1] -= 100000000;
1713 if (tdigit
[1] >= 100000000)
1714 tdigit
[1] -= 100000000;
1718 if (!digit
&& !tdigit
[1])
1726 PD
= (UINT64
) digit
*0x068DB8BBull
;
1727 digit_h
= (UINT32
) (PD
>> 40);
1728 digit_low
= digit
- digit_h
* 10000;
1733 digit_h
= digit_low
;
1737 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
1742 // get P*(2^M[extra_digits])/10^extra_digits
1743 __mul_128x128_full (Qh
, Ql
, CQ
, reciprocals10_128
[nzeros
]);
1745 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
1746 amount
= recip_scale
[nzeros
];
1747 __shr_128 (CQ
, Qh
, amount
);
1749 diff_expon
+= nzeros
;
1753 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
,pfpsf
);
1754 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1755 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1761 if (diff_expon
>= 0) {
1762 #ifdef IEEE_ROUND_NEAREST
1765 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
1766 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
1767 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
1768 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
1770 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
1771 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
1774 if (CQ
.w
[0] < carry64
)
1777 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
1780 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
1781 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
1782 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
1783 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
1785 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
1786 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
1789 if (CQ
.w
[0] < carry64
)
1793 if (sign_x
^ sign_y
&& (unsigned) (rmode
- 1) < 2)
1796 case ROUNDING_TO_NEAREST
: // round to nearest code
1799 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
1800 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
1801 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
1802 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
1803 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
1804 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
1806 if (CQ
.w
[0] < carry64
)
1809 case ROUNDING_TIES_AWAY
:
1812 CA4r
.w
[1] = (CA4
.w
[1] + CA4
.w
[1]) | (CA4
.w
[0] >> 63);
1813 CA4r
.w
[0] = CA4
.w
[0] + CA4
.w
[0];
1814 __sub_borrow_out (CA4r
.w
[0], carry64
, CA4r
.w
[0], CY
.w
[0]);
1815 CA4r
.w
[1] = CA4r
.w
[1] - CY
.w
[1] - carry64
;
1816 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
1817 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
1819 if (CQ
.w
[0] < carry64
)
1823 case ROUNDING_TO_ZERO
:
1825 default: // rounding up
1835 #ifdef SET_STATUS_FLAGS
1836 if (CA4
.w
[0] || CA4
.w
[1]) {
1838 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
1841 handle_UF_128_rem (&res
, sign_x
^ sign_y
, diff_expon
, CQ
,
1842 CA4
.w
[1] | CA4
.w
[0], &rnd_mode
, pfpsf
);
1843 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1844 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1850 get_BID128 (&res
, sign_x
^ sign_y
, diff_expon
, CQ
, &rnd_mode
, pfpsf
);
1851 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1852 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);