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 #include "bid_internal.h"
25
26 /*****************************************************************************
27 * BID64_round_integral_exact
28 ****************************************************************************/
29
30 #if DECIMAL_CALL_BY_REFERENCE
31 void
32 bid64_round_integral_exact (UINT64 * pres,
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
42 bid64_round_integral_exact (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
43 _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
44 #endif
45
46 UINT64 res = 0xbaddbaddbaddbaddull;
47 UINT64 x_sign;
48 int exp; // unbiased exponent
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
55 UINT128 fstar = { {0x0ull, 0x0ull} };
56 UINT128 P128;
57
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
67 // set invalid flag
68 *pfpsf |= INVALID_EXCEPTION;
69 // return quiet (SNaN)
70 res = x & 0xfdffffffffffffffull;
71 } else { // QNaN
72 res = x;
73 }
74 BID_RETURN (res);
75 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity
76 res = x_sign | 0x7800000000000000ull;
77 BID_RETURN (res);
78 }
79 // unpack x
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;
85 if (C1 > 9999999999999999ull) { // non-canonical
86 C1 = 0;
87 }
88 } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS)
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;
145 } // end switch ()
146
147 // q = nr. of decimal digits in x (1 <= q <= 54)
148 // determine first the nr. of bits in x
149 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
150 q = 16;
151 } else { // if x < 2^53
152 tmp1.d = (double) C1; // exact conversion
153 x_nr_bits =
154 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
155 q = nr_digits[x_nr_bits - 1].digits;
156 if (q == 0) {
157 q = nr_digits[x_nr_bits - 1].digits1;
158 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
159 q++;
160 }
161 }
162
163 if (exp >= 0) { // -exp <= 0
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:
171 if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
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
177 C1 = C1 + midpoint64[ind - 1];
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
181 // shiftright128[] and maskhigh128[]
182 // 1 <= x <= 16
183 // kx = 10^(-x) = ten2mk64[ind - 1]
184 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
185 // the approximation of 10^(-x) was rounded up to 64 bits
186 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
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
198 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
199 res = P128.w[1];
200 fstar.w[1] = 0;
201 fstar.w[0] = P128.w[0];
202 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
203 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
204 res = (P128.w[1] >> shift);
205 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
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)
211 && (fstar.w[0] < ten2mk64[ind - 1])) {
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
223 if ((fstar.w[0] - 0x8000000000000000ull) > ten2mk64[ind - 1]) {
224 // set the inexact flag
225 *pfpsf |= INEXACT_EXCEPTION;
226 } // else the result is exact
227 } else { // the result is inexact; f2* <= 1/2
228 // set the inexact flag
229 *pfpsf |= INEXACT_EXCEPTION;
230 }
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])) {
234 // f2* > 1/2 and the result may be exact
235 // Calculate f2* - 1/2
236 if (fstar.w[1] > onehalf128[ind - 1]
237 || fstar.w[0] > ten2mk64[ind - 1]) {
238 // set the inexact flag
239 *pfpsf |= INEXACT_EXCEPTION;
240 } // else the result is exact
241 } else { // the result is inexact; f2* <= 1/2
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);
249 } else { // if exp < 0 and q + exp < 0
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:
257 if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
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
263 C1 = C1 + midpoint64[ind - 1];
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
267 // shiftright128[] and maskhigh128[]
268 // 1 <= x <= 16
269 // kx = 10^(-x) = ten2mk64[ind - 1]
270 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
271 // the approximation of 10^(-x) was rounded up to 64 bits
272 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
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
282 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
283 res = P128.w[1];
284 fstar.w[1] = 0;
285 fstar.w[0] = P128.w[0];
286 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
287 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
288 res = (P128.w[1] >> shift);
289 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
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
302 if ((fstar.w[0] - 0x8000000000000000ull) > ten2mk64[ind - 1]) {
303 // set the inexact flag
304 *pfpsf |= INEXACT_EXCEPTION;
305 } // else the result is exact
306 } else { // the result is inexact; f2* <= 1/2
307 // set the inexact flag
308 *pfpsf |= INEXACT_EXCEPTION;
309 }
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])) {
313 // f2* > 1/2 and the result may be exact
314 // Calculate f2* - 1/2
315 if (fstar.w[1] > onehalf128[ind - 1]
316 || fstar.w[0] > ten2mk64[ind - 1]) {
317 // set the inexact flag
318 *pfpsf |= INEXACT_EXCEPTION;
319 } // else the result is exact
320 } else { // the result is inexact; f2* <= 1/2
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);
328 } else { // if exp < 0 and q + exp < 0
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:
336 if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
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
344 // shiftright128[] and maskhigh128[]
345 // 1 <= x <= 16
346 // kx = 10^(-x) = ten2mk64[ind - 1]
347 // C* = C1 * 10^(-x)
348 // the approximation of 10^(-x) was rounded up to 64 bits
349 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
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
356 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
357 res = P128.w[1];
358 fstar.w[1] = 0;
359 fstar.w[0] = P128.w[0];
360 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
361 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
362 res = (P128.w[1] >> shift);
363 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
364 fstar.w[0] = P128.w[0];
365 }
366 // if (f* > 10^(-x)) then the result is inexact
367 if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1])) {
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);
377 } else { // if exp < 0 and q + exp <= 0
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:
389 if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
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
397 // shiftright128[] and maskhigh128[]
398 // 1 <= x <= 16
399 // kx = 10^(-x) = ten2mk64[ind - 1]
400 // C* = C1 * 10^(-x)
401 // the approximation of 10^(-x) was rounded up to 64 bits
402 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
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
409 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
410 res = P128.w[1];
411 fstar.w[1] = 0;
412 fstar.w[0] = P128.w[0];
413 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
414 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
415 res = (P128.w[1] >> shift);
416 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
417 fstar.w[0] = P128.w[0];
418 }
419 // if (f* > 10^(-x)) then the result is inexact
420 if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1])) {
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);
430 } else { // if exp < 0 and q + exp <= 0
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:
442 if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
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
450 // shiftright128[] and maskhigh128[]
451 // 1 <= x <= 16
452 // kx = 10^(-x) = ten2mk64[ind - 1]
453 // C* = C1 * 10^(-x)
454 // the approximation of 10^(-x) was rounded up to 64 bits
455 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
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
462 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
463 res = P128.w[1];
464 fstar.w[1] = 0;
465 fstar.w[0] = P128.w[0];
466 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
467 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
468 res = (P128.w[1] >> shift);
469 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
470 fstar.w[0] = P128.w[0];
471 }
472 // if (f* > 10^(-x)) then the result is inexact
473 if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1])) {
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);
479 } else { // if exp < 0 and q + exp < 0
480 // the result is +0 or -0
481 res = x_sign | 0x31c0000000000000ull;
482 *pfpsf |= INEXACT_EXCEPTION;
483 BID_RETURN (res);
484 }
485 break;
486 } // end switch ()
487 BID_RETURN (res);
488 }
489
490 /*****************************************************************************
491 * BID64_round_integral_nearest_even
492 ****************************************************************************/
493
494 #if DECIMAL_CALL_BY_REFERENCE
495 void
496 bid64_round_integral_nearest_even (UINT64 * pres,
497 UINT64 *
498 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
499 _EXC_INFO_PARAM) {
500 UINT64 x = *px;
501 #else
502 UINT64
503 bid64_round_integral_nearest_even (UINT64 x _EXC_FLAGS_PARAM
504 _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
505 #endif
506
507 UINT64 res = 0xbaddbaddbaddbaddull;
508 UINT64 x_sign;
509 int exp; // unbiased exponent
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
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
528 *pfpsf |= INVALID_EXCEPTION;
529 // return quiet (SNaN)
530 res = x & 0xfdffffffffffffffull;
531 } else { // QNaN
532 res = x;
533 }
534 BID_RETURN (res);
535 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity
536 res = x_sign | 0x7800000000000000ull;
537 BID_RETURN (res);
538 }
539 // unpack x
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;
545 if (C1 > 9999999999999999ull) { // non-canonical
546 C1 = 0;
547 }
548 } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS)
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
569 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
570 q = 16;
571 } else { // if x < 2^53
572 tmp1.d = (double) C1; // exact conversion
573 x_nr_bits =
574 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
575 q = nr_digits[x_nr_bits - 1].digits;
576 if (q == 0) {
577 q = nr_digits[x_nr_bits - 1].digits1;
578 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
579 q++;
580 }
581 }
582
583 if (exp >= 0) { // -exp <= 0
584 // the argument is an integer already
585 res = x;
586 BID_RETURN (res);
587 } else if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
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
593 C1 = C1 + midpoint64[ind - 1];
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
597 // shiftright128[] and maskhigh128[]
598 // 1 <= x <= 16
599 // kx = 10^(-x) = ten2mk64[ind - 1]
600 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
601 // the approximation of 10^(-x) was rounded up to 64 bits
602 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
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
614 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
615 res = P128.w[1];
616 fstar.w[1] = 0;
617 fstar.w[0] = P128.w[0];
618 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
619 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
620 res = (P128.w[1] >> shift);
621 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
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)
627 && (fstar.w[0] < ten2mk64[ind - 1])) {
628 res--;
629 }
630 // set exponent to zero as it was negative before.
631 res = x_sign | 0x31c0000000000000ull | res;
632 BID_RETURN (res);
633 } else { // if exp < 0 and q + exp < 0
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
646 bid64_round_integral_negative (UINT64 * pres,
647 UINT64 *
648 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
649 _EXC_INFO_PARAM) {
650 UINT64 x = *px;
651 #else
652 UINT64
653 bid64_round_integral_negative (UINT64 x _EXC_FLAGS_PARAM
654 _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
655 #endif
656
657 UINT64 res = 0xbaddbaddbaddbaddull;
658 UINT64 x_sign;
659 int exp; // unbiased exponent
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
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
679 *pfpsf |= INVALID_EXCEPTION;
680 // return quiet (SNaN)
681 res = x & 0xfdffffffffffffffull;
682 } else { // QNaN
683 res = x;
684 }
685 BID_RETURN (res);
686 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity
687 res = x_sign | 0x7800000000000000ull;
688 BID_RETURN (res);
689 }
690 // unpack x
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;
696 if (C1 > 9999999999999999ull) { // non-canonical
697 C1 = 0;
698 }
699 } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS)
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
724 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
725 q = 16;
726 } else { // if x < 2^53
727 tmp1.d = (double) C1; // exact conversion
728 x_nr_bits =
729 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
730 q = nr_digits[x_nr_bits - 1].digits;
731 if (q == 0) {
732 q = nr_digits[x_nr_bits - 1].digits1;
733 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
734 q++;
735 }
736 }
737
738 if (exp >= 0) { // -exp <= 0
739 // the argument is an integer already
740 res = x;
741 BID_RETURN (res);
742 } else if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
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
750 // shiftright128[] and maskhigh128[]
751 // 1 <= x <= 16
752 // kx = 10^(-x) = ten2mk64[ind - 1]
753 // C* = C1 * 10^(-x)
754 // the approximation of 10^(-x) was rounded up to 64 bits
755 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
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
762 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
763 res = P128.w[1];
764 fstar.w[1] = 0;
765 fstar.w[0] = P128.w[0];
766 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
767 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
768 res = (P128.w[1] >> shift);
769 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
770 fstar.w[0] = P128.w[0];
771 }
772 // if (f* > 10^(-x)) then the result is inexact
773 if (x_sign
774 && ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1]))) {
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);
781 } else { // if exp < 0 and q + exp <= 0
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
798 bid64_round_integral_positive (UINT64 * pres,
799 UINT64 *
800 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
801 _EXC_INFO_PARAM) {
802 UINT64 x = *px;
803 #else
804 UINT64
805 bid64_round_integral_positive (UINT64 x _EXC_FLAGS_PARAM
806 _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
807 #endif
808
809 UINT64 res = 0xbaddbaddbaddbaddull;
810 UINT64 x_sign;
811 int exp; // unbiased exponent
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
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
831 *pfpsf |= INVALID_EXCEPTION;
832 // return quiet (SNaN)
833 res = x & 0xfdffffffffffffffull;
834 } else { // QNaN
835 res = x;
836 }
837 BID_RETURN (res);
838 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity
839 res = x_sign | 0x7800000000000000ull;
840 BID_RETURN (res);
841 }
842 // unpack x
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;
848 if (C1 > 9999999999999999ull) { // non-canonical
849 C1 = 0;
850 }
851 } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS)
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
876 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
877 q = 16;
878 } else { // if x < 2^53
879 tmp1.d = (double) C1; // exact conversion
880 x_nr_bits =
881 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
882 q = nr_digits[x_nr_bits - 1].digits;
883 if (q == 0) {
884 q = nr_digits[x_nr_bits - 1].digits1;
885 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
886 q++;
887 }
888 }
889
890 if (exp >= 0) { // -exp <= 0
891 // the argument is an integer already
892 res = x;
893 BID_RETURN (res);
894 } else if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
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
902 // shiftright128[] and maskhigh128[]
903 // 1 <= x <= 16
904 // kx = 10^(-x) = ten2mk64[ind - 1]
905 // C* = C1 * 10^(-x)
906 // the approximation of 10^(-x) was rounded up to 64 bits
907 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
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
914 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
915 res = P128.w[1];
916 fstar.w[1] = 0;
917 fstar.w[0] = P128.w[0];
918 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
919 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
920 res = (P128.w[1] >> shift);
921 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
922 fstar.w[0] = P128.w[0];
923 }
924 // if (f* > 10^(-x)) then the result is inexact
925 if (!x_sign
926 && ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1]))) {
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);
933 } else { // if exp < 0 and q + exp <= 0
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
950 bid64_round_integral_zero (UINT64 * pres,
951 UINT64 *
952 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
953 _EXC_INFO_PARAM) {
954 UINT64 x = *px;
955 #else
956 UINT64
957 bid64_round_integral_zero (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
958 _EXC_INFO_PARAM) {
959 #endif
960
961 UINT64 res = 0xbaddbaddbaddbaddull;
962 UINT64 x_sign;
963 int exp; // unbiased exponent
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
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
982 *pfpsf |= INVALID_EXCEPTION;
983 // return quiet (SNaN)
984 res = x & 0xfdffffffffffffffull;
985 } else { // QNaN
986 res = x;
987 }
988 BID_RETURN (res);
989 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity
990 res = x_sign | 0x7800000000000000ull;
991 BID_RETURN (res);
992 }
993 // unpack x
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;
999 if (C1 > 9999999999999999ull) { // non-canonical
1000 C1 = 0;
1001 }
1002 } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS)
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
1023 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
1024 q = 16;
1025 } else { // if x < 2^53
1026 tmp1.d = (double) C1; // exact conversion
1027 x_nr_bits =
1028 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1029 q = nr_digits[x_nr_bits - 1].digits;
1030 if (q == 0) {
1031 q = nr_digits[x_nr_bits - 1].digits1;
1032 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
1033 q++;
1034 }
1035 }
1036
1037 if (exp >= 0) { // -exp <= 0
1038 // the argument is an integer already
1039 res = x;
1040 BID_RETURN (res);
1041 } else if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
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
1049 // shiftright128[] and maskhigh128[]
1050 // 1 <= x <= 16
1051 // kx = 10^(-x) = ten2mk64[ind - 1]
1052 // C* = C1 * 10^(-x)
1053 // the approximation of 10^(-x) was rounded up to 64 bits
1054 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
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
1061 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
1062 res = P128.w[1];
1063 // redundant fstar.w[1] = 0;
1064 // redundant fstar.w[0] = P128.w[0];
1065 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1066 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
1067 res = (P128.w[1] >> shift);
1068 // redundant fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
1069 // redundant fstar.w[0] = P128.w[0];
1070 }
1071 // if (f* > 10^(-x)) then the result is inexact
1072 // if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind-1])){
1073 // // redundant
1074 // }
1075 // set exponent to zero as it was negative before.
1076 res = x_sign | 0x31c0000000000000ull | res;
1077 BID_RETURN (res);
1078 } else { // if exp < 0 and q + exp < 0
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
1091 bid64_round_integral_nearest_away (UINT64 * pres,
1092 UINT64 *
1093 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1094 _EXC_INFO_PARAM) {
1095 UINT64 x = *px;
1096 #else
1097 UINT64
1098 bid64_round_integral_nearest_away (UINT64 x _EXC_FLAGS_PARAM
1099 _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
1100 #endif
1101
1102 UINT64 res = 0xbaddbaddbaddbaddull;
1103 UINT64 x_sign;
1104 int exp; // unbiased exponent
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
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
1122 *pfpsf |= INVALID_EXCEPTION;
1123 // return quiet (SNaN)
1124 res = x & 0xfdffffffffffffffull;
1125 } else { // QNaN
1126 res = x;
1127 }
1128 BID_RETURN (res);
1129 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity
1130 res = x_sign | 0x7800000000000000ull;
1131 BID_RETURN (res);
1132 }
1133 // unpack x
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;
1139 if (C1 > 9999999999999999ull) { // non-canonical
1140 C1 = 0;
1141 }
1142 } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS)
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
1163 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
1164 q = 16;
1165 } else { // if x < 2^53
1166 tmp1.d = (double) C1; // exact conversion
1167 x_nr_bits =
1168 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1169 q = nr_digits[x_nr_bits - 1].digits;
1170 if (q == 0) {
1171 q = nr_digits[x_nr_bits - 1].digits1;
1172 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
1173 q++;
1174 }
1175 }
1176
1177 if (exp >= 0) { // -exp <= 0
1178 // the argument is an integer already
1179 res = x;
1180 BID_RETURN (res);
1181 } else if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
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
1187 C1 = C1 + midpoint64[ind - 1];
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
1191 // shiftright128[] and maskhigh128[]
1192 // 1 <= x <= 16
1193 // kx = 10^(-x) = ten2mk64[ind - 1]
1194 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1195 // the approximation of 10^(-x) was rounded up to 64 bits
1196 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]);
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
1206 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
1207 res = P128.w[1];
1208 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1209 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
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);
1216 } else { // if exp < 0 and q + exp < 0
1217 // the result is +0 or -0
1218 res = x_sign | 0x31c0000000000000ull;
1219 BID_RETURN (res);
1220 }
1221 }