Mercurial > hg > CbC > CbC_gcc

comparison libgcc/config/libbid/bid128_to_int32.c @ 0:a06113de4d67

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first commit
author kent <kent@cr.ie.u-ryukyu.ac.jp>
date Fri, 17 Jul 2009 14:47:48 +0900
parents
children 04ced10e8804
comparison
equal deleted inserted replaced
-1:000000000000 0:a06113de4d67
1 /* Copyright (C) 2007, 2009 Free Software Foundation, Inc.
2
3 This file is part of GCC.
4
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
8 version.
9
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
13 for more details.
14
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.
18
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/>. */
23
24 #include "bid_internal.h"
25
26 /*****************************************************************************
27 * BID128_to_int32_rnint
28 ****************************************************************************/
29
30 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_rnint, x)
31
32 int res;
33 UINT64 x_sign;
34 UINT64 x_exp;
35 int exp; // unbiased exponent
36 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
37 UINT64 tmp64;
38 BID_UI64DOUBLE tmp1;
39 unsigned int x_nr_bits;
40 int q, ind, shift;
41 UINT128 C1, C;
42 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
43 UINT256 fstar;
44 UINT256 P256;
45
46 // unpack x
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;
50 C1.w[0] = x.w[0];
51
52 // check for NaN or Infinity
53 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
54 // x is 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
57 // set invalid flag
58 *pfpsf |= INVALID_EXCEPTION;
59 // return Integer Indefinite
60 res = 0x80000000;
61 } else { // x is QNaN
62 // set invalid flag
63 *pfpsf |= INVALID_EXCEPTION;
64 // return Integer Indefinite
65 res = 0x80000000;
66 }
67 BID_RETURN (res);
68 } else { // x is not a NaN, so it must be infinity
69 if (!x_sign) { // x is +inf
70 // set invalid flag
71 *pfpsf |= INVALID_EXCEPTION;
72 // return Integer Indefinite
73 res = 0x80000000;
74 } else { // x is -inf
75 // set invalid flag
76 *pfpsf |= INVALID_EXCEPTION;
77 // return Integer Indefinite
78 res = 0x80000000;
79 }
80 BID_RETURN (res);
81 }
82 }
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)) {
88 res = 0x00000000;
89 BID_RETURN (res);
90 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
91 // x is 0
92 res = 0x00000000;
93 BID_RETURN (res);
94 } else { // x is not special and is not zero
95
96 // q = nr. of decimal digits in x
97 // determine first the nr. of bits in x
98 if (C1.w[1] == 0) {
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
103 x_nr_bits =
104 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
105 } else { // x < 2^32
106 tmp1.d = (double) (C1.w[0]); // exact conversion
107 x_nr_bits =
108 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
109 }
110 } else { // if x < 2^53
111 tmp1.d = (double) C1.w[0]; // exact conversion
112 x_nr_bits =
113 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
114 }
115 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
116 tmp1.d = (double) C1.w[1]; // exact conversion
117 x_nr_bits =
118 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
119 }
120 q = nr_digits[x_nr_bits - 1].digits;
121 if (q == 0) {
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))
126 q++;
127 }
128 exp = (x_exp >> 49) - 6176;
129 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
130 // set invalid flag
131 *pfpsf |= INVALID_EXCEPTION;
132 // return Integer Indefinite
133 res = 0x80000000;
134 BID_RETURN (res);
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
145 if (q <= 11) {
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) {
149 // set invalid flag
150 *pfpsf |= INVALID_EXCEPTION;
151 // return Integer Indefinite
152 res = 0x80000000;
153 BID_RETURN (res);
154 }
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]);
166 }
167 if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
168 // set invalid flag
169 *pfpsf |= INVALID_EXCEPTION;
170 // return Integer Indefinite
171 res = 0x80000000;
172 BID_RETURN (res);
173 }
174 // else cases that can be rounded to a 32-bit int fall through
175 // to '1 <= q + exp <= 10'
176 }
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
181 if (q <= 11) {
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) {
185 // set invalid flag
186 *pfpsf |= INVALID_EXCEPTION;
187 // return Integer Indefinite
188 res = 0x80000000;
189 BID_RETURN (res);
190 }
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]);
202 }
203 if (C1.w[1] > C.w[1]
204 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
205 // set invalid flag
206 *pfpsf |= INVALID_EXCEPTION;
207 // return Integer Indefinite
208 res = 0x80000000;
209 BID_RETURN (res);
210 }
211 // else cases that can be rounded to a 32-bit int fall through
212 // to '1 <= q + exp <= 10'
213 }
214 }
215 }
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)
219 // return 0
220 res = 0x00000000;
221 BID_RETURN (res);
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)
224 // res = 0
225 // else
226 // res = +/-1
227 ind = 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
233 } else { // n > 0
234 res = 0x00000001; // return +1
235 }
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
243 } else { // n > 0
244 res = 0x00000001; // return +1
245 }
246 }
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
254 tmp64 = C1.w[0];
255 if (ind <= 19) {
256 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
257 } else {
258 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
259 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
260 }
261 if (C1.w[0] < tmp64)
262 C1.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[]
267 // 1 <= x <= 33
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];
275 fstar.w[3] = 0;
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
280 Cstar.w[1] = 0;
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];
286 }
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)
294 // else
295 // C* = floor(C*) (logical right shift; C has p decimal digits,
296 // correct by Property 1)
297 // n = C* * 10^(e+x)
298
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
302 Cstar.w[0] =
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
307 }
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]
321 }
322 if (x_sign)
323 res = -Cstar.w[0];
324 else
325 res = Cstar.w[0];
326 } else if (exp == 0) {
327 // 1 <= q <= 10
328 // res = +/-C (exact)
329 if (x_sign)
330 res = -C1.w[0];
331 else
332 res = C1.w[0];
333 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
334 // res = +/-C * 10^exp (exact)
335 if (x_sign)
336 res = -C1.w[0] * ten2k64[exp];
337 else
338 res = C1.w[0] * ten2k64[exp];
339 }
340 }
341 }
342
343 BID_RETURN (res);
344 }
345
346 /*****************************************************************************
347 * BID128_to_int32_xrnint
348 ****************************************************************************/
349
350 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_xrnint,
351 x)
352
353 int res;
354 UINT64 x_sign;
355 UINT64 x_exp;
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;
359 BID_UI64DOUBLE tmp1;
360 unsigned int x_nr_bits;
361 int q, ind, shift;
362 UINT128 C1, C;
363 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
364 UINT256 fstar;
365 UINT256 P256;
366
367 // unpack x
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;
371 C1.w[0] = x.w[0];
372
373 // check for NaN or Infinity
374 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
375 // x is 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
378 // set invalid flag
379 *pfpsf |= INVALID_EXCEPTION;
380 // return Integer Indefinite
381 res = 0x80000000;
382 } else { // x is QNaN
383 // set invalid flag
384 *pfpsf |= INVALID_EXCEPTION;
385 // return Integer Indefinite
386 res = 0x80000000;
387 }
388 BID_RETURN (res);
389 } else { // x is not a NaN, so it must be infinity
390 if (!x_sign) { // x is +inf
391 // set invalid flag
392 *pfpsf |= INVALID_EXCEPTION;
393 // return Integer Indefinite
394 res = 0x80000000;
395 } else { // x is -inf
396 // set invalid flag
397 *pfpsf |= INVALID_EXCEPTION;
398 // return Integer Indefinite
399 res = 0x80000000;
400 }
401 BID_RETURN (res);
402 }
403 }
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)) {
409 res = 0x00000000;
410 BID_RETURN (res);
411 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
412 // x is 0
413 res = 0x00000000;
414 BID_RETURN (res);
415 } else { // x is not special and is not zero
416
417 // q = nr. of decimal digits in x
418 // determine first the nr. of bits in x
419 if (C1.w[1] == 0) {
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
424 x_nr_bits =
425 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
426 } else { // x < 2^32
427 tmp1.d = (double) (C1.w[0]); // exact conversion
428 x_nr_bits =
429 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
430 }
431 } else { // if x < 2^53
432 tmp1.d = (double) C1.w[0]; // exact conversion
433 x_nr_bits =
434 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
435 }
436 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
437 tmp1.d = (double) C1.w[1]; // exact conversion
438 x_nr_bits =
439 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
440 }
441 q = nr_digits[x_nr_bits - 1].digits;
442 if (q == 0) {
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))
447 q++;
448 }
449 exp = (x_exp >> 49) - 6176;
450 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
451 // set invalid flag
452 *pfpsf |= INVALID_EXCEPTION;
453 // return Integer Indefinite
454 res = 0x80000000;
455 BID_RETURN (res);
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
466 if (q <= 11) {
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) {
470 // set invalid flag
471 *pfpsf |= INVALID_EXCEPTION;
472 // return Integer Indefinite
473 res = 0x80000000;
474 BID_RETURN (res);
475 }
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]);
487 }
488 if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
489 // set invalid flag
490 *pfpsf |= INVALID_EXCEPTION;
491 // return Integer Indefinite
492 res = 0x80000000;
493 BID_RETURN (res);
494 }
495 // else cases that can be rounded to a 32-bit int fall through
496 // to '1 <= q + exp <= 10'
497 }
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
502 if (q <= 11) {
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) {
506 // set invalid flag
507 *pfpsf |= INVALID_EXCEPTION;
508 // return Integer Indefinite
509 res = 0x80000000;
510 BID_RETURN (res);
511 }
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]);
523 }
524 if (C1.w[1] > C.w[1]
525 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
526 // set invalid flag
527 *pfpsf |= INVALID_EXCEPTION;
528 // return Integer Indefinite
529 res = 0x80000000;
530 BID_RETURN (res);
531 }
532 // else cases that can be rounded to a 32-bit int fall through
533 // to '1 <= q + exp <= 10'
534 }
535 }
536 }
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)
540 // set inexact flag
541 *pfpsf |= INEXACT_EXCEPTION;
542 // return 0
543 res = 0x00000000;
544 BID_RETURN (res);
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)
547 // res = 0
548 // else
549 // res = +/-1
550 ind = 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
556 } else { // n > 0
557 res = 0x00000001; // return +1
558 }
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
566 } else { // n > 0
567 res = 0x00000001; // return +1
568 }
569 }
570 // set inexact flag
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
579 tmp64 = C1.w[0];
580 if (ind <= 19) {
581 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
582 } else {
583 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
584 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
585 }
586 if (C1.w[0] < tmp64)
587 C1.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[]
592 // 1 <= x <= 33
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];
600 fstar.w[3] = 0;
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
605 Cstar.w[1] = 0;
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];
611 }
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)
619 // else
620 // C* = floor(C*) (logical right shift; C has p decimal digits,
621 // correct by Property 1)
622 // n = C* * 10^(e+x)
623
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
627 Cstar.w[0] =
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
632 }
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
638 if (ind - 1 <= 2) {
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;
650 }
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];
659 tmp64A = fstar.w[3];
660 if (tmp64 > fstar.w[2])
661 tmp64A--;
662 if (tmp64A || tmp64
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;
672 }
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;
690 }
691 }
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]
705 }
706 if (x_sign)
707 res = -Cstar.w[0];
708 else
709 res = Cstar.w[0];
710 } else if (exp == 0) {
711 // 1 <= q <= 10
712 // res = +/-C (exact)
713 if (x_sign)
714 res = -C1.w[0];
715 else
716 res = C1.w[0];
717 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
718 // res = +/-C * 10^exp (exact)
719 if (x_sign)
720 res = -C1.w[0] * ten2k64[exp];
721 else
722 res = C1.w[0] * ten2k64[exp];
723 }
724 }
725 }
726
727 BID_RETURN (res);
728 }
729
730 /*****************************************************************************
731 * BID128_to_int32_floor
732 ****************************************************************************/
733
734 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_floor, x)
735
736 int res;
737 UINT64 x_sign;
738 UINT64 x_exp;
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;
742 BID_UI64DOUBLE tmp1;
743 unsigned int x_nr_bits;
744 int q, ind, shift;
745 UINT128 C1, C;
746 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
747 UINT256 fstar;
748 UINT256 P256;
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;
753
754 // unpack x
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;
758 C1.w[0] = x.w[0];
759
760 // check for NaN or Infinity
761 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
762 // x is 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
765 // set invalid flag
766 *pfpsf |= INVALID_EXCEPTION;
767 // return Integer Indefinite
768 res = 0x80000000;
769 } else { // x is QNaN
770 // set invalid flag
771 *pfpsf |= INVALID_EXCEPTION;
772 // return Integer Indefinite
773 res = 0x80000000;
774 }
775 BID_RETURN (res);
776 } else { // x is not a NaN, so it must be infinity
777 if (!x_sign) { // x is +inf
778 // set invalid flag
779 *pfpsf |= INVALID_EXCEPTION;
780 // return Integer Indefinite
781 res = 0x80000000;
782 } else { // x is -inf
783 // set invalid flag
784 *pfpsf |= INVALID_EXCEPTION;
785 // return Integer Indefinite
786 res = 0x80000000;
787 }
788 BID_RETURN (res);
789 }
790 }
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)) {
796 res = 0x00000000;
797 BID_RETURN (res);
798 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
799 // x is 0
800 res = 0x00000000;
801 BID_RETURN (res);
802 } else { // x is not special and is not zero
803
804 // q = nr. of decimal digits in x
805 // determine first the nr. of bits in x
806 if (C1.w[1] == 0) {
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
811 x_nr_bits =
812 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
813 } else { // x < 2^32
814 tmp1.d = (double) (C1.w[0]); // exact conversion
815 x_nr_bits =
816 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
817 }
818 } else { // if x < 2^53
819 tmp1.d = (double) C1.w[0]; // exact conversion
820 x_nr_bits =
821 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
822 }
823 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
824 tmp1.d = (double) C1.w[1]; // exact conversion
825 x_nr_bits =
826 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
827 }
828 q = nr_digits[x_nr_bits - 1].digits;
829 if (q == 0) {
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))
834 q++;
835 }
836 exp = (x_exp >> 49) - 6176;
837 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
838 // set invalid flag
839 *pfpsf |= INVALID_EXCEPTION;
840 // return Integer Indefinite
841 res = 0x80000000;
842 BID_RETURN (res);
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
853 if (q <= 11) {
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) {
857 // set invalid flag
858 *pfpsf |= INVALID_EXCEPTION;
859 // return Integer Indefinite
860 res = 0x80000000;
861 BID_RETURN (res);
862 }
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
868 // (scale 2^31 up)
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]);
874 }
875 if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
876 // set invalid flag
877 *pfpsf |= INVALID_EXCEPTION;
878 // return Integer Indefinite
879 res = 0x80000000;
880 BID_RETURN (res);
881 }
882 // else cases that can be rounded to a 32-bit int fall through
883 // to '1 <= q + exp <= 10'
884 }
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
889 if (q <= 11) {
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) {
893 // set invalid flag
894 *pfpsf |= INVALID_EXCEPTION;
895 // return Integer Indefinite
896 res = 0x80000000;
897 BID_RETURN (res);
898 }
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
904 // (scale 2^31 up)
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]);
910 }
911 if (C1.w[1] > C.w[1]
912 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
913 // set invalid flag
914 *pfpsf |= INVALID_EXCEPTION;
915 // return Integer Indefinite
916 res = 0x80000000;
917 BID_RETURN (res);
918 }
919 // else cases that can be rounded to a 32-bit int fall through
920 // to '1 <= q + exp <= 10'
921 }
922 }
923 }
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)
928 // return 0
929 if (x_sign)
930 res = 0xffffffff;
931 else
932 res = 0x00000000;
933 BID_RETURN (res);
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
941 tmp64 = C1.w[0];
942 if (ind <= 19) {
943 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
944 } else {
945 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
946 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
947 }
948 if (C1.w[0] < tmp64)
949 C1.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[]
954 // 1 <= x <= 33
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];
962 fstar.w[3] = 0;
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
967 Cstar.w[1] = 0;
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];
973 }
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)
981 // else
982 // C* = floor(C*) (logical right shift; C has p decimal digits,
983 // correct by Property 1)
984 // n = C* * 10^(e+x)
985
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
989 Cstar.w[0] =
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
994 }
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
1000 if (ind - 1 <= 2) {
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;
1010 }
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])
1021 tmp64A--;
1022 if (tmp64A || tmp64
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;
1030 }
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;
1046 }
1047 }
1048
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;
1068 }
1069 }
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;
1073 } else if (!x_sign
1074 && (is_midpoint_lt_even || is_inexact_gt_midpoint)) {
1075 Cstar.w[0] = Cstar.w[0] - 1;
1076 } else {
1077 ; // the result is already correct
1078 }
1079 if (x_sign)
1080 res = -Cstar.w[0];
1081 else
1082 res = Cstar.w[0];
1083 } else if (exp == 0) {
1084 // 1 <= q <= 10
1085 // res = +/-C (exact)
1086 if (x_sign)
1087 res = -C1.w[0];
1088 else
1089 res = C1.w[0];
1090 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1091 // res = +/-C * 10^exp (exact)
1092 if (x_sign)
1093 res = -C1.w[0] * ten2k64[exp];
1094 else
1095 res = C1.w[0] * ten2k64[exp];
1096 }
1097 }
1098 }
1099
1100 BID_RETURN (res);
1101 }
1102
1103
1104 /*****************************************************************************
1105 * BID128_to_int32_xfloor
1106 ****************************************************************************/
1107
1108 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_xfloor,
1109 x)
1110
1111 int res;
1112 UINT64 x_sign;
1113 UINT64 x_exp;
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;
1119 int q, ind, shift;
1120 UINT128 C1, C;
1121 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
1122 UINT256 fstar;
1123 UINT256 P256;
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;
1128
1129 // unpack x
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;
1133 C1.w[0] = x.w[0];
1134
1135 // check for NaN or Infinity
1136 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1137 // x is 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
1140 // set invalid flag
1141 *pfpsf |= INVALID_EXCEPTION;
1142 // return Integer Indefinite
1143 res = 0x80000000;
1144 } else { // x is QNaN
1145 // set invalid flag
1146 *pfpsf |= INVALID_EXCEPTION;
1147 // return Integer Indefinite
1148 res = 0x80000000;
1149 }
1150 BID_RETURN (res);
1151 } else { // x is not a NaN, so it must be infinity
1152 if (!x_sign) { // x is +inf
1153 // set invalid flag
1154 *pfpsf |= INVALID_EXCEPTION;
1155 // return Integer Indefinite
1156 res = 0x80000000;
1157 } else { // x is -inf
1158 // set invalid flag
1159 *pfpsf |= INVALID_EXCEPTION;
1160 // return Integer Indefinite
1161 res = 0x80000000;
1162 }
1163 BID_RETURN (res);
1164 }
1165 }
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)) {
1171 res = 0x00000000;
1172 BID_RETURN (res);
1173 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1174 // x is 0
1175 res = 0x00000000;
1176 BID_RETURN (res);
1177 } else { // x is not special and is not zero
1178
1179 // q = nr. of decimal digits in x
1180 // determine first the nr. of bits in x
1181 if (C1.w[1] == 0) {
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
1186 x_nr_bits =
1187 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1188 } else { // x < 2^32
1189 tmp1.d = (double) (C1.w[0]); // exact conversion
1190 x_nr_bits =
1191 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1192 }
1193 } else { // if x < 2^53
1194 tmp1.d = (double) C1.w[0]; // exact conversion
1195 x_nr_bits =
1196 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1197 }
1198 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1199 tmp1.d = (double) C1.w[1]; // exact conversion
1200 x_nr_bits =
1201 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1202 }
1203 q = nr_digits[x_nr_bits - 1].digits;
1204 if (q == 0) {
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))
1209 q++;
1210 }
1211 exp = (x_exp >> 49) - 6176;
1212 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1213 // set invalid flag
1214 *pfpsf |= INVALID_EXCEPTION;
1215 // return Integer Indefinite
1216 res = 0x80000000;
1217 BID_RETURN (res);
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
1228 if (q <= 11) {
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) {
1232 // set invalid flag
1233 *pfpsf |= INVALID_EXCEPTION;
1234 // return Integer Indefinite
1235 res = 0x80000000;
1236 BID_RETURN (res);
1237 }
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
1243 // (scale 2^31 up)
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]);
1249 }
1250 if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
1251 // set invalid flag
1252 *pfpsf |= INVALID_EXCEPTION;
1253 // return Integer Indefinite
1254 res = 0x80000000;
1255 BID_RETURN (res);
1256 }
1257 // else cases that can be rounded to a 32-bit int fall through
1258 // to '1 <= q + exp <= 10'
1259 }
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
1264 if (q <= 11) {
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) {
1268 // set invalid flag
1269 *pfpsf |= INVALID_EXCEPTION;
1270 // return Integer Indefinite
1271 res = 0x80000000;
1272 BID_RETURN (res);
1273 }
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
1279 // (scale 2^31 up)
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]);
1285 }
1286 if (C1.w[1] > C.w[1]
1287 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
1288 // set invalid flag
1289 *pfpsf |= INVALID_EXCEPTION;
1290 // return Integer Indefinite
1291 res = 0x80000000;
1292 BID_RETURN (res);
1293 }
1294 // else cases that can be rounded to a 32-bit int fall through
1295 // to '1 <= q + exp <= 10'
1296 }
1297 }
1298 }
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)
1303 // set inexact flag
1304 *pfpsf |= INEXACT_EXCEPTION;
1305 // return 0
1306 if (x_sign)
1307 res = 0xffffffff;
1308 else
1309 res = 0x00000000;
1310 BID_RETURN (res);
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
1318 tmp64 = C1.w[0];
1319 if (ind <= 19) {
1320 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1321 } else {
1322 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1323 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1324 }
1325 if (C1.w[0] < tmp64)
1326 C1.w[1]++;
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[]
1331 // 1 <= x <= 33
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];
1339 fstar.w[3] = 0;
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
1344 Cstar.w[1] = 0;
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];
1350 }
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)
1358 // else
1359 // C* = floor(C*) (logical right shift; C has p decimal digits,
1360 // correct by Property 1)
1361 // n = C* * 10^(e+x)
1362
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
1366 Cstar.w[0] =
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
1371 }
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
1377 if (ind - 1 <= 2) {
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;
1391 }
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])
1402 tmp64A--;
1403 if (tmp64A || tmp64
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;
1415 }
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;
1435 }
1436 }
1437
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;
1457 }
1458 }
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;
1462 } else if (!x_sign
1463 && (is_midpoint_lt_even || is_inexact_gt_midpoint)) {
1464 Cstar.w[0] = Cstar.w[0] - 1;
1465 } else {
1466 ; // the result is already correct
1467 }
1468 if (x_sign)
1469 res = -Cstar.w[0];
1470 else
1471 res = Cstar.w[0];
1472 } else if (exp == 0) {
1473 // 1 <= q <= 10
1474 // res = +/-C (exact)
1475 if (x_sign)
1476 res = -C1.w[0];
1477 else
1478 res = C1.w[0];
1479 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1480 // res = +/-C * 10^exp (exact)
1481 if (x_sign)
1482 res = -C1.w[0] * ten2k64[exp];
1483 else
1484 res = C1.w[0] * ten2k64[exp];
1485 }
1486 }
1487 }
1488
1489 BID_RETURN (res);
1490 }
1491
1492 /*****************************************************************************
1493 * BID128_to_int32_ceil
1494 ****************************************************************************/
1495
1496 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_ceil, x)
1497
1498 int res;
1499 UINT64 x_sign;
1500 UINT64 x_exp;
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;
1506 int q, ind, shift;
1507 UINT128 C1, C;
1508 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
1509 UINT256 fstar;
1510 UINT256 P256;
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;
1515
1516 // unpack x
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;
1520 C1.w[0] = x.w[0];
1521
1522 // check for NaN or Infinity
1523 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1524 // x is 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
1527 // set invalid flag
1528 *pfpsf |= INVALID_EXCEPTION;
1529 // return Integer Indefinite
1530 res = 0x80000000;
1531 } else { // x is QNaN
1532 // set invalid flag
1533 *pfpsf |= INVALID_EXCEPTION;
1534 // return Integer Indefinite
1535 res = 0x80000000;
1536 }
1537 BID_RETURN (res);
1538 } else { // x is not a NaN, so it must be infinity
1539 if (!x_sign) { // x is +inf
1540 // set invalid flag
1541 *pfpsf |= INVALID_EXCEPTION;
1542 // return Integer Indefinite
1543 res = 0x80000000;
1544 } else { // x is -inf
1545 // set invalid flag
1546 *pfpsf |= INVALID_EXCEPTION;
1547 // return Integer Indefinite
1548 res = 0x80000000;
1549 }
1550 BID_RETURN (res);
1551 }
1552 }
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)) {
1558 res = 0x00000000;
1559 BID_RETURN (res);
1560 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1561 // x is 0
1562 res = 0x00000000;
1563 BID_RETURN (res);
1564 } else { // x is not special and is not zero
1565
1566 // q = nr. of decimal digits in x
1567 // determine first the nr. of bits in x
1568 if (C1.w[1] == 0) {
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
1573 x_nr_bits =
1574 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1575 } else { // x < 2^32
1576 tmp1.d = (double) (C1.w[0]); // exact conversion
1577 x_nr_bits =
1578 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1579 }
1580 } else { // if x < 2^53
1581 tmp1.d = (double) C1.w[0]; // exact conversion
1582 x_nr_bits =
1583 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1584 }
1585 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1586 tmp1.d = (double) C1.w[1]; // exact conversion
1587 x_nr_bits =
1588 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1589 }
1590 q = nr_digits[x_nr_bits - 1].digits;
1591 if (q == 0) {
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))
1596 q++;
1597 }
1598 exp = (x_exp >> 49) - 6176;
1599 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1600 // set invalid flag
1601 *pfpsf |= INVALID_EXCEPTION;
1602 // return Integer Indefinite
1603 res = 0x80000000;
1604 BID_RETURN (res);
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
1615 if (q <= 11) {
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) {
1619 // set invalid flag
1620 *pfpsf |= INVALID_EXCEPTION;
1621 // return Integer Indefinite
1622 res = 0x80000000;
1623 BID_RETURN (res);
1624 }
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]);
1636 }
1637 if (C1.w[1] > C.w[1]
1638 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
1639 // set invalid flag
1640 *pfpsf |= INVALID_EXCEPTION;
1641 // return Integer Indefinite
1642 res = 0x80000000;
1643 BID_RETURN (res);
1644 }
1645 // else cases that can be rounded to a 32-bit int fall through
1646 // to '1 <= q + exp <= 10'
1647 }
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
1652 if (q <= 11) {
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) {
1656 // set invalid flag
1657 *pfpsf |= INVALID_EXCEPTION;
1658 // return Integer Indefinite
1659 res = 0x80000000;
1660 BID_RETURN (res);
1661 }
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
1667 // (scale 2^31 up)
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]);
1673 }
1674 if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
1675 // set invalid flag
1676 *pfpsf |= INVALID_EXCEPTION;
1677 // return Integer Indefinite
1678 res = 0x80000000;
1679 BID_RETURN (res);
1680 }
1681 // else cases that can be rounded to a 32-bit int fall through
1682 // to '1 <= q + exp <= 10'
1683 }
1684 }
1685 }
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)
1690 // return 0
1691 if (x_sign)
1692 res = 0x00000000;
1693 else
1694 res = 0x00000001;
1695 BID_RETURN (res);
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
1703 tmp64 = C1.w[0];
1704 if (ind <= 19) {
1705 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1706 } else {
1707 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1708 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1709 }
1710 if (C1.w[0] < tmp64)
1711 C1.w[1]++;
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[]
1716 // 1 <= x <= 33
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];
1724 fstar.w[3] = 0;
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
1729 Cstar.w[1] = 0;
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];
1735 }
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)
1743 // else
1744 // C* = floor(C*) (logical right shift; C has p decimal digits,
1745 // correct by Property 1)
1746 // n = C* * 10^(e+x)
1747
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
1751 Cstar.w[0] =
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
1756 }
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
1762 if (ind - 1 <= 2) {
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;
1772 }
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])
1783 tmp64A--;
1784 if (tmp64A || tmp64
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;
1792 }
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;
1808 }
1809 }
1810
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;
1830 }
1831 }
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;
1835 } else if (!x_sign
1836 && (is_midpoint_gt_even || is_inexact_lt_midpoint)) {
1837 Cstar.w[0] = Cstar.w[0] + 1;
1838 } else {
1839 ; // the result is already correct
1840 }
1841 if (x_sign)
1842 res = -Cstar.w[0];
1843 else
1844 res = Cstar.w[0];
1845 } else if (exp == 0) {
1846 // 1 <= q <= 10
1847 // res = +/-C (exact)
1848 if (x_sign)
1849 res = -C1.w[0];
1850 else
1851 res = C1.w[0];
1852 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1853 // res = +/-C * 10^exp (exact)
1854 if (x_sign)
1855 res = -C1.w[0] * ten2k64[exp];
1856 else
1857 res = C1.w[0] * ten2k64[exp];
1858 }
1859 }
1860 }
1861
1862 BID_RETURN (res);
1863 }
1864
1865 /*****************************************************************************
1866 * BID128_to_int32_xceil
1867 ****************************************************************************/
1868
1869 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_xceil, x)
1870
1871 int res;
1872 UINT64 x_sign;
1873 UINT64 x_exp;
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;
1879 int q, ind, shift;
1880 UINT128 C1, C;
1881 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
1882 UINT256 fstar;
1883 UINT256 P256;
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;
1888
1889 // unpack x
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;
1893 C1.w[0] = x.w[0];
1894
1895 // check for NaN or Infinity
1896 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1897 // x is 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
1900 // set invalid flag
1901 *pfpsf |= INVALID_EXCEPTION;
1902 // return Integer Indefinite
1903 res = 0x80000000;
1904 } else { // x is QNaN
1905 // set invalid flag
1906 *pfpsf |= INVALID_EXCEPTION;
1907 // return Integer Indefinite
1908 res = 0x80000000;
1909 }
1910 BID_RETURN (res);
1911 } else { // x is not a NaN, so it must be infinity
1912 if (!x_sign) { // x is +inf
1913 // set invalid flag
1914 *pfpsf |= INVALID_EXCEPTION;
1915 // return Integer Indefinite
1916 res = 0x80000000;
1917 } else { // x is -inf
1918 // set invalid flag
1919 *pfpsf |= INVALID_EXCEPTION;
1920 // return Integer Indefinite
1921 res = 0x80000000;
1922 }
1923 BID_RETURN (res);
1924 }
1925 }
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)) {
1931 res = 0x00000000;
1932 BID_RETURN (res);
1933 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1934 // x is 0
1935 res = 0x00000000;
1936 BID_RETURN (res);
1937 } else { // x is not special and is not zero
1938
1939 // q = nr. of decimal digits in x
1940 // determine first the nr. of bits in x
1941 if (C1.w[1] == 0) {
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
1946 x_nr_bits =
1947 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1948 } else { // x < 2^32
1949 tmp1.d = (double) (C1.w[0]); // exact conversion
1950 x_nr_bits =
1951 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1952 }
1953 } else { // if x < 2^53
1954 tmp1.d = (double) C1.w[0]; // exact conversion
1955 x_nr_bits =
1956 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1957 }
1958 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1959 tmp1.d = (double) C1.w[1]; // exact conversion
1960 x_nr_bits =
1961 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1962 }
1963 q = nr_digits[x_nr_bits - 1].digits;
1964 if (q == 0) {
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))
1969 q++;
1970 }
1971 exp = (x_exp >> 49) - 6176;
1972 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1973 // set invalid flag
1974 *pfpsf |= INVALID_EXCEPTION;
1975 // return Integer Indefinite
1976 res = 0x80000000;
1977 BID_RETURN (res);
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
1988 if (q <= 11) {
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) {
1992 // set invalid flag
1993 *pfpsf |= INVALID_EXCEPTION;
1994 // return Integer Indefinite
1995 res = 0x80000000;
1996 BID_RETURN (res);
1997 }
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]);
2009 }
2010 if (C1.w[1] > C.w[1]
2011 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2012 // set invalid flag
2013 *pfpsf |= INVALID_EXCEPTION;
2014 // return Integer Indefinite
2015 res = 0x80000000;
2016 BID_RETURN (res);
2017 }
2018 // else cases that can be rounded to a 32-bit int fall through
2019 // to '1 <= q + exp <= 10'
2020 }
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
2025 if (q <= 11) {
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) {
2029 // set invalid flag
2030 *pfpsf |= INVALID_EXCEPTION;
2031 // return Integer Indefinite
2032 res = 0x80000000;
2033 BID_RETURN (res);
2034 }
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
2040 // (scale 2^31 up)
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]);
2046 }
2047 if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
2048 // set invalid flag
2049 *pfpsf |= INVALID_EXCEPTION;
2050 // return Integer Indefinite
2051 res = 0x80000000;
2052 BID_RETURN (res);
2053 }
2054 // else cases that can be rounded to a 32-bit int fall through
2055 // to '1 <= q + exp <= 10'
2056 }
2057 }
2058 }
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)
2063 // set inexact flag
2064 *pfpsf |= INEXACT_EXCEPTION;
2065 // return 0
2066 if (x_sign)
2067 res = 0x00000000;
2068 else
2069 res = 0x00000001;
2070 BID_RETURN (res);
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
2078 tmp64 = C1.w[0];
2079 if (ind <= 19) {
2080 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
2081 } else {
2082 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
2083 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
2084 }
2085 if (C1.w[0] < tmp64)
2086 C1.w[1]++;
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[]
2091 // 1 <= x <= 33
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];
2099 fstar.w[3] = 0;
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
2104 Cstar.w[1] = 0;
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];
2110 }
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)
2118 // else
2119 // C* = floor(C*) (logical right shift; C has p decimal digits,
2120 // correct by Property 1)
2121 // n = C* * 10^(e+x)
2122
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
2126 Cstar.w[0] =
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
2131 }
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
2137 if (ind - 1 <= 2) {
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;
2151 }
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])
2162 tmp64A--;
2163 if (tmp64A || tmp64
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;
2175 }
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;
2195 }
2196 }
2197
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;
2217 }
2218 }
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;
2222 } else if (!x_sign
2223 && (is_midpoint_gt_even || is_inexact_lt_midpoint)) {
2224 Cstar.w[0] = Cstar.w[0] + 1;
2225 } else {
2226 ; // the result is already correct
2227 }
2228 if (x_sign)
2229 res = -Cstar.w[0];
2230 else
2231 res = Cstar.w[0];
2232 } else if (exp == 0) {
2233 // 1 <= q <= 10
2234 // res = +/-C (exact)
2235 if (x_sign)
2236 res = -C1.w[0];
2237 else
2238 res = C1.w[0];
2239 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2240 // res = +/-C * 10^exp (exact)
2241 if (x_sign)
2242 res = -C1.w[0] * ten2k64[exp];
2243 else
2244 res = C1.w[0] * ten2k64[exp];
2245 }
2246 }
2247 }
2248
2249 BID_RETURN (res);
2250 }
2251
2252 /*****************************************************************************
2253 * BID128_to_int32_int
2254 ****************************************************************************/
2255
2256 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_int, x)
2257
2258 int res;
2259 UINT64 x_sign;
2260 UINT64 x_exp;
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;
2266 int q, ind, shift;
2267 UINT128 C1, C;
2268 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
2269 UINT256 fstar;
2270 UINT256 P256;
2271 int is_inexact_gt_midpoint = 0;
2272 int is_midpoint_lt_even = 0;
2273
2274 // unpack x
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;
2278 C1.w[0] = x.w[0];
2279
2280 // check for NaN or Infinity
2281 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
2282 // x is 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
2285 // set invalid flag
2286 *pfpsf |= INVALID_EXCEPTION;
2287 // return Integer Indefinite
2288 res = 0x80000000;
2289 } else { // x is QNaN
2290 // set invalid flag
2291 *pfpsf |= INVALID_EXCEPTION;
2292 // return Integer Indefinite
2293 res = 0x80000000;
2294 }
2295 BID_RETURN (res);
2296 } else { // x is not a NaN, so it must be infinity
2297 if (!x_sign) { // x is +inf
2298 // set invalid flag
2299 *pfpsf |= INVALID_EXCEPTION;
2300 // return Integer Indefinite
2301 res = 0x80000000;
2302 } else { // x is -inf
2303 // set invalid flag
2304 *pfpsf |= INVALID_EXCEPTION;
2305 // return Integer Indefinite
2306 res = 0x80000000;
2307 }
2308 BID_RETURN (res);
2309 }
2310 }
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)) {
2316 res = 0x00000000;
2317 BID_RETURN (res);
2318 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
2319 // x is 0
2320 res = 0x00000000;
2321 BID_RETURN (res);
2322 } else { // x is not special and is not zero
2323
2324 // q = nr. of decimal digits in x
2325 // determine first the nr. of bits in x
2326 if (C1.w[1] == 0) {
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
2331 x_nr_bits =
2332 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2333 } else { // x < 2^32
2334 tmp1.d = (double) (C1.w[0]); // exact conversion
2335 x_nr_bits =
2336 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2337 }
2338 } else { // if x < 2^53
2339 tmp1.d = (double) C1.w[0]; // exact conversion
2340 x_nr_bits =
2341 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2342 }
2343 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2344 tmp1.d = (double) C1.w[1]; // exact conversion
2345 x_nr_bits =
2346 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2347 }
2348 q = nr_digits[x_nr_bits - 1].digits;
2349 if (q == 0) {
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))
2354 q++;
2355 }
2356 exp = (x_exp >> 49) - 6176;
2357 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
2358 // set invalid flag
2359 *pfpsf |= INVALID_EXCEPTION;
2360 // return Integer Indefinite
2361 res = 0x80000000;
2362 BID_RETURN (res);
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
2373 if (q <= 11) {
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) {
2377 // set invalid flag
2378 *pfpsf |= INVALID_EXCEPTION;
2379 // return Integer Indefinite
2380 res = 0x80000000;
2381 BID_RETURN (res);
2382 }
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]);
2394 }
2395 if (C1.w[1] > C.w[1]
2396 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2397 // set invalid flag
2398 *pfpsf |= INVALID_EXCEPTION;
2399 // return Integer Indefinite
2400 res = 0x80000000;
2401 BID_RETURN (res);
2402 }
2403 // else cases that can be rounded to a 32-bit int fall through
2404 // to '1 <= q + exp <= 10'
2405 }
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
2410 if (q <= 11) {
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) {
2414 // set invalid flag
2415 *pfpsf |= INVALID_EXCEPTION;
2416 // return Integer Indefinite
2417 res = 0x80000000;
2418 BID_RETURN (res);
2419 }
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]);
2431 }
2432 if (C1.w[1] > C.w[1]
2433 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2434 // set invalid flag
2435 *pfpsf |= INVALID_EXCEPTION;
2436 // return Integer Indefinite
2437 res = 0x80000000;
2438 BID_RETURN (res);
2439 }
2440 // else cases that can be rounded to a 32-bit int fall through
2441 // to '1 <= q + exp <= 10'
2442 }
2443 }
2444 }
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)
2449 // return 0
2450 res = 0x00000000;
2451 BID_RETURN (res);
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
2459 tmp64 = C1.w[0];
2460 if (ind <= 19) {
2461 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
2462 } else {
2463 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
2464 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
2465 }
2466 if (C1.w[0] < tmp64)
2467 C1.w[1]++;
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[]
2472 // 1 <= x <= 33
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];
2480 fstar.w[3] = 0;
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
2485 Cstar.w[1] = 0;
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];
2491 }
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)
2499 // else
2500 // C* = floor(C*) (logical right shift; C has p decimal digits,
2501 // correct by Property 1)
2502 // n = C* * 10^(e+x)
2503
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
2507 Cstar.w[0] =
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
2512 }
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
2518 if (ind - 1 <= 2) {
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;
2527 }
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])
2538 tmp64A--;
2539 if (tmp64A || tmp64
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;
2546 }
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;
2561 }
2562 }
2563
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;
2580 }
2581 }
2582 // general correction for RZ
2583 if (is_midpoint_lt_even || is_inexact_gt_midpoint) {
2584 Cstar.w[0] = Cstar.w[0] - 1;
2585 } else {
2586 ; // exact, the result is already correct
2587 }
2588 if (x_sign)
2589 res = -Cstar.w[0];
2590 else
2591 res = Cstar.w[0];
2592 } else if (exp == 0) {
2593 // 1 <= q <= 10
2594 // res = +/-C (exact)
2595 if (x_sign)
2596 res = -C1.w[0];
2597 else
2598 res = C1.w[0];
2599 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2600 // res = +/-C * 10^exp (exact)
2601 if (x_sign)
2602 res = -C1.w[0] * ten2k64[exp];
2603 else
2604 res = C1.w[0] * ten2k64[exp];
2605 }
2606 }
2607 }
2608
2609 BID_RETURN (res);
2610 }
2611
2612 /*****************************************************************************
2613 * BID128_to_int32_xint
2614 ****************************************************************************/
2615
2616 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_xint, x)
2617
2618 int res;
2619 UINT64 x_sign;
2620 UINT64 x_exp;
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;
2626 int q, ind, shift;
2627 UINT128 C1, C;
2628 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
2629 UINT256 fstar;
2630 UINT256 P256;
2631 int is_inexact_gt_midpoint = 0;
2632 int is_midpoint_lt_even = 0;
2633
2634 // unpack x
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;
2638 C1.w[0] = x.w[0];
2639
2640 // check for NaN or Infinity
2641 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
2642 // x is 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
2645 // set invalid flag
2646 *pfpsf |= INVALID_EXCEPTION;
2647 // return Integer Indefinite
2648 res = 0x80000000;
2649 } else { // x is QNaN
2650 // set invalid flag
2651 *pfpsf |= INVALID_EXCEPTION;
2652 // return Integer Indefinite
2653 res = 0x80000000;
2654 }
2655 BID_RETURN (res);
2656 } else { // x is not a NaN, so it must be infinity
2657 if (!x_sign) { // x is +inf
2658 // set invalid flag
2659 *pfpsf |= INVALID_EXCEPTION;
2660 // return Integer Indefinite
2661 res = 0x80000000;
2662 } else { // x is -inf
2663 // set invalid flag
2664 *pfpsf |= INVALID_EXCEPTION;
2665 // return Integer Indefinite
2666 res = 0x80000000;
2667 }
2668 BID_RETURN (res);
2669 }
2670 }
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)) {
2676 res = 0x00000000;
2677 BID_RETURN (res);
2678 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
2679 // x is 0
2680 res = 0x00000000;
2681 BID_RETURN (res);
2682 } else { // x is not special and is not zero
2683
2684 // q = nr. of decimal digits in x
2685 // determine first the nr. of bits in x
2686 if (C1.w[1] == 0) {
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
2691 x_nr_bits =
2692 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2693 } else { // x < 2^32
2694 tmp1.d = (double) (C1.w[0]); // exact conversion
2695 x_nr_bits =
2696 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2697 }
2698 } else { // if x < 2^53
2699 tmp1.d = (double) C1.w[0]; // exact conversion
2700 x_nr_bits =
2701 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2702 }
2703 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2704 tmp1.d = (double) C1.w[1]; // exact conversion
2705 x_nr_bits =
2706 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2707 }
2708 q = nr_digits[x_nr_bits - 1].digits;
2709 if (q == 0) {
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))
2714 q++;
2715 }
2716 exp = (x_exp >> 49) - 6176;
2717 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
2718 // set invalid flag
2719 *pfpsf |= INVALID_EXCEPTION;
2720 // return Integer Indefinite
2721 res = 0x80000000;
2722 BID_RETURN (res);
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
2733 if (q <= 11) {
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) {
2737 // set invalid flag
2738 *pfpsf |= INVALID_EXCEPTION;
2739 // return Integer Indefinite
2740 res = 0x80000000;
2741 BID_RETURN (res);
2742 }
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]);
2754 }
2755 if (C1.w[1] > C.w[1]
2756 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2757 // set invalid flag
2758 *pfpsf |= INVALID_EXCEPTION;
2759 // return Integer Indefinite
2760 res = 0x80000000;
2761 BID_RETURN (res);
2762 }
2763 // else cases that can be rounded to a 32-bit int fall through
2764 // to '1 <= q + exp <= 10'
2765 }
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
2770 if (q <= 11) {
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) {
2774 // set invalid flag
2775 *pfpsf |= INVALID_EXCEPTION;
2776 // return Integer Indefinite
2777 res = 0x80000000;
2778 BID_RETURN (res);
2779 }
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]);
2791 }
2792 if (C1.w[1] > C.w[1]
2793 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2794 // set invalid flag
2795 *pfpsf |= INVALID_EXCEPTION;
2796 // return Integer Indefinite
2797 res = 0x80000000;
2798 BID_RETURN (res);
2799 }
2800 // else cases that can be rounded to a 32-bit int fall through
2801 // to '1 <= q + exp <= 10'
2802 }
2803 }
2804 }
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)
2809 // set inexact flag
2810 *pfpsf |= INEXACT_EXCEPTION;
2811 // return 0
2812 res = 0x00000000;
2813 BID_RETURN (res);
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
2821 tmp64 = C1.w[0];
2822 if (ind <= 19) {
2823 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
2824 } else {
2825 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
2826 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
2827 }
2828 if (C1.w[0] < tmp64)
2829 C1.w[1]++;
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[]
2834 // 1 <= x <= 33
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];
2842 fstar.w[3] = 0;
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
2847 Cstar.w[1] = 0;
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];
2853 }
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)
2861 // else
2862 // C* = floor(C*) (logical right shift; C has p decimal digits,
2863 // correct by Property 1)
2864 // n = C* * 10^(e+x)
2865
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
2869 Cstar.w[0] =
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
2874 }
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
2880 if (ind - 1 <= 2) {
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;
2893 }
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])
2904 tmp64A--;
2905 if (tmp64A || tmp64
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;
2916 }
2917 } else { // if 22 <= ind <= 33
2918 if (fstar.w[3] > onehalf128[ind - 1] ||
2919 (fstar.w[3] == onehalf128[ind - 1] && (fstar.w[2] ||
2920 fstar.w[1]
2921 || fstar.w[0]))) {
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;
2936 }
2937 }
2938
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;
2955 }
2956 }
2957 // general correction for RZ
2958 if (is_midpoint_lt_even || is_inexact_gt_midpoint) {
2959 Cstar.w[0] = Cstar.w[0] - 1;
2960 } else {
2961 ; // exact, the result is already correct
2962 }
2963 if (x_sign)
2964 res = -Cstar.w[0];
2965 else
2966 res = Cstar.w[0];
2967 } else if (exp == 0) {
2968 // 1 <= q <= 10
2969 // res = +/-C (exact)
2970 if (x_sign)
2971 res = -C1.w[0];
2972 else
2973 res = C1.w[0];
2974 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2975 // res = +/-C * 10^exp (exact)
2976 if (x_sign)
2977 res = -C1.w[0] * ten2k64[exp];
2978 else
2979 res = C1.w[0] * ten2k64[exp];
2980 }
2981 }
2982 }
2983
2984 BID_RETURN (res);
2985 }
2986
2987 /*****************************************************************************
2988 * BID128_to_int32_rninta
2989 ****************************************************************************/
2990
2991 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_rninta,
2992 x)
2993
2994 int res;
2995 UINT64 x_sign;
2996 UINT64 x_exp;
2997 int exp; // unbiased exponent
2998 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2999 UINT64 tmp64;
3000 BID_UI64DOUBLE tmp1;
3001 unsigned int x_nr_bits;
3002 int q, ind, shift;
3003 UINT128 C1, C;
3004 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
3005 UINT256 P256;
3006
3007 // unpack x
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;
3011 C1.w[0] = x.w[0];
3012
3013 // check for NaN or Infinity
3014 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
3015 // x is 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
3018 // set invalid flag
3019 *pfpsf |= INVALID_EXCEPTION;
3020 // return Integer Indefinite
3021 res = 0x80000000;
3022 } else { // x is QNaN
3023 // set invalid flag
3024 *pfpsf |= INVALID_EXCEPTION;
3025 // return Integer Indefinite
3026 res = 0x80000000;
3027 }
3028 BID_RETURN (res);
3029 } else { // x is not a NaN, so it must be infinity
3030 if (!x_sign) { // x is +inf
3031 // set invalid flag
3032 *pfpsf |= INVALID_EXCEPTION;
3033 // return Integer Indefinite
3034 res = 0x80000000;
3035 } else { // x is -inf
3036 // set invalid flag
3037 *pfpsf |= INVALID_EXCEPTION;
3038 // return Integer Indefinite
3039 res = 0x80000000;
3040 }
3041 BID_RETURN (res);
3042 }
3043 }
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)) {
3049 res = 0x00000000;
3050 BID_RETURN (res);
3051 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
3052 // x is 0
3053 res = 0x00000000;
3054 BID_RETURN (res);
3055 } else { // x is not special and is not zero
3056
3057 // q = nr. of decimal digits in x
3058 // determine first the nr. of bits in x
3059 if (C1.w[1] == 0) {
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
3064 x_nr_bits =
3065 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3066 } else { // x < 2^32
3067 tmp1.d = (double) (C1.w[0]); // exact conversion
3068 x_nr_bits =
3069 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3070 }
3071 } else { // if x < 2^53
3072 tmp1.d = (double) C1.w[0]; // exact conversion
3073 x_nr_bits =
3074 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3075 }
3076 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
3077 tmp1.d = (double) C1.w[1]; // exact conversion
3078 x_nr_bits =
3079 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3080 }
3081 q = nr_digits[x_nr_bits - 1].digits;
3082 if (q == 0) {
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))
3087 q++;
3088 }
3089 exp = (x_exp >> 49) - 6176;
3090 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
3091 // set invalid flag
3092 *pfpsf |= INVALID_EXCEPTION;
3093 // return Integer Indefinite
3094 res = 0x80000000;
3095 BID_RETURN (res);
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
3106 if (q <= 11) {
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) {
3110 // set invalid flag
3111 *pfpsf |= INVALID_EXCEPTION;
3112 // return Integer Indefinite
3113 res = 0x80000000;
3114 BID_RETURN (res);
3115 }
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]);
3127 }
3128 if (C1.w[1] > C.w[1]
3129 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
3130 // set invalid flag
3131 *pfpsf |= INVALID_EXCEPTION;
3132 // return Integer Indefinite
3133 res = 0x80000000;
3134 BID_RETURN (res);
3135 }
3136 // else cases that can be rounded to a 32-bit int fall through
3137 // to '1 <= q + exp <= 10'
3138 }
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
3143 if (q <= 11) {
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) {
3147 // set invalid flag
3148 *pfpsf |= INVALID_EXCEPTION;
3149 // return Integer Indefinite
3150 res = 0x80000000;
3151 BID_RETURN (res);
3152 }
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]);
3164 }
3165 if (C1.w[1] > C.w[1]
3166 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
3167 // set invalid flag
3168 *pfpsf |= INVALID_EXCEPTION;
3169 // return Integer Indefinite
3170 res = 0x80000000;
3171 BID_RETURN (res);
3172 }
3173 // else cases that can be rounded to a 32-bit int fall through
3174 // to '1 <= q + exp <= 10'
3175 }
3176 }
3177 }
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)
3181 // return 0
3182 res = 0x00000000;
3183 BID_RETURN (res);
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)
3186 // res = 0
3187 // else
3188 // res = +/-1
3189 ind = 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
3195 } else { // n > 0
3196 res = 0x00000001; // return +1
3197 }
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
3205 } else { // n > 0
3206 res = 0x00000001; // return +1
3207 }
3208 }
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
3216 tmp64 = C1.w[0];
3217 if (ind <= 19) {
3218 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
3219 } else {
3220 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
3221 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
3222 }
3223 if (C1.w[0] < tmp64)
3224 C1.w[1]++;
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[]
3229 // 1 <= x <= 33
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
3238 Cstar.w[1] = 0;
3239 Cstar.w[0] = P256.w[3];
3240 }
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)
3248 // else
3249 // C* = floor(C*) (logical right shift; C has p decimal digits,
3250 // correct by Property 1)
3251 // n = C* * 10^(e+x)
3252
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
3256 Cstar.w[0] =
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
3261 }
3262 // if the result was a midpoint, it was already rounded away from zero
3263 if (x_sign)
3264 res = -Cstar.w[0];
3265 else
3266 res = Cstar.w[0];
3267 // no need to check for midpoints - already rounded away from zero!
3268 } else if (exp == 0) {
3269 // 1 <= q <= 10
3270 // res = +/-C (exact)
3271 if (x_sign)
3272 res = -C1.w[0];
3273 else
3274 res = C1.w[0];
3275 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
3276 // res = +/-C * 10^exp (exact)
3277 if (x_sign)
3278 res = -C1.w[0] * ten2k64[exp];
3279 else
3280 res = C1.w[0] * ten2k64[exp];
3281 }
3282 }
3283 }
3284
3285 BID_RETURN (res);
3286 }
3287
3288 /*****************************************************************************
3289 * BID128_to_int32_xrninta
3290 ****************************************************************************/
3291
3292 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_xrninta,
3293 x)
3294
3295 int res;
3296 UINT64 x_sign;
3297 UINT64 x_exp;
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;
3303 int q, ind, shift;
3304 UINT128 C1, C;
3305 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
3306 UINT256 fstar;
3307 UINT256 P256;
3308
3309 // unpack x
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;
3313 C1.w[0] = x.w[0];
3314
3315 // check for NaN or Infinity
3316 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
3317 // x is 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
3320 // set invalid flag
3321 *pfpsf |= INVALID_EXCEPTION;
3322 // return Integer Indefinite
3323 res = 0x80000000;
3324 } else { // x is QNaN
3325 // set invalid flag
3326 *pfpsf |= INVALID_EXCEPTION;
3327 // return Integer Indefinite
3328 res = 0x80000000;
3329 }
3330 BID_RETURN (res);
3331 } else { // x is not a NaN, so it must be infinity
3332 if (!x_sign) { // x is +inf
3333 // set invalid flag
3334 *pfpsf |= INVALID_EXCEPTION;
3335 // return Integer Indefinite
3336 res = 0x80000000;
3337 } else { // x is -inf
3338 // set invalid flag
3339 *pfpsf |= INVALID_EXCEPTION;
3340 // return Integer Indefinite
3341 res = 0x80000000;
3342 }
3343 BID_RETURN (res);
3344 }
3345 }
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)) {
3351 res = 0x00000000;
3352 BID_RETURN (res);
3353 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
3354 // x is 0
3355 res = 0x00000000;
3356 BID_RETURN (res);
3357 } else { // x is not special and is not zero
3358
3359 // q = nr. of decimal digits in x
3360 // determine first the nr. of bits in x
3361 if (C1.w[1] == 0) {
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
3366 x_nr_bits =
3367 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3368 } else { // x < 2^32
3369 tmp1.d = (double) (C1.w[0]); // exact conversion
3370 x_nr_bits =
3371 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3372 }
3373 } else { // if x < 2^53
3374 tmp1.d = (double) C1.w[0]; // exact conversion
3375 x_nr_bits =
3376 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3377 }
3378 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
3379 tmp1.d = (double) C1.w[1]; // exact conversion
3380 x_nr_bits =
3381 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3382 }
3383 q = nr_digits[x_nr_bits - 1].digits;
3384 if (q == 0) {
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))
3389 q++;
3390 }
3391 exp = (x_exp >> 49) - 6176;
3392 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
3393 // set invalid flag
3394 *pfpsf |= INVALID_EXCEPTION;
3395 // return Integer Indefinite
3396 res = 0x80000000;
3397 BID_RETURN (res);
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
3408 if (q <= 11) {
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) {
3412 // set invalid flag
3413 *pfpsf |= INVALID_EXCEPTION;
3414 // return Integer Indefinite
3415 res = 0x80000000;
3416 BID_RETURN (res);
3417 }
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]);
3429 }
3430 if (C1.w[1] > C.w[1]
3431 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
3432 // set invalid flag
3433 *pfpsf |= INVALID_EXCEPTION;
3434 // return Integer Indefinite
3435 res = 0x80000000;
3436 BID_RETURN (res);
3437 }
3438 // else cases that can be rounded to a 32-bit int fall through
3439 // to '1 <= q + exp <= 10'
3440 }
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
3445 if (q <= 11) {
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) {
3449 // set invalid flag
3450 *pfpsf |= INVALID_EXCEPTION;
3451 // return Integer Indefinite
3452 res = 0x80000000;
3453 BID_RETURN (res);
3454 }
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]);
3466 }
3467 if (C1.w[1] > C.w[1]
3468 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
3469 // set invalid flag
3470 *pfpsf |= INVALID_EXCEPTION;
3471 // return Integer Indefinite
3472 res = 0x80000000;
3473 BID_RETURN (res);
3474 }
3475 // else cases that can be rounded to a 32-bit int fall through
3476 // to '1 <= q + exp <= 10'
3477 }
3478 }
3479 }
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)
3483 // set inexact flag
3484 *pfpsf |= INEXACT_EXCEPTION;
3485 // return 0
3486 res = 0x00000000;
3487 BID_RETURN (res);
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)
3490 // res = 0
3491 // else
3492 // res = +/-1
3493 ind = 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
3499 } else { // n > 0
3500 res = 0x00000001; // return +1
3501 }
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
3509 } else { // n > 0
3510 res = 0x00000001; // return +1
3511 }
3512 }
3513 // set inexact flag
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
3522 tmp64 = C1.w[0];
3523 if (ind <= 19) {
3524 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
3525 } else {
3526 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
3527 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
3528 }
3529 if (C1.w[0] < tmp64)
3530 C1.w[1]++;
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[]
3535 // 1 <= x <= 33
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];
3543 fstar.w[3] = 0;
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
3548 Cstar.w[1] = 0;
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];
3554 }
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)
3562 // else
3563 // C* = floor(C*) (logical right shift; C has p decimal digits,
3564 // correct by Property 1)
3565 // n = C* * 10^(e+x)
3566
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
3570 Cstar.w[0] =
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
3575 }
3576 // if the result was a midpoint, it was already rounded away from zero
3577 if (x_sign)
3578 res = -Cstar.w[0];
3579 else
3580 res = Cstar.w[0];
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
3586 if (ind - 1 <= 2) {
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;
3598 }
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])
3609 tmp64A--;
3610 if (tmp64A || tmp64
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;
3620 }
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;
3638 }
3639 }
3640 // no need to check for midpoints - already rounded away from zero!
3641 } else if (exp == 0) {
3642 // 1 <= q <= 10
3643 // res = +/-C (exact)
3644 if (x_sign)
3645 res = -C1.w[0];
3646 else
3647 res = C1.w[0];
3648 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
3649 // res = +/-C * 10^exp (exact)
3650 if (x_sign)
3651 res = -C1.w[0] * ten2k64[exp];
3652 else
3653 res = C1.w[0] * ten2k64[exp];
3654 }
3655 }
3656 }
3657
3658 BID_RETURN (res);
3659 }