1 /* Copyright (C) 2007-2018 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/>. */
24 /*****************************************************************************
26 *****************************************************************************
28 * Algorithm description:
30 * if(coefficient_x<coefficient_y)
31 * p = number_digits(coefficient_y) - number_digits(coefficient_x)
32 * A = coefficient_x*10^p
34 * CA= A*10^(15+j), j=0 for A>=B, 1 otherwise
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)
46 * Q += CA/B (64-bit unsigned divide)
48 * get final Q using double precision divide, followed by 3 integer
50 * if exact result, eliminate trailing zeros
52 * round coefficient to nearest
54 ****************************************************************************/
56 #include "bid_internal.h"
57 #include "bid_div_macros.h"
58 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
61 #define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT
64 extern UINT32 convert_table
[5][128][2];
65 extern SINT8 factors
[][2];
66 extern UINT8 packed_10000_zeros
[];
69 #if DECIMAL_CALL_BY_REFERENCE
72 bid64_div (UINT64
* pres
, UINT64
* px
,
74 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
81 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
82 _EXC_MASKS_PARAM _EXC_INFO_PARAM
) {
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
;
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
;
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;
101 #if DECIMAL_CALL_BY_REFERENCE
102 #if !DECIMAL_GLOBAL_ROUNDING
103 _IDEC_round rnd_mode
= *prnd_mode
;
109 valid_x
= unpack_BID64 (&sign_x
, &exponent_x
, &coefficient_x
, x
);
110 valid_y
= unpack_BID64 (&sign_y
, &exponent_y
, &coefficient_y
, y
);
112 // unpack arguments, check for NaN or Infinity
115 #ifdef SET_STATUS_FLAGS
116 if ((y
& SNAN_MASK64
) == SNAN_MASK64
) // y is sNaN
117 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
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
);
126 BID_RETURN (coefficient_x
& QUIET_MASK64
);
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
);
137 BID_RETURN (NAN_MASK64
);
140 // otherwise return +/-Inf
141 BID_RETURN (((x
^ y
) & 0x8000000000000000ull
) |
146 if (((y
& INFINITY_MASK64
) != INFINITY_MASK64
)
147 && !(coefficient_y
)) {
149 #ifdef SET_STATUS_FLAGS
150 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
152 BID_RETURN (NAN_MASK64
);
154 if (((y
& INFINITY_MASK64
) != INFINITY_MASK64
)) {
155 if ((y
& SPECIAL_ENCODING_MASK64
) == SPECIAL_ENCODING_MASK64
)
156 exponent_y
= ((UINT32
) (y
>> 51)) & 0x3ff;
158 exponent_y
= ((UINT32
) (y
>> 53)) & 0x3ff;
159 sign_y
= y
& 0x8000000000000000ull
;
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)
166 BID_RETURN ((sign_x
^ sign_y
) | (((UINT64
) exponent_x
) << 53));
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
);
179 BID_RETURN (coefficient_y
& QUIET_MASK64
);
182 if ((y
& INFINITY_MASK64
) == INFINITY_MASK64
) {
184 BID_RETURN (((x
^ y
) & 0x8000000000000000ull
));
187 #ifdef SET_STATUS_FLAGS
188 __set_status_flags (pfpsf
, ZERO_DIVIDE_EXCEPTION
);
190 BID_RETURN ((sign_x
^ sign_y
) | INFINITY_MASK64
);
192 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
193 (void) fegetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
195 diff_expon
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS
;
197 if (coefficient_x
< coefficient_y
) {
198 // get number of decimal digits for c_x, c_y
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;
205 A
= coefficient_x
* power10_index_binexp
[bin_index
];
208 temp_b
.d
= (double) B
;
213 ed2
= estimate_decimal_digits
[bin_index
] + ed1
;
214 T
= power10_table_128
[ed1
].w
[0];
215 __mul_64x64_to_128 (CA
, A
, T
);
218 diff_expon
= diff_expon
- ed2
;
220 // adjust double precision db, to ensure that later A/B - (int)(da/db) > -1
221 if (coefficient_y
< 0x0020000000000000ull
) {
225 db
= (double) (B
+ 2 + (B
& 1));
230 // set last bit before conversion to DP
231 A2
= coefficient_x
| 1;
234 db
= (double) coefficient_y
;
237 Q
= (UINT64
) tempq
.d
;
239 R
= coefficient_x
- coefficient_y
* Q
;
241 // will use to get number of dec. digits of Q
242 bin_expon_cx
= (tempq
.i
>> 52) - 0x3ff;
245 D
= ((SINT64
) R
) >> 63;
247 R
+= (coefficient_y
& D
);
250 if (((SINT64
) R
) <= 0) {
251 // can have R==-1 for coeff_y==1
253 get_BID64 (sign_x
^ sign_y
, diff_expon
, (Q
+ R
), rnd_mode
,
255 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
256 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
260 // get decimal digits of Q
261 DU
= power10_index_binexp
[bin_expon_cx
] - Q
- 1;
264 ed2
= 16 - estimate_decimal_digits
[bin_expon_cx
] - (int) DU
;
266 T
= power10_table_128
[ed2
].w
[0];
267 __mul_64x64_to_128 (CA
, R
, T
);
270 Q
*= power10_table_128
[ed2
].w
[0];
279 R
= CA
.w
[0] - Q2
* B
;
284 t_scale
.i
= 0x43f0000000000000ull
;
288 da
= da_h
* t_scale
.d
+ da_l
;
294 // get w[0] remainder
295 R
= CA
.w
[0] - Q2
* B
;
298 D
= ((SINT64
) R
) >> 63;
312 D
= ((SINT64
) R
) >> 63;
313 // restore R if negative
319 D
= ((SINT64
) R
) >> 63;
320 // restore R if negative
326 D
= ((SINT64
) R
) >> 63;
327 // restore R if negative
334 #ifdef SET_STATUS_FLAGS
337 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
339 #ifndef LEAVE_TRAILING_ZEROS
343 #ifndef LEAVE_TRAILING_ZEROS
347 #ifndef LEAVE_TRAILING_ZEROS
349 // eliminate trailing zeros
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];
362 __mul_64x64_to_128 (CT
, Q
, reciprocals10_64
[nzeros
]);
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
;
368 diff_expon
+= nzeros
;
370 tdigit
[0] = Q
& 0x3ffffff;
376 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
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;
387 if (!digit
&& !tdigit
[1])
395 PD
= (UINT64
) digit
*0x068DB8BBull
;
396 digit_h
= (UINT32
) (PD
>> 40);
397 digit_low
= digit
- digit_h
* 10000;
406 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
411 __mul_64x64_to_128 (CT
, Q
, reciprocals10_64
[nzeros
]);
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
;
417 diff_expon
+= nzeros
;
420 if (diff_expon
>= 0) {
422 fast_get_BID64_check_OF (sign_x
^ sign_y
, diff_expon
, Q
,
424 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
425 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
432 if (diff_expon
>= 0) {
433 #ifdef IEEE_ROUND_NEAREST
434 // round to nearest code
440 // compare 10*R to 5*B
442 // correction for (R==0 && (Q&1))
445 D
= ((UINT64
) R
) >> 63;
448 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
449 // round to nearest code
455 // compare 10*R to 5*B
457 // correction for (R==0 && (Q&1))
460 D
= ((UINT64
) R
) >> 63;
464 if (sign_x
^ sign_y
&& (unsigned) (rmode
- 1) < 2)
467 case 0: // round to nearest code
468 case ROUNDING_TIES_AWAY
:
473 // compare 10*R to 5*B
475 // correction for (R==0 && (Q&1))
476 R
-= ((Q
| (rmode
>> 2)) & 1);
478 D
= ((UINT64
) R
) >> 63;
482 case ROUNDING_TO_ZERO
:
484 default: // rounding up
492 fast_get_BID64_check_OF (sign_x
^ sign_y
, diff_expon
, Q
, rnd_mode
,
494 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
495 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
501 #ifdef SET_STATUS_FLAGS
502 if ((diff_expon
+ 16 < 0)) {
504 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
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
);
520 TYPE0_FUNCTION_ARGTYPE1_ARG128 (UINT64
, bid64dq_div
, UINT64
, x
, y
)
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
,
529 int nzeros
, i
, j
, k
, d5
, done
= 0;
531 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
532 fexcept_t binaryflags
= 0;
535 valid_y
= unpack_BID128_value (&sign_y
, &exponent_y
, &CY
, y
);
537 // unpack arguments, check for NaN or Infinity
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
);
546 if (((x
) & 0x7c00000000000000ull
) == 0x7c00000000000000ull
) {
548 BID_RETURN (res
& QUIET_MASK64
);
551 if (((x
) & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
552 // check if y is Inf.
553 if (((y
.w
[1] & 0x7c00000000000000ull
) == 0x7800000000000000ull
))
556 #ifdef SET_STATUS_FLAGS
557 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
559 res
= 0x7c00000000000000ull
;
562 if (((y
.w
[1] & 0x7c00000000000000ull
) != 0x7c00000000000000ull
)) {
563 // otherwise return +/-Inf
565 (((x
) ^ y
.w
[1]) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
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
);
576 res
= 0x7c00000000000000ull
;
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)
586 res
|= (((UINT64
) exponent_x
) << 53);
590 exponent_x
+= (DECIMAL_EXPONENT_BIAS_128
- DECIMAL_EXPONENT_BIAS
);
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
);
600 Tmp
.w
[1] = (CY
.w
[1] & 0x00003fffffffffffull
);
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];
610 if ((y
.w
[1] & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
612 res
= sign_x
^ sign_y
;
615 // y is 0, return +/-Inf
617 (((x
) ^ y
.w
[1]) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
618 #ifdef SET_STATUS_FLAGS
619 __set_status_flags (pfpsf
, ZERO_DIVIDE_EXCEPTION
);
623 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
624 (void) fegetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
626 diff_expon
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS
;
628 if (__unsigned_compare_gt_128 (CY
, CX
)) {
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;
641 T
= power10_index_binexp_128
[bin_index
].w
[0];
642 __mul_64x128_short (CA
, T
, CX
);
644 T128
= power10_index_binexp_128
[bin_index
];
645 __mul_64x128_short (CA
, CX
.w
[0], T128
);
649 if (__unsigned_compare_gt_128 (CY
, CA
))
652 T128
= power10_table_128
[ed2
];
653 __mul_128x128_to_256 (CA4
, CA
, T128
);
655 ed2
+= estimate_decimal_digits
[bin_index
];
656 CQ
.w
[0] = CQ
.w
[1] = 0;
657 diff_expon
= diff_expon
- ed2
;
661 __div_128_by_128 (&CQ
, &CR
, CX
, CY
);
663 // get number of decimal digits in CQ
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;
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
))
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
);
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];
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
));
700 __mul_64x64_to_128 (CQ2
, CQ
.w
[0], (power10_table_128
[ed2
].w
[0]));
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];
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])
709 if (CR
.w
[0] || CR
.w
[1])
717 __div_256_by_128 (&CQ
, &CA4
, CY
);
722 #ifdef SET_STATUS_FLAGS
723 if (CA4
.w
[0] || CA4
.w
[1]) {
725 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
727 #ifndef LEAVE_TRAILING_ZEROS
731 #ifndef LEAVE_TRAILING_ZEROS
732 if (!CA4
.w
[0] && !CA4
.w
[1])
735 #ifndef LEAVE_TRAILING_ZEROS
736 // check whether result is exact
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];
748 // get P*(2^M[extra_digits])/10^extra_digits
749 __mul_128x128_high (Qh
, CQ
, reciprocals10_128
[nzeros
]);
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
);
755 diff_expon
+= nzeros
;
757 // decompose Q as Qh*10^17 + Ql
761 tdigit
[0] = Q_low
& 0x3ffffff;
767 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
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;
777 if (tdigit
[1] >= 100000000) {
778 tdigit
[1] -= 100000000;
779 if (tdigit
[1] >= 100000000)
780 tdigit
[1] -= 100000000;
784 if (!digit
&& !tdigit
[1])
792 PD
= (UINT64
) digit
*0x068DB8BBull
;
793 digit_h
= (UINT32
) (PD
>> 40);
794 digit_low
= digit
- digit_h
* 10000;
803 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
808 // get P*(2^M[extra_digits])/10^extra_digits
809 __mul_128x128_high (Qh
, CQ
, reciprocals10_128
[nzeros
]);
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
);
815 diff_expon
+= nzeros
;
821 fast_get_BID64_check_OF (sign_x
^ sign_y
, diff_expon
, CQ
.w
[0],
823 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
824 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
831 if (diff_expon
>= 0) {
832 #ifdef IEEE_ROUND_NEAREST
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
;
840 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
841 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
845 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
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
;
853 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
854 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
857 if (CQ
.w
[0] < carry64
)
861 if (sign_x
^ sign_y
&& (unsigned) (rmode
- 1) < 2)
864 case ROUNDING_TO_NEAREST
: // round to nearest code
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
);
874 if (CQ
.w
[0] < carry64
)
877 case ROUNDING_TIES_AWAY
:
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
;
887 if (CQ
.w
[0] < carry64
)
891 case ROUNDING_TO_ZERO
:
893 default: // rounding up
903 fast_get_BID64_check_OF (sign_x
^ sign_y
, diff_expon
, CQ
.w
[0], rnd_mode
,
905 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
906 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
912 #ifdef SET_STATUS_FLAGS
913 if ((diff_expon
+ 16 < 0)) {
915 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
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
);
931 //#define LEAVE_TRAILING_ZEROS
933 TYPE0_FUNCTION_ARG128_ARGTYPE2 (UINT64
, bid64qd_div
, x
, UINT64
, y
)
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
,
943 int nzeros
, i
, j
, k
, d5
, done
= 0;
945 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
946 fexcept_t binaryflags
= 0;
949 valid_y
= unpack_BID64 (&sign_y
, &exponent_y
, &CY
.w
[0], (y
));
951 // unpack arguments, check for NaN or Infinity
952 if (!unpack_BID128_value (&sign_x
, &exponent_x
, &CX
, x
)) {
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
);
960 Tmp
.w
[1] = (CX
.w
[1] & 0x00003fffffffffffull
);
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];
970 if ((x
.w
[1] & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
971 // check if y is Inf.
972 if (((y
& 0x7c00000000000000ull
) == 0x7800000000000000ull
))
975 #ifdef SET_STATUS_FLAGS
976 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
978 res
= 0x7c00000000000000ull
;
981 if (((y
& 0x7c00000000000000ull
) != 0x7c00000000000000ull
)) {
982 // otherwise return +/-Inf
984 ((x
.w
[1] ^ (y
)) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
989 if (((y
& INFINITY_MASK64
) != INFINITY_MASK64
) &&
991 #ifdef SET_STATUS_FLAGS
992 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
995 res
= 0x7c00000000000000ull
;
999 if (((y
& 0x7800000000000000ull
) != 0x7800000000000000ull
)) {
1001 #ifdef SET_STATUS_FLAGS
1002 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
1004 res
= 0x7c00000000000000ull
;
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)
1014 res
= (sign_x
^ sign_y
) | (((UINT64
) exponent_x
) << 53);
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
);
1028 BID_RETURN (CY
.w
[0] & QUIET_MASK64
);
1031 if (((y
) & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
1033 res
= sign_x
^ sign_y
;
1036 // y is 0, return +/-Inf
1038 ((x
.w
[1] ^ (y
)) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
1039 #ifdef SET_STATUS_FLAGS
1040 __set_status_flags (pfpsf
, ZERO_DIVIDE_EXCEPTION
);
1044 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1045 (void) fegetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1048 exponent_x
- exponent_y
- DECIMAL_EXPONENT_BIAS_128
+
1049 (DECIMAL_EXPONENT_BIAS
<< 1);
1051 if (__unsigned_compare_gt_128 (CY
, CX
)) {
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;
1064 T
= power10_index_binexp_128
[bin_index
].w
[0];
1065 __mul_64x128_short (CA
, T
, CX
);
1067 T128
= power10_index_binexp_128
[bin_index
];
1068 __mul_64x128_short (CA
, CX
.w
[0], T128
);
1072 if (__unsigned_compare_gt_128 (CY
, CA
))
1075 T128
= power10_table_128
[ed2
];
1076 __mul_128x128_to_256 (CA4
, CA
, T128
);
1078 ed2
+= estimate_decimal_digits
[bin_index
];
1079 CQ
.w
[0] = CQ
.w
[1] = 0;
1080 diff_expon
= diff_expon
- ed2
;
1084 __div_128_by_128 (&CQ
, &CR
, CX
, CY
);
1086 // get number of decimal digits in CQ
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;
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
))
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
);
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];
1115 ed2
= digits_q
- 16;
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
));
1123 __mul_64x64_to_128 (CQ2
, CQ
.w
[0], (power10_table_128
[ed2
].w
[0]));
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];
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])
1132 if (CR
.w
[0] || CR
.w
[1])
1135 if(CA4
.w
[1]|CA4
.w
[0]) {
1136 __mul_64x128_low(CY
, (power10_table_128
[ed2
].w
[0]),CY
);
1144 __div_256_by_128 (&CQ
, &CA4
, CY
);
1147 #ifdef SET_STATUS_FLAGS
1148 if (CA4
.w
[0] || CA4
.w
[1]) {
1150 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
1152 #ifndef LEAVE_TRAILING_ZEROS
1156 #ifndef LEAVE_TRAILING_ZEROS
1157 if (!CA4
.w
[0] && !CA4
.w
[1])
1160 #ifndef LEAVE_TRAILING_ZEROS
1161 // check whether result is exact
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];
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];
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
);
1182 diff_expon
+= nzeros
;
1184 // decompose Q as Qh*10^17 + Ql
1185 //T128 = reciprocals10_128[17];
1189 tdigit
[0] = Q_low
& 0x3ffffff;
1195 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
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;
1205 if (tdigit
[1] >= 100000000) {
1206 tdigit
[1] -= 100000000;
1207 if (tdigit
[1] >= 100000000)
1208 tdigit
[1] -= 100000000;
1212 if (!digit
&& !tdigit
[1])
1220 PD
= (UINT64
) digit
*0x068DB8BBull
;
1221 digit_h
= (UINT32
) (PD
>> 40);
1222 digit_low
= digit
- digit_h
* 10000;
1227 digit_h
= digit_low
;
1231 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
1236 // get P*(2^M[extra_digits])/10^extra_digits
1237 __mul_128x128_high (Qh
, CQ
, reciprocals10_128
[nzeros
]);
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
);
1243 diff_expon
+= nzeros
;
1250 fast_get_BID64_check_OF (sign_x
^ sign_y
, diff_expon
, CQ
.w
[0],
1252 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1253 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1260 if (diff_expon
>= 0) {
1261 #ifdef IEEE_ROUND_NEAREST
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
;
1269 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
1270 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
1273 //if(CQ.w[0]<carry64)
1276 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
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
;
1284 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
1285 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
1288 if (CQ
.w
[0] < carry64
)
1292 if (sign_x
^ sign_y
&& (unsigned) (rmode
- 1) < 2)
1295 case ROUNDING_TO_NEAREST
: // round to nearest code
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
);
1305 if (CQ
.w
[0] < carry64
)
1308 case ROUNDING_TIES_AWAY
:
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
;
1318 if (CQ
.w
[0] < carry64
)
1322 case ROUNDING_TO_ZERO
:
1324 default: // rounding up
1335 fast_get_BID64_check_OF (sign_x
^ sign_y
, diff_expon
, CQ
.w
[0], rnd_mode
,
1337 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1338 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1344 #ifdef SET_STATUS_FLAGS
1345 if ((diff_expon
+ 16 < 0)) {
1347 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
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
);
1362 //#define LEAVE_TRAILING_ZEROS
1364 extern UINT32 convert_table
[5][128][2];
1365 extern SINT8 factors
[][2];
1366 extern UINT8 packed_10000_zeros
[];
1369 //UINT64* bid64_div128x128(UINT64 res, UINT128 *px, UINT128 *py, unsigned rnd_mode, unsigned *pfpsf)
1371 TYPE0_FUNCTION_ARG128_ARG128 (UINT64
, bid64qq_div
, x
, y
)
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
,
1380 int nzeros
, i
, j
, k
, d5
, done
= 0;
1382 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1383 fexcept_t binaryflags
= 0;
1386 valid_y
= unpack_BID128_value (&sign_y
, &exponent_y
, &CY
, y
);
1388 // unpack arguments, check for NaN or Infinity
1389 if (!unpack_BID128_value (&sign_x
, &exponent_x
, &CX
, x
)) {
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
);
1397 Tmp
.w
[1] = (CX
.w
[1] & 0x00003fffffffffffull
);
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];
1407 if ((x
.w
[1] & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
1408 // check if y is Inf.
1409 if (((y
.w
[1] & 0x7c00000000000000ull
) == 0x7800000000000000ull
))
1412 #ifdef SET_STATUS_FLAGS
1413 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
1415 res
= 0x7c00000000000000ull
;
1418 if (((y
.w
[1] & 0x7c00000000000000ull
) != 0x7c00000000000000ull
)) {
1419 // otherwise return +/-Inf
1422 w
[1]) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
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
);
1432 // x=y=0, return NaN
1433 res
= 0x7c00000000000000ull
;
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)
1443 res
|= (((UINT64
) exponent_x
) << 53);
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
);
1456 Tmp
.w
[1] = (CY
.w
[1] & 0x00003fffffffffffull
);
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];
1466 if ((y
.w
[1] & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
1468 res
= sign_x
^ sign_y
;
1471 // y is 0, return +/-Inf
1473 ((x
.w
[1] ^ y
.w
[1]) & 0x8000000000000000ull
) | 0x7800000000000000ull
;
1474 #ifdef SET_STATUS_FLAGS
1475 __set_status_flags (pfpsf
, ZERO_DIVIDE_EXCEPTION
);
1479 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1480 (void) fegetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1482 diff_expon
= exponent_x
- exponent_y
+ DECIMAL_EXPONENT_BIAS
;
1484 if (__unsigned_compare_gt_128 (CY
, CX
)) {
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;
1497 T
= power10_index_binexp_128
[bin_index
].w
[0];
1498 __mul_64x128_short (CA
, T
, CX
);
1500 T128
= power10_index_binexp_128
[bin_index
];
1501 __mul_64x128_short (CA
, CX
.w
[0], T128
);
1505 if (__unsigned_compare_gt_128 (CY
, CA
))
1508 T128
= power10_table_128
[ed2
];
1509 __mul_128x128_to_256 (CA4
, CA
, T128
);
1511 ed2
+= estimate_decimal_digits
[bin_index
];
1512 CQ
.w
[0] = CQ
.w
[1] = 0;
1513 diff_expon
= diff_expon
- ed2
;
1517 __div_128_by_128 (&CQ
, &CR
, CX
, CY
);
1519 // get number of decimal digits in CQ
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;
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
))
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
);
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];
1548 ed2
= digits_q
- 16;
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
));
1556 __mul_64x64_to_128 (CQ2
, CQ
.w
[0], (power10_table_128
[ed2
].w
[0]));
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];
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])
1565 if (CR
.w
[0] || CR
.w
[1])
1568 if(CA4
.w
[1]|CA4
.w
[0]) {
1569 __mul_64x128_low(CY
, (power10_table_128
[ed2
].w
[0]),CY
);
1576 __div_256_by_128 (&CQ
, &CA4
, CY
);
1581 #ifdef SET_STATUS_FLAGS
1582 if (CA4
.w
[0] || CA4
.w
[1]) {
1584 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
1586 #ifndef LEAVE_TRAILING_ZEROS
1590 #ifndef LEAVE_TRAILING_ZEROS
1591 if (!CA4
.w
[0] && !CA4
.w
[1])
1594 #ifndef LEAVE_TRAILING_ZEROS
1595 // check whether result is exact
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];
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];
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
);
1616 diff_expon
+= nzeros
;
1618 // decompose Q as Qh*10^17 + Ql
1619 //T128 = reciprocals10_128[17];
1623 tdigit
[0] = Q_low
& 0x3ffffff;
1629 for (j
= 0; QX32
; j
++, QX32
>>= 7) {
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;
1639 if (tdigit
[1] >= 100000000) {
1640 tdigit
[1] -= 100000000;
1641 if (tdigit
[1] >= 100000000)
1642 tdigit
[1] -= 100000000;
1646 if (!digit
&& !tdigit
[1])
1654 PD
= (UINT64
) digit
*0x068DB8BBull
;
1655 digit_h
= (UINT32
) (PD
>> 40);
1656 digit_low
= digit
- digit_h
* 10000;
1661 digit_h
= digit_low
;
1665 3 & (UINT32
) (packed_10000_zeros
[digit_h
>> 3] >>
1670 // get P*(2^M[extra_digits])/10^extra_digits
1671 __mul_128x128_high (Qh
, CQ
, reciprocals10_128
[nzeros
]);
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
);
1677 diff_expon
+= nzeros
;
1684 fast_get_BID64_check_OF (sign_x
^ sign_y
, diff_expon
, CQ
.w
[0],
1686 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1687 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1696 #ifdef IEEE_ROUND_NEAREST
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
;
1704 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 1 : 0;
1705 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) & ((CQ
.w
[0]) | D
);
1708 //if(CQ.w[0]<carry64)
1711 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
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
;
1719 D
= (CA4r
.w
[1] | CA4r
.w
[0]) ? 0 : 1;
1720 carry64
= (1 + (((SINT64
) CA4r
.w
[1]) >> 63)) | D
;
1723 if (CQ
.w
[0] < carry64
)
1727 if (sign_x
^ sign_y
&& (unsigned) (rmode
- 1) < 2)
1730 case ROUNDING_TO_NEAREST
: // round to nearest code
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
);
1740 if (CQ
.w
[0] < carry64
)
1743 case ROUNDING_TIES_AWAY
:
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
;
1753 if (CQ
.w
[0] < carry64
)
1757 case ROUNDING_TO_ZERO
:
1759 default: // rounding up
1770 fast_get_BID64_check_OF (sign_x
^ sign_y
, diff_expon
, CQ
.w
[0], rnd_mode
,
1772 #ifdef UNCHANGED_BINARY_STATUS_FLAGS
1773 (void) fesetexceptflag (&binaryflags
, FE_ALL_FLAGS
);
1779 #ifdef SET_STATUS_FLAGS
1780 if ((diff_expon
+ 16 < 0)) {
1782 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
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
);