]>
Commit | Line | Data |
---|---|---|
748086b7 | 1 | /* Copyright (C) 2007, 2009 Free Software Foundation, Inc. |
200359e8 L |
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 | |
748086b7 | 7 | Software Foundation; either version 3, or (at your option) any later |
200359e8 L |
8 | version. |
9 | ||
200359e8 L |
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 | ||
748086b7 JJ |
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/>. */ | |
200359e8 L |
23 | |
24 | #include "bid_internal.h" | |
25 | ||
26 | /***************************************************************************** | |
27 | * BID64_round_integral_exact | |
28 | ****************************************************************************/ | |
29 | ||
30 | #if DECIMAL_CALL_BY_REFERENCE | |
31 | void | |
b2a00c89 | 32 | bid64_round_integral_exact (UINT64 * pres, |
200359e8 L |
33 | UINT64 * |
34 | px _RND_MODE_PARAM _EXC_FLAGS_PARAM | |
35 | _EXC_MASKS_PARAM _EXC_INFO_PARAM) { | |
36 | UINT64 x = *px; | |
37 | #if !DECIMAL_GLOBAL_ROUNDING | |
38 | unsigned int rnd_mode = *prnd_mode; | |
39 | #endif | |
40 | #else | |
41 | UINT64 | |
b2a00c89 | 42 | bid64_round_integral_exact (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM |
200359e8 L |
43 | _EXC_MASKS_PARAM _EXC_INFO_PARAM) { |
44 | #endif | |
45 | ||
46 | UINT64 res = 0xbaddbaddbaddbaddull; | |
47 | UINT64 x_sign; | |
b2a00c89 | 48 | int exp; // unbiased exponent |
200359e8 L |
49 | // Note: C1 represents the significand (UINT64) |
50 | BID_UI64DOUBLE tmp1; | |
51 | int x_nr_bits; | |
52 | int q, ind, shift; | |
53 | UINT64 C1; | |
54 | // UINT64 res is C* at first - represents up to 16 decimal digits <= 54 bits | |
b2a00c89 | 55 | UINT128 fstar = { {0x0ull, 0x0ull} }; |
200359e8 L |
56 | UINT128 P128; |
57 | ||
b2a00c89 L |
58 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
59 | ||
60 | // check for NaNs and infinities | |
61 | if ((x & MASK_NAN) == MASK_NAN) { // check for NaN | |
62 | if ((x & 0x0003ffffffffffffull) > 999999999999999ull) | |
63 | x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits | |
64 | else | |
65 | x = x & 0xfe03ffffffffffffull; // clear G6-G12 | |
66 | if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN | |
200359e8 L |
67 | // set invalid flag |
68 | *pfpsf |= INVALID_EXCEPTION; | |
b2a00c89 | 69 | // return quiet (SNaN) |
200359e8 | 70 | res = x & 0xfdffffffffffffffull; |
b2a00c89 L |
71 | } else { // QNaN |
72 | res = x; | |
200359e8 | 73 | } |
b2a00c89 L |
74 | BID_RETURN (res); |
75 | } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity | |
76 | res = x_sign | 0x7800000000000000ull; | |
200359e8 L |
77 | BID_RETURN (res); |
78 | } | |
79 | // unpack x | |
200359e8 L |
80 | if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { |
81 | // if the steering bits are 11 (condition will be 0), then | |
82 | // the exponent is G[0:w+1] | |
83 | exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398; | |
84 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; | |
b2a00c89 | 85 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
86 | C1 = 0; |
87 | } | |
b2a00c89 | 88 | } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) |
200359e8 L |
89 | exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398; |
90 | C1 = (x & MASK_BINARY_SIG1); | |
91 | } | |
92 | ||
93 | // if x is 0 or non-canonical return 0 preserving the sign bit and | |
94 | // the preferred exponent of MAX(Q(x), 0) | |
95 | if (C1 == 0) { | |
96 | if (exp < 0) | |
97 | exp = 0; | |
98 | res = x_sign | (((UINT64) exp + 398) << 53); | |
99 | BID_RETURN (res); | |
100 | } | |
101 | // x is a finite non-zero number (not 0, non-canonical, or special) | |
102 | ||
103 | switch (rnd_mode) { | |
104 | case ROUNDING_TO_NEAREST: | |
105 | case ROUNDING_TIES_AWAY: | |
106 | // return 0 if (exp <= -(p+1)) | |
107 | if (exp <= -17) { | |
108 | res = x_sign | 0x31c0000000000000ull; | |
109 | *pfpsf |= INEXACT_EXCEPTION; | |
110 | BID_RETURN (res); | |
111 | } | |
112 | break; | |
113 | case ROUNDING_DOWN: | |
114 | // return 0 if (exp <= -p) | |
115 | if (exp <= -16) { | |
116 | if (x_sign) { | |
117 | res = 0xb1c0000000000001ull; | |
118 | } else { | |
119 | res = 0x31c0000000000000ull; | |
120 | } | |
121 | *pfpsf |= INEXACT_EXCEPTION; | |
122 | BID_RETURN (res); | |
123 | } | |
124 | break; | |
125 | case ROUNDING_UP: | |
126 | // return 0 if (exp <= -p) | |
127 | if (exp <= -16) { | |
128 | if (x_sign) { | |
129 | res = 0xb1c0000000000000ull; | |
130 | } else { | |
131 | res = 0x31c0000000000001ull; | |
132 | } | |
133 | *pfpsf |= INEXACT_EXCEPTION; | |
134 | BID_RETURN (res); | |
135 | } | |
136 | break; | |
137 | case ROUNDING_TO_ZERO: | |
138 | // return 0 if (exp <= -p) | |
139 | if (exp <= -16) { | |
140 | res = x_sign | 0x31c0000000000000ull; | |
141 | *pfpsf |= INEXACT_EXCEPTION; | |
142 | BID_RETURN (res); | |
143 | } | |
144 | break; | |
b2a00c89 | 145 | } // end switch () |
200359e8 L |
146 | |
147 | // q = nr. of decimal digits in x (1 <= q <= 54) | |
148 | // determine first the nr. of bits in x | |
b2a00c89 | 149 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 150 | q = 16; |
b2a00c89 | 151 | } else { // if x < 2^53 |
200359e8 L |
152 | tmp1.d = (double) C1; // exact conversion |
153 | x_nr_bits = | |
154 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 | 155 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 156 | if (q == 0) { |
b2a00c89 L |
157 | q = nr_digits[x_nr_bits - 1].digits1; |
158 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
159 | q++; |
160 | } | |
161 | } | |
162 | ||
b2a00c89 | 163 | if (exp >= 0) { // -exp <= 0 |
200359e8 L |
164 | // the argument is an integer already |
165 | res = x; | |
166 | BID_RETURN (res); | |
167 | } | |
168 | ||
169 | switch (rnd_mode) { | |
170 | case ROUNDING_TO_NEAREST: | |
b2a00c89 | 171 | if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q |
200359e8 L |
172 | // need to shift right -exp digits from the coefficient; exp will be 0 |
173 | ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' | |
174 | // chop off ind digits from the lower part of C1 | |
175 | // C1 = C1 + 1/2 * 10^x where the result C1 fits in 64 bits | |
176 | // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate | |
b2a00c89 | 177 | C1 = C1 + midpoint64[ind - 1]; |
200359e8 L |
178 | // calculate C* and f* |
179 | // C* is actually floor(C*) in this case | |
180 | // C* and f* need shifting and masking, as shown by | |
b2a00c89 | 181 | // shiftright128[] and maskhigh128[] |
200359e8 | 182 | // 1 <= x <= 16 |
b2a00c89 | 183 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
184 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) |
185 | // the approximation of 10^(-x) was rounded up to 64 bits | |
b2a00c89 | 186 | __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); |
200359e8 L |
187 | |
188 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
189 | // if floor(C*) is even then C* = floor(C*) - logical right | |
190 | // shift; C* has p decimal digits, correct by Prop. 1) | |
191 | // else if floor(C*) is odd C* = floor(C*)-1 (logical right | |
192 | // shift; C* has p decimal digits, correct by Pr. 1) | |
193 | // else | |
194 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
195 | // correct by Property 1) | |
196 | // n = C* * 10^(e+x) | |
197 | ||
b2a00c89 | 198 | if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 |
200359e8 L |
199 | res = P128.w[1]; |
200 | fstar.w[1] = 0; | |
201 | fstar.w[0] = P128.w[0]; | |
b2a00c89 L |
202 | } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 |
203 | shift = shiftright128[ind - 1]; // 3 <= shift <= 63 | |
200359e8 | 204 | res = (P128.w[1] >> shift); |
b2a00c89 | 205 | fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 L |
206 | fstar.w[0] = P128.w[0]; |
207 | } | |
208 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
209 | // since round_to_even, subtract 1 if current result is odd | |
210 | if ((res & 0x0000000000000001ull) && (fstar.w[1] == 0) | |
b2a00c89 | 211 | && (fstar.w[0] < ten2mk64[ind - 1])) { |
200359e8 L |
212 | res--; |
213 | } | |
214 | // determine inexactness of the rounding of C* | |
215 | // if (0 < f* - 1/2 < 10^(-x)) then | |
216 | // the result is exact | |
217 | // else // if (f* - 1/2 > T*) then | |
218 | // the result is inexact | |
219 | if (ind - 1 <= 2) { | |
220 | if (fstar.w[0] > 0x8000000000000000ull) { | |
221 | // f* > 1/2 and the result may be exact | |
222 | // fstar.w[0] - 0x8000000000000000ull is f* - 1/2 | |
b2a00c89 | 223 | if ((fstar.w[0] - 0x8000000000000000ull) > ten2mk64[ind - 1]) { |
200359e8 L |
224 | // set the inexact flag |
225 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
226 | } // else the result is exact |
227 | } else { // the result is inexact; f2* <= 1/2 | |
200359e8 L |
228 | // set the inexact flag |
229 | *pfpsf |= INEXACT_EXCEPTION; | |
230 | } | |
b2a00c89 L |
231 | } else { // if 3 <= ind - 1 <= 21 |
232 | if (fstar.w[1] > onehalf128[ind - 1] || | |
233 | (fstar.w[1] == onehalf128[ind - 1] && fstar.w[0])) { | |
200359e8 L |
234 | // f2* > 1/2 and the result may be exact |
235 | // Calculate f2* - 1/2 | |
b2a00c89 L |
236 | if (fstar.w[1] > onehalf128[ind - 1] |
237 | || fstar.w[0] > ten2mk64[ind - 1]) { | |
200359e8 L |
238 | // set the inexact flag |
239 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
240 | } // else the result is exact |
241 | } else { // the result is inexact; f2* <= 1/2 | |
200359e8 L |
242 | // set the inexact flag |
243 | *pfpsf |= INEXACT_EXCEPTION; | |
244 | } | |
245 | } | |
246 | // set exponent to zero as it was negative before. | |
247 | res = x_sign | 0x31c0000000000000ull | res; | |
248 | BID_RETURN (res); | |
b2a00c89 | 249 | } else { // if exp < 0 and q + exp < 0 |
200359e8 L |
250 | // the result is +0 or -0 |
251 | res = x_sign | 0x31c0000000000000ull; | |
252 | *pfpsf |= INEXACT_EXCEPTION; | |
253 | BID_RETURN (res); | |
254 | } | |
255 | break; | |
256 | case ROUNDING_TIES_AWAY: | |
b2a00c89 | 257 | if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q |
200359e8 L |
258 | // need to shift right -exp digits from the coefficient; exp will be 0 |
259 | ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' | |
260 | // chop off ind digits from the lower part of C1 | |
261 | // C1 = C1 + 1/2 * 10^x where the result C1 fits in 64 bits | |
262 | // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate | |
b2a00c89 | 263 | C1 = C1 + midpoint64[ind - 1]; |
200359e8 L |
264 | // calculate C* and f* |
265 | // C* is actually floor(C*) in this case | |
266 | // C* and f* need shifting and masking, as shown by | |
b2a00c89 | 267 | // shiftright128[] and maskhigh128[] |
200359e8 | 268 | // 1 <= x <= 16 |
b2a00c89 | 269 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
270 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) |
271 | // the approximation of 10^(-x) was rounded up to 64 bits | |
b2a00c89 | 272 | __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); |
200359e8 L |
273 | |
274 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
275 | // C* = floor(C*) - logical right shift; C* has p decimal digits, | |
276 | // correct by Prop. 1) | |
277 | // else | |
278 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
279 | // correct by Property 1) | |
280 | // n = C* * 10^(e+x) | |
281 | ||
b2a00c89 | 282 | if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 |
200359e8 L |
283 | res = P128.w[1]; |
284 | fstar.w[1] = 0; | |
285 | fstar.w[0] = P128.w[0]; | |
b2a00c89 L |
286 | } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 |
287 | shift = shiftright128[ind - 1]; // 3 <= shift <= 63 | |
200359e8 | 288 | res = (P128.w[1] >> shift); |
b2a00c89 | 289 | fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 L |
290 | fstar.w[0] = P128.w[0]; |
291 | } | |
292 | // midpoints are already rounded correctly | |
293 | // determine inexactness of the rounding of C* | |
294 | // if (0 < f* - 1/2 < 10^(-x)) then | |
295 | // the result is exact | |
296 | // else // if (f* - 1/2 > T*) then | |
297 | // the result is inexact | |
298 | if (ind - 1 <= 2) { | |
299 | if (fstar.w[0] > 0x8000000000000000ull) { | |
300 | // f* > 1/2 and the result may be exact | |
301 | // fstar.w[0] - 0x8000000000000000ull is f* - 1/2 | |
b2a00c89 | 302 | if ((fstar.w[0] - 0x8000000000000000ull) > ten2mk64[ind - 1]) { |
200359e8 L |
303 | // set the inexact flag |
304 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
305 | } // else the result is exact |
306 | } else { // the result is inexact; f2* <= 1/2 | |
200359e8 L |
307 | // set the inexact flag |
308 | *pfpsf |= INEXACT_EXCEPTION; | |
309 | } | |
b2a00c89 L |
310 | } else { // if 3 <= ind - 1 <= 21 |
311 | if (fstar.w[1] > onehalf128[ind - 1] || | |
312 | (fstar.w[1] == onehalf128[ind - 1] && fstar.w[0])) { | |
200359e8 L |
313 | // f2* > 1/2 and the result may be exact |
314 | // Calculate f2* - 1/2 | |
b2a00c89 L |
315 | if (fstar.w[1] > onehalf128[ind - 1] |
316 | || fstar.w[0] > ten2mk64[ind - 1]) { | |
200359e8 L |
317 | // set the inexact flag |
318 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
319 | } // else the result is exact |
320 | } else { // the result is inexact; f2* <= 1/2 | |
200359e8 L |
321 | // set the inexact flag |
322 | *pfpsf |= INEXACT_EXCEPTION; | |
323 | } | |
324 | } | |
325 | // set exponent to zero as it was negative before. | |
326 | res = x_sign | 0x31c0000000000000ull | res; | |
327 | BID_RETURN (res); | |
b2a00c89 | 328 | } else { // if exp < 0 and q + exp < 0 |
200359e8 L |
329 | // the result is +0 or -0 |
330 | res = x_sign | 0x31c0000000000000ull; | |
331 | *pfpsf |= INEXACT_EXCEPTION; | |
332 | BID_RETURN (res); | |
333 | } | |
334 | break; | |
335 | case ROUNDING_DOWN: | |
b2a00c89 | 336 | if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q |
200359e8 L |
337 | // need to shift right -exp digits from the coefficient; exp will be 0 |
338 | ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' | |
339 | // chop off ind digits from the lower part of C1 | |
340 | // C1 fits in 64 bits | |
341 | // calculate C* and f* | |
342 | // C* is actually floor(C*) in this case | |
343 | // C* and f* need shifting and masking, as shown by | |
b2a00c89 | 344 | // shiftright128[] and maskhigh128[] |
200359e8 | 345 | // 1 <= x <= 16 |
b2a00c89 | 346 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
347 | // C* = C1 * 10^(-x) |
348 | // the approximation of 10^(-x) was rounded up to 64 bits | |
b2a00c89 | 349 | __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); |
200359e8 L |
350 | |
351 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
352 | // correct by Property 1) | |
353 | // if (0 < f* < 10^(-x)) then the result is exact | |
354 | // n = C* * 10^(e+x) | |
355 | ||
b2a00c89 | 356 | if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 |
200359e8 L |
357 | res = P128.w[1]; |
358 | fstar.w[1] = 0; | |
359 | fstar.w[0] = P128.w[0]; | |
b2a00c89 L |
360 | } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 |
361 | shift = shiftright128[ind - 1]; // 3 <= shift <= 63 | |
200359e8 | 362 | res = (P128.w[1] >> shift); |
b2a00c89 | 363 | fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 L |
364 | fstar.w[0] = P128.w[0]; |
365 | } | |
366 | // if (f* > 10^(-x)) then the result is inexact | |
b2a00c89 | 367 | if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1])) { |
200359e8 L |
368 | if (x_sign) { |
369 | // if negative and not exact, increment magnitude | |
370 | res++; | |
371 | } | |
372 | *pfpsf |= INEXACT_EXCEPTION; | |
373 | } | |
374 | // set exponent to zero as it was negative before. | |
375 | res = x_sign | 0x31c0000000000000ull | res; | |
376 | BID_RETURN (res); | |
b2a00c89 | 377 | } else { // if exp < 0 and q + exp <= 0 |
200359e8 L |
378 | // the result is +0 or -1 |
379 | if (x_sign) { | |
380 | res = 0xb1c0000000000001ull; | |
381 | } else { | |
382 | res = 0x31c0000000000000ull; | |
383 | } | |
384 | *pfpsf |= INEXACT_EXCEPTION; | |
385 | BID_RETURN (res); | |
386 | } | |
387 | break; | |
388 | case ROUNDING_UP: | |
b2a00c89 | 389 | if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q |
200359e8 L |
390 | // need to shift right -exp digits from the coefficient; exp will be 0 |
391 | ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' | |
392 | // chop off ind digits from the lower part of C1 | |
393 | // C1 fits in 64 bits | |
394 | // calculate C* and f* | |
395 | // C* is actually floor(C*) in this case | |
396 | // C* and f* need shifting and masking, as shown by | |
b2a00c89 | 397 | // shiftright128[] and maskhigh128[] |
200359e8 | 398 | // 1 <= x <= 16 |
b2a00c89 | 399 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
400 | // C* = C1 * 10^(-x) |
401 | // the approximation of 10^(-x) was rounded up to 64 bits | |
b2a00c89 | 402 | __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); |
200359e8 L |
403 | |
404 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
405 | // correct by Property 1) | |
406 | // if (0 < f* < 10^(-x)) then the result is exact | |
407 | // n = C* * 10^(e+x) | |
408 | ||
b2a00c89 | 409 | if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 |
200359e8 L |
410 | res = P128.w[1]; |
411 | fstar.w[1] = 0; | |
412 | fstar.w[0] = P128.w[0]; | |
b2a00c89 L |
413 | } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 |
414 | shift = shiftright128[ind - 1]; // 3 <= shift <= 63 | |
200359e8 | 415 | res = (P128.w[1] >> shift); |
b2a00c89 | 416 | fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 L |
417 | fstar.w[0] = P128.w[0]; |
418 | } | |
419 | // if (f* > 10^(-x)) then the result is inexact | |
b2a00c89 | 420 | if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1])) { |
200359e8 L |
421 | if (!x_sign) { |
422 | // if positive and not exact, increment magnitude | |
423 | res++; | |
424 | } | |
425 | *pfpsf |= INEXACT_EXCEPTION; | |
426 | } | |
427 | // set exponent to zero as it was negative before. | |
428 | res = x_sign | 0x31c0000000000000ull | res; | |
429 | BID_RETURN (res); | |
b2a00c89 | 430 | } else { // if exp < 0 and q + exp <= 0 |
200359e8 L |
431 | // the result is -0 or +1 |
432 | if (x_sign) { | |
433 | res = 0xb1c0000000000000ull; | |
434 | } else { | |
435 | res = 0x31c0000000000001ull; | |
436 | } | |
437 | *pfpsf |= INEXACT_EXCEPTION; | |
438 | BID_RETURN (res); | |
439 | } | |
440 | break; | |
441 | case ROUNDING_TO_ZERO: | |
b2a00c89 | 442 | if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q |
200359e8 L |
443 | // need to shift right -exp digits from the coefficient; exp will be 0 |
444 | ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' | |
445 | // chop off ind digits from the lower part of C1 | |
446 | // C1 fits in 127 bits | |
447 | // calculate C* and f* | |
448 | // C* is actually floor(C*) in this case | |
449 | // C* and f* need shifting and masking, as shown by | |
b2a00c89 | 450 | // shiftright128[] and maskhigh128[] |
200359e8 | 451 | // 1 <= x <= 16 |
b2a00c89 | 452 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
453 | // C* = C1 * 10^(-x) |
454 | // the approximation of 10^(-x) was rounded up to 64 bits | |
b2a00c89 | 455 | __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); |
200359e8 L |
456 | |
457 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
458 | // correct by Property 1) | |
459 | // if (0 < f* < 10^(-x)) then the result is exact | |
460 | // n = C* * 10^(e+x) | |
461 | ||
b2a00c89 | 462 | if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 |
200359e8 L |
463 | res = P128.w[1]; |
464 | fstar.w[1] = 0; | |
465 | fstar.w[0] = P128.w[0]; | |
b2a00c89 L |
466 | } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 |
467 | shift = shiftright128[ind - 1]; // 3 <= shift <= 63 | |
200359e8 | 468 | res = (P128.w[1] >> shift); |
b2a00c89 | 469 | fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 L |
470 | fstar.w[0] = P128.w[0]; |
471 | } | |
472 | // if (f* > 10^(-x)) then the result is inexact | |
b2a00c89 | 473 | if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1])) { |
200359e8 L |
474 | *pfpsf |= INEXACT_EXCEPTION; |
475 | } | |
476 | // set exponent to zero as it was negative before. | |
477 | res = x_sign | 0x31c0000000000000ull | res; | |
478 | BID_RETURN (res); | |
b2a00c89 | 479 | } else { // if exp < 0 and q + exp < 0 |
200359e8 L |
480 | // the result is +0 or -0 |
481 | res = x_sign | 0x31c0000000000000ull; | |
482 | *pfpsf |= INEXACT_EXCEPTION; | |
483 | BID_RETURN (res); | |
484 | } | |
485 | break; | |
b2a00c89 | 486 | } // end switch () |
200359e8 L |
487 | BID_RETURN (res); |
488 | } | |
489 | ||
490 | /***************************************************************************** | |
491 | * BID64_round_integral_nearest_even | |
492 | ****************************************************************************/ | |
493 | ||
494 | #if DECIMAL_CALL_BY_REFERENCE | |
495 | void | |
b2a00c89 | 496 | bid64_round_integral_nearest_even (UINT64 * pres, |
200359e8 L |
497 | UINT64 * |
498 | px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
499 | _EXC_INFO_PARAM) { | |
500 | UINT64 x = *px; | |
501 | #else | |
502 | UINT64 | |
b2a00c89 | 503 | bid64_round_integral_nearest_even (UINT64 x _EXC_FLAGS_PARAM |
200359e8 L |
504 | _EXC_MASKS_PARAM _EXC_INFO_PARAM) { |
505 | #endif | |
506 | ||
b2a00c89 | 507 | UINT64 res = 0xbaddbaddbaddbaddull; |
200359e8 | 508 | UINT64 x_sign; |
b2a00c89 | 509 | int exp; // unbiased exponent |
200359e8 L |
510 | // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) |
511 | BID_UI64DOUBLE tmp1; | |
512 | int x_nr_bits; | |
513 | int q, ind, shift; | |
514 | UINT64 C1; | |
515 | UINT128 fstar; | |
516 | UINT128 P128; | |
517 | ||
b2a00c89 L |
518 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
519 | ||
520 | // check for NaNs and infinities | |
521 | if ((x & MASK_NAN) == MASK_NAN) { // check for NaN | |
522 | if ((x & 0x0003ffffffffffffull) > 999999999999999ull) | |
523 | x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits | |
524 | else | |
525 | x = x & 0xfe03ffffffffffffull; // clear G6-G12 | |
526 | if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN | |
527 | // set invalid flag | |
200359e8 | 528 | *pfpsf |= INVALID_EXCEPTION; |
b2a00c89 | 529 | // return quiet (SNaN) |
200359e8 | 530 | res = x & 0xfdffffffffffffffull; |
b2a00c89 L |
531 | } else { // QNaN |
532 | res = x; | |
200359e8 | 533 | } |
b2a00c89 L |
534 | BID_RETURN (res); |
535 | } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity | |
536 | res = x_sign | 0x7800000000000000ull; | |
200359e8 L |
537 | BID_RETURN (res); |
538 | } | |
539 | // unpack x | |
200359e8 L |
540 | if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { |
541 | // if the steering bits are 11 (condition will be 0), then | |
542 | // the exponent is G[0:w+1] | |
543 | exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398; | |
544 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; | |
b2a00c89 | 545 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
546 | C1 = 0; |
547 | } | |
b2a00c89 | 548 | } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) |
200359e8 L |
549 | exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398; |
550 | C1 = (x & MASK_BINARY_SIG1); | |
551 | } | |
552 | ||
553 | // if x is 0 or non-canonical | |
554 | if (C1 == 0) { | |
555 | if (exp < 0) | |
556 | exp = 0; | |
557 | res = x_sign | (((UINT64) exp + 398) << 53); | |
558 | BID_RETURN (res); | |
559 | } | |
560 | // x is a finite non-zero number (not 0, non-canonical, or special) | |
561 | ||
562 | // return 0 if (exp <= -(p+1)) | |
563 | if (exp <= -17) { | |
564 | res = x_sign | 0x31c0000000000000ull; | |
565 | BID_RETURN (res); | |
566 | } | |
567 | // q = nr. of decimal digits in x (1 <= q <= 54) | |
568 | // determine first the nr. of bits in x | |
b2a00c89 | 569 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 570 | q = 16; |
b2a00c89 | 571 | } else { // if x < 2^53 |
200359e8 L |
572 | tmp1.d = (double) C1; // exact conversion |
573 | x_nr_bits = | |
574 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 | 575 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 576 | if (q == 0) { |
b2a00c89 L |
577 | q = nr_digits[x_nr_bits - 1].digits1; |
578 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
579 | q++; |
580 | } | |
581 | } | |
582 | ||
b2a00c89 | 583 | if (exp >= 0) { // -exp <= 0 |
200359e8 L |
584 | // the argument is an integer already |
585 | res = x; | |
586 | BID_RETURN (res); | |
b2a00c89 | 587 | } else if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q |
200359e8 L |
588 | // need to shift right -exp digits from the coefficient; the exp will be 0 |
589 | ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' | |
590 | // chop off ind digits from the lower part of C1 | |
591 | // C1 = C1 + 1/2 * 10^x where the result C1 fits in 64 bits | |
592 | // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate | |
b2a00c89 | 593 | C1 = C1 + midpoint64[ind - 1]; |
200359e8 L |
594 | // calculate C* and f* |
595 | // C* is actually floor(C*) in this case | |
596 | // C* and f* need shifting and masking, as shown by | |
b2a00c89 | 597 | // shiftright128[] and maskhigh128[] |
200359e8 | 598 | // 1 <= x <= 16 |
b2a00c89 | 599 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
600 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) |
601 | // the approximation of 10^(-x) was rounded up to 64 bits | |
b2a00c89 | 602 | __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); |
200359e8 L |
603 | |
604 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
605 | // if floor(C*) is even then C* = floor(C*) - logical right | |
606 | // shift; C* has p decimal digits, correct by Prop. 1) | |
607 | // else if floor(C*) is odd C* = floor(C*)-1 (logical right | |
608 | // shift; C* has p decimal digits, correct by Pr. 1) | |
609 | // else | |
610 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
611 | // correct by Property 1) | |
612 | // n = C* * 10^(e+x) | |
613 | ||
b2a00c89 | 614 | if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 |
200359e8 L |
615 | res = P128.w[1]; |
616 | fstar.w[1] = 0; | |
617 | fstar.w[0] = P128.w[0]; | |
b2a00c89 L |
618 | } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 |
619 | shift = shiftright128[ind - 1]; // 3 <= shift <= 63 | |
200359e8 | 620 | res = (P128.w[1] >> shift); |
b2a00c89 | 621 | fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 L |
622 | fstar.w[0] = P128.w[0]; |
623 | } | |
624 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
625 | // since round_to_even, subtract 1 if current result is odd | |
626 | if ((res & 0x0000000000000001ull) && (fstar.w[1] == 0) | |
b2a00c89 | 627 | && (fstar.w[0] < ten2mk64[ind - 1])) { |
200359e8 L |
628 | res--; |
629 | } | |
630 | // set exponent to zero as it was negative before. | |
631 | res = x_sign | 0x31c0000000000000ull | res; | |
632 | BID_RETURN (res); | |
b2a00c89 | 633 | } else { // if exp < 0 and q + exp < 0 |
200359e8 L |
634 | // the result is +0 or -0 |
635 | res = x_sign | 0x31c0000000000000ull; | |
636 | BID_RETURN (res); | |
637 | } | |
638 | } | |
639 | ||
640 | /***************************************************************************** | |
641 | * BID64_round_integral_negative | |
642 | *****************************************************************************/ | |
643 | ||
644 | #if DECIMAL_CALL_BY_REFERENCE | |
645 | void | |
b2a00c89 | 646 | bid64_round_integral_negative (UINT64 * pres, |
200359e8 L |
647 | UINT64 * |
648 | px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
649 | _EXC_INFO_PARAM) { | |
650 | UINT64 x = *px; | |
651 | #else | |
652 | UINT64 | |
b2a00c89 | 653 | bid64_round_integral_negative (UINT64 x _EXC_FLAGS_PARAM |
200359e8 L |
654 | _EXC_MASKS_PARAM _EXC_INFO_PARAM) { |
655 | #endif | |
656 | ||
b2a00c89 | 657 | UINT64 res = 0xbaddbaddbaddbaddull; |
200359e8 | 658 | UINT64 x_sign; |
b2a00c89 | 659 | int exp; // unbiased exponent |
200359e8 L |
660 | // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) |
661 | BID_UI64DOUBLE tmp1; | |
662 | int x_nr_bits; | |
663 | int q, ind, shift; | |
664 | UINT64 C1; | |
665 | // UINT64 res is C* at first - represents up to 34 decimal digits ~ 113 bits | |
666 | UINT128 fstar; | |
667 | UINT128 P128; | |
668 | ||
b2a00c89 L |
669 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
670 | ||
671 | // check for NaNs and infinities | |
672 | if ((x & MASK_NAN) == MASK_NAN) { // check for NaN | |
673 | if ((x & 0x0003ffffffffffffull) > 999999999999999ull) | |
674 | x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits | |
675 | else | |
676 | x = x & 0xfe03ffffffffffffull; // clear G6-G12 | |
677 | if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN | |
678 | // set invalid flag | |
200359e8 | 679 | *pfpsf |= INVALID_EXCEPTION; |
b2a00c89 | 680 | // return quiet (SNaN) |
200359e8 | 681 | res = x & 0xfdffffffffffffffull; |
b2a00c89 L |
682 | } else { // QNaN |
683 | res = x; | |
200359e8 | 684 | } |
b2a00c89 L |
685 | BID_RETURN (res); |
686 | } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity | |
687 | res = x_sign | 0x7800000000000000ull; | |
200359e8 L |
688 | BID_RETURN (res); |
689 | } | |
690 | // unpack x | |
200359e8 L |
691 | if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { |
692 | // if the steering bits are 11 (condition will be 0), then | |
693 | // the exponent is G[0:w+1] | |
694 | exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398; | |
695 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; | |
b2a00c89 | 696 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
697 | C1 = 0; |
698 | } | |
b2a00c89 | 699 | } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) |
200359e8 L |
700 | exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398; |
701 | C1 = (x & MASK_BINARY_SIG1); | |
702 | } | |
703 | ||
704 | // if x is 0 or non-canonical | |
705 | if (C1 == 0) { | |
706 | if (exp < 0) | |
707 | exp = 0; | |
708 | res = x_sign | (((UINT64) exp + 398) << 53); | |
709 | BID_RETURN (res); | |
710 | } | |
711 | // x is a finite non-zero number (not 0, non-canonical, or special) | |
712 | ||
713 | // return 0 if (exp <= -p) | |
714 | if (exp <= -16) { | |
715 | if (x_sign) { | |
716 | res = 0xb1c0000000000001ull; | |
717 | } else { | |
718 | res = 0x31c0000000000000ull; | |
719 | } | |
720 | BID_RETURN (res); | |
721 | } | |
722 | // q = nr. of decimal digits in x (1 <= q <= 54) | |
723 | // determine first the nr. of bits in x | |
b2a00c89 | 724 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 725 | q = 16; |
b2a00c89 | 726 | } else { // if x < 2^53 |
200359e8 L |
727 | tmp1.d = (double) C1; // exact conversion |
728 | x_nr_bits = | |
729 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 | 730 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 731 | if (q == 0) { |
b2a00c89 L |
732 | q = nr_digits[x_nr_bits - 1].digits1; |
733 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
734 | q++; |
735 | } | |
736 | } | |
737 | ||
b2a00c89 | 738 | if (exp >= 0) { // -exp <= 0 |
200359e8 L |
739 | // the argument is an integer already |
740 | res = x; | |
741 | BID_RETURN (res); | |
b2a00c89 | 742 | } else if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q |
200359e8 L |
743 | // need to shift right -exp digits from the coefficient; the exp will be 0 |
744 | ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' | |
745 | // chop off ind digits from the lower part of C1 | |
746 | // C1 fits in 64 bits | |
747 | // calculate C* and f* | |
748 | // C* is actually floor(C*) in this case | |
749 | // C* and f* need shifting and masking, as shown by | |
b2a00c89 | 750 | // shiftright128[] and maskhigh128[] |
200359e8 | 751 | // 1 <= x <= 16 |
b2a00c89 | 752 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
753 | // C* = C1 * 10^(-x) |
754 | // the approximation of 10^(-x) was rounded up to 64 bits | |
b2a00c89 | 755 | __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); |
200359e8 L |
756 | |
757 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
758 | // correct by Property 1) | |
759 | // if (0 < f* < 10^(-x)) then the result is exact | |
760 | // n = C* * 10^(e+x) | |
761 | ||
b2a00c89 | 762 | if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 |
200359e8 L |
763 | res = P128.w[1]; |
764 | fstar.w[1] = 0; | |
765 | fstar.w[0] = P128.w[0]; | |
b2a00c89 L |
766 | } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 |
767 | shift = shiftright128[ind - 1]; // 3 <= shift <= 63 | |
200359e8 | 768 | res = (P128.w[1] >> shift); |
b2a00c89 | 769 | fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 L |
770 | fstar.w[0] = P128.w[0]; |
771 | } | |
772 | // if (f* > 10^(-x)) then the result is inexact | |
773 | if (x_sign | |
b2a00c89 | 774 | && ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1]))) { |
200359e8 L |
775 | // if negative and not exact, increment magnitude |
776 | res++; | |
777 | } | |
778 | // set exponent to zero as it was negative before. | |
779 | res = x_sign | 0x31c0000000000000ull | res; | |
780 | BID_RETURN (res); | |
b2a00c89 | 781 | } else { // if exp < 0 and q + exp <= 0 |
200359e8 L |
782 | // the result is +0 or -1 |
783 | if (x_sign) { | |
784 | res = 0xb1c0000000000001ull; | |
785 | } else { | |
786 | res = 0x31c0000000000000ull; | |
787 | } | |
788 | BID_RETURN (res); | |
789 | } | |
790 | } | |
791 | ||
792 | /***************************************************************************** | |
793 | * BID64_round_integral_positive | |
794 | ****************************************************************************/ | |
795 | ||
796 | #if DECIMAL_CALL_BY_REFERENCE | |
797 | void | |
b2a00c89 | 798 | bid64_round_integral_positive (UINT64 * pres, |
200359e8 L |
799 | UINT64 * |
800 | px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
801 | _EXC_INFO_PARAM) { | |
802 | UINT64 x = *px; | |
803 | #else | |
804 | UINT64 | |
b2a00c89 | 805 | bid64_round_integral_positive (UINT64 x _EXC_FLAGS_PARAM |
200359e8 L |
806 | _EXC_MASKS_PARAM _EXC_INFO_PARAM) { |
807 | #endif | |
808 | ||
b2a00c89 | 809 | UINT64 res = 0xbaddbaddbaddbaddull; |
200359e8 | 810 | UINT64 x_sign; |
b2a00c89 | 811 | int exp; // unbiased exponent |
200359e8 L |
812 | // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) |
813 | BID_UI64DOUBLE tmp1; | |
814 | int x_nr_bits; | |
815 | int q, ind, shift; | |
816 | UINT64 C1; | |
817 | // UINT64 res is C* at first - represents up to 34 decimal digits ~ 113 bits | |
818 | UINT128 fstar; | |
819 | UINT128 P128; | |
820 | ||
b2a00c89 L |
821 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
822 | ||
823 | // check for NaNs and infinities | |
824 | if ((x & MASK_NAN) == MASK_NAN) { // check for NaN | |
825 | if ((x & 0x0003ffffffffffffull) > 999999999999999ull) | |
826 | x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits | |
827 | else | |
828 | x = x & 0xfe03ffffffffffffull; // clear G6-G12 | |
829 | if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN | |
830 | // set invalid flag | |
200359e8 | 831 | *pfpsf |= INVALID_EXCEPTION; |
b2a00c89 | 832 | // return quiet (SNaN) |
200359e8 | 833 | res = x & 0xfdffffffffffffffull; |
b2a00c89 L |
834 | } else { // QNaN |
835 | res = x; | |
200359e8 | 836 | } |
b2a00c89 L |
837 | BID_RETURN (res); |
838 | } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity | |
839 | res = x_sign | 0x7800000000000000ull; | |
200359e8 L |
840 | BID_RETURN (res); |
841 | } | |
842 | // unpack x | |
200359e8 L |
843 | if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { |
844 | // if the steering bits are 11 (condition will be 0), then | |
845 | // the exponent is G[0:w+1] | |
846 | exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398; | |
847 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; | |
b2a00c89 | 848 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
849 | C1 = 0; |
850 | } | |
b2a00c89 | 851 | } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) |
200359e8 L |
852 | exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398; |
853 | C1 = (x & MASK_BINARY_SIG1); | |
854 | } | |
855 | ||
856 | // if x is 0 or non-canonical | |
857 | if (C1 == 0) { | |
858 | if (exp < 0) | |
859 | exp = 0; | |
860 | res = x_sign | (((UINT64) exp + 398) << 53); | |
861 | BID_RETURN (res); | |
862 | } | |
863 | // x is a finite non-zero number (not 0, non-canonical, or special) | |
864 | ||
865 | // return 0 if (exp <= -p) | |
866 | if (exp <= -16) { | |
867 | if (x_sign) { | |
868 | res = 0xb1c0000000000000ull; | |
869 | } else { | |
870 | res = 0x31c0000000000001ull; | |
871 | } | |
872 | BID_RETURN (res); | |
873 | } | |
874 | // q = nr. of decimal digits in x (1 <= q <= 54) | |
875 | // determine first the nr. of bits in x | |
b2a00c89 | 876 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 877 | q = 16; |
b2a00c89 | 878 | } else { // if x < 2^53 |
200359e8 L |
879 | tmp1.d = (double) C1; // exact conversion |
880 | x_nr_bits = | |
881 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 | 882 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 883 | if (q == 0) { |
b2a00c89 L |
884 | q = nr_digits[x_nr_bits - 1].digits1; |
885 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
886 | q++; |
887 | } | |
888 | } | |
889 | ||
b2a00c89 | 890 | if (exp >= 0) { // -exp <= 0 |
200359e8 L |
891 | // the argument is an integer already |
892 | res = x; | |
893 | BID_RETURN (res); | |
b2a00c89 | 894 | } else if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q |
200359e8 L |
895 | // need to shift right -exp digits from the coefficient; the exp will be 0 |
896 | ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' | |
897 | // chop off ind digits from the lower part of C1 | |
898 | // C1 fits in 64 bits | |
899 | // calculate C* and f* | |
900 | // C* is actually floor(C*) in this case | |
901 | // C* and f* need shifting and masking, as shown by | |
b2a00c89 | 902 | // shiftright128[] and maskhigh128[] |
200359e8 | 903 | // 1 <= x <= 16 |
b2a00c89 | 904 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
905 | // C* = C1 * 10^(-x) |
906 | // the approximation of 10^(-x) was rounded up to 64 bits | |
b2a00c89 | 907 | __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); |
200359e8 L |
908 | |
909 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
910 | // correct by Property 1) | |
911 | // if (0 < f* < 10^(-x)) then the result is exact | |
912 | // n = C* * 10^(e+x) | |
913 | ||
b2a00c89 | 914 | if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 |
200359e8 L |
915 | res = P128.w[1]; |
916 | fstar.w[1] = 0; | |
917 | fstar.w[0] = P128.w[0]; | |
b2a00c89 L |
918 | } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 |
919 | shift = shiftright128[ind - 1]; // 3 <= shift <= 63 | |
200359e8 | 920 | res = (P128.w[1] >> shift); |
b2a00c89 | 921 | fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 L |
922 | fstar.w[0] = P128.w[0]; |
923 | } | |
924 | // if (f* > 10^(-x)) then the result is inexact | |
925 | if (!x_sign | |
b2a00c89 | 926 | && ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1]))) { |
200359e8 L |
927 | // if positive and not exact, increment magnitude |
928 | res++; | |
929 | } | |
930 | // set exponent to zero as it was negative before. | |
931 | res = x_sign | 0x31c0000000000000ull | res; | |
932 | BID_RETURN (res); | |
b2a00c89 | 933 | } else { // if exp < 0 and q + exp <= 0 |
200359e8 L |
934 | // the result is -0 or +1 |
935 | if (x_sign) { | |
936 | res = 0xb1c0000000000000ull; | |
937 | } else { | |
938 | res = 0x31c0000000000001ull; | |
939 | } | |
940 | BID_RETURN (res); | |
941 | } | |
942 | } | |
943 | ||
944 | /***************************************************************************** | |
945 | * BID64_round_integral_zero | |
946 | ****************************************************************************/ | |
947 | ||
948 | #if DECIMAL_CALL_BY_REFERENCE | |
949 | void | |
b2a00c89 | 950 | bid64_round_integral_zero (UINT64 * pres, |
200359e8 L |
951 | UINT64 * |
952 | px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
953 | _EXC_INFO_PARAM) { | |
954 | UINT64 x = *px; | |
955 | #else | |
956 | UINT64 | |
b2a00c89 | 957 | bid64_round_integral_zero (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
200359e8 L |
958 | _EXC_INFO_PARAM) { |
959 | #endif | |
960 | ||
b2a00c89 | 961 | UINT64 res = 0xbaddbaddbaddbaddull; |
200359e8 | 962 | UINT64 x_sign; |
b2a00c89 | 963 | int exp; // unbiased exponent |
200359e8 L |
964 | // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) |
965 | BID_UI64DOUBLE tmp1; | |
966 | int x_nr_bits; | |
967 | int q, ind, shift; | |
968 | UINT64 C1; | |
969 | // UINT64 res is C* at first - represents up to 34 decimal digits ~ 113 bits | |
970 | UINT128 P128; | |
971 | ||
b2a00c89 L |
972 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
973 | ||
974 | // check for NaNs and infinities | |
975 | if ((x & MASK_NAN) == MASK_NAN) { // check for NaN | |
976 | if ((x & 0x0003ffffffffffffull) > 999999999999999ull) | |
977 | x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits | |
978 | else | |
979 | x = x & 0xfe03ffffffffffffull; // clear G6-G12 | |
980 | if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN | |
981 | // set invalid flag | |
200359e8 | 982 | *pfpsf |= INVALID_EXCEPTION; |
b2a00c89 | 983 | // return quiet (SNaN) |
200359e8 | 984 | res = x & 0xfdffffffffffffffull; |
b2a00c89 L |
985 | } else { // QNaN |
986 | res = x; | |
200359e8 | 987 | } |
b2a00c89 L |
988 | BID_RETURN (res); |
989 | } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity | |
990 | res = x_sign | 0x7800000000000000ull; | |
200359e8 L |
991 | BID_RETURN (res); |
992 | } | |
993 | // unpack x | |
200359e8 L |
994 | if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { |
995 | // if the steering bits are 11 (condition will be 0), then | |
996 | // the exponent is G[0:w+1] | |
997 | exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398; | |
998 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; | |
b2a00c89 | 999 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
1000 | C1 = 0; |
1001 | } | |
b2a00c89 | 1002 | } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) |
200359e8 L |
1003 | exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398; |
1004 | C1 = (x & MASK_BINARY_SIG1); | |
1005 | } | |
1006 | ||
1007 | // if x is 0 or non-canonical | |
1008 | if (C1 == 0) { | |
1009 | if (exp < 0) | |
1010 | exp = 0; | |
1011 | res = x_sign | (((UINT64) exp + 398) << 53); | |
1012 | BID_RETURN (res); | |
1013 | } | |
1014 | // x is a finite non-zero number (not 0, non-canonical, or special) | |
1015 | ||
1016 | // return 0 if (exp <= -p) | |
1017 | if (exp <= -16) { | |
1018 | res = x_sign | 0x31c0000000000000ull; | |
1019 | BID_RETURN (res); | |
1020 | } | |
1021 | // q = nr. of decimal digits in x (1 <= q <= 54) | |
1022 | // determine first the nr. of bits in x | |
b2a00c89 | 1023 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 1024 | q = 16; |
b2a00c89 | 1025 | } else { // if x < 2^53 |
200359e8 L |
1026 | tmp1.d = (double) C1; // exact conversion |
1027 | x_nr_bits = | |
1028 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 | 1029 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 1030 | if (q == 0) { |
b2a00c89 L |
1031 | q = nr_digits[x_nr_bits - 1].digits1; |
1032 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
1033 | q++; |
1034 | } | |
1035 | } | |
1036 | ||
b2a00c89 | 1037 | if (exp >= 0) { // -exp <= 0 |
200359e8 L |
1038 | // the argument is an integer already |
1039 | res = x; | |
1040 | BID_RETURN (res); | |
b2a00c89 | 1041 | } else if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q |
200359e8 L |
1042 | // need to shift right -exp digits from the coefficient; the exp will be 0 |
1043 | ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' | |
1044 | // chop off ind digits from the lower part of C1 | |
1045 | // C1 fits in 127 bits | |
1046 | // calculate C* and f* | |
1047 | // C* is actually floor(C*) in this case | |
1048 | // C* and f* need shifting and masking, as shown by | |
b2a00c89 | 1049 | // shiftright128[] and maskhigh128[] |
200359e8 | 1050 | // 1 <= x <= 16 |
b2a00c89 | 1051 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
1052 | // C* = C1 * 10^(-x) |
1053 | // the approximation of 10^(-x) was rounded up to 64 bits | |
b2a00c89 | 1054 | __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); |
200359e8 L |
1055 | |
1056 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
1057 | // correct by Property 1) | |
1058 | // if (0 < f* < 10^(-x)) then the result is exact | |
1059 | // n = C* * 10^(e+x) | |
1060 | ||
b2a00c89 | 1061 | if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 |
200359e8 L |
1062 | res = P128.w[1]; |
1063 | // redundant fstar.w[1] = 0; | |
1064 | // redundant fstar.w[0] = P128.w[0]; | |
b2a00c89 L |
1065 | } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 |
1066 | shift = shiftright128[ind - 1]; // 3 <= shift <= 63 | |
200359e8 | 1067 | res = (P128.w[1] >> shift); |
b2a00c89 | 1068 | // redundant fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 L |
1069 | // redundant fstar.w[0] = P128.w[0]; |
1070 | } | |
1071 | // if (f* > 10^(-x)) then the result is inexact | |
b2a00c89 | 1072 | // if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind-1])){ |
200359e8 L |
1073 | // // redundant |
1074 | // } | |
1075 | // set exponent to zero as it was negative before. | |
1076 | res = x_sign | 0x31c0000000000000ull | res; | |
1077 | BID_RETURN (res); | |
b2a00c89 | 1078 | } else { // if exp < 0 and q + exp < 0 |
200359e8 L |
1079 | // the result is +0 or -0 |
1080 | res = x_sign | 0x31c0000000000000ull; | |
1081 | BID_RETURN (res); | |
1082 | } | |
1083 | } | |
1084 | ||
1085 | /***************************************************************************** | |
1086 | * BID64_round_integral_nearest_away | |
1087 | ****************************************************************************/ | |
1088 | ||
1089 | #if DECIMAL_CALL_BY_REFERENCE | |
1090 | void | |
b2a00c89 | 1091 | bid64_round_integral_nearest_away (UINT64 * pres, |
200359e8 L |
1092 | UINT64 * |
1093 | px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM | |
1094 | _EXC_INFO_PARAM) { | |
1095 | UINT64 x = *px; | |
1096 | #else | |
1097 | UINT64 | |
b2a00c89 | 1098 | bid64_round_integral_nearest_away (UINT64 x _EXC_FLAGS_PARAM |
200359e8 L |
1099 | _EXC_MASKS_PARAM _EXC_INFO_PARAM) { |
1100 | #endif | |
1101 | ||
b2a00c89 | 1102 | UINT64 res = 0xbaddbaddbaddbaddull; |
200359e8 | 1103 | UINT64 x_sign; |
b2a00c89 | 1104 | int exp; // unbiased exponent |
200359e8 L |
1105 | // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) |
1106 | BID_UI64DOUBLE tmp1; | |
1107 | int x_nr_bits; | |
1108 | int q, ind, shift; | |
1109 | UINT64 C1; | |
1110 | UINT128 P128; | |
1111 | ||
b2a00c89 L |
1112 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
1113 | ||
1114 | // check for NaNs and infinities | |
1115 | if ((x & MASK_NAN) == MASK_NAN) { // check for NaN | |
1116 | if ((x & 0x0003ffffffffffffull) > 999999999999999ull) | |
1117 | x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits | |
1118 | else | |
1119 | x = x & 0xfe03ffffffffffffull; // clear G6-G12 | |
1120 | if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN | |
1121 | // set invalid flag | |
200359e8 | 1122 | *pfpsf |= INVALID_EXCEPTION; |
b2a00c89 | 1123 | // return quiet (SNaN) |
200359e8 | 1124 | res = x & 0xfdffffffffffffffull; |
b2a00c89 L |
1125 | } else { // QNaN |
1126 | res = x; | |
200359e8 | 1127 | } |
b2a00c89 L |
1128 | BID_RETURN (res); |
1129 | } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity | |
1130 | res = x_sign | 0x7800000000000000ull; | |
200359e8 L |
1131 | BID_RETURN (res); |
1132 | } | |
1133 | // unpack x | |
200359e8 L |
1134 | if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { |
1135 | // if the steering bits are 11 (condition will be 0), then | |
1136 | // the exponent is G[0:w+1] | |
1137 | exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398; | |
1138 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; | |
b2a00c89 | 1139 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
1140 | C1 = 0; |
1141 | } | |
b2a00c89 | 1142 | } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) |
200359e8 L |
1143 | exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398; |
1144 | C1 = (x & MASK_BINARY_SIG1); | |
1145 | } | |
1146 | ||
1147 | // if x is 0 or non-canonical | |
1148 | if (C1 == 0) { | |
1149 | if (exp < 0) | |
1150 | exp = 0; | |
1151 | res = x_sign | (((UINT64) exp + 398) << 53); | |
1152 | BID_RETURN (res); | |
1153 | } | |
1154 | // x is a finite non-zero number (not 0, non-canonical, or special) | |
1155 | ||
1156 | // return 0 if (exp <= -(p+1)) | |
1157 | if (exp <= -17) { | |
1158 | res = x_sign | 0x31c0000000000000ull; | |
1159 | BID_RETURN (res); | |
1160 | } | |
1161 | // q = nr. of decimal digits in x (1 <= q <= 54) | |
1162 | // determine first the nr. of bits in x | |
b2a00c89 | 1163 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 1164 | q = 16; |
b2a00c89 | 1165 | } else { // if x < 2^53 |
200359e8 L |
1166 | tmp1.d = (double) C1; // exact conversion |
1167 | x_nr_bits = | |
1168 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 | 1169 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 1170 | if (q == 0) { |
b2a00c89 L |
1171 | q = nr_digits[x_nr_bits - 1].digits1; |
1172 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
1173 | q++; |
1174 | } | |
1175 | } | |
1176 | ||
b2a00c89 | 1177 | if (exp >= 0) { // -exp <= 0 |
200359e8 L |
1178 | // the argument is an integer already |
1179 | res = x; | |
1180 | BID_RETURN (res); | |
b2a00c89 | 1181 | } else if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q |
200359e8 L |
1182 | // need to shift right -exp digits from the coefficient; the exp will be 0 |
1183 | ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' | |
1184 | // chop off ind digits from the lower part of C1 | |
1185 | // C1 = C1 + 1/2 * 10^x where the result C1 fits in 64 bits | |
1186 | // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate | |
b2a00c89 | 1187 | C1 = C1 + midpoint64[ind - 1]; |
200359e8 L |
1188 | // calculate C* and f* |
1189 | // C* is actually floor(C*) in this case | |
1190 | // C* and f* need shifting and masking, as shown by | |
b2a00c89 | 1191 | // shiftright128[] and maskhigh128[] |
200359e8 | 1192 | // 1 <= x <= 16 |
b2a00c89 | 1193 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
1194 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) |
1195 | // the approximation of 10^(-x) was rounded up to 64 bits | |
b2a00c89 | 1196 | __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); |
200359e8 L |
1197 | |
1198 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
1199 | // C* = floor(C*) - logical right shift; C* has p decimal digits, | |
1200 | // correct by Prop. 1) | |
1201 | // else | |
1202 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
1203 | // correct by Property 1) | |
1204 | // n = C* * 10^(e+x) | |
1205 | ||
b2a00c89 | 1206 | if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 |
200359e8 | 1207 | res = P128.w[1]; |
b2a00c89 L |
1208 | } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 |
1209 | shift = shiftright128[ind - 1]; // 3 <= shift <= 63 | |
200359e8 L |
1210 | res = (P128.w[1] >> shift); |
1211 | } | |
1212 | // midpoints are already rounded correctly | |
1213 | // set exponent to zero as it was negative before. | |
1214 | res = x_sign | 0x31c0000000000000ull | res; | |
1215 | BID_RETURN (res); | |
b2a00c89 | 1216 | } else { // if exp < 0 and q + exp < 0 |
200359e8 L |
1217 | // the result is +0 or -0 |
1218 | res = x_sign | 0x31c0000000000000ull; | |
1219 | BID_RETURN (res); | |
1220 | } | |
1221 | } |