1 /* Copyright (C) 2007-2020 Free Software Foundation, Inc.
3 This file is part of GCC.
5 GCC is free software; you can redistribute it and/or modify it under
6 the terms of the GNU General Public License as published by the Free
7 Software Foundation; either version 3, or (at your option) any later
10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11 WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 Under Section 7 of GPL version 3, you are granted additional
16 permissions described in the GCC Runtime Library Exception, version
17 3.1, as published by the Free Software Foundation.
19 You should have received a copy of the GNU General Public License and
20 a copy of the GCC Runtime Library Exception along with this program;
21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
22 <http://www.gnu.org/licenses/>. */
24 #include "bid_internal.h"
26 /*****************************************************************************
27 * BID64_to_uint32_rnint
28 ****************************************************************************/
30 #if DECIMAL_CALL_BY_REFERENCE
32 bid64_to_uint32_rnint (unsigned int *pres
, UINT64
* px
33 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
38 bid64_to_uint32_rnint (UINT64 x
39 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
45 int exp
; // unbiased exponent
46 // Note: C1 represents x_significand (UINT64)
49 unsigned int x_nr_bits
;
52 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
56 // check for NaN or Infinity
57 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
59 *pfpsf
|= INVALID_EXCEPTION
;
60 // return Integer Indefinite
65 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
66 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
67 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
68 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
69 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
70 if (C1
> 9999999999999999ull) { // non-canonical
75 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
76 C1
= x
& MASK_BINARY_SIG1
;
79 // check for zeros (possibly from non-canonical values)
85 // x is not special and is not zero
87 // q = nr. of decimal digits in x (1 <= q <= 54)
88 // determine first the nr. of bits in x
89 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
90 // split the 64-bit value in two 32-bit halves to avoid rounding errors
91 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
92 tmp1
.d
= (double) (C1
>> 32); // exact conversion
94 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
96 tmp1
.d
= (double) C1
; // exact conversion
98 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
100 } else { // if x < 2^53
101 tmp1
.d
= (double) C1
; // exact conversion
103 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
105 q
= nr_digits
[x_nr_bits
- 1].digits
;
107 q
= nr_digits
[x_nr_bits
- 1].digits1
;
108 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
111 exp
= x_exp
- 398; // unbiased exponent
113 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
115 *pfpsf
|= INVALID_EXCEPTION
;
116 // return Integer Indefinite
119 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
120 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
121 // so x rounded to an integer may or may not fit in an unsigned 32-bit int
122 // the cases that do not fit are identified here; the ones that fit
123 // fall through and will be handled with other cases further,
124 // under '1 <= q + exp <= 10'
125 if (x_sign
) { // if n < 0 and q + exp = 10 then x is much less than -1/2
126 // => set invalid flag
127 *pfpsf
|= INVALID_EXCEPTION
;
128 // return Integer Indefinite
131 } else { // if n > 0 and q + exp = 10
132 // if n >= 2^32 - 1/2 then n is too large
133 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
134 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=16
135 // <=> C * 10^(11-q) >= 0x9fffffffb, 1<=q<=16
137 // Note: C * 10^(11-q) has 10 or 11 digits; 0x9fffffffb has 11 digits
138 tmp64
= C1
* ten2k64
[11 - q
]; // C scaled up to 11-digit int
139 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
140 if (tmp64
>= 0x9fffffffbull
) {
142 *pfpsf
|= INVALID_EXCEPTION
;
143 // return Integer Indefinite
147 // else cases that can be rounded to a 32-bit unsigned int fall through
148 // to '1 <= q + exp <= 10'
149 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
150 // C * 10^(11-q) >= 0x9fffffffb <=>
151 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 5
152 // (scale 2^32-1/2 up)
153 // Note: 0x9fffffffb*10^(q-11) has q-1 or q digits, where q <= 16
154 tmp64
= 0x9fffffffbull
* ten2k64
[q
- 11];
157 *pfpsf
|= INVALID_EXCEPTION
;
158 // return Integer Indefinite
162 // else cases that can be rounded to a 32-bit int fall through
163 // to '1 <= q + exp <= 10'
167 // n is not too large to be converted to int32 if -1/2 <= n < 2^32 - 1/2
168 // Note: some of the cases tested for above fall through to this point
169 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
173 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
174 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
181 if (C1
<= midpoint64
[ind
]) {
182 res
= 0x00000000; // return 0
183 } else if (x_sign
) { // n < 0
185 *pfpsf
|= INVALID_EXCEPTION
;
186 // return Integer Indefinite
190 res
= 0x00000001; // return +1
192 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
193 // -2^32-1/2 <= x <= -1 or 1 <= x < 2^32-1/2 so if positive, x can be
194 // rounded to nearest to a 32-bit unsigned integer
195 if (x_sign
) { // x <= -1
197 *pfpsf
|= INVALID_EXCEPTION
;
198 // return Integer Indefinite
202 // 1 <= x < 2^32-1/2 so x can be rounded
203 // to nearest to a 32-bit unsigned integer
204 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
205 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
206 // chop off ind digits from the lower part of C1
207 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
208 C1
= C1
+ midpoint64
[ind
- 1];
209 // calculate C* and f*
210 // C* is actually floor(C*) in this case
211 // C* and f* need shifting and masking, as shown by
212 // shiftright128[] and maskhigh128[]
214 // kx = 10^(-x) = ten2mk64[ind - 1]
215 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
216 // the approximation of 10^(-x) was rounded up to 54 bits
217 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
219 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
220 fstar
.w
[0] = P128
.w
[0];
221 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
222 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
223 // if (0 < f* < 10^(-x)) then the result is a midpoint
224 // if floor(C*) is even then C* = floor(C*) - logical right
225 // shift; C* has p decimal digits, correct by Prop. 1)
226 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
227 // shift; C* has p decimal digits, correct by Pr. 1)
229 // C* = floor(C*) (logical right shift; C has p decimal digits,
230 // correct by Property 1)
233 // shift right C* by Ex-64 = shiftright128[ind]
234 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
235 Cstar
= Cstar
>> shift
;
237 // if the result was a midpoint it was rounded away from zero, so
238 // it will need a correction
239 // check for midpoints
240 if ((fstar
.w
[1] == 0) && fstar
.w
[0] &&
241 (fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[1])) {
242 // ten2mk128trunc[ind -1].w[1] is identical to
243 // ten2mk128[ind -1].w[1]
244 // the result is a midpoint; round to nearest
245 if (Cstar
& 0x01) { // Cstar is odd; MP in [EVEN, ODD]
246 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
247 Cstar
--; // Cstar is now even
248 } // else MP in [ODD, EVEN]
250 res
= Cstar
; // the result is positive
251 } else if (exp
== 0) {
254 res
= C1
; // the result is positive
255 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
256 // res = +C * 10^exp (exact)
257 res
= C1
* ten2k64
[exp
]; // the result is positive
263 /*****************************************************************************
264 * BID64_to_uint32_xrnint
265 ****************************************************************************/
267 #if DECIMAL_CALL_BY_REFERENCE
269 bid64_to_uint32_xrnint (unsigned int *pres
, UINT64
* px
270 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
275 bid64_to_uint32_xrnint (UINT64 x
276 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
282 int exp
; // unbiased exponent
283 // Note: C1 represents x_significand (UINT64)
286 unsigned int x_nr_bits
;
289 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
293 // check for NaN or Infinity
294 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
296 *pfpsf
|= INVALID_EXCEPTION
;
297 // return Integer Indefinite
302 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
303 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
304 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
305 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
306 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
307 if (C1
> 9999999999999999ull) { // non-canonical
312 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
313 C1
= x
& MASK_BINARY_SIG1
;
316 // check for zeros (possibly from non-canonical values)
322 // x is not special and is not zero
324 // q = nr. of decimal digits in x (1 <= q <= 54)
325 // determine first the nr. of bits in x
326 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
327 // split the 64-bit value in two 32-bit halves to avoid rounding errors
328 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
329 tmp1
.d
= (double) (C1
>> 32); // exact conversion
331 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
333 tmp1
.d
= (double) C1
; // exact conversion
335 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
337 } else { // if x < 2^53
338 tmp1
.d
= (double) C1
; // exact conversion
340 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
342 q
= nr_digits
[x_nr_bits
- 1].digits
;
344 q
= nr_digits
[x_nr_bits
- 1].digits1
;
345 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
348 exp
= x_exp
- 398; // unbiased exponent
350 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
352 *pfpsf
|= INVALID_EXCEPTION
;
353 // return Integer Indefinite
356 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
357 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
358 // so x rounded to an integer may or may not fit in an unsigned 32-bit int
359 // the cases that do not fit are identified here; the ones that fit
360 // fall through and will be handled with other cases further,
361 // under '1 <= q + exp <= 10'
362 if (x_sign
) { // if n < 0 and q + exp = 10 then x is much less than -1/2
363 // => set invalid flag
364 *pfpsf
|= INVALID_EXCEPTION
;
365 // return Integer Indefinite
368 } else { // if n > 0 and q + exp = 10
369 // if n >= 2^32 - 1/2 then n is too large
370 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
371 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=16
372 // <=> C * 10^(11-q) >= 0x9fffffffb, 1<=q<=16
374 // Note: C * 10^(11-q) has 10 or 11 digits; 0x9fffffffb has 11 digits
375 tmp64
= C1
* ten2k64
[11 - q
]; // C scaled up to 11-digit int
376 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
377 if (tmp64
>= 0x9fffffffbull
) {
379 *pfpsf
|= INVALID_EXCEPTION
;
380 // return Integer Indefinite
384 // else cases that can be rounded to a 32-bit unsigned int fall through
385 // to '1 <= q + exp <= 10'
386 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
387 // C * 10^(11-q) >= 0x9fffffffb <=>
388 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 5
389 // (scale 2^32-1/2 up)
390 // Note: 0x9fffffffb*10^(q-11) has q-1 or q digits, where q <= 16
391 tmp64
= 0x9fffffffbull
* ten2k64
[q
- 11];
394 *pfpsf
|= INVALID_EXCEPTION
;
395 // return Integer Indefinite
399 // else cases that can be rounded to a 32-bit int fall through
400 // to '1 <= q + exp <= 10'
404 // n is not too large to be converted to int32 if -1/2 <= n < 2^32 - 1/2
405 // Note: some of the cases tested for above fall through to this point
406 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
408 *pfpsf
|= INEXACT_EXCEPTION
;
412 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
413 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
420 if (C1
<= midpoint64
[ind
]) {
421 res
= 0x00000000; // return 0
422 } else if (x_sign
) { // n < 0
424 *pfpsf
|= INVALID_EXCEPTION
;
425 // return Integer Indefinite
429 res
= 0x00000001; // return +1
432 *pfpsf
|= INEXACT_EXCEPTION
;
433 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
434 // -2^32-1/2 <= x <= -1 or 1 <= x < 2^32-1/2 so if positive, x can be
435 // rounded to nearest to a 32-bit unsigned integer
436 if (x_sign
) { // x <= -1
438 *pfpsf
|= INVALID_EXCEPTION
;
439 // return Integer Indefinite
443 // 1 <= x < 2^32-1/2 so x can be rounded
444 // to nearest to a 32-bit unsigned integer
445 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
446 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
447 // chop off ind digits from the lower part of C1
448 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
449 C1
= C1
+ midpoint64
[ind
- 1];
450 // calculate C* and f*
451 // C* is actually floor(C*) in this case
452 // C* and f* need shifting and masking, as shown by
453 // shiftright128[] and maskhigh128[]
455 // kx = 10^(-x) = ten2mk64[ind - 1]
456 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
457 // the approximation of 10^(-x) was rounded up to 54 bits
458 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
460 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
461 fstar
.w
[0] = P128
.w
[0];
462 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
463 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
464 // if (0 < f* < 10^(-x)) then the result is a midpoint
465 // if floor(C*) is even then C* = floor(C*) - logical right
466 // shift; C* has p decimal digits, correct by Prop. 1)
467 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
468 // shift; C* has p decimal digits, correct by Pr. 1)
470 // C* = floor(C*) (logical right shift; C has p decimal digits,
471 // correct by Property 1)
474 // shift right C* by Ex-64 = shiftright128[ind]
475 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
476 Cstar
= Cstar
>> shift
;
477 // determine inexactness of the rounding of C*
478 // if (0 < f* - 1/2 < 10^(-x)) then
479 // the result is exact
480 // else // if (f* - 1/2 > T*) then
481 // the result is inexact
482 if (ind
- 1 <= 2) { // fstar.w[1] is 0
483 if (fstar
.w
[0] > 0x8000000000000000ull
) {
484 // f* > 1/2 and the result may be exact
485 tmp64
= fstar
.w
[0] - 0x8000000000000000ull
; // f* - 1/2
486 if ((tmp64
> ten2mk128trunc
[ind
- 1].w
[1])) {
487 // ten2mk128trunc[ind -1].w[1] is identical to
488 // ten2mk128[ind -1].w[1]
489 // set the inexact flag
490 *pfpsf
|= INEXACT_EXCEPTION
;
491 } // else the result is exact
492 } else { // the result is inexact; f2* <= 1/2
493 // set the inexact flag
494 *pfpsf
|= INEXACT_EXCEPTION
;
496 } else { // if 3 <= ind - 1 <= 14
497 if (fstar
.w
[1] > onehalf128
[ind
- 1] ||
498 (fstar
.w
[1] == onehalf128
[ind
- 1] && fstar
.w
[0])) {
499 // f2* > 1/2 and the result may be exact
500 // Calculate f2* - 1/2
501 tmp64
= fstar
.w
[1] - onehalf128
[ind
- 1];
502 if (tmp64
|| fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
503 // ten2mk128trunc[ind -1].w[1] is identical to
504 // ten2mk128[ind -1].w[1]
505 // set the inexact flag
506 *pfpsf
|= INEXACT_EXCEPTION
;
507 } // else the result is exact
508 } else { // the result is inexact; f2* <= 1/2
509 // set the inexact flag
510 *pfpsf
|= INEXACT_EXCEPTION
;
514 // if the result was a midpoint it was rounded away from zero, so
515 // it will need a correction
516 // check for midpoints
517 if ((fstar
.w
[1] == 0) && fstar
.w
[0] &&
518 (fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[1])) {
519 // ten2mk128trunc[ind -1].w[1] is identical to
520 // ten2mk128[ind -1].w[1]
521 // the result is a midpoint; round to nearest
522 if (Cstar
& 0x01) { // Cstar is odd; MP in [EVEN, ODD]
523 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
524 Cstar
--; // Cstar is now even
525 } // else MP in [ODD, EVEN]
527 res
= Cstar
; // the result is positive
528 } else if (exp
== 0) {
531 res
= C1
; // the result is positive
532 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
533 // res = +C * 10^exp (exact)
534 res
= C1
* ten2k64
[exp
]; // the result is positive
540 /*****************************************************************************
541 * BID64_to_uint32_floor
542 ****************************************************************************/
544 #if DECIMAL_CALL_BY_REFERENCE
546 bid64_to_uint32_floor (unsigned int *pres
, UINT64
* px
547 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
552 bid64_to_uint32_floor (UINT64 x
553 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
559 int exp
; // unbiased exponent
560 // Note: C1 represents x_significand (UINT64)
563 unsigned int x_nr_bits
;
566 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
569 // check for NaN or Infinity
570 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
572 *pfpsf
|= INVALID_EXCEPTION
;
573 // return Integer Indefinite
578 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
579 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
580 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
581 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
582 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
583 if (C1
> 9999999999999999ull) { // non-canonical
588 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
589 C1
= x
& MASK_BINARY_SIG1
;
592 // check for zeros (possibly from non-canonical values)
598 // x is not special and is not zero
600 if (x_sign
) { // if n < 0 the conversion is invalid
602 *pfpsf
|= INVALID_EXCEPTION
;
603 // return Integer Indefinite
607 // q = nr. of decimal digits in x (1 <= q <= 54)
608 // determine first the nr. of bits in x
609 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
610 // split the 64-bit value in two 32-bit halves to avoid rounding errors
611 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
612 tmp1
.d
= (double) (C1
>> 32); // exact conversion
614 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
616 tmp1
.d
= (double) C1
; // exact conversion
618 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
620 } else { // if x < 2^53
621 tmp1
.d
= (double) C1
; // exact conversion
623 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
625 q
= nr_digits
[x_nr_bits
- 1].digits
;
627 q
= nr_digits
[x_nr_bits
- 1].digits1
;
628 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
631 exp
= x_exp
- 398; // unbiased exponent
633 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
635 *pfpsf
|= INVALID_EXCEPTION
;
636 // return Integer Indefinite
639 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
640 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
641 // so x rounded to an integer may or may not fit in an unsigned 32-bit int
642 // the cases that do not fit are identified here; the ones that fit
643 // fall through and will be handled with other cases further,
644 // under '1 <= q + exp <= 10'
645 // n > 0 and q + exp = 10
646 // if n >= 2^32 then n is too large
647 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
648 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=16
649 // <=> C * 10^(11-q) >= 0xa00000000, 1<=q<=16
651 // Note: C * 10^(11-q) has 10 or 11 digits; 0xa00000000 has 11 digits
652 tmp64
= C1
* ten2k64
[11 - q
]; // C scaled up to 11-digit int
653 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
654 if (tmp64
>= 0xa00000000ull
) {
656 *pfpsf
|= INVALID_EXCEPTION
;
657 // return Integer Indefinite
661 // else cases that can be rounded to a 32-bit unsigned int fall through
662 // to '1 <= q + exp <= 10'
663 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
664 // C * 10^(11-q) >= 0xa00000000 <=>
665 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 5
666 // (scale 2^32-1/2 up)
667 // Note: 0xa00000000*10^(q-11) has q-1 or q digits, where q <= 16
668 tmp64
= 0xa00000000ull
* ten2k64
[q
- 11];
671 *pfpsf
|= INVALID_EXCEPTION
;
672 // return Integer Indefinite
676 // else cases that can be rounded to a 32-bit int fall through
677 // to '1 <= q + exp <= 10'
680 // n is not too large to be converted to int32 if -1 < n < 2^32
681 // Note: some of the cases tested for above fall through to this point
682 if ((q
+ exp
) <= 0) { // n = +0.[0...0]c(0)c(1)...c(q-1)
686 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
687 // 1 <= x < 2^32 so x can be rounded
688 // to nearest to a 32-bit unsigned integer
689 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
690 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
691 // chop off ind digits from the lower part of C1
692 // C1 fits in 64 bits
693 // calculate C* and f*
694 // C* is actually floor(C*) in this case
695 // C* and f* need shifting and masking, as shown by
696 // shiftright128[] and maskhigh128[]
698 // kx = 10^(-x) = ten2mk64[ind - 1]
700 // the approximation of 10^(-x) was rounded up to 54 bits
701 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
703 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
704 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
705 // C* = floor(C*) (logical right shift; C has p decimal digits,
706 // correct by Property 1)
709 // shift right C* by Ex-64 = shiftright128[ind]
710 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
711 Cstar
= Cstar
>> shift
;
713 res
= Cstar
; // the result is positive
714 } else if (exp
== 0) {
717 res
= C1
; // the result is positive
718 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
719 // res = +C * 10^exp (exact)
720 res
= C1
* ten2k64
[exp
]; // the result is positive
726 /*****************************************************************************
727 * BID64_to_uint32_xfloor
728 ****************************************************************************/
730 #if DECIMAL_CALL_BY_REFERENCE
732 bid64_to_uint32_xfloor (unsigned int *pres
, UINT64
* px
733 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
738 bid64_to_uint32_xfloor (UINT64 x
739 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
745 int exp
; // unbiased exponent
746 // Note: C1 represents x_significand (UINT64)
749 unsigned int x_nr_bits
;
752 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
756 // check for NaN or Infinity
757 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
759 *pfpsf
|= INVALID_EXCEPTION
;
760 // return Integer Indefinite
765 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
766 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
767 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
768 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
769 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
770 if (C1
> 9999999999999999ull) { // non-canonical
775 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
776 C1
= x
& MASK_BINARY_SIG1
;
779 // check for zeros (possibly from non-canonical values)
785 // x is not special and is not zero
787 if (x_sign
) { // if n < 0 the conversion is invalid
789 *pfpsf
|= INVALID_EXCEPTION
;
790 // return Integer Indefinite
794 // q = nr. of decimal digits in x (1 <= q <= 54)
795 // determine first the nr. of bits in x
796 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
797 // split the 64-bit value in two 32-bit halves to avoid rounding errors
798 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
799 tmp1
.d
= (double) (C1
>> 32); // exact conversion
801 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
803 tmp1
.d
= (double) C1
; // exact conversion
805 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
807 } else { // if x < 2^53
808 tmp1
.d
= (double) C1
; // exact conversion
810 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
812 q
= nr_digits
[x_nr_bits
- 1].digits
;
814 q
= nr_digits
[x_nr_bits
- 1].digits1
;
815 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
818 exp
= x_exp
- 398; // unbiased exponent
820 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
822 *pfpsf
|= INVALID_EXCEPTION
;
823 // return Integer Indefinite
826 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
827 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
828 // so x rounded to an integer may or may not fit in an unsigned 32-bit int
829 // the cases that do not fit are identified here; the ones that fit
830 // fall through and will be handled with other cases further,
831 // under '1 <= q + exp <= 10'
832 // if n > 0 and q + exp = 10
833 // if n >= 2^32 then n is too large
834 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
835 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=16
836 // <=> C * 10^(11-q) >= 0xa00000000, 1<=q<=16
838 // Note: C * 10^(11-q) has 10 or 11 digits; 0xa00000000 has 11 digits
839 tmp64
= C1
* ten2k64
[11 - q
]; // C scaled up to 11-digit int
840 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
841 if (tmp64
>= 0xa00000000ull
) {
843 *pfpsf
|= INVALID_EXCEPTION
;
844 // return Integer Indefinite
848 // else cases that can be rounded to a 32-bit unsigned int fall through
849 // to '1 <= q + exp <= 10'
850 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
851 // C * 10^(11-q) >= 0xa00000000 <=>
852 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 5
853 // (scale 2^32-1/2 up)
854 // Note: 0xa00000000*10^(q-11) has q-1 or q digits, where q <= 16
855 tmp64
= 0xa00000000ull
* ten2k64
[q
- 11];
858 *pfpsf
|= INVALID_EXCEPTION
;
859 // return Integer Indefinite
863 // else cases that can be rounded to a 32-bit int fall through
864 // to '1 <= q + exp <= 10'
867 // n is not too large to be converted to int32 if -1 < n < 2^32
868 // Note: some of the cases tested for above fall through to this point
869 if ((q
+ exp
) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
871 *pfpsf
|= INEXACT_EXCEPTION
;
875 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
876 // 1 <= x < 2^32 so x can be rounded
877 // to nearest to a 32-bit unsigned integer
878 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
879 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
880 // chop off ind digits from the lower part of C1
881 // C1 fits in 64 bits
882 // calculate C* and f*
883 // C* is actually floor(C*) in this case
884 // C* and f* need shifting and masking, as shown by
885 // shiftright128[] and maskhigh128[]
887 // kx = 10^(-x) = ten2mk64[ind - 1]
889 // the approximation of 10^(-x) was rounded up to 54 bits
890 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
892 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
893 fstar
.w
[0] = P128
.w
[0];
894 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
895 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
896 // C* = floor(C*) (logical right shift; C has p decimal digits,
897 // correct by Property 1)
900 // shift right C* by Ex-64 = shiftright128[ind]
901 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
902 Cstar
= Cstar
>> shift
;
903 // determine inexactness of the rounding of C*
904 // if (0 < f* < 10^(-x)) then
905 // the result is exact
906 // else // if (f* > T*) then
907 // the result is inexact
909 if (fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
910 // ten2mk128trunc[ind -1].w[1] is identical to
911 // ten2mk128[ind -1].w[1]
912 // set the inexact flag
913 *pfpsf
|= INEXACT_EXCEPTION
;
914 } // else the result is exact
915 } else { // if 3 <= ind - 1 <= 14
916 if (fstar
.w
[1] || fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
917 // ten2mk128trunc[ind -1].w[1] is identical to
918 // ten2mk128[ind -1].w[1]
919 // set the inexact flag
920 *pfpsf
|= INEXACT_EXCEPTION
;
921 } // else the result is exact
924 res
= Cstar
; // the result is positive
925 } else if (exp
== 0) {
928 res
= C1
; // the result is positive
929 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
930 // res = +C * 10^exp (exact)
931 res
= C1
* ten2k64
[exp
]; // the result is positive
937 /*****************************************************************************
938 * BID64_to_uint32_ceil
939 ****************************************************************************/
941 #if DECIMAL_CALL_BY_REFERENCE
943 bid64_to_uint32_ceil (unsigned int *pres
, UINT64
* px
944 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
949 bid64_to_uint32_ceil (UINT64 x
950 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
956 int exp
; // unbiased exponent
957 // Note: C1 represents x_significand (UINT64)
960 unsigned int x_nr_bits
;
963 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
967 // check for NaN or Infinity
968 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
970 *pfpsf
|= INVALID_EXCEPTION
;
971 // return Integer Indefinite
976 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
977 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
978 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
979 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
980 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
981 if (C1
> 9999999999999999ull) { // non-canonical
986 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
987 C1
= x
& MASK_BINARY_SIG1
;
990 // check for zeros (possibly from non-canonical values)
996 // x is not special and is not zero
998 // q = nr. of decimal digits in x (1 <= q <= 54)
999 // determine first the nr. of bits in x
1000 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
1001 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1002 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
1003 tmp1
.d
= (double) (C1
>> 32); // exact conversion
1005 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1006 } else { // x < 2^32
1007 tmp1
.d
= (double) C1
; // exact conversion
1009 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1011 } else { // if x < 2^53
1012 tmp1
.d
= (double) C1
; // exact conversion
1014 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1016 q
= nr_digits
[x_nr_bits
- 1].digits
;
1018 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1019 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
1022 exp
= x_exp
- 398; // unbiased exponent
1024 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1026 *pfpsf
|= INVALID_EXCEPTION
;
1027 // return Integer Indefinite
1030 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1031 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1032 // so x rounded to an integer may or may not fit in an unsigned 32-bit int
1033 // the cases that do not fit are identified here; the ones that fit
1034 // fall through and will be handled with other cases further,
1035 // under '1 <= q + exp <= 10'
1036 if (x_sign
) { // if n < 0 and q + exp = 10 then x is much less than -1
1037 // => set invalid flag
1038 *pfpsf
|= INVALID_EXCEPTION
;
1039 // return Integer Indefinite
1042 } else { // if n > 0 and q + exp = 10
1043 // if n > 2^32 - 1 then n is too large
1044 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^32 - 1
1045 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6, 1<=q<=16
1046 // <=> C * 10^(11-q) > 0x9fffffff6, 1<=q<=16
1048 // Note: C * 10^(11-q) has 10 or 11 digits; 0x9fffffff6 has 11 digits
1049 tmp64
= C1
* ten2k64
[11 - q
]; // C scaled up to 11-digit int
1050 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1051 if (tmp64
> 0x9fffffff6ull
) {
1053 *pfpsf
|= INVALID_EXCEPTION
;
1054 // return Integer Indefinite
1058 // else cases that can be rounded to a 32-bit unsigned int fall through
1059 // to '1 <= q + exp <= 10'
1060 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1061 // C * 10^(11-q) > 0x9fffffff6 <=>
1062 // C > 0x9fffffff6 * 10^(q-11) where 1 <= q - 11 <= 5
1063 // (scale 2^32-1 up)
1064 // Note: 0x9fffffff6*10^(q-11) has q-1 or q digits, where q <= 16
1065 tmp64
= 0x9fffffff6ull
* ten2k64
[q
- 11];
1068 *pfpsf
|= INVALID_EXCEPTION
;
1069 // return Integer Indefinite
1073 // else cases that can be rounded to a 32-bit int fall through
1074 // to '1 <= q + exp <= 10'
1078 // n is not too large to be converted to int32 if -1 < n < 2^32
1079 // Note: some of the cases tested for above fall through to this point
1080 if ((q
+ exp
) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1087 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
1088 // x <= -1 or 1 <= x <= 2^32 - 1 so if positive, x can be
1089 // rounded to nearest to a 32-bit unsigned integer
1090 if (x_sign
) { // x <= -1
1092 *pfpsf
|= INVALID_EXCEPTION
;
1093 // return Integer Indefinite
1097 // 1 <= x <= 2^32 - 1 so x can be rounded
1098 // to nearest to a 32-bit unsigned integer
1099 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
1100 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
1101 // chop off ind digits from the lower part of C1
1102 // C1 fits in 64 bits
1103 // calculate C* and f*
1104 // C* is actually floor(C*) in this case
1105 // C* and f* need shifting and masking, as shown by
1106 // shiftright128[] and maskhigh128[]
1108 // kx = 10^(-x) = ten2mk64[ind - 1]
1109 // C* = C1 * 10^(-x)
1110 // the approximation of 10^(-x) was rounded up to 54 bits
1111 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
1113 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
1114 fstar
.w
[0] = P128
.w
[0];
1115 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1116 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1117 // C* = floor(C*) (logical right shift; C has p decimal digits,
1118 // correct by Property 1)
1119 // n = C* * 10^(e+x)
1121 // shift right C* by Ex-64 = shiftright128[ind]
1122 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
1123 Cstar
= Cstar
>> shift
;
1124 // determine inexactness of the rounding of C*
1125 // if (0 < f* < 10^(-x)) then
1126 // the result is exact
1127 // else // if (f* > T*) then
1128 // the result is inexact
1129 if (ind
- 1 <= 2) { // fstar.w[1] is 0
1130 if (fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1131 // ten2mk128trunc[ind -1].w[1] is identical to
1132 // ten2mk128[ind -1].w[1]
1134 } // else the result is exact
1135 } else { // if 3 <= ind - 1 <= 14
1136 if (fstar
.w
[1] || fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1137 // ten2mk128trunc[ind -1].w[1] is identical to
1138 // ten2mk128[ind -1].w[1]
1140 } // else the result is exact
1143 res
= Cstar
; // the result is positive
1144 } else if (exp
== 0) {
1147 res
= C1
; // the result is positive
1148 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1149 // res = +C * 10^exp (exact)
1150 res
= C1
* ten2k64
[exp
]; // the result is positive
1156 /*****************************************************************************
1157 * BID64_to_uint32_xceil
1158 ****************************************************************************/
1160 #if DECIMAL_CALL_BY_REFERENCE
1162 bid64_to_uint32_xceil (unsigned int *pres
, UINT64
* px
1163 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1168 bid64_to_uint32_xceil (UINT64 x
1169 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1175 int exp
; // unbiased exponent
1176 // Note: C1 represents x_significand (UINT64)
1178 BID_UI64DOUBLE tmp1
;
1179 unsigned int x_nr_bits
;
1182 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
1186 // check for NaN or Infinity
1187 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
1189 *pfpsf
|= INVALID_EXCEPTION
;
1190 // return Integer Indefinite
1195 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1196 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1197 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
1198 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
1199 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
1200 if (C1
> 9999999999999999ull) { // non-canonical
1205 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
1206 C1
= x
& MASK_BINARY_SIG1
;
1209 // check for zeros (possibly from non-canonical values)
1215 // x is not special and is not zero
1217 // q = nr. of decimal digits in x (1 <= q <= 54)
1218 // determine first the nr. of bits in x
1219 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
1220 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1221 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
1222 tmp1
.d
= (double) (C1
>> 32); // exact conversion
1224 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1225 } else { // x < 2^32
1226 tmp1
.d
= (double) C1
; // exact conversion
1228 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1230 } else { // if x < 2^53
1231 tmp1
.d
= (double) C1
; // exact conversion
1233 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1235 q
= nr_digits
[x_nr_bits
- 1].digits
;
1237 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1238 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
1241 exp
= x_exp
- 398; // unbiased exponent
1243 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1245 *pfpsf
|= INVALID_EXCEPTION
;
1246 // return Integer Indefinite
1249 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1250 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1251 // so x rounded to an integer may or may not fit in an unsigned 32-bit int
1252 // the cases that do not fit are identified here; the ones that fit
1253 // fall through and will be handled with other cases further,
1254 // under '1 <= q + exp <= 10'
1255 if (x_sign
) { // if n < 0 and q + exp = 10 then x is much less than -1
1256 // => set invalid flag
1257 *pfpsf
|= INVALID_EXCEPTION
;
1258 // return Integer Indefinite
1261 } else { // if n > 0 and q + exp = 10
1262 // if n > 2^32 - 1 then n is too large
1263 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^32 - 1
1264 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6, 1<=q<=16
1265 // <=> C * 10^(11-q) > 0x9fffffff6, 1<=q<=16
1267 // Note: C * 10^(11-q) has 10 or 11 digits; 0x9fffffff6 has 11 digits
1268 tmp64
= C1
* ten2k64
[11 - q
]; // C scaled up to 11-digit int
1269 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1270 if (tmp64
> 0x9fffffff6ull
) {
1272 *pfpsf
|= INVALID_EXCEPTION
;
1273 // return Integer Indefinite
1277 // else cases that can be rounded to a 32-bit unsigned int fall through
1278 // to '1 <= q + exp <= 10'
1279 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1280 // C * 10^(11-q) > 0x9fffffff6 <=>
1281 // C > 0x9fffffff6 * 10^(q-11) where 1 <= q - 11 <= 5
1282 // (scale 2^32-1 up)
1283 // Note: 0x9fffffff6*10^(q-11) has q-1 or q digits, where q <= 16
1284 tmp64
= 0x9fffffff6ull
* ten2k64
[q
- 11];
1287 *pfpsf
|= INVALID_EXCEPTION
;
1288 // return Integer Indefinite
1292 // else cases that can be rounded to a 32-bit int fall through
1293 // to '1 <= q + exp <= 10'
1297 // n is not too large to be converted to int32 if -1 < n < 2^32
1298 // Note: some of the cases tested for above fall through to this point
1299 if ((q
+ exp
) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1301 *pfpsf
|= INEXACT_EXCEPTION
;
1308 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
1309 // x <= -1 or 1 <= x < 2^32 so if positive, x can be
1310 // rounded to nearest to a 32-bit unsigned integer
1311 if (x_sign
) { // x <= -1
1313 *pfpsf
|= INVALID_EXCEPTION
;
1314 // return Integer Indefinite
1318 // 1 <= x < 2^32 so x can be rounded
1319 // to nearest to a 32-bit unsigned integer
1320 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
1321 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
1322 // chop off ind digits from the lower part of C1
1323 // C1 fits in 64 bits
1324 // calculate C* and f*
1325 // C* is actually floor(C*) in this case
1326 // C* and f* need shifting and masking, as shown by
1327 // shiftright128[] and maskhigh128[]
1329 // kx = 10^(-x) = ten2mk64[ind - 1]
1330 // C* = C1 * 10^(-x)
1331 // the approximation of 10^(-x) was rounded up to 54 bits
1332 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
1334 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
1335 fstar
.w
[0] = P128
.w
[0];
1336 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1337 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1338 // C* = floor(C*) (logical right shift; C has p decimal digits,
1339 // correct by Property 1)
1340 // n = C* * 10^(e+x)
1342 // shift right C* by Ex-64 = shiftright128[ind]
1343 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
1344 Cstar
= Cstar
>> shift
;
1345 // determine inexactness of the rounding of C*
1346 // if (0 < f* < 10^(-x)) then
1347 // the result is exact
1348 // else // if (f* > T*) then
1349 // the result is inexact
1350 if (ind
- 1 <= 2) { // fstar.w[1] is 0
1351 if (fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1352 // ten2mk128trunc[ind -1].w[1] is identical to
1353 // ten2mk128[ind -1].w[1]
1355 // set the inexact flag
1356 *pfpsf
|= INEXACT_EXCEPTION
;
1357 } // else the result is exact
1358 } else { // if 3 <= ind - 1 <= 14
1359 if (fstar
.w
[1] || fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1360 // ten2mk128trunc[ind -1].w[1] is identical to
1361 // ten2mk128[ind -1].w[1]
1363 // set the inexact flag
1364 *pfpsf
|= INEXACT_EXCEPTION
;
1365 } // else the result is exact
1368 res
= Cstar
; // the result is positive
1369 } else if (exp
== 0) {
1372 res
= C1
; // the result is positive
1373 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1374 // res = +C * 10^exp (exact)
1375 res
= C1
* ten2k64
[exp
]; // the result is positive
1381 /*****************************************************************************
1382 * BID64_to_uint32_int
1383 ****************************************************************************/
1385 #if DECIMAL_CALL_BY_REFERENCE
1387 bid64_to_uint32_int (unsigned int *pres
, UINT64
* px
1388 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM
)
1393 bid64_to_uint32_int (UINT64 x
1394 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM
)
1400 int exp
; // unbiased exponent
1401 // Note: C1 represents x_significand (UINT64)
1403 BID_UI64DOUBLE tmp1
;
1404 unsigned int x_nr_bits
;
1407 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
1410 // check for NaN or Infinity
1411 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
1413 *pfpsf
|= INVALID_EXCEPTION
;
1414 // return Integer Indefinite
1419 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1420 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1421 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
1422 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
1423 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
1424 if (C1
> 9999999999999999ull) { // non-canonical
1429 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
1430 C1
= x
& MASK_BINARY_SIG1
;
1433 // check for zeros (possibly from non-canonical values)
1439 // x is not special and is not zero
1441 // q = nr. of decimal digits in x (1 <= q <= 54)
1442 // determine first the nr. of bits in x
1443 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
1444 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1445 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
1446 tmp1
.d
= (double) (C1
>> 32); // exact conversion
1448 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1449 } else { // x < 2^32
1450 tmp1
.d
= (double) C1
; // exact conversion
1452 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1454 } else { // if x < 2^53
1455 tmp1
.d
= (double) C1
; // exact conversion
1457 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1459 q
= nr_digits
[x_nr_bits
- 1].digits
;
1461 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1462 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
1465 exp
= x_exp
- 398; // unbiased exponent
1467 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1469 *pfpsf
|= INVALID_EXCEPTION
;
1470 // return Integer Indefinite
1473 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1474 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1475 // so x rounded to an integer may or may not fit in an unsigned 32-bit int
1476 // the cases that do not fit are identified here; the ones that fit
1477 // fall through and will be handled with other cases further,
1478 // under '1 <= q + exp <= 10'
1479 if (x_sign
) { // if n < 0 and q + exp = 10 then x is much less than -1
1480 // => set invalid flag
1481 *pfpsf
|= INVALID_EXCEPTION
;
1482 // return Integer Indefinite
1485 } else { // if n > 0 and q + exp = 10
1486 // if n >= 2^32 then n is too large
1487 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
1488 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=16
1489 // <=> C * 10^(11-q) >= 0xa00000000, 1<=q<=16
1491 // Note: C * 10^(11-q) has 10 or 11 digits; 0xa00000000 has 11 digits
1492 tmp64
= C1
* ten2k64
[11 - q
]; // C scaled up to 11-digit int
1493 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1494 if (tmp64
>= 0xa00000000ull
) {
1496 *pfpsf
|= INVALID_EXCEPTION
;
1497 // return Integer Indefinite
1501 // else cases that can be rounded to a 32-bit unsigned int fall through
1502 // to '1 <= q + exp <= 10'
1503 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1504 // C * 10^(11-q) >= 0xa00000000 <=>
1505 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 5
1506 // (scale 2^32-1/2 up)
1507 // Note: 0xa00000000*10^(q-11) has q-1 or q digits, where q <= 16
1508 tmp64
= 0xa00000000ull
* ten2k64
[q
- 11];
1511 *pfpsf
|= INVALID_EXCEPTION
;
1512 // return Integer Indefinite
1516 // else cases that can be rounded to a 32-bit int fall through
1517 // to '1 <= q + exp <= 10'
1521 // n is not too large to be converted to int32 if -1 < n < 2^32
1522 // Note: some of the cases tested for above fall through to this point
1523 if ((q
+ exp
) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1527 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
1528 // x <= -1 or 1 <= x < 2^32 so if positive, x can be
1529 // rounded to nearest to a 32-bit unsigned integer
1530 if (x_sign
) { // x <= -1
1532 *pfpsf
|= INVALID_EXCEPTION
;
1533 // return Integer Indefinite
1537 // 1 <= x < 2^32 so x can be rounded
1538 // to nearest to a 32-bit unsigned integer
1539 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
1540 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
1541 // chop off ind digits from the lower part of C1
1542 // C1 fits in 64 bits
1543 // calculate C* and f*
1544 // C* is actually floor(C*) in this case
1545 // C* and f* need shifting and masking, as shown by
1546 // shiftright128[] and maskhigh128[]
1548 // kx = 10^(-x) = ten2mk64[ind - 1]
1549 // C* = C1 * 10^(-x)
1550 // the approximation of 10^(-x) was rounded up to 54 bits
1551 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
1553 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1554 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1555 // C* = floor(C*) (logical right shift; C has p decimal digits,
1556 // correct by Property 1)
1557 // n = C* * 10^(e+x)
1559 // shift right C* by Ex-64 = shiftright128[ind]
1560 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
1561 Cstar
= Cstar
>> shift
;
1563 res
= Cstar
; // the result is positive
1564 } else if (exp
== 0) {
1567 res
= C1
; // the result is positive
1568 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1569 // res = +C * 10^exp (exact)
1570 res
= C1
* ten2k64
[exp
]; // the result is positive
1576 /*****************************************************************************
1577 * BID64_to_uint32_xint
1578 ****************************************************************************/
1580 #if DECIMAL_CALL_BY_REFERENCE
1582 bid64_to_uint32_xint (unsigned int *pres
, UINT64
* px
1583 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1588 bid64_to_uint32_xint (UINT64 x
1589 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1595 int exp
; // unbiased exponent
1596 // Note: C1 represents x_significand (UINT64)
1598 BID_UI64DOUBLE tmp1
;
1599 unsigned int x_nr_bits
;
1602 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
1606 // check for NaN or Infinity
1607 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
1609 *pfpsf
|= INVALID_EXCEPTION
;
1610 // return Integer Indefinite
1615 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1616 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1617 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
1618 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
1619 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
1620 if (C1
> 9999999999999999ull) { // non-canonical
1625 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
1626 C1
= x
& MASK_BINARY_SIG1
;
1629 // check for zeros (possibly from non-canonical values)
1635 // x is not special and is not zero
1637 // q = nr. of decimal digits in x (1 <= q <= 54)
1638 // determine first the nr. of bits in x
1639 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
1640 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1641 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
1642 tmp1
.d
= (double) (C1
>> 32); // exact conversion
1644 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1645 } else { // x < 2^32
1646 tmp1
.d
= (double) C1
; // exact conversion
1648 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1650 } else { // if x < 2^53
1651 tmp1
.d
= (double) C1
; // exact conversion
1653 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1655 q
= nr_digits
[x_nr_bits
- 1].digits
;
1657 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1658 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
1661 exp
= x_exp
- 398; // unbiased exponent
1663 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1665 *pfpsf
|= INVALID_EXCEPTION
;
1666 // return Integer Indefinite
1669 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1670 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1671 // so x rounded to an integer may or may not fit in an unsigned 32-bit int
1672 // the cases that do not fit are identified here; the ones that fit
1673 // fall through and will be handled with other cases further,
1674 // under '1 <= q + exp <= 10'
1675 if (x_sign
) { // if n < 0 and q + exp = 10 then x is much less than -1
1676 // => set invalid flag
1677 *pfpsf
|= INVALID_EXCEPTION
;
1678 // return Integer Indefinite
1681 } else { // if n > 0 and q + exp = 10
1682 // if n >= 2^32 then n is too large
1683 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
1684 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=16
1685 // <=> C * 10^(11-q) >= 0xa00000000, 1<=q<=16
1687 // Note: C * 10^(11-q) has 10 or 11 digits; 0xa00000000 has 11 digits
1688 tmp64
= C1
* ten2k64
[11 - q
]; // C scaled up to 11-digit int
1689 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1690 if (tmp64
>= 0xa00000000ull
) {
1692 *pfpsf
|= INVALID_EXCEPTION
;
1693 // return Integer Indefinite
1697 // else cases that can be rounded to a 32-bit unsigned int fall through
1698 // to '1 <= q + exp <= 10'
1699 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1700 // C * 10^(11-q) >= 0xa00000000 <=>
1701 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 5
1702 // (scale 2^32-1/2 up)
1703 // Note: 0xa00000000*10^(q-11) has q-1 or q digits, where q <= 16
1704 tmp64
= 0xa00000000ull
* ten2k64
[q
- 11];
1707 *pfpsf
|= INVALID_EXCEPTION
;
1708 // return Integer Indefinite
1712 // else cases that can be rounded to a 32-bit int fall through
1713 // to '1 <= q + exp <= 10'
1717 // n is not too large to be converted to int32 if -1 < n < 2^32
1718 // Note: some of the cases tested for above fall through to this point
1719 if ((q
+ exp
) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1721 *pfpsf
|= INEXACT_EXCEPTION
;
1725 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
1726 // x <= -1 or 1 <= x < 2^32 so if positive, x can be
1727 // rounded to nearest to a 32-bit unsigned integer
1728 if (x_sign
) { // x <= -1
1730 *pfpsf
|= INVALID_EXCEPTION
;
1731 // return Integer Indefinite
1735 // 1 <= x < 2^32 so x can be rounded
1736 // to nearest to a 32-bit unsigned integer
1737 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
1738 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
1739 // chop off ind digits from the lower part of C1
1740 // C1 fits in 64 bits
1741 // calculate C* and f*
1742 // C* is actually floor(C*) in this case
1743 // C* and f* need shifting and masking, as shown by
1744 // shiftright128[] and maskhigh128[]
1746 // kx = 10^(-x) = ten2mk64[ind - 1]
1747 // C* = C1 * 10^(-x)
1748 // the approximation of 10^(-x) was rounded up to 54 bits
1749 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
1751 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
1752 fstar
.w
[0] = P128
.w
[0];
1753 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1754 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1755 // C* = floor(C*) (logical right shift; C has p decimal digits,
1756 // correct by Property 1)
1757 // n = C* * 10^(e+x)
1759 // shift right C* by Ex-64 = shiftright128[ind]
1760 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
1761 Cstar
= Cstar
>> shift
;
1762 // determine inexactness of the rounding of C*
1763 // if (0 < f* < 10^(-x)) then
1764 // the result is exact
1765 // else // if (f* > T*) then
1766 // the result is inexact
1767 if (ind
- 1 <= 2) { // fstar.w[1] is 0
1768 if (fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1769 // ten2mk128trunc[ind -1].w[1] is identical to
1770 // ten2mk128[ind -1].w[1]
1771 // set the inexact flag
1772 *pfpsf
|= INEXACT_EXCEPTION
;
1773 } // else the result is exact
1774 } else { // if 3 <= ind - 1 <= 14
1775 if (fstar
.w
[1] || fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1776 // ten2mk128trunc[ind -1].w[1] is identical to
1777 // ten2mk128[ind -1].w[1]
1778 // set the inexact flag
1779 *pfpsf
|= INEXACT_EXCEPTION
;
1780 } // else the result is exact
1783 res
= Cstar
; // the result is positive
1784 } else if (exp
== 0) {
1787 res
= C1
; // the result is positive
1788 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1789 // res = +C * 10^exp (exact)
1790 res
= C1
* ten2k64
[exp
]; // the result is positive
1796 /*****************************************************************************
1797 * BID64_to_uint32_rninta
1798 ****************************************************************************/
1800 #if DECIMAL_CALL_BY_REFERENCE
1802 bid64_to_uint32_rninta (unsigned int *pres
, UINT64
* px
1803 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1808 bid64_to_uint32_rninta (UINT64 x
1809 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1815 int exp
; // unbiased exponent
1816 // Note: C1 represents x_significand (UINT64)
1818 BID_UI64DOUBLE tmp1
;
1819 unsigned int x_nr_bits
;
1822 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
1825 // check for NaN or Infinity
1826 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
1828 *pfpsf
|= INVALID_EXCEPTION
;
1829 // return Integer Indefinite
1834 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1835 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1836 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
1837 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
1838 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
1839 if (C1
> 9999999999999999ull) { // non-canonical
1844 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
1845 C1
= x
& MASK_BINARY_SIG1
;
1848 // check for zeros (possibly from non-canonical values)
1854 // x is not special and is not zero
1856 // q = nr. of decimal digits in x (1 <= q <= 54)
1857 // determine first the nr. of bits in x
1858 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
1859 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1860 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
1861 tmp1
.d
= (double) (C1
>> 32); // exact conversion
1863 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1864 } else { // x < 2^32
1865 tmp1
.d
= (double) C1
; // exact conversion
1867 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1869 } else { // if x < 2^53
1870 tmp1
.d
= (double) C1
; // exact conversion
1872 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1874 q
= nr_digits
[x_nr_bits
- 1].digits
;
1876 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1877 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
1880 exp
= x_exp
- 398; // unbiased exponent
1882 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1884 *pfpsf
|= INVALID_EXCEPTION
;
1885 // return Integer Indefinite
1888 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1889 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1890 // so x rounded to an integer may or may not fit in an unsigned 32-bit int
1891 // the cases that do not fit are identified here; the ones that fit
1892 // fall through and will be handled with other cases further,
1893 // under '1 <= q + exp <= 10'
1894 if (x_sign
) { // if n < 0 and q + exp = 10 then x is much less than -1/2
1895 // => set invalid flag
1896 *pfpsf
|= INVALID_EXCEPTION
;
1897 // return Integer Indefinite
1900 } else { // if n > 0 and q + exp = 10
1901 // if n >= 2^32 - 1/2 then n is too large
1902 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
1903 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=16
1904 // <=> C * 10^(11-q) >= 0x9fffffffb, 1<=q<=16
1906 // Note: C * 10^(11-q) has 10 or 11 digits; 0x9fffffffb has 11 digits
1907 tmp64
= C1
* ten2k64
[11 - q
]; // C scaled up to 11-digit int
1908 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1909 if (tmp64
>= 0x9fffffffbull
) {
1911 *pfpsf
|= INVALID_EXCEPTION
;
1912 // return Integer Indefinite
1916 // else cases that can be rounded to a 32-bit unsigned int fall through
1917 // to '1 <= q + exp <= 10'
1918 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1919 // C * 10^(11-q) >= 0x9fffffffb <=>
1920 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 5
1921 // (scale 2^32-1/2 up)
1922 // Note: 0x9fffffffb*10^(q-11) has q-1 or q digits, where q <= 16
1923 tmp64
= 0x9fffffffbull
* ten2k64
[q
- 11];
1926 *pfpsf
|= INVALID_EXCEPTION
;
1927 // return Integer Indefinite
1931 // else cases that can be rounded to a 32-bit int fall through
1932 // to '1 <= q + exp <= 10'
1936 // n is not too large to be converted to int32 if -1/2 < n < 2^32 - 1/2
1937 // Note: some of the cases tested for above fall through to this point
1938 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
1942 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
1943 // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
1950 if (C1
< midpoint64
[ind
]) {
1951 res
= 0x00000000; // return 0
1952 } else if (x_sign
) { // n < 0
1954 *pfpsf
|= INVALID_EXCEPTION
;
1955 // return Integer Indefinite
1959 res
= 0x00000001; // return +1
1961 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
1962 // -2^32-1/2 <= x <= -1 or 1 <= x < 2^32-1/2 so if positive, x can be
1963 // rounded to nearest to a 32-bit unsigned integer
1964 if (x_sign
) { // x <= -1
1966 *pfpsf
|= INVALID_EXCEPTION
;
1967 // return Integer Indefinite
1971 // 1 <= x < 2^32-1/2 so x can be rounded
1972 // to nearest to a 32-bit unsigned integer
1973 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
1974 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
1975 // chop off ind digits from the lower part of C1
1976 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
1977 C1
= C1
+ midpoint64
[ind
- 1];
1978 // calculate C* and f*
1979 // C* is actually floor(C*) in this case
1980 // C* and f* need shifting and masking, as shown by
1981 // shiftright128[] and maskhigh128[]
1983 // kx = 10^(-x) = ten2mk64[ind - 1]
1984 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1985 // the approximation of 10^(-x) was rounded up to 54 bits
1986 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
1988 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1989 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1990 // C* = floor(C*) (logical right shift; C has p decimal digits,
1991 // correct by Property 1)
1992 // n = C* * 10^(e+x)
1994 // shift right C* by Ex-64 = shiftright128[ind]
1995 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
1996 Cstar
= Cstar
>> shift
;
1998 // if the result was a midpoint it was rounded away from zero
1999 res
= Cstar
; // the result is positive
2000 } else if (exp
== 0) {
2003 res
= C1
; // the result is positive
2004 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2005 // res = +C * 10^exp (exact)
2006 res
= C1
* ten2k64
[exp
]; // the result is positive
2012 /*****************************************************************************
2013 * BID64_to_uint32_xrninta
2014 ****************************************************************************/
2016 #if DECIMAL_CALL_BY_REFERENCE
2018 bid64_to_uint32_xrninta (unsigned int *pres
, UINT64
* px
2019 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
2024 bid64_to_uint32_xrninta (UINT64 x
2025 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
2031 int exp
; // unbiased exponent
2032 // Note: C1 represents x_significand (UINT64)
2034 BID_UI64DOUBLE tmp1
;
2035 unsigned int x_nr_bits
;
2038 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
2042 // check for NaN or Infinity
2043 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
2045 *pfpsf
|= INVALID_EXCEPTION
;
2046 // return Integer Indefinite
2051 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
2052 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
2053 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
2054 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
2055 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
2056 if (C1
> 9999999999999999ull) { // non-canonical
2061 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
2062 C1
= x
& MASK_BINARY_SIG1
;
2065 // check for zeros (possibly from non-canonical values)
2071 // x is not special and is not zero
2073 // q = nr. of decimal digits in x (1 <= q <= 54)
2074 // determine first the nr. of bits in x
2075 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
2076 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2077 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
2078 tmp1
.d
= (double) (C1
>> 32); // exact conversion
2080 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2081 } else { // x < 2^32
2082 tmp1
.d
= (double) C1
; // exact conversion
2084 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2086 } else { // if x < 2^53
2087 tmp1
.d
= (double) C1
; // exact conversion
2089 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2091 q
= nr_digits
[x_nr_bits
- 1].digits
;
2093 q
= nr_digits
[x_nr_bits
- 1].digits1
;
2094 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
2097 exp
= x_exp
- 398; // unbiased exponent
2099 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
2101 *pfpsf
|= INVALID_EXCEPTION
;
2102 // return Integer Indefinite
2105 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
2106 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
2107 // so x rounded to an integer may or may not fit in an unsigned 32-bit int
2108 // the cases that do not fit are identified here; the ones that fit
2109 // fall through and will be handled with other cases further,
2110 // under '1 <= q + exp <= 10'
2111 if (x_sign
) { // if n < 0 and q + exp = 10 then x is much less than -1/2
2112 // => set invalid flag
2113 *pfpsf
|= INVALID_EXCEPTION
;
2114 // return Integer Indefinite
2117 } else { // if n > 0 and q + exp = 10
2118 // if n >= 2^32 - 1/2 then n is too large
2119 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
2120 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=16
2121 // <=> C * 10^(11-q) >= 0x9fffffffb, 1<=q<=16
2123 // Note: C * 10^(11-q) has 10 or 11 digits; 0x9fffffffb has 11 digits
2124 tmp64
= C1
* ten2k64
[11 - q
]; // C scaled up to 11-digit int
2125 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2126 if (tmp64
>= 0x9fffffffbull
) {
2128 *pfpsf
|= INVALID_EXCEPTION
;
2129 // return Integer Indefinite
2133 // else cases that can be rounded to a 32-bit unsigned int fall through
2134 // to '1 <= q + exp <= 10'
2135 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
2136 // C * 10^(11-q) >= 0x9fffffffb <=>
2137 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 5
2138 // (scale 2^32-1/2 up)
2139 // Note: 0x9fffffffb*10^(q-11) has q-1 or q digits, where q <= 16
2140 tmp64
= 0x9fffffffbull
* ten2k64
[q
- 11];
2143 *pfpsf
|= INVALID_EXCEPTION
;
2144 // return Integer Indefinite
2148 // else cases that can be rounded to a 32-bit int fall through
2149 // to '1 <= q + exp <= 10'
2153 // n is not too large to be converted to int32 if -1/2 < n < 2^32 - 1/2
2154 // Note: some of the cases tested for above fall through to this point
2155 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
2157 *pfpsf
|= INEXACT_EXCEPTION
;
2161 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
2162 // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
2169 if (C1
< midpoint64
[ind
]) {
2170 res
= 0x00000000; // return 0
2171 } else if (x_sign
) { // n < 0
2173 *pfpsf
|= INVALID_EXCEPTION
;
2174 // return Integer Indefinite
2178 res
= 0x00000001; // return +1
2181 *pfpsf
|= INEXACT_EXCEPTION
;
2182 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
2183 // -2^32-1/2 <= x <= -1 or 1 <= x < 2^32-1/2 so if positive, x can be
2184 // rounded to nearest to a 32-bit unsigned integer
2185 if (x_sign
) { // x <= -1
2187 *pfpsf
|= INVALID_EXCEPTION
;
2188 // return Integer Indefinite
2192 // 1 <= x < 2^32-1/2 so x can be rounded
2193 // to nearest to a 32-bit unsigned integer
2194 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
2195 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
2196 // chop off ind digits from the lower part of C1
2197 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
2198 C1
= C1
+ midpoint64
[ind
- 1];
2199 // calculate C* and f*
2200 // C* is actually floor(C*) in this case
2201 // C* and f* need shifting and masking, as shown by
2202 // shiftright128[] and maskhigh128[]
2204 // kx = 10^(-x) = ten2mk64[ind - 1]
2205 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2206 // the approximation of 10^(-x) was rounded up to 54 bits
2207 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
2209 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
2210 fstar
.w
[0] = P128
.w
[0];
2211 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
2212 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
2213 // C* = floor(C*) (logical right shift; C has p decimal digits,
2214 // correct by Property 1)
2215 // n = C* * 10^(e+x)
2217 // shift right C* by Ex-64 = shiftright128[ind]
2218 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
2219 Cstar
= Cstar
>> shift
;
2221 // determine inexactness of the rounding of C*
2222 // if (0 < f* - 1/2 < 10^(-x)) then
2223 // the result is exact
2224 // else // if (f* - 1/2 > T*) then
2225 // the result is inexact
2226 if (ind
- 1 <= 2) { // fstar.w[1] is 0
2227 if (fstar
.w
[0] > 0x8000000000000000ull
) {
2228 // f* > 1/2 and the result may be exact
2229 tmp64
= fstar
.w
[0] - 0x8000000000000000ull
; // f* - 1/2
2230 if ((tmp64
> ten2mk128trunc
[ind
- 1].w
[1])) {
2231 // ten2mk128trunc[ind -1].w[1] is identical to
2232 // ten2mk128[ind -1].w[1]
2233 // set the inexact flag
2234 *pfpsf
|= INEXACT_EXCEPTION
;
2235 } // else the result is exact
2236 } else { // the result is inexact; f2* <= 1/2
2237 // set the inexact flag
2238 *pfpsf
|= INEXACT_EXCEPTION
;
2240 } else { // if 3 <= ind - 1 <= 14
2241 if (fstar
.w
[1] > onehalf128
[ind
- 1] ||
2242 (fstar
.w
[1] == onehalf128
[ind
- 1] && fstar
.w
[0])) {
2243 // f2* > 1/2 and the result may be exact
2244 // Calculate f2* - 1/2
2245 tmp64
= fstar
.w
[1] - onehalf128
[ind
- 1];
2246 if (tmp64
|| fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
2247 // ten2mk128trunc[ind -1].w[1] is identical to
2248 // ten2mk128[ind -1].w[1]
2249 // set the inexact flag
2250 *pfpsf
|= INEXACT_EXCEPTION
;
2251 } // else the result is exact
2252 } else { // the result is inexact; f2* <= 1/2
2253 // set the inexact flag
2254 *pfpsf
|= INEXACT_EXCEPTION
;
2258 // if the result was a midpoint it was rounded away from zero
2259 res
= Cstar
; // the result is positive
2260 } else if (exp
== 0) {
2263 res
= C1
; // the result is positive
2264 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2265 // res = +C * 10^exp (exact)
2266 res
= C1
* ten2k64
[exp
]; // the result is positive