1 /* Copyright (C) 2007 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 2, or (at your option) any later
10 In addition to the permissions in the GNU General Public License, the
11 Free Software Foundation gives you unlimited permission to link the
12 compiled version of this file into combinations with other programs,
13 and to distribute those combinations without any restriction coming
14 from the use of this file. (The General Public License restrictions
15 do apply in other respects; for example, they cover modification of
16 the file, and distribution when not linked into a combine
19 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
20 WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
24 You should have received a copy of the GNU General Public License
25 along with GCC; see the file COPYING. If not, write to the Free
26 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
29 #include "bid_internal.h"
31 /*****************************************************************************
32 * BID64_to_uint64_rnint
33 ****************************************************************************/
35 #if DECIMAL_CALL_BY_REFERENCE
37 bid64_to_uint64_rnint (UINT64
* pres
, UINT64
* px
38 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
43 bid64_to_uint64_rnint (UINT64 x
44 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
50 int exp
; // unbiased exponent
51 // Note: C1 represents x_significand (UINT64)
53 unsigned int x_nr_bits
;
57 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
61 // check for NaN or Infinity
62 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
64 *pfpsf
|= INVALID_EXCEPTION
;
65 // return Integer Indefinite
66 res
= 0x8000000000000000ull
;
70 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
71 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
72 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
73 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
74 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
75 if (C1
> 9999999999999999ull) { // non-canonical
80 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
81 C1
= x
& MASK_BINARY_SIG1
;
84 // check for zeros (possibly from non-canonical values)
87 res
= 0x0000000000000000ull
;
90 // x is not special and is not zero
92 // q = nr. of decimal digits in x (1 <= q <= 54)
93 // determine first the nr. of bits in x
94 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
95 // split the 64-bit value in two 32-bit halves to avoid rounding errors
96 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
97 tmp1
.d
= (double) (C1
>> 32); // exact conversion
99 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
101 tmp1
.d
= (double) C1
; // exact conversion
103 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
105 } else { // if x < 2^53
106 tmp1
.d
= (double) C1
; // exact conversion
108 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
110 q
= nr_digits
[x_nr_bits
- 1].digits
;
112 q
= nr_digits
[x_nr_bits
- 1].digits1
;
113 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
116 exp
= x_exp
- 398; // unbiased exponent
118 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
120 *pfpsf
|= INVALID_EXCEPTION
;
121 // return Integer Indefinite
122 res
= 0x8000000000000000ull
;
124 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
125 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
126 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
127 // the cases that do not fit are identified here; the ones that fit
128 // fall through and will be handled with other cases further,
129 // under '1 <= q + exp <= 20'
130 if (x_sign
) { // if n < 0 and q + exp = 20 then x is much less than -1/2
131 // => set invalid flag
132 *pfpsf
|= INVALID_EXCEPTION
;
133 // return Integer Indefinite
134 res
= 0x8000000000000000ull
;
136 } else { // if n > 0 and q + exp = 20
137 // if n >= 2^64 - 1/2 then n is too large
138 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
139 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
140 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
141 // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=16
143 // C * 10^20 >= 0x9fffffffffffffffb
144 __mul_128x64_to_128 (C
, C1
, ten2k128
[0]); // 10^20 * C
146 (C
.w
[1] == 0x09 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
148 *pfpsf
|= INVALID_EXCEPTION
;
149 // return Integer Indefinite
150 res
= 0x8000000000000000ull
;
153 // else cases that can be rounded to a 64-bit int fall through
154 // to '1 <= q + exp <= 20'
155 } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
156 // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffffb
158 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[21 - q
]);
160 (C
.w
[1] == 0x09 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
162 *pfpsf
|= INVALID_EXCEPTION
;
163 // return Integer Indefinite
164 res
= 0x8000000000000000ull
;
167 // else cases that can be rounded to a 64-bit int fall through
168 // to '1 <= q + exp <= 20'
172 // n is not too large to be converted to int64 if -1/2 <= n < 2^64 - 1/2
173 // Note: some of the cases tested for above fall through to this point
174 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
176 res
= 0x0000000000000000ull
;
178 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
179 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
185 ind
= q
- 1; // 0 <= ind <= 15
186 if (C1
<= midpoint64
[ind
]) {
187 res
= 0x0000000000000000ull
; // return 0
188 } else if (!x_sign
) { // n > 0
189 res
= 0x0000000000000001ull
; // return +1
191 res
= 0x8000000000000000ull
;
192 *pfpsf
|= INVALID_EXCEPTION
;
195 } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
196 // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
197 // to nearest to a 64-bit unsigned signed integer
198 if (x_sign
) { // x <= -1
200 *pfpsf
|= INVALID_EXCEPTION
;
201 // return Integer Indefinite
202 res
= 0x8000000000000000ull
;
205 // 1 <= x < 2^64-1/2 so x can be rounded
206 // to nearest to a 64-bit unsigned integer
207 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
208 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
209 // chop off ind digits from the lower part of C1
210 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
211 C1
= C1
+ midpoint64
[ind
- 1];
212 // calculate C* and f*
213 // C* is actually floor(C*) in this case
214 // C* and f* need shifting and masking, as shown by
215 // shiftright128[] and maskhigh128[]
217 // kx = 10^(-x) = ten2mk64[ind - 1]
218 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
219 // the approximation of 10^(-x) was rounded up to 54 bits
220 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
222 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
223 fstar
.w
[0] = P128
.w
[0];
224 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
225 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
226 // if (0 < f* < 10^(-x)) then the result is a midpoint
227 // if floor(C*) is even then C* = floor(C*) - logical right
228 // shift; C* has p decimal digits, correct by Prop. 1)
229 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
230 // shift; C* has p decimal digits, correct by Pr. 1)
232 // C* = floor(C*) (logical right shift; C has p decimal digits,
233 // correct by Property 1)
236 // shift right C* by Ex-64 = shiftright128[ind]
237 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
238 Cstar
= Cstar
>> shift
;
240 // if the result was a midpoint it was rounded away from zero, so
241 // it will need a correction
242 // check for midpoints
243 if ((fstar
.w
[1] == 0) && fstar
.w
[0] &&
244 (fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[1])) {
245 // ten2mk128trunc[ind -1].w[1] is identical to
246 // ten2mk128[ind -1].w[1]
247 // the result is a midpoint; round to nearest
248 if (Cstar
& 0x01) { // Cstar is odd; MP in [EVEN, ODD]
249 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
250 Cstar
--; // Cstar is now even
251 } // else MP in [ODD, EVEN]
253 res
= Cstar
; // the result is positive
254 } else if (exp
== 0) {
257 res
= C1
; // the result is positive
258 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
259 // res = +C * 10^exp (exact)
260 res
= C1
* ten2k64
[exp
]; // the result is positive
266 /*****************************************************************************
267 * BID64_to_uint64_xrnint
268 ****************************************************************************/
270 #if DECIMAL_CALL_BY_REFERENCE
272 bid64_to_uint64_xrnint (UINT64
* pres
, UINT64
* px
273 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
278 bid64_to_uint64_xrnint (UINT64 x
279 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
285 int exp
; // unbiased exponent
286 // Note: C1 represents x_significand (UINT64)
289 unsigned int x_nr_bits
;
293 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
297 // check for NaN or Infinity
298 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
300 *pfpsf
|= INVALID_EXCEPTION
;
301 // return Integer Indefinite
302 res
= 0x8000000000000000ull
;
306 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
307 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
308 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
309 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
310 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
311 if (C1
> 9999999999999999ull) { // non-canonical
316 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
317 C1
= x
& MASK_BINARY_SIG1
;
320 // check for zeros (possibly from non-canonical values)
323 res
= 0x0000000000000000ull
;
326 // x is not special and is not zero
328 // q = nr. of decimal digits in x (1 <= q <= 54)
329 // determine first the nr. of bits in x
330 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
331 // split the 64-bit value in two 32-bit halves to avoid rounding errors
332 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
333 tmp1
.d
= (double) (C1
>> 32); // exact conversion
335 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
337 tmp1
.d
= (double) C1
; // exact conversion
339 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
341 } else { // if x < 2^53
342 tmp1
.d
= (double) C1
; // exact conversion
344 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
346 q
= nr_digits
[x_nr_bits
- 1].digits
;
348 q
= nr_digits
[x_nr_bits
- 1].digits1
;
349 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
352 exp
= x_exp
- 398; // unbiased exponent
354 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
356 *pfpsf
|= INVALID_EXCEPTION
;
357 // return Integer Indefinite
358 res
= 0x8000000000000000ull
;
360 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
361 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
362 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
363 // the cases that do not fit are identified here; the ones that fit
364 // fall through and will be handled with other cases further,
365 // under '1 <= q + exp <= 20'
366 if (x_sign
) { // if n < 0 and q + exp = 20 then x is much less than -1/2
367 // => set invalid flag
368 *pfpsf
|= INVALID_EXCEPTION
;
369 // return Integer Indefinite
370 res
= 0x8000000000000000ull
;
372 } else { // if n > 0 and q + exp = 20
373 // if n >= 2^64 - 1/2 then n is too large
374 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
375 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
376 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
377 // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=16
379 // C * 10^20 >= 0x9fffffffffffffffb
380 __mul_128x64_to_128 (C
, C1
, ten2k128
[0]); // 10^20 * C
382 (C
.w
[1] == 0x09 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
384 *pfpsf
|= INVALID_EXCEPTION
;
385 // return Integer Indefinite
386 res
= 0x8000000000000000ull
;
389 // else cases that can be rounded to a 64-bit int fall through
390 // to '1 <= q + exp <= 20'
391 } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
392 // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffffb
394 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[21 - q
]);
396 (C
.w
[1] == 0x09 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
398 *pfpsf
|= INVALID_EXCEPTION
;
399 // return Integer Indefinite
400 res
= 0x8000000000000000ull
;
403 // else cases that can be rounded to a 64-bit int fall through
404 // to '1 <= q + exp <= 20'
408 // n is not too large to be converted to int64 if -1/2 <= n < 2^64 - 1/2
409 // Note: some of the cases tested for above fall through to this point
410 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
412 *pfpsf
|= INEXACT_EXCEPTION
;
414 res
= 0x0000000000000000ull
;
416 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
417 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
423 ind
= q
- 1; // 0 <= ind <= 15
424 if (C1
<= midpoint64
[ind
]) {
425 res
= 0x0000000000000000ull
; // return 0
426 } else if (!x_sign
) { // n > 0
427 res
= 0x0000000000000001ull
; // return +1
429 res
= 0x8000000000000000ull
;
430 *pfpsf
|= INVALID_EXCEPTION
;
434 *pfpsf
|= INEXACT_EXCEPTION
;
435 } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
436 // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
437 // to nearest to a 64-bit unsigned signed integer
438 if (x_sign
) { // x <= -1
440 *pfpsf
|= INVALID_EXCEPTION
;
441 // return Integer Indefinite
442 res
= 0x8000000000000000ull
;
445 // 1 <= x < 2^64-1/2 so x can be rounded
446 // to nearest to a 64-bit unsigned integer
447 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
448 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
449 // chop off ind digits from the lower part of C1
450 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
451 C1
= C1
+ midpoint64
[ind
- 1];
452 // calculate C* and f*
453 // C* is actually floor(C*) in this case
454 // C* and f* need shifting and masking, as shown by
455 // shiftright128[] and maskhigh128[]
457 // kx = 10^(-x) = ten2mk64[ind - 1]
458 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
459 // the approximation of 10^(-x) was rounded up to 54 bits
460 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
462 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
463 fstar
.w
[0] = P128
.w
[0];
464 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
465 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
466 // if (0 < f* < 10^(-x)) then the result is a midpoint
467 // if floor(C*) is even then C* = floor(C*) - logical right
468 // shift; C* has p decimal digits, correct by Prop. 1)
469 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
470 // shift; C* has p decimal digits, correct by Pr. 1)
472 // C* = floor(C*) (logical right shift; C has p decimal digits,
473 // correct by Property 1)
476 // shift right C* by Ex-64 = shiftright128[ind]
477 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
478 Cstar
= Cstar
>> shift
;
479 // determine inexactness of the rounding of C*
480 // if (0 < f* - 1/2 < 10^(-x)) then
481 // the result is exact
482 // else // if (f* - 1/2 > T*) then
483 // the result is inexact
484 if (ind
- 1 <= 2) { // fstar.w[1] is 0
485 if (fstar
.w
[0] > 0x8000000000000000ull
) {
486 // f* > 1/2 and the result may be exact
487 tmp64
= fstar
.w
[0] - 0x8000000000000000ull
; // f* - 1/2
488 if ((tmp64
> ten2mk128trunc
[ind
- 1].w
[1])) {
489 // ten2mk128trunc[ind -1].w[1] is identical to
490 // ten2mk128[ind -1].w[1]
491 // set the inexact flag
492 *pfpsf
|= INEXACT_EXCEPTION
;
493 } // else the result is exact
494 } else { // the result is inexact; f2* <= 1/2
495 // set the inexact flag
496 *pfpsf
|= INEXACT_EXCEPTION
;
498 } else { // if 3 <= ind - 1 <= 14
499 if (fstar
.w
[1] > onehalf128
[ind
- 1] ||
500 (fstar
.w
[1] == onehalf128
[ind
- 1] && fstar
.w
[0])) {
501 // f2* > 1/2 and the result may be exact
502 // Calculate f2* - 1/2
503 tmp64
= fstar
.w
[1] - onehalf128
[ind
- 1];
504 if (tmp64
|| fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
505 // ten2mk128trunc[ind -1].w[1] is identical to
506 // ten2mk128[ind -1].w[1]
507 // set the inexact flag
508 *pfpsf
|= INEXACT_EXCEPTION
;
509 } // else the result is exact
510 } else { // the result is inexact; f2* <= 1/2
511 // set the inexact flag
512 *pfpsf
|= INEXACT_EXCEPTION
;
516 // if the result was a midpoint it was rounded away from zero, so
517 // it will need a correction
518 // check for midpoints
519 if ((fstar
.w
[1] == 0) && fstar
.w
[0] &&
520 (fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[1])) {
521 // ten2mk128trunc[ind -1].w[1] is identical to
522 // ten2mk128[ind -1].w[1]
523 // the result is a midpoint; round to nearest
524 if (Cstar
& 0x01) { // Cstar is odd; MP in [EVEN, ODD]
525 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
526 Cstar
--; // Cstar is now even
527 } // else MP in [ODD, EVEN]
529 res
= Cstar
; // the result is positive
530 } else if (exp
== 0) {
533 res
= C1
; // the result is positive
534 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
535 // res = +C * 10^exp (exact)
536 res
= C1
* ten2k64
[exp
]; // the result is positive
542 /*****************************************************************************
543 * BID64_to_uint64_floor
544 ****************************************************************************/
546 #if DECIMAL_CALL_BY_REFERENCE
548 bid64_to_uint64_floor (UINT64
* pres
, UINT64
* px
549 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
554 bid64_to_uint64_floor (UINT64 x
555 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
561 int exp
; // unbiased exponent
562 // Note: C1 represents x_significand (UINT64)
564 unsigned int x_nr_bits
;
568 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
571 // check for NaN or Infinity
572 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
574 *pfpsf
|= INVALID_EXCEPTION
;
575 // return Integer Indefinite
576 res
= 0x8000000000000000ull
;
580 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
581 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
582 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
583 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
584 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
585 if (C1
> 9999999999999999ull) { // non-canonical
590 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
591 C1
= x
& MASK_BINARY_SIG1
;
594 // check for zeros (possibly from non-canonical values)
597 res
= 0x0000000000000000ull
;
600 // x is not special and is not zero
602 if (x_sign
) { // if n < 0 the conversion is invalid
604 *pfpsf
|= INVALID_EXCEPTION
;
605 // return Integer Indefinite
606 res
= 0x8000000000000000ull
;
609 // q = nr. of decimal digits in x (1 <= q <= 54)
610 // determine first the nr. of bits in x
611 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
612 // split the 64-bit value in two 32-bit halves to avoid rounding errors
613 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
614 tmp1
.d
= (double) (C1
>> 32); // exact conversion
616 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
618 tmp1
.d
= (double) C1
; // exact conversion
620 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
622 } else { // if x < 2^53
623 tmp1
.d
= (double) C1
; // exact conversion
625 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
627 q
= nr_digits
[x_nr_bits
- 1].digits
;
629 q
= nr_digits
[x_nr_bits
- 1].digits1
;
630 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
633 exp
= x_exp
- 398; // unbiased exponent
635 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
637 *pfpsf
|= INVALID_EXCEPTION
;
638 // return Integer Indefinite
639 res
= 0x8000000000000000ull
;
641 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
642 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
643 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
644 // the cases that do not fit are identified here; the ones that fit
645 // fall through and will be handled with other cases further,
646 // under '1 <= q + exp <= 20'
647 // n > 0 and q + exp = 20
648 // if n >= 2^64 then n is too large
649 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
650 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
651 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65)
652 // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=16
654 // C * 10^20 >= 0xa0000000000000000
655 __mul_128x64_to_128 (C
, C1
, ten2k128
[0]); // 10^20 * C
656 if (C
.w
[1] >= 0x0a) {
657 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
659 *pfpsf
|= INVALID_EXCEPTION
;
660 // return Integer Indefinite
661 res
= 0x8000000000000000ull
;
664 // else cases that can be rounded to a 64-bit int fall through
665 // to '1 <= q + exp <= 20'
666 } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
667 // Note: C * 10^(21-q) has 20 or 21 digits; 0xa0000000000000000
669 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[21 - q
]);
670 if (C
.w
[1] >= 0x0a) {
671 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
673 *pfpsf
|= INVALID_EXCEPTION
;
674 // return Integer Indefinite
675 res
= 0x8000000000000000ull
;
678 // else cases that can be rounded to a 64-bit int fall through
679 // to '1 <= q + exp <= 20'
682 // n is not too large to be converted to int64 if -1 < n < 2^64
683 // Note: some of the cases tested for above fall through to this point
684 if ((q
+ exp
) <= 0) { // n = +0.[0...0]c(0)c(1)...c(q-1)
686 res
= 0x0000000000000000ull
;
688 } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
689 // 1 <= x < 2^64 so x can be rounded
690 // to nearest to a 64-bit unsigned signed integer
691 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
692 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
693 // chop off ind digits from the lower part of C1
694 // C1 fits in 64 bits
695 // calculate C* and f*
696 // C* is actually floor(C*) in this case
697 // C* and f* need shifting and masking, as shown by
698 // shiftright128[] and maskhigh128[]
700 // kx = 10^(-x) = ten2mk64[ind - 1]
702 // the approximation of 10^(-x) was rounded up to 54 bits
703 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
705 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
706 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
707 // C* = floor(C*) (logical right shift; C has p decimal digits,
708 // correct by Property 1)
711 // shift right C* by Ex-64 = shiftright128[ind]
712 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
713 Cstar
= Cstar
>> shift
;
714 res
= Cstar
; // the result is positive
715 } else if (exp
== 0) {
718 res
= C1
; // the result is positive
719 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
720 // res = +C * 10^exp (exact)
721 res
= C1
* ten2k64
[exp
]; // the result is positive
727 /*****************************************************************************
728 * BID64_to_uint64_xfloor
729 ****************************************************************************/
731 #if DECIMAL_CALL_BY_REFERENCE
733 bid64_to_uint64_xfloor (UINT64
* pres
, UINT64
* px
734 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
739 bid64_to_uint64_xfloor (UINT64 x
740 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
746 int exp
; // unbiased exponent
747 // Note: C1 represents x_significand (UINT64)
749 unsigned int x_nr_bits
;
753 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
757 // check for NaN or Infinity
758 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
760 *pfpsf
|= INVALID_EXCEPTION
;
761 // return Integer Indefinite
762 res
= 0x8000000000000000ull
;
766 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
767 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
768 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
769 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
770 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
771 if (C1
> 9999999999999999ull) { // non-canonical
776 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
777 C1
= x
& MASK_BINARY_SIG1
;
780 // check for zeros (possibly from non-canonical values)
783 res
= 0x0000000000000000ull
;
786 // x is not special and is not zero
788 if (x_sign
) { // if n < 0 the conversion is invalid
790 *pfpsf
|= INVALID_EXCEPTION
;
791 // return Integer Indefinite
792 res
= 0x8000000000000000ull
;
795 // q = nr. of decimal digits in x (1 <= q <= 54)
796 // determine first the nr. of bits in x
797 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
798 // split the 64-bit value in two 32-bit halves to avoid rounding errors
799 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
800 tmp1
.d
= (double) (C1
>> 32); // exact conversion
802 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
804 tmp1
.d
= (double) C1
; // exact conversion
806 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
808 } else { // if x < 2^53
809 tmp1
.d
= (double) C1
; // exact conversion
811 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
813 q
= nr_digits
[x_nr_bits
- 1].digits
;
815 q
= nr_digits
[x_nr_bits
- 1].digits1
;
816 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
819 exp
= x_exp
- 398; // unbiased exponent
821 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
823 *pfpsf
|= INVALID_EXCEPTION
;
824 // return Integer Indefinite
825 res
= 0x8000000000000000ull
;
827 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
828 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
829 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
830 // the cases that do not fit are identified here; the ones that fit
831 // fall through and will be handled with other cases further,
832 // under '1 <= q + exp <= 20'
833 // n > 0 and q + exp = 20
834 // if n >= 2^64 then n is too large
835 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
836 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
837 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65)
838 // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=16
840 // C * 10^20 >= 0xa0000000000000000
841 __mul_128x64_to_128 (C
, C1
, ten2k128
[0]); // 10^20 * C
842 if (C
.w
[1] >= 0x0a) {
843 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
845 *pfpsf
|= INVALID_EXCEPTION
;
846 // return Integer Indefinite
847 res
= 0x8000000000000000ull
;
850 // else cases that can be rounded to a 64-bit int fall through
851 // to '1 <= q + exp <= 20'
852 } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
853 // Note: C * 10^(21-q) has 20 or 21 digits; 0xa0000000000000000
855 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[21 - q
]);
856 if (C
.w
[1] >= 0x0a) {
857 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
859 *pfpsf
|= INVALID_EXCEPTION
;
860 // return Integer Indefinite
861 res
= 0x8000000000000000ull
;
864 // else cases that can be rounded to a 64-bit int fall through
865 // to '1 <= q + exp <= 20'
868 // n is not too large to be converted to int64 if -1 < n < 2^64
869 // Note: some of the cases tested for above fall through to this point
870 if ((q
+ exp
) <= 0) { // n = +0.[0...0]c(0)c(1)...c(q-1)
872 *pfpsf
|= INEXACT_EXCEPTION
;
874 res
= 0x0000000000000000ull
;
876 } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
877 // 1 <= x < 2^64 so x can be rounded
878 // to nearest to a 64-bit unsigned signed integer
879 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
880 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
881 // chop off ind digits from the lower part of C1
882 // C1 fits in 64 bits
883 // calculate C* and f*
884 // C* is actually floor(C*) in this case
885 // C* and f* need shifting and masking, as shown by
886 // shiftright128[] and maskhigh128[]
888 // kx = 10^(-x) = ten2mk64[ind - 1]
890 // the approximation of 10^(-x) was rounded up to 54 bits
891 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
893 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
894 fstar
.w
[0] = P128
.w
[0];
895 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
896 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
897 // C* = floor(C*) (logical right shift; C has p decimal digits,
898 // correct by Property 1)
901 // shift right C* by Ex-64 = shiftright128[ind]
902 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
903 Cstar
= Cstar
>> shift
;
904 // determine inexactness of the rounding of C*
905 // if (0 < f* < 10^(-x)) then
906 // the result is exact
907 // else // if (f* > T*) then
908 // the result is inexact
909 if (ind
- 1 <= 2) { // fstar.w[1] is 0
910 if (fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
911 // ten2mk128trunc[ind -1].w[1] is identical to
912 // ten2mk128[ind -1].w[1]
913 // set the inexact flag
914 *pfpsf
|= INEXACT_EXCEPTION
;
915 } // else the result is exact
916 } else { // if 3 <= ind - 1 <= 14
917 if (fstar
.w
[1] || fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
918 // ten2mk128trunc[ind -1].w[1] is identical to
919 // ten2mk128[ind -1].w[1]
920 // set the inexact flag
921 *pfpsf
|= INEXACT_EXCEPTION
;
922 } // else the result is exact
925 res
= Cstar
; // the result is positive
926 } else if (exp
== 0) {
929 res
= C1
; // the result is positive
930 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
931 // res = +C * 10^exp (exact)
932 res
= C1
* ten2k64
[exp
]; // the result is positive
938 /*****************************************************************************
939 * BID64_to_uint64_ceil
940 ****************************************************************************/
942 #if DECIMAL_CALL_BY_REFERENCE
944 bid64_to_uint64_ceil (UINT64
* pres
, UINT64
* px
945 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
950 bid64_to_uint64_ceil (UINT64 x
951 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
957 int exp
; // unbiased exponent
958 // Note: C1 represents x_significand (UINT64)
960 unsigned int x_nr_bits
;
964 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
968 // check for NaN or Infinity
969 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
971 *pfpsf
|= INVALID_EXCEPTION
;
972 // return Integer Indefinite
973 res
= 0x8000000000000000ull
;
977 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
978 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
979 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
980 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
981 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
982 if (C1
> 9999999999999999ull) { // non-canonical
987 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
988 C1
= x
& MASK_BINARY_SIG1
;
991 // check for zeros (possibly from non-canonical values)
994 res
= 0x0000000000000000ull
;
997 // x is not special and is not zero
999 // q = nr. of decimal digits in x (1 <= q <= 54)
1000 // determine first the nr. of bits in x
1001 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
1002 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1003 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
1004 tmp1
.d
= (double) (C1
>> 32); // exact conversion
1006 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1007 } else { // x < 2^32
1008 tmp1
.d
= (double) C1
; // exact conversion
1010 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1012 } else { // if x < 2^53
1013 tmp1
.d
= (double) C1
; // exact conversion
1015 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1017 q
= nr_digits
[x_nr_bits
- 1].digits
;
1019 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1020 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
1023 exp
= x_exp
- 398; // unbiased exponent
1025 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
1027 *pfpsf
|= INVALID_EXCEPTION
;
1028 // return Integer Indefinite
1029 res
= 0x8000000000000000ull
;
1031 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
1032 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1033 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
1034 // the cases that do not fit are identified here; the ones that fit
1035 // fall through and will be handled with other cases further,
1036 // under '1 <= q + exp <= 20'
1037 if (x_sign
) { // if n < 0 and q + exp = 20 then x is much less than -1
1038 // => set invalid flag
1039 *pfpsf
|= INVALID_EXCEPTION
;
1040 // return Integer Indefinite
1041 res
= 0x8000000000000000ull
;
1043 } else { // if n > 0 and q + exp = 20
1044 // if n > 2^64 - 1 then n is too large
1045 // <=> c(0)c(1)...c(19).c(20)...c(q-1) > 2^64 - 1
1046 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 > 2^64 - 1
1047 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65 - 2)
1048 // <=> C * 10^(21-q) > 0x9fffffffffffffff6, 1<=q<=16
1050 // C * 10^20 > 0x9fffffffffffffff6
1051 __mul_128x64_to_128 (C
, C1
, ten2k128
[0]); // 10^20 * C
1052 if (C
.w
[1] > 0x09 ||
1053 (C
.w
[1] == 0x09 && C
.w
[0] > 0xfffffffffffffff6ull
)) {
1055 *pfpsf
|= INVALID_EXCEPTION
;
1056 // return Integer Indefinite
1057 res
= 0x8000000000000000ull
;
1060 // else cases that can be rounded to a 64-bit int fall through
1061 // to '1 <= q + exp <= 20'
1062 } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
1063 // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffff6
1065 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[21 - q
]);
1066 if (C
.w
[1] > 0x09 ||
1067 (C
.w
[1] == 0x09 && C
.w
[0] > 0xfffffffffffffff6ull
)) {
1069 *pfpsf
|= INVALID_EXCEPTION
;
1070 // return Integer Indefinite
1071 res
= 0x8000000000000000ull
;
1074 // else cases that can be rounded to a 64-bit int fall through
1075 // to '1 <= q + exp <= 20'
1079 // n is not too large to be converted to int64 if -1 < n < 2^64
1080 // Note: some of the cases tested for above fall through to this point
1081 if ((q
+ exp
) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1084 res
= 0x0000000000000000ull
;
1086 res
= 0x0000000000000001ull
;
1088 } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
1089 // x <= -1 or 1 <= x <= 2^64 - 1 so if positive x can be rounded
1090 // to nearest to a 64-bit unsigned signed integer
1091 if (x_sign
) { // x <= -1
1093 *pfpsf
|= INVALID_EXCEPTION
;
1094 // return Integer Indefinite
1095 res
= 0x8000000000000000ull
;
1098 // 1 <= x <= 2^64 - 1 so x can be rounded
1099 // to nearest to a 64-bit unsigned integer
1100 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
1101 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
1102 // chop off ind digits from the lower part of C1
1103 // C1 fits in 64 bits
1104 // calculate C* and f*
1105 // C* is actually floor(C*) in this case
1106 // C* and f* need shifting and masking, as shown by
1107 // shiftright128[] and maskhigh128[]
1109 // kx = 10^(-x) = ten2mk64[ind - 1]
1110 // C* = C1 * 10^(-x)
1111 // the approximation of 10^(-x) was rounded up to 54 bits
1112 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
1114 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
1115 fstar
.w
[0] = P128
.w
[0];
1116 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1117 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1118 // C* = floor(C*) (logical right shift; C has p decimal digits,
1119 // correct by Property 1)
1120 // n = C* * 10^(e+x)
1122 // shift right C* by Ex-64 = shiftright128[ind]
1123 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
1124 Cstar
= Cstar
>> shift
;
1125 // determine inexactness of the rounding of C*
1126 // if (0 < f* < 10^(-x)) then
1127 // the result is exact
1128 // else // if (f* > T*) then
1129 // the result is inexact
1130 if (ind
- 1 <= 2) { // fstar.w[1] is 0
1131 if (fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1132 // ten2mk128trunc[ind -1].w[1] is identical to
1133 // ten2mk128[ind -1].w[1]
1135 } // else the result is exact
1136 } else { // if 3 <= ind - 1 <= 14
1137 if (fstar
.w
[1] || fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1138 // ten2mk128trunc[ind -1].w[1] is identical to
1139 // ten2mk128[ind -1].w[1]
1141 } // else the result is exact
1144 res
= Cstar
; // the result is positive
1145 } else if (exp
== 0) {
1148 res
= C1
; // the result is positive
1149 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1150 // res = +C * 10^exp (exact)
1151 res
= C1
* ten2k64
[exp
]; // the result is positive
1157 /*****************************************************************************
1158 * BID64_to_uint64_xceil
1159 ****************************************************************************/
1161 #if DECIMAL_CALL_BY_REFERENCE
1163 bid64_to_uint64_xceil (UINT64
* pres
, UINT64
* px
1164 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1169 bid64_to_uint64_xceil (UINT64 x
1170 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1176 int exp
; // unbiased exponent
1177 // Note: C1 represents x_significand (UINT64)
1178 BID_UI64DOUBLE tmp1
;
1179 unsigned int x_nr_bits
;
1183 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
1187 // check for NaN or Infinity
1188 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
1190 *pfpsf
|= INVALID_EXCEPTION
;
1191 // return Integer Indefinite
1192 res
= 0x8000000000000000ull
;
1196 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1197 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1198 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
1199 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
1200 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
1201 if (C1
> 9999999999999999ull) { // non-canonical
1206 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
1207 C1
= x
& MASK_BINARY_SIG1
;
1210 // check for zeros (possibly from non-canonical values)
1213 res
= 0x0000000000000000ull
;
1216 // x is not special and is not zero
1218 // q = nr. of decimal digits in x (1 <= q <= 54)
1219 // determine first the nr. of bits in x
1220 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
1221 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1222 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
1223 tmp1
.d
= (double) (C1
>> 32); // exact conversion
1225 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1226 } else { // x < 2^32
1227 tmp1
.d
= (double) C1
; // exact conversion
1229 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1231 } else { // if x < 2^53
1232 tmp1
.d
= (double) C1
; // exact conversion
1234 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1236 q
= nr_digits
[x_nr_bits
- 1].digits
;
1238 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1239 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
1242 exp
= x_exp
- 398; // unbiased exponent
1244 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
1246 *pfpsf
|= INVALID_EXCEPTION
;
1247 // return Integer Indefinite
1248 res
= 0x8000000000000000ull
;
1250 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
1251 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1252 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
1253 // the cases that do not fit are identified here; the ones that fit
1254 // fall through and will be handled with other cases further,
1255 // under '1 <= q + exp <= 20'
1256 if (x_sign
) { // if n < 0 and q + exp = 20 then x is much less than -1
1257 // => set invalid flag
1258 *pfpsf
|= INVALID_EXCEPTION
;
1259 // return Integer Indefinite
1260 res
= 0x8000000000000000ull
;
1262 } else { // if n > 0 and q + exp = 20
1263 // if n > 2^64 - 1 then n is too large
1264 // <=> c(0)c(1)...c(19).c(20)...c(q-1) > 2^64 - 1
1265 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 > 2^64 - 1
1266 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65 - 2)
1267 // <=> C * 10^(21-q) > 0x9fffffffffffffff6, 1<=q<=16
1269 // C * 10^20 > 0x9fffffffffffffff6
1270 __mul_128x64_to_128 (C
, C1
, ten2k128
[0]); // 10^20 * C
1271 if (C
.w
[1] > 0x09 ||
1272 (C
.w
[1] == 0x09 && C
.w
[0] > 0xfffffffffffffff6ull
)) {
1274 *pfpsf
|= INVALID_EXCEPTION
;
1275 // return Integer Indefinite
1276 res
= 0x8000000000000000ull
;
1279 // else cases that can be rounded to a 64-bit int fall through
1280 // to '1 <= q + exp <= 20'
1281 } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
1282 // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffff6
1284 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[21 - q
]);
1285 if (C
.w
[1] > 0x09 ||
1286 (C
.w
[1] == 0x09 && C
.w
[0] > 0xfffffffffffffff6ull
)) {
1288 *pfpsf
|= INVALID_EXCEPTION
;
1289 // return Integer Indefinite
1290 res
= 0x8000000000000000ull
;
1293 // else cases that can be rounded to a 64-bit int fall through
1294 // to '1 <= q + exp <= 20'
1298 // n is not too large to be converted to int64 if -1 < n < 2^64
1299 // Note: some of the cases tested for above fall through to this point
1300 if ((q
+ exp
) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1302 *pfpsf
|= INEXACT_EXCEPTION
;
1305 res
= 0x0000000000000000ull
;
1307 res
= 0x0000000000000001ull
;
1309 } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
1310 // x <= -1 or 1 <= x <= 2^64 - 1 so if positive x can be rounded
1311 // to nearest to a 64-bit unsigned signed integer
1312 if (x_sign
) { // x <= -1
1314 *pfpsf
|= INVALID_EXCEPTION
;
1315 // return Integer Indefinite
1316 res
= 0x8000000000000000ull
;
1319 // 1 <= x <= 2^64 - 1 so x can be rounded
1320 // to nearest to a 64-bit unsigned integer
1321 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
1322 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
1323 // chop off ind digits from the lower part of C1
1324 // C1 fits in 64 bits
1325 // calculate C* and f*
1326 // C* is actually floor(C*) in this case
1327 // C* and f* need shifting and masking, as shown by
1328 // shiftright128[] and maskhigh128[]
1330 // kx = 10^(-x) = ten2mk64[ind - 1]
1331 // C* = C1 * 10^(-x)
1332 // the approximation of 10^(-x) was rounded up to 54 bits
1333 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
1335 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
1336 fstar
.w
[0] = P128
.w
[0];
1337 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1338 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1339 // C* = floor(C*) (logical right shift; C has p decimal digits,
1340 // correct by Property 1)
1341 // n = C* * 10^(e+x)
1343 // shift right C* by Ex-64 = shiftright128[ind]
1344 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
1345 Cstar
= Cstar
>> shift
;
1346 // determine inexactness of the rounding of C*
1347 // if (0 < f* < 10^(-x)) then
1348 // the result is exact
1349 // else // if (f* > T*) then
1350 // the result is inexact
1351 if (ind
- 1 <= 2) { // fstar.w[1] is 0
1352 if (fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1353 // ten2mk128trunc[ind -1].w[1] is identical to
1354 // ten2mk128[ind -1].w[1]
1356 // set the inexact flag
1357 *pfpsf
|= INEXACT_EXCEPTION
;
1358 } // else the result is exact
1359 } else { // if 3 <= ind - 1 <= 14
1360 if (fstar
.w
[1] || fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1361 // ten2mk128trunc[ind -1].w[1] is identical to
1362 // ten2mk128[ind -1].w[1]
1364 // set the inexact flag
1365 *pfpsf
|= INEXACT_EXCEPTION
;
1366 } // else the result is exact
1369 res
= Cstar
; // the result is positive
1370 } else if (exp
== 0) {
1373 res
= C1
; // the result is positive
1374 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1375 // res = +C * 10^exp (exact)
1376 res
= C1
* ten2k64
[exp
]; // the result is positive
1382 /*****************************************************************************
1383 * BID64_to_uint64_int
1384 ****************************************************************************/
1386 #if DECIMAL_CALL_BY_REFERENCE
1388 bid64_to_uint64_int (UINT64
* pres
, UINT64
* px
1389 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM
)
1394 bid64_to_uint64_int (UINT64 x
1395 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM
)
1401 int exp
; // unbiased exponent
1402 // Note: C1 represents x_significand (UINT64)
1403 BID_UI64DOUBLE tmp1
;
1404 unsigned int x_nr_bits
;
1408 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
1411 // check for NaN or Infinity
1412 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
1414 *pfpsf
|= INVALID_EXCEPTION
;
1415 // return Integer Indefinite
1416 res
= 0x8000000000000000ull
;
1420 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1421 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1422 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
1423 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
1424 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
1425 if (C1
> 9999999999999999ull) { // non-canonical
1430 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
1431 C1
= x
& MASK_BINARY_SIG1
;
1434 // check for zeros (possibly from non-canonical values)
1437 res
= 0x0000000000000000ull
;
1440 // x is not special and is not zero
1442 // q = nr. of decimal digits in x (1 <= q <= 54)
1443 // determine first the nr. of bits in x
1444 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
1445 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1446 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
1447 tmp1
.d
= (double) (C1
>> 32); // exact conversion
1449 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1450 } else { // x < 2^32
1451 tmp1
.d
= (double) C1
; // exact conversion
1453 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1455 } else { // if x < 2^53
1456 tmp1
.d
= (double) C1
; // exact conversion
1458 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1460 q
= nr_digits
[x_nr_bits
- 1].digits
;
1462 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1463 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
1466 exp
= x_exp
- 398; // unbiased exponent
1468 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
1470 *pfpsf
|= INVALID_EXCEPTION
;
1471 // return Integer Indefinite
1472 res
= 0x8000000000000000ull
;
1474 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
1475 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1476 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
1477 // the cases that do not fit are identified here; the ones that fit
1478 // fall through and will be handled with other cases further,
1479 // under '1 <= q + exp <= 20'
1480 if (x_sign
) { // if n < 0 and q + exp = 20 then x is much less than -1
1481 // => set invalid flag
1482 *pfpsf
|= INVALID_EXCEPTION
;
1483 // return Integer Indefinite
1484 res
= 0x8000000000000000ull
;
1486 } else { // if n > 0 and q + exp = 20
1487 // if n >= 2^64 then n is too large
1488 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
1489 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
1490 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65)
1491 // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=16
1493 // C * 10^20 >= 0xa0000000000000000
1494 __mul_128x64_to_128 (C
, C1
, ten2k128
[0]); // 10^20 * C
1495 if (C
.w
[1] >= 0x0a) {
1496 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
1498 *pfpsf
|= INVALID_EXCEPTION
;
1499 // return Integer Indefinite
1500 res
= 0x8000000000000000ull
;
1503 // else cases that can be rounded to a 64-bit int fall through
1504 // to '1 <= q + exp <= 20'
1505 } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
1506 // Note: C * 10^(21-q) has 20 or 21 digits; 0xa0000000000000000
1508 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[21 - q
]);
1509 if (C
.w
[1] >= 0x0a) {
1510 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
1512 *pfpsf
|= INVALID_EXCEPTION
;
1513 // return Integer Indefinite
1514 res
= 0x8000000000000000ull
;
1517 // else cases that can be rounded to a 64-bit int fall through
1518 // to '1 <= q + exp <= 20'
1522 // n is not too large to be converted to int64 if -1 < n < 2^64
1523 // Note: some of the cases tested for above fall through to this point
1524 if ((q
+ exp
) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1526 res
= 0x0000000000000000ull
;
1528 } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
1529 // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
1530 // to nearest to a 64-bit unsigned signed integer
1531 if (x_sign
) { // x <= -1
1533 *pfpsf
|= INVALID_EXCEPTION
;
1534 // return Integer Indefinite
1535 res
= 0x8000000000000000ull
;
1538 // 1 <= x < 2^64 so x can be rounded
1539 // to nearest to a 64-bit unsigned integer
1540 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
1541 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
1542 // chop off ind digits from the lower part of C1
1543 // C1 fits in 64 bits
1544 // calculate C* and f*
1545 // C* is actually floor(C*) in this case
1546 // C* and f* need shifting and masking, as shown by
1547 // shiftright128[] and maskhigh128[]
1549 // kx = 10^(-x) = ten2mk64[ind - 1]
1550 // C* = C1 * 10^(-x)
1551 // the approximation of 10^(-x) was rounded up to 54 bits
1552 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
1554 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1555 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1556 // C* = floor(C*) (logical right shift; C has p decimal digits,
1557 // correct by Property 1)
1558 // n = C* * 10^(e+x)
1560 // shift right C* by Ex-64 = shiftright128[ind]
1561 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
1562 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_uint64_xint
1578 ****************************************************************************/
1580 #if DECIMAL_CALL_BY_REFERENCE
1582 bid64_to_uint64_xint (UINT64
* pres
, UINT64
* px
1583 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1588 bid64_to_uint64_xint (UINT64 x
1589 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1595 int exp
; // unbiased exponent
1596 // Note: C1 represents x_significand (UINT64)
1597 BID_UI64DOUBLE tmp1
;
1598 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
1611 res
= 0x8000000000000000ull
;
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)
1632 res
= 0x0000000000000000ull
;
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
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
1665 *pfpsf
|= INVALID_EXCEPTION
;
1666 // return Integer Indefinite
1667 res
= 0x8000000000000000ull
;
1669 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
1670 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1671 // so x rounded to an integer may or may not fit in an unsigned 64-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 <= 20'
1675 if (x_sign
) { // if n < 0 and q + exp = 20 then x is much less than -1
1676 // => set invalid flag
1677 *pfpsf
|= INVALID_EXCEPTION
;
1678 // return Integer Indefinite
1679 res
= 0x8000000000000000ull
;
1681 } else { // if n > 0 and q + exp = 20
1682 // if n >= 2^64 then n is too large
1683 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
1684 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
1685 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65)
1686 // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=16
1688 // C * 10^20 >= 0xa0000000000000000
1689 __mul_128x64_to_128 (C
, C1
, ten2k128
[0]); // 10^20 * C
1690 if (C
.w
[1] >= 0x0a) {
1691 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
1693 *pfpsf
|= INVALID_EXCEPTION
;
1694 // return Integer Indefinite
1695 res
= 0x8000000000000000ull
;
1698 // else cases that can be rounded to a 64-bit int fall through
1699 // to '1 <= q + exp <= 20'
1700 } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
1701 // Note: C * 10^(21-q) has 20 or 21 digits; 0xa0000000000000000
1703 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[21 - q
]);
1704 if (C
.w
[1] >= 0x0a) {
1705 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
1707 *pfpsf
|= INVALID_EXCEPTION
;
1708 // return Integer Indefinite
1709 res
= 0x8000000000000000ull
;
1712 // else cases that can be rounded to a 64-bit int fall through
1713 // to '1 <= q + exp <= 20'
1717 // n is not too large to be converted to int64 if -1 < n < 2^64
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
;
1723 res
= 0x0000000000000000ull
;
1725 } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
1726 // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
1727 // to nearest to a 64-bit unsigned signed integer
1728 if (x_sign
) { // x <= -1
1730 *pfpsf
|= INVALID_EXCEPTION
;
1731 // return Integer Indefinite
1732 res
= 0x8000000000000000ull
;
1735 // 1 <= x < 2^64 so x can be rounded
1736 // to nearest to a 64-bit unsigned integer
1737 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
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_uint64_rninta
1798 ****************************************************************************/
1800 #if DECIMAL_CALL_BY_REFERENCE
1802 bid64_to_uint64_rninta (UINT64
* pres
, UINT64
* px
1803 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1808 bid64_to_uint64_rninta (UINT64 x
1809 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1815 int exp
; // unbiased exponent
1816 // Note: C1 represents x_significand (UINT64)
1817 BID_UI64DOUBLE tmp1
;
1818 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
1830 res
= 0x8000000000000000ull
;
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)
1851 res
= 0x0000000000000000ull
;
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
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
1884 *pfpsf
|= INVALID_EXCEPTION
;
1885 // return Integer Indefinite
1886 res
= 0x8000000000000000ull
;
1888 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
1889 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1890 // so x rounded to an integer may or may not fit in an unsigned 64-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 <= 20'
1894 if (x_sign
) { // if n < 0 and q + exp = 20 then x is much less than -1/2
1895 // => set invalid flag
1896 *pfpsf
|= INVALID_EXCEPTION
;
1897 // return Integer Indefinite
1898 res
= 0x8000000000000000ull
;
1900 } else { // if n > 0 and q + exp = 20
1901 // if n >= 2^64 - 1/2 then n is too large
1902 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
1903 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
1904 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
1905 // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=16
1907 // C * 10^20 >= 0x9fffffffffffffffb
1908 __mul_128x64_to_128 (C
, C1
, ten2k128
[0]); // 10^20 * C
1909 if (C
.w
[1] > 0x09 ||
1910 (C
.w
[1] == 0x09 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
1912 *pfpsf
|= INVALID_EXCEPTION
;
1913 // return Integer Indefinite
1914 res
= 0x8000000000000000ull
;
1917 // else cases that can be rounded to a 64-bit int fall through
1918 // to '1 <= q + exp <= 20'
1919 } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
1920 // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffffb
1922 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[21 - q
]);
1923 if (C
.w
[1] > 0x09 ||
1924 (C
.w
[1] == 0x09 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
1926 *pfpsf
|= INVALID_EXCEPTION
;
1927 // return Integer Indefinite
1928 res
= 0x8000000000000000ull
;
1931 // else cases that can be rounded to a 64-bit int fall through
1932 // to '1 <= q + exp <= 20'
1936 // n is not too large to be converted to int64 if -1/2 <= n < 2^64 - 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)
1940 res
= 0x0000000000000000ull
;
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)
1949 ind
= q
- 1; // 0 <= ind <= 15
1950 if (C1
< midpoint64
[ind
]) {
1951 res
= 0x0000000000000000ull
; // return 0
1952 } else if (!x_sign
) { // n > 0
1953 res
= 0x0000000000000001ull
; // return +1
1954 } else { // if n < 0
1955 res
= 0x8000000000000000ull
;
1956 *pfpsf
|= INVALID_EXCEPTION
;
1959 } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
1960 // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
1961 // to nearest to a 64-bit unsigned signed integer
1962 if (x_sign
) { // x <= -1
1964 *pfpsf
|= INVALID_EXCEPTION
;
1965 // return Integer Indefinite
1966 res
= 0x8000000000000000ull
;
1969 // 1 <= x < 2^64-1/2 so x can be rounded
1970 // to nearest to a 64-bit unsigned integer
1971 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
1972 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
1973 // chop off ind digits from the lower part of C1
1974 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
1975 C1
= C1
+ midpoint64
[ind
- 1];
1976 // calculate C* and f*
1977 // C* is actually floor(C*) in this case
1978 // C* and f* need shifting and masking, as shown by
1979 // shiftright128[] and maskhigh128[]
1981 // kx = 10^(-x) = ten2mk64[ind - 1]
1982 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1983 // the approximation of 10^(-x) was rounded up to 54 bits
1984 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
1986 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1987 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1988 // if (0 < f* < 10^(-x)) then the result is a midpoint
1989 // if floor(C*) is even then C* = floor(C*) - logical right
1990 // shift; C* has p decimal digits, correct by Prop. 1)
1991 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
1992 // shift; C* has p decimal digits, correct by Pr. 1)
1994 // C* = floor(C*) (logical right shift; C has p decimal digits,
1995 // correct by Property 1)
1996 // n = C* * 10^(e+x)
1998 // shift right C* by Ex-64 = shiftright128[ind]
1999 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
2000 Cstar
= Cstar
>> shift
;
2002 // if the result was a midpoint it was rounded away from zero
2003 res
= Cstar
; // the result is positive
2004 } else if (exp
== 0) {
2007 res
= C1
; // the result is positive
2008 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2009 // res = +C * 10^exp (exact)
2010 res
= C1
* ten2k64
[exp
]; // the result is positive
2016 /*****************************************************************************
2017 * BID64_to_uint64_xrninta
2018 ****************************************************************************/
2020 #if DECIMAL_CALL_BY_REFERENCE
2022 bid64_to_uint64_xrninta (UINT64
* pres
, UINT64
* px
2023 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
2028 bid64_to_uint64_xrninta (UINT64 x
2029 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
2035 int exp
; // unbiased exponent
2036 // Note: C1 represents x_significand (UINT64)
2038 BID_UI64DOUBLE tmp1
;
2039 unsigned int x_nr_bits
;
2043 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
2047 // check for NaN or Infinity
2048 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
2050 *pfpsf
|= INVALID_EXCEPTION
;
2051 // return Integer Indefinite
2052 res
= 0x8000000000000000ull
;
2056 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
2057 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
2058 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
2059 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
2060 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
2061 if (C1
> 9999999999999999ull) { // non-canonical
2066 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
2067 C1
= x
& MASK_BINARY_SIG1
;
2070 // check for zeros (possibly from non-canonical values)
2073 res
= 0x0000000000000000ull
;
2076 // x is not special and is not zero
2078 // q = nr. of decimal digits in x (1 <= q <= 54)
2079 // determine first the nr. of bits in x
2080 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
2081 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2082 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
2083 tmp1
.d
= (double) (C1
>> 32); // exact conversion
2085 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2086 } else { // x < 2^32
2087 tmp1
.d
= (double) C1
; // exact conversion
2089 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2091 } else { // if x < 2^53
2092 tmp1
.d
= (double) C1
; // exact conversion
2094 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2096 q
= nr_digits
[x_nr_bits
- 1].digits
;
2098 q
= nr_digits
[x_nr_bits
- 1].digits1
;
2099 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
2102 exp
= x_exp
- 398; // unbiased exponent
2104 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
2106 *pfpsf
|= INVALID_EXCEPTION
;
2107 // return Integer Indefinite
2108 res
= 0x8000000000000000ull
;
2110 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
2111 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
2112 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
2113 // the cases that do not fit are identified here; the ones that fit
2114 // fall through and will be handled with other cases further,
2115 // under '1 <= q + exp <= 20'
2116 if (x_sign
) { // if n < 0 and q + exp = 20 then x is much less than -1/2
2117 // => set invalid flag
2118 *pfpsf
|= INVALID_EXCEPTION
;
2119 // return Integer Indefinite
2120 res
= 0x8000000000000000ull
;
2122 } else { // if n > 0 and q + exp = 20
2123 // if n >= 2^64 - 1/2 then n is too large
2124 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
2125 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
2126 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
2127 // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=16
2129 // C * 10^20 >= 0x9fffffffffffffffb
2130 __mul_128x64_to_128 (C
, C1
, ten2k128
[0]); // 10^20 * C
2131 if (C
.w
[1] > 0x09 ||
2132 (C
.w
[1] == 0x09 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
2134 *pfpsf
|= INVALID_EXCEPTION
;
2135 // return Integer Indefinite
2136 res
= 0x8000000000000000ull
;
2139 // else cases that can be rounded to a 64-bit int fall through
2140 // to '1 <= q + exp <= 20'
2141 } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
2142 // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffffb
2144 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[21 - q
]);
2145 if (C
.w
[1] > 0x09 ||
2146 (C
.w
[1] == 0x09 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
2148 *pfpsf
|= INVALID_EXCEPTION
;
2149 // return Integer Indefinite
2150 res
= 0x8000000000000000ull
;
2153 // else cases that can be rounded to a 64-bit int fall through
2154 // to '1 <= q + exp <= 20'
2158 // n is not too large to be converted to int64 if -1/2 <= n < 2^64 - 1/2
2159 // Note: some of the cases tested for above fall through to this point
2160 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
2162 *pfpsf
|= INEXACT_EXCEPTION
;
2164 res
= 0x0000000000000000ull
;
2166 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
2167 // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
2173 ind
= q
- 1; // 0 <= ind <= 15
2174 if (C1
< midpoint64
[ind
]) {
2175 res
= 0x0000000000000000ull
; // return 0
2176 } else if (!x_sign
) { // n > 0
2177 res
= 0x0000000000000001ull
; // return +1
2178 } else { // if n < 0
2179 res
= 0x8000000000000000ull
;
2180 *pfpsf
|= INVALID_EXCEPTION
;
2184 *pfpsf
|= INEXACT_EXCEPTION
;
2185 } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
2186 // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
2187 // to nearest to a 64-bit unsigned signed integer
2188 if (x_sign
) { // x <= -1
2190 *pfpsf
|= INVALID_EXCEPTION
;
2191 // return Integer Indefinite
2192 res
= 0x8000000000000000ull
;
2195 // 1 <= x < 2^64-1/2 so x can be rounded
2196 // to nearest to a 64-bit unsigned integer
2197 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
2198 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
2199 // chop off ind digits from the lower part of C1
2200 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
2201 C1
= C1
+ midpoint64
[ind
- 1];
2202 // calculate C* and f*
2203 // C* is actually floor(C*) in this case
2204 // C* and f* need shifting and masking, as shown by
2205 // shiftright128[] and maskhigh128[]
2207 // kx = 10^(-x) = ten2mk64[ind - 1]
2208 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2209 // the approximation of 10^(-x) was rounded up to 54 bits
2210 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
2212 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
2213 fstar
.w
[0] = P128
.w
[0];
2214 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
2215 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
2216 // if (0 < f* < 10^(-x)) then the result is a midpoint
2217 // if floor(C*) is even then C* = floor(C*) - logical right
2218 // shift; C* has p decimal digits, correct by Prop. 1)
2219 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
2220 // shift; C* has p decimal digits, correct by Pr. 1)
2222 // C* = floor(C*) (logical right shift; C has p decimal digits,
2223 // correct by Property 1)
2224 // n = C* * 10^(e+x)
2226 // shift right C* by Ex-64 = shiftright128[ind]
2227 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
2228 Cstar
= Cstar
>> shift
;
2229 // determine inexactness of the rounding of C*
2230 // if (0 < f* - 1/2 < 10^(-x)) then
2231 // the result is exact
2232 // else // if (f* - 1/2 > T*) then
2233 // the result is inexact
2234 if (ind
- 1 <= 2) { // fstar.w[1] is 0
2235 if (fstar
.w
[0] > 0x8000000000000000ull
) {
2236 // f* > 1/2 and the result may be exact
2237 tmp64
= fstar
.w
[0] - 0x8000000000000000ull
; // f* - 1/2
2238 if ((tmp64
> ten2mk128trunc
[ind
- 1].w
[1])) {
2239 // ten2mk128trunc[ind -1].w[1] is identical to
2240 // ten2mk128[ind -1].w[1]
2241 // set the inexact flag
2242 *pfpsf
|= INEXACT_EXCEPTION
;
2243 } // else the result is exact
2244 } else { // the result is inexact; f2* <= 1/2
2245 // set the inexact flag
2246 *pfpsf
|= INEXACT_EXCEPTION
;
2248 } else { // if 3 <= ind - 1 <= 14
2249 if (fstar
.w
[1] > onehalf128
[ind
- 1] ||
2250 (fstar
.w
[1] == onehalf128
[ind
- 1] && fstar
.w
[0])) {
2251 // f2* > 1/2 and the result may be exact
2252 // Calculate f2* - 1/2
2253 tmp64
= fstar
.w
[1] - onehalf128
[ind
- 1];
2254 if (tmp64
|| fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
2255 // ten2mk128trunc[ind -1].w[1] is identical to
2256 // ten2mk128[ind -1].w[1]
2257 // set the inexact flag
2258 *pfpsf
|= INEXACT_EXCEPTION
;
2259 } // else the result is exact
2260 } else { // the result is inexact; f2* <= 1/2
2261 // set the inexact flag
2262 *pfpsf
|= INEXACT_EXCEPTION
;
2266 // if the result was a midpoint it was rounded away from zero
2267 res
= Cstar
; // the result is positive
2268 } else if (exp
== 0) {
2271 res
= C1
; // the result is positive
2272 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2273 // res = +C * 10^exp (exact)
2274 res
= C1
* ten2k64
[exp
]; // the result is positive