1 /* Copyright (C) 2007-2023 Free Software Foundation, Inc.
2
3 This file is part of GCC.
4
5 GCC is free software; you can redistribute it and/or modify it under
6 the terms of the GNU General Public License as published by the Free
7 Software Foundation; either version 3, or (at your option) any later
8 version.
9
10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11 WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13 for more details.
14
15 Under Section 7 of GPL version 3, you are granted additional
16 permissions described in the GCC Runtime Library Exception, version
17 3.1, as published by the Free Software Foundation.
18
19 You should have received a copy of the GNU General Public License and
20 a copy of the GCC Runtime Library Exception along with this program;
21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
22 <http://www.gnu.org/licenses/>. */
23
24 /*****************************************************************************
25 * BID64 remainder
26 *****************************************************************************
27 *
28 * Algorithm description:
29 *
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
34 * than coefficient_x)
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)
40 *
41 ****************************************************************************/
42
43 #include "bid_internal.h"
44
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
50
51 #if DECIMAL_CALL_BY_REFERENCE
52
53 void
54 bid64_rem (UINT64 * pres, UINT64 * px,
55 UINT64 *
56 py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
57 UINT64 x, y;
58 #else
59
60 UINT64
61 bid64_rem (UINT64 x,
62 UINT64 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
63 #endif
64 UINT128 CY;
65 UINT64 sign_x, sign_y, coefficient_x, coefficient_y, res;
66 UINT64 Q, R, R2, T, valid_y, valid_x;
67 int_float tempx;
68 int exponent_x, exponent_y, bin_expon, e_scale;
69 int digits_x, diff_expon;
70
71 #if DECIMAL_CALL_BY_REFERENCE
72 x = *px;
73 y = *py;
74 #endif
75
76 valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
77 valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
78
79 // unpack arguments, check for NaN or Infinity
80 if (!valid_x) {
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);
85 #endif
86
87 // test if x is NaN
88 if ((x & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
89 #ifdef SET_STATUS_FLAGS
90 if (((x & SNAN_MASK64) == SNAN_MASK64))
91 __set_status_flags (pfpsf, INVALID_EXCEPTION);
92 #endif
93 res = coefficient_x & QUIET_MASK64;;
94 BID_RETURN (res);
95 }
96 // x is Infinity?
97 if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {
98 if (((y & NAN_MASK64) != NAN_MASK64)) {
99 #ifdef SET_STATUS_FLAGS
100 __set_status_flags (pfpsf, INVALID_EXCEPTION);
101 #endif
102 // return NaN
103 res = 0x7c00000000000000ull;
104 BID_RETURN (res);
105 }
106 }
107 // x is 0
108 // return x if y != 0
109 if (((y & 0x7800000000000000ull) < 0x7800000000000000ull) &&
110 coefficient_y) {
111 if ((y & 0x6000000000000000ull) == 0x6000000000000000ull)
112 exponent_y = (y >> 51) & 0x3ff;
113 else
114 exponent_y = (y >> 53) & 0x3ff;
115
116 if (exponent_y < exponent_x)
117 exponent_x = exponent_y;
118
119 x = exponent_x;
120 x <<= 53;
121
122 res = x | sign_x;
123 BID_RETURN (res);
124 }
125
126 }
127 if (!valid_y) {
128 // y is Inf. or NaN
129
130 // test if y is NaN
131 if ((y & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
132 #ifdef SET_STATUS_FLAGS
133 if (((y & SNAN_MASK64) == SNAN_MASK64))
134 __set_status_flags (pfpsf, INVALID_EXCEPTION);
135 #endif
136 res = coefficient_y & QUIET_MASK64;;
137 BID_RETURN (res);
138 }
139 // y is Infinity?
140 if ((y & 0x7800000000000000ull) == 0x7800000000000000ull) {
141 res = very_fast_get_BID64 (sign_x, exponent_x, coefficient_x);
142 BID_RETURN (res);
143 }
144 // y is 0, return NaN
145 {
146 #ifdef SET_STATUS_FLAGS
147 __set_status_flags (pfpsf, INVALID_EXCEPTION);
148 #endif
149 res = 0x7c00000000000000ull;
150 BID_RETURN (res);
151 }
152 }
153
154
155 diff_expon = exponent_x - exponent_y;
156 if (diff_expon <= 0) {
157 diff_expon = -diff_expon;
158
159 if (diff_expon > 16) {
160 // |x|<|y| in this case
161 res = x;
162 BID_RETURN (res);
163 }
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);
167
168 if (CY.w[1] || CY.w[0] > (coefficient_x << 1)) {
169 res = x;
170 BID_RETURN (res);
171 }
172
173 Q = coefficient_x / CY.w[0];
174 R = coefficient_x - Q * CY.w[0];
175
176 R2 = R + R;
177 if (R2 > CY.w[0] || (R2 == CY.w[0] && (Q & 1))) {
178 R = CY.w[0] - R;
179 sign_x ^= 0x8000000000000000ull;
180 }
181
182 res = very_fast_get_BID64 (sign_x, exponent_x, R);
183 BID_RETURN (res);
184 }
185
186
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])
194 // digits_x++;
195
196 e_scale = 18 - digits_x;
197 if (diff_expon >= e_scale) {
198 diff_expon -= e_scale;
199 } else {
200 e_scale = diff_expon;
201 diff_expon = 0;
202 }
203
204 // scale dividend to 18 or 19 digits
205 coefficient_x *= power10_table_128[e_scale].w[0];
206
207 // quotient
208 Q = coefficient_x / coefficient_y;
209 // remainder
210 coefficient_x -= Q * coefficient_y;
211
212 // check for remainder == 0
213 if (!coefficient_x) {
214 res = very_fast_get_BID64_small_mantissa (sign_x, exponent_y, 0);
215 BID_RETURN (res);
216 }
217 }
218
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;
223 }
224
225 res = very_fast_get_BID64 (sign_x, exponent_y, coefficient_x);
226 BID_RETURN (res);
227
228 }