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 * BID128_to_uint64_rnint
33 ****************************************************************************/
35 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64
,
36 bid128_to_uint64_rnint
, x
)
41 int exp
; // unbiased exponent
42 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
45 unsigned int x_nr_bits
;
48 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
53 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
54 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
55 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
58 // check for NaN or Infinity
59 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
61 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
62 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
64 *pfpsf
|= INVALID_EXCEPTION
;
65 // return Integer Indefinite
66 res
= 0x8000000000000000ull
;
69 *pfpsf
|= INVALID_EXCEPTION
;
70 // return Integer Indefinite
71 res
= 0x8000000000000000ull
;
74 } else { // x is not a NaN, so it must be infinity
75 if (!x_sign
) { // x is +inf
77 *pfpsf
|= INVALID_EXCEPTION
;
78 // return Integer Indefinite
79 res
= 0x8000000000000000ull
;
82 *pfpsf
|= INVALID_EXCEPTION
;
83 // return Integer Indefinite
84 res
= 0x8000000000000000ull
;
89 // check for non-canonical values (after the check for special values)
90 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
91 || (C1
.w
[1] == 0x0001ed09bead87c0ull
92 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
93 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
94 res
= 0x0000000000000000ull
;
96 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
98 res
= 0x0000000000000000ull
;
100 } else { // x is not special and is not zero
102 // q = nr. of decimal digits in x
103 // determine first the nr. of bits in x
105 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
106 // split the 64-bit value in two 32-bit halves to avoid rounding errors
107 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
108 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
110 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
112 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
114 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
116 } else { // if x < 2^53
117 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
119 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
121 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
122 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
124 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
126 q
= nr_digits
[x_nr_bits
- 1].digits
;
128 q
= nr_digits
[x_nr_bits
- 1].digits1
;
129 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
130 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
131 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
134 exp
= (x_exp
>> 49) - 6176;
136 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
138 *pfpsf
|= INVALID_EXCEPTION
;
139 // return Integer Indefinite
140 res
= 0x8000000000000000ull
;
142 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
143 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
144 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
145 // the cases that do not fit are identified here; the ones that fit
146 // fall through and will be handled with other cases further,
147 // under '1 <= q + exp <= 20'
148 if (x_sign
) { // if n < 0 and q + exp = 20
149 // if n < -1/2 then n cannot be converted to uint64 with RN
150 // too large if c(0)c(1)...c(19).c(20)...c(q-1) > 1/2
151 // <=> 0.c(0)c(1)...c(q-1) * 10^21 > 0x05, 1<=q<=34
152 // <=> C * 10^(21-q) > 0x05, 1<=q<=34
155 if (C1
.w
[1] != 0 || C1
.w
[0] > 0x05ull
) {
157 *pfpsf
|= INVALID_EXCEPTION
;
158 // return Integer Indefinite
159 res
= 0x8000000000000000ull
;
162 // else cases that can be rounded to 64-bit unsigned int fall through
163 // to '1 <= q + exp <= 20'
166 // C * 10^(21-q) > 5 is true because C >= 1 and 10^(21-q) >= 10
167 // if 22 <= q <= 34 => 1 <= q - 21 <= 13
168 // C > 5 * 10^(q-21) is true because C > 2^64 and 5*10^(q-21) < 2^64
170 *pfpsf
|= INVALID_EXCEPTION
;
171 // return Integer Indefinite
172 res
= 0x8000000000000000ull
;
175 } else { // if n > 0 and q + exp = 20
176 // if n >= 2^64 - 1/2 then n is too large
177 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
178 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
179 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
180 // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=34
182 // C * 10^20 >= 0x9fffffffffffffffb
183 __mul_128x64_to_128 (C
, C1
.w
[0], ten2k128
[0]); // 10^20 * C
184 if (C
.w
[1] > 0x09 || (C
.w
[1] == 0x09
185 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
187 *pfpsf
|= INVALID_EXCEPTION
;
188 // return Integer Indefinite
189 res
= 0x8000000000000000ull
;
192 // else cases that can be rounded to a 64-bit int fall through
193 // to '1 <= q + exp <= 20'
194 } else if (q
<= 19) {
195 // C * 10^(21-q) >= 0x9fffffffffffffffb
196 __mul_64x64_to_128MACH (C
, C1
.w
[0], ten2k64
[21 - q
]);
197 if (C
.w
[1] > 0x09 || (C
.w
[1] == 0x09
198 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
200 *pfpsf
|= INVALID_EXCEPTION
;
201 // return Integer Indefinite
202 res
= 0x8000000000000000ull
;
205 // else cases that can be rounded to a 64-bit int fall through
206 // to '1 <= q + exp <= 20'
207 } else if (q
== 20) {
208 // C * 10 >= 0x9fffffffffffffffb <=> C * 2 > 1ffffffffffffffff
209 C
.w
[0] = C1
.w
[0] + C1
.w
[0];
210 C
.w
[1] = C1
.w
[1] + C1
.w
[1];
211 if (C
.w
[0] < C1
.w
[0])
213 if (C
.w
[1] > 0x01 || (C
.w
[1] == 0x01
214 && C
.w
[0] >= 0xffffffffffffffffull
)) {
216 *pfpsf
|= INVALID_EXCEPTION
;
217 // return Integer Indefinite
218 res
= 0x8000000000000000ull
;
221 // else cases that can be rounded to a 64-bit int fall through
222 // to '1 <= q + exp <= 20'
223 } else if (q
== 21) {
224 // C >= 0x9fffffffffffffffb
225 if (C1
.w
[1] > 0x09 || (C1
.w
[1] == 0x09
226 && C1
.w
[0] >= 0xfffffffffffffffbull
)) {
228 *pfpsf
|= INVALID_EXCEPTION
;
229 // return Integer Indefinite
230 res
= 0x8000000000000000ull
;
233 // else cases that can be rounded to a 64-bit int fall through
234 // to '1 <= q + exp <= 20'
235 } else { // if 22 <= q <= 34 => 1 <= q - 21 <= 13
236 // C >= 10^(q-21) * 0x9fffffffffffffffb max 44 bits x 68 bits
238 C
.w
[0] = 0xfffffffffffffffbull
;
239 __mul_128x64_to_128 (C
, ten2k64
[q
- 21], C
);
241 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
243 *pfpsf
|= INVALID_EXCEPTION
;
244 // return Integer Indefinite
245 res
= 0x8000000000000000ull
;
248 // else cases that can be rounded to a 64-bit int fall through
249 // to '1 <= q + exp <= 20'
253 // n is not too large to be converted to int64 if -1/2 <= n < 2^64 - 1/2
254 // Note: some of the cases tested for above fall through to this point
255 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
257 res
= 0x0000000000000000ull
;
259 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
260 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
267 if (ind
<= 18) { // 0 <= ind <= 18
268 if ((C1
.w
[1] == 0) && (C1
.w
[0] <= midpoint64
[ind
])) {
269 res
= 0x0000000000000000ull
; // return 0
270 } else if (!x_sign
) { // n > 0
271 res
= 0x00000001; // return +1
273 res
= 0x8000000000000000ull
;
274 *pfpsf
|= INVALID_EXCEPTION
;
276 } else { // 19 <= ind <= 33
277 if ((C1
.w
[1] < midpoint128
[ind
- 19].w
[1])
278 || ((C1
.w
[1] == midpoint128
[ind
- 19].w
[1])
279 && (C1
.w
[0] <= midpoint128
[ind
- 19].w
[0]))) {
280 res
= 0x0000000000000000ull
; // return 0
281 } else if (!x_sign
) { // n > 0
282 res
= 0x00000001; // return +1
284 res
= 0x8000000000000000ull
;
285 *pfpsf
|= INVALID_EXCEPTION
;
288 } else { // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
289 // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
290 // to nearest to a 64-bit unsigned signed integer
291 if (x_sign
) { // x <= -1
293 *pfpsf
|= INVALID_EXCEPTION
;
294 // return Integer Indefinite
295 res
= 0x8000000000000000ull
;
298 // 1 <= x < 2^64-1/2 so x can be rounded
299 // to nearest to a 64-bit unsigned integer
300 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
301 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
302 // chop off ind digits from the lower part of C1
303 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
306 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
308 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
309 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
313 // calculate C* and f*
314 // C* is actually floor(C*) in this case
315 // C* and f* need shifting and masking, as shown by
316 // shiftright128[] and maskhigh128[]
318 // kx = 10^(-x) = ten2mk128[ind - 1]
319 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
320 // the approximation of 10^(-x) was rounded up to 118 bits
321 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
322 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
323 Cstar
.w
[1] = P256
.w
[3];
324 Cstar
.w
[0] = P256
.w
[2];
326 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
327 fstar
.w
[1] = P256
.w
[1];
328 fstar
.w
[0] = P256
.w
[0];
329 } else { // 22 <= ind - 1 <= 33
331 Cstar
.w
[0] = P256
.w
[3];
332 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
333 fstar
.w
[2] = P256
.w
[2];
334 fstar
.w
[1] = P256
.w
[1];
335 fstar
.w
[0] = P256
.w
[0];
337 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
338 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
339 // if (0 < f* < 10^(-x)) then the result is a midpoint
340 // if floor(C*) is even then C* = floor(C*) - logical right
341 // shift; C* has p decimal digits, correct by Prop. 1)
342 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
343 // shift; C* has p decimal digits, correct by Pr. 1)
345 // C* = floor(C*) (logical right shift; C has p decimal digits,
346 // correct by Property 1)
349 // shift right C* by Ex-128 = shiftright128[ind]
350 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
351 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
353 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
354 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
355 } else { // 22 <= ind - 1 <= 33
356 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
358 // if the result was a midpoint it was rounded away from zero, so
359 // it will need a correction
360 // check for midpoints
361 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
362 && (fstar
.w
[1] || fstar
.w
[0])
363 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
364 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
365 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
366 // the result is a midpoint; round to nearest
367 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
368 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
369 Cstar
.w
[0]--; // Cstar.w[0] is now even
370 } // else MP in [ODD, EVEN]
372 res
= Cstar
.w
[0]; // the result is positive
373 } else if (exp
== 0) {
374 // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
378 // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
379 // res = C * 10^exp (exact) - must fit in 64 bits
380 res
= C1
.w
[0] * ten2k64
[exp
];
388 /*****************************************************************************
389 * BID128_to_uint64_xrnint
390 ****************************************************************************/
392 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64
,
393 bid128_to_uint64_xrnint
, x
)
398 int exp
; // unbiased exponent
399 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
400 UINT64 tmp64
, tmp64A
;
402 unsigned int x_nr_bits
;
405 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
410 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
411 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
412 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
415 // check for NaN or Infinity
416 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
418 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
419 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
421 *pfpsf
|= INVALID_EXCEPTION
;
422 // return Integer Indefinite
423 res
= 0x8000000000000000ull
;
424 } else { // x is QNaN
426 *pfpsf
|= INVALID_EXCEPTION
;
427 // return Integer Indefinite
428 res
= 0x8000000000000000ull
;
431 } else { // x is not a NaN, so it must be infinity
432 if (!x_sign
) { // x is +inf
434 *pfpsf
|= INVALID_EXCEPTION
;
435 // return Integer Indefinite
436 res
= 0x8000000000000000ull
;
437 } else { // x is -inf
439 *pfpsf
|= INVALID_EXCEPTION
;
440 // return Integer Indefinite
441 res
= 0x8000000000000000ull
;
446 // check for non-canonical values (after the check for special values)
447 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
448 || (C1
.w
[1] == 0x0001ed09bead87c0ull
449 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
450 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
451 res
= 0x0000000000000000ull
;
453 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
455 res
= 0x0000000000000000ull
;
457 } else { // x is not special and is not zero
459 // q = nr. of decimal digits in x
460 // determine first the nr. of bits in x
462 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
463 // split the 64-bit value in two 32-bit halves to avoid rounding errors
464 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
465 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
467 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
469 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
471 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
473 } else { // if x < 2^53
474 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
476 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
478 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
479 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
481 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
483 q
= nr_digits
[x_nr_bits
- 1].digits
;
485 q
= nr_digits
[x_nr_bits
- 1].digits1
;
486 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
487 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
488 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
491 exp
= (x_exp
>> 49) - 6176;
493 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
495 *pfpsf
|= INVALID_EXCEPTION
;
496 // return Integer Indefinite
497 res
= 0x8000000000000000ull
;
499 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
500 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
501 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
502 // the cases that do not fit are identified here; the ones that fit
503 // fall through and will be handled with other cases further,
504 // under '1 <= q + exp <= 20'
505 if (x_sign
) { // if n < 0 and q + exp = 20
506 // if n < -1/2 then n cannot be converted to uint64 with RN
507 // too large if c(0)c(1)...c(19).c(20)...c(q-1) > 1/2
508 // <=> 0.c(0)c(1)...c(q-1) * 10^21 > 0x05, 1<=q<=34
509 // <=> C * 10^(21-q) > 0x05, 1<=q<=34
512 if (C1
.w
[1] != 0 || C1
.w
[0] > 0x05ull
) {
514 *pfpsf
|= INVALID_EXCEPTION
;
515 // return Integer Indefinite
516 res
= 0x8000000000000000ull
;
519 // else cases that can be rounded to 64-bit unsigned int fall through
520 // to '1 <= q + exp <= 20'
523 // C * 10^(21-q) > 5 is true because C >= 1 and 10^(21-q) >= 10
524 // if 22 <= q <= 34 => 1 <= q - 21 <= 13
525 // C > 5 * 10^(q-21) is true because C > 2^64 and 5*10^(q-21) < 2^64
527 *pfpsf
|= INVALID_EXCEPTION
;
528 // return Integer Indefinite
529 res
= 0x8000000000000000ull
;
532 } else { // if n > 0 and q + exp = 20
533 // if n >= 2^64 - 1/2 then n is too large
534 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
535 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
536 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
537 // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=34
539 // C * 10^20 >= 0x9fffffffffffffffb
540 __mul_128x64_to_128 (C
, C1
.w
[0], ten2k128
[0]); // 10^20 * C
541 if (C
.w
[1] > 0x09 || (C
.w
[1] == 0x09
542 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
544 *pfpsf
|= INVALID_EXCEPTION
;
545 // return Integer Indefinite
546 res
= 0x8000000000000000ull
;
549 // else cases that can be rounded to a 64-bit int fall through
550 // to '1 <= q + exp <= 20'
551 } else if (q
<= 19) {
552 // C * 10^(21-q) >= 0x9fffffffffffffffb
553 __mul_64x64_to_128MACH (C
, C1
.w
[0], ten2k64
[21 - q
]);
554 if (C
.w
[1] > 0x09 || (C
.w
[1] == 0x09
555 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
557 *pfpsf
|= INVALID_EXCEPTION
;
558 // return Integer Indefinite
559 res
= 0x8000000000000000ull
;
562 // else cases that can be rounded to a 64-bit int fall through
563 // to '1 <= q + exp <= 20'
564 } else if (q
== 20) {
565 // C * 10 >= 0x9fffffffffffffffb <=> C * 2 > 1ffffffffffffffff
566 C
.w
[0] = C1
.w
[0] + C1
.w
[0];
567 C
.w
[1] = C1
.w
[1] + C1
.w
[1];
568 if (C
.w
[0] < C1
.w
[0])
570 if (C
.w
[1] > 0x01 || (C
.w
[1] == 0x01
571 && C
.w
[0] >= 0xffffffffffffffffull
)) {
573 *pfpsf
|= INVALID_EXCEPTION
;
574 // return Integer Indefinite
575 res
= 0x8000000000000000ull
;
578 // else cases that can be rounded to a 64-bit int fall through
579 // to '1 <= q + exp <= 20'
580 } else if (q
== 21) {
581 // C >= 0x9fffffffffffffffb
582 if (C1
.w
[1] > 0x09 || (C1
.w
[1] == 0x09
583 && C1
.w
[0] >= 0xfffffffffffffffbull
)) {
585 *pfpsf
|= INVALID_EXCEPTION
;
586 // return Integer Indefinite
587 res
= 0x8000000000000000ull
;
590 // else cases that can be rounded to a 64-bit int fall through
591 // to '1 <= q + exp <= 20'
592 } else { // if 22 <= q <= 34 => 1 <= q - 21 <= 13
593 // C >= 10^(q-21) * 0x9fffffffffffffffb max 44 bits x 68 bits
595 C
.w
[0] = 0xfffffffffffffffbull
;
596 __mul_128x64_to_128 (C
, ten2k64
[q
- 21], C
);
598 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
600 *pfpsf
|= INVALID_EXCEPTION
;
601 // return Integer Indefinite
602 res
= 0x8000000000000000ull
;
605 // else cases that can be rounded to a 64-bit int fall through
606 // to '1 <= q + exp <= 20'
610 // n is not too large to be converted to int64 if -1/2 <= n < 2^64 - 1/2
611 // Note: some of the cases tested for above fall through to this point
612 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
614 *pfpsf
|= INEXACT_EXCEPTION
;
616 res
= 0x0000000000000000ull
;
618 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
619 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
626 if (ind
<= 18) { // 0 <= ind <= 18
627 if ((C1
.w
[1] == 0) && (C1
.w
[0] <= midpoint64
[ind
])) {
628 res
= 0x0000000000000000ull
; // return 0
629 } else if (!x_sign
) { // n > 0
630 res
= 0x00000001; // return +1
632 res
= 0x8000000000000000ull
;
633 *pfpsf
|= INVALID_EXCEPTION
;
636 } else { // 19 <= ind <= 33
637 if ((C1
.w
[1] < midpoint128
[ind
- 19].w
[1])
638 || ((C1
.w
[1] == midpoint128
[ind
- 19].w
[1])
639 && (C1
.w
[0] <= midpoint128
[ind
- 19].w
[0]))) {
640 res
= 0x0000000000000000ull
; // return 0
641 } else if (!x_sign
) { // n > 0
642 res
= 0x00000001; // return +1
644 res
= 0x8000000000000000ull
;
645 *pfpsf
|= INVALID_EXCEPTION
;
650 *pfpsf
|= INEXACT_EXCEPTION
;
651 } else { // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
652 // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
653 // to nearest to a 64-bit unsigned signed integer
654 if (x_sign
) { // x <= -1
656 *pfpsf
|= INVALID_EXCEPTION
;
657 // return Integer Indefinite
658 res
= 0x8000000000000000ull
;
661 // 1 <= x < 2^64-1/2 so x can be rounded
662 // to nearest to a 64-bit unsigned integer
663 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
664 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
665 // chop off ind digits from the lower part of C1
666 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
669 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
671 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
672 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
676 // calculate C* and f*
677 // C* is actually floor(C*) in this case
678 // C* and f* need shifting and masking, as shown by
679 // shiftright128[] and maskhigh128[]
681 // kx = 10^(-x) = ten2mk128[ind - 1]
682 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
683 // the approximation of 10^(-x) was rounded up to 118 bits
684 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
685 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
686 Cstar
.w
[1] = P256
.w
[3];
687 Cstar
.w
[0] = P256
.w
[2];
689 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
690 fstar
.w
[1] = P256
.w
[1];
691 fstar
.w
[0] = P256
.w
[0];
692 } else { // 22 <= ind - 1 <= 33
694 Cstar
.w
[0] = P256
.w
[3];
695 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
696 fstar
.w
[2] = P256
.w
[2];
697 fstar
.w
[1] = P256
.w
[1];
698 fstar
.w
[0] = P256
.w
[0];
700 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
701 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
702 // if (0 < f* < 10^(-x)) then the result is a midpoint
703 // if floor(C*) is even then C* = floor(C*) - logical right
704 // shift; C* has p decimal digits, correct by Prop. 1)
705 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
706 // shift; C* has p decimal digits, correct by Pr. 1)
708 // C* = floor(C*) (logical right shift; C has p decimal digits,
709 // correct by Property 1)
712 // shift right C* by Ex-128 = shiftright128[ind]
713 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
714 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
716 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
717 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
718 } else { // 22 <= ind - 1 <= 33
719 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
721 // determine inexactness of the rounding of C*
722 // if (0 < f* - 1/2 < 10^(-x)) then
723 // the result is exact
724 // else // if (f* - 1/2 > T*) then
725 // the result is inexact
727 if (fstar
.w
[1] > 0x8000000000000000ull
||
728 (fstar
.w
[1] == 0x8000000000000000ull
729 && fstar
.w
[0] > 0x0ull
)) {
730 // f* > 1/2 and the result may be exact
731 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
732 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
733 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
734 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
735 // set the inexact flag
736 *pfpsf
|= INEXACT_EXCEPTION
;
737 } // else the result is exact
738 } else { // the result is inexact; f2* <= 1/2
739 // set the inexact flag
740 *pfpsf
|= INEXACT_EXCEPTION
;
742 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
743 if (fstar
.w
[3] > 0x0 ||
744 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
745 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
746 (fstar
.w
[1] || fstar
.w
[0]))) {
747 // f2* > 1/2 and the result may be exact
748 // Calculate f2* - 1/2
749 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
751 if (tmp64
> fstar
.w
[2])
754 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
755 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
756 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
757 // set the inexact flag
758 *pfpsf
|= INEXACT_EXCEPTION
;
759 } // else the result is exact
760 } else { // the result is inexact; f2* <= 1/2
761 // set the inexact flag
762 *pfpsf
|= INEXACT_EXCEPTION
;
764 } else { // if 22 <= ind <= 33
765 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
766 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
767 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
768 // f2* > 1/2 and the result may be exact
769 // Calculate f2* - 1/2
770 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
771 if (tmp64
|| fstar
.w
[2]
772 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
773 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
774 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
775 // set the inexact flag
776 *pfpsf
|= INEXACT_EXCEPTION
;
777 } // else the result is exact
778 } else { // the result is inexact; f2* <= 1/2
779 // set the inexact flag
780 *pfpsf
|= INEXACT_EXCEPTION
;
784 // if the result was a midpoint it was rounded away from zero, so
785 // it will need a correction
786 // check for midpoints
787 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
788 && (fstar
.w
[1] || fstar
.w
[0])
789 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
790 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
791 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
792 // the result is a midpoint; round to nearest
793 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
794 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
795 Cstar
.w
[0]--; // Cstar.w[0] is now even
796 } // else MP in [ODD, EVEN]
798 res
= Cstar
.w
[0]; // the result is positive
799 } else if (exp
== 0) {
800 // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
804 // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
805 // res = C * 10^exp (exact) - must fit in 64 bits
806 res
= C1
.w
[0] * ten2k64
[exp
];
814 /*****************************************************************************
815 * BID128_to_uint64_floor
816 ****************************************************************************/
818 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64
,
819 bid128_to_uint64_floor
, x
)
824 int exp
; // unbiased exponent
825 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
827 unsigned int x_nr_bits
;
830 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
834 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
835 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
836 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
839 // check for NaN or Infinity
840 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
842 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
843 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
845 *pfpsf
|= INVALID_EXCEPTION
;
846 // return Integer Indefinite
847 res
= 0x8000000000000000ull
;
848 } else { // x is QNaN
850 *pfpsf
|= INVALID_EXCEPTION
;
851 // return Integer Indefinite
852 res
= 0x8000000000000000ull
;
855 } else { // x is not a NaN, so it must be infinity
856 if (!x_sign
) { // x is +inf
858 *pfpsf
|= INVALID_EXCEPTION
;
859 // return Integer Indefinite
860 res
= 0x8000000000000000ull
;
861 } else { // x is -inf
863 *pfpsf
|= INVALID_EXCEPTION
;
864 // return Integer Indefinite
865 res
= 0x8000000000000000ull
;
870 // check for non-canonical values (after the check for special values)
871 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
872 || (C1
.w
[1] == 0x0001ed09bead87c0ull
873 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
874 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
875 res
= 0x0000000000000000ull
;
877 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
879 res
= 0x0000000000000000ull
;
881 } else { // x is not special and is not zero
883 // if n < 0 then n cannot be converted to uint64 with RM
884 if (x_sign
) { // if n < 0 and q + exp = 20
885 // too large if c(0)c(1)...c(19).c(20)...c(q-1) > 0
887 *pfpsf
|= INVALID_EXCEPTION
;
888 // return Integer Indefinite
889 res
= 0x8000000000000000ull
;
892 // q = nr. of decimal digits in x
893 // determine first the nr. of bits in x
895 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
896 // split the 64-bit value in two 32-bit halves to avoid rounding errors
897 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
898 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
900 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
902 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
904 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
906 } else { // if x < 2^53
907 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
909 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
911 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
912 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
914 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
916 q
= nr_digits
[x_nr_bits
- 1].digits
;
918 q
= nr_digits
[x_nr_bits
- 1].digits1
;
919 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
920 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
921 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
924 exp
= (x_exp
>> 49) - 6176;
925 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
927 *pfpsf
|= INVALID_EXCEPTION
;
928 // return Integer Indefinite
929 res
= 0x8000000000000000ull
;
931 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
932 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
933 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
934 // the cases that do not fit are identified here; the ones that fit
935 // fall through and will be handled with other cases further,
936 // under '1 <= q + exp <= 20'
937 // if n > 0 and q + exp = 20
938 // if n >= 2^64 then n is too large
939 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
940 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
941 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*2^65
942 // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=34
944 // C * 10^20 >= 0xa0000000000000000
945 __mul_128x64_to_128 (C
, C1
.w
[0], ten2k128
[0]); // 10^20 * C
946 if (C
.w
[1] >= 0x0a) {
947 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
949 *pfpsf
|= INVALID_EXCEPTION
;
950 // return Integer Indefinite
951 res
= 0x8000000000000000ull
;
954 // else cases that can be rounded to a 64-bit int fall through
955 // to '1 <= q + exp <= 20'
956 } else if (q
<= 19) {
957 // C * 10^(21-q) >= 0xa0000000000000000
958 __mul_64x64_to_128MACH (C
, C1
.w
[0], ten2k64
[21 - q
]);
959 if (C
.w
[1] >= 0x0a) {
960 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
962 *pfpsf
|= INVALID_EXCEPTION
;
963 // return Integer Indefinite
964 res
= 0x8000000000000000ull
;
967 // else cases that can be rounded to a 64-bit int fall through
968 // to '1 <= q + exp <= 20'
969 } else if (q
== 20) {
970 // C >= 0x10000000000000000
971 if (C1
.w
[1] >= 0x01) {
972 // actually C1.w[1] == 0x01 && C1.w[0] >= 0x0000000000000000ull) {
974 *pfpsf
|= INVALID_EXCEPTION
;
975 // return Integer Indefinite
976 res
= 0x8000000000000000ull
;
979 // else cases that can be rounded to a 64-bit int fall through
980 // to '1 <= q + exp <= 20'
981 } else if (q
== 21) {
982 // C >= 0xa0000000000000000
983 if (C1
.w
[1] >= 0x0a) {
984 // actually C1.w[1] == 0x0a && C1.w[0] >= 0x0000000000000000ull) {
986 *pfpsf
|= INVALID_EXCEPTION
;
987 // return Integer Indefinite
988 res
= 0x8000000000000000ull
;
991 // else cases that can be rounded to a 64-bit int fall through
992 // to '1 <= q + exp <= 20'
993 } else { // if 22 <= q <= 34 => 1 <= q - 21 <= 13
994 // C >= 10^(q-21) * 0xa0000000000000000 max 44 bits x 68 bits
996 C
.w
[0] = 0x0000000000000000ull
;
997 __mul_128x64_to_128 (C
, ten2k64
[q
- 21], C
);
998 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
1000 *pfpsf
|= INVALID_EXCEPTION
;
1001 // return Integer Indefinite
1002 res
= 0x8000000000000000ull
;
1005 // else cases that can be rounded to a 64-bit int fall through
1006 // to '1 <= q + exp <= 20'
1009 // n is not too large to be converted to int64 if 0 <= n < 2^64
1010 // Note: some of the cases tested for above fall through to this point
1011 if ((q
+ exp
) <= 0) { // n = +0.[0...0]c(0)c(1)...c(q-1)
1013 res
= 0x0000000000000000ull
;
1015 } else { // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
1016 // 1 <= x < 2^64 so x can be rounded
1017 // down to a 64-bit unsigned signed integer
1018 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
1019 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
1020 // chop off ind digits from the lower part of C1
1021 // C1 fits in 127 bits
1022 // calculate C* and f*
1023 // C* is actually floor(C*) in this case
1024 // C* and f* need shifting and masking, as shown by
1025 // shiftright128[] and maskhigh128[]
1027 // kx = 10^(-x) = ten2mk128[ind - 1]
1028 // C* = C1 * 10^(-x)
1029 // the approximation of 10^(-x) was rounded up to 118 bits
1030 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
1031 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1032 Cstar
.w
[1] = P256
.w
[3];
1033 Cstar
.w
[0] = P256
.w
[2];
1034 } else { // 22 <= ind - 1 <= 33
1036 Cstar
.w
[0] = P256
.w
[3];
1038 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1039 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1040 // C* = floor(C*) (logical right shift; C has p decimal digits,
1041 // correct by Property 1)
1042 // n = C* * 10^(e+x)
1044 // shift right C* by Ex-128 = shiftright128[ind]
1045 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
1046 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1048 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
1049 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1050 } else { // 22 <= ind - 1 <= 33
1051 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
1053 res
= Cstar
.w
[0]; // the result is positive
1054 } else if (exp
== 0) {
1055 // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
1059 // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
1060 // res = C * 10^exp (exact) - must fit in 64 bits
1061 res
= C1
.w
[0] * ten2k64
[exp
];
1069 /*****************************************************************************
1070 * BID128_to_uint64_xfloor
1071 ****************************************************************************/
1073 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64
,
1074 bid128_to_uint64_xfloor
, x
)
1079 int exp
; // unbiased exponent
1080 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1081 BID_UI64DOUBLE tmp1
;
1082 unsigned int x_nr_bits
;
1085 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
1090 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1091 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
1092 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
1095 // check for NaN or Infinity
1096 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
1098 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
1099 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
1101 *pfpsf
|= INVALID_EXCEPTION
;
1102 // return Integer Indefinite
1103 res
= 0x8000000000000000ull
;
1104 } else { // x is QNaN
1106 *pfpsf
|= INVALID_EXCEPTION
;
1107 // return Integer Indefinite
1108 res
= 0x8000000000000000ull
;
1111 } else { // x is not a NaN, so it must be infinity
1112 if (!x_sign
) { // x is +inf
1114 *pfpsf
|= INVALID_EXCEPTION
;
1115 // return Integer Indefinite
1116 res
= 0x8000000000000000ull
;
1117 } else { // x is -inf
1119 *pfpsf
|= INVALID_EXCEPTION
;
1120 // return Integer Indefinite
1121 res
= 0x8000000000000000ull
;
1126 // check for non-canonical values (after the check for special values)
1127 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
1128 || (C1
.w
[1] == 0x0001ed09bead87c0ull
1129 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
1130 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
1131 res
= 0x0000000000000000ull
;
1133 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
1135 res
= 0x0000000000000000ull
;
1137 } else { // x is not special and is not zero
1139 // if n < 0 then n cannot be converted to uint64 with RM
1140 if (x_sign
) { // if n < 0 and q + exp = 20
1141 // too large if c(0)c(1)...c(19).c(20)...c(q-1) > 0
1143 *pfpsf
|= INVALID_EXCEPTION
;
1144 // return Integer Indefinite
1145 res
= 0x8000000000000000ull
;
1148 // q = nr. of decimal digits in x
1149 // determine first the nr. of bits in x
1151 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
1152 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1153 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
1154 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
1156 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1157 } else { // x < 2^32
1158 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
1160 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1162 } else { // if x < 2^53
1163 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
1165 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1167 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1168 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
1170 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1172 q
= nr_digits
[x_nr_bits
- 1].digits
;
1174 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1175 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
1176 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
1177 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
1180 exp
= (x_exp
>> 49) - 6176;
1181 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
1183 *pfpsf
|= INVALID_EXCEPTION
;
1184 // return Integer Indefinite
1185 res
= 0x8000000000000000ull
;
1187 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
1188 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1189 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
1190 // the cases that do not fit are identified here; the ones that fit
1191 // fall through and will be handled with other cases further,
1192 // under '1 <= q + exp <= 20'
1193 // if n > 0 and q + exp = 20
1194 // if n >= 2^64 then n is too large
1195 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
1196 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
1197 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*2^65
1198 // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=34
1200 // C * 10^20 >= 0xa0000000000000000
1201 __mul_128x64_to_128 (C
, C1
.w
[0], ten2k128
[0]); // 10^20 * C
1202 if (C
.w
[1] >= 0x0a) {
1203 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
1205 *pfpsf
|= INVALID_EXCEPTION
;
1206 // return Integer Indefinite
1207 res
= 0x8000000000000000ull
;
1210 // else cases that can be rounded to a 64-bit int fall through
1211 // to '1 <= q + exp <= 20'
1212 } else if (q
<= 19) {
1213 // C * 10^(21-q) >= 0xa0000000000000000
1214 __mul_64x64_to_128MACH (C
, C1
.w
[0], ten2k64
[21 - q
]);
1215 if (C
.w
[1] >= 0x0a) {
1216 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
1218 *pfpsf
|= INVALID_EXCEPTION
;
1219 // return Integer Indefinite
1220 res
= 0x8000000000000000ull
;
1223 // else cases that can be rounded to a 64-bit int fall through
1224 // to '1 <= q + exp <= 20'
1225 } else if (q
== 20) {
1226 // C >= 0x10000000000000000
1227 if (C1
.w
[1] >= 0x01) {
1228 // actually C1.w[1] == 0x01 && C1.w[0] >= 0x0000000000000000ull) {
1230 *pfpsf
|= INVALID_EXCEPTION
;
1231 // return Integer Indefinite
1232 res
= 0x8000000000000000ull
;
1235 // else cases that can be rounded to a 64-bit int fall through
1236 // to '1 <= q + exp <= 20'
1237 } else if (q
== 21) {
1238 // C >= 0xa0000000000000000
1239 if (C1
.w
[1] >= 0x0a) {
1240 // actually C1.w[1] == 0x0a && C1.w[0] >= 0x0000000000000000ull) {
1242 *pfpsf
|= INVALID_EXCEPTION
;
1243 // return Integer Indefinite
1244 res
= 0x8000000000000000ull
;
1247 // else cases that can be rounded to a 64-bit int fall through
1248 // to '1 <= q + exp <= 20'
1249 } else { // if 22 <= q <= 34 => 1 <= q - 21 <= 13
1250 // C >= 10^(q-21) * 0xa0000000000000000 max 44 bits x 68 bits
1252 C
.w
[0] = 0x0000000000000000ull
;
1253 __mul_128x64_to_128 (C
, ten2k64
[q
- 21], C
);
1254 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
1256 *pfpsf
|= INVALID_EXCEPTION
;
1257 // return Integer Indefinite
1258 res
= 0x8000000000000000ull
;
1261 // else cases that can be rounded to a 64-bit int fall through
1262 // to '1 <= q + exp <= 20'
1265 // n is not too large to be converted to int64 if 0 <= n < 2^64
1266 // Note: some of the cases tested for above fall through to this point
1267 if ((q
+ exp
) <= 0) { // n = +0.[0...0]c(0)c(1)...c(q-1)
1269 *pfpsf
|= INEXACT_EXCEPTION
;
1271 res
= 0x0000000000000000ull
;
1273 } else { // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
1274 // 1 <= x < 2^64 so x can be rounded
1275 // down to a 64-bit unsigned signed integer
1276 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
1277 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
1278 // chop off ind digits from the lower part of C1
1279 // C1 fits in 127 bits
1280 // calculate C* and f*
1281 // C* is actually floor(C*) in this case
1282 // C* and f* need shifting and masking, as shown by
1283 // shiftright128[] and maskhigh128[]
1285 // kx = 10^(-x) = ten2mk128[ind - 1]
1286 // C* = C1 * 10^(-x)
1287 // the approximation of 10^(-x) was rounded up to 118 bits
1288 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
1289 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1290 Cstar
.w
[1] = P256
.w
[3];
1291 Cstar
.w
[0] = P256
.w
[2];
1293 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
1294 fstar
.w
[1] = P256
.w
[1];
1295 fstar
.w
[0] = P256
.w
[0];
1296 } else { // 22 <= ind - 1 <= 33
1298 Cstar
.w
[0] = P256
.w
[3];
1299 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
1300 fstar
.w
[2] = P256
.w
[2];
1301 fstar
.w
[1] = P256
.w
[1];
1302 fstar
.w
[0] = P256
.w
[0];
1304 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1305 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1306 // C* = floor(C*) (logical right shift; C has p decimal digits,
1307 // correct by Property 1)
1308 // n = C* * 10^(e+x)
1310 // shift right C* by Ex-128 = shiftright128[ind]
1311 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
1312 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1314 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
1315 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1316 } else { // 22 <= ind - 1 <= 33
1317 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
1319 // determine inexactness of the rounding of C*
1320 // if (0 < f* < 10^(-x)) then
1321 // the result is exact
1322 // else // if (f* > T*) then
1323 // the result is inexact
1325 if (fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1] ||
1326 (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1] &&
1327 fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1328 // set the inexact flag
1329 *pfpsf
|= INEXACT_EXCEPTION
;
1330 } // else the result is exact
1331 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
1332 if (fstar
.w
[2] || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1333 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1334 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1335 // set the inexact flag
1336 *pfpsf
|= INEXACT_EXCEPTION
;
1337 } // else the result is exact
1338 } else { // if 22 <= ind <= 33
1339 if (fstar
.w
[3] || fstar
.w
[2]
1340 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1341 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1342 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1343 // set the inexact flag
1344 *pfpsf
|= INEXACT_EXCEPTION
;
1345 } // else the result is exact
1348 res
= Cstar
.w
[0]; // the result is positive
1349 } else if (exp
== 0) {
1350 // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
1354 // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
1355 // res = C * 10^exp (exact) - must fit in 64 bits
1356 res
= C1
.w
[0] * ten2k64
[exp
];
1364 /*****************************************************************************
1365 * BID128_to_uint64_ceil
1366 ****************************************************************************/
1368 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64
, bid128_to_uint64_ceil
,
1374 int exp
; // unbiased exponent
1375 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1376 BID_UI64DOUBLE tmp1
;
1377 unsigned int x_nr_bits
;
1380 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
1385 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1386 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
1387 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
1390 // check for NaN or Infinity
1391 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
1393 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
1394 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
1396 *pfpsf
|= INVALID_EXCEPTION
;
1397 // return Integer Indefinite
1398 res
= 0x8000000000000000ull
;
1399 } else { // x is QNaN
1401 *pfpsf
|= INVALID_EXCEPTION
;
1402 // return Integer Indefinite
1403 res
= 0x8000000000000000ull
;
1406 } else { // x is not a NaN, so it must be infinity
1407 if (!x_sign
) { // x is +inf
1409 *pfpsf
|= INVALID_EXCEPTION
;
1410 // return Integer Indefinite
1411 res
= 0x8000000000000000ull
;
1412 } else { // x is -inf
1414 *pfpsf
|= INVALID_EXCEPTION
;
1415 // return Integer Indefinite
1416 res
= 0x8000000000000000ull
;
1421 // check for non-canonical values (after the check for special values)
1422 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
1423 || (C1
.w
[1] == 0x0001ed09bead87c0ull
1424 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
1425 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
1426 res
= 0x0000000000000000ull
;
1428 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
1430 res
= 0x0000000000000000ull
;
1432 } else { // x is not special and is not zero
1434 // q = nr. of decimal digits in x
1435 // determine first the nr. of bits in x
1437 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
1438 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1439 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
1440 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
1442 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1443 } else { // x < 2^32
1444 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
1446 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1448 } else { // if x < 2^53
1449 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
1451 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1453 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1454 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
1456 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1458 q
= nr_digits
[x_nr_bits
- 1].digits
;
1460 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1461 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
1462 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
1463 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
1466 exp
= (x_exp
>> 49) - 6176;
1467 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
1469 *pfpsf
|= INVALID_EXCEPTION
;
1470 // return Integer Indefinite
1471 res
= 0x8000000000000000ull
;
1473 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
1474 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1475 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
1476 // the cases that do not fit are identified here; the ones that fit
1477 // fall through and will be handled with other cases further,
1478 // under '1 <= q + exp <= 20'
1479 if (x_sign
) { // if n < 0 and q + exp = 20
1480 // if n <= -1 then n cannot be converted to uint64 with RZ
1481 // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1
1482 // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x0a, 1<=q<=34
1483 // <=> C * 10^(21-q) >= 0x0a, 1<=q<=34
1486 if (C1
.w
[1] != 0 || C1
.w
[0] >= 0x0aull
) {
1488 *pfpsf
|= INVALID_EXCEPTION
;
1489 // return Integer Indefinite
1490 res
= 0x8000000000000000ull
;
1493 // else cases that can be rounded to 64-bit unsigned int fall through
1494 // to '1 <= q + exp <= 20'
1497 // C * 10^(21-q) >= a is true because C >= 1 and 10^(21-q) >= 10
1498 // if 22 <= q <= 34 => 1 <= q - 21 <= 13
1499 // C >= a * 10^(q-21) is true because C > 2^64 and a*10^(q-21) < 2^64
1501 *pfpsf
|= INVALID_EXCEPTION
;
1502 // return Integer Indefinite
1503 res
= 0x8000000000000000ull
;
1506 } else { // if n > 0 and q + exp = 20
1507 // if n > 2^64 - 1 then n is too large
1508 // <=> c(0)c(1)...c(19).c(20)...c(q-1) > 2^64 - 1
1509 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 > 2^64 - 1
1510 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 > 10 * (2^64 - 1)
1511 // <=> C * 10^(21-q) > 0x9fffffffffffffff6, 1<=q<=34
1513 // C * 10^20 > 0x9fffffffffffffff6
1514 __mul_128x64_to_128 (C
, C1
.w
[0], ten2k128
[0]); // 10^20 * C
1515 if (C
.w
[1] > 0x09 || (C
.w
[1] == 0x09
1516 && C
.w
[0] > 0xfffffffffffffff6ull
)) {
1518 *pfpsf
|= INVALID_EXCEPTION
;
1519 // return Integer Indefinite
1520 res
= 0x8000000000000000ull
;
1523 // else cases that can be rounded to a 64-bit int fall through
1524 // to '1 <= q + exp <= 20'
1525 } else if (q
<= 19) {
1526 // C * 10^(21-q) > 0x9fffffffffffffff6
1527 __mul_64x64_to_128MACH (C
, C1
.w
[0], ten2k64
[21 - q
]);
1528 if (C
.w
[1] > 0x09 || (C
.w
[1] == 0x09
1529 && C
.w
[0] > 0xfffffffffffffff6ull
)) {
1531 *pfpsf
|= INVALID_EXCEPTION
;
1532 // return Integer Indefinite
1533 res
= 0x8000000000000000ull
;
1536 // else cases that can be rounded to a 64-bit int fall through
1537 // to '1 <= q + exp <= 20'
1538 } else if (q
== 20) {
1539 // C > 0xffffffffffffffff
1542 *pfpsf
|= INVALID_EXCEPTION
;
1543 // return Integer Indefinite
1544 res
= 0x8000000000000000ull
;
1547 // else cases that can be rounded to a 64-bit int fall through
1548 // to '1 <= q + exp <= 20'
1549 } else if (q
== 21) {
1550 // C > 0x9fffffffffffffff6
1551 if (C1
.w
[1] > 0x09 || (C1
.w
[1] == 0x09
1552 && C1
.w
[0] > 0xfffffffffffffff6ull
)) {
1554 *pfpsf
|= INVALID_EXCEPTION
;
1555 // return Integer Indefinite
1556 res
= 0x8000000000000000ull
;
1559 // else cases that can be rounded to a 64-bit int fall through
1560 // to '1 <= q + exp <= 20'
1561 } else { // if 22 <= q <= 34 => 1 <= q - 21 <= 13
1562 // C > 10^(q-21) * 0x9fffffffffffffff6 max 44 bits x 68 bits
1564 C
.w
[0] = 0xfffffffffffffff6ull
;
1565 __mul_128x64_to_128 (C
, ten2k64
[q
- 21], C
);
1566 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] > C
.w
[0])) {
1568 *pfpsf
|= INVALID_EXCEPTION
;
1569 // return Integer Indefinite
1570 res
= 0x8000000000000000ull
;
1573 // else cases that can be rounded to a 64-bit int fall through
1574 // to '1 <= q + exp <= 20'
1578 // n is not too large to be converted to int64 if -1 < n <= 2^64 - 1
1579 // Note: some of the cases tested for above fall through to this point
1580 if ((q
+ exp
) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1583 res
= 0x0000000000000000ull
;
1585 res
= 0x0000000000000001ull
;
1587 } else { // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
1588 // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
1589 // to zero to a 64-bit unsigned signed integer
1590 if (x_sign
) { // x <= -1
1592 *pfpsf
|= INVALID_EXCEPTION
;
1593 // return Integer Indefinite
1594 res
= 0x8000000000000000ull
;
1597 // 1 <= x <= 2^64 - 1 so x can be rounded
1598 // to zero to a 64-bit unsigned integer
1599 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
1600 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
1601 // chop off ind digits from the lower part of C1
1602 // C1 fits in 127 bits
1603 // calculate C* and f*
1604 // C* is actually floor(C*) in this case
1605 // C* and f* need shifting and masking, as shown by
1606 // shiftright128[] and maskhigh128[]
1608 // kx = 10^(-x) = ten2mk128[ind - 1]
1609 // C* = C1 * 10^(-x)
1610 // the approximation of 10^(-x) was rounded up to 118 bits
1611 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
1612 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1613 Cstar
.w
[1] = P256
.w
[3];
1614 Cstar
.w
[0] = P256
.w
[2];
1616 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
1617 fstar
.w
[1] = P256
.w
[1];
1618 fstar
.w
[0] = P256
.w
[0];
1619 } else { // 22 <= ind - 1 <= 33
1621 Cstar
.w
[0] = P256
.w
[3];
1622 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
1623 fstar
.w
[2] = P256
.w
[2];
1624 fstar
.w
[1] = P256
.w
[1];
1625 fstar
.w
[0] = P256
.w
[0];
1627 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1628 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1629 // C* = floor(C*) (logical right shift; C has p decimal digits,
1630 // correct by Property 1)
1631 // n = C* * 10^(e+x)
1633 // shift right C* by Ex-128 = shiftright128[ind]
1634 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
1635 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1637 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
1638 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1639 } else { // 22 <= ind - 1 <= 33
1640 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
1642 // if the result is positive and inexact, need to add 1 to it
1644 // determine inexactness of the rounding of C*
1645 // if (0 < f* < 10^(-x)) then
1646 // the result is exact
1647 // else // if (f* > T*) then
1648 // the result is inexact
1650 if (fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1651 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1652 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1653 if (!x_sign
) { // positive and inexact
1655 if (Cstar
.w
[0] == 0x0)
1658 } // else the result is exact
1659 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
1660 if (fstar
.w
[2] || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1661 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1662 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1663 if (!x_sign
) { // positive and inexact
1665 if (Cstar
.w
[0] == 0x0)
1668 } // else the result is exact
1669 } else { // if 22 <= ind <= 33
1670 if (fstar
.w
[3] || fstar
.w
[2]
1671 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1672 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1673 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1674 if (!x_sign
) { // positive and inexact
1676 if (Cstar
.w
[0] == 0x0)
1679 } // else the result is exact
1682 res
= Cstar
.w
[0]; // the result is positive
1683 } else if (exp
== 0) {
1684 // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
1688 // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
1689 // res = C * 10^exp (exact) - must fit in 64 bits
1690 res
= C1
.w
[0] * ten2k64
[exp
];
1698 /*****************************************************************************
1699 * BID128_to_uint64_xceil
1700 ****************************************************************************/
1702 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64
,
1703 bid128_to_uint64_xceil
, x
)
1708 int exp
; // unbiased exponent
1709 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1710 BID_UI64DOUBLE tmp1
;
1711 unsigned int x_nr_bits
;
1714 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
1719 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1720 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
1721 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
1724 // check for NaN or Infinity
1725 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
1727 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
1728 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
1730 *pfpsf
|= INVALID_EXCEPTION
;
1731 // return Integer Indefinite
1732 res
= 0x8000000000000000ull
;
1733 } else { // x is QNaN
1735 *pfpsf
|= INVALID_EXCEPTION
;
1736 // return Integer Indefinite
1737 res
= 0x8000000000000000ull
;
1740 } else { // x is not a NaN, so it must be infinity
1741 if (!x_sign
) { // x is +inf
1743 *pfpsf
|= INVALID_EXCEPTION
;
1744 // return Integer Indefinite
1745 res
= 0x8000000000000000ull
;
1746 } else { // x is -inf
1748 *pfpsf
|= INVALID_EXCEPTION
;
1749 // return Integer Indefinite
1750 res
= 0x8000000000000000ull
;
1755 // check for non-canonical values (after the check for special values)
1756 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
1757 || (C1
.w
[1] == 0x0001ed09bead87c0ull
1758 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
1759 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
1760 res
= 0x0000000000000000ull
;
1762 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
1764 res
= 0x0000000000000000ull
;
1766 } else { // x is not special and is not zero
1768 // q = nr. of decimal digits in x
1769 // determine first the nr. of bits in x
1771 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
1772 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1773 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
1774 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
1776 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1777 } else { // x < 2^32
1778 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
1780 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1782 } else { // if x < 2^53
1783 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
1785 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1787 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1788 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
1790 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1792 q
= nr_digits
[x_nr_bits
- 1].digits
;
1794 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1795 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
1796 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
1797 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
1800 exp
= (x_exp
>> 49) - 6176;
1801 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
1803 *pfpsf
|= INVALID_EXCEPTION
;
1804 // return Integer Indefinite
1805 res
= 0x8000000000000000ull
;
1807 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
1808 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
1809 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
1810 // the cases that do not fit are identified here; the ones that fit
1811 // fall through and will be handled with other cases further,
1812 // under '1 <= q + exp <= 20'
1813 if (x_sign
) { // if n < 0 and q + exp = 20
1814 // if n <= -1 then n cannot be converted to uint64 with RZ
1815 // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1
1816 // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x0a, 1<=q<=34
1817 // <=> C * 10^(21-q) >= 0x0a, 1<=q<=34
1820 if (C1
.w
[1] != 0 || C1
.w
[0] >= 0x0aull
) {
1822 *pfpsf
|= INVALID_EXCEPTION
;
1823 // return Integer Indefinite
1824 res
= 0x8000000000000000ull
;
1827 // else cases that can be rounded to 64-bit unsigned int fall through
1828 // to '1 <= q + exp <= 20'
1831 // C * 10^(21-q) >= a is true because C >= 1 and 10^(21-q) >= 10
1832 // if 22 <= q <= 34 => 1 <= q - 21 <= 13
1833 // C >= a * 10^(q-21) is true because C > 2^64 and a*10^(q-21) < 2^64
1835 *pfpsf
|= INVALID_EXCEPTION
;
1836 // return Integer Indefinite
1837 res
= 0x8000000000000000ull
;
1840 } else { // if n > 0 and q + exp = 20
1841 // if n > 2^64 - 1 then n is too large
1842 // <=> c(0)c(1)...c(19).c(20)...c(q-1) > 2^64 - 1
1843 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 > 2^64 - 1
1844 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 > 10 * (2^64 - 1)
1845 // <=> C * 10^(21-q) > 0x9fffffffffffffff6, 1<=q<=34
1847 // C * 10^20 > 0x9fffffffffffffff6
1848 __mul_128x64_to_128 (C
, C1
.w
[0], ten2k128
[0]); // 10^20 * C
1849 if (C
.w
[1] > 0x09 || (C
.w
[1] == 0x09
1850 && C
.w
[0] > 0xfffffffffffffff6ull
)) {
1852 *pfpsf
|= INVALID_EXCEPTION
;
1853 // return Integer Indefinite
1854 res
= 0x8000000000000000ull
;
1857 // else cases that can be rounded to a 64-bit int fall through
1858 // to '1 <= q + exp <= 20'
1859 } else if (q
<= 19) {
1860 // C * 10^(21-q) > 0x9fffffffffffffff6
1861 __mul_64x64_to_128MACH (C
, C1
.w
[0], ten2k64
[21 - q
]);
1862 if (C
.w
[1] > 0x09 || (C
.w
[1] == 0x09
1863 && C
.w
[0] > 0xfffffffffffffff6ull
)) {
1865 *pfpsf
|= INVALID_EXCEPTION
;
1866 // return Integer Indefinite
1867 res
= 0x8000000000000000ull
;
1870 // else cases that can be rounded to a 64-bit int fall through
1871 // to '1 <= q + exp <= 20'
1872 } else if (q
== 20) {
1873 // C > 0xffffffffffffffff
1876 *pfpsf
|= INVALID_EXCEPTION
;
1877 // return Integer Indefinite
1878 res
= 0x8000000000000000ull
;
1881 // else cases that can be rounded to a 64-bit int fall through
1882 // to '1 <= q + exp <= 20'
1883 } else if (q
== 21) {
1884 // C > 0x9fffffffffffffff6
1885 if (C1
.w
[1] > 0x09 || (C1
.w
[1] == 0x09
1886 && C1
.w
[0] > 0xfffffffffffffff6ull
)) {
1888 *pfpsf
|= INVALID_EXCEPTION
;
1889 // return Integer Indefinite
1890 res
= 0x8000000000000000ull
;
1893 // else cases that can be rounded to a 64-bit int fall through
1894 // to '1 <= q + exp <= 20'
1895 } else { // if 22 <= q <= 34 => 1 <= q - 21 <= 13
1896 // C > 10^(q-21) * 0x9fffffffffffffff6 max 44 bits x 68 bits
1898 C
.w
[0] = 0xfffffffffffffff6ull
;
1899 __mul_128x64_to_128 (C
, ten2k64
[q
- 21], C
);
1900 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] > C
.w
[0])) {
1902 *pfpsf
|= INVALID_EXCEPTION
;
1903 // return Integer Indefinite
1904 res
= 0x8000000000000000ull
;
1907 // else cases that can be rounded to a 64-bit int fall through
1908 // to '1 <= q + exp <= 20'
1912 // n is not too large to be converted to int64 if -1 < n <= 2^64 - 1
1913 // Note: some of the cases tested for above fall through to this point
1914 if ((q
+ exp
) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1916 *pfpsf
|= INEXACT_EXCEPTION
;
1919 res
= 0x0000000000000000ull
;
1921 res
= 0x0000000000000001ull
;
1923 } else { // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
1924 // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
1925 // to zero to a 64-bit unsigned signed integer
1926 if (x_sign
) { // x <= -1
1928 *pfpsf
|= INVALID_EXCEPTION
;
1929 // return Integer Indefinite
1930 res
= 0x8000000000000000ull
;
1933 // 1 <= x <= 2^64 - 1 so x can be rounded
1934 // to zero to a 64-bit unsigned integer
1935 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
1936 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
1937 // chop off ind digits from the lower part of C1
1938 // C1 fits in 127 bits
1939 // calculate C* and f*
1940 // C* is actually floor(C*) in this case
1941 // C* and f* need shifting and masking, as shown by
1942 // shiftright128[] and maskhigh128[]
1944 // kx = 10^(-x) = ten2mk128[ind - 1]
1945 // C* = C1 * 10^(-x)
1946 // the approximation of 10^(-x) was rounded up to 118 bits
1947 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
1948 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1949 Cstar
.w
[1] = P256
.w
[3];
1950 Cstar
.w
[0] = P256
.w
[2];
1952 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
1953 fstar
.w
[1] = P256
.w
[1];
1954 fstar
.w
[0] = P256
.w
[0];
1955 } else { // 22 <= ind - 1 <= 33
1957 Cstar
.w
[0] = P256
.w
[3];
1958 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
1959 fstar
.w
[2] = P256
.w
[2];
1960 fstar
.w
[1] = P256
.w
[1];
1961 fstar
.w
[0] = P256
.w
[0];
1963 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1964 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1965 // C* = floor(C*) (logical right shift; C has p decimal digits,
1966 // correct by Property 1)
1967 // n = C* * 10^(e+x)
1969 // shift right C* by Ex-128 = shiftright128[ind]
1970 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
1971 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1973 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
1974 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1975 } else { // 22 <= ind - 1 <= 33
1976 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
1978 // if the result is positive and inexact, need to add 1 to it
1980 // determine inexactness of the rounding of C*
1981 // if (0 < f* < 10^(-x)) then
1982 // the result is exact
1983 // else // if (f* > T*) then
1984 // the result is inexact
1986 if (fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1987 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1988 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1989 if (!x_sign
) { // positive and inexact
1991 if (Cstar
.w
[0] == 0x0)
1994 // set the inexact flag
1995 *pfpsf
|= INEXACT_EXCEPTION
;
1996 } // else the result is exact
1997 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
1998 if (fstar
.w
[2] || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1999 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2000 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2001 if (!x_sign
) { // positive and inexact
2003 if (Cstar
.w
[0] == 0x0)
2006 // set the inexact flag
2007 *pfpsf
|= INEXACT_EXCEPTION
;
2008 } // else the result is exact
2009 } else { // if 22 <= ind <= 33
2010 if (fstar
.w
[3] || fstar
.w
[2]
2011 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2012 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2013 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2014 if (!x_sign
) { // positive and inexact
2016 if (Cstar
.w
[0] == 0x0)
2019 // set the inexact flag
2020 *pfpsf
|= INEXACT_EXCEPTION
;
2021 } // else the result is exact
2024 res
= Cstar
.w
[0]; // the result is positive
2025 } else if (exp
== 0) {
2026 // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
2030 // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
2031 // res = C * 10^exp (exact) - must fit in 64 bits
2032 res
= C1
.w
[0] * ten2k64
[exp
];
2040 /*****************************************************************************
2041 * BID128_to_uint64_int
2042 ****************************************************************************/
2044 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64
, bid128_to_uint64_int
,
2050 int exp
; // unbiased exponent
2051 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2052 BID_UI64DOUBLE tmp1
;
2053 unsigned int x_nr_bits
;
2056 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
2060 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
2061 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
2062 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
2065 // check for NaN or Infinity
2066 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
2068 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
2069 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
2071 *pfpsf
|= INVALID_EXCEPTION
;
2072 // return Integer Indefinite
2073 res
= 0x8000000000000000ull
;
2074 } else { // x is QNaN
2076 *pfpsf
|= INVALID_EXCEPTION
;
2077 // return Integer Indefinite
2078 res
= 0x8000000000000000ull
;
2081 } else { // x is not a NaN, so it must be infinity
2082 if (!x_sign
) { // x is +inf
2084 *pfpsf
|= INVALID_EXCEPTION
;
2085 // return Integer Indefinite
2086 res
= 0x8000000000000000ull
;
2087 } else { // x is -inf
2089 *pfpsf
|= INVALID_EXCEPTION
;
2090 // return Integer Indefinite
2091 res
= 0x8000000000000000ull
;
2096 // check for non-canonical values (after the check for special values)
2097 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
2098 || (C1
.w
[1] == 0x0001ed09bead87c0ull
2099 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
2100 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
2101 res
= 0x0000000000000000ull
;
2103 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
2105 res
= 0x0000000000000000ull
;
2107 } else { // x is not special and is not zero
2109 // q = nr. of decimal digits in x
2110 // determine first the nr. of bits in x
2112 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
2113 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2114 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
2115 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
2117 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2118 } else { // x < 2^32
2119 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
2121 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2123 } else { // if x < 2^53
2124 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
2126 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2128 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2129 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
2131 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2133 q
= nr_digits
[x_nr_bits
- 1].digits
;
2135 q
= nr_digits
[x_nr_bits
- 1].digits1
;
2136 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
2137 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
2138 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
2141 exp
= (x_exp
>> 49) - 6176;
2143 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
2145 *pfpsf
|= INVALID_EXCEPTION
;
2146 // return Integer Indefinite
2147 res
= 0x8000000000000000ull
;
2149 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
2150 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
2151 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
2152 // the cases that do not fit are identified here; the ones that fit
2153 // fall through and will be handled with other cases further,
2154 // under '1 <= q + exp <= 20'
2155 if (x_sign
) { // if n < 0 and q + exp = 20
2156 // if n <= -1 then n cannot be converted to uint64 with RZ
2157 // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1
2158 // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x0a, 1<=q<=34
2159 // <=> C * 10^(21-q) >= 0x0a, 1<=q<=34
2162 if (C1
.w
[1] != 0 || C1
.w
[0] >= 0x0aull
) {
2164 *pfpsf
|= INVALID_EXCEPTION
;
2165 // return Integer Indefinite
2166 res
= 0x8000000000000000ull
;
2169 // else cases that can be rounded to 64-bit unsigned int fall through
2170 // to '1 <= q + exp <= 20'
2173 // C * 10^(21-q) >= a is true because C >= 1 and 10^(21-q) >= 10
2174 // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2175 // C >= a * 10^(q-21) is true because C > 2^64 and a*10^(q-21) < 2^64
2177 *pfpsf
|= INVALID_EXCEPTION
;
2178 // return Integer Indefinite
2179 res
= 0x8000000000000000ull
;
2182 } else { // if n > 0 and q + exp = 20
2183 // if n >= 2^64 then n is too large
2184 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
2185 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
2186 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*2^65
2187 // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=34
2189 // C * 10^20 >= 0xa0000000000000000
2190 __mul_128x64_to_128 (C
, C1
.w
[0], ten2k128
[0]); // 10^20 * C
2191 if (C
.w
[1] >= 0x0a) {
2192 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
2194 *pfpsf
|= INVALID_EXCEPTION
;
2195 // return Integer Indefinite
2196 res
= 0x8000000000000000ull
;
2199 // else cases that can be rounded to a 64-bit int fall through
2200 // to '1 <= q + exp <= 20'
2201 } else if (q
<= 19) {
2202 // C * 10^(21-q) >= 0xa0000000000000000
2203 __mul_64x64_to_128MACH (C
, C1
.w
[0], ten2k64
[21 - q
]);
2204 if (C
.w
[1] >= 0x0a) {
2205 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
2207 *pfpsf
|= INVALID_EXCEPTION
;
2208 // return Integer Indefinite
2209 res
= 0x8000000000000000ull
;
2212 // else cases that can be rounded to a 64-bit int fall through
2213 // to '1 <= q + exp <= 20'
2214 } else if (q
== 20) {
2215 // C >= 0x10000000000000000
2216 if (C1
.w
[1] >= 0x01) {
2217 // actually C1.w[1] == 0x01 && C1.w[0] >= 0x0000000000000000ull) {
2219 *pfpsf
|= INVALID_EXCEPTION
;
2220 // return Integer Indefinite
2221 res
= 0x8000000000000000ull
;
2224 // else cases that can be rounded to a 64-bit int fall through
2225 // to '1 <= q + exp <= 20'
2226 } else if (q
== 21) {
2227 // C >= 0xa0000000000000000
2228 if (C1
.w
[1] >= 0x0a) {
2229 // actually C1.w[1] == 0x0a && C1.w[0] >= 0x0000000000000000ull) {
2231 *pfpsf
|= INVALID_EXCEPTION
;
2232 // return Integer Indefinite
2233 res
= 0x8000000000000000ull
;
2236 // else cases that can be rounded to a 64-bit int fall through
2237 // to '1 <= q + exp <= 20'
2238 } else { // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2239 // C >= 10^(q-21) * 0xa0000000000000000 max 44 bits x 68 bits
2241 C
.w
[0] = 0x0000000000000000ull
;
2242 __mul_128x64_to_128 (C
, ten2k64
[q
- 21], C
);
2243 if (C1
.w
[1] > C
.w
[1]
2244 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
2246 *pfpsf
|= INVALID_EXCEPTION
;
2247 // return Integer Indefinite
2248 res
= 0x8000000000000000ull
;
2251 // else cases that can be rounded to a 64-bit int fall through
2252 // to '1 <= q + exp <= 20'
2256 // n is not too large to be converted to int64 if -1 < n < 2^64
2257 // Note: some of the cases tested for above fall through to this point
2258 if ((q
+ exp
) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
2260 res
= 0x0000000000000000ull
;
2262 } else { // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
2263 // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
2264 // to zero to a 64-bit unsigned signed integer
2265 if (x_sign
) { // x <= -1
2267 *pfpsf
|= INVALID_EXCEPTION
;
2268 // return Integer Indefinite
2269 res
= 0x8000000000000000ull
;
2272 // 1 <= x < 2^64 so x can be rounded
2273 // to zero to a 64-bit unsigned integer
2274 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
2275 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
2276 // chop off ind digits from the lower part of C1
2277 // C1 fits in 127 bits
2278 // calculate C* and f*
2279 // C* is actually floor(C*) in this case
2280 // C* and f* need shifting and masking, as shown by
2281 // shiftright128[] and maskhigh128[]
2283 // kx = 10^(-x) = ten2mk128[ind - 1]
2284 // C* = C1 * 10^(-x)
2285 // the approximation of 10^(-x) was rounded up to 118 bits
2286 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
2287 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2288 Cstar
.w
[1] = P256
.w
[3];
2289 Cstar
.w
[0] = P256
.w
[2];
2290 } else { // 22 <= ind - 1 <= 33
2292 Cstar
.w
[0] = P256
.w
[3];
2294 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2295 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2296 // C* = floor(C*) (logical right shift; C has p decimal digits,
2297 // correct by Property 1)
2298 // n = C* * 10^(e+x)
2300 // shift right C* by Ex-128 = shiftright128[ind]
2301 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
2302 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2304 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
2305 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2306 } else { // 22 <= ind - 1 <= 33
2307 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
2309 res
= Cstar
.w
[0]; // the result is positive
2310 } else if (exp
== 0) {
2311 // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
2315 // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
2316 // res = C * 10^exp (exact) - must fit in 64 bits
2317 res
= C1
.w
[0] * ten2k64
[exp
];
2325 /*****************************************************************************
2326 * BID128_to_uint64_xint
2327 ****************************************************************************/
2329 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64
, bid128_to_uint64_xint
,
2335 int exp
; // unbiased exponent
2336 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2337 BID_UI64DOUBLE tmp1
;
2338 unsigned int x_nr_bits
;
2341 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
2346 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
2347 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
2348 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
2351 // check for NaN or Infinity
2352 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
2354 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
2355 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
2357 *pfpsf
|= INVALID_EXCEPTION
;
2358 // return Integer Indefinite
2359 res
= 0x8000000000000000ull
;
2360 } else { // x is QNaN
2362 *pfpsf
|= INVALID_EXCEPTION
;
2363 // return Integer Indefinite
2364 res
= 0x8000000000000000ull
;
2367 } else { // x is not a NaN, so it must be infinity
2368 if (!x_sign
) { // x is +inf
2370 *pfpsf
|= INVALID_EXCEPTION
;
2371 // return Integer Indefinite
2372 res
= 0x8000000000000000ull
;
2373 } else { // x is -inf
2375 *pfpsf
|= INVALID_EXCEPTION
;
2376 // return Integer Indefinite
2377 res
= 0x8000000000000000ull
;
2382 // check for non-canonical values (after the check for special values)
2383 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
2384 || (C1
.w
[1] == 0x0001ed09bead87c0ull
2385 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
2386 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
2387 res
= 0x0000000000000000ull
;
2389 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
2391 res
= 0x0000000000000000ull
;
2393 } else { // x is not special and is not zero
2395 // q = nr. of decimal digits in x
2396 // determine first the nr. of bits in x
2398 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
2399 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2400 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
2401 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
2403 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2404 } else { // x < 2^32
2405 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
2407 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2409 } else { // if x < 2^53
2410 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
2412 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2414 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2415 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
2417 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2419 q
= nr_digits
[x_nr_bits
- 1].digits
;
2421 q
= nr_digits
[x_nr_bits
- 1].digits1
;
2422 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
2423 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
2424 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
2427 exp
= (x_exp
>> 49) - 6176;
2428 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
2430 *pfpsf
|= INVALID_EXCEPTION
;
2431 // return Integer Indefinite
2432 res
= 0x8000000000000000ull
;
2434 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
2435 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
2436 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
2437 // the cases that do not fit are identified here; the ones that fit
2438 // fall through and will be handled with other cases further,
2439 // under '1 <= q + exp <= 20'
2440 if (x_sign
) { // if n < 0 and q + exp = 20
2441 // if n <= -1 then n cannot be converted to uint64 with RZ
2442 // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1
2443 // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x0a, 1<=q<=34
2444 // <=> C * 10^(21-q) >= 0x0a, 1<=q<=34
2447 if (C1
.w
[1] != 0 || C1
.w
[0] >= 0x0aull
) {
2449 *pfpsf
|= INVALID_EXCEPTION
;
2450 // return Integer Indefinite
2451 res
= 0x8000000000000000ull
;
2454 // else cases that can be rounded to 64-bit unsigned int fall through
2455 // to '1 <= q + exp <= 20'
2458 // C * 10^(21-q) >= a is true because C >= 1 and 10^(21-q) >= 10
2459 // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2460 // C >= a * 10^(q-21) is true because C > 2^64 and a*10^(q-21) < 2^64
2462 *pfpsf
|= INVALID_EXCEPTION
;
2463 // return Integer Indefinite
2464 res
= 0x8000000000000000ull
;
2467 } else { // if n > 0 and q + exp = 20
2468 // if n >= 2^64 then n is too large
2469 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
2470 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
2471 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*2^65
2472 // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=34
2474 // C * 10^20 >= 0xa0000000000000000
2475 __mul_128x64_to_128 (C
, C1
.w
[0], ten2k128
[0]); // 10^20 * C
2476 if (C
.w
[1] >= 0x0a) {
2477 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
2479 *pfpsf
|= INVALID_EXCEPTION
;
2480 // return Integer Indefinite
2481 res
= 0x8000000000000000ull
;
2484 // else cases that can be rounded to a 64-bit int fall through
2485 // to '1 <= q + exp <= 20'
2486 } else if (q
<= 19) {
2487 // C * 10^(21-q) >= 0xa0000000000000000
2488 __mul_64x64_to_128MACH (C
, C1
.w
[0], ten2k64
[21 - q
]);
2489 if (C
.w
[1] >= 0x0a) {
2490 // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
2492 *pfpsf
|= INVALID_EXCEPTION
;
2493 // return Integer Indefinite
2494 res
= 0x8000000000000000ull
;
2497 // else cases that can be rounded to a 64-bit int fall through
2498 // to '1 <= q + exp <= 20'
2499 } else if (q
== 20) {
2500 // C >= 0x10000000000000000
2501 if (C1
.w
[1] >= 0x01) {
2502 // actually C1.w[1] == 0x01 && C1.w[0] >= 0x0000000000000000ull) {
2504 *pfpsf
|= INVALID_EXCEPTION
;
2505 // return Integer Indefinite
2506 res
= 0x8000000000000000ull
;
2509 // else cases that can be rounded to a 64-bit int fall through
2510 // to '1 <= q + exp <= 20'
2511 } else if (q
== 21) {
2512 // C >= 0xa0000000000000000
2513 if (C1
.w
[1] >= 0x0a) {
2514 // actually C1.w[1] == 0x0a && C1.w[0] >= 0x0000000000000000ull) {
2516 *pfpsf
|= INVALID_EXCEPTION
;
2517 // return Integer Indefinite
2518 res
= 0x8000000000000000ull
;
2521 // else cases that can be rounded to a 64-bit int fall through
2522 // to '1 <= q + exp <= 20'
2523 } else { // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2524 // C >= 10^(q-21) * 0xa0000000000000000 max 44 bits x 68 bits
2526 C
.w
[0] = 0x0000000000000000ull
;
2527 __mul_128x64_to_128 (C
, ten2k64
[q
- 21], C
);
2528 if (C1
.w
[1] > C
.w
[1]
2529 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
2531 *pfpsf
|= INVALID_EXCEPTION
;
2532 // return Integer Indefinite
2533 res
= 0x8000000000000000ull
;
2536 // else cases that can be rounded to a 64-bit int fall through
2537 // to '1 <= q + exp <= 20'
2541 // n is not too large to be converted to int64 if -1 < n < 2^64
2542 // Note: some of the cases tested for above fall through to this point
2543 if ((q
+ exp
) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
2545 *pfpsf
|= INEXACT_EXCEPTION
;
2547 res
= 0x0000000000000000ull
;
2549 } else { // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
2550 // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
2551 // to zero to a 64-bit unsigned signed integer
2552 if (x_sign
) { // x <= -1
2554 *pfpsf
|= INVALID_EXCEPTION
;
2555 // return Integer Indefinite
2556 res
= 0x8000000000000000ull
;
2559 // 1 <= x < 2^64 so x can be rounded
2560 // to zero to a 64-bit unsigned integer
2561 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
2562 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
2563 // chop off ind digits from the lower part of C1
2564 // C1 fits in 127 bits
2565 // calculate C* and f*
2566 // C* is actually floor(C*) in this case
2567 // C* and f* need shifting and masking, as shown by
2568 // shiftright128[] and maskhigh128[]
2570 // kx = 10^(-x) = ten2mk128[ind - 1]
2571 // C* = C1 * 10^(-x)
2572 // the approximation of 10^(-x) was rounded up to 118 bits
2573 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
2574 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2575 Cstar
.w
[1] = P256
.w
[3];
2576 Cstar
.w
[0] = P256
.w
[2];
2578 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
2579 fstar
.w
[1] = P256
.w
[1];
2580 fstar
.w
[0] = P256
.w
[0];
2581 } else { // 22 <= ind - 1 <= 33
2583 Cstar
.w
[0] = P256
.w
[3];
2584 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
2585 fstar
.w
[2] = P256
.w
[2];
2586 fstar
.w
[1] = P256
.w
[1];
2587 fstar
.w
[0] = P256
.w
[0];
2589 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2590 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2591 // C* = floor(C*) (logical right shift; C has p decimal digits,
2592 // correct by Property 1)
2593 // n = C* * 10^(e+x)
2595 // shift right C* by Ex-128 = shiftright128[ind]
2596 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
2597 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2599 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
2600 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2601 } else { // 22 <= ind - 1 <= 33
2602 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
2604 // determine inexactness of the rounding of C*
2605 // if (0 < f* < 10^(-x)) then
2606 // the result is exact
2607 // else // if (f* > T*) then
2608 // the result is inexact
2610 if (fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2611 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2612 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2613 // set the inexact flag
2614 *pfpsf
|= INEXACT_EXCEPTION
;
2615 } // else the result is exact
2616 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
2617 if (fstar
.w
[2] || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2618 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2619 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2620 // set the inexact flag
2621 *pfpsf
|= INEXACT_EXCEPTION
;
2622 } // else the result is exact
2623 } else { // if 22 <= ind <= 33
2624 if (fstar
.w
[3] || fstar
.w
[2]
2625 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2626 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2627 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2628 // set the inexact flag
2629 *pfpsf
|= INEXACT_EXCEPTION
;
2630 } // else the result is exact
2633 res
= Cstar
.w
[0]; // the result is positive
2634 } else if (exp
== 0) {
2635 // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
2639 // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
2640 // res = C * 10^exp (exact) - must fit in 64 bits
2641 res
= C1
.w
[0] * ten2k64
[exp
];
2649 /*****************************************************************************
2650 * BID128_to_uint64_rninta
2651 ****************************************************************************/
2653 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64
,
2654 bid128_to_uint64_rninta
, x
)
2659 int exp
; // unbiased exponent
2660 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2662 BID_UI64DOUBLE tmp1
;
2663 unsigned int x_nr_bits
;
2666 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
2670 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
2671 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
2672 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
2675 // check for NaN or Infinity
2676 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
2678 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
2679 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
2681 *pfpsf
|= INVALID_EXCEPTION
;
2682 // return Integer Indefinite
2683 res
= 0x8000000000000000ull
;
2684 } else { // x is QNaN
2686 *pfpsf
|= INVALID_EXCEPTION
;
2687 // return Integer Indefinite
2688 res
= 0x8000000000000000ull
;
2691 } else { // x is not a NaN, so it must be infinity
2692 if (!x_sign
) { // x is +inf
2694 *pfpsf
|= INVALID_EXCEPTION
;
2695 // return Integer Indefinite
2696 res
= 0x8000000000000000ull
;
2697 } else { // x is -inf
2699 *pfpsf
|= INVALID_EXCEPTION
;
2700 // return Integer Indefinite
2701 res
= 0x8000000000000000ull
;
2706 // check for non-canonical values (after the check for special values)
2707 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
2708 || (C1
.w
[1] == 0x0001ed09bead87c0ull
2709 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
2710 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
2711 res
= 0x0000000000000000ull
;
2713 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
2715 res
= 0x0000000000000000ull
;
2717 } else { // x is not special and is not zero
2719 // q = nr. of decimal digits in x
2720 // determine first the nr. of bits in x
2722 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
2723 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2724 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
2725 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
2727 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2728 } else { // x < 2^32
2729 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
2731 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2733 } else { // if x < 2^53
2734 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
2736 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2738 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2739 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
2741 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2743 q
= nr_digits
[x_nr_bits
- 1].digits
;
2745 q
= nr_digits
[x_nr_bits
- 1].digits1
;
2746 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
2747 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
2748 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
2751 exp
= (x_exp
>> 49) - 6176;
2752 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
2754 *pfpsf
|= INVALID_EXCEPTION
;
2755 // return Integer Indefinite
2756 res
= 0x8000000000000000ull
;
2758 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
2759 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
2760 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
2761 // the cases that do not fit are identified here; the ones that fit
2762 // fall through and will be handled with other cases further,
2763 // under '1 <= q + exp <= 20'
2764 if (x_sign
) { // if n < 0 and q + exp = 20
2765 // if n <= -1/2 then n cannot be converted to uint64 with RN
2766 // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1/2
2767 // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x05, 1<=q<=34
2768 // <=> C * 10^(21-q) >= 0x05, 1<=q<=34
2771 if (C1
.w
[1] != 0 || C1
.w
[0] >= 0x05ull
) {
2773 *pfpsf
|= INVALID_EXCEPTION
;
2774 // return Integer Indefinite
2775 res
= 0x8000000000000000ull
;
2778 // else cases that can be rounded to 64-bit unsigned int fall through
2779 // to '1 <= q + exp <= 20'
2782 // C * 10^(21-q) >= 5 is true because C >= 1 and 10^(21-q) >= 10
2783 // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2784 // C >= 5 * 10^(q-21) is true because C > 2^64 and 5*10^(q-21) < 2^64
2786 *pfpsf
|= INVALID_EXCEPTION
;
2787 // return Integer Indefinite
2788 res
= 0x8000000000000000ull
;
2791 } else { // if n > 0 and q + exp = 20
2792 // if n >= 2^64 - 1/2 then n is too large
2793 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
2794 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
2795 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
2796 // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=34
2798 // C * 10^20 >= 0x9fffffffffffffffb
2799 __mul_128x64_to_128 (C
, C1
.w
[0], ten2k128
[0]); // 10^20 * C
2800 if (C
.w
[1] > 0x09 || (C
.w
[1] == 0x09
2801 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
2803 *pfpsf
|= INVALID_EXCEPTION
;
2804 // return Integer Indefinite
2805 res
= 0x8000000000000000ull
;
2808 // else cases that can be rounded to a 64-bit int fall through
2809 // to '1 <= q + exp <= 20'
2810 } else if (q
<= 19) {
2811 // C * 10^(21-q) >= 0x9fffffffffffffffb
2812 __mul_64x64_to_128MACH (C
, C1
.w
[0], ten2k64
[21 - q
]);
2813 if (C
.w
[1] > 0x09 || (C
.w
[1] == 0x09
2814 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
2816 *pfpsf
|= INVALID_EXCEPTION
;
2817 // return Integer Indefinite
2818 res
= 0x8000000000000000ull
;
2821 // else cases that can be rounded to a 64-bit int fall through
2822 // to '1 <= q + exp <= 20'
2823 } else if (q
== 20) {
2824 // C * 10 >= 0x9fffffffffffffffb <=> C * 2 > 1ffffffffffffffff
2825 C
.w
[0] = C1
.w
[0] + C1
.w
[0];
2826 C
.w
[1] = C1
.w
[1] + C1
.w
[1];
2827 if (C
.w
[0] < C1
.w
[0])
2829 if (C
.w
[1] > 0x01 || (C
.w
[1] == 0x01
2830 && C
.w
[0] >= 0xffffffffffffffffull
)) {
2832 *pfpsf
|= INVALID_EXCEPTION
;
2833 // return Integer Indefinite
2834 res
= 0x8000000000000000ull
;
2837 // else cases that can be rounded to a 64-bit int fall through
2838 // to '1 <= q + exp <= 20'
2839 } else if (q
== 21) {
2840 // C >= 0x9fffffffffffffffb
2841 if (C1
.w
[1] > 0x09 || (C1
.w
[1] == 0x09
2842 && C1
.w
[0] >= 0xfffffffffffffffbull
)) {
2844 *pfpsf
|= INVALID_EXCEPTION
;
2845 // return Integer Indefinite
2846 res
= 0x8000000000000000ull
;
2849 // else cases that can be rounded to a 64-bit int fall through
2850 // to '1 <= q + exp <= 20'
2851 } else { // if 22 <= q <= 34 => 1 <= q - 21 <= 13
2852 // C >= 10^(q-21) * 0x9fffffffffffffffb max 44 bits x 68 bits
2854 C
.w
[0] = 0xfffffffffffffffbull
;
2855 __mul_128x64_to_128 (C
, ten2k64
[q
- 21], C
);
2856 if (C1
.w
[1] > C
.w
[1]
2857 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
2859 *pfpsf
|= INVALID_EXCEPTION
;
2860 // return Integer Indefinite
2861 res
= 0x8000000000000000ull
;
2864 // else cases that can be rounded to a 64-bit int fall through
2865 // to '1 <= q + exp <= 20'
2869 // n is not too large to be converted to int64 if -1/2 < n < 2^64 - 1/2
2870 // Note: some of the cases tested for above fall through to this point
2871 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
2873 res
= 0x0000000000000000ull
;
2875 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
2876 // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
2883 if (ind
<= 18) { // 0 <= ind <= 18
2884 if ((C1
.w
[1] == 0) && (C1
.w
[0] < midpoint64
[ind
])) {
2885 res
= 0x0000000000000000ull
; // return 0
2886 } else if (!x_sign
) { // n > 0
2887 res
= 0x00000001; // return +1
2890 *pfpsf
|= INVALID_EXCEPTION
;
2891 // return Integer Indefinite
2892 res
= 0x8000000000000000ull
;
2895 } else { // 19 <= ind <= 33
2896 if ((C1
.w
[1] < midpoint128
[ind
- 19].w
[1])
2897 || ((C1
.w
[1] == midpoint128
[ind
- 19].w
[1])
2898 && (C1
.w
[0] < midpoint128
[ind
- 19].w
[0]))) {
2899 res
= 0x0000000000000000ull
; // return 0
2900 } else if (!x_sign
) { // n > 0
2901 res
= 0x00000001; // return +1
2904 *pfpsf
|= INVALID_EXCEPTION
;
2905 // return Integer Indefinite
2906 res
= 0x8000000000000000ull
;
2910 } else { // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
2911 // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
2912 // to nearest to a 64-bit unsigned signed integer
2913 if (x_sign
) { // x <= -1
2915 *pfpsf
|= INVALID_EXCEPTION
;
2916 // return Integer Indefinite
2917 res
= 0x8000000000000000ull
;
2920 // 1 <= x < 2^64-1/2 so x can be rounded
2921 // to nearest to a 64-bit unsigned integer
2922 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
2923 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
2924 // chop off ind digits from the lower part of C1
2925 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
2928 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
2930 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
2931 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
2933 if (C1
.w
[0] < tmp64
)
2935 // calculate C* and f*
2936 // C* is actually floor(C*) in this case
2937 // C* and f* need shifting and masking, as shown by
2938 // shiftright128[] and maskhigh128[]
2940 // kx = 10^(-x) = ten2mk128[ind - 1]
2941 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2942 // the approximation of 10^(-x) was rounded up to 118 bits
2943 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
2944 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2945 Cstar
.w
[1] = P256
.w
[3];
2946 Cstar
.w
[0] = P256
.w
[2];
2947 } else { // 22 <= ind - 1 <= 33
2949 Cstar
.w
[0] = P256
.w
[3];
2951 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2952 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2953 // if (0 < f* < 10^(-x)) then the result is a midpoint
2954 // if floor(C*) is even then C* = floor(C*) - logical right
2955 // shift; C* has p decimal digits, correct by Prop. 1)
2956 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
2957 // shift; C* has p decimal digits, correct by Pr. 1)
2959 // C* = floor(C*) (logical right shift; C has p decimal digits,
2960 // correct by Property 1)
2961 // n = C* * 10^(e+x)
2963 // shift right C* by Ex-128 = shiftright128[ind]
2964 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
2965 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2967 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
2968 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2969 } else { // 22 <= ind - 1 <= 33
2970 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
2973 // if the result was a midpoint it was rounded away from zero
2974 res
= Cstar
.w
[0]; // the result is positive
2975 } else if (exp
== 0) {
2976 // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
2980 // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
2981 // res = C * 10^exp (exact) - must fit in 64 bits
2982 res
= C1
.w
[0] * ten2k64
[exp
];
2990 /*****************************************************************************
2991 * BID128_to_uint64_xrninta
2992 ****************************************************************************/
2994 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64
,
2995 bid128_to_uint64_xrninta
, x
)
3000 int exp
; // unbiased exponent
3001 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
3002 UINT64 tmp64
, tmp64A
;
3003 BID_UI64DOUBLE tmp1
;
3004 unsigned int x_nr_bits
;
3007 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
3012 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
3013 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
3014 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
3017 // check for NaN or Infinity
3018 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
3020 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
3021 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
3023 *pfpsf
|= INVALID_EXCEPTION
;
3024 // return Integer Indefinite
3025 res
= 0x8000000000000000ull
;
3026 } else { // x is QNaN
3028 *pfpsf
|= INVALID_EXCEPTION
;
3029 // return Integer Indefinite
3030 res
= 0x8000000000000000ull
;
3033 } else { // x is not a NaN, so it must be infinity
3034 if (!x_sign
) { // x is +inf
3036 *pfpsf
|= INVALID_EXCEPTION
;
3037 // return Integer Indefinite
3038 res
= 0x8000000000000000ull
;
3039 } else { // x is -inf
3041 *pfpsf
|= INVALID_EXCEPTION
;
3042 // return Integer Indefinite
3043 res
= 0x8000000000000000ull
;
3048 // check for non-canonical values (after the check for special values)
3049 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
3050 || (C1
.w
[1] == 0x0001ed09bead87c0ull
3051 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
3052 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
3053 res
= 0x0000000000000000ull
;
3055 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
3057 res
= 0x0000000000000000ull
;
3059 } else { // x is not special and is not zero
3061 // q = nr. of decimal digits in x
3062 // determine first the nr. of bits in x
3064 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
3065 // split the 64-bit value in two 32-bit halves to avoid rounding errors
3066 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
3067 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
3069 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3070 } else { // x < 2^32
3071 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
3073 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3075 } else { // if x < 2^53
3076 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
3078 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3080 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
3081 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
3083 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3085 q
= nr_digits
[x_nr_bits
- 1].digits
;
3087 q
= nr_digits
[x_nr_bits
- 1].digits1
;
3088 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
3089 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
3090 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
3093 exp
= (x_exp
>> 49) - 6176;
3095 if ((q
+ exp
) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
3097 *pfpsf
|= INVALID_EXCEPTION
;
3098 // return Integer Indefinite
3099 res
= 0x8000000000000000ull
;
3101 } else if ((q
+ exp
) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
3102 // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
3103 // so x rounded to an integer may or may not fit in an unsigned 64-bit int
3104 // the cases that do not fit are identified here; the ones that fit
3105 // fall through and will be handled with other cases further,
3106 // under '1 <= q + exp <= 20'
3107 if (x_sign
) { // if n < 0 and q + exp = 20
3108 // if n <= -1/2 then n cannot be converted to uint64 with RN
3109 // too large if c(0)c(1)...c(19).c(20)...c(q-1) >= 1/2
3110 // <=> 0.c(0)c(1)...c(q-1) * 10^21 >= 0x05, 1<=q<=34
3111 // <=> C * 10^(21-q) >= 0x05, 1<=q<=34
3114 if (C1
.w
[1] != 0 || C1
.w
[0] >= 0x05ull
) {
3116 *pfpsf
|= INVALID_EXCEPTION
;
3117 // return Integer Indefinite
3118 res
= 0x8000000000000000ull
;
3121 // else cases that can be rounded to 64-bit unsigned int fall through
3122 // to '1 <= q + exp <= 20'
3125 // C * 10^(21-q) >= 5 is true because C >= 1 and 10^(21-q) >= 10
3126 // if 22 <= q <= 34 => 1 <= q - 21 <= 13
3127 // C >= 5 * 10^(q-21) is true because C > 2^64 and 5*10^(q-21) < 2^64
3129 *pfpsf
|= INVALID_EXCEPTION
;
3130 // return Integer Indefinite
3131 res
= 0x8000000000000000ull
;
3134 } else { // if n > 0 and q + exp = 20
3135 // if n >= 2^64 - 1/2 then n is too large
3136 // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
3137 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
3138 // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
3139 // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=34
3141 // C * 10^20 >= 0x9fffffffffffffffb
3142 __mul_128x64_to_128 (C
, C1
.w
[0], ten2k128
[0]); // 10^20 * C
3143 if (C
.w
[1] > 0x09 || (C
.w
[1] == 0x09
3144 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
3146 *pfpsf
|= INVALID_EXCEPTION
;
3147 // return Integer Indefinite
3148 res
= 0x8000000000000000ull
;
3151 // else cases that can be rounded to a 64-bit int fall through
3152 // to '1 <= q + exp <= 20'
3153 } else if (q
<= 19) {
3154 // C * 10^(21-q) >= 0x9fffffffffffffffb
3155 __mul_64x64_to_128MACH (C
, C1
.w
[0], ten2k64
[21 - q
]);
3156 if (C
.w
[1] > 0x09 || (C
.w
[1] == 0x09
3157 && C
.w
[0] >= 0xfffffffffffffffbull
)) {
3159 *pfpsf
|= INVALID_EXCEPTION
;
3160 // return Integer Indefinite
3161 res
= 0x8000000000000000ull
;
3164 // else cases that can be rounded to a 64-bit int fall through
3165 // to '1 <= q + exp <= 20'
3166 } else if (q
== 20) {
3167 // C * 10 >= 0x9fffffffffffffffb <=> C * 2 > 1ffffffffffffffff
3168 C
.w
[0] = C1
.w
[0] + C1
.w
[0];
3169 C
.w
[1] = C1
.w
[1] + C1
.w
[1];
3170 if (C
.w
[0] < C1
.w
[0])
3172 if (C
.w
[1] > 0x01 || (C
.w
[1] == 0x01
3173 && C
.w
[0] >= 0xffffffffffffffffull
)) {
3175 *pfpsf
|= INVALID_EXCEPTION
;
3176 // return Integer Indefinite
3177 res
= 0x8000000000000000ull
;
3180 // else cases that can be rounded to a 64-bit int fall through
3181 // to '1 <= q + exp <= 20'
3182 } else if (q
== 21) {
3183 // C >= 0x9fffffffffffffffb
3184 if (C1
.w
[1] > 0x09 || (C1
.w
[1] == 0x09
3185 && C1
.w
[0] >= 0xfffffffffffffffbull
)) {
3187 *pfpsf
|= INVALID_EXCEPTION
;
3188 // return Integer Indefinite
3189 res
= 0x8000000000000000ull
;
3192 // else cases that can be rounded to a 64-bit int fall through
3193 // to '1 <= q + exp <= 20'
3194 } else { // if 22 <= q <= 34 => 1 <= q - 21 <= 13
3195 // C >= 10^(q-21) * 0x9fffffffffffffffb max 44 bits x 68 bits
3197 C
.w
[0] = 0xfffffffffffffffbull
;
3198 __mul_128x64_to_128 (C
, ten2k64
[q
- 21], C
);
3199 if (C1
.w
[1] > C
.w
[1]
3200 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
3202 *pfpsf
|= INVALID_EXCEPTION
;
3203 // return Integer Indefinite
3204 res
= 0x8000000000000000ull
;
3207 // else cases that can be rounded to a 64-bit int fall through
3208 // to '1 <= q + exp <= 20'
3212 // n is not too large to be converted to int64 if -1/2 < n < 2^64 - 1/2
3213 // Note: some of the cases tested for above fall through to this point
3214 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
3216 *pfpsf
|= INEXACT_EXCEPTION
;
3218 res
= 0x0000000000000000ull
;
3220 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
3221 // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
3228 if (ind
<= 18) { // 0 <= ind <= 18
3229 if ((C1
.w
[1] == 0) && (C1
.w
[0] < midpoint64
[ind
])) {
3230 res
= 0x0000000000000000ull
; // return 0
3231 } else if (!x_sign
) { // n > 0
3232 res
= 0x00000001; // return +1
3234 res
= 0x8000000000000000ull
;
3236 *pfpsf
|= INVALID_EXCEPTION
;
3237 // return Integer Indefinite
3238 res
= 0x8000000000000000ull
;
3241 } else { // 19 <= ind <= 33
3242 if ((C1
.w
[1] < midpoint128
[ind
- 19].w
[1])
3243 || ((C1
.w
[1] == midpoint128
[ind
- 19].w
[1])
3244 && (C1
.w
[0] < midpoint128
[ind
- 19].w
[0]))) {
3245 res
= 0x0000000000000000ull
; // return 0
3246 } else if (!x_sign
) { // n > 0
3247 res
= 0x00000001; // return +1
3249 res
= 0x8000000000000000ull
;
3250 *pfpsf
|= INVALID_EXCEPTION
;
3251 // return Integer Indefinite
3252 res
= 0x8000000000000000ull
;
3257 *pfpsf
|= INEXACT_EXCEPTION
;
3258 } else { // if (1 <= q + exp <= 20, 1 <= q <= 34, -33 <= exp <= 19)
3259 // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
3260 // to nearest to a 64-bit unsigned signed integer
3261 if (x_sign
) { // x <= -1
3263 *pfpsf
|= INVALID_EXCEPTION
;
3264 // return Integer Indefinite
3265 res
= 0x8000000000000000ull
;
3268 // 1 <= x < 2^64-1/2 so x can be rounded
3269 // to nearest to a 64-bit unsigned integer
3270 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 20
3271 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
3272 // chop off ind digits from the lower part of C1
3273 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
3276 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
3278 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
3279 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
3281 if (C1
.w
[0] < tmp64
)
3283 // calculate C* and f*
3284 // C* is actually floor(C*) in this case
3285 // C* and f* need shifting and masking, as shown by
3286 // shiftright128[] and maskhigh128[]
3288 // kx = 10^(-x) = ten2mk128[ind - 1]
3289 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
3290 // the approximation of 10^(-x) was rounded up to 118 bits
3291 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
3292 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
3293 Cstar
.w
[1] = P256
.w
[3];
3294 Cstar
.w
[0] = P256
.w
[2];
3296 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
3297 fstar
.w
[1] = P256
.w
[1];
3298 fstar
.w
[0] = P256
.w
[0];
3299 } else { // 22 <= ind - 1 <= 33
3301 Cstar
.w
[0] = P256
.w
[3];
3302 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
3303 fstar
.w
[2] = P256
.w
[2];
3304 fstar
.w
[1] = P256
.w
[1];
3305 fstar
.w
[0] = P256
.w
[0];
3307 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
3308 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
3309 // if (0 < f* < 10^(-x)) then the result is a midpoint
3310 // if floor(C*) is even then C* = floor(C*) - logical right
3311 // shift; C* has p decimal digits, correct by Prop. 1)
3312 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
3313 // shift; C* has p decimal digits, correct by Pr. 1)
3315 // C* = floor(C*) (logical right shift; C has p decimal digits,
3316 // correct by Property 1)
3317 // n = C* * 10^(e+x)
3319 // shift right C* by Ex-128 = shiftright128[ind]
3320 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
3321 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
3323 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
3324 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
3325 } else { // 22 <= ind - 1 <= 33
3326 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
3328 // determine inexactness of the rounding of C*
3329 // if (0 < f* - 1/2 < 10^(-x)) then
3330 // the result is exact
3331 // else // if (f* - 1/2 > T*) then
3332 // the result is inexact
3334 if (fstar
.w
[1] > 0x8000000000000000ull
||
3335 (fstar
.w
[1] == 0x8000000000000000ull
3336 && fstar
.w
[0] > 0x0ull
)) {
3337 // f* > 1/2 and the result may be exact
3338 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
3339 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
3340 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
3341 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
3342 // set the inexact flag
3343 *pfpsf
|= INEXACT_EXCEPTION
;
3344 } // else the result is exact
3345 } else { // the result is inexact; f2* <= 1/2
3346 // set the inexact flag
3347 *pfpsf
|= INEXACT_EXCEPTION
;
3349 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
3350 if (fstar
.w
[3] > 0x0 ||
3351 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
3352 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
3353 (fstar
.w
[1] || fstar
.w
[0]))) {
3354 // f2* > 1/2 and the result may be exact
3355 // Calculate f2* - 1/2
3356 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
3357 tmp64A
= fstar
.w
[3];
3358 if (tmp64
> fstar
.w
[2])
3361 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
3362 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
3363 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
3364 // set the inexact flag
3365 *pfpsf
|= INEXACT_EXCEPTION
;
3366 } // else the result is exact
3367 } else { // the result is inexact; f2* <= 1/2
3368 // set the inexact flag
3369 *pfpsf
|= INEXACT_EXCEPTION
;
3371 } else { // if 22 <= ind <= 33
3372 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
3373 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
3374 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
3375 // f2* > 1/2 and the result may be exact
3376 // Calculate f2* - 1/2
3377 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
3378 if (tmp64
|| fstar
.w
[2]
3379 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
3380 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
3381 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
3382 // set the inexact flag
3383 *pfpsf
|= INEXACT_EXCEPTION
;
3384 } // else the result is exact
3385 } else { // the result is inexact; f2* <= 1/2
3386 // set the inexact flag
3387 *pfpsf
|= INEXACT_EXCEPTION
;
3391 // if the result was a midpoint it was rounded away from zero
3392 res
= Cstar
.w
[0]; // the result is positive
3393 } else if (exp
== 0) {
3394 // 1 <= q <= 20, but x < 2^64 - 1/2 so in this case C1.w[1] has to be 0
3398 // if (exp > 0) => 1 <= exp <= 19, 1 <= q < 19, 2 <= q + exp <= 20
3399 // res = C * 10^exp (exact) - must fit in 64 bits
3400 res
= C1
.w
[0] * ten2k64
[exp
];