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_int64_rnint
33 ****************************************************************************/
35 #if DECIMAL_CALL_BY_REFERENCE
37 bid64_to_int64_rnint (SINT64
* pres
, UINT64
* px
38 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
43 bid64_to_int64_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)
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
) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
120 *pfpsf
|= INVALID_EXCEPTION
;
121 // return Integer Indefinite
122 res
= 0x8000000000000000ull
;
124 } else if ((q
+ exp
) == 19) { // x = c(0)c(1)...c(18).c(19)...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 a signed 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 <= 19'
130 if (x_sign
) { // if n < 0 and q + exp = 19
131 // if n < -2^63 - 1/2 then n is too large
132 // too large if c(0)c(1)...c(18).c(19)...c(q-1) > 2^63+1/2
133 // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 5*(2^64+1), 1<=q<=16
134 // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 0x50000000000000005, 1<=q<=16
135 // <=> C * 10^(20-q) > 0x50000000000000005, 1<=q<=16
136 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
137 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
138 // Note: C1 * 10^(11-q) has 19 or 20 digits; 0x50000000000000005, has 20
139 if (C
.w
[1] > 0x05ull
|| (C
.w
[1] == 0x05ull
&& C
.w
[0] > 0x05ull
)) {
141 *pfpsf
|= INVALID_EXCEPTION
;
142 // return Integer Indefinite
143 res
= 0x8000000000000000ull
;
146 // else cases that can be rounded to a 64-bit int fall through
147 // to '1 <= q + exp <= 19'
148 } else { // if n > 0 and q + exp = 19
149 // if n >= 2^63 - 1/2 then n is too large
150 // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63-1/2
151 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64-1), 1<=q<=16
152 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x4fffffffffffffffb, 1<=q<=16
153 // <=> if C * 10^(20-q) >= 0x4fffffffffffffffb, 1<=q<=16
154 C
.w
[1] = 0x0000000000000004ull
;
155 C
.w
[0] = 0xfffffffffffffffbull
;
156 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
157 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
158 if (C
.w
[1] > 0x04ull
||
159 (C
.w
[1] == 0x04ull
&& C
.w
[0] >= 0xfffffffffffffffbull
)) {
161 *pfpsf
|= INVALID_EXCEPTION
;
162 // return Integer Indefinite
163 res
= 0x8000000000000000ull
;
166 // else cases that can be rounded to a 64-bit int fall through
167 // to '1 <= q + exp <= 19'
168 } // end else if n > 0 and q + exp = 19
169 } // end else if ((q + exp) == 19)
171 // n is not too large to be converted to int64: -2^63-1/2 <= n < 2^63-1/2
172 // Note: some of the cases tested for above fall through to this point
173 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
175 res
= 0x0000000000000000ull
;
177 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
178 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
182 ind
= q
- 1; // 0 <= ind <= 15
183 if (C1
<= midpoint64
[ind
]) {
184 res
= 0x0000000000000000ull
; // return 0
185 } else if (x_sign
) { // n < 0
186 res
= 0xffffffffffffffffull
; // return -1
188 res
= 0x0000000000000001ull
; // return +1
190 } else { // if (1 <= q + exp <= 19, 1 <= q <= 16, -15 <= exp <= 18)
191 // -2^63-1/2 <= x <= -1 or 1 <= x < 2^63-1/2 so x can be rounded
192 // to nearest to a 64-bit signed integer
193 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 19
194 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
195 // chop off ind digits from the lower part of C1
196 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
197 C1
= C1
+ midpoint64
[ind
- 1];
198 // calculate C* and f*
199 // C* is actually floor(C*) in this case
200 // C* and f* need shifting and masking, as shown by
201 // shiftright128[] and maskhigh128[]
203 // kx = 10^(-x) = ten2mk64[ind - 1]
204 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
205 // the approximation of 10^(-x) was rounded up to 54 bits
206 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
208 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
209 fstar
.w
[0] = P128
.w
[0];
210 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
211 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
212 // if (0 < f* < 10^(-x)) then the result is a midpoint
213 // if floor(C*) is even then C* = floor(C*) - logical right
214 // shift; C* has p decimal digits, correct by Prop. 1)
215 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
216 // shift; C* has p decimal digits, correct by Pr. 1)
218 // C* = floor(C*) (logical right shift; C has p decimal digits,
219 // correct by Property 1)
222 // shift right C* by Ex-64 = shiftright128[ind]
223 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
224 Cstar
= Cstar
>> shift
;
226 // if the result was a midpoint it was rounded away from zero, so
227 // it will need a correction
228 // check for midpoints
229 if ((fstar
.w
[1] == 0) && fstar
.w
[0] &&
230 (fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[1])) {
231 // ten2mk128trunc[ind -1].w[1] is identical to
232 // ten2mk128[ind -1].w[1]
233 // the result is a midpoint; round to nearest
234 if (Cstar
& 0x01) { // Cstar is odd; MP in [EVEN, ODD]
235 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
236 Cstar
--; // Cstar is now even
237 } // else MP in [ODD, EVEN]
243 } else if (exp
== 0) {
245 // res = +/-C (exact)
250 } else { // if (exp > 0) => 1 <= exp <= 18, 1 <= q <= 16, 2 <= q + exp <= 20
251 // (the upper limit of 20 on q + exp is due to the fact that
252 // +/-C * 10^exp is guaranteed to fit in 64 bits)
253 // res = +/-C * 10^exp (exact)
255 res
= -C1
* ten2k64
[exp
];
257 res
= C1
* ten2k64
[exp
];
263 /*****************************************************************************
264 * BID64_to_int64_xrnint
265 ****************************************************************************/
267 #if DECIMAL_CALL_BY_REFERENCE
269 bid64_to_int64_xrnint (SINT64
* pres
, UINT64
* px
270 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
275 bid64_to_int64_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
;
290 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
294 // check for NaN or Infinity
295 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
297 *pfpsf
|= INVALID_EXCEPTION
;
298 // return Integer Indefinite
299 res
= 0x8000000000000000ull
;
303 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
304 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
305 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
306 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
307 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
308 if (C1
> 9999999999999999ull) { // non-canonical
313 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
314 C1
= x
& MASK_BINARY_SIG1
;
317 // check for zeros (possibly from non-canonical values)
323 // x is not special and is not zero
325 // q = nr. of decimal digits in x (1 <= q <= 54)
326 // determine first the nr. of bits in x
327 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
328 // split the 64-bit value in two 32-bit halves to avoid rounding errors
329 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
330 tmp1
.d
= (double) (C1
>> 32); // exact conversion
332 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
334 tmp1
.d
= (double) C1
; // exact conversion
336 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
338 } else { // if x < 2^53
339 tmp1
.d
= (double) C1
; // exact conversion
341 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
343 q
= nr_digits
[x_nr_bits
- 1].digits
;
345 q
= nr_digits
[x_nr_bits
- 1].digits1
;
346 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
349 exp
= x_exp
- 398; // unbiased exponent
351 if ((q
+ exp
) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
353 *pfpsf
|= INVALID_EXCEPTION
;
354 // return Integer Indefinite
355 res
= 0x8000000000000000ull
;
357 } else if ((q
+ exp
) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
358 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
359 // so x rounded to an integer may or may not fit in a signed 64-bit int
360 // the cases that do not fit are identified here; the ones that fit
361 // fall through and will be handled with other cases further,
362 // under '1 <= q + exp <= 19'
363 if (x_sign
) { // if n < 0 and q + exp = 19
364 // if n < -2^63 - 1/2 then n is too large
365 // too large if c(0)c(1)...c(18).c(19)...c(q-1) > 2^63+1/2
366 // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 5*(2^64+1), 1<=q<=16
367 // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 0x50000000000000005, 1<=q<=16
368 // <=> C * 10^(20-q) > 0x50000000000000005, 1<=q<=16
369 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
370 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
371 // Note: C1 * 10^(11-q) has 19 or 20 digits; 0x50000000000000005, has 20
372 if (C
.w
[1] > 0x05ull
|| (C
.w
[1] == 0x05ull
&& C
.w
[0] > 0x05ull
)) {
374 *pfpsf
|= INVALID_EXCEPTION
;
375 // return Integer Indefinite
376 res
= 0x8000000000000000ull
;
379 // else cases that can be rounded to a 64-bit int fall through
380 // to '1 <= q + exp <= 19'
381 } else { // if n > 0 and q + exp = 19
382 // if n >= 2^63 - 1/2 then n is too large
383 // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63-1/2
384 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64-1), 1<=q<=16
385 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x4fffffffffffffffb, 1<=q<=16
386 // <=> if C * 10^(20-q) >= 0x4fffffffffffffffb, 1<=q<=16
387 C
.w
[1] = 0x0000000000000004ull
;
388 C
.w
[0] = 0xfffffffffffffffbull
;
389 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
390 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
391 if (C
.w
[1] > 0x04ull
||
392 (C
.w
[1] == 0x04ull
&& C
.w
[0] >= 0xfffffffffffffffbull
)) {
394 *pfpsf
|= INVALID_EXCEPTION
;
395 // return Integer Indefinite
396 res
= 0x8000000000000000ull
;
399 // else cases that can be rounded to a 64-bit int fall through
400 // to '1 <= q + exp <= 19'
401 } // end else if n > 0 and q + exp = 19
402 } // end else if ((q + exp) == 19)
404 // n is not too large to be converted to int64: -2^63-1/2 <= n < 2^63-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
;
410 res
= 0x0000000000000000ull
;
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)
417 ind
= q
- 1; // 0 <= ind <= 15
418 if (C1
<= midpoint64
[ind
]) {
419 res
= 0x0000000000000000ull
; // return 0
420 } else if (x_sign
) { // n < 0
421 res
= 0xffffffffffffffffull
; // return -1
423 res
= 0x0000000000000001ull
; // return +1
426 *pfpsf
|= INEXACT_EXCEPTION
;
427 } else { // if (1 <= q + exp <= 19, 1 <= q <= 16, -15 <= exp <= 18)
428 // -2^63-1/2 <= x <= -1 or 1 <= x < 2^63-1/2 so x can be rounded
429 // to nearest to a 64-bit signed integer
430 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 19
431 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
432 // chop off ind digits from the lower part of C1
433 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
434 C1
= C1
+ midpoint64
[ind
- 1];
435 // calculate C* and f*
436 // C* is actually floor(C*) in this case
437 // C* and f* need shifting and masking, as shown by
438 // shiftright128[] and maskhigh128[]
440 // kx = 10^(-x) = ten2mk64[ind - 1]
441 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
442 // the approximation of 10^(-x) was rounded up to 54 bits
443 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
445 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
446 fstar
.w
[0] = P128
.w
[0];
447 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
448 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
449 // if (0 < f* < 10^(-x)) then the result is a midpoint
450 // if floor(C*) is even then C* = floor(C*) - logical right
451 // shift; C* has p decimal digits, correct by Prop. 1)
452 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
453 // shift; C* has p decimal digits, correct by Pr. 1)
455 // C* = floor(C*) (logical right shift; C has p decimal digits,
456 // correct by Property 1)
459 // shift right C* by Ex-64 = shiftright128[ind]
460 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
461 Cstar
= Cstar
>> shift
;
462 // determine inexactness of the rounding of C*
463 // if (0 < f* - 1/2 < 10^(-x)) then
464 // the result is exact
465 // else // if (f* - 1/2 > T*) then
466 // the result is inexact
468 if (fstar
.w
[0] > 0x8000000000000000ull
) {
469 // f* > 1/2 and the result may be exact
470 tmp64
= fstar
.w
[0] - 0x8000000000000000ull
; // f* - 1/2
471 if ((tmp64
> ten2mk128trunc
[ind
- 1].w
[1])) {
472 // ten2mk128trunc[ind -1].w[1] is identical to
473 // ten2mk128[ind -1].w[1]
474 // set the inexact flag
475 *pfpsf
|= INEXACT_EXCEPTION
;
476 } // else the result is exact
477 } else { // the result is inexact; f2* <= 1/2
478 // set the inexact flag
479 *pfpsf
|= INEXACT_EXCEPTION
;
481 } else { // if 3 <= ind - 1 <= 14
482 if (fstar
.w
[1] > onehalf128
[ind
- 1] ||
483 (fstar
.w
[1] == onehalf128
[ind
- 1] && fstar
.w
[0])) {
484 // f2* > 1/2 and the result may be exact
485 // Calculate f2* - 1/2
486 tmp64
= fstar
.w
[1] - onehalf128
[ind
- 1];
487 if (tmp64
|| fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
488 // ten2mk128trunc[ind -1].w[1] is identical to
489 // ten2mk128[ind -1].w[1]
490 // set the inexact flag
491 *pfpsf
|= INEXACT_EXCEPTION
;
492 } // else the result is exact
493 } else { // the result is inexact; f2* <= 1/2
494 // set the inexact flag
495 *pfpsf
|= INEXACT_EXCEPTION
;
499 // if the result was a midpoint it was rounded away from zero, so
500 // it will need a correction
501 // check for midpoints
502 if ((fstar
.w
[1] == 0) && fstar
.w
[0] &&
503 (fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[1])) {
504 // ten2mk128trunc[ind -1].w[1] is identical to
505 // ten2mk128[ind -1].w[1]
506 // the result is a midpoint; round to nearest
507 if (Cstar
& 0x01) { // Cstar is odd; MP in [EVEN, ODD]
508 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
509 Cstar
--; // Cstar is now even
510 } // else MP in [ODD, EVEN]
516 } else if (exp
== 0) {
518 // res = +/-C (exact)
523 } else { // if (exp > 0) => 1 <= exp <= 18, 1 <= q <= 16, 2 <= q + exp <= 20
524 // (the upper limit of 20 on q + exp is due to the fact that
525 // +/-C * 10^exp is guaranteed to fit in 64 bits)
526 // res = +/-C * 10^exp (exact)
528 res
= -C1
* ten2k64
[exp
];
530 res
= C1
* ten2k64
[exp
];
536 /*****************************************************************************
537 * BID64_to_int64_floor
538 ****************************************************************************/
540 #if DECIMAL_CALL_BY_REFERENCE
542 bid64_to_int64_floor (SINT64
* pres
, UINT64
* px
543 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
548 bid64_to_int64_floor (UINT64 x
549 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
555 int exp
; // unbiased exponent
556 // Note: C1 represents x_significand (UINT64)
558 unsigned int x_nr_bits
;
562 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
566 // check for NaN or Infinity
567 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
569 *pfpsf
|= INVALID_EXCEPTION
;
570 // return Integer Indefinite
571 res
= 0x8000000000000000ull
;
575 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
576 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
577 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
578 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
579 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
580 if (C1
> 9999999999999999ull) { // non-canonical
585 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
586 C1
= x
& MASK_BINARY_SIG1
;
589 // check for zeros (possibly from non-canonical values)
592 res
= 0x0000000000000000ull
;
595 // x is not special and is not zero
597 // q = nr. of decimal digits in x (1 <= q <= 54)
598 // determine first the nr. of bits in x
599 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
600 // split the 64-bit value in two 32-bit halves to avoid rounding errors
601 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
602 tmp1
.d
= (double) (C1
>> 32); // exact conversion
604 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
606 tmp1
.d
= (double) C1
; // exact conversion
608 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
610 } else { // if x < 2^53
611 tmp1
.d
= (double) C1
; // exact conversion
613 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
615 q
= nr_digits
[x_nr_bits
- 1].digits
;
617 q
= nr_digits
[x_nr_bits
- 1].digits1
;
618 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
621 exp
= x_exp
- 398; // unbiased exponent
623 if ((q
+ exp
) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
625 *pfpsf
|= INVALID_EXCEPTION
;
626 // return Integer Indefinite
627 res
= 0x8000000000000000ull
;
629 } else if ((q
+ exp
) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
630 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
631 // so x rounded to an integer may or may not fit in a signed 64-bit int
632 // the cases that do not fit are identified here; the ones that fit
633 // fall through and will be handled with other cases further,
634 // under '1 <= q + exp <= 19'
635 if (x_sign
) { // if n < 0 and q + exp = 19
636 // if n < -2^63 then n is too large
637 // too large if c(0)c(1)...c(18).c(19)...c(q-1) > 2^63
638 // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 5*2^64, 1<=q<=16
639 // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 0x50000000000000000, 1<=q<=16
640 // <=> C * 10^(20-q) > 0x50000000000000000, 1<=q<=16
641 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
642 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
643 // Note: C1 * 10^(11-q) has 19 or 20 digits; 0x5000000000000000a, has 20
644 if (C
.w
[1] > 0x05ull
|| (C
.w
[1] == 0x05ull
&& C
.w
[0] != 0)) {
646 *pfpsf
|= INVALID_EXCEPTION
;
647 // return Integer Indefinite
648 res
= 0x8000000000000000ull
;
651 // else cases that can be rounded to a 64-bit int fall through
652 // to '1 <= q + exp <= 19'
653 } else { // if n > 0 and q + exp = 19
654 // if n >= 2^63 then n is too large
655 // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63
656 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*2^64, 1<=q<=16
657 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x50000000000000000, 1<=q<=16
658 // <=> if C * 10^(20-q) >= 0x50000000000000000, 1<=q<=16
659 C
.w
[1] = 0x0000000000000005ull
;
660 C
.w
[0] = 0x0000000000000000ull
;
661 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
662 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
663 if (C
.w
[1] >= 0x05ull
) {
664 // actually C.w[1] == 0x05ull && C.w[0] >= 0x0000000000000000ull) {
666 *pfpsf
|= INVALID_EXCEPTION
;
667 // return Integer Indefinite
668 res
= 0x8000000000000000ull
;
671 // else cases that can be rounded to a 64-bit int fall through
672 // to '1 <= q + exp <= 19'
673 } // end else if n > 0 and q + exp = 19
674 } // end else if ((q + exp) == 19)
676 // n is not too large to be converted to int64: -2^63 <= n < 2^63
677 // Note: some of the cases tested for above fall through to this point
678 if ((q
+ exp
) <= 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
681 res
= 0xffffffffffffffffull
;
683 res
= 0x0000000000000000ull
;
685 } else { // if (1 <= q + exp <= 19, 1 <= q <= 16, -15 <= exp <= 18)
686 // -2^63 <= x <= -1 or 1 <= x < 2^63 so x can be rounded
687 // to nearest to a 64-bit signed integer
688 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 19
689 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
690 // chop off ind digits from the lower part of C1
691 // C1 fits in 64 bits
692 // calculate C* and f*
693 // C* is actually floor(C*) in this case
694 // C* and f* need shifting and masking, as shown by
695 // shiftright128[] and maskhigh128[]
697 // kx = 10^(-x) = ten2mk64[ind - 1]
699 // the approximation of 10^(-x) was rounded up to 54 bits
700 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
702 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
703 fstar
.w
[0] = P128
.w
[0];
704 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
705 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
706 // C* = floor(C*) (logical right shift; C has p decimal digits,
707 // correct by Property 1)
710 // shift right C* by Ex-64 = shiftright128[ind]
711 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
712 Cstar
= Cstar
>> shift
;
713 // determine inexactness of the rounding of C*
714 // if (0 < f* < 10^(-x)) then
715 // the result is exact
716 // else // if (f* > T*) then
717 // the result is inexact
718 if (ind
- 1 <= 2) { // fstar.w[1] is 0
719 if (fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
720 // ten2mk128trunc[ind -1].w[1] is identical to
721 // ten2mk128[ind -1].w[1]
722 if (x_sign
) { // negative and inexact
725 } // else the result is exact
726 } else { // if 3 <= ind - 1 <= 14
727 if (fstar
.w
[1] || fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
728 // ten2mk128trunc[ind -1].w[1] is identical to
729 // ten2mk128[ind -1].w[1]
730 if (x_sign
) { // negative and inexact
733 } // else the result is exact
740 } else if (exp
== 0) {
742 // res = +/-C (exact)
747 } else { // if (exp > 0) => 1 <= exp <= 18, 1 <= q <= 16, 2 <= q + exp <= 20
748 // (the upper limit of 20 on q + exp is due to the fact that
749 // +/-C * 10^exp is guaranteed to fit in 64 bits)
750 // res = +/-C * 10^exp (exact)
752 res
= -C1
* ten2k64
[exp
];
754 res
= C1
* ten2k64
[exp
];
760 /*****************************************************************************
761 * BID64_to_int64_xfloor
762 ****************************************************************************/
764 #if DECIMAL_CALL_BY_REFERENCE
766 bid64_to_int64_xfloor (SINT64
* pres
, UINT64
* px
767 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
772 bid64_to_int64_xfloor (UINT64 x
773 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
779 int exp
; // unbiased exponent
780 // Note: C1 represents x_significand (UINT64)
782 unsigned int x_nr_bits
;
786 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
790 // check for NaN or Infinity
791 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
793 *pfpsf
|= INVALID_EXCEPTION
;
794 // return Integer Indefinite
795 res
= 0x8000000000000000ull
;
799 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
800 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
801 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
802 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
803 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
804 if (C1
> 9999999999999999ull) { // non-canonical
809 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
810 C1
= x
& MASK_BINARY_SIG1
;
813 // check for zeros (possibly from non-canonical values)
816 res
= 0x0000000000000000ull
;
819 // x is not special and is not zero
821 // q = nr. of decimal digits in x (1 <= q <= 54)
822 // determine first the nr. of bits in x
823 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
824 // split the 64-bit value in two 32-bit halves to avoid rounding errors
825 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
826 tmp1
.d
= (double) (C1
>> 32); // exact conversion
828 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
830 tmp1
.d
= (double) C1
; // exact conversion
832 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
834 } else { // if x < 2^53
835 tmp1
.d
= (double) C1
; // exact conversion
837 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
839 q
= nr_digits
[x_nr_bits
- 1].digits
;
841 q
= nr_digits
[x_nr_bits
- 1].digits1
;
842 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
845 exp
= x_exp
- 398; // unbiased exponent
847 if ((q
+ exp
) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
849 *pfpsf
|= INVALID_EXCEPTION
;
850 // return Integer Indefinite
851 res
= 0x8000000000000000ull
;
853 } else if ((q
+ exp
) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
854 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
855 // so x rounded to an integer may or may not fit in a signed 64-bit int
856 // the cases that do not fit are identified here; the ones that fit
857 // fall through and will be handled with other cases further,
858 // under '1 <= q + exp <= 19'
859 if (x_sign
) { // if n < 0 and q + exp = 19
860 // if n < -2^63 then n is too large
861 // too large if c(0)c(1)...c(18).c(19)...c(q-1) > 2^63
862 // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 5*2^64, 1<=q<=16
863 // <=> 0.c(0)c(1)...c(q-1) * 10^20 > 0x50000000000000000, 1<=q<=16
864 // <=> C * 10^(20-q) > 0x50000000000000000, 1<=q<=16
865 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
866 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
867 // Note: C1 * 10^(11-q) has 19 or 20 digits; 0x5000000000000000a, has 20
868 if (C
.w
[1] > 0x05ull
|| (C
.w
[1] == 0x05ull
&& C
.w
[0] != 0)) {
870 *pfpsf
|= INVALID_EXCEPTION
;
871 // return Integer Indefinite
872 res
= 0x8000000000000000ull
;
875 // else cases that can be rounded to a 64-bit int fall through
876 // to '1 <= q + exp <= 19'
877 } else { // if n > 0 and q + exp = 19
878 // if n >= 2^63 then n is too large
879 // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63
880 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*2^64, 1<=q<=16
881 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x50000000000000000, 1<=q<=16
882 // <=> if C * 10^(20-q) >= 0x50000000000000000, 1<=q<=16
883 C
.w
[1] = 0x0000000000000005ull
;
884 C
.w
[0] = 0x0000000000000000ull
;
885 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
886 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
887 if (C
.w
[1] >= 0x05ull
) {
888 // actually C.w[1] == 0x05ull && C.w[0] >= 0x0000000000000000ull) {
890 *pfpsf
|= INVALID_EXCEPTION
;
891 // return Integer Indefinite
892 res
= 0x8000000000000000ull
;
895 // else cases that can be rounded to a 64-bit int fall through
896 // to '1 <= q + exp <= 19'
897 } // end else if n > 0 and q + exp = 19
898 } // end else if ((q + exp) == 19)
900 // n is not too large to be converted to int64: -2^63 <= n < 2^63
901 // Note: some of the cases tested for above fall through to this point
902 if ((q
+ exp
) <= 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
904 *pfpsf
|= INEXACT_EXCEPTION
;
907 res
= 0xffffffffffffffffull
;
909 res
= 0x0000000000000000ull
;
911 } else { // if (1 <= q + exp <= 19, 1 <= q <= 16, -15 <= exp <= 18)
912 // -2^63 <= x <= -1 or 1 <= x < 2^63 so x can be rounded
913 // to nearest to a 64-bit signed integer
914 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 19
915 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
916 // chop off ind digits from the lower part of C1
917 // C1 fits in 64 bits
918 // calculate C* and f*
919 // C* is actually floor(C*) in this case
920 // C* and f* need shifting and masking, as shown by
921 // shiftright128[] and maskhigh128[]
923 // kx = 10^(-x) = ten2mk64[ind - 1]
925 // the approximation of 10^(-x) was rounded up to 54 bits
926 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
928 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
929 fstar
.w
[0] = P128
.w
[0];
930 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
931 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
932 // C* = floor(C*) (logical right shift; C has p decimal digits,
933 // correct by Property 1)
936 // shift right C* by Ex-64 = shiftright128[ind]
937 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
938 Cstar
= Cstar
>> shift
;
939 // determine inexactness of the rounding of C*
940 // if (0 < f* < 10^(-x)) then
941 // the result is exact
942 // else // if (f* > T*) then
943 // the result is inexact
944 if (ind
- 1 <= 2) { // fstar.w[1] is 0
945 if (fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
946 // ten2mk128trunc[ind -1].w[1] is identical to
947 // ten2mk128[ind -1].w[1]
948 if (x_sign
) { // negative and inexact
951 // set the inexact flag
952 *pfpsf
|= INEXACT_EXCEPTION
;
953 } // else the result is exact
954 } else { // if 3 <= ind - 1 <= 14
955 if (fstar
.w
[1] || fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
956 // ten2mk128trunc[ind -1].w[1] is identical to
957 // ten2mk128[ind -1].w[1]
958 if (x_sign
) { // negative and inexact
961 // set the inexact flag
962 *pfpsf
|= INEXACT_EXCEPTION
;
963 } // else the result is exact
970 } else if (exp
== 0) {
972 // res = +/-C (exact)
977 } else { // if (exp > 0) => 1 <= exp <= 18, 1 <= q <= 16, 2 <= q + exp <= 20
978 // (the upper limit of 20 on q + exp is due to the fact that
979 // +/-C * 10^exp is guaranteed to fit in 64 bits)
980 // res = +/-C * 10^exp (exact)
982 res
= -C1
* ten2k64
[exp
];
984 res
= C1
* ten2k64
[exp
];
990 /*****************************************************************************
991 * BID64_to_int64_ceil
992 ****************************************************************************/
994 #if DECIMAL_CALL_BY_REFERENCE
996 bid64_to_int64_ceil (SINT64
* pres
, UINT64
* px
997 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM
)
1002 bid64_to_int64_ceil (UINT64 x
1003 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM
)
1009 int exp
; // unbiased exponent
1010 // Note: C1 represents x_significand (UINT64)
1011 BID_UI64DOUBLE tmp1
;
1012 unsigned int x_nr_bits
;
1016 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
1020 // check for NaN or Infinity
1021 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
1023 *pfpsf
|= INVALID_EXCEPTION
;
1024 // return Integer Indefinite
1025 res
= 0x8000000000000000ull
;
1029 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1030 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1031 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
1032 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
1033 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
1034 if (C1
> 9999999999999999ull) { // non-canonical
1039 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
1040 C1
= x
& MASK_BINARY_SIG1
;
1043 // check for zeros (possibly from non-canonical values)
1049 // x is not special and is not zero
1051 // q = nr. of decimal digits in x (1 <= q <= 54)
1052 // determine first the nr. of bits in x
1053 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
1054 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1055 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
1056 tmp1
.d
= (double) (C1
>> 32); // exact conversion
1058 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1059 } else { // x < 2^32
1060 tmp1
.d
= (double) C1
; // exact conversion
1062 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1064 } else { // if x < 2^53
1065 tmp1
.d
= (double) C1
; // exact conversion
1067 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1069 q
= nr_digits
[x_nr_bits
- 1].digits
;
1071 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1072 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
1075 exp
= x_exp
- 398; // unbiased exponent
1077 if ((q
+ exp
) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
1079 *pfpsf
|= INVALID_EXCEPTION
;
1080 // return Integer Indefinite
1081 res
= 0x8000000000000000ull
;
1083 } else if ((q
+ exp
) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
1084 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1085 // so x rounded to an integer may or may not fit in a signed 64-bit int
1086 // the cases that do not fit are identified here; the ones that fit
1087 // fall through and will be handled with other cases further,
1088 // under '1 <= q + exp <= 19'
1089 if (x_sign
) { // if n < 0 and q + exp = 19
1090 // if n <= -2^63 - 1 then n is too large
1091 // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63+1
1092 // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64+2), 1<=q<=16
1093 // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 0x5000000000000000a, 1<=q<=16
1094 // <=> C * 10^(20-q) >= 0x5000000000000000a, 1<=q<=16
1095 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
1096 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
1097 // Note: C1 * 10^(11-q) has 19 or 20 digits; 0x5000000000000000a, has 20
1098 if (C
.w
[1] > 0x05ull
|| (C
.w
[1] == 0x05ull
&& C
.w
[0] >= 0x0aull
)) {
1100 *pfpsf
|= INVALID_EXCEPTION
;
1101 // return Integer Indefinite
1102 res
= 0x8000000000000000ull
;
1105 // else cases that can be rounded to a 64-bit int fall through
1106 // to '1 <= q + exp <= 19'
1107 } else { // if n > 0 and q + exp = 19
1108 // if n > 2^63 - 1 then n is too large
1109 // too large if c(0)c(1)...c(18).c(19)...c(q-1) > 2^63 - 1
1110 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 > 5*(2^64-2), 1<=q<=16
1111 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 > 0x4fffffffffffffff6, 1<=q<=16
1112 // <=> if C * 10^(20-q) > 0x4fffffffffffffff6, 1<=q<=16
1113 C
.w
[1] = 0x0000000000000004ull
;
1114 C
.w
[0] = 0xfffffffffffffff6ull
;
1115 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
1116 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
1117 if (C
.w
[1] > 0x04ull
||
1118 (C
.w
[1] == 0x04ull
&& C
.w
[0] > 0xfffffffffffffff6ull
)) {
1120 *pfpsf
|= INVALID_EXCEPTION
;
1121 // return Integer Indefinite
1122 res
= 0x8000000000000000ull
;
1125 // else cases that can be rounded to a 64-bit int fall through
1126 // to '1 <= q + exp <= 19'
1127 } // end else if n > 0 and q + exp = 19
1128 } // end else if ((q + exp) == 19)
1130 // n is not too large to be converted to int64: -2^63-1 < n < 2^63
1131 // Note: some of the cases tested for above fall through to this point
1132 if ((q
+ exp
) <= 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
1139 } else { // if (1 <= q + exp <= 19, 1 <= q <= 16, -15 <= exp <= 18)
1140 // -2^63-1 < x <= -1 or 1 <= x <= 2^63 - 1 so x can be rounded
1141 // to nearest to a 64-bit signed integer
1142 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 19
1143 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
1144 // chop off ind digits from the lower part of C1
1145 // C1 fits in 64 bits
1146 // calculate C* and f*
1147 // C* is actually floor(C*) in this case
1148 // C* and f* need shifting and masking, as shown by
1149 // shiftright128[] and maskhigh128[]
1151 // kx = 10^(-x) = ten2mk64[ind - 1]
1152 // C* = C1 * 10^(-x)
1153 // the approximation of 10^(-x) was rounded up to 54 bits
1154 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
1156 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
1157 fstar
.w
[0] = P128
.w
[0];
1158 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1159 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1160 // C* = floor(C*) (logical right shift; C has p decimal digits,
1161 // correct by Property 1)
1162 // n = C* * 10^(e+x)
1164 // shift right C* by Ex-64 = shiftright128[ind]
1165 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
1166 Cstar
= Cstar
>> shift
;
1167 // determine inexactness of the rounding of C*
1168 // if (0 < f* < 10^(-x)) then
1169 // the result is exact
1170 // else // if (f* > T*) then
1171 // the result is inexact
1172 if (ind
- 1 <= 2) { // fstar.w[1] is 0
1173 if (fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1174 // ten2mk128trunc[ind -1].w[1] is identical to
1175 // ten2mk128[ind -1].w[1]
1176 if (!x_sign
) { // positive and inexact
1179 } // else the result is exact
1180 } else { // if 3 <= ind - 1 <= 14
1181 if (fstar
.w
[1] || fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1182 // ten2mk128trunc[ind -1].w[1] is identical to
1183 // ten2mk128[ind -1].w[1]
1184 if (!x_sign
) { // positive and inexact
1187 } // else the result is exact
1194 } else if (exp
== 0) {
1196 // res = +/-C (exact)
1201 } else { // if (exp > 0) => 1 <= exp <= 18, 1 <= q <= 16, 2 <= q + exp <= 20
1202 // (the upper limit of 20 on q + exp is due to the fact that
1203 // +/-C * 10^exp is guaranteed to fit in 64 bits)
1204 // res = +/-C * 10^exp (exact)
1206 res
= -C1
* ten2k64
[exp
];
1208 res
= C1
* ten2k64
[exp
];
1214 /*****************************************************************************
1215 * BID64_to_int64_xceil
1216 ****************************************************************************/
1218 #if DECIMAL_CALL_BY_REFERENCE
1220 bid64_to_int64_xceil (SINT64
* pres
, UINT64
* px
1221 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1226 bid64_to_int64_xceil (UINT64 x
1227 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1233 int exp
; // unbiased exponent
1234 // Note: C1 represents x_significand (UINT64)
1235 BID_UI64DOUBLE tmp1
;
1236 unsigned int x_nr_bits
;
1240 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
1244 // check for NaN or Infinity
1245 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
1247 *pfpsf
|= INVALID_EXCEPTION
;
1248 // return Integer Indefinite
1249 res
= 0x8000000000000000ull
;
1253 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1254 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1255 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
1256 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
1257 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
1258 if (C1
> 9999999999999999ull) { // non-canonical
1263 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
1264 C1
= x
& MASK_BINARY_SIG1
;
1267 // check for zeros (possibly from non-canonical values)
1273 // x is not special and is not zero
1275 // q = nr. of decimal digits in x (1 <= q <= 54)
1276 // determine first the nr. of bits in x
1277 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
1278 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1279 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
1280 tmp1
.d
= (double) (C1
>> 32); // exact conversion
1282 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1283 } else { // x < 2^32
1284 tmp1
.d
= (double) C1
; // exact conversion
1286 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1288 } else { // if x < 2^53
1289 tmp1
.d
= (double) C1
; // exact conversion
1291 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1293 q
= nr_digits
[x_nr_bits
- 1].digits
;
1295 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1296 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
1299 exp
= x_exp
- 398; // unbiased exponent
1301 if ((q
+ exp
) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
1303 *pfpsf
|= INVALID_EXCEPTION
;
1304 // return Integer Indefinite
1305 res
= 0x8000000000000000ull
;
1307 } else if ((q
+ exp
) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
1308 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1309 // so x rounded to an integer may or may not fit in a signed 64-bit int
1310 // the cases that do not fit are identified here; the ones that fit
1311 // fall through and will be handled with other cases further,
1312 // under '1 <= q + exp <= 19'
1313 if (x_sign
) { // if n < 0 and q + exp = 19
1314 // if n <= -2^63 - 1 then n is too large
1315 // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63+1
1316 // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64+2), 1<=q<=16
1317 // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 0x5000000000000000a, 1<=q<=16
1318 // <=> C * 10^(20-q) >= 0x5000000000000000a, 1<=q<=16
1319 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
1320 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
1321 // Note: C1 * 10^(11-q) has 19 or 20 digits; 0x5000000000000000a, has 20
1322 if (C
.w
[1] > 0x05ull
|| (C
.w
[1] == 0x05ull
&& C
.w
[0] >= 0x0aull
)) {
1324 *pfpsf
|= INVALID_EXCEPTION
;
1325 // return Integer Indefinite
1326 res
= 0x8000000000000000ull
;
1329 // else cases that can be rounded to a 64-bit int fall through
1330 // to '1 <= q + exp <= 19'
1331 } else { // if n > 0 and q + exp = 19
1332 // if n > 2^63 - 1 then n is too large
1333 // too large if c(0)c(1)...c(18).c(19)...c(q-1) > 2^63 - 1
1334 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 > 5*(2^64-2), 1<=q<=16
1335 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 > 0x4fffffffffffffff6, 1<=q<=16
1336 // <=> if C * 10^(20-q) > 0x4fffffffffffffff6, 1<=q<=16
1337 C
.w
[1] = 0x0000000000000004ull
;
1338 C
.w
[0] = 0xfffffffffffffff6ull
;
1339 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
1340 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
1341 if (C
.w
[1] > 0x04ull
||
1342 (C
.w
[1] == 0x04ull
&& C
.w
[0] > 0xfffffffffffffff6ull
)) {
1344 *pfpsf
|= INVALID_EXCEPTION
;
1345 // return Integer Indefinite
1346 res
= 0x8000000000000000ull
;
1349 // else cases that can be rounded to a 64-bit int fall through
1350 // to '1 <= q + exp <= 19'
1351 } // end else if n > 0 and q + exp = 19
1352 } // end else if ((q + exp) == 19)
1354 // n is not too large to be converted to int64: -2^63-1 < n < 2^63
1355 // Note: some of the cases tested for above fall through to this point
1356 if ((q
+ exp
) <= 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
1358 *pfpsf
|= INEXACT_EXCEPTION
;
1365 } else { // if (1 <= q + exp <= 19, 1 <= q <= 16, -15 <= exp <= 18)
1366 // -2^63-1 < x <= -1 or 1 <= x <= 2^63 - 1 so x can be rounded
1367 // to nearest to a 64-bit signed integer
1368 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 19
1369 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
1370 // chop off ind digits from the lower part of C1
1371 // C1 fits in 64 bits
1372 // calculate C* and f*
1373 // C* is actually floor(C*) in this case
1374 // C* and f* need shifting and masking, as shown by
1375 // shiftright128[] and maskhigh128[]
1377 // kx = 10^(-x) = ten2mk64[ind - 1]
1378 // C* = C1 * 10^(-x)
1379 // the approximation of 10^(-x) was rounded up to 54 bits
1380 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
1382 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
1383 fstar
.w
[0] = P128
.w
[0];
1384 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1385 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1386 // C* = floor(C*) (logical right shift; C has p decimal digits,
1387 // correct by Property 1)
1388 // n = C* * 10^(e+x)
1390 // shift right C* by Ex-64 = shiftright128[ind]
1391 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
1392 Cstar
= Cstar
>> shift
;
1393 // determine inexactness of the rounding of C*
1394 // if (0 < f* < 10^(-x)) then
1395 // the result is exact
1396 // else // if (f* > T*) then
1397 // the result is inexact
1398 if (ind
- 1 <= 2) { // fstar.w[1] is 0
1399 if (fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1400 // ten2mk128trunc[ind -1].w[1] is identical to
1401 // ten2mk128[ind -1].w[1]
1402 if (!x_sign
) { // positive and inexact
1405 // set the inexact flag
1406 *pfpsf
|= INEXACT_EXCEPTION
;
1407 } // else the result is exact
1408 } else { // if 3 <= ind - 1 <= 14
1409 if (fstar
.w
[1] || fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1410 // ten2mk128trunc[ind -1].w[1] is identical to
1411 // ten2mk128[ind -1].w[1]
1412 if (!x_sign
) { // positive and inexact
1415 // set the inexact flag
1416 *pfpsf
|= INEXACT_EXCEPTION
;
1417 } // else the result is exact
1424 } else if (exp
== 0) {
1426 // res = +/-C (exact)
1431 } else { // if (exp > 0) => 1 <= exp <= 18, 1 <= q <= 16, 2 <= q + exp <= 20
1432 // (the upper limit of 20 on q + exp is due to the fact that
1433 // +/-C * 10^exp is guaranteed to fit in 64 bits)
1434 // res = +/-C * 10^exp (exact)
1436 res
= -C1
* ten2k64
[exp
];
1438 res
= C1
* ten2k64
[exp
];
1444 /*****************************************************************************
1445 * BID64_to_int64_int
1446 ****************************************************************************/
1448 #if DECIMAL_CALL_BY_REFERENCE
1450 bid64_to_int64_int (SINT64
* pres
, UINT64
* px
1451 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM
) {
1455 bid64_to_int64_int (UINT64 x
1456 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM
) {
1461 int exp
; // unbiased exponent
1462 // Note: C1 represents x_significand (UINT64)
1463 BID_UI64DOUBLE tmp1
;
1464 unsigned int x_nr_bits
;
1468 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
1471 // check for NaN or Infinity
1472 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
1474 *pfpsf
|= INVALID_EXCEPTION
;
1475 // return Integer Indefinite
1476 res
= 0x8000000000000000ull
;
1480 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1481 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1482 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
1483 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
1484 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
1485 if (C1
> 9999999999999999ull) { // non-canonical
1490 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
1491 C1
= x
& MASK_BINARY_SIG1
;
1494 // check for zeros (possibly from non-canonical values)
1500 // x is not special and is not zero
1502 // q = nr. of decimal digits in x (1 <= q <= 54)
1503 // determine first the nr. of bits in x
1504 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
1505 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1506 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
1507 tmp1
.d
= (double) (C1
>> 32); // exact conversion
1509 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1510 } else { // x < 2^32
1511 tmp1
.d
= (double) C1
; // exact conversion
1513 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1515 } else { // if x < 2^53
1516 tmp1
.d
= (double) C1
; // exact conversion
1518 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1520 q
= nr_digits
[x_nr_bits
- 1].digits
;
1522 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1523 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
1526 exp
= x_exp
- 398; // unbiased exponent
1528 if ((q
+ exp
) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
1530 *pfpsf
|= INVALID_EXCEPTION
;
1531 // return Integer Indefinite
1532 res
= 0x8000000000000000ull
;
1534 } else if ((q
+ exp
) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
1535 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1536 // so x rounded to an integer may or may not fit in a signed 64-bit int
1537 // the cases that do not fit are identified here; the ones that fit
1538 // fall through and will be handled with other cases further,
1539 // under '1 <= q + exp <= 19'
1540 if (x_sign
) { // if n < 0 and q + exp = 19
1541 // if n <= -2^63 - 1 then n is too large
1542 // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63+1
1543 // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64+2), 1<=q<=16
1544 // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 0x5000000000000000a, 1<=q<=16
1545 // <=> C * 10^(20-q) >= 0x5000000000000000a, 1<=q<=16
1546 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
1547 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
1548 // Note: C1 * 10^(11-q) has 19 or 20 digits; 0x5000000000000000a, has 20
1549 if (C
.w
[1] > 0x05ull
|| (C
.w
[1] == 0x05ull
&& C
.w
[0] >= 0x0aull
)) {
1551 *pfpsf
|= INVALID_EXCEPTION
;
1552 // return Integer Indefinite
1553 res
= 0x8000000000000000ull
;
1556 // else cases that can be rounded to a 64-bit int fall through
1557 // to '1 <= q + exp <= 19'
1558 } else { // if n > 0 and q + exp = 19
1559 // if n >= 2^63 then n is too large
1560 // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63
1561 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*2^64, 1<=q<=16
1562 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x50000000000000000, 1<=q<=16
1563 // <=> if C * 10^(20-q) >= 0x50000000000000000, 1<=q<=16
1564 C
.w
[1] = 0x0000000000000005ull
;
1565 C
.w
[0] = 0x0000000000000000ull
;
1566 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
1567 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
1568 if (C
.w
[1] >= 0x05ull
) {
1569 // actually C.w[1] == 0x05ull && C.w[0] >= 0x0000000000000000ull) {
1571 *pfpsf
|= INVALID_EXCEPTION
;
1572 // return Integer Indefinite
1573 res
= 0x8000000000000000ull
;
1576 // else cases that can be rounded to a 64-bit int fall through
1577 // to '1 <= q + exp <= 19'
1578 } // end else if n > 0 and q + exp = 19
1579 } // end else if ((q + exp) == 19)
1581 // n is not too large to be converted to int64: -2^63-1 < n < 2^63
1582 // Note: some of the cases tested for above fall through to this point
1583 if ((q
+ exp
) <= 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
1585 res
= 0x0000000000000000ull
;
1587 } else { // if (1 <= q + exp <= 19, 1 <= q <= 16, -15 <= exp <= 18)
1588 // -2^63-1 < x <= -1 or 1 <= x < 2^63 so x can be rounded
1589 // to nearest to a 64-bit signed integer
1590 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 19
1591 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
1592 // chop off ind digits from the lower part of C1
1593 // C1 fits in 64 bits
1594 // calculate C* and f*
1595 // C* is actually floor(C*) in this case
1596 // C* and f* need shifting and masking, as shown by
1597 // shiftright128[] and maskhigh128[]
1599 // kx = 10^(-x) = ten2mk64[ind - 1]
1600 // C* = C1 * 10^(-x)
1601 // the approximation of 10^(-x) was rounded up to 54 bits
1602 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
1604 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1605 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1606 // C* = floor(C*) (logical right shift; C has p decimal digits,
1607 // correct by Property 1)
1608 // n = C* * 10^(e+x)
1610 // shift right C* by Ex-64 = shiftright128[ind]
1611 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
1612 Cstar
= Cstar
>> shift
;
1618 } else if (exp
== 0) {
1620 // res = +/-C (exact)
1625 } else { // if (exp > 0) => 1 <= exp <= 18, 1 <= q <= 16, 2 <= q + exp <= 20
1626 // (the upper limit of 20 on q + exp is due to the fact that
1627 // +/-C * 10^exp is guaranteed to fit in 64 bits)
1628 // res = +/-C * 10^exp (exact)
1630 res
= -C1
* ten2k64
[exp
];
1632 res
= C1
* ten2k64
[exp
];
1638 /*****************************************************************************
1639 * BID64_to_int64_xint
1640 ****************************************************************************/
1642 #if DECIMAL_CALL_BY_REFERENCE
1644 bid64_to_int64_xint (SINT64
* pres
, UINT64
* px
1645 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM
)
1650 bid64_to_int64_xint (UINT64 x
1651 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM
)
1657 int exp
; // unbiased exponent
1658 // Note: C1 represents x_significand (UINT64)
1659 BID_UI64DOUBLE tmp1
;
1660 unsigned int x_nr_bits
;
1664 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
1668 // check for NaN or Infinity
1669 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
1671 *pfpsf
|= INVALID_EXCEPTION
;
1672 // return Integer Indefinite
1673 res
= 0x8000000000000000ull
;
1677 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1678 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1679 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
1680 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
1681 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
1682 if (C1
> 9999999999999999ull) { // non-canonical
1687 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
1688 C1
= x
& MASK_BINARY_SIG1
;
1691 // check for zeros (possibly from non-canonical values)
1697 // x is not special and is not zero
1699 // q = nr. of decimal digits in x (1 <= q <= 54)
1700 // determine first the nr. of bits in x
1701 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
1702 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1703 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
1704 tmp1
.d
= (double) (C1
>> 32); // exact conversion
1706 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1707 } else { // x < 2^32
1708 tmp1
.d
= (double) C1
; // exact conversion
1710 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1712 } else { // if x < 2^53
1713 tmp1
.d
= (double) C1
; // exact conversion
1715 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1717 q
= nr_digits
[x_nr_bits
- 1].digits
;
1719 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1720 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
1723 exp
= x_exp
- 398; // unbiased exponent
1725 if ((q
+ exp
) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
1727 *pfpsf
|= INVALID_EXCEPTION
;
1728 // return Integer Indefinite
1729 res
= 0x8000000000000000ull
;
1731 } else if ((q
+ exp
) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
1732 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1733 // so x rounded to an integer may or may not fit in a signed 64-bit int
1734 // the cases that do not fit are identified here; the ones that fit
1735 // fall through and will be handled with other cases further,
1736 // under '1 <= q + exp <= 19'
1737 if (x_sign
) { // if n < 0 and q + exp = 19
1738 // if n <= -2^63 - 1 then n is too large
1739 // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63+1
1740 // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64+2), 1<=q<=16
1741 // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 0x5000000000000000a, 1<=q<=16
1742 // <=> C * 10^(20-q) >= 0x5000000000000000a, 1<=q<=16
1743 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
1744 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
1745 // Note: C1 * 10^(11-q) has 19 or 20 digits; 0x5000000000000000a, has 20
1746 if (C
.w
[1] > 0x05ull
|| (C
.w
[1] == 0x05ull
&& C
.w
[0] >= 0x0aull
)) {
1748 *pfpsf
|= INVALID_EXCEPTION
;
1749 // return Integer Indefinite
1750 res
= 0x8000000000000000ull
;
1753 // else cases that can be rounded to a 64-bit int fall through
1754 // to '1 <= q + exp <= 19'
1755 } else { // if n > 0 and q + exp = 19
1756 // if n >= 2^63 then n is too large
1757 // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63
1758 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*2^64, 1<=q<=16
1759 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x50000000000000000, 1<=q<=16
1760 // <=> if C * 10^(20-q) >= 0x50000000000000000, 1<=q<=16
1761 C
.w
[1] = 0x0000000000000005ull
;
1762 C
.w
[0] = 0x0000000000000000ull
;
1763 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
1764 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
1765 if (C
.w
[1] >= 0x05ull
) {
1766 // actually C.w[1] == 0x05ull && C.w[0] >= 0x0000000000000000ull) {
1768 *pfpsf
|= INVALID_EXCEPTION
;
1769 // return Integer Indefinite
1770 res
= 0x8000000000000000ull
;
1773 // else cases that can be rounded to a 64-bit int fall through
1774 // to '1 <= q + exp <= 19'
1775 } // end else if n > 0 and q + exp = 19
1776 } // end else if ((q + exp) == 19)
1778 // n is not too large to be converted to int64: -2^63-1 < n < 2^63
1779 // Note: some of the cases tested for above fall through to this point
1780 if ((q
+ exp
) <= 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
1782 *pfpsf
|= INEXACT_EXCEPTION
;
1784 res
= 0x0000000000000000ull
;
1786 } else { // if (1 <= q + exp <= 19, 1 <= q <= 16, -15 <= exp <= 18)
1787 // -2^63-1 < x <= -1 or 1 <= x < 2^63 so x can be rounded
1788 // to nearest to a 64-bit signed integer
1789 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 19
1790 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
1791 // chop off ind digits from the lower part of C1
1792 // C1 fits in 64 bits
1793 // calculate C* and f*
1794 // C* is actually floor(C*) in this case
1795 // C* and f* need shifting and masking, as shown by
1796 // shiftright128[] and maskhigh128[]
1798 // kx = 10^(-x) = ten2mk64[ind - 1]
1799 // C* = C1 * 10^(-x)
1800 // the approximation of 10^(-x) was rounded up to 54 bits
1801 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
1803 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
1804 fstar
.w
[0] = P128
.w
[0];
1805 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1806 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1807 // C* = floor(C*) (logical right shift; C has p decimal digits,
1808 // correct by Property 1)
1809 // n = C* * 10^(e+x)
1811 // shift right C* by Ex-64 = shiftright128[ind]
1812 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
1813 Cstar
= Cstar
>> shift
;
1814 // determine inexactness of the rounding of C*
1815 // if (0 < f* < 10^(-x)) then
1816 // the result is exact
1817 // else // if (f* > T*) then
1818 // the result is inexact
1819 if (ind
- 1 <= 2) { // fstar.w[1] is 0
1820 if (fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1821 // ten2mk128trunc[ind -1].w[1] is identical to
1822 // ten2mk128[ind -1].w[1]
1823 // set the inexact flag
1824 *pfpsf
|= INEXACT_EXCEPTION
;
1825 } // else the result is exact
1826 } else { // if 3 <= ind - 1 <= 14
1827 if (fstar
.w
[1] || fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
1828 // ten2mk128trunc[ind -1].w[1] is identical to
1829 // ten2mk128[ind -1].w[1]
1830 // set the inexact flag
1831 *pfpsf
|= INEXACT_EXCEPTION
;
1832 } // else the result is exact
1838 } else if (exp
== 0) {
1840 // res = +/-C (exact)
1845 } else { // if (exp > 0) => 1 <= exp <= 18, 1 <= q <= 16, 2 <= q + exp <= 20
1846 // (the upper limit of 20 on q + exp is due to the fact that
1847 // +/-C * 10^exp is guaranteed to fit in 64 bits)
1848 // res = +/-C * 10^exp (exact)
1850 res
= -C1
* ten2k64
[exp
];
1852 res
= C1
* ten2k64
[exp
];
1858 /*****************************************************************************
1859 * BID64_to_int64_rninta
1860 ****************************************************************************/
1862 #if DECIMAL_CALL_BY_REFERENCE
1864 bid64_to_int64_rninta (SINT64
* pres
, UINT64
* px
1865 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1870 bid64_to_int64_rninta (UINT64 x
1871 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1877 int exp
; // unbiased exponent
1878 // Note: C1 represents x_significand (UINT64)
1879 BID_UI64DOUBLE tmp1
;
1880 unsigned int x_nr_bits
;
1884 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
1887 // check for NaN or Infinity
1888 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
1890 *pfpsf
|= INVALID_EXCEPTION
;
1891 // return Integer Indefinite
1892 res
= 0x8000000000000000ull
;
1896 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1897 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1898 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
1899 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
1900 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
1901 if (C1
> 9999999999999999ull) { // non-canonical
1906 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
1907 C1
= x
& MASK_BINARY_SIG1
;
1910 // check for zeros (possibly from non-canonical values)
1916 // x is not special and is not zero
1918 // q = nr. of decimal digits in x (1 <= q <= 54)
1919 // determine first the nr. of bits in x
1920 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
1921 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1922 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
1923 tmp1
.d
= (double) (C1
>> 32); // exact conversion
1925 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1926 } else { // x < 2^32
1927 tmp1
.d
= (double) C1
; // exact conversion
1929 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1931 } else { // if x < 2^53
1932 tmp1
.d
= (double) C1
; // exact conversion
1934 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1936 q
= nr_digits
[x_nr_bits
- 1].digits
;
1938 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1939 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
1942 exp
= x_exp
- 398; // unbiased exponent
1944 if ((q
+ exp
) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
1946 *pfpsf
|= INVALID_EXCEPTION
;
1947 // return Integer Indefinite
1948 res
= 0x8000000000000000ull
;
1950 } else if ((q
+ exp
) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
1951 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1952 // so x rounded to an integer may or may not fit in a signed 64-bit int
1953 // the cases that do not fit are identified here; the ones that fit
1954 // fall through and will be handled with other cases further,
1955 // under '1 <= q + exp <= 19'
1956 if (x_sign
) { // if n < 0 and q + exp = 19
1957 // if n <= -2^63 - 1/2 then n is too large
1958 // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63+1/2
1959 // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64+1), 1<=q<=16
1960 // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 0x50000000000000005, 1<=q<=16
1961 // <=> C * 10^(20-q) >= 0x50000000000000005, 1<=q<=16
1962 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
1963 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
1964 // Note: C1 * 10^(11-q) has 19 or 20 digits; 0x50000000000000005, has 20
1965 if (C
.w
[1] > 0x05ull
|| (C
.w
[1] == 0x05ull
&& C
.w
[0] >= 0x05ull
)) {
1967 *pfpsf
|= INVALID_EXCEPTION
;
1968 // return Integer Indefinite
1969 res
= 0x8000000000000000ull
;
1972 // else cases that can be rounded to a 64-bit int fall through
1973 // to '1 <= q + exp <= 19'
1974 } else { // if n > 0 and q + exp = 19
1975 // if n >= 2^63 - 1/2 then n is too large
1976 // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63-1/2
1977 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64-1), 1<=q<=16
1978 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x4fffffffffffffffb, 1<=q<=16
1979 // <=> if C * 10^(20-q) >= 0x4fffffffffffffffb, 1<=q<=16
1980 C
.w
[1] = 0x0000000000000004ull
;
1981 C
.w
[0] = 0xfffffffffffffffbull
;
1982 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
1983 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
1984 if (C
.w
[1] > 0x04ull
||
1985 (C
.w
[1] == 0x04ull
&& C
.w
[0] >= 0xfffffffffffffffbull
)) {
1987 *pfpsf
|= INVALID_EXCEPTION
;
1988 // return Integer Indefinite
1989 res
= 0x8000000000000000ull
;
1992 // else cases that can be rounded to a 64-bit int fall through
1993 // to '1 <= q + exp <= 19'
1994 } // end else if n > 0 and q + exp = 19
1995 } // end else if ((q + exp) == 19)
1997 // n is not too large to be converted to int64: -2^63-1/2 < n < 2^63-1/2
1998 // Note: some of the cases tested for above fall through to this point
1999 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
2001 res
= 0x0000000000000000ull
;
2003 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
2004 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
2008 ind
= q
- 1; // 0 <= ind <= 15
2009 if (C1
< midpoint64
[ind
]) {
2010 res
= 0x0000000000000000ull
; // return 0
2011 } else if (x_sign
) { // n < 0
2012 res
= 0xffffffffffffffffull
; // return -1
2014 res
= 0x0000000000000001ull
; // return +1
2016 } else { // if (1 <= q + exp <= 19, 1 <= q <= 16, -15 <= exp <= 18)
2017 // -2^63-1/2 < x <= -1 or 1 <= x < 2^63-1/2 so x can be rounded
2018 // to nearest to a 64-bit signed integer
2019 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 19
2020 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
2021 // chop off ind digits from the lower part of C1
2022 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
2023 C1
= C1
+ midpoint64
[ind
- 1];
2024 // calculate C* and f*
2025 // C* is actually floor(C*) in this case
2026 // C* and f* need shifting and masking, as shown by
2027 // shiftright128[] and maskhigh128[]
2029 // kx = 10^(-x) = ten2mk64[ind - 1]
2030 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2031 // the approximation of 10^(-x) was rounded up to 54 bits
2032 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
2034 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
2035 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
2036 // if (0 < f* < 10^(-x)) then the result is a midpoint
2037 // if floor(C*) is even then C* = floor(C*) - logical right
2038 // shift; C* has p decimal digits, correct by Prop. 1)
2039 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
2040 // shift; C* has p decimal digits, correct by Pr. 1)
2042 // C* = floor(C*) (logical right shift; C has p decimal digits,
2043 // correct by Property 1)
2044 // n = C* * 10^(e+x)
2046 // shift right C* by Ex-64 = shiftright128[ind]
2047 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
2048 Cstar
= Cstar
>> shift
;
2050 // if the result was a midpoint it was rounded away from zero
2055 } else if (exp
== 0) {
2057 // res = +/-C (exact)
2062 } else { // if (exp > 0) => 1 <= exp <= 18, 1 <= q <= 16, 2 <= q + exp <= 20
2063 // (the upper limit of 20 on q + exp is due to the fact that
2064 // +/-C * 10^exp is guaranteed to fit in 64 bits)
2065 // res = +/-C * 10^exp (exact)
2067 res
= -C1
* ten2k64
[exp
];
2069 res
= C1
* ten2k64
[exp
];
2075 /*****************************************************************************
2076 * BID64_to_int64_xrninta
2077 ****************************************************************************/
2079 #if DECIMAL_CALL_BY_REFERENCE
2081 bid64_to_int64_xrninta (SINT64
* pres
, UINT64
* px
2082 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
2087 bid64_to_int64_xrninta (UINT64 x
2088 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
2094 int exp
; // unbiased exponent
2095 // Note: C1 represents x_significand (UINT64)
2097 BID_UI64DOUBLE tmp1
;
2098 unsigned int x_nr_bits
;
2102 UINT64 Cstar
; // C* represents up to 16 decimal digits ~ 54 bits
2106 // check for NaN or Infinity
2107 if ((x
& MASK_NAN
) == MASK_NAN
|| (x
& MASK_INF
) == MASK_INF
) {
2109 *pfpsf
|= INVALID_EXCEPTION
;
2110 // return Integer Indefinite
2111 res
= 0x8000000000000000ull
;
2115 x_sign
= x
& MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
2116 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
2117 if ((x
& MASK_STEERING_BITS
) == MASK_STEERING_BITS
) {
2118 x_exp
= (x
& MASK_BINARY_EXPONENT2
) >> 51; // biased
2119 C1
= (x
& MASK_BINARY_SIG2
) | MASK_BINARY_OR2
;
2120 if (C1
> 9999999999999999ull) { // non-canonical
2125 x_exp
= (x
& MASK_BINARY_EXPONENT1
) >> 53; // biased
2126 C1
= x
& MASK_BINARY_SIG1
;
2129 // check for zeros (possibly from non-canonical values)
2135 // x is not special and is not zero
2137 // q = nr. of decimal digits in x (1 <= q <= 54)
2138 // determine first the nr. of bits in x
2139 if (C1
>= 0x0020000000000000ull
) { // x >= 2^53
2140 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2141 if (C1
>= 0x0000000100000000ull
) { // x >= 2^32
2142 tmp1
.d
= (double) (C1
>> 32); // exact conversion
2144 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2145 } else { // x < 2^32
2146 tmp1
.d
= (double) C1
; // exact conversion
2148 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2150 } else { // if x < 2^53
2151 tmp1
.d
= (double) C1
; // exact conversion
2153 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2155 q
= nr_digits
[x_nr_bits
- 1].digits
;
2157 q
= nr_digits
[x_nr_bits
- 1].digits1
;
2158 if (C1
>= nr_digits
[x_nr_bits
- 1].threshold_lo
)
2161 exp
= x_exp
- 398; // unbiased exponent
2163 if ((q
+ exp
) > 19) { // x >= 10^19 ~= 2^63.11... (cannot fit in SINT64)
2165 *pfpsf
|= INVALID_EXCEPTION
;
2166 // return Integer Indefinite
2167 res
= 0x8000000000000000ull
;
2169 } else if ((q
+ exp
) == 19) { // x = c(0)c(1)...c(18).c(19)...c(q-1)
2170 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
2171 // so x rounded to an integer may or may not fit in a signed 64-bit int
2172 // the cases that do not fit are identified here; the ones that fit
2173 // fall through and will be handled with other cases further,
2174 // under '1 <= q + exp <= 19'
2175 if (x_sign
) { // if n < 0 and q + exp = 19
2176 // if n <= -2^63 - 1/2 then n is too large
2177 // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63+1/2
2178 // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64+1), 1<=q<=16
2179 // <=> 0.c(0)c(1)...c(q-1) * 10^20 >= 0x50000000000000005, 1<=q<=16
2180 // <=> C * 10^(20-q) >= 0x50000000000000005, 1<=q<=16
2181 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
2182 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
2183 // Note: C1 * 10^(11-q) has 19 or 20 digits; 0x50000000000000005, has 20
2184 if (C
.w
[1] > 0x05ull
|| (C
.w
[1] == 0x05ull
&& C
.w
[0] >= 0x05ull
)) {
2186 *pfpsf
|= INVALID_EXCEPTION
;
2187 // return Integer Indefinite
2188 res
= 0x8000000000000000ull
;
2191 // else cases that can be rounded to a 64-bit int fall through
2192 // to '1 <= q + exp <= 19'
2193 } else { // if n > 0 and q + exp = 19
2194 // if n >= 2^63 - 1/2 then n is too large
2195 // too large if c(0)c(1)...c(18).c(19)...c(q-1) >= 2^63-1/2
2196 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 5*(2^64-1), 1<=q<=16
2197 // <=> if 0.c(0)c(1)...c(q-1) * 10^20 >= 0x4fffffffffffffffb, 1<=q<=16
2198 // <=> if C * 10^(20-q) >= 0x4fffffffffffffffb, 1<=q<=16
2199 C
.w
[1] = 0x0000000000000004ull
;
2200 C
.w
[0] = 0xfffffffffffffffbull
;
2201 // 1 <= q <= 16 => 4 <= 20-q <= 19 => 10^(20-q) is 64-bit, and so is C1
2202 __mul_64x64_to_128MACH (C
, C1
, ten2k64
[20 - q
]);
2203 if (C
.w
[1] > 0x04ull
||
2204 (C
.w
[1] == 0x04ull
&& C
.w
[0] >= 0xfffffffffffffffbull
)) {
2206 *pfpsf
|= INVALID_EXCEPTION
;
2207 // return Integer Indefinite
2208 res
= 0x8000000000000000ull
;
2211 // else cases that can be rounded to a 64-bit int fall through
2212 // to '1 <= q + exp <= 19'
2213 } // end else if n > 0 and q + exp = 19
2214 } // end else if ((q + exp) == 19)
2216 // n is not too large to be converted to int64: -2^63-1/2 < n < 2^63-1/2
2217 // Note: some of the cases tested for above fall through to this point
2218 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
2220 *pfpsf
|= INEXACT_EXCEPTION
;
2222 res
= 0x0000000000000000ull
;
2224 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
2225 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
2229 ind
= q
- 1; // 0 <= ind <= 15
2230 if (C1
< midpoint64
[ind
]) {
2231 res
= 0x0000000000000000ull
; // return 0
2232 } else if (x_sign
) { // n < 0
2233 res
= 0xffffffffffffffffull
; // return -1
2235 res
= 0x0000000000000001ull
; // return +1
2238 *pfpsf
|= INEXACT_EXCEPTION
;
2239 } else { // if (1 <= q + exp <= 19, 1 <= q <= 16, -15 <= exp <= 18)
2240 // -2^63-1/2 < x <= -1 or 1 <= x < 2^63-1/2 so x can be rounded
2241 // to nearest to a 64-bit signed integer
2242 if (exp
< 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 19
2243 ind
= -exp
; // 1 <= ind <= 15; ind is a synonym for 'x'
2244 // chop off ind digits from the lower part of C1
2245 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
2246 C1
= C1
+ midpoint64
[ind
- 1];
2247 // calculate C* and f*
2248 // C* is actually floor(C*) in this case
2249 // C* and f* need shifting and masking, as shown by
2250 // shiftright128[] and maskhigh128[]
2252 // kx = 10^(-x) = ten2mk64[ind - 1]
2253 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2254 // the approximation of 10^(-x) was rounded up to 54 bits
2255 __mul_64x64_to_128MACH (P128
, C1
, ten2mk64
[ind
- 1]);
2257 fstar
.w
[1] = P128
.w
[1] & maskhigh128
[ind
- 1];
2258 fstar
.w
[0] = P128
.w
[0];
2259 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
2260 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
2261 // if (0 < f* < 10^(-x)) then the result is a midpoint
2262 // if floor(C*) is even then C* = floor(C*) - logical right
2263 // shift; C* has p decimal digits, correct by Prop. 1)
2264 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
2265 // shift; C* has p decimal digits, correct by Pr. 1)
2267 // C* = floor(C*) (logical right shift; C has p decimal digits,
2268 // correct by Property 1)
2269 // n = C* * 10^(e+x)
2271 // shift right C* by Ex-64 = shiftright128[ind]
2272 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 39
2273 Cstar
= Cstar
>> shift
;
2274 // determine inexactness of the rounding of C*
2275 // if (0 < f* - 1/2 < 10^(-x)) then
2276 // the result is exact
2277 // else // if (f* - 1/2 > T*) then
2278 // the result is inexact
2280 if (fstar
.w
[0] > 0x8000000000000000ull
) {
2281 // f* > 1/2 and the result may be exact
2282 tmp64
= fstar
.w
[0] - 0x8000000000000000ull
; // f* - 1/2
2283 if ((tmp64
> ten2mk128trunc
[ind
- 1].w
[1])) {
2284 // ten2mk128trunc[ind -1].w[1] is identical to
2285 // ten2mk128[ind -1].w[1]
2286 // set the inexact flag
2287 *pfpsf
|= INEXACT_EXCEPTION
;
2288 } // else the result is exact
2289 } else { // the result is inexact; f2* <= 1/2
2290 // set the inexact flag
2291 *pfpsf
|= INEXACT_EXCEPTION
;
2293 } else { // if 3 <= ind - 1 <= 14
2294 if (fstar
.w
[1] > onehalf128
[ind
- 1] ||
2295 (fstar
.w
[1] == onehalf128
[ind
- 1] && fstar
.w
[0])) {
2296 // f2* > 1/2 and the result may be exact
2297 // Calculate f2* - 1/2
2298 tmp64
= fstar
.w
[1] - onehalf128
[ind
- 1];
2299 if (tmp64
|| fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[1]) {
2300 // ten2mk128trunc[ind -1].w[1] is identical to
2301 // ten2mk128[ind -1].w[1]
2302 // set the inexact flag
2303 *pfpsf
|= INEXACT_EXCEPTION
;
2304 } // else the result is exact
2305 } else { // the result is inexact; f2* <= 1/2
2306 // set the inexact flag
2307 *pfpsf
|= INEXACT_EXCEPTION
;
2311 // if the result was a midpoint it was rounded away from zero
2316 } else if (exp
== 0) {
2318 // res = +/-C (exact)
2323 } else { // if (exp > 0) => 1 <= exp <= 18, 1 <= q <= 16, 2 <= q + exp <= 20
2324 // (the upper limit of 20 on q + exp is due to the fact that
2325 // +/-C * 10^exp is guaranteed to fit in 64 bits)
2326 // res = +/-C * 10^exp (exact)
2328 res
= -C1
* ten2k64
[exp
];
2330 res
= C1
* ten2k64
[exp
];