1 /* Copyright (C) 2007, 2009 Free Software Foundation, Inc.
3 This file is part of GCC.
5 GCC is free software; you can redistribute it and/or modify it under
6 the terms of the GNU General Public License as published by the Free
7 Software Foundation; either version 3, or (at your option) any later
10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11 WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 Under Section 7 of GPL version 3, you are granted additional
16 permissions described in the GCC Runtime Library Exception, version
17 3.1, as published by the Free Software Foundation.
19 You should have received a copy of the GNU General Public License and
20 a copy of the GCC Runtime Library Exception along with this program;
21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
22 <http://www.gnu.org/licenses/>. */
24 #include "bid_internal.h"
26 /*****************************************************************************
27 * BID128_to_uint32_rnint
28 ****************************************************************************/
30 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
31 bid128_to_uint32_rnint
, x
)
36 int exp
; // unbiased exponent
37 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
40 unsigned int x_nr_bits
;
43 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
48 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
49 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
50 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
53 // check for NaN or Infinity
54 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
56 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
57 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
59 *pfpsf
|= INVALID_EXCEPTION
;
60 // return Integer Indefinite
64 *pfpsf
|= INVALID_EXCEPTION
;
65 // return Integer Indefinite
69 } else { // x is not a NaN, so it must be infinity
70 if (!x_sign
) { // x is +inf
72 *pfpsf
|= INVALID_EXCEPTION
;
73 // return Integer Indefinite
77 *pfpsf
|= INVALID_EXCEPTION
;
78 // return Integer Indefinite
84 // check for non-canonical values (after the check for special values)
85 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
86 || (C1
.w
[1] == 0x0001ed09bead87c0ull
87 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
88 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
91 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
95 } else { // x is not special and is not zero
97 // q = nr. of decimal digits in x
98 // determine first the nr. of bits in x
100 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
101 // split the 64-bit value in two 32-bit halves to avoid rounding errors
102 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
103 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
105 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
107 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
109 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
111 } else { // if x < 2^53
112 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
114 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
116 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
117 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
119 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
121 q
= nr_digits
[x_nr_bits
- 1].digits
;
123 q
= nr_digits
[x_nr_bits
- 1].digits1
;
124 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
125 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
126 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
129 exp
= (x_exp
>> 49) - 6176;
130 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
132 *pfpsf
|= INVALID_EXCEPTION
;
133 // return Integer Indefinite
136 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
137 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
138 // so x rounded to an integer may or may not fit in an unsigned 32-bit int
139 // the cases that do not fit are identified here; the ones that fit
140 // fall through and will be handled with other cases further,
141 // under '1 <= q + exp <= 10'
142 if (x_sign
) { // if n < 0 and q + exp = 10
143 // if n < -1/2 then n cannot be converted to uint32 with RN
144 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 1/2
145 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x05, 1<=q<=34
147 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
148 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
149 if (tmp64
> 0x05ull
) {
151 *pfpsf
|= INVALID_EXCEPTION
;
152 // return Integer Indefinite
156 // else cases that can be rounded to a 32-bit int fall through
157 // to '1 <= q + exp <= 10'
158 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
159 // 0.c(0)c(1)...c(q-1) * 10^11 > 0x05 <=>
160 // C > 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
163 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
164 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
165 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
166 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
168 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] > C
.w
[0])) {
170 *pfpsf
|= INVALID_EXCEPTION
;
171 // return Integer Indefinite
175 // else cases that can be rounded to a 32-bit int fall through
176 // to '1 <= q + exp <= 10'
178 } else { // if n > 0 and q + exp = 10
179 // if n >= 2^32 - 1/2 then n is too large
180 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
181 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
183 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
184 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
185 if (tmp64
>= 0x9fffffffbull
) {
187 *pfpsf
|= INVALID_EXCEPTION
;
188 // return Integer Indefinite
192 // else cases that can be rounded to a 32-bit int fall through
193 // to '1 <= q + exp <= 10'
194 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
195 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
196 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
197 // (scale 2^32-1/2 up)
198 tmp64
= 0x9fffffffbull
;
199 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
200 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
201 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
202 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
205 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
207 *pfpsf
|= INVALID_EXCEPTION
;
208 // return Integer Indefinite
212 // else cases that can be rounded to a 32-bit int fall through
213 // to '1 <= q + exp <= 10'
217 // n is not too large to be converted to int32: -1/2 <= n < 2^32 - 1/2
218 // Note: some of the cases tested for above fall through to this point
219 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
223 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
224 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
231 if (ind
<= 18) { // 0 <= ind <= 18
232 if ((C1
.w
[1] == 0) && (C1
.w
[0] <= midpoint64
[ind
])) {
233 res
= 0x00000000; // return 0
234 } else if (!x_sign
) { // n > 0
235 res
= 0x00000001; // return +1
238 *pfpsf
|= INVALID_EXCEPTION
;
240 } else { // 19 <= ind <= 33
241 if ((C1
.w
[1] < midpoint128
[ind
- 19].w
[1])
242 || ((C1
.w
[1] == midpoint128
[ind
- 19].w
[1])
243 && (C1
.w
[0] <= midpoint128
[ind
- 19].w
[0]))) {
244 res
= 0x00000000; // return 0
245 } else if (!x_sign
) { // n > 0
246 res
= 0x00000001; // return +1
249 *pfpsf
|= INVALID_EXCEPTION
;
252 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
253 if (x_sign
) { // x <= -1
255 *pfpsf
|= INVALID_EXCEPTION
;
256 // return Integer Indefinite
260 // 1 <= x < 2^32-1/2 so x can be rounded
261 // to nearest to a 32-bit unsigned integer
262 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
263 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
264 // chop off ind digits from the lower part of C1
265 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
268 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
270 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
271 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
275 // calculate C* and f*
276 // C* is actually floor(C*) in this case
277 // C* and f* need shifting and masking, as shown by
278 // shiftright128[] and maskhigh128[]
280 // kx = 10^(-x) = ten2mk128[ind - 1]
281 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
282 // the approximation of 10^(-x) was rounded up to 118 bits
283 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
284 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
285 Cstar
.w
[1] = P256
.w
[3];
286 Cstar
.w
[0] = P256
.w
[2];
288 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
289 fstar
.w
[1] = P256
.w
[1];
290 fstar
.w
[0] = P256
.w
[0];
291 } else { // 22 <= ind - 1 <= 33
293 Cstar
.w
[0] = P256
.w
[3];
294 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
295 fstar
.w
[2] = P256
.w
[2];
296 fstar
.w
[1] = P256
.w
[1];
297 fstar
.w
[0] = P256
.w
[0];
299 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
300 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
301 // if (0 < f* < 10^(-x)) then the result is a midpoint
302 // if floor(C*) is even then C* = floor(C*) - logical right
303 // shift; C* has p decimal digits, correct by Prop. 1)
304 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
305 // shift; C* has p decimal digits, correct by Pr. 1)
307 // C* = floor(C*) (logical right shift; C has p decimal digits,
308 // correct by Property 1)
311 // shift right C* by Ex-128 = shiftright128[ind]
312 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
313 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
315 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
316 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
317 } else { // 22 <= ind - 1 <= 33
318 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
320 // if the result was a midpoint it was rounded away from zero, so
321 // it will need a correction
322 // check for midpoints
323 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
324 && (fstar
.w
[1] || fstar
.w
[0])
325 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
326 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
327 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
328 // the result is a midpoint; round to nearest
329 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
330 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
331 Cstar
.w
[0]--; // Cstar.w[0] is now even
332 } // else MP in [ODD, EVEN]
334 res
= Cstar
.w
[0]; // the result is positive
335 } else if (exp
== 0) {
339 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
340 // res = C * 10^exp (exact)
341 res
= C1
.w
[0] * ten2k64
[exp
];
349 /*****************************************************************************
350 * BID128_to_uint32_xrnint
351 ****************************************************************************/
353 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
354 bid128_to_uint32_xrnint
, x
)
359 int exp
; // unbiased exponent
360 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
361 UINT64 tmp64
, tmp64A
;
363 unsigned int x_nr_bits
;
366 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
369 unsigned int tmp_inexact
= 0;
372 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
373 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
374 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
377 // check for NaN or Infinity
378 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
380 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
381 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
383 *pfpsf
|= INVALID_EXCEPTION
;
384 // return Integer Indefinite
386 } else { // x is QNaN
388 *pfpsf
|= INVALID_EXCEPTION
;
389 // return Integer Indefinite
393 } else { // x is not a NaN, so it must be infinity
394 if (!x_sign
) { // x is +inf
396 *pfpsf
|= INVALID_EXCEPTION
;
397 // return Integer Indefinite
399 } else { // x is -inf
401 *pfpsf
|= INVALID_EXCEPTION
;
402 // return Integer Indefinite
408 // check for non-canonical values (after the check for special values)
409 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
410 || (C1
.w
[1] == 0x0001ed09bead87c0ull
411 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
412 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
415 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
419 } else { // x is not special and is not zero
421 // q = nr. of decimal digits in x
422 // determine first the nr. of bits in x
424 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
425 // split the 64-bit value in two 32-bit halves to avoid rounding errors
426 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
427 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
429 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
431 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
433 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
435 } else { // if x < 2^53
436 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
438 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
440 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
441 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
443 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
445 q
= nr_digits
[x_nr_bits
- 1].digits
;
447 q
= nr_digits
[x_nr_bits
- 1].digits1
;
448 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
449 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
450 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
453 exp
= (x_exp
>> 49) - 6176;
454 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
456 *pfpsf
|= INVALID_EXCEPTION
;
457 // return Integer Indefinite
460 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
461 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
462 // so x rounded to an integer may or may not fit in an unsigned 32-bit int
463 // the cases that do not fit are identified here; the ones that fit
464 // fall through and will be handled with other cases further,
465 // under '1 <= q + exp <= 10'
466 if (x_sign
) { // if n < 0 and q + exp = 10
467 // if n < -1/2 then n cannot be converted to uint32 with RN
468 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 1/2
469 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x05, 1<=q<=34
471 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
472 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
473 if (tmp64
> 0x05ull
) {
475 *pfpsf
|= INVALID_EXCEPTION
;
476 // return Integer Indefinite
480 // else cases that can be rounded to a 32-bit int fall through
481 // to '1 <= q + exp <= 10'
482 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
483 // 0.c(0)c(1)...c(q-1) * 10^11 > 0x05 <=>
484 // C > 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
487 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
488 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
489 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
490 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
492 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] > C
.w
[0])) {
494 *pfpsf
|= INVALID_EXCEPTION
;
495 // return Integer Indefinite
499 // else cases that can be rounded to a 32-bit int fall through
500 // to '1 <= q + exp <= 10'
502 } else { // if n > 0 and q + exp = 10
503 // if n >= 2^32 - 1/2 then n is too large
504 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
505 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
507 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
508 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
509 if (tmp64
>= 0x9fffffffbull
) {
511 *pfpsf
|= INVALID_EXCEPTION
;
512 // return Integer Indefinite
516 // else cases that can be rounded to a 32-bit int fall through
517 // to '1 <= q + exp <= 10'
518 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
519 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
520 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
521 // (scale 2^32-1/2 up)
522 tmp64
= 0x9fffffffbull
;
523 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
524 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
525 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
526 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
529 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
531 *pfpsf
|= INVALID_EXCEPTION
;
532 // return Integer Indefinite
536 // else cases that can be rounded to a 32-bit int fall through
537 // to '1 <= q + exp <= 10'
541 // n is not too large to be converted to int32: -1/2 <= n < 2^32 - 1/2
542 // Note: some of the cases tested for above fall through to this point
543 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
545 *pfpsf
|= INEXACT_EXCEPTION
;
549 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
550 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
557 if (ind
<= 18) { // 0 <= ind <= 18
558 if ((C1
.w
[1] == 0) && (C1
.w
[0] <= midpoint64
[ind
])) {
559 res
= 0x00000000; // return 0
560 } else if (!x_sign
) { // n > 0
561 res
= 0x00000001; // return +1
564 *pfpsf
|= INVALID_EXCEPTION
;
567 } else { // 19 <= ind <= 33
568 if ((C1
.w
[1] < midpoint128
[ind
- 19].w
[1])
569 || ((C1
.w
[1] == midpoint128
[ind
- 19].w
[1])
570 && (C1
.w
[0] <= midpoint128
[ind
- 19].w
[0]))) {
571 res
= 0x00000000; // return 0
572 } else if (!x_sign
) { // n > 0
573 res
= 0x00000001; // return +1
576 *pfpsf
|= INVALID_EXCEPTION
;
581 *pfpsf
|= INEXACT_EXCEPTION
;
582 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
583 if (x_sign
) { // x <= -1
585 *pfpsf
|= INVALID_EXCEPTION
;
586 // return Integer Indefinite
590 // 1 <= x < 2^32-1/2 so x can be rounded
591 // to nearest to a 32-bit unsigned integer
592 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
593 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
594 // chop off ind digits from the lower part of C1
595 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
598 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
600 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
601 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
605 // calculate C* and f*
606 // C* is actually floor(C*) in this case
607 // C* and f* need shifting and masking, as shown by
608 // shiftright128[] and maskhigh128[]
610 // kx = 10^(-x) = ten2mk128[ind - 1]
611 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
612 // the approximation of 10^(-x) was rounded up to 118 bits
613 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
614 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
615 Cstar
.w
[1] = P256
.w
[3];
616 Cstar
.w
[0] = P256
.w
[2];
618 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
619 fstar
.w
[1] = P256
.w
[1];
620 fstar
.w
[0] = P256
.w
[0];
621 } else { // 22 <= ind - 1 <= 33
623 Cstar
.w
[0] = P256
.w
[3];
624 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
625 fstar
.w
[2] = P256
.w
[2];
626 fstar
.w
[1] = P256
.w
[1];
627 fstar
.w
[0] = P256
.w
[0];
629 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
630 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
631 // if (0 < f* < 10^(-x)) then the result is a midpoint
632 // if floor(C*) is even then C* = floor(C*) - logical right
633 // shift; C* has p decimal digits, correct by Prop. 1)
634 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
635 // shift; C* has p decimal digits, correct by Pr. 1)
637 // C* = floor(C*) (logical right shift; C has p decimal digits,
638 // correct by Property 1)
641 // shift right C* by Ex-128 = shiftright128[ind]
642 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
643 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
645 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
646 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
647 } else { // 22 <= ind - 1 <= 33
648 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
650 // determine inexactness of the rounding of C*
651 // if (0 < f* - 1/2 < 10^(-x)) then
652 // the result is exact
653 // else // if (f* - 1/2 > T*) then
654 // the result is inexact
656 if (fstar
.w
[1] > 0x8000000000000000ull
||
657 (fstar
.w
[1] == 0x8000000000000000ull
658 && fstar
.w
[0] > 0x0ull
)) {
659 // f* > 1/2 and the result may be exact
660 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
661 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
662 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
663 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
664 // set the inexact flag
665 // *pfpsf |= INEXACT_EXCEPTION;
667 } // else the result is exact
668 } else { // the result is inexact; f2* <= 1/2
669 // set the inexact flag
670 // *pfpsf |= INEXACT_EXCEPTION;
673 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
674 if (fstar
.w
[3] > 0x0 ||
675 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
676 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
677 (fstar
.w
[1] || fstar
.w
[0]))) {
678 // f2* > 1/2 and the result may be exact
679 // Calculate f2* - 1/2
680 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
682 if (tmp64
> fstar
.w
[2])
685 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
686 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
687 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
688 // set the inexact flag
689 // *pfpsf |= INEXACT_EXCEPTION;
691 } // else the result is exact
692 } else { // the result is inexact; f2* <= 1/2
693 // set the inexact flag
694 // *pfpsf |= INEXACT_EXCEPTION;
697 } else { // if 22 <= ind <= 33
698 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
699 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
700 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
701 // f2* > 1/2 and the result may be exact
702 // Calculate f2* - 1/2
703 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
704 if (tmp64
|| fstar
.w
[2]
705 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
706 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
707 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
708 // set the inexact flag
709 // *pfpsf |= INEXACT_EXCEPTION;
711 } // else the result is exact
712 } else { // the result is inexact; f2* <= 1/2
713 // set the inexact flag
714 // *pfpsf |= INEXACT_EXCEPTION;
719 // if the result was a midpoint it was rounded away from zero, so
720 // it will need a correction
721 // check for midpoints
722 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
723 && (fstar
.w
[1] || fstar
.w
[0])
724 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
725 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
726 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
727 // the result is a midpoint; round to nearest
728 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
729 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
730 Cstar
.w
[0]--; // Cstar.w[0] is now even
731 } // else MP in [ODD, EVEN]
733 res
= Cstar
.w
[0]; // the result is positive
735 *pfpsf
|= INEXACT_EXCEPTION
;
736 } else if (exp
== 0) {
740 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
741 // res = C * 10^exp (exact)
742 res
= C1
.w
[0] * ten2k64
[exp
];
750 /*****************************************************************************
751 * BID128_to_uint32_floor
752 ****************************************************************************/
754 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
755 bid128_to_uint32_floor
, x
)
760 int exp
; // unbiased exponent
761 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
762 UINT64 tmp64
, tmp64A
;
764 unsigned int x_nr_bits
;
767 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
770 int is_inexact_gt_midpoint
= 0;
771 int is_midpoint_lt_even
= 0;
774 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
775 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
776 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
779 // check for NaN or Infinity
780 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
782 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
783 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
785 *pfpsf
|= INVALID_EXCEPTION
;
786 // return Integer Indefinite
788 } else { // x is QNaN
790 *pfpsf
|= INVALID_EXCEPTION
;
791 // return Integer Indefinite
795 } else { // x is not a NaN, so it must be infinity
796 if (!x_sign
) { // x is +inf
798 *pfpsf
|= INVALID_EXCEPTION
;
799 // return Integer Indefinite
801 } else { // x is -inf
803 *pfpsf
|= INVALID_EXCEPTION
;
804 // return Integer Indefinite
810 // check for non-canonical values (after the check for special values)
811 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
812 || (C1
.w
[1] == 0x0001ed09bead87c0ull
813 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
814 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
817 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
821 } else { // x is not special and is not zero
825 *pfpsf
|= INVALID_EXCEPTION
;
826 // return Integer Indefinite
830 // x > 0 from this point on
831 // q = nr. of decimal digits in x
832 // determine first the nr. of bits in x
834 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
835 // split the 64-bit value in two 32-bit halves to avoid rounding errors
836 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
837 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
839 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
841 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
843 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
845 } else { // if x < 2^53
846 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
848 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
850 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
851 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
853 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
855 q
= nr_digits
[x_nr_bits
- 1].digits
;
857 q
= nr_digits
[x_nr_bits
- 1].digits1
;
858 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
859 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
860 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
863 exp
= (x_exp
>> 49) - 6176;
865 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
867 *pfpsf
|= INVALID_EXCEPTION
;
868 // return Integer Indefinite
871 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
872 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
873 // so x rounded to an integer may or may not fit in a signed 32-bit int
874 // the cases that do not fit are identified here; the ones that fit
875 // fall through and will be handled with other cases further,
876 // under '1 <= q + exp <= 10'
877 // n > 0 and q + exp = 10
878 // if n >= 2^32 then n is too large
879 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
880 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
882 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
883 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
884 if (tmp64
>= 0xa00000000ull
) {
886 *pfpsf
|= INVALID_EXCEPTION
;
887 // return Integer Indefinite
891 // else cases that can be rounded to a 32-bit int fall through
892 // to '1 <= q + exp <= 10'
893 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
894 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
895 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
897 tmp64
= 0xa00000000ull
;
898 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
899 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
900 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
901 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
903 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
905 *pfpsf
|= INVALID_EXCEPTION
;
906 // return Integer Indefinite
910 // else cases that can be rounded to a 32-bit int fall through
911 // to '1 <= q + exp <= 10'
914 // n is not too large to be converted to int32: 0 <= n < 2^31
915 // Note: some of the cases tested for above fall through to this point
916 if ((q
+ exp
) <= 0) {
917 // n = +0.0...c(0)c(1)...c(q-1) or n = +0.c(0)c(1)...c(q-1)
921 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
922 // 1 <= x < 2^32 so x can be rounded down to a 32-bit unsigned integer
923 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
924 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
925 // chop off ind digits from the lower part of C1
926 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
929 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
931 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
932 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
936 // calculate C* and f*
937 // C* is actually floor(C*) in this case
938 // C* and f* need shifting and masking, as shown by
939 // shiftright128[] and maskhigh128[]
941 // kx = 10^(-x) = ten2mk128[ind - 1]
942 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
943 // the approximation of 10^(-x) was rounded up to 118 bits
944 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
945 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
946 Cstar
.w
[1] = P256
.w
[3];
947 Cstar
.w
[0] = P256
.w
[2];
949 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
950 fstar
.w
[1] = P256
.w
[1];
951 fstar
.w
[0] = P256
.w
[0];
952 } else { // 22 <= ind - 1 <= 33
954 Cstar
.w
[0] = P256
.w
[3];
955 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
956 fstar
.w
[2] = P256
.w
[2];
957 fstar
.w
[1] = P256
.w
[1];
958 fstar
.w
[0] = P256
.w
[0];
960 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
961 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
962 // if (0 < f* < 10^(-x)) then the result is a midpoint
963 // if floor(C*) is even then C* = floor(C*) - logical right
964 // shift; C* has p decimal digits, correct by Prop. 1)
965 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
966 // shift; C* has p decimal digits, correct by Pr. 1)
968 // C* = floor(C*) (logical right shift; C has p decimal digits,
969 // correct by Property 1)
972 // shift right C* by Ex-128 = shiftright128[ind]
973 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
974 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
976 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
977 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
978 } else { // 22 <= ind - 1 <= 33
979 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
981 // determine inexactness of the rounding of C*
982 // if (0 < f* - 1/2 < 10^(-x)) then
983 // the result is exact
984 // else // if (f* - 1/2 > T*) then
985 // the result is inexact
987 if (fstar
.w
[1] > 0x8000000000000000ull
||
988 (fstar
.w
[1] == 0x8000000000000000ull
989 && fstar
.w
[0] > 0x0ull
)) {
990 // f* > 1/2 and the result may be exact
991 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
992 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
993 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
994 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
995 } // else the result is exact
996 } else { // the result is inexact; f2* <= 1/2
997 is_inexact_gt_midpoint
= 1;
999 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
1000 if (fstar
.w
[3] > 0x0 ||
1001 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
1002 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
1003 (fstar
.w
[1] || fstar
.w
[0]))) {
1004 // f2* > 1/2 and the result may be exact
1005 // Calculate f2* - 1/2
1006 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
1007 tmp64A
= fstar
.w
[3];
1008 if (tmp64
> fstar
.w
[2])
1011 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1012 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1013 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1014 } // else the result is exact
1015 } else { // the result is inexact; f2* <= 1/2
1016 is_inexact_gt_midpoint
= 1;
1018 } else { // if 22 <= ind <= 33
1019 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
1020 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
1021 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
1022 // f2* > 1/2 and the result may be exact
1023 // Calculate f2* - 1/2
1024 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
1025 if (tmp64
|| fstar
.w
[2]
1026 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1027 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1028 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1029 } // else the result is exact
1030 } else { // the result is inexact; f2* <= 1/2
1031 is_inexact_gt_midpoint
= 1;
1035 // if the result was a midpoint it was rounded away from zero, so
1036 // it will need a correction
1037 // check for midpoints
1038 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
1039 && (fstar
.w
[1] || fstar
.w
[0])
1040 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
1041 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1042 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
1043 // the result is a midpoint; round to nearest
1044 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
1045 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
1046 Cstar
.w
[0]--; // Cstar.w[0] is now even
1047 is_inexact_gt_midpoint
= 0;
1048 } else { // else MP in [ODD, EVEN]
1049 is_midpoint_lt_even
= 1;
1050 is_inexact_gt_midpoint
= 0;
1053 // general correction for RM
1054 if (is_midpoint_lt_even
|| is_inexact_gt_midpoint
) {
1055 Cstar
.w
[0] = Cstar
.w
[0] - 1;
1057 ; // the result is already correct
1060 } else if (exp
== 0) {
1064 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1065 // res = +C * 10^exp (exact)
1066 res
= C1
.w
[0] * ten2k64
[exp
];
1074 /*****************************************************************************
1075 * BID128_to_uint32_xfloor
1076 ****************************************************************************/
1078 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
1079 bid128_to_uint32_xfloor
, x
)
1084 int exp
; // unbiased exponent
1085 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1086 UINT64 tmp64
, tmp64A
;
1087 BID_UI64DOUBLE tmp1
;
1088 unsigned int x_nr_bits
;
1091 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
1094 int is_inexact_gt_midpoint
= 0;
1095 int is_midpoint_lt_even
= 0;
1098 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1099 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
1100 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
1103 // check for NaN or Infinity
1104 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
1106 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
1107 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
1109 *pfpsf
|= INVALID_EXCEPTION
;
1110 // return Integer Indefinite
1112 } else { // x is QNaN
1114 *pfpsf
|= INVALID_EXCEPTION
;
1115 // return Integer Indefinite
1119 } else { // x is not a NaN, so it must be infinity
1120 if (!x_sign
) { // x is +inf
1122 *pfpsf
|= INVALID_EXCEPTION
;
1123 // return Integer Indefinite
1125 } else { // x is -inf
1127 *pfpsf
|= INVALID_EXCEPTION
;
1128 // return Integer Indefinite
1134 // check for non-canonical values (after the check for special values)
1135 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
1136 || (C1
.w
[1] == 0x0001ed09bead87c0ull
1137 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
1138 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
1141 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
1145 } else { // x is not special and is not zero
1149 *pfpsf
|= INVALID_EXCEPTION
;
1150 // return Integer Indefinite
1154 // x > 0 from this point on
1155 // q = nr. of decimal digits in x
1156 // determine first the nr. of bits in x
1158 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
1159 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1160 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
1161 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
1163 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1164 } else { // x < 2^32
1165 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
1167 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1169 } else { // if x < 2^53
1170 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
1172 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1174 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1175 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
1177 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1179 q
= nr_digits
[x_nr_bits
- 1].digits
;
1181 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1182 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
1183 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
1184 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
1187 exp
= (x_exp
>> 49) - 6176;
1188 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1190 *pfpsf
|= INVALID_EXCEPTION
;
1191 // return Integer Indefinite
1194 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1195 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1196 // so x rounded to an integer may or may not fit in a signed 32-bit int
1197 // the cases that do not fit are identified here; the ones that fit
1198 // fall through and will be handled with other cases further,
1199 // under '1 <= q + exp <= 10'
1200 // n > 0 and q + exp = 10
1201 // if n >= 2^32 then n is too large
1202 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
1203 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
1205 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
1206 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1207 if (tmp64
>= 0xa00000000ull
) {
1209 *pfpsf
|= INVALID_EXCEPTION
;
1210 // return Integer Indefinite
1214 // else cases that can be rounded to a 32-bit int fall through
1215 // to '1 <= q + exp <= 10'
1216 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1217 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
1218 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
1220 tmp64
= 0xa00000000ull
;
1221 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1222 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
1223 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1224 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
1226 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
1228 *pfpsf
|= INVALID_EXCEPTION
;
1229 // return Integer Indefinite
1233 // else cases that can be rounded to a 32-bit int fall through
1234 // to '1 <= q + exp <= 10'
1237 // n is not too large to be converted to int32: 0 <= n < 2^31
1238 // Note: some of the cases tested for above fall through to this point
1239 if ((q
+ exp
) <= 0) {
1240 // n = +0.0...c(0)c(1)...c(q-1) or n = +0.c(0)c(1)...c(q-1)
1242 *pfpsf
|= INEXACT_EXCEPTION
;
1246 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
1247 // 1 <= x < 2^32 so x can be rounded down to a 32-bit unsigned integer
1248 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
1249 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
1250 // chop off ind digits from the lower part of C1
1251 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
1254 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
1256 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
1257 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
1259 if (C1
.w
[0] < tmp64
)
1261 // calculate C* and f*
1262 // C* is actually floor(C*) in this case
1263 // C* and f* need shifting and masking, as shown by
1264 // shiftright128[] and maskhigh128[]
1266 // kx = 10^(-x) = ten2mk128[ind - 1]
1267 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1268 // the approximation of 10^(-x) was rounded up to 118 bits
1269 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
1270 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1271 Cstar
.w
[1] = P256
.w
[3];
1272 Cstar
.w
[0] = P256
.w
[2];
1274 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
1275 fstar
.w
[1] = P256
.w
[1];
1276 fstar
.w
[0] = P256
.w
[0];
1277 } else { // 22 <= ind - 1 <= 33
1279 Cstar
.w
[0] = P256
.w
[3];
1280 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
1281 fstar
.w
[2] = P256
.w
[2];
1282 fstar
.w
[1] = P256
.w
[1];
1283 fstar
.w
[0] = P256
.w
[0];
1285 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1286 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1287 // if (0 < f* < 10^(-x)) then the result is a midpoint
1288 // if floor(C*) is even then C* = floor(C*) - logical right
1289 // shift; C* has p decimal digits, correct by Prop. 1)
1290 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
1291 // shift; C* has p decimal digits, correct by Pr. 1)
1293 // C* = floor(C*) (logical right shift; C has p decimal digits,
1294 // correct by Property 1)
1295 // n = C* * 10^(e+x)
1297 // shift right C* by Ex-128 = shiftright128[ind]
1298 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
1299 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1301 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
1302 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1303 } else { // 22 <= ind - 1 <= 33
1304 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
1306 // determine inexactness of the rounding of C*
1307 // if (0 < f* - 1/2 < 10^(-x)) then
1308 // the result is exact
1309 // else // if (f* - 1/2 > T*) then
1310 // the result is inexact
1312 if (fstar
.w
[1] > 0x8000000000000000ull
||
1313 (fstar
.w
[1] == 0x8000000000000000ull
1314 && fstar
.w
[0] > 0x0ull
)) {
1315 // f* > 1/2 and the result may be exact
1316 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
1317 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
1318 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
1319 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
1320 // set the inexact flag
1321 *pfpsf
|= INEXACT_EXCEPTION
;
1322 } // else the result is exact
1323 } else { // the result is inexact; f2* <= 1/2
1324 // set the inexact flag
1325 *pfpsf
|= INEXACT_EXCEPTION
;
1326 is_inexact_gt_midpoint
= 1;
1328 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
1329 if (fstar
.w
[3] > 0x0 ||
1330 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
1331 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
1332 (fstar
.w
[1] || fstar
.w
[0]))) {
1333 // f2* > 1/2 and the result may be exact
1334 // Calculate f2* - 1/2
1335 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
1336 tmp64A
= fstar
.w
[3];
1337 if (tmp64
> 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
1346 } else { // the result is inexact; f2* <= 1/2
1347 // set the inexact flag
1348 *pfpsf
|= INEXACT_EXCEPTION
;
1349 is_inexact_gt_midpoint
= 1;
1351 } else { // if 22 <= ind <= 33
1352 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
1353 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
1354 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
1355 // f2* > 1/2 and the result may be exact
1356 // Calculate f2* - 1/2
1357 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
1358 if (tmp64
|| fstar
.w
[2]
1359 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1360 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1361 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1362 // set the inexact flag
1363 *pfpsf
|= INEXACT_EXCEPTION
;
1364 } // else the result is exact
1365 } else { // the result is inexact; f2* <= 1/2
1366 // set the inexact flag
1367 *pfpsf
|= INEXACT_EXCEPTION
;
1368 is_inexact_gt_midpoint
= 1;
1372 // if the result was a midpoint it was rounded away from zero, so
1373 // it will need a correction
1374 // check for midpoints
1375 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
1376 && (fstar
.w
[1] || fstar
.w
[0])
1377 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
1378 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1379 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
1380 // the result is a midpoint; round to nearest
1381 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
1382 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
1383 Cstar
.w
[0]--; // Cstar.w[0] is now even
1384 is_inexact_gt_midpoint
= 0;
1385 } else { // else MP in [ODD, EVEN]
1386 is_midpoint_lt_even
= 1;
1387 is_inexact_gt_midpoint
= 0;
1390 // general correction for RM
1391 if (is_midpoint_lt_even
|| is_inexact_gt_midpoint
) {
1392 Cstar
.w
[0] = Cstar
.w
[0] - 1;
1394 ; // the result is already correct
1397 } else if (exp
== 0) {
1401 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1402 // res = +C * 10^exp (exact)
1403 res
= C1
.w
[0] * ten2k64
[exp
];
1411 /*****************************************************************************
1412 * BID128_to_uint32_ceil
1413 ****************************************************************************/
1415 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
1416 bid128_to_uint32_ceil
, x
)
1421 int exp
; // unbiased exponent
1422 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1423 UINT64 tmp64
, tmp64A
;
1424 BID_UI64DOUBLE tmp1
;
1425 unsigned int x_nr_bits
;
1428 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
1431 int is_inexact_lt_midpoint
= 0;
1432 int is_midpoint_gt_even
= 0;
1435 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1436 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
1437 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
1440 // check for NaN or Infinity
1441 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
1443 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
1444 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
1446 *pfpsf
|= INVALID_EXCEPTION
;
1447 // return Integer Indefinite
1449 } else { // x is QNaN
1451 *pfpsf
|= INVALID_EXCEPTION
;
1452 // return Integer Indefinite
1456 } else { // x is not a NaN, so it must be infinity
1457 if (!x_sign
) { // x is +inf
1459 *pfpsf
|= INVALID_EXCEPTION
;
1460 // return Integer Indefinite
1462 } else { // x is -inf
1464 *pfpsf
|= INVALID_EXCEPTION
;
1465 // return Integer Indefinite
1471 // check for non-canonical values (after the check for special values)
1472 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
1473 || (C1
.w
[1] == 0x0001ed09bead87c0ull
1474 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
1475 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
1478 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
1482 } else { // x is not special and is not zero
1484 // q = nr. of decimal digits in x
1485 // determine first the nr. of bits in x
1487 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
1488 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1489 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
1490 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
1492 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1493 } else { // x < 2^32
1494 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
1496 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1498 } else { // if x < 2^53
1499 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
1501 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1503 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1504 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
1506 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1508 q
= nr_digits
[x_nr_bits
- 1].digits
;
1510 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1511 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
1512 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
1513 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
1516 exp
= (x_exp
>> 49) - 6176;
1518 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1520 *pfpsf
|= INVALID_EXCEPTION
;
1521 // return Integer Indefinite
1524 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1525 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1526 // so x rounded to an integer may or may not fit in a signed 32-bit int
1527 // the cases that do not fit are identified here; the ones that fit
1528 // fall through and will be handled with other cases further,
1529 // under '1 <= q + exp <= 10'
1530 if (x_sign
) { // if n < 0 and q + exp = 10
1531 // if n <= -1 then n is too large
1532 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
1533 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34
1535 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
1536 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1537 if (tmp64
>= 0x0aull
) {
1539 *pfpsf
|= INVALID_EXCEPTION
;
1540 // return Integer Indefinite
1544 // else cases that can be rounded to a 32-bit int fall through
1545 // to '1 <= q + exp <= 10'
1546 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1547 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
1548 // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
1551 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1552 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
1553 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1554 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
1556 if (C1
.w
[1] > C
.w
[1]
1557 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
1559 *pfpsf
|= INVALID_EXCEPTION
;
1560 // return Integer Indefinite
1564 // else cases that can be rounded to a 32-bit int fall through
1565 // to '1 <= q + exp <= 10'
1567 } else { // if n > 0 and q + exp = 10
1568 // if n > 2^32 - 1 then n is too large
1569 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^32 - 1
1570 // too large if 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6, 1<=q<=34
1572 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
1573 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1574 if (tmp64
> 0x9fffffff6ull
) {
1576 *pfpsf
|= INVALID_EXCEPTION
;
1577 // return Integer Indefinite
1581 // else cases that can be rounded to a 32-bit int fall through
1582 // to '1 <= q + exp <= 10'
1583 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1584 // 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6 <=>
1585 // C > 0x9fffffff6 * 10^(q-11) where 1 <= q - 11 <= 23
1587 tmp64
= 0x9fffffff6ull
;
1588 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1589 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
1590 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1591 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
1593 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] > C
.w
[0])) {
1595 *pfpsf
|= INVALID_EXCEPTION
;
1596 // return Integer Indefinite
1600 // else cases that can be rounded to a 32-bit int fall through
1601 // to '1 <= q + exp <= 10'
1605 // n is not too large to be converted to int32: -2^32-1 < n <= 2^32-1
1606 // Note: some of the cases tested for above fall through to this point
1607 if ((q
+ exp
) <= 0) {
1608 // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
1615 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
1616 // -2^32-1 < x <= -1 or 1 <= x <= 2^32-1 so x can be rounded
1617 // toward positive infinity to a 32-bit signed integer
1618 if (x_sign
) { // x <= -1 is invalid
1620 *pfpsf
|= INVALID_EXCEPTION
;
1621 // return Integer Indefinite
1625 // x > 0 from this point on
1626 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
1627 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
1628 // chop off ind digits from the lower part of C1
1629 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
1632 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
1634 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
1635 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
1637 if (C1
.w
[0] < tmp64
)
1639 // calculate C* and f*
1640 // C* is actually floor(C*) in this case
1641 // C* and f* need shifting and masking, as shown by
1642 // shiftright128[] and maskhigh128[]
1644 // kx = 10^(-x) = ten2mk128[ind - 1]
1645 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1646 // the approximation of 10^(-x) was rounded up to 118 bits
1647 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
1648 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1649 Cstar
.w
[1] = P256
.w
[3];
1650 Cstar
.w
[0] = P256
.w
[2];
1652 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
1653 fstar
.w
[1] = P256
.w
[1];
1654 fstar
.w
[0] = P256
.w
[0];
1655 } else { // 22 <= ind - 1 <= 33
1657 Cstar
.w
[0] = P256
.w
[3];
1658 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
1659 fstar
.w
[2] = P256
.w
[2];
1660 fstar
.w
[1] = P256
.w
[1];
1661 fstar
.w
[0] = P256
.w
[0];
1663 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1664 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1665 // if (0 < f* < 10^(-x)) then the result is a midpoint
1666 // if floor(C*) is even then C* = floor(C*) - logical right
1667 // shift; C* has p decimal digits, correct by Prop. 1)
1668 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
1669 // shift; C* has p decimal digits, correct by Pr. 1)
1671 // C* = floor(C*) (logical right shift; C has p decimal digits,
1672 // correct by Property 1)
1673 // n = C* * 10^(e+x)
1675 // shift right C* by Ex-128 = shiftright128[ind]
1676 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
1677 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1679 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
1680 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1681 } else { // 22 <= ind - 1 <= 33
1682 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
1684 // determine inexactness of the rounding of C*
1685 // if (0 < f* - 1/2 < 10^(-x)) then
1686 // the result is exact
1687 // else // if (f* - 1/2 > T*) then
1688 // the result is inexact
1690 if (fstar
.w
[1] > 0x8000000000000000ull
||
1691 (fstar
.w
[1] == 0x8000000000000000ull
1692 && fstar
.w
[0] > 0x0ull
)) {
1693 // f* > 1/2 and the result may be exact
1694 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
1695 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
1696 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
1697 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
1698 is_inexact_lt_midpoint
= 1;
1699 } // else the result is exact
1700 } else { // the result is inexact; f2* <= 1/2
1703 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
1704 if (fstar
.w
[3] > 0x0 ||
1705 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
1706 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
1707 (fstar
.w
[1] || fstar
.w
[0]))) {
1708 // f2* > 1/2 and the result may be exact
1709 // Calculate f2* - 1/2
1710 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
1711 tmp64A
= fstar
.w
[3];
1712 if (tmp64
> fstar
.w
[2])
1715 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1716 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1717 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1718 is_inexact_lt_midpoint
= 1;
1719 } // else the result is exact
1720 } else { // the result is inexact; f2* <= 1/2
1722 } else { // if 22 <= ind <= 33
1723 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
1724 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
1725 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
1726 // f2* > 1/2 and the result may be exact
1727 // Calculate f2* - 1/2
1728 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
1729 if (tmp64
|| fstar
.w
[2]
1730 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1731 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1732 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1733 is_inexact_lt_midpoint
= 1;
1734 } // else the result is exact
1735 } else { // the result is inexact; f2* <= 1/2
1740 // if the result was a midpoint it was rounded away from zero, so
1741 // it will need a correction
1742 // check for midpoints
1743 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
1744 && (fstar
.w
[1] || fstar
.w
[0])
1745 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
1746 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1747 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
1748 // the result is a midpoint; round to nearest
1749 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
1750 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
1751 Cstar
.w
[0]--; // Cstar.w[0] is now even
1752 is_midpoint_gt_even
= 1;
1753 is_inexact_lt_midpoint
= 0;
1754 } else { // else MP in [ODD, EVEN]
1755 is_inexact_lt_midpoint
= 0;
1758 // general correction for RM
1759 if (is_midpoint_gt_even
|| is_inexact_lt_midpoint
) {
1760 Cstar
.w
[0] = Cstar
.w
[0] + 1;
1762 ; // the result is already correct
1765 } else if (exp
== 0) {
1769 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1770 // res = +C * 10^exp (exact)
1771 res
= C1
.w
[0] * ten2k64
[exp
];
1779 /*****************************************************************************
1780 * BID128_to_uint32_xceil
1781 ****************************************************************************/
1783 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
1784 bid128_to_uint32_xceil
, x
)
1789 int exp
; // unbiased exponent
1790 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1791 UINT64 tmp64
, tmp64A
;
1792 BID_UI64DOUBLE tmp1
;
1793 unsigned int x_nr_bits
;
1796 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
1799 int is_inexact_lt_midpoint
= 0;
1800 int is_midpoint_gt_even
= 0;
1803 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1804 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
1805 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
1808 // check for NaN or Infinity
1809 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
1811 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
1812 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
1814 *pfpsf
|= INVALID_EXCEPTION
;
1815 // return Integer Indefinite
1817 } else { // x is QNaN
1819 *pfpsf
|= INVALID_EXCEPTION
;
1820 // return Integer Indefinite
1824 } else { // x is not a NaN, so it must be infinity
1825 if (!x_sign
) { // x is +inf
1827 *pfpsf
|= INVALID_EXCEPTION
;
1828 // return Integer Indefinite
1830 } else { // x is -inf
1832 *pfpsf
|= INVALID_EXCEPTION
;
1833 // return Integer Indefinite
1839 // check for non-canonical values (after the check for special values)
1840 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
1841 || (C1
.w
[1] == 0x0001ed09bead87c0ull
1842 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
1843 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
1846 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
1850 } else { // x is not special and is not zero
1852 // q = nr. of decimal digits in x
1853 // determine first the nr. of bits in x
1855 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
1856 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1857 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
1858 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
1860 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1861 } else { // x < 2^32
1862 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
1864 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1866 } else { // if x < 2^53
1867 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
1869 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1871 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1872 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
1874 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1876 q
= nr_digits
[x_nr_bits
- 1].digits
;
1878 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1879 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
1880 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
1881 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
1884 exp
= (x_exp
>> 49) - 6176;
1885 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1887 *pfpsf
|= INVALID_EXCEPTION
;
1888 // return Integer Indefinite
1891 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1892 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1893 // so x rounded to an integer may or may not fit in a signed 32-bit int
1894 // the cases that do not fit are identified here; the ones that fit
1895 // fall through and will be handled with other cases further,
1896 // under '1 <= q + exp <= 10'
1897 if (x_sign
) { // if n < 0 and q + exp = 10
1898 // if n <= -1 then n is too large
1899 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
1900 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34
1902 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
1903 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1904 if (tmp64
>= 0x0aull
) {
1906 *pfpsf
|= INVALID_EXCEPTION
;
1907 // return Integer Indefinite
1911 // else cases that can be rounded to a 32-bit int fall through
1912 // to '1 <= q + exp <= 10'
1913 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1914 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
1915 // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
1918 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1919 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
1920 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1921 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
1923 if (C1
.w
[1] > C
.w
[1]
1924 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
1926 *pfpsf
|= INVALID_EXCEPTION
;
1927 // return Integer Indefinite
1931 // else cases that can be rounded to a 32-bit int fall through
1932 // to '1 <= q + exp <= 10'
1934 } else { // if n > 0 and q + exp = 10
1935 // if n > 2^32 - 1 then n is too large
1936 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^32 - 1
1937 // too large if 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6, 1<=q<=34
1939 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
1940 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1941 if (tmp64
> 0x9fffffff6ull
) {
1943 *pfpsf
|= INVALID_EXCEPTION
;
1944 // return Integer Indefinite
1948 // else cases that can be rounded to a 32-bit int fall through
1949 // to '1 <= q + exp <= 10'
1950 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1951 // 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6 <=>
1952 // C > 0x9fffffff6 * 10^(q-11) where 1 <= q - 11 <= 23
1954 tmp64
= 0x9fffffff6ull
;
1955 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1956 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
1957 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1958 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
1960 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] > C
.w
[0])) {
1962 *pfpsf
|= INVALID_EXCEPTION
;
1963 // return Integer Indefinite
1967 // else cases that can be rounded to a 32-bit int fall through
1968 // to '1 <= q + exp <= 10'
1972 // n is not too large to be converted to int32: -2^32-1 < n <= 2^32-1
1973 // Note: some of the cases tested for above fall through to this point
1974 if ((q
+ exp
) <= 0) {
1975 // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
1977 *pfpsf
|= INEXACT_EXCEPTION
;
1984 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
1985 // -2^32-1 < x <= -1 or 1 <= x <= 2^32-1 so x can be rounded
1986 // toward positive infinity to a 32-bit signed integer
1987 if (x_sign
) { // x <= -1 is invalid
1989 *pfpsf
|= INVALID_EXCEPTION
;
1990 // return Integer Indefinite
1994 // x > 0 from this point on
1995 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
1996 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
1997 // chop off ind digits from the lower part of C1
1998 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
2001 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
2003 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
2004 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
2006 if (C1
.w
[0] < tmp64
)
2008 // calculate C* and f*
2009 // C* is actually floor(C*) in this case
2010 // C* and f* need shifting and masking, as shown by
2011 // shiftright128[] and maskhigh128[]
2013 // kx = 10^(-x) = ten2mk128[ind - 1]
2014 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2015 // the approximation of 10^(-x) was rounded up to 118 bits
2016 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
2017 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2018 Cstar
.w
[1] = P256
.w
[3];
2019 Cstar
.w
[0] = P256
.w
[2];
2021 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
2022 fstar
.w
[1] = P256
.w
[1];
2023 fstar
.w
[0] = P256
.w
[0];
2024 } else { // 22 <= ind - 1 <= 33
2026 Cstar
.w
[0] = P256
.w
[3];
2027 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
2028 fstar
.w
[2] = P256
.w
[2];
2029 fstar
.w
[1] = P256
.w
[1];
2030 fstar
.w
[0] = P256
.w
[0];
2032 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2033 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2034 // if (0 < f* < 10^(-x)) then the result is a midpoint
2035 // if floor(C*) is even then C* = floor(C*) - logical right
2036 // shift; C* has p decimal digits, correct by Prop. 1)
2037 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
2038 // shift; C* has p decimal digits, correct by Pr. 1)
2040 // C* = floor(C*) (logical right shift; C has p decimal digits,
2041 // correct by Property 1)
2042 // n = C* * 10^(e+x)
2044 // shift right C* by Ex-128 = shiftright128[ind]
2045 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
2046 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2048 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
2049 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2050 } else { // 22 <= ind - 1 <= 33
2051 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
2053 // determine inexactness of the rounding of C*
2054 // if (0 < f* - 1/2 < 10^(-x)) then
2055 // the result is exact
2056 // else // if (f* - 1/2 > T*) then
2057 // the result is inexact
2059 if (fstar
.w
[1] > 0x8000000000000000ull
||
2060 (fstar
.w
[1] == 0x8000000000000000ull
2061 && fstar
.w
[0] > 0x0ull
)) {
2062 // f* > 1/2 and the result may be exact
2063 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
2064 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
2065 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
2066 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
2067 // set the inexact flag
2068 *pfpsf
|= INEXACT_EXCEPTION
;
2069 is_inexact_lt_midpoint
= 1;
2070 } // else the result is exact
2071 } else { // the result is inexact; f2* <= 1/2
2072 // set the inexact flag
2073 *pfpsf
|= INEXACT_EXCEPTION
;
2075 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
2076 if (fstar
.w
[3] > 0x0 ||
2077 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
2078 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
2079 (fstar
.w
[1] || fstar
.w
[0]))) {
2080 // f2* > 1/2 and the result may be exact
2081 // Calculate f2* - 1/2
2082 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
2083 tmp64A
= fstar
.w
[3];
2084 if (tmp64
> fstar
.w
[2])
2087 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2088 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2089 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2090 // set the inexact flag
2091 *pfpsf
|= INEXACT_EXCEPTION
;
2092 is_inexact_lt_midpoint
= 1;
2093 } // else the result is exact
2094 } else { // the result is inexact; f2* <= 1/2
2095 // set the inexact flag
2096 *pfpsf
|= INEXACT_EXCEPTION
;
2098 } else { // if 22 <= ind <= 33
2099 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
2100 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
2101 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
2102 // f2* > 1/2 and the result may be exact
2103 // Calculate f2* - 1/2
2104 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
2105 if (tmp64
|| fstar
.w
[2]
2106 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2107 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2108 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2109 // set the inexact flag
2110 *pfpsf
|= INEXACT_EXCEPTION
;
2111 is_inexact_lt_midpoint
= 1;
2112 } // else the result is exact
2113 } else { // the result is inexact; f2* <= 1/2
2114 // set the inexact flag
2115 *pfpsf
|= INEXACT_EXCEPTION
;
2119 // if the result was a midpoint it was rounded away from zero, so
2120 // it will need a correction
2121 // check for midpoints
2122 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
2123 && (fstar
.w
[1] || fstar
.w
[0])
2124 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
2125 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2126 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
2127 // the result is a midpoint; round to nearest
2128 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
2129 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
2130 Cstar
.w
[0]--; // Cstar.w[0] is now even
2131 is_midpoint_gt_even
= 1;
2132 is_inexact_lt_midpoint
= 0;
2133 } else { // else MP in [ODD, EVEN]
2134 is_inexact_lt_midpoint
= 0;
2137 // general correction for RM
2138 if (is_midpoint_gt_even
|| is_inexact_lt_midpoint
) {
2139 Cstar
.w
[0] = Cstar
.w
[0] + 1;
2141 ; // the result is already correct
2144 } else if (exp
== 0) {
2148 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2149 // res = +C * 10^exp (exact)
2150 res
= C1
.w
[0] * ten2k64
[exp
];
2158 /*****************************************************************************
2159 * BID128_to_uint32_int
2160 ****************************************************************************/
2162 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
2163 bid128_to_uint32_int
, x
)
2168 int exp
; // unbiased exponent
2169 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2170 UINT64 tmp64
, tmp64A
;
2171 BID_UI64DOUBLE tmp1
;
2172 unsigned int x_nr_bits
;
2175 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
2178 int is_inexact_gt_midpoint
= 0;
2179 int is_midpoint_lt_even
= 0;
2182 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
2183 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
2184 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
2187 // check for NaN or Infinity
2188 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
2190 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
2191 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
2193 *pfpsf
|= INVALID_EXCEPTION
;
2194 // return Integer Indefinite
2196 } else { // x is QNaN
2198 *pfpsf
|= INVALID_EXCEPTION
;
2199 // return Integer Indefinite
2203 } else { // x is not a NaN, so it must be infinity
2204 if (!x_sign
) { // x is +inf
2206 *pfpsf
|= INVALID_EXCEPTION
;
2207 // return Integer Indefinite
2209 } else { // x is -inf
2211 *pfpsf
|= INVALID_EXCEPTION
;
2212 // return Integer Indefinite
2218 // check for non-canonical values (after the check for special values)
2219 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
2220 || (C1
.w
[1] == 0x0001ed09bead87c0ull
2221 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
2222 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
2225 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
2229 } else { // x is not special and is not zero
2231 // q = nr. of decimal digits in x
2232 // determine first the nr. of bits in x
2234 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
2235 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2236 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
2237 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
2239 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2240 } else { // x < 2^32
2241 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
2243 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2245 } else { // if x < 2^53
2246 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
2248 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2250 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2251 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
2253 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2255 q
= nr_digits
[x_nr_bits
- 1].digits
;
2257 q
= nr_digits
[x_nr_bits
- 1].digits1
;
2258 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
2259 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
2260 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
2263 exp
= (x_exp
>> 49) - 6176;
2264 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
2266 *pfpsf
|= INVALID_EXCEPTION
;
2267 // return Integer Indefinite
2270 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
2271 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
2272 // so x rounded to an integer may or may not fit in a signed 32-bit int
2273 // the cases that do not fit are identified here; the ones that fit
2274 // fall through and will be handled with other cases further,
2275 // under '1 <= q + exp <= 10'
2276 if (x_sign
) { // if n < 0 and q + exp = 10
2277 // if n <= -1 then n is too large
2278 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
2279 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a, 1<=q<=34
2281 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
2282 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2283 if (tmp64
>= 0x0aull
) {
2285 *pfpsf
|= INVALID_EXCEPTION
;
2286 // return Integer Indefinite
2290 // else cases that can be rounded to a 32-bit uint fall through
2291 // to '1 <= q + exp <= 10'
2292 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2293 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
2294 // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
2297 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2298 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
2299 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2300 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
2302 if (C1
.w
[1] > C
.w
[1]
2303 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
2305 *pfpsf
|= INVALID_EXCEPTION
;
2306 // return Integer Indefinite
2310 // else cases that can be rounded to a 32-bit int fall through
2311 // to '1 <= q + exp <= 10'
2313 } else { // if n > 0 and q + exp = 10
2314 // if n >= 2^32 then n is too large
2315 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
2316 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
2318 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
2319 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2320 if (tmp64
>= 0xa00000000ull
) {
2322 *pfpsf
|= INVALID_EXCEPTION
;
2323 // return Integer Indefinite
2327 // else cases that can be rounded to a 32-bit uint fall through
2328 // to '1 <= q + exp <= 10'
2329 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2330 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
2331 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
2333 tmp64
= 0xa00000000ull
;
2334 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2335 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
2336 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2337 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
2339 if (C1
.w
[1] > C
.w
[1]
2340 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
2342 *pfpsf
|= INVALID_EXCEPTION
;
2343 // return Integer Indefinite
2347 // else cases that can be rounded to a 32-bit int fall through
2348 // to '1 <= q + exp <= 10'
2352 // n is not too large to be converted to uint32: -2^32 < n < 2^32
2353 // Note: some of the cases tested for above fall through to this point
2354 if ((q
+ exp
) <= 0) {
2355 // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
2359 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
2360 // x = d(0)...d(k).d(k+1)..., k >= 0, d(0) != 0
2361 if (x_sign
) { // x <= -1
2363 *pfpsf
|= INVALID_EXCEPTION
;
2364 // return Integer Indefinite
2368 // x > 0 from this point on
2369 // 1 <= x < 2^32 so x can be rounded to zero to a 32-bit unsigned integer
2370 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
2371 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
2372 // chop off ind digits from the lower part of C1
2373 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
2376 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
2378 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
2379 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
2381 if (C1
.w
[0] < tmp64
)
2383 // calculate C* and f*
2384 // C* is actually floor(C*) in this case
2385 // C* and f* need shifting and masking, as shown by
2386 // shiftright128[] and maskhigh128[]
2388 // kx = 10^(-x) = ten2mk128[ind - 1]
2389 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2390 // the approximation of 10^(-x) was rounded up to 118 bits
2391 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
2392 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2393 Cstar
.w
[1] = P256
.w
[3];
2394 Cstar
.w
[0] = P256
.w
[2];
2396 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
2397 fstar
.w
[1] = P256
.w
[1];
2398 fstar
.w
[0] = P256
.w
[0];
2399 } else { // 22 <= ind - 1 <= 33
2401 Cstar
.w
[0] = P256
.w
[3];
2402 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
2403 fstar
.w
[2] = P256
.w
[2];
2404 fstar
.w
[1] = P256
.w
[1];
2405 fstar
.w
[0] = P256
.w
[0];
2407 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2408 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2409 // if (0 < f* < 10^(-x)) then the result is a midpoint
2410 // if floor(C*) is even then C* = floor(C*) - logical right
2411 // shift; C* has p decimal digits, correct by Prop. 1)
2412 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
2413 // shift; C* has p decimal digits, correct by Pr. 1)
2415 // C* = floor(C*) (logical right shift; C has p decimal digits,
2416 // correct by Property 1)
2417 // n = C* * 10^(e+x)
2419 // shift right C* by Ex-128 = shiftright128[ind]
2420 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
2421 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2423 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
2424 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2425 } else { // 22 <= ind - 1 <= 33
2426 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
2428 // determine inexactness of the rounding of C*
2429 // if (0 < f* - 1/2 < 10^(-x)) then
2430 // the result is exact
2431 // else // if (f* - 1/2 > T*) then
2432 // the result is inexact
2434 if (fstar
.w
[1] > 0x8000000000000000ull
||
2435 (fstar
.w
[1] == 0x8000000000000000ull
2436 && fstar
.w
[0] > 0x0ull
)) {
2437 // f* > 1/2 and the result may be exact
2438 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
2439 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
2440 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
2441 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
2442 } // else the result is exact
2443 } else { // the result is inexact; f2* <= 1/2
2444 is_inexact_gt_midpoint
= 1;
2446 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
2447 if (fstar
.w
[3] > 0x0 ||
2448 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
2449 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
2450 (fstar
.w
[1] || fstar
.w
[0]))) {
2451 // f2* > 1/2 and the result may be exact
2452 // Calculate f2* - 1/2
2453 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
2454 tmp64A
= fstar
.w
[3];
2455 if (tmp64
> fstar
.w
[2])
2458 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2459 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2460 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2461 } // else the result is exact
2462 } else { // the result is inexact; f2* <= 1/2
2463 is_inexact_gt_midpoint
= 1;
2465 } else { // if 22 <= ind <= 33
2466 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
2467 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
2468 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
2469 // f2* > 1/2 and the result may be exact
2470 // Calculate f2* - 1/2
2471 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
2472 if (tmp64
|| fstar
.w
[2]
2473 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2474 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2475 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2476 } // else the result is exact
2477 } else { // the result is inexact; f2* <= 1/2
2478 is_inexact_gt_midpoint
= 1;
2482 // if the result was a midpoint it was rounded away from zero, so
2483 // it will need a correction
2484 // check for midpoints
2485 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
2486 && (fstar
.w
[1] || fstar
.w
[0])
2487 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
2488 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2489 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
2490 // the result is a midpoint; round to nearest
2491 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
2492 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
2493 Cstar
.w
[0]--; // Cstar.w[0] is now even
2494 is_inexact_gt_midpoint
= 0;
2495 } else { // else MP in [ODD, EVEN]
2496 is_midpoint_lt_even
= 1;
2497 is_inexact_gt_midpoint
= 0;
2500 // general correction for RZ
2501 if (is_midpoint_lt_even
|| is_inexact_gt_midpoint
) {
2502 Cstar
.w
[0] = Cstar
.w
[0] - 1;
2504 ; // exact, the result is already correct
2507 } else if (exp
== 0) {
2511 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2512 // res = +C * 10^exp (exact)
2513 res
= C1
.w
[0] * ten2k64
[exp
];
2521 /*****************************************************************************
2522 * BID128_to_uint32_xint
2523 ****************************************************************************/
2525 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
2526 bid128_to_uint32_xint
, x
)
2531 int exp
; // unbiased exponent
2532 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2533 UINT64 tmp64
, tmp64A
;
2534 BID_UI64DOUBLE tmp1
;
2535 unsigned int x_nr_bits
;
2538 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
2541 int is_inexact_gt_midpoint
= 0;
2542 int is_midpoint_lt_even
= 0;
2545 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
2546 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
2547 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
2550 // check for NaN or Infinity
2551 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
2553 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
2554 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
2556 *pfpsf
|= INVALID_EXCEPTION
;
2557 // return Integer Indefinite
2559 } else { // x is QNaN
2561 *pfpsf
|= INVALID_EXCEPTION
;
2562 // return Integer Indefinite
2566 } else { // x is not a NaN, so it must be infinity
2567 if (!x_sign
) { // x is +inf
2569 *pfpsf
|= INVALID_EXCEPTION
;
2570 // return Integer Indefinite
2572 } else { // x is -inf
2574 *pfpsf
|= INVALID_EXCEPTION
;
2575 // return Integer Indefinite
2581 // check for non-canonical values (after the check for special values)
2582 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
2583 || (C1
.w
[1] == 0x0001ed09bead87c0ull
2584 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
2585 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
2588 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
2592 } else { // x is not special and is not zero
2594 // q = nr. of decimal digits in x
2595 // determine first the nr. of bits in x
2597 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
2598 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2599 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
2600 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
2602 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2603 } else { // x < 2^32
2604 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
2606 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2608 } else { // if x < 2^53
2609 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
2611 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2613 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2614 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
2616 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2618 q
= nr_digits
[x_nr_bits
- 1].digits
;
2620 q
= nr_digits
[x_nr_bits
- 1].digits1
;
2621 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
2622 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
2623 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
2626 exp
= (x_exp
>> 49) - 6176;
2627 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
2629 *pfpsf
|= INVALID_EXCEPTION
;
2630 // return Integer Indefinite
2633 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
2634 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
2635 // so x rounded to an integer may or may not fit in a signed 32-bit int
2636 // the cases that do not fit are identified here; the ones that fit
2637 // fall through and will be handled with other cases further,
2638 // under '1 <= q + exp <= 10'
2639 if (x_sign
) { // if n < 0 and q + exp = 10
2640 // if n <= -1 then n is too large
2641 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
2642 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a, 1<=q<=34
2644 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
2645 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2646 if (tmp64
>= 0x0aull
) {
2648 *pfpsf
|= INVALID_EXCEPTION
;
2649 // return Integer Indefinite
2653 // else cases that can be rounded to a 32-bit uint fall through
2654 // to '1 <= q + exp <= 10'
2655 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2656 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
2657 // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
2660 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2661 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
2662 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2663 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
2665 if (C1
.w
[1] > C
.w
[1]
2666 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
2668 *pfpsf
|= INVALID_EXCEPTION
;
2669 // return Integer Indefinite
2673 // else cases that can be rounded to a 32-bit int fall through
2674 // to '1 <= q + exp <= 10'
2676 } else { // if n > 0 and q + exp = 10
2677 // if n >= 2^32 then n is too large
2678 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
2679 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
2681 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
2682 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2683 if (tmp64
>= 0xa00000000ull
) {
2685 *pfpsf
|= INVALID_EXCEPTION
;
2686 // return Integer Indefinite
2690 // else cases that can be rounded to a 32-bit uint fall through
2691 // to '1 <= q + exp <= 10'
2692 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2693 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
2694 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
2696 tmp64
= 0xa00000000ull
;
2697 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2698 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
2699 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2700 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
2702 if (C1
.w
[1] > C
.w
[1]
2703 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
2705 *pfpsf
|= INVALID_EXCEPTION
;
2706 // return Integer Indefinite
2710 // else cases that can be rounded to a 32-bit int fall through
2711 // to '1 <= q + exp <= 10'
2715 // n is not too large to be converted to uint32: -2^32 < n < 2^32
2716 // Note: some of the cases tested for above fall through to this point
2717 if ((q
+ exp
) <= 0) {
2718 // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
2720 *pfpsf
|= INEXACT_EXCEPTION
;
2724 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
2725 // x = d(0)...d(k).d(k+1)..., k >= 0, d(0) != 0
2726 if (x_sign
) { // x <= -1
2728 *pfpsf
|= INVALID_EXCEPTION
;
2729 // return Integer Indefinite
2733 // x > 0 from this point on
2734 // 1 <= x < 2^32 so x can be rounded to zero to a 32-bit unsigned integer
2735 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
2736 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
2737 // chop off ind digits from the lower part of C1
2738 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
2741 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
2743 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
2744 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
2746 if (C1
.w
[0] < tmp64
)
2748 // calculate C* and f*
2749 // C* is actually floor(C*) in this case
2750 // C* and f* need shifting and masking, as shown by
2751 // shiftright128[] and maskhigh128[]
2753 // kx = 10^(-x) = ten2mk128[ind - 1]
2754 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2755 // the approximation of 10^(-x) was rounded up to 118 bits
2756 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
2757 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2758 Cstar
.w
[1] = P256
.w
[3];
2759 Cstar
.w
[0] = P256
.w
[2];
2761 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
2762 fstar
.w
[1] = P256
.w
[1];
2763 fstar
.w
[0] = P256
.w
[0];
2764 } else { // 22 <= ind - 1 <= 33
2766 Cstar
.w
[0] = P256
.w
[3];
2767 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
2768 fstar
.w
[2] = P256
.w
[2];
2769 fstar
.w
[1] = P256
.w
[1];
2770 fstar
.w
[0] = P256
.w
[0];
2772 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2773 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2774 // if (0 < f* < 10^(-x)) then the result is a midpoint
2775 // if floor(C*) is even then C* = floor(C*) - logical right
2776 // shift; C* has p decimal digits, correct by Prop. 1)
2777 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
2778 // shift; C* has p decimal digits, correct by Pr. 1)
2780 // C* = floor(C*) (logical right shift; C has p decimal digits,
2781 // correct by Property 1)
2782 // n = C* * 10^(e+x)
2784 // shift right C* by Ex-128 = shiftright128[ind]
2785 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
2786 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2788 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
2789 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2790 } else { // 22 <= ind - 1 <= 33
2791 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
2793 // determine inexactness of the rounding of C*
2794 // if (0 < f* - 1/2 < 10^(-x)) then
2795 // the result is exact
2796 // else // if (f* - 1/2 > T*) then
2797 // the result is inexact
2799 if (fstar
.w
[1] > 0x8000000000000000ull
||
2800 (fstar
.w
[1] == 0x8000000000000000ull
2801 && fstar
.w
[0] > 0x0ull
)) {
2802 // f* > 1/2 and the result may be exact
2803 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
2804 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
2805 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
2806 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
2807 // set the inexact flag
2808 *pfpsf
|= INEXACT_EXCEPTION
;
2809 } // else the result is exact
2810 } else { // the result is inexact; f2* <= 1/2
2811 // set the inexact flag
2812 *pfpsf
|= INEXACT_EXCEPTION
;
2813 is_inexact_gt_midpoint
= 1;
2815 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
2816 if (fstar
.w
[3] > 0x0 ||
2817 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
2818 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
2819 (fstar
.w
[1] || fstar
.w
[0]))) {
2820 // f2* > 1/2 and the result may be exact
2821 // Calculate f2* - 1/2
2822 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
2823 tmp64A
= fstar
.w
[3];
2824 if (tmp64
> fstar
.w
[2])
2827 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2828 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2829 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2830 // set the inexact flag
2831 *pfpsf
|= INEXACT_EXCEPTION
;
2832 } // else the result is exact
2833 } else { // the result is inexact; f2* <= 1/2
2834 // set the inexact flag
2835 *pfpsf
|= INEXACT_EXCEPTION
;
2836 is_inexact_gt_midpoint
= 1;
2838 } else { // if 22 <= ind <= 33
2839 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
2840 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
2841 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
2842 // f2* > 1/2 and the result may be exact
2843 // Calculate f2* - 1/2
2844 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
2845 if (tmp64
|| fstar
.w
[2]
2846 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2847 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2848 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2849 // set the inexact flag
2850 *pfpsf
|= INEXACT_EXCEPTION
;
2851 } // else the result is exact
2852 } else { // the result is inexact; f2* <= 1/2
2853 // set the inexact flag
2854 *pfpsf
|= INEXACT_EXCEPTION
;
2855 is_inexact_gt_midpoint
= 1;
2859 // if the result was a midpoint it was rounded away from zero, so
2860 // it will need a correction
2861 // check for midpoints
2862 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
2863 && (fstar
.w
[1] || fstar
.w
[0])
2864 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
2865 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2866 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
2867 // the result is a midpoint; round to nearest
2868 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
2869 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
2870 Cstar
.w
[0]--; // Cstar.w[0] is now even
2871 is_inexact_gt_midpoint
= 0;
2872 } else { // else MP in [ODD, EVEN]
2873 is_midpoint_lt_even
= 1;
2874 is_inexact_gt_midpoint
= 0;
2877 // general correction for RZ
2878 if (is_midpoint_lt_even
|| is_inexact_gt_midpoint
) {
2879 Cstar
.w
[0] = Cstar
.w
[0] - 1;
2881 ; // exact, the result is already correct
2884 } else if (exp
== 0) {
2888 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2889 // res = +C * 10^exp (exact)
2890 res
= C1
.w
[0] * ten2k64
[exp
];
2898 /*****************************************************************************
2899 * BID128_to_uint32_rninta
2900 ****************************************************************************/
2902 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
2903 bid128_to_uint32_rninta
, x
)
2908 int exp
; // unbiased exponent
2909 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2911 BID_UI64DOUBLE tmp1
;
2912 unsigned int x_nr_bits
;
2915 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
2919 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
2920 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
2921 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
2924 // check for NaN or Infinity
2925 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
2927 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
2928 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
2930 *pfpsf
|= INVALID_EXCEPTION
;
2931 // return Integer Indefinite
2933 } else { // x is QNaN
2935 *pfpsf
|= INVALID_EXCEPTION
;
2936 // return Integer Indefinite
2940 } else { // x is not a NaN, so it must be infinity
2941 if (!x_sign
) { // x is +inf
2943 *pfpsf
|= INVALID_EXCEPTION
;
2944 // return Integer Indefinite
2946 } else { // x is -inf
2948 *pfpsf
|= INVALID_EXCEPTION
;
2949 // return Integer Indefinite
2955 // check for non-canonical values (after the check for special values)
2956 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
2957 || (C1
.w
[1] == 0x0001ed09bead87c0ull
2958 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
2959 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
2962 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
2966 } else { // x is not special and is not zero
2968 // q = nr. of decimal digits in x
2969 // determine first the nr. of bits in x
2971 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
2972 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2973 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
2974 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
2976 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2977 } else { // x < 2^32
2978 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
2980 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2982 } else { // if x < 2^53
2983 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
2985 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2987 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2988 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
2990 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2992 q
= nr_digits
[x_nr_bits
- 1].digits
;
2994 q
= nr_digits
[x_nr_bits
- 1].digits1
;
2995 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
2996 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
2997 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
3000 exp
= (x_exp
>> 49) - 6176;
3001 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
3003 *pfpsf
|= INVALID_EXCEPTION
;
3004 // return Integer Indefinite
3007 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
3008 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
3009 // so x rounded to an integer may or may not fit in a signed 32-bit int
3010 // the cases that do not fit are identified here; the ones that fit
3011 // fall through and will be handled with other cases further,
3012 // under '1 <= q + exp <= 10'
3013 if (x_sign
) { // if n < 0 and q + exp = 10
3014 // if n <= -1/2 then n is too large
3015 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1/2
3016 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05, 1<=q<=34
3018 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
3019 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3020 if (tmp64
>= 0x05ull
) {
3022 *pfpsf
|= INVALID_EXCEPTION
;
3023 // return Integer Indefinite
3027 // else cases that can be rounded to a 32-bit int fall through
3028 // to '1 <= q + exp <= 10'
3029 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3030 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05 <=>
3031 // C >= 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
3034 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3035 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
3036 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3037 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
3039 if (C1
.w
[1] > C
.w
[1]
3040 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
3042 *pfpsf
|= INVALID_EXCEPTION
;
3043 // return Integer Indefinite
3047 // else cases that can be rounded to a 32-bit int fall through
3048 // to '1 <= q + exp <= 10'
3050 } else { // if n > 0 and q + exp = 10
3051 // if n >= 2^32 - 1/2 then n is too large
3052 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
3053 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
3055 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
3056 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3057 if (tmp64
>= 0x9fffffffbull
) {
3059 *pfpsf
|= INVALID_EXCEPTION
;
3060 // return Integer Indefinite
3064 // else cases that can be rounded to a 32-bit int fall through
3065 // to '1 <= q + exp <= 10'
3066 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3067 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
3068 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
3069 // (scale 2^32-1/2 up)
3070 tmp64
= 0x9fffffffbull
;
3071 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3072 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
3073 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3074 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
3076 if (C1
.w
[1] > C
.w
[1]
3077 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
3079 *pfpsf
|= INVALID_EXCEPTION
;
3080 // return Integer Indefinite
3084 // else cases that can be rounded to a 32-bit int fall through
3085 // to '1 <= q + exp <= 10'
3089 // n is not too large to be converted to int32: -1/2 < n < 2^32 - 1/2
3090 // Note: some of the cases tested for above fall through to this point
3091 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
3095 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
3096 // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
3103 if (ind
<= 18) { // 0 <= ind <= 18
3104 if ((C1
.w
[1] == 0) && (C1
.w
[0] < midpoint64
[ind
])) {
3105 res
= 0x00000000; // return 0
3106 } else if (!x_sign
) { // n > 0
3107 res
= 0x00000001; // return +1
3110 *pfpsf
|= INVALID_EXCEPTION
;
3112 } else { // 19 <= ind <= 33
3113 if ((C1
.w
[1] < midpoint128
[ind
- 19].w
[1])
3114 || ((C1
.w
[1] == midpoint128
[ind
- 19].w
[1])
3115 && (C1
.w
[0] < midpoint128
[ind
- 19].w
[0]))) {
3116 res
= 0x00000000; // return 0
3117 } else if (!x_sign
) { // n > 0
3118 res
= 0x00000001; // return +1
3121 *pfpsf
|= INVALID_EXCEPTION
;
3124 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
3125 if (x_sign
) { // x <= -1
3127 *pfpsf
|= INVALID_EXCEPTION
;
3128 // return Integer Indefinite
3132 // 1 <= x < 2^31-1/2 so x can be rounded
3133 // to nearest-away to a 32-bit signed integer
3134 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
3135 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
3136 // chop off ind digits from the lower part of C1
3137 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
3140 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
3142 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
3143 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
3145 if (C1
.w
[0] < tmp64
)
3147 // calculate C* and f*
3148 // C* is actually floor(C*) in this case
3149 // C* and f* need shifting and masking, as shown by
3150 // shiftright128[] and maskhigh128[]
3152 // kx = 10^(-x) = ten2mk128[ind - 1]
3153 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
3154 // the approximation of 10^(-x) was rounded up to 118 bits
3155 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
3156 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
3157 Cstar
.w
[1] = P256
.w
[3];
3158 Cstar
.w
[0] = P256
.w
[2];
3159 } else { // 22 <= ind - 1 <= 33
3161 Cstar
.w
[0] = P256
.w
[3];
3163 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
3164 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
3165 // if (0 < f* < 10^(-x)) then the result is a midpoint
3166 // if floor(C*) is even then C* = floor(C*) - logical right
3167 // shift; C* has p decimal digits, correct by Prop. 1)
3168 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
3169 // shift; C* has p decimal digits, correct by Pr. 1)
3171 // C* = floor(C*) (logical right shift; C has p decimal digits,
3172 // correct by Property 1)
3173 // n = C* * 10^(e+x)
3175 // shift right C* by Ex-128 = shiftright128[ind]
3176 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
3177 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
3179 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
3180 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
3181 } else { // 22 <= ind - 1 <= 33
3182 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
3184 // if the result was a midpoint, it was already rounded away from zero
3185 res
= Cstar
.w
[0]; // always positive
3186 // no need to check for midpoints - already rounded away from zero!
3187 } else if (exp
== 0) {
3191 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
3192 // res = +C * 10^exp (exact)
3193 res
= C1
.w
[0] * ten2k64
[exp
];
3201 /*****************************************************************************
3202 * BID128_to_uint32_xrninta
3203 ****************************************************************************/
3205 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
3206 bid128_to_uint32_xrninta
, x
)
3211 int exp
; // unbiased exponent
3212 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
3213 UINT64 tmp64
, tmp64A
;
3214 BID_UI64DOUBLE tmp1
;
3215 unsigned int x_nr_bits
;
3218 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
3221 unsigned int tmp_inexact
= 0;
3224 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
3225 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
3226 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
3229 // check for NaN or Infinity
3230 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
3232 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
3233 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
3235 *pfpsf
|= INVALID_EXCEPTION
;
3236 // return Integer Indefinite
3238 } else { // x is QNaN
3240 *pfpsf
|= INVALID_EXCEPTION
;
3241 // return Integer Indefinite
3245 } else { // x is not a NaN, so it must be infinity
3246 if (!x_sign
) { // x is +inf
3248 *pfpsf
|= INVALID_EXCEPTION
;
3249 // return Integer Indefinite
3251 } else { // x is -inf
3253 *pfpsf
|= INVALID_EXCEPTION
;
3254 // return Integer Indefinite
3260 // check for non-canonical values (after the check for special values)
3261 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
3262 || (C1
.w
[1] == 0x0001ed09bead87c0ull
3263 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
3264 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
3267 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
3271 } else { // x is not special and is not zero
3273 // q = nr. of decimal digits in x
3274 // determine first the nr. of bits in x
3276 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
3277 // split the 64-bit value in two 32-bit halves to avoid rounding errors
3278 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
3279 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
3281 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3282 } else { // x < 2^32
3283 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
3285 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3287 } else { // if x < 2^53
3288 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
3290 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3292 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
3293 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
3295 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3297 q
= nr_digits
[x_nr_bits
- 1].digits
;
3299 q
= nr_digits
[x_nr_bits
- 1].digits1
;
3300 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
3301 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
3302 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
3305 exp
= (x_exp
>> 49) - 6176;
3306 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
3308 *pfpsf
|= INVALID_EXCEPTION
;
3309 // return Integer Indefinite
3312 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
3313 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
3314 // so x rounded to an integer may or may not fit in a signed 32-bit int
3315 // the cases that do not fit are identified here; the ones that fit
3316 // fall through and will be handled with other cases further,
3317 // under '1 <= q + exp <= 10'
3318 if (x_sign
) { // if n < 0 and q + exp = 10
3319 // if n <= -1/2 then n is too large
3320 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1/2
3321 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05, 1<=q<=34
3323 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
3324 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3325 if (tmp64
>= 0x05ull
) {
3327 *pfpsf
|= INVALID_EXCEPTION
;
3328 // return Integer Indefinite
3332 // else cases that can be rounded to a 32-bit int fall through
3333 // to '1 <= q + exp <= 10'
3334 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3335 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05 <=>
3336 // C >= 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
3339 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3340 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
3341 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3342 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
3344 if (C1
.w
[1] > C
.w
[1]
3345 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
3347 *pfpsf
|= INVALID_EXCEPTION
;
3348 // return Integer Indefinite
3352 // else cases that can be rounded to a 32-bit int fall through
3353 // to '1 <= q + exp <= 10'
3355 } else { // if n > 0 and q + exp = 10
3356 // if n >= 2^32 - 1/2 then n is too large
3357 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
3358 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
3360 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
3361 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3362 if (tmp64
>= 0x9fffffffbull
) {
3364 *pfpsf
|= INVALID_EXCEPTION
;
3365 // return Integer Indefinite
3369 // else cases that can be rounded to a 32-bit int fall through
3370 // to '1 <= q + exp <= 10'
3371 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3372 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
3373 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
3374 // (scale 2^32-1/2 up)
3375 tmp64
= 0x9fffffffbull
;
3376 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3377 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
3378 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3379 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
3381 if (C1
.w
[1] > C
.w
[1]
3382 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
3384 *pfpsf
|= INVALID_EXCEPTION
;
3385 // return Integer Indefinite
3389 // else cases that can be rounded to a 32-bit int fall through
3390 // to '1 <= q + exp <= 10'
3394 // n is not too large to be converted to int32: -1/2 < n < 2^32 - 1/2
3395 // Note: some of the cases tested for above fall through to this point
3396 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
3398 *pfpsf
|= INEXACT_EXCEPTION
;
3402 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
3403 // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
3410 if (ind
<= 18) { // 0 <= ind <= 18
3411 if ((C1
.w
[1] == 0) && (C1
.w
[0] < midpoint64
[ind
])) {
3412 res
= 0x00000000; // return 0
3413 } else if (!x_sign
) { // n > 0
3414 res
= 0x00000001; // return +1
3417 *pfpsf
|= INVALID_EXCEPTION
;
3420 } else { // 19 <= ind <= 33
3421 if ((C1
.w
[1] < midpoint128
[ind
- 19].w
[1])
3422 || ((C1
.w
[1] == midpoint128
[ind
- 19].w
[1])
3423 && (C1
.w
[0] < midpoint128
[ind
- 19].w
[0]))) {
3424 res
= 0x00000000; // return 0
3425 } else if (!x_sign
) { // n > 0
3426 res
= 0x00000001; // return +1
3429 *pfpsf
|= INVALID_EXCEPTION
;
3434 *pfpsf
|= INEXACT_EXCEPTION
;
3435 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
3436 if (x_sign
) { // x <= -1
3438 *pfpsf
|= INVALID_EXCEPTION
;
3439 // return Integer Indefinite
3443 // 1 <= x < 2^31-1/2 so x can be rounded
3444 // to nearest-away to a 32-bit signed integer
3445 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
3446 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
3447 // chop off ind digits from the lower part of C1
3448 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
3451 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
3453 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
3454 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
3456 if (C1
.w
[0] < tmp64
)
3458 // calculate C* and f*
3459 // C* is actually floor(C*) in this case
3460 // C* and f* need shifting and masking, as shown by
3461 // shiftright128[] and maskhigh128[]
3463 // kx = 10^(-x) = ten2mk128[ind - 1]
3464 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
3465 // the approximation of 10^(-x) was rounded up to 118 bits
3466 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
3467 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
3468 Cstar
.w
[1] = P256
.w
[3];
3469 Cstar
.w
[0] = P256
.w
[2];
3471 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
3472 fstar
.w
[1] = P256
.w
[1];
3473 fstar
.w
[0] = P256
.w
[0];
3474 } else { // 22 <= ind - 1 <= 33
3476 Cstar
.w
[0] = P256
.w
[3];
3477 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
3478 fstar
.w
[2] = P256
.w
[2];
3479 fstar
.w
[1] = P256
.w
[1];
3480 fstar
.w
[0] = P256
.w
[0];
3482 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
3483 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
3484 // if (0 < f* < 10^(-x)) then the result is a midpoint
3485 // if floor(C*) is even then C* = floor(C*) - logical right
3486 // shift; C* has p decimal digits, correct by Prop. 1)
3487 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
3488 // shift; C* has p decimal digits, correct by Pr. 1)
3490 // C* = floor(C*) (logical right shift; C has p decimal digits,
3491 // correct by Property 1)
3492 // n = C* * 10^(e+x)
3494 // shift right C* by Ex-128 = shiftright128[ind]
3495 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
3496 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
3498 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
3499 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
3500 } else { // 22 <= ind - 1 <= 33
3501 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
3503 // if the result was a midpoint, it was already rounded away from zero
3504 // determine inexactness of the rounding of C*
3505 // if (0 < f* - 1/2 < 10^(-x)) then
3506 // the result is exact
3507 // else // if (f* - 1/2 > T*) then
3508 // the result is inexact
3510 if (fstar
.w
[1] > 0x8000000000000000ull
||
3511 (fstar
.w
[1] == 0x8000000000000000ull
3512 && fstar
.w
[0] > 0x0ull
)) {
3513 // f* > 1/2 and the result may be exact
3514 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
3515 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
3516 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
3517 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
3518 // set the inexact flag
3519 // *pfpsf |= INEXACT_EXCEPTION;
3521 } // else the result is exact
3522 } else { // the result is inexact; f2* <= 1/2
3523 // set the inexact flag
3524 // *pfpsf |= INEXACT_EXCEPTION;
3527 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
3528 if (fstar
.w
[3] > 0x0 ||
3529 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
3530 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
3531 (fstar
.w
[1] || fstar
.w
[0]))) {
3532 // f2* > 1/2 and the result may be exact
3533 // Calculate f2* - 1/2
3534 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
3535 tmp64A
= fstar
.w
[3];
3536 if (tmp64
> fstar
.w
[2])
3539 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
3540 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
3541 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
3542 // set the inexact flag
3543 // *pfpsf |= INEXACT_EXCEPTION;
3545 } // else the result is exact
3546 } else { // the result is inexact; f2* <= 1/2
3547 // set the inexact flag
3548 // *pfpsf |= INEXACT_EXCEPTION;
3551 } else { // if 22 <= ind <= 33
3552 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
3553 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
3554 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
3555 // f2* > 1/2 and the result may be exact
3556 // Calculate f2* - 1/2
3557 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
3558 if (tmp64
|| fstar
.w
[2]
3559 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
3560 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
3561 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
3562 // set the inexact flag
3563 // *pfpsf |= INEXACT_EXCEPTION;
3565 } // else the result is exact
3566 } else { // the result is inexact; f2* <= 1/2
3567 // set the inexact flag
3568 // *pfpsf |= INEXACT_EXCEPTION;
3572 // no need to check for midpoints - already rounded away from zero!
3573 res
= Cstar
.w
[0]; // the result is positive
3575 *pfpsf
|= INEXACT_EXCEPTION
;
3576 } else if (exp
== 0) {
3580 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
3581 // res = +C * 10^exp (exact)
3582 res
= C1
.w
[0] * ten2k64
[exp
];