1 /* Copyright (C) 2007-2018 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_int32_rnint
28 ****************************************************************************/
30 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_rnint
, x
)
35 int exp
; // unbiased exponent
36 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
39 unsigned int x_nr_bits
;
42 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
47 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
48 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
49 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
52 // check for NaN or Infinity
53 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
55 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
56 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
58 *pfpsf
|= INVALID_EXCEPTION
;
59 // return Integer Indefinite
63 *pfpsf
|= INVALID_EXCEPTION
;
64 // return Integer Indefinite
68 } else { // x is not a NaN, so it must be infinity
69 if (!x_sign
) { // x is +inf
71 *pfpsf
|= INVALID_EXCEPTION
;
72 // return Integer Indefinite
76 *pfpsf
|= INVALID_EXCEPTION
;
77 // return Integer Indefinite
83 // check for non-canonical values (after the check for special values)
84 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
85 || (C1
.w
[1] == 0x0001ed09bead87c0ull
86 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
87 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
90 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
94 } else { // x is not special and is not zero
96 // q = nr. of decimal digits in x
97 // determine first the nr. of bits in x
99 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
100 // split the 64-bit value in two 32-bit halves to avoid rounding errors
101 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
102 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
104 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
106 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
108 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
110 } else { // if x < 2^53
111 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
113 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
115 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
116 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
118 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
120 q
= nr_digits
[x_nr_bits
- 1].digits
;
122 q
= nr_digits
[x_nr_bits
- 1].digits1
;
123 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
124 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
125 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
128 exp
= (x_exp
>> 49) - 6176;
129 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
131 *pfpsf
|= INVALID_EXCEPTION
;
132 // return Integer Indefinite
135 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
136 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
137 // so x rounded to an integer may or may not fit in a signed 32-bit int
138 // the cases that do not fit are identified here; the ones that fit
139 // fall through and will be handled with other cases further,
140 // under '1 <= q + exp <= 10'
141 if (x_sign
) { // if n < 0 and q + exp = 10
142 // if n < -2^31 - 1/2 then n is too large
143 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31+1/2
144 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000005, 1<=q<=34
146 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
147 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
148 if (tmp64
> 0x500000005ull
) {
150 *pfpsf
|= INVALID_EXCEPTION
;
151 // return Integer Indefinite
155 // else cases that can be rounded to a 32-bit int fall through
156 // to '1 <= q + exp <= 10'
157 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
158 // 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000005 <=>
159 // C > 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 23
160 // (scale 2^31+1/2 up)
161 tmp64
= 0x500000005ull
;
162 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
163 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
164 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
165 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
167 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] > C
.w
[0])) {
169 *pfpsf
|= INVALID_EXCEPTION
;
170 // return Integer Indefinite
174 // else cases that can be rounded to a 32-bit int fall through
175 // to '1 <= q + exp <= 10'
177 } else { // if n > 0 and q + exp = 10
178 // if n >= 2^31 - 1/2 then n is too large
179 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2
180 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=34
182 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
183 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
184 if (tmp64
>= 0x4fffffffbull
) {
186 *pfpsf
|= INVALID_EXCEPTION
;
187 // return Integer Indefinite
191 // else cases that can be rounded to a 32-bit int fall through
192 // to '1 <= q + exp <= 10'
193 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
194 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb <=>
195 // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
196 // (scale 2^31-1/2 up)
197 tmp64
= 0x4fffffffbull
;
198 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
199 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
200 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
201 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
204 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
206 *pfpsf
|= INVALID_EXCEPTION
;
207 // return Integer Indefinite
211 // else cases that can be rounded to a 32-bit int fall through
212 // to '1 <= q + exp <= 10'
216 // n is not too large to be converted to int32: -2^31 - 1/2 < n < 2^31 - 1/2
217 // Note: some of the cases tested for above fall through to this point
218 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
222 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
223 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
228 if (ind
<= 18) { // 0 <= ind <= 18
229 if ((C1
.w
[1] == 0) && (C1
.w
[0] <= midpoint64
[ind
])) {
230 res
= 0x00000000; // return 0
231 } else if (x_sign
) { // n < 0
232 res
= 0xffffffff; // return -1
234 res
= 0x00000001; // return +1
236 } else { // 19 <= ind <= 33
237 if ((C1
.w
[1] < midpoint128
[ind
- 19].w
[1])
238 || ((C1
.w
[1] == midpoint128
[ind
- 19].w
[1])
239 && (C1
.w
[0] <= midpoint128
[ind
- 19].w
[0]))) {
240 res
= 0x00000000; // return 0
241 } else if (x_sign
) { // n < 0
242 res
= 0xffffffff; // return -1
244 res
= 0x00000001; // return +1
247 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
248 // -2^31-1/2 <= x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded
249 // to nearest to a 32-bit signed integer
250 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
251 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
252 // chop off ind digits from the lower part of C1
253 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
256 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
258 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
259 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
263 // calculate C* and f*
264 // C* is actually floor(C*) in this case
265 // C* and f* need shifting and masking, as shown by
266 // shiftright128[] and maskhigh128[]
268 // kx = 10^(-x) = ten2mk128[ind - 1]
269 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
270 // the approximation of 10^(-x) was rounded up to 118 bits
271 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
272 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
273 Cstar
.w
[1] = P256
.w
[3];
274 Cstar
.w
[0] = P256
.w
[2];
276 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
277 fstar
.w
[1] = P256
.w
[1];
278 fstar
.w
[0] = P256
.w
[0];
279 } else { // 22 <= ind - 1 <= 33
281 Cstar
.w
[0] = P256
.w
[3];
282 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
283 fstar
.w
[2] = P256
.w
[2];
284 fstar
.w
[1] = P256
.w
[1];
285 fstar
.w
[0] = P256
.w
[0];
287 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
288 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
289 // if (0 < f* < 10^(-x)) then the result is a midpoint
290 // if floor(C*) is even then C* = floor(C*) - logical right
291 // shift; C* has p decimal digits, correct by Prop. 1)
292 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
293 // shift; C* has p decimal digits, correct by Pr. 1)
295 // C* = floor(C*) (logical right shift; C has p decimal digits,
296 // correct by Property 1)
299 // shift right C* by Ex-128 = shiftright128[ind]
300 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
301 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
303 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
304 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
305 } else { // 22 <= ind - 1 <= 33
306 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
308 // if the result was a midpoint it was rounded away from zero, so
309 // it will need a correction
310 // check for midpoints
311 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
312 && (fstar
.w
[1] || fstar
.w
[0])
313 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
314 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
315 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
316 // the result is a midpoint; round to nearest
317 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
318 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
319 Cstar
.w
[0]--; // Cstar.w[0] is now even
320 } // else MP in [ODD, EVEN]
326 } else if (exp
== 0) {
328 // res = +/-C (exact)
333 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
334 // res = +/-C * 10^exp (exact)
336 res
= -C1
.w
[0] * ten2k64
[exp
];
338 res
= C1
.w
[0] * ten2k64
[exp
];
346 /*****************************************************************************
347 * BID128_to_int32_xrnint
348 ****************************************************************************/
350 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_xrnint
,
356 int exp
; // unbiased exponent
357 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
358 UINT64 tmp64
, tmp64A
;
360 unsigned int x_nr_bits
;
363 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
368 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
369 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
370 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
373 // check for NaN or Infinity
374 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
376 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
377 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
379 *pfpsf
|= INVALID_EXCEPTION
;
380 // return Integer Indefinite
382 } else { // x is QNaN
384 *pfpsf
|= INVALID_EXCEPTION
;
385 // return Integer Indefinite
389 } else { // x is not a NaN, so it must be infinity
390 if (!x_sign
) { // x is +inf
392 *pfpsf
|= INVALID_EXCEPTION
;
393 // return Integer Indefinite
395 } else { // x is -inf
397 *pfpsf
|= INVALID_EXCEPTION
;
398 // return Integer Indefinite
404 // check for non-canonical values (after the check for special values)
405 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
406 || (C1
.w
[1] == 0x0001ed09bead87c0ull
407 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
408 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
411 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
415 } else { // x is not special and is not zero
417 // q = nr. of decimal digits in x
418 // determine first the nr. of bits in x
420 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
421 // split the 64-bit value in two 32-bit halves to avoid rounding errors
422 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
423 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
425 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
427 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
429 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
431 } else { // if x < 2^53
432 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
434 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
436 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
437 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
439 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
441 q
= nr_digits
[x_nr_bits
- 1].digits
;
443 q
= nr_digits
[x_nr_bits
- 1].digits1
;
444 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
445 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
446 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
449 exp
= (x_exp
>> 49) - 6176;
450 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
452 *pfpsf
|= INVALID_EXCEPTION
;
453 // return Integer Indefinite
456 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
457 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
458 // so x rounded to an integer may or may not fit in a signed 32-bit int
459 // the cases that do not fit are identified here; the ones that fit
460 // fall through and will be handled with other cases further,
461 // under '1 <= q + exp <= 10'
462 if (x_sign
) { // if n < 0 and q + exp = 10
463 // if n < -2^31 - 1/2 then n is too large
464 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31+1/2
465 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000005, 1<=q<=34
467 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
468 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
469 if (tmp64
> 0x500000005ull
) {
471 *pfpsf
|= INVALID_EXCEPTION
;
472 // return Integer Indefinite
476 // else cases that can be rounded to a 32-bit int fall through
477 // to '1 <= q + exp <= 10'
478 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
479 // 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000005 <=>
480 // C > 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 23
481 // (scale 2^31+1/2 up)
482 tmp64
= 0x500000005ull
;
483 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
484 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
485 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
486 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
488 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] > C
.w
[0])) {
490 *pfpsf
|= INVALID_EXCEPTION
;
491 // return Integer Indefinite
495 // else cases that can be rounded to a 32-bit int fall through
496 // to '1 <= q + exp <= 10'
498 } else { // if n > 0 and q + exp = 10
499 // if n >= 2^31 - 1/2 then n is too large
500 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2
501 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=34
503 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
504 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
505 if (tmp64
>= 0x4fffffffbull
) {
507 *pfpsf
|= INVALID_EXCEPTION
;
508 // return Integer Indefinite
512 // else cases that can be rounded to a 32-bit int fall through
513 // to '1 <= q + exp <= 10'
514 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
515 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb <=>
516 // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
517 // (scale 2^31-1/2 up)
518 tmp64
= 0x4fffffffbull
;
519 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
520 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
521 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
522 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
525 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
527 *pfpsf
|= INVALID_EXCEPTION
;
528 // return Integer Indefinite
532 // else cases that can be rounded to a 32-bit int fall through
533 // to '1 <= q + exp <= 10'
537 // n is not too large to be converted to int32: -2^31 - 1/2 < n < 2^31 - 1/2
538 // Note: some of the cases tested for above fall through to this point
539 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
541 *pfpsf
|= INEXACT_EXCEPTION
;
545 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
546 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
551 if (ind
<= 18) { // 0 <= ind <= 18
552 if ((C1
.w
[1] == 0) && (C1
.w
[0] <= midpoint64
[ind
])) {
553 res
= 0x00000000; // return 0
554 } else if (x_sign
) { // n < 0
555 res
= 0xffffffff; // return -1
557 res
= 0x00000001; // return +1
559 } else { // 19 <= ind <= 33
560 if ((C1
.w
[1] < midpoint128
[ind
- 19].w
[1])
561 || ((C1
.w
[1] == midpoint128
[ind
- 19].w
[1])
562 && (C1
.w
[0] <= midpoint128
[ind
- 19].w
[0]))) {
563 res
= 0x00000000; // return 0
564 } else if (x_sign
) { // n < 0
565 res
= 0xffffffff; // return -1
567 res
= 0x00000001; // return +1
571 *pfpsf
|= INEXACT_EXCEPTION
;
572 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
573 // -2^31-1/2 <= x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded
574 // to nearest to a 32-bit signed integer
575 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
576 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
577 // chop off ind digits from the lower part of C1
578 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
581 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
583 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
584 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
588 // calculate C* and f*
589 // C* is actually floor(C*) in this case
590 // C* and f* need shifting and masking, as shown by
591 // shiftright128[] and maskhigh128[]
593 // kx = 10^(-x) = ten2mk128[ind - 1]
594 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
595 // the approximation of 10^(-x) was rounded up to 118 bits
596 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
597 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
598 Cstar
.w
[1] = P256
.w
[3];
599 Cstar
.w
[0] = P256
.w
[2];
601 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
602 fstar
.w
[1] = P256
.w
[1];
603 fstar
.w
[0] = P256
.w
[0];
604 } else { // 22 <= ind - 1 <= 33
606 Cstar
.w
[0] = P256
.w
[3];
607 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
608 fstar
.w
[2] = P256
.w
[2];
609 fstar
.w
[1] = P256
.w
[1];
610 fstar
.w
[0] = P256
.w
[0];
612 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
613 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
614 // if (0 < f* < 10^(-x)) then the result is a midpoint
615 // if floor(C*) is even then C* = floor(C*) - logical right
616 // shift; C* has p decimal digits, correct by Prop. 1)
617 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
618 // shift; C* has p decimal digits, correct by Pr. 1)
620 // C* = floor(C*) (logical right shift; C has p decimal digits,
621 // correct by Property 1)
624 // shift right C* by Ex-128 = shiftright128[ind]
625 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
626 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
628 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
629 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
630 } else { // 22 <= ind - 1 <= 33
631 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
633 // determine inexactness of the rounding of C*
634 // if (0 < f* - 1/2 < 10^(-x)) then
635 // the result is exact
636 // else // if (f* - 1/2 > T*) then
637 // the result is inexact
639 if (fstar
.w
[1] > 0x8000000000000000ull
|| (fstar
.w
[1] == 0x8000000000000000ull
&& fstar
.w
[0] > 0x0ull
)) { // f* > 1/2 and the result may be exact
640 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
641 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
642 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
643 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
644 // set the inexact flag
645 *pfpsf
|= INEXACT_EXCEPTION
;
646 } // else the result is exact
647 } else { // the result is inexact; f2* <= 1/2
648 // set the inexact flag
649 *pfpsf
|= INEXACT_EXCEPTION
;
651 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
652 if (fstar
.w
[3] > 0x0 ||
653 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
654 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
655 (fstar
.w
[1] || fstar
.w
[0]))) {
656 // f2* > 1/2 and the result may be exact
657 // Calculate f2* - 1/2
658 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
660 if (tmp64
> fstar
.w
[2])
663 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
664 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
665 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
666 // set the inexact flag
667 *pfpsf
|= INEXACT_EXCEPTION
;
668 } // else the result is exact
669 } else { // the result is inexact; f2* <= 1/2
670 // set the inexact flag
671 *pfpsf
|= INEXACT_EXCEPTION
;
673 } else { // if 22 <= ind <= 33
674 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
675 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
676 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
677 // f2* > 1/2 and the result may be exact
678 // Calculate f2* - 1/2
679 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
680 if (tmp64
|| fstar
.w
[2]
681 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
682 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
683 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
684 // set the inexact flag
685 *pfpsf
|= INEXACT_EXCEPTION
;
686 } // else the result is exact
687 } else { // the result is inexact; f2* <= 1/2
688 // set the inexact flag
689 *pfpsf
|= INEXACT_EXCEPTION
;
692 // if the result was a midpoint it was rounded away from zero, so
693 // it will need a correction
694 // check for midpoints
695 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
696 && (fstar
.w
[1] || fstar
.w
[0])
697 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
698 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
699 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
700 // the result is a midpoint; round to nearest
701 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
702 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
703 Cstar
.w
[0]--; // Cstar.w[0] is now even
704 } // else MP in [ODD, EVEN]
710 } else if (exp
== 0) {
712 // res = +/-C (exact)
717 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
718 // res = +/-C * 10^exp (exact)
720 res
= -C1
.w
[0] * ten2k64
[exp
];
722 res
= C1
.w
[0] * ten2k64
[exp
];
730 /*****************************************************************************
731 * BID128_to_int32_floor
732 ****************************************************************************/
734 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_floor
, x
)
739 int exp
; // unbiased exponent
740 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
741 UINT64 tmp64
, tmp64A
;
743 unsigned int x_nr_bits
;
746 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
749 int is_inexact_lt_midpoint
= 0;
750 int is_inexact_gt_midpoint
= 0;
751 int is_midpoint_lt_even
= 0;
752 int is_midpoint_gt_even
= 0;
755 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
756 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
757 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
760 // check for NaN or Infinity
761 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
763 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
764 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
766 *pfpsf
|= INVALID_EXCEPTION
;
767 // return Integer Indefinite
769 } else { // x is QNaN
771 *pfpsf
|= INVALID_EXCEPTION
;
772 // return Integer Indefinite
776 } else { // x is not a NaN, so it must be infinity
777 if (!x_sign
) { // x is +inf
779 *pfpsf
|= INVALID_EXCEPTION
;
780 // return Integer Indefinite
782 } else { // x is -inf
784 *pfpsf
|= INVALID_EXCEPTION
;
785 // return Integer Indefinite
791 // check for non-canonical values (after the check for special values)
792 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
793 || (C1
.w
[1] == 0x0001ed09bead87c0ull
794 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
795 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
798 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
802 } else { // x is not special and is not zero
804 // q = nr. of decimal digits in x
805 // determine first the nr. of bits in x
807 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
808 // split the 64-bit value in two 32-bit halves to avoid rounding errors
809 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
810 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
812 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
814 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
816 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
818 } else { // if x < 2^53
819 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
821 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
823 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
824 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
826 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
828 q
= nr_digits
[x_nr_bits
- 1].digits
;
830 q
= nr_digits
[x_nr_bits
- 1].digits1
;
831 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
832 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
833 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
836 exp
= (x_exp
>> 49) - 6176;
837 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
839 *pfpsf
|= INVALID_EXCEPTION
;
840 // return Integer Indefinite
843 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
844 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
845 // so x rounded to an integer may or may not fit in a signed 32-bit int
846 // the cases that do not fit are identified here; the ones that fit
847 // fall through and will be handled with other cases further,
848 // under '1 <= q + exp <= 10'
849 if (x_sign
) { // if n < 0 and q + exp = 10
850 // if n < -2^31 then n is too large
851 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31
852 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000000, 1<=q<=34
854 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
855 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
856 if (tmp64
> 0x500000000ull
) {
858 *pfpsf
|= INVALID_EXCEPTION
;
859 // return Integer Indefinite
863 // else cases that can be rounded to a 32-bit int fall through
864 // to '1 <= q + exp <= 10'
865 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
866 // 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000000 <=>
867 // C > 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 23
869 tmp64
= 0x500000000ull
;
870 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
871 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
872 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
873 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
875 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] > C
.w
[0])) {
877 *pfpsf
|= INVALID_EXCEPTION
;
878 // return Integer Indefinite
882 // else cases that can be rounded to a 32-bit int fall through
883 // to '1 <= q + exp <= 10'
885 } else { // if n > 0 and q + exp = 10
886 // if n >= 2^31 then n is too large
887 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31
888 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=34
890 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
891 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
892 if (tmp64
>= 0x500000000ull
) {
894 *pfpsf
|= INVALID_EXCEPTION
;
895 // return Integer Indefinite
899 // else cases that can be rounded to a 32-bit int fall through
900 // to '1 <= q + exp <= 10'
901 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
902 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000 <=>
903 // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 23
905 tmp64
= 0x500000000ull
;
906 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
907 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
908 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
909 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
912 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
914 *pfpsf
|= INVALID_EXCEPTION
;
915 // return Integer Indefinite
919 // else cases that can be rounded to a 32-bit int fall through
920 // to '1 <= q + exp <= 10'
924 // n is not too large to be converted to int32: -2^31 <= n < 2^31
925 // Note: some of the cases tested for above fall through to this point
926 if ((q
+ exp
) <= 0) {
927 // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
934 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
935 // -2^31 <= x <= -1 or 1 <= x < 2^31 so x can be rounded
936 // toward negative infinity to a 32-bit signed integer
937 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
938 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
939 // chop off ind digits from the lower part of C1
940 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
943 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
945 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
946 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
950 // calculate C* and f*
951 // C* is actually floor(C*) in this case
952 // C* and f* need shifting and masking, as shown by
953 // shiftright128[] and maskhigh128[]
955 // kx = 10^(-x) = ten2mk128[ind - 1]
956 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
957 // the approximation of 10^(-x) was rounded up to 118 bits
958 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
959 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
960 Cstar
.w
[1] = P256
.w
[3];
961 Cstar
.w
[0] = P256
.w
[2];
963 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
964 fstar
.w
[1] = P256
.w
[1];
965 fstar
.w
[0] = P256
.w
[0];
966 } else { // 22 <= ind - 1 <= 33
968 Cstar
.w
[0] = P256
.w
[3];
969 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
970 fstar
.w
[2] = P256
.w
[2];
971 fstar
.w
[1] = P256
.w
[1];
972 fstar
.w
[0] = P256
.w
[0];
974 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
975 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
976 // if (0 < f* < 10^(-x)) then the result is a midpoint
977 // if floor(C*) is even then C* = floor(C*) - logical right
978 // shift; C* has p decimal digits, correct by Prop. 1)
979 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
980 // shift; C* has p decimal digits, correct by Pr. 1)
982 // C* = floor(C*) (logical right shift; C has p decimal digits,
983 // correct by Property 1)
986 // shift right C* by Ex-128 = shiftright128[ind]
987 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
988 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
990 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
991 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
992 } else { // 22 <= ind - 1 <= 33
993 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
995 // determine inexactness of the rounding of C*
996 // if (0 < f* - 1/2 < 10^(-x)) then
997 // the result is exact
998 // else // if (f* - 1/2 > T*) then
999 // the result is inexact
1001 if (fstar
.w
[1] > 0x8000000000000000ull
|| (fstar
.w
[1] == 0x8000000000000000ull
&& fstar
.w
[0] > 0x0ull
)) { // f* > 1/2 and the result may be exact
1002 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
1003 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
1004 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
1005 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
1006 is_inexact_lt_midpoint
= 1;
1007 } // else the result is exact
1008 } else { // the result is inexact; f2* <= 1/2
1009 is_inexact_gt_midpoint
= 1;
1011 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
1012 if (fstar
.w
[3] > 0x0 ||
1013 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
1014 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
1015 (fstar
.w
[1] || fstar
.w
[0]))) {
1016 // f2* > 1/2 and the result may be exact
1017 // Calculate f2* - 1/2
1018 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
1019 tmp64A
= fstar
.w
[3];
1020 if (tmp64
> fstar
.w
[2])
1023 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1024 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1025 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1026 is_inexact_lt_midpoint
= 1;
1027 } // else the result is exact
1028 } else { // the result is inexact; f2* <= 1/2
1029 is_inexact_gt_midpoint
= 1;
1031 } else { // if 22 <= ind <= 33
1032 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
1033 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
1034 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
1035 // f2* > 1/2 and the result may be exact
1036 // Calculate f2* - 1/2
1037 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
1038 if (tmp64
|| fstar
.w
[2]
1039 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1040 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1041 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1042 is_inexact_lt_midpoint
= 1;
1043 } // else the result is exact
1044 } else { // the result is inexact; f2* <= 1/2
1045 is_inexact_gt_midpoint
= 1;
1049 // if the result was a midpoint it was rounded away from zero, so
1050 // it will need a correction
1051 // check for midpoints
1052 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
1053 && (fstar
.w
[1] || fstar
.w
[0])
1054 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
1055 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1056 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
1057 // the result is a midpoint; round to nearest
1058 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
1059 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
1060 Cstar
.w
[0]--; // Cstar.w[0] is now even
1061 is_midpoint_gt_even
= 1;
1062 is_inexact_lt_midpoint
= 0;
1063 is_inexact_gt_midpoint
= 0;
1064 } else { // else MP in [ODD, EVEN]
1065 is_midpoint_lt_even
= 1;
1066 is_inexact_lt_midpoint
= 0;
1067 is_inexact_gt_midpoint
= 0;
1070 // general correction for RM
1071 if (x_sign
&& (is_midpoint_gt_even
|| is_inexact_lt_midpoint
)) {
1072 Cstar
.w
[0] = Cstar
.w
[0] + 1;
1074 && (is_midpoint_lt_even
|| is_inexact_gt_midpoint
)) {
1075 Cstar
.w
[0] = Cstar
.w
[0] - 1;
1077 ; // the result is already correct
1083 } else if (exp
== 0) {
1085 // res = +/-C (exact)
1090 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1091 // res = +/-C * 10^exp (exact)
1093 res
= -C1
.w
[0] * ten2k64
[exp
];
1095 res
= C1
.w
[0] * ten2k64
[exp
];
1104 /*****************************************************************************
1105 * BID128_to_int32_xfloor
1106 ****************************************************************************/
1108 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_xfloor
,
1114 int exp
; // unbiased exponent
1115 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1116 UINT64 tmp64
, tmp64A
;
1117 BID_UI64DOUBLE tmp1
;
1118 unsigned int x_nr_bits
;
1121 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
1124 int is_inexact_lt_midpoint
= 0;
1125 int is_inexact_gt_midpoint
= 0;
1126 int is_midpoint_lt_even
= 0;
1127 int is_midpoint_gt_even
= 0;
1130 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1131 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
1132 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
1135 // check for NaN or Infinity
1136 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
1138 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
1139 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
1141 *pfpsf
|= INVALID_EXCEPTION
;
1142 // return Integer Indefinite
1144 } else { // x is QNaN
1146 *pfpsf
|= INVALID_EXCEPTION
;
1147 // return Integer Indefinite
1151 } else { // x is not a NaN, so it must be infinity
1152 if (!x_sign
) { // x is +inf
1154 *pfpsf
|= INVALID_EXCEPTION
;
1155 // return Integer Indefinite
1157 } else { // x is -inf
1159 *pfpsf
|= INVALID_EXCEPTION
;
1160 // return Integer Indefinite
1166 // check for non-canonical values (after the check for special values)
1167 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
1168 || (C1
.w
[1] == 0x0001ed09bead87c0ull
1169 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
1170 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
1173 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
1177 } else { // x is not special and is not zero
1179 // q = nr. of decimal digits in x
1180 // determine first the nr. of bits in x
1182 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
1183 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1184 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
1185 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
1187 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1188 } else { // x < 2^32
1189 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
1191 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1193 } else { // if x < 2^53
1194 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
1196 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1198 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1199 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
1201 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1203 q
= nr_digits
[x_nr_bits
- 1].digits
;
1205 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1206 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
1207 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
1208 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
1211 exp
= (x_exp
>> 49) - 6176;
1212 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1214 *pfpsf
|= INVALID_EXCEPTION
;
1215 // return Integer Indefinite
1218 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1219 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1220 // so x rounded to an integer may or may not fit in a signed 32-bit int
1221 // the cases that do not fit are identified here; the ones that fit
1222 // fall through and will be handled with other cases further,
1223 // under '1 <= q + exp <= 10'
1224 if (x_sign
) { // if n < 0 and q + exp = 10
1225 // if n < -2^31 then n is too large
1226 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31
1227 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000000, 1<=q<=34
1229 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
1230 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1231 if (tmp64
> 0x500000000ull
) {
1233 *pfpsf
|= INVALID_EXCEPTION
;
1234 // return Integer Indefinite
1238 // else cases that can be rounded to a 32-bit int fall through
1239 // to '1 <= q + exp <= 10'
1240 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1241 // 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000000 <=>
1242 // C > 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 23
1244 tmp64
= 0x500000000ull
;
1245 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1246 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
1247 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1248 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
1250 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] > C
.w
[0])) {
1252 *pfpsf
|= INVALID_EXCEPTION
;
1253 // return Integer Indefinite
1257 // else cases that can be rounded to a 32-bit int fall through
1258 // to '1 <= q + exp <= 10'
1260 } else { // if n > 0 and q + exp = 10
1261 // if n >= 2^31 then n is too large
1262 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31
1263 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=34
1265 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
1266 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1267 if (tmp64
>= 0x500000000ull
) {
1269 *pfpsf
|= INVALID_EXCEPTION
;
1270 // return Integer Indefinite
1274 // else cases that can be rounded to a 32-bit int fall through
1275 // to '1 <= q + exp <= 10'
1276 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1277 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000 <=>
1278 // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 23
1280 tmp64
= 0x500000000ull
;
1281 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1282 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
1283 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1284 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
1286 if (C1
.w
[1] > C
.w
[1]
1287 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
1289 *pfpsf
|= INVALID_EXCEPTION
;
1290 // return Integer Indefinite
1294 // else cases that can be rounded to a 32-bit int fall through
1295 // to '1 <= q + exp <= 10'
1299 // n is not too large to be converted to int32: -2^31 <= n < 2^31
1300 // Note: some of the cases tested for above fall through to this point
1301 if ((q
+ exp
) <= 0) {
1302 // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
1304 *pfpsf
|= INEXACT_EXCEPTION
;
1311 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
1312 // -2^31 <= x <= -1 or 1 <= x < 2^31 so x can be rounded
1313 // toward negative infinity to a 32-bit signed integer
1314 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
1315 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
1316 // chop off ind digits from the lower part of C1
1317 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
1320 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
1322 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
1323 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
1325 if (C1
.w
[0] < tmp64
)
1327 // calculate C* and f*
1328 // C* is actually floor(C*) in this case
1329 // C* and f* need shifting and masking, as shown by
1330 // shiftright128[] and maskhigh128[]
1332 // kx = 10^(-x) = ten2mk128[ind - 1]
1333 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1334 // the approximation of 10^(-x) was rounded up to 118 bits
1335 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
1336 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1337 Cstar
.w
[1] = P256
.w
[3];
1338 Cstar
.w
[0] = P256
.w
[2];
1340 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
1341 fstar
.w
[1] = P256
.w
[1];
1342 fstar
.w
[0] = P256
.w
[0];
1343 } else { // 22 <= ind - 1 <= 33
1345 Cstar
.w
[0] = P256
.w
[3];
1346 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
1347 fstar
.w
[2] = P256
.w
[2];
1348 fstar
.w
[1] = P256
.w
[1];
1349 fstar
.w
[0] = P256
.w
[0];
1351 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1352 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1353 // if (0 < f* < 10^(-x)) then the result is a midpoint
1354 // if floor(C*) is even then C* = floor(C*) - logical right
1355 // shift; C* has p decimal digits, correct by Prop. 1)
1356 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
1357 // shift; C* has p decimal digits, correct by Pr. 1)
1359 // C* = floor(C*) (logical right shift; C has p decimal digits,
1360 // correct by Property 1)
1361 // n = C* * 10^(e+x)
1363 // shift right C* by Ex-128 = shiftright128[ind]
1364 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
1365 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1367 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
1368 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1369 } else { // 22 <= ind - 1 <= 33
1370 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
1372 // determine inexactness of the rounding of C*
1373 // if (0 < f* - 1/2 < 10^(-x)) then
1374 // the result is exact
1375 // else // if (f* - 1/2 > T*) then
1376 // the result is inexact
1378 if (fstar
.w
[1] > 0x8000000000000000ull
|| (fstar
.w
[1] == 0x8000000000000000ull
&& fstar
.w
[0] > 0x0ull
)) { // f* > 1/2 and the result may be exact
1379 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
1380 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
1381 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
1382 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
1383 // set the inexact flag
1384 *pfpsf
|= INEXACT_EXCEPTION
;
1385 is_inexact_lt_midpoint
= 1;
1386 } // else the result is exact
1387 } else { // the result is inexact; f2* <= 1/2
1388 // set the inexact flag
1389 *pfpsf
|= INEXACT_EXCEPTION
;
1390 is_inexact_gt_midpoint
= 1;
1392 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
1393 if (fstar
.w
[3] > 0x0 ||
1394 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
1395 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
1396 (fstar
.w
[1] || fstar
.w
[0]))) {
1397 // f2* > 1/2 and the result may be exact
1398 // Calculate f2* - 1/2
1399 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
1400 tmp64A
= fstar
.w
[3];
1401 if (tmp64
> fstar
.w
[2])
1404 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1405 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1406 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1407 // set the inexact flag
1408 *pfpsf
|= INEXACT_EXCEPTION
;
1409 is_inexact_lt_midpoint
= 1;
1410 } // else the result is exact
1411 } else { // the result is inexact; f2* <= 1/2
1412 // set the inexact flag
1413 *pfpsf
|= INEXACT_EXCEPTION
;
1414 is_inexact_gt_midpoint
= 1;
1416 } else { // if 22 <= ind <= 33
1417 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
1418 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
1419 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
1420 // f2* > 1/2 and the result may be exact
1421 // Calculate f2* - 1/2
1422 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
1423 if (tmp64
|| fstar
.w
[2]
1424 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1425 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1426 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1427 // set the inexact flag
1428 *pfpsf
|= INEXACT_EXCEPTION
;
1429 is_inexact_lt_midpoint
= 1;
1430 } // else the result is exact
1431 } else { // the result is inexact; f2* <= 1/2
1432 // set the inexact flag
1433 *pfpsf
|= INEXACT_EXCEPTION
;
1434 is_inexact_gt_midpoint
= 1;
1438 // if the result was a midpoint it was rounded away from zero, so
1439 // it will need a correction
1440 // check for midpoints
1441 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
1442 && (fstar
.w
[1] || fstar
.w
[0])
1443 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
1444 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1445 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
1446 // the result is a midpoint; round to nearest
1447 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
1448 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
1449 Cstar
.w
[0]--; // Cstar.w[0] is now even
1450 is_midpoint_gt_even
= 1;
1451 is_inexact_lt_midpoint
= 0;
1452 is_inexact_gt_midpoint
= 0;
1453 } else { // else MP in [ODD, EVEN]
1454 is_midpoint_lt_even
= 1;
1455 is_inexact_lt_midpoint
= 0;
1456 is_inexact_gt_midpoint
= 0;
1459 // general correction for RM
1460 if (x_sign
&& (is_midpoint_gt_even
|| is_inexact_lt_midpoint
)) {
1461 Cstar
.w
[0] = Cstar
.w
[0] + 1;
1463 && (is_midpoint_lt_even
|| is_inexact_gt_midpoint
)) {
1464 Cstar
.w
[0] = Cstar
.w
[0] - 1;
1466 ; // the result is already correct
1472 } else if (exp
== 0) {
1474 // res = +/-C (exact)
1479 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1480 // res = +/-C * 10^exp (exact)
1482 res
= -C1
.w
[0] * ten2k64
[exp
];
1484 res
= C1
.w
[0] * ten2k64
[exp
];
1492 /*****************************************************************************
1493 * BID128_to_int32_ceil
1494 ****************************************************************************/
1496 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_ceil
, x
)
1501 int exp
; // unbiased exponent
1502 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1503 UINT64 tmp64
, tmp64A
;
1504 BID_UI64DOUBLE tmp1
;
1505 unsigned int x_nr_bits
;
1508 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
1511 int is_inexact_lt_midpoint
= 0;
1512 int is_inexact_gt_midpoint
= 0;
1513 int is_midpoint_lt_even
= 0;
1514 int is_midpoint_gt_even
= 0;
1517 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1518 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
1519 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
1522 // check for NaN or Infinity
1523 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
1525 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
1526 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
1528 *pfpsf
|= INVALID_EXCEPTION
;
1529 // return Integer Indefinite
1531 } else { // x is QNaN
1533 *pfpsf
|= INVALID_EXCEPTION
;
1534 // return Integer Indefinite
1538 } else { // x is not a NaN, so it must be infinity
1539 if (!x_sign
) { // x is +inf
1541 *pfpsf
|= INVALID_EXCEPTION
;
1542 // return Integer Indefinite
1544 } else { // x is -inf
1546 *pfpsf
|= INVALID_EXCEPTION
;
1547 // return Integer Indefinite
1553 // check for non-canonical values (after the check for special values)
1554 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
1555 || (C1
.w
[1] == 0x0001ed09bead87c0ull
1556 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
1557 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
1560 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
1564 } else { // x is not special and is not zero
1566 // q = nr. of decimal digits in x
1567 // determine first the nr. of bits in x
1569 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
1570 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1571 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
1572 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
1574 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1575 } else { // x < 2^32
1576 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
1578 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1580 } else { // if x < 2^53
1581 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
1583 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1585 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1586 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
1588 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1590 q
= nr_digits
[x_nr_bits
- 1].digits
;
1592 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1593 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
1594 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
1595 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
1598 exp
= (x_exp
>> 49) - 6176;
1599 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1601 *pfpsf
|= INVALID_EXCEPTION
;
1602 // return Integer Indefinite
1605 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1606 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1607 // so x rounded to an integer may or may not fit in a signed 32-bit int
1608 // the cases that do not fit are identified here; the ones that fit
1609 // fall through and will be handled with other cases further,
1610 // under '1 <= q + exp <= 10'
1611 if (x_sign
) { // if n < 0 and q + exp = 10
1612 // if n <= -2^31-1 then n is too large
1613 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1
1614 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34
1616 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
1617 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1618 if (tmp64
>= 0x50000000aull
) {
1620 *pfpsf
|= INVALID_EXCEPTION
;
1621 // return Integer Indefinite
1625 // else cases that can be rounded to a 32-bit int fall through
1626 // to '1 <= q + exp <= 10'
1627 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1628 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a <=>
1629 // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 23
1630 // (scale 2^31+1 up)
1631 tmp64
= 0x50000000aull
;
1632 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1633 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
1634 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1635 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
1637 if (C1
.w
[1] > C
.w
[1]
1638 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
1640 *pfpsf
|= INVALID_EXCEPTION
;
1641 // return Integer Indefinite
1645 // else cases that can be rounded to a 32-bit int fall through
1646 // to '1 <= q + exp <= 10'
1648 } else { // if n > 0 and q + exp = 10
1649 // if n > 2^31 - 1 then n is too large
1650 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31 - 1
1651 // too large if 0.c(0)c(1)...c(q-1) * 10^11 > 0x4fffffff6, 1<=q<=34
1653 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
1654 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1655 if (tmp64
> 0x4fffffff6ull
) {
1657 *pfpsf
|= INVALID_EXCEPTION
;
1658 // return Integer Indefinite
1662 // else cases that can be rounded to a 32-bit int fall through
1663 // to '1 <= q + exp <= 10'
1664 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1665 // 0.c(0)c(1)...c(q-1) * 10^11 > 0x4fffffff6 <=>
1666 // C > 0x4fffffff6 * 10^(q-11) where 1 <= q - 11 <= 23
1668 tmp64
= 0x4fffffff6ull
;
1669 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1670 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
1671 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1672 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
1674 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] > C
.w
[0])) {
1676 *pfpsf
|= INVALID_EXCEPTION
;
1677 // return Integer Indefinite
1681 // else cases that can be rounded to a 32-bit int fall through
1682 // to '1 <= q + exp <= 10'
1686 // n is not too large to be converted to int32: -2^31-1 < n <= 2^31-1
1687 // Note: some of the cases tested for above fall through to this point
1688 if ((q
+ exp
) <= 0) {
1689 // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
1696 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
1697 // -2^31-1 < x <= -1 or 1 <= x <= 2^31-1 so x can be rounded
1698 // toward positive infinity to a 32-bit signed integer
1699 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
1700 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
1701 // chop off ind digits from the lower part of C1
1702 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
1705 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
1707 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
1708 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
1710 if (C1
.w
[0] < tmp64
)
1712 // calculate C* and f*
1713 // C* is actually floor(C*) in this case
1714 // C* and f* need shifting and masking, as shown by
1715 // shiftright128[] and maskhigh128[]
1717 // kx = 10^(-x) = ten2mk128[ind - 1]
1718 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1719 // the approximation of 10^(-x) was rounded up to 118 bits
1720 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
1721 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1722 Cstar
.w
[1] = P256
.w
[3];
1723 Cstar
.w
[0] = P256
.w
[2];
1725 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
1726 fstar
.w
[1] = P256
.w
[1];
1727 fstar
.w
[0] = P256
.w
[0];
1728 } else { // 22 <= ind - 1 <= 33
1730 Cstar
.w
[0] = P256
.w
[3];
1731 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
1732 fstar
.w
[2] = P256
.w
[2];
1733 fstar
.w
[1] = P256
.w
[1];
1734 fstar
.w
[0] = P256
.w
[0];
1736 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1737 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1738 // if (0 < f* < 10^(-x)) then the result is a midpoint
1739 // if floor(C*) is even then C* = floor(C*) - logical right
1740 // shift; C* has p decimal digits, correct by Prop. 1)
1741 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
1742 // shift; C* has p decimal digits, correct by Pr. 1)
1744 // C* = floor(C*) (logical right shift; C has p decimal digits,
1745 // correct by Property 1)
1746 // n = C* * 10^(e+x)
1748 // shift right C* by Ex-128 = shiftright128[ind]
1749 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
1750 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
1752 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
1753 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1754 } else { // 22 <= ind - 1 <= 33
1755 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
1757 // determine inexactness of the rounding of C*
1758 // if (0 < f* - 1/2 < 10^(-x)) then
1759 // the result is exact
1760 // else // if (f* - 1/2 > T*) then
1761 // the result is inexact
1763 if (fstar
.w
[1] > 0x8000000000000000ull
|| (fstar
.w
[1] == 0x8000000000000000ull
&& fstar
.w
[0] > 0x0ull
)) { // f* > 1/2 and the result may be exact
1764 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
1765 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
1766 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
1767 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
1768 is_inexact_lt_midpoint
= 1;
1769 } // else the result is exact
1770 } else { // the result is inexact; f2* <= 1/2
1771 is_inexact_gt_midpoint
= 1;
1773 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
1774 if (fstar
.w
[3] > 0x0 ||
1775 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
1776 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
1777 (fstar
.w
[1] || fstar
.w
[0]))) {
1778 // f2* > 1/2 and the result may be exact
1779 // Calculate f2* - 1/2
1780 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
1781 tmp64A
= fstar
.w
[3];
1782 if (tmp64
> fstar
.w
[2])
1785 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1786 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1787 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1788 is_inexact_lt_midpoint
= 1;
1789 } // else the result is exact
1790 } else { // the result is inexact; f2* <= 1/2
1791 is_inexact_gt_midpoint
= 1;
1793 } else { // if 22 <= ind <= 33
1794 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
1795 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
1796 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
1797 // f2* > 1/2 and the result may be exact
1798 // Calculate f2* - 1/2
1799 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
1800 if (tmp64
|| fstar
.w
[2]
1801 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
1802 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1803 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
1804 is_inexact_lt_midpoint
= 1;
1805 } // else the result is exact
1806 } else { // the result is inexact; f2* <= 1/2
1807 is_inexact_gt_midpoint
= 1;
1811 // if the result was a midpoint it was rounded away from zero, so
1812 // it will need a correction
1813 // check for midpoints
1814 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
1815 && (fstar
.w
[1] || fstar
.w
[0])
1816 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
1817 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
1818 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
1819 // the result is a midpoint; round to nearest
1820 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
1821 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
1822 Cstar
.w
[0]--; // Cstar.w[0] is now even
1823 is_midpoint_gt_even
= 1;
1824 is_inexact_lt_midpoint
= 0;
1825 is_inexact_gt_midpoint
= 0;
1826 } else { // else MP in [ODD, EVEN]
1827 is_midpoint_lt_even
= 1;
1828 is_inexact_lt_midpoint
= 0;
1829 is_inexact_gt_midpoint
= 0;
1832 // general correction for RM
1833 if (x_sign
&& (is_midpoint_lt_even
|| is_inexact_gt_midpoint
)) {
1834 Cstar
.w
[0] = Cstar
.w
[0] - 1;
1836 && (is_midpoint_gt_even
|| is_inexact_lt_midpoint
)) {
1837 Cstar
.w
[0] = Cstar
.w
[0] + 1;
1839 ; // the result is already correct
1845 } else if (exp
== 0) {
1847 // res = +/-C (exact)
1852 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1853 // res = +/-C * 10^exp (exact)
1855 res
= -C1
.w
[0] * ten2k64
[exp
];
1857 res
= C1
.w
[0] * ten2k64
[exp
];
1865 /*****************************************************************************
1866 * BID128_to_int32_xceil
1867 ****************************************************************************/
1869 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_xceil
, x
)
1874 int exp
; // unbiased exponent
1875 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1876 UINT64 tmp64
, tmp64A
;
1877 BID_UI64DOUBLE tmp1
;
1878 unsigned int x_nr_bits
;
1881 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
1884 int is_inexact_lt_midpoint
= 0;
1885 int is_inexact_gt_midpoint
= 0;
1886 int is_midpoint_lt_even
= 0;
1887 int is_midpoint_gt_even
= 0;
1890 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
1891 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
1892 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
1895 // check for NaN or Infinity
1896 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
1898 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
1899 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
1901 *pfpsf
|= INVALID_EXCEPTION
;
1902 // return Integer Indefinite
1904 } else { // x is QNaN
1906 *pfpsf
|= INVALID_EXCEPTION
;
1907 // return Integer Indefinite
1911 } else { // x is not a NaN, so it must be infinity
1912 if (!x_sign
) { // x is +inf
1914 *pfpsf
|= INVALID_EXCEPTION
;
1915 // return Integer Indefinite
1917 } else { // x is -inf
1919 *pfpsf
|= INVALID_EXCEPTION
;
1920 // return Integer Indefinite
1926 // check for non-canonical values (after the check for special values)
1927 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
1928 || (C1
.w
[1] == 0x0001ed09bead87c0ull
1929 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
1930 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
1933 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
1937 } else { // x is not special and is not zero
1939 // q = nr. of decimal digits in x
1940 // determine first the nr. of bits in x
1942 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
1943 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1944 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
1945 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
1947 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1948 } else { // x < 2^32
1949 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
1951 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1953 } else { // if x < 2^53
1954 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
1956 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1958 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1959 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
1961 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
1963 q
= nr_digits
[x_nr_bits
- 1].digits
;
1965 q
= nr_digits
[x_nr_bits
- 1].digits1
;
1966 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
1967 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
1968 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
1971 exp
= (x_exp
>> 49) - 6176;
1972 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1974 *pfpsf
|= INVALID_EXCEPTION
;
1975 // return Integer Indefinite
1978 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1979 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1980 // so x rounded to an integer may or may not fit in a signed 32-bit int
1981 // the cases that do not fit are identified here; the ones that fit
1982 // fall through and will be handled with other cases further,
1983 // under '1 <= q + exp <= 10'
1984 if (x_sign
) { // if n < 0 and q + exp = 10
1985 // if n <= -2^31-1 then n is too large
1986 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1
1987 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34
1989 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
1990 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1991 if (tmp64
>= 0x50000000aull
) {
1993 *pfpsf
|= INVALID_EXCEPTION
;
1994 // return Integer Indefinite
1998 // else cases that can be rounded to a 32-bit int fall through
1999 // to '1 <= q + exp <= 10'
2000 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2001 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a <=>
2002 // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 23
2003 // (scale 2^31+1 up)
2004 tmp64
= 0x50000000aull
;
2005 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2006 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
2007 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2008 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
2010 if (C1
.w
[1] > C
.w
[1]
2011 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
2013 *pfpsf
|= INVALID_EXCEPTION
;
2014 // return Integer Indefinite
2018 // else cases that can be rounded to a 32-bit int fall through
2019 // to '1 <= q + exp <= 10'
2021 } else { // if n > 0 and q + exp = 10
2022 // if n > 2^31 - 1 then n is too large
2023 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31 - 1
2024 // too large if 0.c(0)c(1)...c(q-1) * 10^11 > 0x4fffffff6, 1<=q<=34
2026 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
2027 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2028 if (tmp64
> 0x4fffffff6ull
) {
2030 *pfpsf
|= INVALID_EXCEPTION
;
2031 // return Integer Indefinite
2035 // else cases that can be rounded to a 32-bit int fall through
2036 // to '1 <= q + exp <= 10'
2037 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2038 // 0.c(0)c(1)...c(q-1) * 10^11 > 0x4fffffff6 <=>
2039 // C > 0x4fffffff6 * 10^(q-11) where 1 <= q - 11 <= 23
2041 tmp64
= 0x4fffffff6ull
;
2042 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2043 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
2044 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2045 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
2047 if (C1
.w
[1] > C
.w
[1] || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] > C
.w
[0])) {
2049 *pfpsf
|= INVALID_EXCEPTION
;
2050 // return Integer Indefinite
2054 // else cases that can be rounded to a 32-bit int fall through
2055 // to '1 <= q + exp <= 10'
2059 // n is not too large to be converted to int32: -2^31-1 < n <= 2^31-1
2060 // Note: some of the cases tested for above fall through to this point
2061 if ((q
+ exp
) <= 0) {
2062 // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
2064 *pfpsf
|= INEXACT_EXCEPTION
;
2071 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
2072 // -2^31-1 < x <= -1 or 1 <= x <= 2^31-1 so x can be rounded
2073 // toward positive infinity to a 32-bit signed integer
2074 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
2075 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
2076 // chop off ind digits from the lower part of C1
2077 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
2080 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
2082 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
2083 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
2085 if (C1
.w
[0] < tmp64
)
2087 // calculate C* and f*
2088 // C* is actually floor(C*) in this case
2089 // C* and f* need shifting and masking, as shown by
2090 // shiftright128[] and maskhigh128[]
2092 // kx = 10^(-x) = ten2mk128[ind - 1]
2093 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2094 // the approximation of 10^(-x) was rounded up to 118 bits
2095 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
2096 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2097 Cstar
.w
[1] = P256
.w
[3];
2098 Cstar
.w
[0] = P256
.w
[2];
2100 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
2101 fstar
.w
[1] = P256
.w
[1];
2102 fstar
.w
[0] = P256
.w
[0];
2103 } else { // 22 <= ind - 1 <= 33
2105 Cstar
.w
[0] = P256
.w
[3];
2106 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
2107 fstar
.w
[2] = P256
.w
[2];
2108 fstar
.w
[1] = P256
.w
[1];
2109 fstar
.w
[0] = P256
.w
[0];
2111 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2112 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2113 // if (0 < f* < 10^(-x)) then the result is a midpoint
2114 // if floor(C*) is even then C* = floor(C*) - logical right
2115 // shift; C* has p decimal digits, correct by Prop. 1)
2116 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
2117 // shift; C* has p decimal digits, correct by Pr. 1)
2119 // C* = floor(C*) (logical right shift; C has p decimal digits,
2120 // correct by Property 1)
2121 // n = C* * 10^(e+x)
2123 // shift right C* by Ex-128 = shiftright128[ind]
2124 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
2125 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2127 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
2128 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2129 } else { // 22 <= ind - 1 <= 33
2130 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
2132 // determine inexactness of the rounding of C*
2133 // if (0 < f* - 1/2 < 10^(-x)) then
2134 // the result is exact
2135 // else // if (f* - 1/2 > T*) then
2136 // the result is inexact
2138 if (fstar
.w
[1] > 0x8000000000000000ull
|| (fstar
.w
[1] == 0x8000000000000000ull
&& fstar
.w
[0] > 0x0ull
)) { // f* > 1/2 and the result may be exact
2139 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
2140 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
2141 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
2142 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
2143 // set the inexact flag
2144 *pfpsf
|= INEXACT_EXCEPTION
;
2145 is_inexact_lt_midpoint
= 1;
2146 } // else the result is exact
2147 } else { // the result is inexact; f2* <= 1/2
2148 // set the inexact flag
2149 *pfpsf
|= INEXACT_EXCEPTION
;
2150 is_inexact_gt_midpoint
= 1;
2152 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
2153 if (fstar
.w
[3] > 0x0 ||
2154 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
2155 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
2156 (fstar
.w
[1] || fstar
.w
[0]))) {
2157 // f2* > 1/2 and the result may be exact
2158 // Calculate f2* - 1/2
2159 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
2160 tmp64A
= fstar
.w
[3];
2161 if (tmp64
> fstar
.w
[2])
2164 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2165 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2166 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2167 // set the inexact flag
2168 *pfpsf
|= INEXACT_EXCEPTION
;
2169 is_inexact_lt_midpoint
= 1;
2170 } // else the result is exact
2171 } else { // the result is inexact; f2* <= 1/2
2172 // set the inexact flag
2173 *pfpsf
|= INEXACT_EXCEPTION
;
2174 is_inexact_gt_midpoint
= 1;
2176 } else { // if 22 <= ind <= 33
2177 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
2178 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
2179 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
2180 // f2* > 1/2 and the result may be exact
2181 // Calculate f2* - 1/2
2182 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
2183 if (tmp64
|| fstar
.w
[2]
2184 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2185 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2186 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2187 // set the inexact flag
2188 *pfpsf
|= INEXACT_EXCEPTION
;
2189 is_inexact_lt_midpoint
= 1;
2190 } // else the result is exact
2191 } else { // the result is inexact; f2* <= 1/2
2192 // set the inexact flag
2193 *pfpsf
|= INEXACT_EXCEPTION
;
2194 is_inexact_gt_midpoint
= 1;
2198 // if the result was a midpoint it was rounded away from zero, so
2199 // it will need a correction
2200 // check for midpoints
2201 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
2202 && (fstar
.w
[1] || fstar
.w
[0])
2203 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
2204 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2205 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
2206 // the result is a midpoint; round to nearest
2207 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
2208 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
2209 Cstar
.w
[0]--; // Cstar.w[0] is now even
2210 is_midpoint_gt_even
= 1;
2211 is_inexact_lt_midpoint
= 0;
2212 is_inexact_gt_midpoint
= 0;
2213 } else { // else MP in [ODD, EVEN]
2214 is_midpoint_lt_even
= 1;
2215 is_inexact_lt_midpoint
= 0;
2216 is_inexact_gt_midpoint
= 0;
2219 // general correction for RM
2220 if (x_sign
&& (is_midpoint_lt_even
|| is_inexact_gt_midpoint
)) {
2221 Cstar
.w
[0] = Cstar
.w
[0] - 1;
2223 && (is_midpoint_gt_even
|| is_inexact_lt_midpoint
)) {
2224 Cstar
.w
[0] = Cstar
.w
[0] + 1;
2226 ; // the result is already correct
2232 } else if (exp
== 0) {
2234 // res = +/-C (exact)
2239 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2240 // res = +/-C * 10^exp (exact)
2242 res
= -C1
.w
[0] * ten2k64
[exp
];
2244 res
= C1
.w
[0] * ten2k64
[exp
];
2252 /*****************************************************************************
2253 * BID128_to_int32_int
2254 ****************************************************************************/
2256 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_int
, x
)
2261 int exp
; // unbiased exponent
2262 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2263 UINT64 tmp64
, tmp64A
;
2264 BID_UI64DOUBLE tmp1
;
2265 unsigned int x_nr_bits
;
2268 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
2271 int is_inexact_gt_midpoint
= 0;
2272 int is_midpoint_lt_even
= 0;
2275 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
2276 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
2277 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
2280 // check for NaN or Infinity
2281 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
2283 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
2284 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
2286 *pfpsf
|= INVALID_EXCEPTION
;
2287 // return Integer Indefinite
2289 } else { // x is QNaN
2291 *pfpsf
|= INVALID_EXCEPTION
;
2292 // return Integer Indefinite
2296 } else { // x is not a NaN, so it must be infinity
2297 if (!x_sign
) { // x is +inf
2299 *pfpsf
|= INVALID_EXCEPTION
;
2300 // return Integer Indefinite
2302 } else { // x is -inf
2304 *pfpsf
|= INVALID_EXCEPTION
;
2305 // return Integer Indefinite
2311 // check for non-canonical values (after the check for special values)
2312 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
2313 || (C1
.w
[1] == 0x0001ed09bead87c0ull
2314 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
2315 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
2318 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
2322 } else { // x is not special and is not zero
2324 // q = nr. of decimal digits in x
2325 // determine first the nr. of bits in x
2327 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
2328 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2329 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
2330 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
2332 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2333 } else { // x < 2^32
2334 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
2336 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2338 } else { // if x < 2^53
2339 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
2341 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2343 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2344 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
2346 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2348 q
= nr_digits
[x_nr_bits
- 1].digits
;
2350 q
= nr_digits
[x_nr_bits
- 1].digits1
;
2351 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
2352 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
2353 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
2356 exp
= (x_exp
>> 49) - 6176;
2357 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
2359 *pfpsf
|= INVALID_EXCEPTION
;
2360 // return Integer Indefinite
2363 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
2364 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
2365 // so x rounded to an integer may or may not fit in a signed 32-bit int
2366 // the cases that do not fit are identified here; the ones that fit
2367 // fall through and will be handled with other cases further,
2368 // under '1 <= q + exp <= 10'
2369 if (x_sign
) { // if n < 0 and q + exp = 10
2370 // if n <= -2^31 - 1 then n is too large
2371 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1
2372 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34
2374 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
2375 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2376 if (tmp64
>= 0x50000000aull
) {
2378 *pfpsf
|= INVALID_EXCEPTION
;
2379 // return Integer Indefinite
2383 // else cases that can be rounded to a 32-bit int fall through
2384 // to '1 <= q + exp <= 10'
2385 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2386 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a <=>
2387 // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 23
2388 // (scale 2^31+1 up)
2389 tmp64
= 0x50000000aull
;
2390 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2391 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
2392 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2393 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
2395 if (C1
.w
[1] > C
.w
[1]
2396 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
2398 *pfpsf
|= INVALID_EXCEPTION
;
2399 // return Integer Indefinite
2403 // else cases that can be rounded to a 32-bit int fall through
2404 // to '1 <= q + exp <= 10'
2406 } else { // if n > 0 and q + exp = 10
2407 // if n >= 2^31 then n is too large
2408 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31
2409 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=34
2411 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
2412 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2413 if (tmp64
>= 0x500000000ull
) {
2415 *pfpsf
|= INVALID_EXCEPTION
;
2416 // return Integer Indefinite
2420 // else cases that can be rounded to a 32-bit int fall through
2421 // to '1 <= q + exp <= 10'
2422 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2423 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000 <=>
2424 // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 23
2425 // (scale 2^31-1/2 up)
2426 tmp64
= 0x500000000ull
;
2427 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2428 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
2429 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2430 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
2432 if (C1
.w
[1] > C
.w
[1]
2433 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
2435 *pfpsf
|= INVALID_EXCEPTION
;
2436 // return Integer Indefinite
2440 // else cases that can be rounded to a 32-bit int fall through
2441 // to '1 <= q + exp <= 10'
2445 // n is not too large to be converted to int32: -2^31 - 1 < n < 2^31
2446 // Note: some of the cases tested for above fall through to this point
2447 if ((q
+ exp
) <= 0) {
2448 // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
2452 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
2453 // -2^31-1 < x <= -1 or 1 <= x < 2^31 so x can be rounded
2454 // toward zero to a 32-bit signed integer
2455 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
2456 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
2457 // chop off ind digits from the lower part of C1
2458 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
2461 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
2463 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
2464 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
2466 if (C1
.w
[0] < tmp64
)
2468 // calculate C* and f*
2469 // C* is actually floor(C*) in this case
2470 // C* and f* need shifting and masking, as shown by
2471 // shiftright128[] and maskhigh128[]
2473 // kx = 10^(-x) = ten2mk128[ind - 1]
2474 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2475 // the approximation of 10^(-x) was rounded up to 118 bits
2476 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
2477 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2478 Cstar
.w
[1] = P256
.w
[3];
2479 Cstar
.w
[0] = P256
.w
[2];
2481 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
2482 fstar
.w
[1] = P256
.w
[1];
2483 fstar
.w
[0] = P256
.w
[0];
2484 } else { // 22 <= ind - 1 <= 33
2486 Cstar
.w
[0] = P256
.w
[3];
2487 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
2488 fstar
.w
[2] = P256
.w
[2];
2489 fstar
.w
[1] = P256
.w
[1];
2490 fstar
.w
[0] = P256
.w
[0];
2492 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2493 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2494 // if (0 < f* < 10^(-x)) then the result is a midpoint
2495 // if floor(C*) is even then C* = floor(C*) - logical right
2496 // shift; C* has p decimal digits, correct by Prop. 1)
2497 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
2498 // shift; C* has p decimal digits, correct by Pr. 1)
2500 // C* = floor(C*) (logical right shift; C has p decimal digits,
2501 // correct by Property 1)
2502 // n = C* * 10^(e+x)
2504 // shift right C* by Ex-128 = shiftright128[ind]
2505 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
2506 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2508 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
2509 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2510 } else { // 22 <= ind - 1 <= 33
2511 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
2513 // determine inexactness of the rounding of C*
2514 // if (0 < f* - 1/2 < 10^(-x)) then
2515 // the result is exact
2516 // else // if (f* - 1/2 > T*) then
2517 // the result is inexact
2519 if (fstar
.w
[1] > 0x8000000000000000ull
|| (fstar
.w
[1] == 0x8000000000000000ull
&& fstar
.w
[0] > 0x0ull
)) { // f* > 1/2 and the result may be exact
2520 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
2521 if ((tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
2522 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
2523 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0]))) {
2524 } // else the result is exact
2525 } else { // the result is inexact; f2* <= 1/2
2526 is_inexact_gt_midpoint
= 1;
2528 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
2529 if (fstar
.w
[3] > 0x0 ||
2530 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
2531 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
2532 (fstar
.w
[1] || fstar
.w
[0]))) {
2533 // f2* > 1/2 and the result may be exact
2534 // Calculate f2* - 1/2
2535 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
2536 tmp64A
= fstar
.w
[3];
2537 if (tmp64
> fstar
.w
[2])
2540 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2541 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2542 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2543 } // else the result is exact
2544 } else { // the result is inexact; f2* <= 1/2
2545 is_inexact_gt_midpoint
= 1;
2547 } else { // if 22 <= ind <= 33
2548 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
2549 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
2550 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
2551 // f2* > 1/2 and the result may be exact
2552 // Calculate f2* - 1/2
2553 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
2554 if (tmp64
|| fstar
.w
[2]
2555 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2556 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2557 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2558 } // else the result is exact
2559 } else { // the result is inexact; f2* <= 1/2
2560 is_inexact_gt_midpoint
= 1;
2564 // if the result was a midpoint it was rounded away from zero, so
2565 // it will need a correction
2566 // check for midpoints
2567 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0) &&
2568 (fstar
.w
[1] || fstar
.w
[0]) &&
2569 (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1] ||
2570 (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1] &&
2571 fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
2572 // the result is a midpoint; round to nearest
2573 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
2574 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
2575 Cstar
.w
[0]--; // Cstar.w[0] is now even
2576 is_inexact_gt_midpoint
= 0;
2577 } else { // else MP in [ODD, EVEN]
2578 is_midpoint_lt_even
= 1;
2579 is_inexact_gt_midpoint
= 0;
2582 // general correction for RZ
2583 if (is_midpoint_lt_even
|| is_inexact_gt_midpoint
) {
2584 Cstar
.w
[0] = Cstar
.w
[0] - 1;
2586 ; // exact, the result is already correct
2592 } else if (exp
== 0) {
2594 // res = +/-C (exact)
2599 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2600 // res = +/-C * 10^exp (exact)
2602 res
= -C1
.w
[0] * ten2k64
[exp
];
2604 res
= C1
.w
[0] * ten2k64
[exp
];
2612 /*****************************************************************************
2613 * BID128_to_int32_xint
2614 ****************************************************************************/
2616 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_xint
, x
)
2621 int exp
; // unbiased exponent
2622 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2623 UINT64 tmp64
, tmp64A
;
2624 BID_UI64DOUBLE tmp1
;
2625 unsigned int x_nr_bits
;
2628 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
2631 int is_inexact_gt_midpoint
= 0;
2632 int is_midpoint_lt_even
= 0;
2635 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
2636 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
2637 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
2640 // check for NaN or Infinity
2641 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
2643 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
2644 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
2646 *pfpsf
|= INVALID_EXCEPTION
;
2647 // return Integer Indefinite
2649 } else { // x is QNaN
2651 *pfpsf
|= INVALID_EXCEPTION
;
2652 // return Integer Indefinite
2656 } else { // x is not a NaN, so it must be infinity
2657 if (!x_sign
) { // x is +inf
2659 *pfpsf
|= INVALID_EXCEPTION
;
2660 // return Integer Indefinite
2662 } else { // x is -inf
2664 *pfpsf
|= INVALID_EXCEPTION
;
2665 // return Integer Indefinite
2671 // check for non-canonical values (after the check for special values)
2672 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
2673 || (C1
.w
[1] == 0x0001ed09bead87c0ull
2674 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
2675 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
2678 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
2682 } else { // x is not special and is not zero
2684 // q = nr. of decimal digits in x
2685 // determine first the nr. of bits in x
2687 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
2688 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2689 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
2690 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
2692 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2693 } else { // x < 2^32
2694 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
2696 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2698 } else { // if x < 2^53
2699 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
2701 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2703 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2704 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
2706 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
2708 q
= nr_digits
[x_nr_bits
- 1].digits
;
2710 q
= nr_digits
[x_nr_bits
- 1].digits1
;
2711 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
2712 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
2713 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
2716 exp
= (x_exp
>> 49) - 6176;
2717 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
2719 *pfpsf
|= INVALID_EXCEPTION
;
2720 // return Integer Indefinite
2723 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
2724 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
2725 // so x rounded to an integer may or may not fit in a signed 32-bit int
2726 // the cases that do not fit are identified here; the ones that fit
2727 // fall through and will be handled with other cases further,
2728 // under '1 <= q + exp <= 10'
2729 if (x_sign
) { // if n < 0 and q + exp = 10
2730 // if n <= -2^31 - 1 then n is too large
2731 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1
2732 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34
2734 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
2735 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2736 if (tmp64
>= 0x50000000aull
) {
2738 *pfpsf
|= INVALID_EXCEPTION
;
2739 // return Integer Indefinite
2743 // else cases that can be rounded to a 32-bit int fall through
2744 // to '1 <= q + exp <= 10'
2745 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2746 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a <=>
2747 // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 23
2748 // (scale 2^31+1 up)
2749 tmp64
= 0x50000000aull
;
2750 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2751 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
2752 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2753 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
2755 if (C1
.w
[1] > C
.w
[1]
2756 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
2758 *pfpsf
|= INVALID_EXCEPTION
;
2759 // return Integer Indefinite
2763 // else cases that can be rounded to a 32-bit int fall through
2764 // to '1 <= q + exp <= 10'
2766 } else { // if n > 0 and q + exp = 10
2767 // if n >= 2^31 then n is too large
2768 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31
2769 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=34
2771 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
2772 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2773 if (tmp64
>= 0x500000000ull
) {
2775 *pfpsf
|= INVALID_EXCEPTION
;
2776 // return Integer Indefinite
2780 // else cases that can be rounded to a 32-bit int fall through
2781 // to '1 <= q + exp <= 10'
2782 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2783 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000 <=>
2784 // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 23
2785 // (scale 2^31-1/2 up)
2786 tmp64
= 0x500000000ull
;
2787 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2788 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
2789 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2790 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
2792 if (C1
.w
[1] > C
.w
[1]
2793 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
2795 *pfpsf
|= INVALID_EXCEPTION
;
2796 // return Integer Indefinite
2800 // else cases that can be rounded to a 32-bit int fall through
2801 // to '1 <= q + exp <= 10'
2805 // n is not too large to be converted to int32: -2^31 - 1 < n < 2^31
2806 // Note: some of the cases tested for above fall through to this point
2807 if ((q
+ exp
) <= 0) {
2808 // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
2810 *pfpsf
|= INEXACT_EXCEPTION
;
2814 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
2815 // -2^31-1 < x <= -1 or 1 <= x < 2^31 so x can be rounded
2816 // toward zero to a 32-bit signed integer
2817 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
2818 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
2819 // chop off ind digits from the lower part of C1
2820 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
2823 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
2825 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
2826 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
2828 if (C1
.w
[0] < tmp64
)
2830 // calculate C* and f*
2831 // C* is actually floor(C*) in this case
2832 // C* and f* need shifting and masking, as shown by
2833 // shiftright128[] and maskhigh128[]
2835 // kx = 10^(-x) = ten2mk128[ind - 1]
2836 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2837 // the approximation of 10^(-x) was rounded up to 118 bits
2838 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
2839 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2840 Cstar
.w
[1] = P256
.w
[3];
2841 Cstar
.w
[0] = P256
.w
[2];
2843 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
2844 fstar
.w
[1] = P256
.w
[1];
2845 fstar
.w
[0] = P256
.w
[0];
2846 } else { // 22 <= ind - 1 <= 33
2848 Cstar
.w
[0] = P256
.w
[3];
2849 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
2850 fstar
.w
[2] = P256
.w
[2];
2851 fstar
.w
[1] = P256
.w
[1];
2852 fstar
.w
[0] = P256
.w
[0];
2854 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2855 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2856 // if (0 < f* < 10^(-x)) then the result is a midpoint
2857 // if floor(C*) is even then C* = floor(C*) - logical right
2858 // shift; C* has p decimal digits, correct by Prop. 1)
2859 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
2860 // shift; C* has p decimal digits, correct by Pr. 1)
2862 // C* = floor(C*) (logical right shift; C has p decimal digits,
2863 // correct by Property 1)
2864 // n = C* * 10^(e+x)
2866 // shift right C* by Ex-128 = shiftright128[ind]
2867 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
2868 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
2870 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
2871 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2872 } else { // 22 <= ind - 1 <= 33
2873 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
2875 // determine inexactness of the rounding of C*
2876 // if (0 < f* - 1/2 < 10^(-x)) then
2877 // the result is exact
2878 // else // if (f* - 1/2 > T*) then
2879 // the result is inexact
2881 if (fstar
.w
[1] > 0x8000000000000000ull
|| (fstar
.w
[1] == 0x8000000000000000ull
&& fstar
.w
[0] > 0x0ull
)) { // f* > 1/2 and the result may be exact
2882 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
2883 if (tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
2884 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
2885 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0])) {
2886 // set the inexact flag
2887 *pfpsf
|= INEXACT_EXCEPTION
;
2888 } // else the result is exact
2889 } else { // the result is inexact; f2* <= 1/2
2890 // set the inexact flag
2891 *pfpsf
|= INEXACT_EXCEPTION
;
2892 is_inexact_gt_midpoint
= 1;
2894 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
2895 if (fstar
.w
[3] > 0x0 ||
2896 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
2897 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
2898 (fstar
.w
[1] || fstar
.w
[0]))) {
2899 // f2* > 1/2 and the result may be exact
2900 // Calculate f2* - 1/2
2901 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
2902 tmp64A
= fstar
.w
[3];
2903 if (tmp64
> fstar
.w
[2])
2906 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2907 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2908 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2909 // set the inexact flag
2910 *pfpsf
|= INEXACT_EXCEPTION
;
2911 } // else the result is exact
2912 } else { // the result is inexact; f2* <= 1/2
2913 // set the inexact flag
2914 *pfpsf
|= INEXACT_EXCEPTION
;
2915 is_inexact_gt_midpoint
= 1;
2917 } else { // if 22 <= ind <= 33
2918 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
2919 (fstar
.w
[3] == onehalf128
[ind
- 1] && (fstar
.w
[2] ||
2922 // f2* > 1/2 and the result may be exact
2923 // Calculate f2* - 1/2
2924 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
2925 if (tmp64
|| fstar
.w
[2]
2926 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
2927 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2928 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
2929 // set the inexact flag
2930 *pfpsf
|= INEXACT_EXCEPTION
;
2931 } // else the result is exact
2932 } else { // the result is inexact; f2* <= 1/2
2933 // set the inexact flag
2934 *pfpsf
|= INEXACT_EXCEPTION
;
2935 is_inexact_gt_midpoint
= 1;
2939 // if the result was a midpoint it was rounded away from zero, so
2940 // it will need a correction
2941 // check for midpoints
2942 if ((fstar
.w
[3] == 0) && (fstar
.w
[2] == 0)
2943 && (fstar
.w
[1] || fstar
.w
[0])
2944 && (fstar
.w
[1] < ten2mk128trunc
[ind
- 1].w
[1]
2945 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
2946 && fstar
.w
[0] <= ten2mk128trunc
[ind
- 1].w
[0]))) {
2947 // the result is a midpoint; round to nearest
2948 if (Cstar
.w
[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
2949 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
2950 Cstar
.w
[0]--; // Cstar.w[0] is now even
2951 is_inexact_gt_midpoint
= 0;
2952 } else { // else MP in [ODD, EVEN]
2953 is_midpoint_lt_even
= 1;
2954 is_inexact_gt_midpoint
= 0;
2957 // general correction for RZ
2958 if (is_midpoint_lt_even
|| is_inexact_gt_midpoint
) {
2959 Cstar
.w
[0] = Cstar
.w
[0] - 1;
2961 ; // exact, the result is already correct
2967 } else if (exp
== 0) {
2969 // res = +/-C (exact)
2974 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2975 // res = +/-C * 10^exp (exact)
2977 res
= -C1
.w
[0] * ten2k64
[exp
];
2979 res
= C1
.w
[0] * ten2k64
[exp
];
2987 /*****************************************************************************
2988 * BID128_to_int32_rninta
2989 ****************************************************************************/
2991 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_rninta
,
2997 int exp
; // unbiased exponent
2998 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
3000 BID_UI64DOUBLE tmp1
;
3001 unsigned int x_nr_bits
;
3004 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
3008 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
3009 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
3010 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
3013 // check for NaN or Infinity
3014 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
3016 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
3017 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
3019 *pfpsf
|= INVALID_EXCEPTION
;
3020 // return Integer Indefinite
3022 } else { // x is QNaN
3024 *pfpsf
|= INVALID_EXCEPTION
;
3025 // return Integer Indefinite
3029 } else { // x is not a NaN, so it must be infinity
3030 if (!x_sign
) { // x is +inf
3032 *pfpsf
|= INVALID_EXCEPTION
;
3033 // return Integer Indefinite
3035 } else { // x is -inf
3037 *pfpsf
|= INVALID_EXCEPTION
;
3038 // return Integer Indefinite
3044 // check for non-canonical values (after the check for special values)
3045 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
3046 || (C1
.w
[1] == 0x0001ed09bead87c0ull
3047 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
3048 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
3051 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
3055 } else { // x is not special and is not zero
3057 // q = nr. of decimal digits in x
3058 // determine first the nr. of bits in x
3060 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
3061 // split the 64-bit value in two 32-bit halves to avoid rounding errors
3062 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
3063 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
3065 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3066 } else { // x < 2^32
3067 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
3069 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3071 } else { // if x < 2^53
3072 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
3074 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3076 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
3077 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
3079 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3081 q
= nr_digits
[x_nr_bits
- 1].digits
;
3083 q
= nr_digits
[x_nr_bits
- 1].digits1
;
3084 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
3085 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
3086 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
3089 exp
= (x_exp
>> 49) - 6176;
3090 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
3092 *pfpsf
|= INVALID_EXCEPTION
;
3093 // return Integer Indefinite
3096 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
3097 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
3098 // so x rounded to an integer may or may not fit in a signed 32-bit int
3099 // the cases that do not fit are identified here; the ones that fit
3100 // fall through and will be handled with other cases further,
3101 // under '1 <= q + exp <= 10'
3102 if (x_sign
) { // if n < 0 and q + exp = 10
3103 // if n <= -2^31 - 1/2 then n is too large
3104 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1/2
3105 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000005, 1<=q<=34
3107 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
3108 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3109 if (tmp64
>= 0x500000005ull
) {
3111 *pfpsf
|= INVALID_EXCEPTION
;
3112 // return Integer Indefinite
3116 // else cases that can be rounded to a 32-bit int fall through
3117 // to '1 <= q + exp <= 10'
3118 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3119 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000005 <=>
3120 // C >= 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 23
3121 // (scale 2^31+1/2 up)
3122 tmp64
= 0x500000005ull
;
3123 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3124 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
3125 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3126 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
3128 if (C1
.w
[1] > C
.w
[1]
3129 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
3131 *pfpsf
|= INVALID_EXCEPTION
;
3132 // return Integer Indefinite
3136 // else cases that can be rounded to a 32-bit int fall through
3137 // to '1 <= q + exp <= 10'
3139 } else { // if n > 0 and q + exp = 10
3140 // if n >= 2^31 - 1/2 then n is too large
3141 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2
3142 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=34
3144 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
3145 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3146 if (tmp64
>= 0x4fffffffbull
) {
3148 *pfpsf
|= INVALID_EXCEPTION
;
3149 // return Integer Indefinite
3153 // else cases that can be rounded to a 32-bit int fall through
3154 // to '1 <= q + exp <= 10'
3155 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3156 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb <=>
3157 // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
3158 // (scale 2^31-1/2 up)
3159 tmp64
= 0x4fffffffbull
;
3160 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3161 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
3162 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3163 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
3165 if (C1
.w
[1] > C
.w
[1]
3166 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
3168 *pfpsf
|= INVALID_EXCEPTION
;
3169 // return Integer Indefinite
3173 // else cases that can be rounded to a 32-bit int fall through
3174 // to '1 <= q + exp <= 10'
3178 // n is not too large to be converted to int32: -2^31 - 1/2 < n < 2^31 - 1/2
3179 // Note: some of the cases tested for above fall through to this point
3180 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
3184 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
3185 // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
3190 if (ind
<= 18) { // 0 <= ind <= 18
3191 if ((C1
.w
[1] == 0) && (C1
.w
[0] < midpoint64
[ind
])) {
3192 res
= 0x00000000; // return 0
3193 } else if (x_sign
) { // n < 0
3194 res
= 0xffffffff; // return -1
3196 res
= 0x00000001; // return +1
3198 } else { // 19 <= ind <= 33
3199 if ((C1
.w
[1] < midpoint128
[ind
- 19].w
[1])
3200 || ((C1
.w
[1] == midpoint128
[ind
- 19].w
[1])
3201 && (C1
.w
[0] < midpoint128
[ind
- 19].w
[0]))) {
3202 res
= 0x00000000; // return 0
3203 } else if (x_sign
) { // n < 0
3204 res
= 0xffffffff; // return -1
3206 res
= 0x00000001; // return +1
3209 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
3210 // -2^31-1/2 < x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded
3211 // to nearest-away to a 32-bit signed integer
3212 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
3213 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
3214 // chop off ind digits from the lower part of C1
3215 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
3218 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
3220 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
3221 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
3223 if (C1
.w
[0] < tmp64
)
3225 // calculate C* and f*
3226 // C* is actually floor(C*) in this case
3227 // C* and f* need shifting and masking, as shown by
3228 // shiftright128[] and maskhigh128[]
3230 // kx = 10^(-x) = ten2mk128[ind - 1]
3231 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
3232 // the approximation of 10^(-x) was rounded up to 118 bits
3233 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
3234 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
3235 Cstar
.w
[1] = P256
.w
[3];
3236 Cstar
.w
[0] = P256
.w
[2];
3237 } else { // 22 <= ind - 1 <= 33
3239 Cstar
.w
[0] = P256
.w
[3];
3241 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
3242 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
3243 // if (0 < f* < 10^(-x)) then the result is a midpoint
3244 // if floor(C*) is even then C* = floor(C*) - logical right
3245 // shift; C* has p decimal digits, correct by Prop. 1)
3246 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
3247 // shift; C* has p decimal digits, correct by Pr. 1)
3249 // C* = floor(C*) (logical right shift; C has p decimal digits,
3250 // correct by Property 1)
3251 // n = C* * 10^(e+x)
3253 // shift right C* by Ex-128 = shiftright128[ind]
3254 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
3255 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
3257 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
3258 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
3259 } else { // 22 <= ind - 1 <= 33
3260 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
3262 // if the result was a midpoint, it was already rounded away from zero
3267 // no need to check for midpoints - already rounded away from zero!
3268 } else if (exp
== 0) {
3270 // res = +/-C (exact)
3275 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
3276 // res = +/-C * 10^exp (exact)
3278 res
= -C1
.w
[0] * ten2k64
[exp
];
3280 res
= C1
.w
[0] * ten2k64
[exp
];
3288 /*****************************************************************************
3289 * BID128_to_int32_xrninta
3290 ****************************************************************************/
3292 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_xrninta
,
3298 int exp
; // unbiased exponent
3299 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
3300 UINT64 tmp64
, tmp64A
;
3301 BID_UI64DOUBLE tmp1
;
3302 unsigned int x_nr_bits
;
3305 UINT128 Cstar
; // C* represents up to 34 decimal digits ~ 113 bits
3310 x_sign
= x
.w
[1] & MASK_SIGN
; // 0 for positive, MASK_SIGN for negative
3311 x_exp
= x
.w
[1] & MASK_EXP
; // biased and shifted left 49 bit positions
3312 C1
.w
[1] = x
.w
[1] & MASK_COEFF
;
3315 // check for NaN or Infinity
3316 if ((x
.w
[1] & MASK_SPECIAL
) == MASK_SPECIAL
) {
3318 if ((x
.w
[1] & MASK_NAN
) == MASK_NAN
) { // x is NAN
3319 if ((x
.w
[1] & MASK_SNAN
) == MASK_SNAN
) { // x is SNAN
3321 *pfpsf
|= INVALID_EXCEPTION
;
3322 // return Integer Indefinite
3324 } else { // x is QNaN
3326 *pfpsf
|= INVALID_EXCEPTION
;
3327 // return Integer Indefinite
3331 } else { // x is not a NaN, so it must be infinity
3332 if (!x_sign
) { // x is +inf
3334 *pfpsf
|= INVALID_EXCEPTION
;
3335 // return Integer Indefinite
3337 } else { // x is -inf
3339 *pfpsf
|= INVALID_EXCEPTION
;
3340 // return Integer Indefinite
3346 // check for non-canonical values (after the check for special values)
3347 if ((C1
.w
[1] > 0x0001ed09bead87c0ull
)
3348 || (C1
.w
[1] == 0x0001ed09bead87c0ull
3349 && (C1
.w
[0] > 0x378d8e63ffffffffull
))
3350 || ((x
.w
[1] & 0x6000000000000000ull
) == 0x6000000000000000ull
)) {
3353 } else if ((C1
.w
[1] == 0x0ull
) && (C1
.w
[0] == 0x0ull
)) {
3357 } else { // x is not special and is not zero
3359 // q = nr. of decimal digits in x
3360 // determine first the nr. of bits in x
3362 if (C1
.w
[0] >= 0x0020000000000000ull
) { // x >= 2^53
3363 // split the 64-bit value in two 32-bit halves to avoid rounding errors
3364 if (C1
.w
[0] >= 0x0000000100000000ull
) { // x >= 2^32
3365 tmp1
.d
= (double) (C1
.w
[0] >> 32); // exact conversion
3367 33 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3368 } else { // x < 2^32
3369 tmp1
.d
= (double) (C1
.w
[0]); // exact conversion
3371 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3373 } else { // if x < 2^53
3374 tmp1
.d
= (double) C1
.w
[0]; // exact conversion
3376 1 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3378 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
3379 tmp1
.d
= (double) C1
.w
[1]; // exact conversion
3381 65 + ((((unsigned int) (tmp1
.ui64
>> 52)) & 0x7ff) - 0x3ff);
3383 q
= nr_digits
[x_nr_bits
- 1].digits
;
3385 q
= nr_digits
[x_nr_bits
- 1].digits1
;
3386 if (C1
.w
[1] > nr_digits
[x_nr_bits
- 1].threshold_hi
3387 || (C1
.w
[1] == nr_digits
[x_nr_bits
- 1].threshold_hi
3388 && C1
.w
[0] >= nr_digits
[x_nr_bits
- 1].threshold_lo
))
3391 exp
= (x_exp
>> 49) - 6176;
3392 if ((q
+ exp
) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
3394 *pfpsf
|= INVALID_EXCEPTION
;
3395 // return Integer Indefinite
3398 } else if ((q
+ exp
) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
3399 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
3400 // so x rounded to an integer may or may not fit in a signed 32-bit int
3401 // the cases that do not fit are identified here; the ones that fit
3402 // fall through and will be handled with other cases further,
3403 // under '1 <= q + exp <= 10'
3404 if (x_sign
) { // if n < 0 and q + exp = 10
3405 // if n <= -2^31 - 1/2 then n is too large
3406 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1/2
3407 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000005, 1<=q<=34
3409 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
3410 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3411 if (tmp64
>= 0x500000005ull
) {
3413 *pfpsf
|= INVALID_EXCEPTION
;
3414 // return Integer Indefinite
3418 // else cases that can be rounded to a 32-bit int fall through
3419 // to '1 <= q + exp <= 10'
3420 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3421 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000005 <=>
3422 // C >= 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 23
3423 // (scale 2^31+1/2 up)
3424 tmp64
= 0x500000005ull
;
3425 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3426 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
3427 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3428 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
3430 if (C1
.w
[1] > C
.w
[1]
3431 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
3433 *pfpsf
|= INVALID_EXCEPTION
;
3434 // return Integer Indefinite
3438 // else cases that can be rounded to a 32-bit int fall through
3439 // to '1 <= q + exp <= 10'
3441 } else { // if n > 0 and q + exp = 10
3442 // if n >= 2^31 - 1/2 then n is too large
3443 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2
3444 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=34
3446 tmp64
= C1
.w
[0] * ten2k64
[11 - q
]; // C scaled up to 11-digit int
3447 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3448 if (tmp64
>= 0x4fffffffbull
) {
3450 *pfpsf
|= INVALID_EXCEPTION
;
3451 // return Integer Indefinite
3455 // else cases that can be rounded to a 32-bit int fall through
3456 // to '1 <= q + exp <= 10'
3457 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3458 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb <=>
3459 // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
3460 // (scale 2^31-1/2 up)
3461 tmp64
= 0x4fffffffbull
;
3462 if (q
- 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3463 __mul_64x64_to_128MACH (C
, tmp64
, ten2k64
[q
- 11]);
3464 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3465 __mul_128x64_to_128 (C
, tmp64
, ten2k128
[q
- 31]);
3467 if (C1
.w
[1] > C
.w
[1]
3468 || (C1
.w
[1] == C
.w
[1] && C1
.w
[0] >= C
.w
[0])) {
3470 *pfpsf
|= INVALID_EXCEPTION
;
3471 // return Integer Indefinite
3475 // else cases that can be rounded to a 32-bit int fall through
3476 // to '1 <= q + exp <= 10'
3480 // n is not too large to be converted to int32: -2^31 - 1/2 < n < 2^31 - 1/2
3481 // Note: some of the cases tested for above fall through to this point
3482 if ((q
+ exp
) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
3484 *pfpsf
|= INEXACT_EXCEPTION
;
3488 } else if ((q
+ exp
) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
3489 // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
3494 if (ind
<= 18) { // 0 <= ind <= 18
3495 if ((C1
.w
[1] == 0) && (C1
.w
[0] < midpoint64
[ind
])) {
3496 res
= 0x00000000; // return 0
3497 } else if (x_sign
) { // n < 0
3498 res
= 0xffffffff; // return -1
3500 res
= 0x00000001; // return +1
3502 } else { // 19 <= ind <= 33
3503 if ((C1
.w
[1] < midpoint128
[ind
- 19].w
[1])
3504 || ((C1
.w
[1] == midpoint128
[ind
- 19].w
[1])
3505 && (C1
.w
[0] < midpoint128
[ind
- 19].w
[0]))) {
3506 res
= 0x00000000; // return 0
3507 } else if (x_sign
) { // n < 0
3508 res
= 0xffffffff; // return -1
3510 res
= 0x00000001; // return +1
3514 *pfpsf
|= INEXACT_EXCEPTION
;
3515 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
3516 // -2^31-1/2 < x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded
3517 // to nearest-away to a 32-bit signed integer
3518 if (exp
< 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
3519 ind
= -exp
; // 1 <= ind <= 33; ind is a synonym for 'x'
3520 // chop off ind digits from the lower part of C1
3521 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
3524 C1
.w
[0] = C1
.w
[0] + midpoint64
[ind
- 1];
3526 C1
.w
[0] = C1
.w
[0] + midpoint128
[ind
- 20].w
[0];
3527 C1
.w
[1] = C1
.w
[1] + midpoint128
[ind
- 20].w
[1];
3529 if (C1
.w
[0] < tmp64
)
3531 // calculate C* and f*
3532 // C* is actually floor(C*) in this case
3533 // C* and f* need shifting and masking, as shown by
3534 // shiftright128[] and maskhigh128[]
3536 // kx = 10^(-x) = ten2mk128[ind - 1]
3537 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
3538 // the approximation of 10^(-x) was rounded up to 118 bits
3539 __mul_128x128_to_256 (P256
, C1
, ten2mk128
[ind
- 1]);
3540 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
3541 Cstar
.w
[1] = P256
.w
[3];
3542 Cstar
.w
[0] = P256
.w
[2];
3544 fstar
.w
[2] = P256
.w
[2] & maskhigh128
[ind
- 1];
3545 fstar
.w
[1] = P256
.w
[1];
3546 fstar
.w
[0] = P256
.w
[0];
3547 } else { // 22 <= ind - 1 <= 33
3549 Cstar
.w
[0] = P256
.w
[3];
3550 fstar
.w
[3] = P256
.w
[3] & maskhigh128
[ind
- 1];
3551 fstar
.w
[2] = P256
.w
[2];
3552 fstar
.w
[1] = P256
.w
[1];
3553 fstar
.w
[0] = P256
.w
[0];
3555 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
3556 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
3557 // if (0 < f* < 10^(-x)) then the result is a midpoint
3558 // if floor(C*) is even then C* = floor(C*) - logical right
3559 // shift; C* has p decimal digits, correct by Prop. 1)
3560 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
3561 // shift; C* has p decimal digits, correct by Pr. 1)
3563 // C* = floor(C*) (logical right shift; C has p decimal digits,
3564 // correct by Property 1)
3565 // n = C* * 10^(e+x)
3567 // shift right C* by Ex-128 = shiftright128[ind]
3568 shift
= shiftright128
[ind
- 1]; // 0 <= shift <= 102
3569 if (ind
- 1 <= 21) { // 0 <= ind - 1 <= 21
3571 (Cstar
.w
[0] >> shift
) | (Cstar
.w
[1] << (64 - shift
));
3572 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
3573 } else { // 22 <= ind - 1 <= 33
3574 Cstar
.w
[0] = (Cstar
.w
[0] >> (shift
- 64)); // 2 <= shift - 64 <= 38
3576 // if the result was a midpoint, it was already rounded away from zero
3581 // determine inexactness of the rounding of C*
3582 // if (0 < f* - 1/2 < 10^(-x)) then
3583 // the result is exact
3584 // else // if (f* - 1/2 > T*) then
3585 // the result is inexact
3587 if (fstar
.w
[1] > 0x8000000000000000ull
|| (fstar
.w
[1] == 0x8000000000000000ull
&& fstar
.w
[0] > 0x0ull
)) { // f* > 1/2 and the result may be exact
3588 tmp64
= fstar
.w
[1] - 0x8000000000000000ull
; // f* - 1/2
3589 if ((tmp64
> ten2mk128trunc
[ind
- 1].w
[1]
3590 || (tmp64
== ten2mk128trunc
[ind
- 1].w
[1]
3591 && fstar
.w
[0] >= ten2mk128trunc
[ind
- 1].w
[0]))) {
3592 // set the inexact flag
3593 *pfpsf
|= INEXACT_EXCEPTION
;
3594 } // else the result is exact
3595 } else { // the result is inexact; f2* <= 1/2
3596 // set the inexact flag
3597 *pfpsf
|= INEXACT_EXCEPTION
;
3599 } else if (ind
- 1 <= 21) { // if 3 <= ind <= 21
3600 if (fstar
.w
[3] > 0x0 ||
3601 (fstar
.w
[3] == 0x0 && fstar
.w
[2] > onehalf128
[ind
- 1]) ||
3602 (fstar
.w
[3] == 0x0 && fstar
.w
[2] == onehalf128
[ind
- 1] &&
3603 (fstar
.w
[1] || fstar
.w
[0]))) {
3604 // f2* > 1/2 and the result may be exact
3605 // Calculate f2* - 1/2
3606 tmp64
= fstar
.w
[2] - onehalf128
[ind
- 1];
3607 tmp64A
= fstar
.w
[3];
3608 if (tmp64
> fstar
.w
[2])
3611 || fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
3612 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
3613 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
3614 // set the inexact flag
3615 *pfpsf
|= INEXACT_EXCEPTION
;
3616 } // else the result is exact
3617 } else { // the result is inexact; f2* <= 1/2
3618 // set the inexact flag
3619 *pfpsf
|= INEXACT_EXCEPTION
;
3621 } else { // if 22 <= ind <= 33
3622 if (fstar
.w
[3] > onehalf128
[ind
- 1] ||
3623 (fstar
.w
[3] == onehalf128
[ind
- 1] &&
3624 (fstar
.w
[2] || fstar
.w
[1] || fstar
.w
[0]))) {
3625 // f2* > 1/2 and the result may be exact
3626 // Calculate f2* - 1/2
3627 tmp64
= fstar
.w
[3] - onehalf128
[ind
- 1];
3628 if (tmp64
|| fstar
.w
[2] ||
3629 fstar
.w
[1] > ten2mk128trunc
[ind
- 1].w
[1]
3630 || (fstar
.w
[1] == ten2mk128trunc
[ind
- 1].w
[1]
3631 && fstar
.w
[0] > ten2mk128trunc
[ind
- 1].w
[0])) {
3632 // set the inexact flag
3633 *pfpsf
|= INEXACT_EXCEPTION
;
3634 } // else the result is exact
3635 } else { // the result is inexact; f2* <= 1/2
3636 // set the inexact flag
3637 *pfpsf
|= INEXACT_EXCEPTION
;
3640 // no need to check for midpoints - already rounded away from zero!
3641 } else if (exp
== 0) {
3643 // res = +/-C (exact)
3648 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
3649 // res = +/-C * 10^exp (exact)
3651 res
= -C1
.w
[0] * ten2k64
[exp
];
3653 res
= C1
.w
[0] * ten2k64
[exp
];