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/>. */
24 /*****************************************************************************
26 *****************************************************************************
28 * Algorithm description:
30 * if(exponent_x < exponent_y)
31 * scale coefficient_y so exponents are aligned
32 * perform coefficient divide (64-bit integer divide), unless
33 * coefficient_y is longer than 64 bits (clearly larger
35 * else // exponent_x > exponent_y
36 * use a loop to scale coefficient_x to 18_digits, divide by
37 * coefficient_y (64-bit integer divide), calculate remainder
38 * as new_coefficient_x and repeat until final remainder is obtained
39 * (when new_exponent_x < exponent_y)
41 ****************************************************************************/
43 #include "bid_internal.h"
45 #define MAX_FORMAT_DIGITS 16
46 #define DECIMAL_EXPONENT_BIAS 398
47 #define MASK_BINARY_EXPONENT 0x7ff0000000000000ull
48 #define BINARY_EXPONENT_BIAS 0x3ff
49 #define UPPER_EXPON_LIMIT 51
51 #if DECIMAL_CALL_BY_REFERENCE
54 bid64_rem (UINT64
* pres
, UINT64
* px
,
56 py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM
) {
62 UINT64 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM
) {
65 UINT64 sign_x
, sign_y
, coefficient_x
, coefficient_y
, res
;
66 UINT64 Q
, R
, R2
, T
, valid_y
, valid_x
;
68 int exponent_x
, exponent_y
, bin_expon
, e_scale
;
69 int digits_x
, diff_expon
;
71 #if DECIMAL_CALL_BY_REFERENCE
76 valid_y
= unpack_BID64 (&sign_y
, &exponent_y
, &coefficient_y
, y
);
77 valid_x
= unpack_BID64 (&sign_x
, &exponent_x
, &coefficient_x
, x
);
79 // unpack arguments, check for NaN or Infinity
81 // x is Inf. or NaN or 0
82 #ifdef SET_STATUS_FLAGS
83 if ((y
& SNAN_MASK64
) == SNAN_MASK64
) // y is sNaN
84 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
88 if ((x
& 0x7c00000000000000ull
) == 0x7c00000000000000ull
) {
89 #ifdef SET_STATUS_FLAGS
90 if (((x
& SNAN_MASK64
) == SNAN_MASK64
))
91 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
93 res
= coefficient_x
& QUIET_MASK64
;;
97 if ((x
& 0x7800000000000000ull
) == 0x7800000000000000ull
) {
98 if (((y
& NAN_MASK64
) != NAN_MASK64
)) {
99 #ifdef SET_STATUS_FLAGS
100 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
103 res
= 0x7c00000000000000ull
;
108 // return x if y != 0
109 if (((y
& 0x7800000000000000ull
) < 0x7800000000000000ull
) &&
111 if ((y
& 0x6000000000000000ull
) == 0x6000000000000000ull
)
112 exponent_y
= (y
>> 51) & 0x3ff;
114 exponent_y
= (y
>> 53) & 0x3ff;
116 if (exponent_y
< exponent_x
)
117 exponent_x
= exponent_y
;
131 if ((y
& 0x7c00000000000000ull
) == 0x7c00000000000000ull
) {
132 #ifdef SET_STATUS_FLAGS
133 if (((y
& SNAN_MASK64
) == SNAN_MASK64
))
134 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
136 res
= coefficient_y
& QUIET_MASK64
;;
140 if ((y
& 0x7800000000000000ull
) == 0x7800000000000000ull
) {
141 res
= very_fast_get_BID64 (sign_x
, exponent_x
, coefficient_x
);
144 // y is 0, return NaN
146 #ifdef SET_STATUS_FLAGS
147 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
149 res
= 0x7c00000000000000ull
;
155 diff_expon
= exponent_x
- exponent_y
;
156 if (diff_expon
<= 0) {
157 diff_expon
= -diff_expon
;
159 if (diff_expon
> 16) {
160 // |x|<|y| in this case
164 // set exponent of y to exponent_x, scale coefficient_y
165 T
= power10_table_128
[diff_expon
].w
[0];
166 __mul_64x64_to_128 (CY
, coefficient_y
, T
);
168 if (CY
.w
[1] || CY
.w
[0] > (coefficient_x
<< 1)) {
173 Q
= coefficient_x
/ CY
.w
[0];
174 R
= coefficient_x
- Q
* CY
.w
[0];
177 if (R2
> CY
.w
[0] || (R2
== CY
.w
[0] && (Q
& 1))) {
179 sign_x
^= 0x8000000000000000ull
;
182 res
= very_fast_get_BID64 (sign_x
, exponent_x
, R
);
187 while (diff_expon
> 0) {
188 // get number of digits in coeff_x
189 tempx
.d
= (float) coefficient_x
;
190 bin_expon
= ((tempx
.i
>> 23) & 0xff) - 0x7f;
191 digits_x
= estimate_decimal_digits
[bin_expon
];
192 // will not use this test, dividend will have 18 or 19 digits
193 //if(coefficient_x >= power10_table_128[digits_x].w[0])
196 e_scale
= 18 - digits_x
;
197 if (diff_expon
>= e_scale
) {
198 diff_expon
-= e_scale
;
200 e_scale
= diff_expon
;
204 // scale dividend to 18 or 19 digits
205 coefficient_x
*= power10_table_128
[e_scale
].w
[0];
208 Q
= coefficient_x
/ coefficient_y
;
210 coefficient_x
-= Q
* coefficient_y
;
212 // check for remainder == 0
213 if (!coefficient_x
) {
214 res
= very_fast_get_BID64_small_mantissa (sign_x
, exponent_y
, 0);
219 R2
= coefficient_x
+ coefficient_x
;
220 if (R2
> coefficient_y
|| (R2
== coefficient_y
&& (Q
& 1))) {
221 coefficient_x
= coefficient_y
- coefficient_x
;
222 sign_x
^= 0x8000000000000000ull
;
225 res
= very_fast_get_BID64 (sign_x
, exponent_y
, coefficient_x
);