Mercurial > hg > CbC > CbC_gcc

comparison libgcc/config/libbid/bid128_to_uint32.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_uint32_rnint
28 ****************************************************************************/
29
30 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
31 bid128_to_uint32_rnint, x)
32
33 unsigned int res;
34 UINT64 x_sign;
35 UINT64 x_exp;
36 int exp; // unbiased exponent
37 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
38 UINT64 tmp64;
39 BID_UI64DOUBLE tmp1;
40 unsigned int x_nr_bits;
41 int q, ind, shift;
42 UINT128 C1, C;
43 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
44 UINT256 fstar;
45 UINT256 P256;
46
47 // unpack x
48 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
49 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
50 C1.w[1] = x.w[1] & MASK_COEFF;
51 C1.w[0] = x.w[0];
52
53 // check for NaN or Infinity
54 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
55 // x is special
56 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
57 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
58 // set invalid flag
59 *pfpsf |= INVALID_EXCEPTION;
60 // return Integer Indefinite
61 res = 0x80000000;
62 } else { // x is QNaN
63 // set invalid flag
64 *pfpsf |= INVALID_EXCEPTION;
65 // return Integer Indefinite
66 res = 0x80000000;
67 }
68 BID_RETURN (res);
69 } else { // x is not a NaN, so it must be infinity
70 if (!x_sign) { // x is +inf
71 // set invalid flag
72 *pfpsf |= INVALID_EXCEPTION;
73 // return Integer Indefinite
74 res = 0x80000000;
75 } else { // x is -inf
76 // set invalid flag
77 *pfpsf |= INVALID_EXCEPTION;
78 // return Integer Indefinite
79 res = 0x80000000;
80 }
81 BID_RETURN (res);
82 }
83 }
84 // check for non-canonical values (after the check for special values)
85 if ((C1.w[1] > 0x0001ed09bead87c0ull)
86 || (C1.w[1] == 0x0001ed09bead87c0ull
87 && (C1.w[0] > 0x378d8e63ffffffffull))
88 || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
89 res = 0x00000000;
90 BID_RETURN (res);
91 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
92 // x is 0
93 res = 0x00000000;
94 BID_RETURN (res);
95 } else { // x is not special and is not zero
96
97 // q = nr. of decimal digits in x
98 // determine first the nr. of bits in x
99 if (C1.w[1] == 0) {
100 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
101 // split the 64-bit value in two 32-bit halves to avoid rounding errors
102 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
103 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
104 x_nr_bits =
105 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
106 } else { // x < 2^32
107 tmp1.d = (double) (C1.w[0]); // exact conversion
108 x_nr_bits =
109 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
110 }
111 } else { // if x < 2^53
112 tmp1.d = (double) C1.w[0]; // exact conversion
113 x_nr_bits =
114 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
115 }
116 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
117 tmp1.d = (double) C1.w[1]; // exact conversion
118 x_nr_bits =
119 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
120 }
121 q = nr_digits[x_nr_bits - 1].digits;
122 if (q == 0) {
123 q = nr_digits[x_nr_bits - 1].digits1;
124 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
125 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
126 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
127 q++;
128 }
129 exp = (x_exp >> 49) - 6176;
130 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
131 // set invalid flag
132 *pfpsf |= INVALID_EXCEPTION;
133 // return Integer Indefinite
134 res = 0x80000000;
135 BID_RETURN (res);
136 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
137 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
138 // so x rounded to an integer may or may not fit in an unsigned 32-bit int
139 // the cases that do not fit are identified here; the ones that fit
140 // fall through and will be handled with other cases further,
141 // under '1 <= q + exp <= 10'
142 if (x_sign) { // if n < 0 and q + exp = 10
143 // if n < -1/2 then n cannot be converted to uint32 with RN
144 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 1/2
145 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x05, 1<=q<=34
146 if (q <= 11) {
147 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
148 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
149 if (tmp64 > 0x05ull) {
150 // set invalid flag
151 *pfpsf |= INVALID_EXCEPTION;
152 // return Integer Indefinite
153 res = 0x80000000;
154 BID_RETURN (res);
155 }
156 // else cases that can be rounded to a 32-bit int fall through
157 // to '1 <= q + exp <= 10'
158 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
159 // 0.c(0)c(1)...c(q-1) * 10^11 > 0x05 <=>
160 // C > 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
161 // (scale 1/2 up)
162 tmp64 = 0x05ull;
163 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
164 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
165 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
166 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
167 }
168 if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
169 // set invalid flag
170 *pfpsf |= INVALID_EXCEPTION;
171 // return Integer Indefinite
172 res = 0x80000000;
173 BID_RETURN (res);
174 }
175 // else cases that can be rounded to a 32-bit int fall through
176 // to '1 <= q + exp <= 10'
177 }
178 } else { // if n > 0 and q + exp = 10
179 // if n >= 2^32 - 1/2 then n is too large
180 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
181 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
182 if (q <= 11) {
183 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
184 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
185 if (tmp64 >= 0x9fffffffbull) {
186 // set invalid flag
187 *pfpsf |= INVALID_EXCEPTION;
188 // return Integer Indefinite
189 res = 0x80000000;
190 BID_RETURN (res);
191 }
192 // else cases that can be rounded to a 32-bit int fall through
193 // to '1 <= q + exp <= 10'
194 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
195 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
196 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
197 // (scale 2^32-1/2 up)
198 tmp64 = 0x9fffffffbull;
199 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
200 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
201 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
202 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
203 }
204 if (C1.w[1] > C.w[1]
205 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
206 // set invalid flag
207 *pfpsf |= INVALID_EXCEPTION;
208 // return Integer Indefinite
209 res = 0x80000000;
210 BID_RETURN (res);
211 }
212 // else cases that can be rounded to a 32-bit int fall through
213 // to '1 <= q + exp <= 10'
214 }
215 }
216 }
217 // n is not too large to be converted to int32: -1/2 <= n < 2^32 - 1/2
218 // Note: some of the cases tested for above fall through to this point
219 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
220 // return 0
221 res = 0x00000000;
222 BID_RETURN (res);
223 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
224 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
225 // res = 0
226 // else if x > 0
227 // res = +1
228 // else // if x < 0
229 // invalid exc
230 ind = q - 1;
231 if (ind <= 18) { // 0 <= ind <= 18
232 if ((C1.w[1] == 0) && (C1.w[0] <= midpoint64[ind])) {
233 res = 0x00000000; // return 0
234 } else if (!x_sign) { // n > 0
235 res = 0x00000001; // return +1
236 } else {
237 res = 0x80000000;
238 *pfpsf |= INVALID_EXCEPTION;
239 }
240 } else { // 19 <= ind <= 33
241 if ((C1.w[1] < midpoint128[ind - 19].w[1])
242 || ((C1.w[1] == midpoint128[ind - 19].w[1])
243 && (C1.w[0] <= midpoint128[ind - 19].w[0]))) {
244 res = 0x00000000; // return 0
245 } else if (!x_sign) { // n > 0
246 res = 0x00000001; // return +1
247 } else {
248 res = 0x80000000;
249 *pfpsf |= INVALID_EXCEPTION;
250 }
251 }
252 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
253 if (x_sign) { // x <= -1
254 // set invalid flag
255 *pfpsf |= INVALID_EXCEPTION;
256 // return Integer Indefinite
257 res = 0x80000000;
258 BID_RETURN (res);
259 }
260 // 1 <= x < 2^32-1/2 so x can be rounded
261 // to nearest to a 32-bit unsigned integer
262 if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
263 ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
264 // chop off ind digits from the lower part of C1
265 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
266 tmp64 = C1.w[0];
267 if (ind <= 19) {
268 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
269 } else {
270 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
271 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
272 }
273 if (C1.w[0] < tmp64)
274 C1.w[1]++;
275 // calculate C* and f*
276 // C* is actually floor(C*) in this case
277 // C* and f* need shifting and masking, as shown by
278 // shiftright128[] and maskhigh128[]
279 // 1 <= x <= 33
280 // kx = 10^(-x) = ten2mk128[ind - 1]
281 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
282 // the approximation of 10^(-x) was rounded up to 118 bits
283 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
284 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
285 Cstar.w[1] = P256.w[3];
286 Cstar.w[0] = P256.w[2];
287 fstar.w[3] = 0;
288 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
289 fstar.w[1] = P256.w[1];
290 fstar.w[0] = P256.w[0];
291 } else { // 22 <= ind - 1 <= 33
292 Cstar.w[1] = 0;
293 Cstar.w[0] = P256.w[3];
294 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
295 fstar.w[2] = P256.w[2];
296 fstar.w[1] = P256.w[1];
297 fstar.w[0] = P256.w[0];
298 }
299 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
300 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
301 // if (0 < f* < 10^(-x)) then the result is a midpoint
302 // if floor(C*) is even then C* = floor(C*) - logical right
303 // shift; C* has p decimal digits, correct by Prop. 1)
304 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
305 // shift; C* has p decimal digits, correct by Pr. 1)
306 // else
307 // C* = floor(C*) (logical right shift; C has p decimal digits,
308 // correct by Property 1)
309 // n = C* * 10^(e+x)
310
311 // shift right C* by Ex-128 = shiftright128[ind]
312 shift = shiftright128[ind - 1]; // 0 <= shift <= 102
313 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
314 Cstar.w[0] =
315 (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
316 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
317 } else { // 22 <= ind - 1 <= 33
318 Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
319 }
320 // if the result was a midpoint it was rounded away from zero, so
321 // it will need a correction
322 // check for midpoints
323 if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
324 && (fstar.w[1] || fstar.w[0])
325 && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
326 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
327 && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
328 // the result is a midpoint; round to nearest
329 if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
330 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
331 Cstar.w[0]--; // Cstar.w[0] is now even
332 } // else MP in [ODD, EVEN]
333 }
334 res = Cstar.w[0]; // the result is positive
335 } else if (exp == 0) {
336 // 1 <= q <= 10
337 // res = C (exact)
338 res = C1.w[0];
339 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
340 // res = C * 10^exp (exact)
341 res = C1.w[0] * ten2k64[exp];
342 }
343 }
344 }
345
346 BID_RETURN (res);
347 }
348
349 /*****************************************************************************
350 * BID128_to_uint32_xrnint
351 ****************************************************************************/
352
353 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
354 bid128_to_uint32_xrnint, x)
355
356 unsigned int res;
357 UINT64 x_sign;
358 UINT64 x_exp;
359 int exp; // unbiased exponent
360 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
361 UINT64 tmp64, tmp64A;
362 BID_UI64DOUBLE tmp1;
363 unsigned int x_nr_bits;
364 int q, ind, shift;
365 UINT128 C1, C;
366 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
367 UINT256 fstar;
368 UINT256 P256;
369 unsigned int tmp_inexact = 0;
370
371 // unpack x
372 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
373 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
374 C1.w[1] = x.w[1] & MASK_COEFF;
375 C1.w[0] = x.w[0];
376
377 // check for NaN or Infinity
378 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
379 // x is special
380 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
381 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
382 // set invalid flag
383 *pfpsf |= INVALID_EXCEPTION;
384 // return Integer Indefinite
385 res = 0x80000000;
386 } else { // x is QNaN
387 // set invalid flag
388 *pfpsf |= INVALID_EXCEPTION;
389 // return Integer Indefinite
390 res = 0x80000000;
391 }
392 BID_RETURN (res);
393 } else { // x is not a NaN, so it must be infinity
394 if (!x_sign) { // x is +inf
395 // set invalid flag
396 *pfpsf |= INVALID_EXCEPTION;
397 // return Integer Indefinite
398 res = 0x80000000;
399 } else { // x is -inf
400 // set invalid flag
401 *pfpsf |= INVALID_EXCEPTION;
402 // return Integer Indefinite
403 res = 0x80000000;
404 }
405 BID_RETURN (res);
406 }
407 }
408 // check for non-canonical values (after the check for special values)
409 if ((C1.w[1] > 0x0001ed09bead87c0ull)
410 || (C1.w[1] == 0x0001ed09bead87c0ull
411 && (C1.w[0] > 0x378d8e63ffffffffull))
412 || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
413 res = 0x00000000;
414 BID_RETURN (res);
415 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
416 // x is 0
417 res = 0x00000000;
418 BID_RETURN (res);
419 } else { // x is not special and is not zero
420
421 // q = nr. of decimal digits in x
422 // determine first the nr. of bits in x
423 if (C1.w[1] == 0) {
424 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
425 // split the 64-bit value in two 32-bit halves to avoid rounding errors
426 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
427 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
428 x_nr_bits =
429 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
430 } else { // x < 2^32
431 tmp1.d = (double) (C1.w[0]); // exact conversion
432 x_nr_bits =
433 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
434 }
435 } else { // if x < 2^53
436 tmp1.d = (double) C1.w[0]; // exact conversion
437 x_nr_bits =
438 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
439 }
440 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
441 tmp1.d = (double) C1.w[1]; // exact conversion
442 x_nr_bits =
443 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
444 }
445 q = nr_digits[x_nr_bits - 1].digits;
446 if (q == 0) {
447 q = nr_digits[x_nr_bits - 1].digits1;
448 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
449 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
450 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
451 q++;
452 }
453 exp = (x_exp >> 49) - 6176;
454 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
455 // set invalid flag
456 *pfpsf |= INVALID_EXCEPTION;
457 // return Integer Indefinite
458 res = 0x80000000;
459 BID_RETURN (res);
460 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
461 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
462 // so x rounded to an integer may or may not fit in an unsigned 32-bit int
463 // the cases that do not fit are identified here; the ones that fit
464 // fall through and will be handled with other cases further,
465 // under '1 <= q + exp <= 10'
466 if (x_sign) { // if n < 0 and q + exp = 10
467 // if n < -1/2 then n cannot be converted to uint32 with RN
468 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 1/2
469 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x05, 1<=q<=34
470 if (q <= 11) {
471 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
472 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
473 if (tmp64 > 0x05ull) {
474 // set invalid flag
475 *pfpsf |= INVALID_EXCEPTION;
476 // return Integer Indefinite
477 res = 0x80000000;
478 BID_RETURN (res);
479 }
480 // else cases that can be rounded to a 32-bit int fall through
481 // to '1 <= q + exp <= 10'
482 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
483 // 0.c(0)c(1)...c(q-1) * 10^11 > 0x05 <=>
484 // C > 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
485 // (scale 1/2 up)
486 tmp64 = 0x05ull;
487 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
488 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
489 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
490 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
491 }
492 if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
493 // set invalid flag
494 *pfpsf |= INVALID_EXCEPTION;
495 // return Integer Indefinite
496 res = 0x80000000;
497 BID_RETURN (res);
498 }
499 // else cases that can be rounded to a 32-bit int fall through
500 // to '1 <= q + exp <= 10'
501 }
502 } else { // if n > 0 and q + exp = 10
503 // if n >= 2^32 - 1/2 then n is too large
504 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
505 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
506 if (q <= 11) {
507 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
508 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
509 if (tmp64 >= 0x9fffffffbull) {
510 // set invalid flag
511 *pfpsf |= INVALID_EXCEPTION;
512 // return Integer Indefinite
513 res = 0x80000000;
514 BID_RETURN (res);
515 }
516 // else cases that can be rounded to a 32-bit int fall through
517 // to '1 <= q + exp <= 10'
518 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
519 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
520 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
521 // (scale 2^32-1/2 up)
522 tmp64 = 0x9fffffffbull;
523 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
524 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
525 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
526 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
527 }
528 if (C1.w[1] > C.w[1]
529 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
530 // set invalid flag
531 *pfpsf |= INVALID_EXCEPTION;
532 // return Integer Indefinite
533 res = 0x80000000;
534 BID_RETURN (res);
535 }
536 // else cases that can be rounded to a 32-bit int fall through
537 // to '1 <= q + exp <= 10'
538 }
539 }
540 }
541 // n is not too large to be converted to int32: -1/2 <= n < 2^32 - 1/2
542 // Note: some of the cases tested for above fall through to this point
543 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
544 // set inexact flag
545 *pfpsf |= INEXACT_EXCEPTION;
546 // return 0
547 res = 0x00000000;
548 BID_RETURN (res);
549 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
550 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
551 // res = 0
552 // else if x > 0
553 // res = +1
554 // else // if x < 0
555 // invalid exc
556 ind = q - 1;
557 if (ind <= 18) { // 0 <= ind <= 18
558 if ((C1.w[1] == 0) && (C1.w[0] <= midpoint64[ind])) {
559 res = 0x00000000; // return 0
560 } else if (!x_sign) { // n > 0
561 res = 0x00000001; // return +1
562 } else {
563 res = 0x80000000;
564 *pfpsf |= INVALID_EXCEPTION;
565 BID_RETURN (res);
566 }
567 } else { // 19 <= ind <= 33
568 if ((C1.w[1] < midpoint128[ind - 19].w[1])
569 || ((C1.w[1] == midpoint128[ind - 19].w[1])
570 && (C1.w[0] <= midpoint128[ind - 19].w[0]))) {
571 res = 0x00000000; // return 0
572 } else if (!x_sign) { // n > 0
573 res = 0x00000001; // return +1
574 } else {
575 res = 0x80000000;
576 *pfpsf |= INVALID_EXCEPTION;
577 BID_RETURN (res);
578 }
579 }
580 // set inexact flag
581 *pfpsf |= INEXACT_EXCEPTION;
582 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
583 if (x_sign) { // x <= -1
584 // set invalid flag
585 *pfpsf |= INVALID_EXCEPTION;
586 // return Integer Indefinite
587 res = 0x80000000;
588 BID_RETURN (res);
589 }
590 // 1 <= x < 2^32-1/2 so x can be rounded
591 // to nearest to a 32-bit unsigned integer
592 if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
593 ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
594 // chop off ind digits from the lower part of C1
595 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
596 tmp64 = C1.w[0];
597 if (ind <= 19) {
598 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
599 } else {
600 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
601 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
602 }
603 if (C1.w[0] < tmp64)
604 C1.w[1]++;
605 // calculate C* and f*
606 // C* is actually floor(C*) in this case
607 // C* and f* need shifting and masking, as shown by
608 // shiftright128[] and maskhigh128[]
609 // 1 <= x <= 33
610 // kx = 10^(-x) = ten2mk128[ind - 1]
611 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
612 // the approximation of 10^(-x) was rounded up to 118 bits
613 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
614 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
615 Cstar.w[1] = P256.w[3];
616 Cstar.w[0] = P256.w[2];
617 fstar.w[3] = 0;
618 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
619 fstar.w[1] = P256.w[1];
620 fstar.w[0] = P256.w[0];
621 } else { // 22 <= ind - 1 <= 33
622 Cstar.w[1] = 0;
623 Cstar.w[0] = P256.w[3];
624 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
625 fstar.w[2] = P256.w[2];
626 fstar.w[1] = P256.w[1];
627 fstar.w[0] = P256.w[0];
628 }
629 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
630 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
631 // if (0 < f* < 10^(-x)) then the result is a midpoint
632 // if floor(C*) is even then C* = floor(C*) - logical right
633 // shift; C* has p decimal digits, correct by Prop. 1)
634 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
635 // shift; C* has p decimal digits, correct by Pr. 1)
636 // else
637 // C* = floor(C*) (logical right shift; C has p decimal digits,
638 // correct by Property 1)
639 // n = C* * 10^(e+x)
640
641 // shift right C* by Ex-128 = shiftright128[ind]
642 shift = shiftright128[ind - 1]; // 0 <= shift <= 102
643 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
644 Cstar.w[0] =
645 (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
646 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
647 } else { // 22 <= ind - 1 <= 33
648 Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
649 }
650 // determine inexactness of the rounding of C*
651 // if (0 < f* - 1/2 < 10^(-x)) then
652 // the result is exact
653 // else // if (f* - 1/2 > T*) then
654 // the result is inexact
655 if (ind - 1 <= 2) {
656 if (fstar.w[1] > 0x8000000000000000ull ||
657 (fstar.w[1] == 0x8000000000000000ull
658 && fstar.w[0] > 0x0ull)) {
659 // f* > 1/2 and the result may be exact
660 tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
661 if (tmp64 > ten2mk128trunc[ind - 1].w[1]
662 || (tmp64 == ten2mk128trunc[ind - 1].w[1]
663 && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
664 // set the inexact flag
665 // *pfpsf |= INEXACT_EXCEPTION;
666 tmp_inexact = 1;
667 } // else the result is exact
668 } else { // the result is inexact; f2* <= 1/2
669 // set the inexact flag
670 // *pfpsf |= INEXACT_EXCEPTION;
671 tmp_inexact = 1;
672 }
673 } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
674 if (fstar.w[3] > 0x0 ||
675 (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
676 (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
677 (fstar.w[1] || fstar.w[0]))) {
678 // f2* > 1/2 and the result may be exact
679 // Calculate f2* - 1/2
680 tmp64 = fstar.w[2] - onehalf128[ind - 1];
681 tmp64A = fstar.w[3];
682 if (tmp64 > fstar.w[2])
683 tmp64A--;
684 if (tmp64A || tmp64
685 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
686 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
687 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
688 // set the inexact flag
689 // *pfpsf |= INEXACT_EXCEPTION;
690 tmp_inexact = 1;
691 } // else the result is exact
692 } else { // the result is inexact; f2* <= 1/2
693 // set the inexact flag
694 // *pfpsf |= INEXACT_EXCEPTION;
695 tmp_inexact = 1;
696 }
697 } else { // if 22 <= ind <= 33
698 if (fstar.w[3] > onehalf128[ind - 1] ||
699 (fstar.w[3] == onehalf128[ind - 1] &&
700 (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
701 // f2* > 1/2 and the result may be exact
702 // Calculate f2* - 1/2
703 tmp64 = fstar.w[3] - onehalf128[ind - 1];
704 if (tmp64 || fstar.w[2]
705 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
706 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
707 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
708 // set the inexact flag
709 // *pfpsf |= INEXACT_EXCEPTION;
710 tmp_inexact = 1;
711 } // else the result is exact
712 } else { // the result is inexact; f2* <= 1/2
713 // set the inexact flag
714 // *pfpsf |= INEXACT_EXCEPTION;
715 tmp_inexact = 1;
716 }
717 }
718
719 // if the result was a midpoint it was rounded away from zero, so
720 // it will need a correction
721 // check for midpoints
722 if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
723 && (fstar.w[1] || fstar.w[0])
724 && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
725 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
726 && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
727 // the result is a midpoint; round to nearest
728 if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
729 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
730 Cstar.w[0]--; // Cstar.w[0] is now even
731 } // else MP in [ODD, EVEN]
732 }
733 res = Cstar.w[0]; // the result is positive
734 if (tmp_inexact)
735 *pfpsf |= INEXACT_EXCEPTION;
736 } else if (exp == 0) {
737 // 1 <= q <= 10
738 // res = C (exact)
739 res = C1.w[0];
740 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
741 // res = C * 10^exp (exact)
742 res = C1.w[0] * ten2k64[exp];
743 }
744 }
745 }
746
747 BID_RETURN (res);
748 }
749
750 /*****************************************************************************
751 * BID128_to_uint32_floor
752 ****************************************************************************/
753
754 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
755 bid128_to_uint32_floor, x)
756
757 unsigned int res;
758 UINT64 x_sign;
759 UINT64 x_exp;
760 int exp; // unbiased exponent
761 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
762 UINT64 tmp64, tmp64A;
763 BID_UI64DOUBLE tmp1;
764 unsigned int x_nr_bits;
765 int q, ind, shift;
766 UINT128 C1, C;
767 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
768 UINT256 fstar;
769 UINT256 P256;
770 int is_inexact_gt_midpoint = 0;
771 int is_midpoint_lt_even = 0;
772
773 // unpack x
774 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
775 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
776 C1.w[1] = x.w[1] & MASK_COEFF;
777 C1.w[0] = x.w[0];
778
779 // check for NaN or Infinity
780 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
781 // x is special
782 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
783 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
784 // set invalid flag
785 *pfpsf |= INVALID_EXCEPTION;
786 // return Integer Indefinite
787 res = 0x80000000;
788 } else { // x is QNaN
789 // set invalid flag
790 *pfpsf |= INVALID_EXCEPTION;
791 // return Integer Indefinite
792 res = 0x80000000;
793 }
794 BID_RETURN (res);
795 } else { // x is not a NaN, so it must be infinity
796 if (!x_sign) { // x is +inf
797 // set invalid flag
798 *pfpsf |= INVALID_EXCEPTION;
799 // return Integer Indefinite
800 res = 0x80000000;
801 } else { // x is -inf
802 // set invalid flag
803 *pfpsf |= INVALID_EXCEPTION;
804 // return Integer Indefinite
805 res = 0x80000000;
806 }
807 BID_RETURN (res);
808 }
809 }
810 // check for non-canonical values (after the check for special values)
811 if ((C1.w[1] > 0x0001ed09bead87c0ull)
812 || (C1.w[1] == 0x0001ed09bead87c0ull
813 && (C1.w[0] > 0x378d8e63ffffffffull))
814 || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
815 res = 0x00000000;
816 BID_RETURN (res);
817 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
818 // x is 0
819 res = 0x00000000;
820 BID_RETURN (res);
821 } else { // x is not special and is not zero
822 // x < 0 is invalid
823 if (x_sign) {
824 // set invalid flag
825 *pfpsf |= INVALID_EXCEPTION;
826 // return Integer Indefinite
827 res = 0x80000000;
828 BID_RETURN (res);
829 }
830 // x > 0 from this point on
831 // q = nr. of decimal digits in x
832 // determine first the nr. of bits in x
833 if (C1.w[1] == 0) {
834 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
835 // split the 64-bit value in two 32-bit halves to avoid rounding errors
836 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
837 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
838 x_nr_bits =
839 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
840 } else { // x < 2^32
841 tmp1.d = (double) (C1.w[0]); // exact conversion
842 x_nr_bits =
843 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
844 }
845 } else { // if x < 2^53
846 tmp1.d = (double) C1.w[0]; // exact conversion
847 x_nr_bits =
848 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
849 }
850 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
851 tmp1.d = (double) C1.w[1]; // exact conversion
852 x_nr_bits =
853 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
854 }
855 q = nr_digits[x_nr_bits - 1].digits;
856 if (q == 0) {
857 q = nr_digits[x_nr_bits - 1].digits1;
858 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
859 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
860 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
861 q++;
862 }
863 exp = (x_exp >> 49) - 6176;
864
865 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
866 // set invalid flag
867 *pfpsf |= INVALID_EXCEPTION;
868 // return Integer Indefinite
869 res = 0x80000000;
870 BID_RETURN (res);
871 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
872 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
873 // so x rounded to an integer may or may not fit in a signed 32-bit int
874 // the cases that do not fit are identified here; the ones that fit
875 // fall through and will be handled with other cases further,
876 // under '1 <= q + exp <= 10'
877 // n > 0 and q + exp = 10
878 // if n >= 2^32 then n is too large
879 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
880 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
881 if (q <= 11) {
882 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
883 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
884 if (tmp64 >= 0xa00000000ull) {
885 // set invalid flag
886 *pfpsf |= INVALID_EXCEPTION;
887 // return Integer Indefinite
888 res = 0x80000000;
889 BID_RETURN (res);
890 }
891 // else cases that can be rounded to a 32-bit int fall through
892 // to '1 <= q + exp <= 10'
893 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
894 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
895 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
896 // (scale 2^32 up)
897 tmp64 = 0xa00000000ull;
898 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
899 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
900 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
901 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
902 }
903 if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
904 // set invalid flag
905 *pfpsf |= INVALID_EXCEPTION;
906 // return Integer Indefinite
907 res = 0x80000000;
908 BID_RETURN (res);
909 }
910 // else cases that can be rounded to a 32-bit int fall through
911 // to '1 <= q + exp <= 10'
912 }
913 }
914 // n is not too large to be converted to int32: 0 <= n < 2^31
915 // Note: some of the cases tested for above fall through to this point
916 if ((q + exp) <= 0) {
917 // n = +0.0...c(0)c(1)...c(q-1) or n = +0.c(0)c(1)...c(q-1)
918 // return 0
919 res = 0x00000000;
920 BID_RETURN (res);
921 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
922 // 1 <= x < 2^32 so x can be rounded down to a 32-bit unsigned integer
923 if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
924 ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
925 // chop off ind digits from the lower part of C1
926 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
927 tmp64 = C1.w[0];
928 if (ind <= 19) {
929 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
930 } else {
931 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
932 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
933 }
934 if (C1.w[0] < tmp64)
935 C1.w[1]++;
936 // calculate C* and f*
937 // C* is actually floor(C*) in this case
938 // C* and f* need shifting and masking, as shown by
939 // shiftright128[] and maskhigh128[]
940 // 1 <= x <= 33
941 // kx = 10^(-x) = ten2mk128[ind - 1]
942 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
943 // the approximation of 10^(-x) was rounded up to 118 bits
944 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
945 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
946 Cstar.w[1] = P256.w[3];
947 Cstar.w[0] = P256.w[2];
948 fstar.w[3] = 0;
949 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
950 fstar.w[1] = P256.w[1];
951 fstar.w[0] = P256.w[0];
952 } else { // 22 <= ind - 1 <= 33
953 Cstar.w[1] = 0;
954 Cstar.w[0] = P256.w[3];
955 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
956 fstar.w[2] = P256.w[2];
957 fstar.w[1] = P256.w[1];
958 fstar.w[0] = P256.w[0];
959 }
960 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
961 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
962 // if (0 < f* < 10^(-x)) then the result is a midpoint
963 // if floor(C*) is even then C* = floor(C*) - logical right
964 // shift; C* has p decimal digits, correct by Prop. 1)
965 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
966 // shift; C* has p decimal digits, correct by Pr. 1)
967 // else
968 // C* = floor(C*) (logical right shift; C has p decimal digits,
969 // correct by Property 1)
970 // n = C* * 10^(e+x)
971
972 // shift right C* by Ex-128 = shiftright128[ind]
973 shift = shiftright128[ind - 1]; // 0 <= shift <= 102
974 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
975 Cstar.w[0] =
976 (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
977 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
978 } else { // 22 <= ind - 1 <= 33
979 Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
980 }
981 // determine inexactness of the rounding of C*
982 // if (0 < f* - 1/2 < 10^(-x)) then
983 // the result is exact
984 // else // if (f* - 1/2 > T*) then
985 // the result is inexact
986 if (ind - 1 <= 2) {
987 if (fstar.w[1] > 0x8000000000000000ull ||
988 (fstar.w[1] == 0x8000000000000000ull
989 && fstar.w[0] > 0x0ull)) {
990 // f* > 1/2 and the result may be exact
991 tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
992 if (tmp64 > ten2mk128trunc[ind - 1].w[1]
993 || (tmp64 == ten2mk128trunc[ind - 1].w[1]
994 && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
995 } // else the result is exact
996 } else { // the result is inexact; f2* <= 1/2
997 is_inexact_gt_midpoint = 1;
998 }
999 } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
1000 if (fstar.w[3] > 0x0 ||
1001 (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
1002 (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
1003 (fstar.w[1] || fstar.w[0]))) {
1004 // f2* > 1/2 and the result may be exact
1005 // Calculate f2* - 1/2
1006 tmp64 = fstar.w[2] - onehalf128[ind - 1];
1007 tmp64A = fstar.w[3];
1008 if (tmp64 > fstar.w[2])
1009 tmp64A--;
1010 if (tmp64A || tmp64
1011 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1012 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1013 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1014 } // else the result is exact
1015 } else { // the result is inexact; f2* <= 1/2
1016 is_inexact_gt_midpoint = 1;
1017 }
1018 } else { // if 22 <= ind <= 33
1019 if (fstar.w[3] > onehalf128[ind - 1] ||
1020 (fstar.w[3] == onehalf128[ind - 1] &&
1021 (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
1022 // f2* > 1/2 and the result may be exact
1023 // Calculate f2* - 1/2
1024 tmp64 = fstar.w[3] - onehalf128[ind - 1];
1025 if (tmp64 || fstar.w[2]
1026 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1027 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1028 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1029 } // else the result is exact
1030 } else { // the result is inexact; f2* <= 1/2
1031 is_inexact_gt_midpoint = 1;
1032 }
1033 }
1034
1035 // if the result was a midpoint it was rounded away from zero, so
1036 // it will need a correction
1037 // check for midpoints
1038 if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
1039 && (fstar.w[1] || fstar.w[0])
1040 && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
1041 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1042 && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
1043 // the result is a midpoint; round to nearest
1044 if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
1045 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
1046 Cstar.w[0]--; // Cstar.w[0] is now even
1047 is_inexact_gt_midpoint = 0;
1048 } else { // else MP in [ODD, EVEN]
1049 is_midpoint_lt_even = 1;
1050 is_inexact_gt_midpoint = 0;
1051 }
1052 }
1053 // general correction for RM
1054 if (is_midpoint_lt_even || is_inexact_gt_midpoint) {
1055 Cstar.w[0] = Cstar.w[0] - 1;
1056 } else {
1057 ; // the result is already correct
1058 }
1059 res = Cstar.w[0];
1060 } else if (exp == 0) {
1061 // 1 <= q <= 10
1062 // res = +C (exact)
1063 res = C1.w[0];
1064 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1065 // res = +C * 10^exp (exact)
1066 res = C1.w[0] * ten2k64[exp];
1067 }
1068 }
1069 }
1070
1071 BID_RETURN (res);
1072 }
1073
1074 /*****************************************************************************
1075 * BID128_to_uint32_xfloor
1076 ****************************************************************************/
1077
1078 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
1079 bid128_to_uint32_xfloor, x)
1080
1081 unsigned int res;
1082 UINT64 x_sign;
1083 UINT64 x_exp;
1084 int exp; // unbiased exponent
1085 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1086 UINT64 tmp64, tmp64A;
1087 BID_UI64DOUBLE tmp1;
1088 unsigned int x_nr_bits;
1089 int q, ind, shift;
1090 UINT128 C1, C;
1091 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
1092 UINT256 fstar;
1093 UINT256 P256;
1094 int is_inexact_gt_midpoint = 0;
1095 int is_midpoint_lt_even = 0;
1096
1097 // unpack x
1098 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1099 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
1100 C1.w[1] = x.w[1] & MASK_COEFF;
1101 C1.w[0] = x.w[0];
1102
1103 // check for NaN or Infinity
1104 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1105 // x is special
1106 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
1107 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
1108 // set invalid flag
1109 *pfpsf |= INVALID_EXCEPTION;
1110 // return Integer Indefinite
1111 res = 0x80000000;
1112 } else { // x is QNaN
1113 // set invalid flag
1114 *pfpsf |= INVALID_EXCEPTION;
1115 // return Integer Indefinite
1116 res = 0x80000000;
1117 }
1118 BID_RETURN (res);
1119 } else { // x is not a NaN, so it must be infinity
1120 if (!x_sign) { // x is +inf
1121 // set invalid flag
1122 *pfpsf |= INVALID_EXCEPTION;
1123 // return Integer Indefinite
1124 res = 0x80000000;
1125 } else { // x is -inf
1126 // set invalid flag
1127 *pfpsf |= INVALID_EXCEPTION;
1128 // return Integer Indefinite
1129 res = 0x80000000;
1130 }
1131 BID_RETURN (res);
1132 }
1133 }
1134 // check for non-canonical values (after the check for special values)
1135 if ((C1.w[1] > 0x0001ed09bead87c0ull)
1136 || (C1.w[1] == 0x0001ed09bead87c0ull
1137 && (C1.w[0] > 0x378d8e63ffffffffull))
1138 || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
1139 res = 0x00000000;
1140 BID_RETURN (res);
1141 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1142 // x is 0
1143 res = 0x00000000;
1144 BID_RETURN (res);
1145 } else { // x is not special and is not zero
1146 // x < 0 is invalid
1147 if (x_sign) {
1148 // set invalid flag
1149 *pfpsf |= INVALID_EXCEPTION;
1150 // return Integer Indefinite
1151 res = 0x80000000;
1152 BID_RETURN (res);
1153 }
1154 // x > 0 from this point on
1155 // q = nr. of decimal digits in x
1156 // determine first the nr. of bits in x
1157 if (C1.w[1] == 0) {
1158 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
1159 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1160 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
1161 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
1162 x_nr_bits =
1163 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1164 } else { // x < 2^32
1165 tmp1.d = (double) (C1.w[0]); // exact conversion
1166 x_nr_bits =
1167 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1168 }
1169 } else { // if x < 2^53
1170 tmp1.d = (double) C1.w[0]; // exact conversion
1171 x_nr_bits =
1172 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1173 }
1174 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1175 tmp1.d = (double) C1.w[1]; // exact conversion
1176 x_nr_bits =
1177 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1178 }
1179 q = nr_digits[x_nr_bits - 1].digits;
1180 if (q == 0) {
1181 q = nr_digits[x_nr_bits - 1].digits1;
1182 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
1183 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
1184 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1185 q++;
1186 }
1187 exp = (x_exp >> 49) - 6176;
1188 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1189 // set invalid flag
1190 *pfpsf |= INVALID_EXCEPTION;
1191 // return Integer Indefinite
1192 res = 0x80000000;
1193 BID_RETURN (res);
1194 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1195 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1196 // so x rounded to an integer may or may not fit in a signed 32-bit int
1197 // the cases that do not fit are identified here; the ones that fit
1198 // fall through and will be handled with other cases further,
1199 // under '1 <= q + exp <= 10'
1200 // n > 0 and q + exp = 10
1201 // if n >= 2^32 then n is too large
1202 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
1203 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
1204 if (q <= 11) {
1205 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
1206 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1207 if (tmp64 >= 0xa00000000ull) {
1208 // set invalid flag
1209 *pfpsf |= INVALID_EXCEPTION;
1210 // return Integer Indefinite
1211 res = 0x80000000;
1212 BID_RETURN (res);
1213 }
1214 // else cases that can be rounded to a 32-bit int fall through
1215 // to '1 <= q + exp <= 10'
1216 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1217 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
1218 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
1219 // (scale 2^32 up)
1220 tmp64 = 0xa00000000ull;
1221 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1222 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
1223 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1224 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
1225 }
1226 if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
1227 // set invalid flag
1228 *pfpsf |= INVALID_EXCEPTION;
1229 // return Integer Indefinite
1230 res = 0x80000000;
1231 BID_RETURN (res);
1232 }
1233 // else cases that can be rounded to a 32-bit int fall through
1234 // to '1 <= q + exp <= 10'
1235 }
1236 }
1237 // n is not too large to be converted to int32: 0 <= n < 2^31
1238 // Note: some of the cases tested for above fall through to this point
1239 if ((q + exp) <= 0) {
1240 // n = +0.0...c(0)c(1)...c(q-1) or n = +0.c(0)c(1)...c(q-1)
1241 // set inexact flag
1242 *pfpsf |= INEXACT_EXCEPTION;
1243 // return 0
1244 res = 0x00000000;
1245 BID_RETURN (res);
1246 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
1247 // 1 <= x < 2^32 so x can be rounded down to a 32-bit unsigned integer
1248 if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
1249 ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
1250 // chop off ind digits from the lower part of C1
1251 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
1252 tmp64 = C1.w[0];
1253 if (ind <= 19) {
1254 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1255 } else {
1256 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1257 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1258 }
1259 if (C1.w[0] < tmp64)
1260 C1.w[1]++;
1261 // calculate C* and f*
1262 // C* is actually floor(C*) in this case
1263 // C* and f* need shifting and masking, as shown by
1264 // shiftright128[] and maskhigh128[]
1265 // 1 <= x <= 33
1266 // kx = 10^(-x) = ten2mk128[ind - 1]
1267 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1268 // the approximation of 10^(-x) was rounded up to 118 bits
1269 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1270 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
1271 Cstar.w[1] = P256.w[3];
1272 Cstar.w[0] = P256.w[2];
1273 fstar.w[3] = 0;
1274 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1275 fstar.w[1] = P256.w[1];
1276 fstar.w[0] = P256.w[0];
1277 } else { // 22 <= ind - 1 <= 33
1278 Cstar.w[1] = 0;
1279 Cstar.w[0] = P256.w[3];
1280 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1281 fstar.w[2] = P256.w[2];
1282 fstar.w[1] = P256.w[1];
1283 fstar.w[0] = P256.w[0];
1284 }
1285 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1286 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1287 // if (0 < f* < 10^(-x)) then the result is a midpoint
1288 // if floor(C*) is even then C* = floor(C*) - logical right
1289 // shift; C* has p decimal digits, correct by Prop. 1)
1290 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
1291 // shift; C* has p decimal digits, correct by Pr. 1)
1292 // else
1293 // C* = floor(C*) (logical right shift; C has p decimal digits,
1294 // correct by Property 1)
1295 // n = C* * 10^(e+x)
1296
1297 // shift right C* by Ex-128 = shiftright128[ind]
1298 shift = shiftright128[ind - 1]; // 0 <= shift <= 102
1299 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
1300 Cstar.w[0] =
1301 (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
1302 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1303 } else { // 22 <= ind - 1 <= 33
1304 Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
1305 }
1306 // determine inexactness of the rounding of C*
1307 // if (0 < f* - 1/2 < 10^(-x)) then
1308 // the result is exact
1309 // else // if (f* - 1/2 > T*) then
1310 // the result is inexact
1311 if (ind - 1 <= 2) {
1312 if (fstar.w[1] > 0x8000000000000000ull ||
1313 (fstar.w[1] == 0x8000000000000000ull
1314 && fstar.w[0] > 0x0ull)) {
1315 // f* > 1/2 and the result may be exact
1316 tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
1317 if (tmp64 > ten2mk128trunc[ind - 1].w[1]
1318 || (tmp64 == ten2mk128trunc[ind - 1].w[1]
1319 && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
1320 // set the inexact flag
1321 *pfpsf |= INEXACT_EXCEPTION;
1322 } // else the result is exact
1323 } else { // the result is inexact; f2* <= 1/2
1324 // set the inexact flag
1325 *pfpsf |= INEXACT_EXCEPTION;
1326 is_inexact_gt_midpoint = 1;
1327 }
1328 } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
1329 if (fstar.w[3] > 0x0 ||
1330 (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
1331 (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
1332 (fstar.w[1] || fstar.w[0]))) {
1333 // f2* > 1/2 and the result may be exact
1334 // Calculate f2* - 1/2
1335 tmp64 = fstar.w[2] - onehalf128[ind - 1];
1336 tmp64A = fstar.w[3];
1337 if (tmp64 > fstar.w[2])
1338 tmp64A--;
1339 if (tmp64A || tmp64
1340 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1341 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1342 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1343 // set the inexact flag
1344 *pfpsf |= INEXACT_EXCEPTION;
1345 } // else the result is exact
1346 } else { // the result is inexact; f2* <= 1/2
1347 // set the inexact flag
1348 *pfpsf |= INEXACT_EXCEPTION;
1349 is_inexact_gt_midpoint = 1;
1350 }
1351 } else { // if 22 <= ind <= 33
1352 if (fstar.w[3] > onehalf128[ind - 1] ||
1353 (fstar.w[3] == onehalf128[ind - 1] &&
1354 (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
1355 // f2* > 1/2 and the result may be exact
1356 // Calculate f2* - 1/2
1357 tmp64 = fstar.w[3] - onehalf128[ind - 1];
1358 if (tmp64 || fstar.w[2]
1359 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1360 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1361 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1362 // set the inexact flag
1363 *pfpsf |= INEXACT_EXCEPTION;
1364 } // else the result is exact
1365 } else { // the result is inexact; f2* <= 1/2
1366 // set the inexact flag
1367 *pfpsf |= INEXACT_EXCEPTION;
1368 is_inexact_gt_midpoint = 1;
1369 }
1370 }
1371
1372 // if the result was a midpoint it was rounded away from zero, so
1373 // it will need a correction
1374 // check for midpoints
1375 if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
1376 && (fstar.w[1] || fstar.w[0])
1377 && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
1378 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1379 && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
1380 // the result is a midpoint; round to nearest
1381 if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
1382 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
1383 Cstar.w[0]--; // Cstar.w[0] is now even
1384 is_inexact_gt_midpoint = 0;
1385 } else { // else MP in [ODD, EVEN]
1386 is_midpoint_lt_even = 1;
1387 is_inexact_gt_midpoint = 0;
1388 }
1389 }
1390 // general correction for RM
1391 if (is_midpoint_lt_even || is_inexact_gt_midpoint) {
1392 Cstar.w[0] = Cstar.w[0] - 1;
1393 } else {
1394 ; // the result is already correct
1395 }
1396 res = Cstar.w[0];
1397 } else if (exp == 0) {
1398 // 1 <= q <= 10
1399 // res = +C (exact)
1400 res = C1.w[0];
1401 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1402 // res = +C * 10^exp (exact)
1403 res = C1.w[0] * ten2k64[exp];
1404 }
1405 }
1406 }
1407
1408 BID_RETURN (res);
1409 }
1410
1411 /*****************************************************************************
1412 * BID128_to_uint32_ceil
1413 ****************************************************************************/
1414
1415 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
1416 bid128_to_uint32_ceil, x)
1417
1418 unsigned int res;
1419 UINT64 x_sign;
1420 UINT64 x_exp;
1421 int exp; // unbiased exponent
1422 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1423 UINT64 tmp64, tmp64A;
1424 BID_UI64DOUBLE tmp1;
1425 unsigned int x_nr_bits;
1426 int q, ind, shift;
1427 UINT128 C1, C;
1428 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
1429 UINT256 fstar;
1430 UINT256 P256;
1431 int is_inexact_lt_midpoint = 0;
1432 int is_midpoint_gt_even = 0;
1433
1434 // unpack x
1435 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1436 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
1437 C1.w[1] = x.w[1] & MASK_COEFF;
1438 C1.w[0] = x.w[0];
1439
1440 // check for NaN or Infinity
1441 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1442 // x is special
1443 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
1444 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
1445 // set invalid flag
1446 *pfpsf |= INVALID_EXCEPTION;
1447 // return Integer Indefinite
1448 res = 0x80000000;
1449 } else { // x is QNaN
1450 // set invalid flag
1451 *pfpsf |= INVALID_EXCEPTION;
1452 // return Integer Indefinite
1453 res = 0x80000000;
1454 }
1455 BID_RETURN (res);
1456 } else { // x is not a NaN, so it must be infinity
1457 if (!x_sign) { // x is +inf
1458 // set invalid flag
1459 *pfpsf |= INVALID_EXCEPTION;
1460 // return Integer Indefinite
1461 res = 0x80000000;
1462 } else { // x is -inf
1463 // set invalid flag
1464 *pfpsf |= INVALID_EXCEPTION;
1465 // return Integer Indefinite
1466 res = 0x80000000;
1467 }
1468 BID_RETURN (res);
1469 }
1470 }
1471 // check for non-canonical values (after the check for special values)
1472 if ((C1.w[1] > 0x0001ed09bead87c0ull)
1473 || (C1.w[1] == 0x0001ed09bead87c0ull
1474 && (C1.w[0] > 0x378d8e63ffffffffull))
1475 || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
1476 res = 0x00000000;
1477 BID_RETURN (res);
1478 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1479 // x is 0
1480 res = 0x00000000;
1481 BID_RETURN (res);
1482 } else { // x is not special and is not zero
1483
1484 // q = nr. of decimal digits in x
1485 // determine first the nr. of bits in x
1486 if (C1.w[1] == 0) {
1487 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
1488 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1489 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
1490 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
1491 x_nr_bits =
1492 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1493 } else { // x < 2^32
1494 tmp1.d = (double) (C1.w[0]); // exact conversion
1495 x_nr_bits =
1496 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1497 }
1498 } else { // if x < 2^53
1499 tmp1.d = (double) C1.w[0]; // exact conversion
1500 x_nr_bits =
1501 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1502 }
1503 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1504 tmp1.d = (double) C1.w[1]; // exact conversion
1505 x_nr_bits =
1506 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1507 }
1508 q = nr_digits[x_nr_bits - 1].digits;
1509 if (q == 0) {
1510 q = nr_digits[x_nr_bits - 1].digits1;
1511 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
1512 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
1513 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1514 q++;
1515 }
1516 exp = (x_exp >> 49) - 6176;
1517
1518 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1519 // set invalid flag
1520 *pfpsf |= INVALID_EXCEPTION;
1521 // return Integer Indefinite
1522 res = 0x80000000;
1523 BID_RETURN (res);
1524 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1525 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1526 // so x rounded to an integer may or may not fit in a signed 32-bit int
1527 // the cases that do not fit are identified here; the ones that fit
1528 // fall through and will be handled with other cases further,
1529 // under '1 <= q + exp <= 10'
1530 if (x_sign) { // if n < 0 and q + exp = 10
1531 // if n <= -1 then n is too large
1532 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
1533 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34
1534 if (q <= 11) {
1535 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
1536 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1537 if (tmp64 >= 0x0aull) {
1538 // set invalid flag
1539 *pfpsf |= INVALID_EXCEPTION;
1540 // return Integer Indefinite
1541 res = 0x80000000;
1542 BID_RETURN (res);
1543 }
1544 // else cases that can be rounded to a 32-bit int fall through
1545 // to '1 <= q + exp <= 10'
1546 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1547 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
1548 // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
1549 // (scale 1 up)
1550 tmp64 = 0x0aull;
1551 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1552 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
1553 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1554 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
1555 }
1556 if (C1.w[1] > C.w[1]
1557 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
1558 // set invalid flag
1559 *pfpsf |= INVALID_EXCEPTION;
1560 // return Integer Indefinite
1561 res = 0x80000000;
1562 BID_RETURN (res);
1563 }
1564 // else cases that can be rounded to a 32-bit int fall through
1565 // to '1 <= q + exp <= 10'
1566 }
1567 } else { // if n > 0 and q + exp = 10
1568 // if n > 2^32 - 1 then n is too large
1569 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^32 - 1
1570 // too large if 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6, 1<=q<=34
1571 if (q <= 11) {
1572 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
1573 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1574 if (tmp64 > 0x9fffffff6ull) {
1575 // set invalid flag
1576 *pfpsf |= INVALID_EXCEPTION;
1577 // return Integer Indefinite
1578 res = 0x80000000;
1579 BID_RETURN (res);
1580 }
1581 // else cases that can be rounded to a 32-bit int fall through
1582 // to '1 <= q + exp <= 10'
1583 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1584 // 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6 <=>
1585 // C > 0x9fffffff6 * 10^(q-11) where 1 <= q - 11 <= 23
1586 // (scale 2^32 up)
1587 tmp64 = 0x9fffffff6ull;
1588 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1589 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
1590 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1591 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
1592 }
1593 if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
1594 // set invalid flag
1595 *pfpsf |= INVALID_EXCEPTION;
1596 // return Integer Indefinite
1597 res = 0x80000000;
1598 BID_RETURN (res);
1599 }
1600 // else cases that can be rounded to a 32-bit int fall through
1601 // to '1 <= q + exp <= 10'
1602 }
1603 }
1604 }
1605 // n is not too large to be converted to int32: -2^32-1 < n <= 2^32-1
1606 // Note: some of the cases tested for above fall through to this point
1607 if ((q + exp) <= 0) {
1608 // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
1609 // return 0
1610 if (x_sign)
1611 res = 0x00000000;
1612 else
1613 res = 0x00000001;
1614 BID_RETURN (res);
1615 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
1616 // -2^32-1 < x <= -1 or 1 <= x <= 2^32-1 so x can be rounded
1617 // toward positive infinity to a 32-bit signed integer
1618 if (x_sign) { // x <= -1 is invalid
1619 // set invalid flag
1620 *pfpsf |= INVALID_EXCEPTION;
1621 // return Integer Indefinite
1622 res = 0x80000000;
1623 BID_RETURN (res);
1624 }
1625 // x > 0 from this point on
1626 if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
1627 ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
1628 // chop off ind digits from the lower part of C1
1629 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
1630 tmp64 = C1.w[0];
1631 if (ind <= 19) {
1632 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1633 } else {
1634 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1635 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1636 }
1637 if (C1.w[0] < tmp64)
1638 C1.w[1]++;
1639 // calculate C* and f*
1640 // C* is actually floor(C*) in this case
1641 // C* and f* need shifting and masking, as shown by
1642 // shiftright128[] and maskhigh128[]
1643 // 1 <= x <= 33
1644 // kx = 10^(-x) = ten2mk128[ind - 1]
1645 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1646 // the approximation of 10^(-x) was rounded up to 118 bits
1647 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1648 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
1649 Cstar.w[1] = P256.w[3];
1650 Cstar.w[0] = P256.w[2];
1651 fstar.w[3] = 0;
1652 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1653 fstar.w[1] = P256.w[1];
1654 fstar.w[0] = P256.w[0];
1655 } else { // 22 <= ind - 1 <= 33
1656 Cstar.w[1] = 0;
1657 Cstar.w[0] = P256.w[3];
1658 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1659 fstar.w[2] = P256.w[2];
1660 fstar.w[1] = P256.w[1];
1661 fstar.w[0] = P256.w[0];
1662 }
1663 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1664 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1665 // if (0 < f* < 10^(-x)) then the result is a midpoint
1666 // if floor(C*) is even then C* = floor(C*) - logical right
1667 // shift; C* has p decimal digits, correct by Prop. 1)
1668 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
1669 // shift; C* has p decimal digits, correct by Pr. 1)
1670 // else
1671 // C* = floor(C*) (logical right shift; C has p decimal digits,
1672 // correct by Property 1)
1673 // n = C* * 10^(e+x)
1674
1675 // shift right C* by Ex-128 = shiftright128[ind]
1676 shift = shiftright128[ind - 1]; // 0 <= shift <= 102
1677 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
1678 Cstar.w[0] =
1679 (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
1680 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1681 } else { // 22 <= ind - 1 <= 33
1682 Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
1683 }
1684 // determine inexactness of the rounding of C*
1685 // if (0 < f* - 1/2 < 10^(-x)) then
1686 // the result is exact
1687 // else // if (f* - 1/2 > T*) then
1688 // the result is inexact
1689 if (ind - 1 <= 2) {
1690 if (fstar.w[1] > 0x8000000000000000ull ||
1691 (fstar.w[1] == 0x8000000000000000ull
1692 && fstar.w[0] > 0x0ull)) {
1693 // f* > 1/2 and the result may be exact
1694 tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
1695 if (tmp64 > ten2mk128trunc[ind - 1].w[1]
1696 || (tmp64 == ten2mk128trunc[ind - 1].w[1]
1697 && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
1698 is_inexact_lt_midpoint = 1;
1699 } // else the result is exact
1700 } else { // the result is inexact; f2* <= 1/2
1701 ;
1702 }
1703 } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
1704 if (fstar.w[3] > 0x0 ||
1705 (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
1706 (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
1707 (fstar.w[1] || fstar.w[0]))) {
1708 // f2* > 1/2 and the result may be exact
1709 // Calculate f2* - 1/2
1710 tmp64 = fstar.w[2] - onehalf128[ind - 1];
1711 tmp64A = fstar.w[3];
1712 if (tmp64 > fstar.w[2])
1713 tmp64A--;
1714 if (tmp64A || tmp64
1715 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1716 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1717 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1718 is_inexact_lt_midpoint = 1;
1719 } // else the result is exact
1720 } else { // the result is inexact; f2* <= 1/2
1721 }
1722 } else { // if 22 <= ind <= 33
1723 if (fstar.w[3] > onehalf128[ind - 1] ||
1724 (fstar.w[3] == onehalf128[ind - 1] &&
1725 (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
1726 // f2* > 1/2 and the result may be exact
1727 // Calculate f2* - 1/2
1728 tmp64 = fstar.w[3] - onehalf128[ind - 1];
1729 if (tmp64 || fstar.w[2]
1730 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1731 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1732 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1733 is_inexact_lt_midpoint = 1;
1734 } // else the result is exact
1735 } else { // the result is inexact; f2* <= 1/2
1736 ;
1737 }
1738 }
1739
1740 // if the result was a midpoint it was rounded away from zero, so
1741 // it will need a correction
1742 // check for midpoints
1743 if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
1744 && (fstar.w[1] || fstar.w[0])
1745 && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
1746 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1747 && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
1748 // the result is a midpoint; round to nearest
1749 if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
1750 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
1751 Cstar.w[0]--; // Cstar.w[0] is now even
1752 is_midpoint_gt_even = 1;
1753 is_inexact_lt_midpoint = 0;
1754 } else { // else MP in [ODD, EVEN]
1755 is_inexact_lt_midpoint = 0;
1756 }
1757 }
1758 // general correction for RM
1759 if (is_midpoint_gt_even || is_inexact_lt_midpoint) {
1760 Cstar.w[0] = Cstar.w[0] + 1;
1761 } else {
1762 ; // the result is already correct
1763 }
1764 res = Cstar.w[0];
1765 } else if (exp == 0) {
1766 // 1 <= q <= 10
1767 // res = +C (exact)
1768 res = C1.w[0];
1769 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1770 // res = +C * 10^exp (exact)
1771 res = C1.w[0] * ten2k64[exp];
1772 }
1773 }
1774 }
1775
1776 BID_RETURN (res);
1777 }
1778
1779 /*****************************************************************************
1780 * BID128_to_uint32_xceil
1781 ****************************************************************************/
1782
1783 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
1784 bid128_to_uint32_xceil, x)
1785
1786 unsigned int res;
1787 UINT64 x_sign;
1788 UINT64 x_exp;
1789 int exp; // unbiased exponent
1790 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1791 UINT64 tmp64, tmp64A;
1792 BID_UI64DOUBLE tmp1;
1793 unsigned int x_nr_bits;
1794 int q, ind, shift;
1795 UINT128 C1, C;
1796 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
1797 UINT256 fstar;
1798 UINT256 P256;
1799 int is_inexact_lt_midpoint = 0;
1800 int is_midpoint_gt_even = 0;
1801
1802 // unpack x
1803 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1804 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
1805 C1.w[1] = x.w[1] & MASK_COEFF;
1806 C1.w[0] = x.w[0];
1807
1808 // check for NaN or Infinity
1809 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1810 // x is special
1811 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
1812 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
1813 // set invalid flag
1814 *pfpsf |= INVALID_EXCEPTION;
1815 // return Integer Indefinite
1816 res = 0x80000000;
1817 } else { // x is QNaN
1818 // set invalid flag
1819 *pfpsf |= INVALID_EXCEPTION;
1820 // return Integer Indefinite
1821 res = 0x80000000;
1822 }
1823 BID_RETURN (res);
1824 } else { // x is not a NaN, so it must be infinity
1825 if (!x_sign) { // x is +inf
1826 // set invalid flag
1827 *pfpsf |= INVALID_EXCEPTION;
1828 // return Integer Indefinite
1829 res = 0x80000000;
1830 } else { // x is -inf
1831 // set invalid flag
1832 *pfpsf |= INVALID_EXCEPTION;
1833 // return Integer Indefinite
1834 res = 0x80000000;
1835 }
1836 BID_RETURN (res);
1837 }
1838 }
1839 // check for non-canonical values (after the check for special values)
1840 if ((C1.w[1] > 0x0001ed09bead87c0ull)
1841 || (C1.w[1] == 0x0001ed09bead87c0ull
1842 && (C1.w[0] > 0x378d8e63ffffffffull))
1843 || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
1844 res = 0x00000000;
1845 BID_RETURN (res);
1846 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1847 // x is 0
1848 res = 0x00000000;
1849 BID_RETURN (res);
1850 } else { // x is not special and is not zero
1851
1852 // q = nr. of decimal digits in x
1853 // determine first the nr. of bits in x
1854 if (C1.w[1] == 0) {
1855 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
1856 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1857 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
1858 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
1859 x_nr_bits =
1860 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1861 } else { // x < 2^32
1862 tmp1.d = (double) (C1.w[0]); // exact conversion
1863 x_nr_bits =
1864 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1865 }
1866 } else { // if x < 2^53
1867 tmp1.d = (double) C1.w[0]; // exact conversion
1868 x_nr_bits =
1869 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1870 }
1871 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1872 tmp1.d = (double) C1.w[1]; // exact conversion
1873 x_nr_bits =
1874 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1875 }
1876 q = nr_digits[x_nr_bits - 1].digits;
1877 if (q == 0) {
1878 q = nr_digits[x_nr_bits - 1].digits1;
1879 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
1880 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
1881 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1882 q++;
1883 }
1884 exp = (x_exp >> 49) - 6176;
1885 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1886 // set invalid flag
1887 *pfpsf |= INVALID_EXCEPTION;
1888 // return Integer Indefinite
1889 res = 0x80000000;
1890 BID_RETURN (res);
1891 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1892 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1893 // so x rounded to an integer may or may not fit in a signed 32-bit int
1894 // the cases that do not fit are identified here; the ones that fit
1895 // fall through and will be handled with other cases further,
1896 // under '1 <= q + exp <= 10'
1897 if (x_sign) { // if n < 0 and q + exp = 10
1898 // if n <= -1 then n is too large
1899 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
1900 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34
1901 if (q <= 11) {
1902 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
1903 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1904 if (tmp64 >= 0x0aull) {
1905 // set invalid flag
1906 *pfpsf |= INVALID_EXCEPTION;
1907 // return Integer Indefinite
1908 res = 0x80000000;
1909 BID_RETURN (res);
1910 }
1911 // else cases that can be rounded to a 32-bit int fall through
1912 // to '1 <= q + exp <= 10'
1913 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1914 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
1915 // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
1916 // (scale 1 up)
1917 tmp64 = 0x0aull;
1918 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1919 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
1920 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1921 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
1922 }
1923 if (C1.w[1] > C.w[1]
1924 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
1925 // set invalid flag
1926 *pfpsf |= INVALID_EXCEPTION;
1927 // return Integer Indefinite
1928 res = 0x80000000;
1929 BID_RETURN (res);
1930 }
1931 // else cases that can be rounded to a 32-bit int fall through
1932 // to '1 <= q + exp <= 10'
1933 }
1934 } else { // if n > 0 and q + exp = 10
1935 // if n > 2^32 - 1 then n is too large
1936 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^32 - 1
1937 // too large if 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6, 1<=q<=34
1938 if (q <= 11) {
1939 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
1940 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1941 if (tmp64 > 0x9fffffff6ull) {
1942 // set invalid flag
1943 *pfpsf |= INVALID_EXCEPTION;
1944 // return Integer Indefinite
1945 res = 0x80000000;
1946 BID_RETURN (res);
1947 }
1948 // else cases that can be rounded to a 32-bit int fall through
1949 // to '1 <= q + exp <= 10'
1950 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1951 // 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6 <=>
1952 // C > 0x9fffffff6 * 10^(q-11) where 1 <= q - 11 <= 23
1953 // (scale 2^32 up)
1954 tmp64 = 0x9fffffff6ull;
1955 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1956 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
1957 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1958 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
1959 }
1960 if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
1961 // set invalid flag
1962 *pfpsf |= INVALID_EXCEPTION;
1963 // return Integer Indefinite
1964 res = 0x80000000;
1965 BID_RETURN (res);
1966 }
1967 // else cases that can be rounded to a 32-bit int fall through
1968 // to '1 <= q + exp <= 10'
1969 }
1970 }
1971 }
1972 // n is not too large to be converted to int32: -2^32-1 < n <= 2^32-1
1973 // Note: some of the cases tested for above fall through to this point
1974 if ((q + exp) <= 0) {
1975 // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
1976 // set inexact flag
1977 *pfpsf |= INEXACT_EXCEPTION;
1978 // return 0
1979 if (x_sign)
1980 res = 0x00000000;
1981 else
1982 res = 0x00000001;
1983 BID_RETURN (res);
1984 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
1985 // -2^32-1 < x <= -1 or 1 <= x <= 2^32-1 so x can be rounded
1986 // toward positive infinity to a 32-bit signed integer
1987 if (x_sign) { // x <= -1 is invalid
1988 // set invalid flag
1989 *pfpsf |= INVALID_EXCEPTION;
1990 // return Integer Indefinite
1991 res = 0x80000000;
1992 BID_RETURN (res);
1993 }
1994 // x > 0 from this point on
1995 if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
1996 ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
1997 // chop off ind digits from the lower part of C1
1998 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
1999 tmp64 = C1.w[0];
2000 if (ind <= 19) {
2001 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
2002 } else {
2003 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
2004 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
2005 }
2006 if (C1.w[0] < tmp64)
2007 C1.w[1]++;
2008 // calculate C* and f*
2009 // C* is actually floor(C*) in this case
2010 // C* and f* need shifting and masking, as shown by
2011 // shiftright128[] and maskhigh128[]
2012 // 1 <= x <= 33
2013 // kx = 10^(-x) = ten2mk128[ind - 1]
2014 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2015 // the approximation of 10^(-x) was rounded up to 118 bits
2016 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
2017 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
2018 Cstar.w[1] = P256.w[3];
2019 Cstar.w[0] = P256.w[2];
2020 fstar.w[3] = 0;
2021 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
2022 fstar.w[1] = P256.w[1];
2023 fstar.w[0] = P256.w[0];
2024 } else { // 22 <= ind - 1 <= 33
2025 Cstar.w[1] = 0;
2026 Cstar.w[0] = P256.w[3];
2027 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
2028 fstar.w[2] = P256.w[2];
2029 fstar.w[1] = P256.w[1];
2030 fstar.w[0] = P256.w[0];
2031 }
2032 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2033 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2034 // if (0 < f* < 10^(-x)) then the result is a midpoint
2035 // if floor(C*) is even then C* = floor(C*) - logical right
2036 // shift; C* has p decimal digits, correct by Prop. 1)
2037 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
2038 // shift; C* has p decimal digits, correct by Pr. 1)
2039 // else
2040 // C* = floor(C*) (logical right shift; C has p decimal digits,
2041 // correct by Property 1)
2042 // n = C* * 10^(e+x)
2043
2044 // shift right C* by Ex-128 = shiftright128[ind]
2045 shift = shiftright128[ind - 1]; // 0 <= shift <= 102
2046 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
2047 Cstar.w[0] =
2048 (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
2049 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2050 } else { // 22 <= ind - 1 <= 33
2051 Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
2052 }
2053 // determine inexactness of the rounding of C*
2054 // if (0 < f* - 1/2 < 10^(-x)) then
2055 // the result is exact
2056 // else // if (f* - 1/2 > T*) then
2057 // the result is inexact
2058 if (ind - 1 <= 2) {
2059 if (fstar.w[1] > 0x8000000000000000ull ||
2060 (fstar.w[1] == 0x8000000000000000ull
2061 && fstar.w[0] > 0x0ull)) {
2062 // f* > 1/2 and the result may be exact
2063 tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
2064 if (tmp64 > ten2mk128trunc[ind - 1].w[1]
2065 || (tmp64 == ten2mk128trunc[ind - 1].w[1]
2066 && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
2067 // set the inexact flag
2068 *pfpsf |= INEXACT_EXCEPTION;
2069 is_inexact_lt_midpoint = 1;
2070 } // else the result is exact
2071 } else { // the result is inexact; f2* <= 1/2
2072 // set the inexact flag
2073 *pfpsf |= INEXACT_EXCEPTION;
2074 }
2075 } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
2076 if (fstar.w[3] > 0x0 ||
2077 (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
2078 (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
2079 (fstar.w[1] || fstar.w[0]))) {
2080 // f2* > 1/2 and the result may be exact
2081 // Calculate f2* - 1/2
2082 tmp64 = fstar.w[2] - onehalf128[ind - 1];
2083 tmp64A = fstar.w[3];
2084 if (tmp64 > fstar.w[2])
2085 tmp64A--;
2086 if (tmp64A || tmp64
2087 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2088 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2089 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2090 // set the inexact flag
2091 *pfpsf |= INEXACT_EXCEPTION;
2092 is_inexact_lt_midpoint = 1;
2093 } // else the result is exact
2094 } else { // the result is inexact; f2* <= 1/2
2095 // set the inexact flag
2096 *pfpsf |= INEXACT_EXCEPTION;
2097 }
2098 } else { // if 22 <= ind <= 33
2099 if (fstar.w[3] > onehalf128[ind - 1] ||
2100 (fstar.w[3] == onehalf128[ind - 1] &&
2101 (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
2102 // f2* > 1/2 and the result may be exact
2103 // Calculate f2* - 1/2
2104 tmp64 = fstar.w[3] - onehalf128[ind - 1];
2105 if (tmp64 || fstar.w[2]
2106 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2107 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2108 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2109 // set the inexact flag
2110 *pfpsf |= INEXACT_EXCEPTION;
2111 is_inexact_lt_midpoint = 1;
2112 } // else the result is exact
2113 } else { // the result is inexact; f2* <= 1/2
2114 // set the inexact flag
2115 *pfpsf |= INEXACT_EXCEPTION;
2116 }
2117 }
2118
2119 // if the result was a midpoint it was rounded away from zero, so
2120 // it will need a correction
2121 // check for midpoints
2122 if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
2123 && (fstar.w[1] || fstar.w[0])
2124 && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
2125 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2126 && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
2127 // the result is a midpoint; round to nearest
2128 if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
2129 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
2130 Cstar.w[0]--; // Cstar.w[0] is now even
2131 is_midpoint_gt_even = 1;
2132 is_inexact_lt_midpoint = 0;
2133 } else { // else MP in [ODD, EVEN]
2134 is_inexact_lt_midpoint = 0;
2135 }
2136 }
2137 // general correction for RM
2138 if (is_midpoint_gt_even || is_inexact_lt_midpoint) {
2139 Cstar.w[0] = Cstar.w[0] + 1;
2140 } else {
2141 ; // the result is already correct
2142 }
2143 res = Cstar.w[0];
2144 } else if (exp == 0) {
2145 // 1 <= q <= 10
2146 // res = +C (exact)
2147 res = C1.w[0];
2148 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2149 // res = +C * 10^exp (exact)
2150 res = C1.w[0] * ten2k64[exp];
2151 }
2152 }
2153 }
2154
2155 BID_RETURN (res);
2156 }
2157
2158 /*****************************************************************************
2159 * BID128_to_uint32_int
2160 ****************************************************************************/
2161
2162 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
2163 bid128_to_uint32_int, x)
2164
2165 int res;
2166 UINT64 x_sign;
2167 UINT64 x_exp;
2168 int exp; // unbiased exponent
2169 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2170 UINT64 tmp64, tmp64A;
2171 BID_UI64DOUBLE tmp1;
2172 unsigned int x_nr_bits;
2173 int q, ind, shift;
2174 UINT128 C1, C;
2175 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
2176 UINT256 fstar;
2177 UINT256 P256;
2178 int is_inexact_gt_midpoint = 0;
2179 int is_midpoint_lt_even = 0;
2180
2181 // unpack x
2182 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
2183 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
2184 C1.w[1] = x.w[1] & MASK_COEFF;
2185 C1.w[0] = x.w[0];
2186
2187 // check for NaN or Infinity
2188 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
2189 // x is special
2190 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
2191 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
2192 // set invalid flag
2193 *pfpsf |= INVALID_EXCEPTION;
2194 // return Integer Indefinite
2195 res = 0x80000000;
2196 } else { // x is QNaN
2197 // set invalid flag
2198 *pfpsf |= INVALID_EXCEPTION;
2199 // return Integer Indefinite
2200 res = 0x80000000;
2201 }
2202 BID_RETURN (res);
2203 } else { // x is not a NaN, so it must be infinity
2204 if (!x_sign) { // x is +inf
2205 // set invalid flag
2206 *pfpsf |= INVALID_EXCEPTION;
2207 // return Integer Indefinite
2208 res = 0x80000000;
2209 } else { // x is -inf
2210 // set invalid flag
2211 *pfpsf |= INVALID_EXCEPTION;
2212 // return Integer Indefinite
2213 res = 0x80000000;
2214 }
2215 BID_RETURN (res);
2216 }
2217 }
2218 // check for non-canonical values (after the check for special values)
2219 if ((C1.w[1] > 0x0001ed09bead87c0ull)
2220 || (C1.w[1] == 0x0001ed09bead87c0ull
2221 && (C1.w[0] > 0x378d8e63ffffffffull))
2222 || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
2223 res = 0x00000000;
2224 BID_RETURN (res);
2225 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
2226 // x is 0
2227 res = 0x00000000;
2228 BID_RETURN (res);
2229 } else { // x is not special and is not zero
2230
2231 // q = nr. of decimal digits in x
2232 // determine first the nr. of bits in x
2233 if (C1.w[1] == 0) {
2234 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
2235 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2236 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
2237 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
2238 x_nr_bits =
2239 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2240 } else { // x < 2^32
2241 tmp1.d = (double) (C1.w[0]); // exact conversion
2242 x_nr_bits =
2243 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2244 }
2245 } else { // if x < 2^53
2246 tmp1.d = (double) C1.w[0]; // exact conversion
2247 x_nr_bits =
2248 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2249 }
2250 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2251 tmp1.d = (double) C1.w[1]; // exact conversion
2252 x_nr_bits =
2253 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2254 }
2255 q = nr_digits[x_nr_bits - 1].digits;
2256 if (q == 0) {
2257 q = nr_digits[x_nr_bits - 1].digits1;
2258 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
2259 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
2260 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
2261 q++;
2262 }
2263 exp = (x_exp >> 49) - 6176;
2264 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
2265 // set invalid flag
2266 *pfpsf |= INVALID_EXCEPTION;
2267 // return Integer Indefinite
2268 res = 0x80000000;
2269 BID_RETURN (res);
2270 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
2271 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
2272 // so x rounded to an integer may or may not fit in a signed 32-bit int
2273 // the cases that do not fit are identified here; the ones that fit
2274 // fall through and will be handled with other cases further,
2275 // under '1 <= q + exp <= 10'
2276 if (x_sign) { // if n < 0 and q + exp = 10
2277 // if n <= -1 then n is too large
2278 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
2279 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a, 1<=q<=34
2280 if (q <= 11) {
2281 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
2282 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2283 if (tmp64 >= 0x0aull) {
2284 // set invalid flag
2285 *pfpsf |= INVALID_EXCEPTION;
2286 // return Integer Indefinite
2287 res = 0x80000000;
2288 BID_RETURN (res);
2289 }
2290 // else cases that can be rounded to a 32-bit uint fall through
2291 // to '1 <= q + exp <= 10'
2292 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2293 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
2294 // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
2295 // (scale 1 up)
2296 tmp64 = 0x0aull;
2297 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2298 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
2299 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2300 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
2301 }
2302 if (C1.w[1] > C.w[1]
2303 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2304 // set invalid flag
2305 *pfpsf |= INVALID_EXCEPTION;
2306 // return Integer Indefinite
2307 res = 0x80000000;
2308 BID_RETURN (res);
2309 }
2310 // else cases that can be rounded to a 32-bit int fall through
2311 // to '1 <= q + exp <= 10'
2312 }
2313 } else { // if n > 0 and q + exp = 10
2314 // if n >= 2^32 then n is too large
2315 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
2316 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
2317 if (q <= 11) {
2318 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
2319 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2320 if (tmp64 >= 0xa00000000ull) {
2321 // set invalid flag
2322 *pfpsf |= INVALID_EXCEPTION;
2323 // return Integer Indefinite
2324 res = 0x80000000;
2325 BID_RETURN (res);
2326 }
2327 // else cases that can be rounded to a 32-bit uint fall through
2328 // to '1 <= q + exp <= 10'
2329 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2330 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
2331 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
2332 // (scale 2^32 up)
2333 tmp64 = 0xa00000000ull;
2334 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2335 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
2336 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2337 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
2338 }
2339 if (C1.w[1] > C.w[1]
2340 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2341 // set invalid flag
2342 *pfpsf |= INVALID_EXCEPTION;
2343 // return Integer Indefinite
2344 res = 0x80000000;
2345 BID_RETURN (res);
2346 }
2347 // else cases that can be rounded to a 32-bit int fall through
2348 // to '1 <= q + exp <= 10'
2349 }
2350 }
2351 }
2352 // n is not too large to be converted to uint32: -2^32 < n < 2^32
2353 // Note: some of the cases tested for above fall through to this point
2354 if ((q + exp) <= 0) {
2355 // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
2356 // return 0
2357 res = 0x00000000;
2358 BID_RETURN (res);
2359 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
2360 // x = d(0)...d(k).d(k+1)..., k >= 0, d(0) != 0
2361 if (x_sign) { // x <= -1
2362 // set invalid flag
2363 *pfpsf |= INVALID_EXCEPTION;
2364 // return Integer Indefinite
2365 res = 0x80000000;
2366 BID_RETURN (res);
2367 }
2368 // x > 0 from this point on
2369 // 1 <= x < 2^32 so x can be rounded to zero to a 32-bit unsigned integer
2370 if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
2371 ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
2372 // chop off ind digits from the lower part of C1
2373 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
2374 tmp64 = C1.w[0];
2375 if (ind <= 19) {
2376 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
2377 } else {
2378 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
2379 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
2380 }
2381 if (C1.w[0] < tmp64)
2382 C1.w[1]++;
2383 // calculate C* and f*
2384 // C* is actually floor(C*) in this case
2385 // C* and f* need shifting and masking, as shown by
2386 // shiftright128[] and maskhigh128[]
2387 // 1 <= x <= 33
2388 // kx = 10^(-x) = ten2mk128[ind - 1]
2389 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2390 // the approximation of 10^(-x) was rounded up to 118 bits
2391 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
2392 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
2393 Cstar.w[1] = P256.w[3];
2394 Cstar.w[0] = P256.w[2];
2395 fstar.w[3] = 0;
2396 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
2397 fstar.w[1] = P256.w[1];
2398 fstar.w[0] = P256.w[0];
2399 } else { // 22 <= ind - 1 <= 33
2400 Cstar.w[1] = 0;
2401 Cstar.w[0] = P256.w[3];
2402 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
2403 fstar.w[2] = P256.w[2];
2404 fstar.w[1] = P256.w[1];
2405 fstar.w[0] = P256.w[0];
2406 }
2407 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2408 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2409 // if (0 < f* < 10^(-x)) then the result is a midpoint
2410 // if floor(C*) is even then C* = floor(C*) - logical right
2411 // shift; C* has p decimal digits, correct by Prop. 1)
2412 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
2413 // shift; C* has p decimal digits, correct by Pr. 1)
2414 // else
2415 // C* = floor(C*) (logical right shift; C has p decimal digits,
2416 // correct by Property 1)
2417 // n = C* * 10^(e+x)
2418
2419 // shift right C* by Ex-128 = shiftright128[ind]
2420 shift = shiftright128[ind - 1]; // 0 <= shift <= 102
2421 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
2422 Cstar.w[0] =
2423 (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
2424 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2425 } else { // 22 <= ind - 1 <= 33
2426 Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
2427 }
2428 // determine inexactness of the rounding of C*
2429 // if (0 < f* - 1/2 < 10^(-x)) then
2430 // the result is exact
2431 // else // if (f* - 1/2 > T*) then
2432 // the result is inexact
2433 if (ind - 1 <= 2) {
2434 if (fstar.w[1] > 0x8000000000000000ull ||
2435 (fstar.w[1] == 0x8000000000000000ull
2436 && fstar.w[0] > 0x0ull)) {
2437 // f* > 1/2 and the result may be exact
2438 tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
2439 if (tmp64 > ten2mk128trunc[ind - 1].w[1]
2440 || (tmp64 == ten2mk128trunc[ind - 1].w[1]
2441 && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
2442 } // else the result is exact
2443 } else { // the result is inexact; f2* <= 1/2
2444 is_inexact_gt_midpoint = 1;
2445 }
2446 } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
2447 if (fstar.w[3] > 0x0 ||
2448 (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
2449 (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
2450 (fstar.w[1] || fstar.w[0]))) {
2451 // f2* > 1/2 and the result may be exact
2452 // Calculate f2* - 1/2
2453 tmp64 = fstar.w[2] - onehalf128[ind - 1];
2454 tmp64A = fstar.w[3];
2455 if (tmp64 > fstar.w[2])
2456 tmp64A--;
2457 if (tmp64A || tmp64
2458 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2459 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2460 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2461 } // else the result is exact
2462 } else { // the result is inexact; f2* <= 1/2
2463 is_inexact_gt_midpoint = 1;
2464 }
2465 } else { // if 22 <= ind <= 33
2466 if (fstar.w[3] > onehalf128[ind - 1] ||
2467 (fstar.w[3] == onehalf128[ind - 1] &&
2468 (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
2469 // f2* > 1/2 and the result may be exact
2470 // Calculate f2* - 1/2
2471 tmp64 = fstar.w[3] - onehalf128[ind - 1];
2472 if (tmp64 || fstar.w[2]
2473 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2474 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2475 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2476 } // else the result is exact
2477 } else { // the result is inexact; f2* <= 1/2
2478 is_inexact_gt_midpoint = 1;
2479 }
2480 }
2481
2482 // if the result was a midpoint it was rounded away from zero, so
2483 // it will need a correction
2484 // check for midpoints
2485 if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
2486 && (fstar.w[1] || fstar.w[0])
2487 && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
2488 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2489 && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
2490 // the result is a midpoint; round to nearest
2491 if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
2492 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
2493 Cstar.w[0]--; // Cstar.w[0] is now even
2494 is_inexact_gt_midpoint = 0;
2495 } else { // else MP in [ODD, EVEN]
2496 is_midpoint_lt_even = 1;
2497 is_inexact_gt_midpoint = 0;
2498 }
2499 }
2500 // general correction for RZ
2501 if (is_midpoint_lt_even || is_inexact_gt_midpoint) {
2502 Cstar.w[0] = Cstar.w[0] - 1;
2503 } else {
2504 ; // exact, the result is already correct
2505 }
2506 res = Cstar.w[0];
2507 } else if (exp == 0) {
2508 // 1 <= q <= 10
2509 // res = +C (exact)
2510 res = C1.w[0];
2511 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2512 // res = +C * 10^exp (exact)
2513 res = C1.w[0] * ten2k64[exp];
2514 }
2515 }
2516 }
2517
2518 BID_RETURN (res);
2519 }
2520
2521 /*****************************************************************************
2522 * BID128_to_uint32_xint
2523 ****************************************************************************/
2524
2525 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
2526 bid128_to_uint32_xint, x)
2527
2528 int res;
2529 UINT64 x_sign;
2530 UINT64 x_exp;
2531 int exp; // unbiased exponent
2532 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2533 UINT64 tmp64, tmp64A;
2534 BID_UI64DOUBLE tmp1;
2535 unsigned int x_nr_bits;
2536 int q, ind, shift;
2537 UINT128 C1, C;
2538 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
2539 UINT256 fstar;
2540 UINT256 P256;
2541 int is_inexact_gt_midpoint = 0;
2542 int is_midpoint_lt_even = 0;
2543
2544 // unpack x
2545 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
2546 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
2547 C1.w[1] = x.w[1] & MASK_COEFF;
2548 C1.w[0] = x.w[0];
2549
2550 // check for NaN or Infinity
2551 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
2552 // x is special
2553 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
2554 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
2555 // set invalid flag
2556 *pfpsf |= INVALID_EXCEPTION;
2557 // return Integer Indefinite
2558 res = 0x80000000;
2559 } else { // x is QNaN
2560 // set invalid flag
2561 *pfpsf |= INVALID_EXCEPTION;
2562 // return Integer Indefinite
2563 res = 0x80000000;
2564 }
2565 BID_RETURN (res);
2566 } else { // x is not a NaN, so it must be infinity
2567 if (!x_sign) { // x is +inf
2568 // set invalid flag
2569 *pfpsf |= INVALID_EXCEPTION;
2570 // return Integer Indefinite
2571 res = 0x80000000;
2572 } else { // x is -inf
2573 // set invalid flag
2574 *pfpsf |= INVALID_EXCEPTION;
2575 // return Integer Indefinite
2576 res = 0x80000000;
2577 }
2578 BID_RETURN (res);
2579 }
2580 }
2581 // check for non-canonical values (after the check for special values)
2582 if ((C1.w[1] > 0x0001ed09bead87c0ull)
2583 || (C1.w[1] == 0x0001ed09bead87c0ull
2584 && (C1.w[0] > 0x378d8e63ffffffffull))
2585 || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
2586 res = 0x00000000;
2587 BID_RETURN (res);
2588 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
2589 // x is 0
2590 res = 0x00000000;
2591 BID_RETURN (res);
2592 } else { // x is not special and is not zero
2593
2594 // q = nr. of decimal digits in x
2595 // determine first the nr. of bits in x
2596 if (C1.w[1] == 0) {
2597 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
2598 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2599 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
2600 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
2601 x_nr_bits =
2602 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2603 } else { // x < 2^32
2604 tmp1.d = (double) (C1.w[0]); // exact conversion
2605 x_nr_bits =
2606 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2607 }
2608 } else { // if x < 2^53
2609 tmp1.d = (double) C1.w[0]; // exact conversion
2610 x_nr_bits =
2611 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2612 }
2613 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2614 tmp1.d = (double) C1.w[1]; // exact conversion
2615 x_nr_bits =
2616 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2617 }
2618 q = nr_digits[x_nr_bits - 1].digits;
2619 if (q == 0) {
2620 q = nr_digits[x_nr_bits - 1].digits1;
2621 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
2622 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
2623 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
2624 q++;
2625 }
2626 exp = (x_exp >> 49) - 6176;
2627 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
2628 // set invalid flag
2629 *pfpsf |= INVALID_EXCEPTION;
2630 // return Integer Indefinite
2631 res = 0x80000000;
2632 BID_RETURN (res);
2633 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
2634 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
2635 // so x rounded to an integer may or may not fit in a signed 32-bit int
2636 // the cases that do not fit are identified here; the ones that fit
2637 // fall through and will be handled with other cases further,
2638 // under '1 <= q + exp <= 10'
2639 if (x_sign) { // if n < 0 and q + exp = 10
2640 // if n <= -1 then n is too large
2641 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
2642 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a, 1<=q<=34
2643 if (q <= 11) {
2644 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
2645 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2646 if (tmp64 >= 0x0aull) {
2647 // set invalid flag
2648 *pfpsf |= INVALID_EXCEPTION;
2649 // return Integer Indefinite
2650 res = 0x80000000;
2651 BID_RETURN (res);
2652 }
2653 // else cases that can be rounded to a 32-bit uint fall through
2654 // to '1 <= q + exp <= 10'
2655 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2656 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
2657 // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
2658 // (scale 1 up)
2659 tmp64 = 0x0aull;
2660 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2661 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
2662 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2663 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
2664 }
2665 if (C1.w[1] > C.w[1]
2666 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2667 // set invalid flag
2668 *pfpsf |= INVALID_EXCEPTION;
2669 // return Integer Indefinite
2670 res = 0x80000000;
2671 BID_RETURN (res);
2672 }
2673 // else cases that can be rounded to a 32-bit int fall through
2674 // to '1 <= q + exp <= 10'
2675 }
2676 } else { // if n > 0 and q + exp = 10
2677 // if n >= 2^32 then n is too large
2678 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
2679 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
2680 if (q <= 11) {
2681 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
2682 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2683 if (tmp64 >= 0xa00000000ull) {
2684 // set invalid flag
2685 *pfpsf |= INVALID_EXCEPTION;
2686 // return Integer Indefinite
2687 res = 0x80000000;
2688 BID_RETURN (res);
2689 }
2690 // else cases that can be rounded to a 32-bit uint fall through
2691 // to '1 <= q + exp <= 10'
2692 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2693 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
2694 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
2695 // (scale 2^32 up)
2696 tmp64 = 0xa00000000ull;
2697 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2698 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
2699 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2700 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
2701 }
2702 if (C1.w[1] > C.w[1]
2703 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2704 // set invalid flag
2705 *pfpsf |= INVALID_EXCEPTION;
2706 // return Integer Indefinite
2707 res = 0x80000000;
2708 BID_RETURN (res);
2709 }
2710 // else cases that can be rounded to a 32-bit int fall through
2711 // to '1 <= q + exp <= 10'
2712 }
2713 }
2714 }
2715 // n is not too large to be converted to uint32: -2^32 < n < 2^32
2716 // Note: some of the cases tested for above fall through to this point
2717 if ((q + exp) <= 0) {
2718 // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
2719 // set inexact flag
2720 *pfpsf |= INEXACT_EXCEPTION;
2721 // return 0
2722 res = 0x00000000;
2723 BID_RETURN (res);
2724 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
2725 // x = d(0)...d(k).d(k+1)..., k >= 0, d(0) != 0
2726 if (x_sign) { // x <= -1
2727 // set invalid flag
2728 *pfpsf |= INVALID_EXCEPTION;
2729 // return Integer Indefinite
2730 res = 0x80000000;
2731 BID_RETURN (res);
2732 }
2733 // x > 0 from this point on
2734 // 1 <= x < 2^32 so x can be rounded to zero to a 32-bit unsigned integer
2735 if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
2736 ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
2737 // chop off ind digits from the lower part of C1
2738 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
2739 tmp64 = C1.w[0];
2740 if (ind <= 19) {
2741 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
2742 } else {
2743 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
2744 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
2745 }
2746 if (C1.w[0] < tmp64)
2747 C1.w[1]++;
2748 // calculate C* and f*
2749 // C* is actually floor(C*) in this case
2750 // C* and f* need shifting and masking, as shown by
2751 // shiftright128[] and maskhigh128[]
2752 // 1 <= x <= 33
2753 // kx = 10^(-x) = ten2mk128[ind - 1]
2754 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2755 // the approximation of 10^(-x) was rounded up to 118 bits
2756 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
2757 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
2758 Cstar.w[1] = P256.w[3];
2759 Cstar.w[0] = P256.w[2];
2760 fstar.w[3] = 0;
2761 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
2762 fstar.w[1] = P256.w[1];
2763 fstar.w[0] = P256.w[0];
2764 } else { // 22 <= ind - 1 <= 33
2765 Cstar.w[1] = 0;
2766 Cstar.w[0] = P256.w[3];
2767 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
2768 fstar.w[2] = P256.w[2];
2769 fstar.w[1] = P256.w[1];
2770 fstar.w[0] = P256.w[0];
2771 }
2772 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2773 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2774 // if (0 < f* < 10^(-x)) then the result is a midpoint
2775 // if floor(C*) is even then C* = floor(C*) - logical right
2776 // shift; C* has p decimal digits, correct by Prop. 1)
2777 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
2778 // shift; C* has p decimal digits, correct by Pr. 1)
2779 // else
2780 // C* = floor(C*) (logical right shift; C has p decimal digits,
2781 // correct by Property 1)
2782 // n = C* * 10^(e+x)
2783
2784 // shift right C* by Ex-128 = shiftright128[ind]
2785 shift = shiftright128[ind - 1]; // 0 <= shift <= 102
2786 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
2787 Cstar.w[0] =
2788 (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
2789 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2790 } else { // 22 <= ind - 1 <= 33
2791 Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
2792 }
2793 // determine inexactness of the rounding of C*
2794 // if (0 < f* - 1/2 < 10^(-x)) then
2795 // the result is exact
2796 // else // if (f* - 1/2 > T*) then
2797 // the result is inexact
2798 if (ind - 1 <= 2) {
2799 if (fstar.w[1] > 0x8000000000000000ull ||
2800 (fstar.w[1] == 0x8000000000000000ull
2801 && fstar.w[0] > 0x0ull)) {
2802 // f* > 1/2 and the result may be exact
2803 tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
2804 if (tmp64 > ten2mk128trunc[ind - 1].w[1]
2805 || (tmp64 == ten2mk128trunc[ind - 1].w[1]
2806 && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
2807 // set the inexact flag
2808 *pfpsf |= INEXACT_EXCEPTION;
2809 } // else the result is exact
2810 } else { // the result is inexact; f2* <= 1/2
2811 // set the inexact flag
2812 *pfpsf |= INEXACT_EXCEPTION;
2813 is_inexact_gt_midpoint = 1;
2814 }
2815 } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
2816 if (fstar.w[3] > 0x0 ||
2817 (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
2818 (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
2819 (fstar.w[1] || fstar.w[0]))) {
2820 // f2* > 1/2 and the result may be exact
2821 // Calculate f2* - 1/2
2822 tmp64 = fstar.w[2] - onehalf128[ind - 1];
2823 tmp64A = fstar.w[3];
2824 if (tmp64 > fstar.w[2])
2825 tmp64A--;
2826 if (tmp64A || tmp64
2827 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2828 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2829 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2830 // set the inexact flag
2831 *pfpsf |= INEXACT_EXCEPTION;
2832 } // else the result is exact
2833 } else { // the result is inexact; f2* <= 1/2
2834 // set the inexact flag
2835 *pfpsf |= INEXACT_EXCEPTION;
2836 is_inexact_gt_midpoint = 1;
2837 }
2838 } else { // if 22 <= ind <= 33
2839 if (fstar.w[3] > onehalf128[ind - 1] ||
2840 (fstar.w[3] == onehalf128[ind - 1] &&
2841 (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
2842 // f2* > 1/2 and the result may be exact
2843 // Calculate f2* - 1/2
2844 tmp64 = fstar.w[3] - onehalf128[ind - 1];
2845 if (tmp64 || fstar.w[2]
2846 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2847 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2848 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2849 // set the inexact flag
2850 *pfpsf |= INEXACT_EXCEPTION;
2851 } // else the result is exact
2852 } else { // the result is inexact; f2* <= 1/2
2853 // set the inexact flag
2854 *pfpsf |= INEXACT_EXCEPTION;
2855 is_inexact_gt_midpoint = 1;
2856 }
2857 }
2858
2859 // if the result was a midpoint it was rounded away from zero, so
2860 // it will need a correction
2861 // check for midpoints
2862 if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
2863 && (fstar.w[1] || fstar.w[0])
2864 && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
2865 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2866 && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
2867 // the result is a midpoint; round to nearest
2868 if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
2869 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
2870 Cstar.w[0]--; // Cstar.w[0] is now even
2871 is_inexact_gt_midpoint = 0;
2872 } else { // else MP in [ODD, EVEN]
2873 is_midpoint_lt_even = 1;
2874 is_inexact_gt_midpoint = 0;
2875 }
2876 }
2877 // general correction for RZ
2878 if (is_midpoint_lt_even || is_inexact_gt_midpoint) {
2879 Cstar.w[0] = Cstar.w[0] - 1;
2880 } else {
2881 ; // exact, the result is already correct
2882 }
2883 res = Cstar.w[0];
2884 } else if (exp == 0) {
2885 // 1 <= q <= 10
2886 // res = +C (exact)
2887 res = C1.w[0];
2888 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2889 // res = +C * 10^exp (exact)
2890 res = C1.w[0] * ten2k64[exp];
2891 }
2892 }
2893 }
2894
2895 BID_RETURN (res);
2896 }
2897
2898 /*****************************************************************************
2899 * BID128_to_uint32_rninta
2900 ****************************************************************************/
2901
2902 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
2903 bid128_to_uint32_rninta, x)
2904
2905 unsigned int res;
2906 UINT64 x_sign;
2907 UINT64 x_exp;
2908 int exp; // unbiased exponent
2909 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2910 UINT64 tmp64;
2911 BID_UI64DOUBLE tmp1;
2912 unsigned int x_nr_bits;
2913 int q, ind, shift;
2914 UINT128 C1, C;
2915 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
2916 UINT256 P256;
2917
2918 // unpack x
2919 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
2920 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
2921 C1.w[1] = x.w[1] & MASK_COEFF;
2922 C1.w[0] = x.w[0];
2923
2924 // check for NaN or Infinity
2925 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
2926 // x is special
2927 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
2928 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
2929 // set invalid flag
2930 *pfpsf |= INVALID_EXCEPTION;
2931 // return Integer Indefinite
2932 res = 0x80000000;
2933 } else { // x is QNaN
2934 // set invalid flag
2935 *pfpsf |= INVALID_EXCEPTION;
2936 // return Integer Indefinite
2937 res = 0x80000000;
2938 }
2939 BID_RETURN (res);
2940 } else { // x is not a NaN, so it must be infinity
2941 if (!x_sign) { // x is +inf
2942 // set invalid flag
2943 *pfpsf |= INVALID_EXCEPTION;
2944 // return Integer Indefinite
2945 res = 0x80000000;
2946 } else { // x is -inf
2947 // set invalid flag
2948 *pfpsf |= INVALID_EXCEPTION;
2949 // return Integer Indefinite
2950 res = 0x80000000;
2951 }
2952 BID_RETURN (res);
2953 }
2954 }
2955 // check for non-canonical values (after the check for special values)
2956 if ((C1.w[1] > 0x0001ed09bead87c0ull)
2957 || (C1.w[1] == 0x0001ed09bead87c0ull
2958 && (C1.w[0] > 0x378d8e63ffffffffull))
2959 || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
2960 res = 0x00000000;
2961 BID_RETURN (res);
2962 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
2963 // x is 0
2964 res = 0x00000000;
2965 BID_RETURN (res);
2966 } else { // x is not special and is not zero
2967
2968 // q = nr. of decimal digits in x
2969 // determine first the nr. of bits in x
2970 if (C1.w[1] == 0) {
2971 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
2972 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2973 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
2974 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
2975 x_nr_bits =
2976 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2977 } else { // x < 2^32
2978 tmp1.d = (double) (C1.w[0]); // exact conversion
2979 x_nr_bits =
2980 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2981 }
2982 } else { // if x < 2^53
2983 tmp1.d = (double) C1.w[0]; // exact conversion
2984 x_nr_bits =
2985 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2986 }
2987 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2988 tmp1.d = (double) C1.w[1]; // exact conversion
2989 x_nr_bits =
2990 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2991 }
2992 q = nr_digits[x_nr_bits - 1].digits;
2993 if (q == 0) {
2994 q = nr_digits[x_nr_bits - 1].digits1;
2995 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
2996 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
2997 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
2998 q++;
2999 }
3000 exp = (x_exp >> 49) - 6176;
3001 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
3002 // set invalid flag
3003 *pfpsf |= INVALID_EXCEPTION;
3004 // return Integer Indefinite
3005 res = 0x80000000;
3006 BID_RETURN (res);
3007 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
3008 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
3009 // so x rounded to an integer may or may not fit in a signed 32-bit int
3010 // the cases that do not fit are identified here; the ones that fit
3011 // fall through and will be handled with other cases further,
3012 // under '1 <= q + exp <= 10'
3013 if (x_sign) { // if n < 0 and q + exp = 10
3014 // if n <= -1/2 then n is too large
3015 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1/2
3016 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05, 1<=q<=34
3017 if (q <= 11) {
3018 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
3019 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3020 if (tmp64 >= 0x05ull) {
3021 // set invalid flag
3022 *pfpsf |= INVALID_EXCEPTION;
3023 // return Integer Indefinite
3024 res = 0x80000000;
3025 BID_RETURN (res);
3026 }
3027 // else cases that can be rounded to a 32-bit int fall through
3028 // to '1 <= q + exp <= 10'
3029 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3030 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05 <=>
3031 // C >= 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
3032 // (scale 1/2 up)
3033 tmp64 = 0x05ull;
3034 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3035 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
3036 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3037 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
3038 }
3039 if (C1.w[1] > C.w[1]
3040 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
3041 // set invalid flag
3042 *pfpsf |= INVALID_EXCEPTION;
3043 // return Integer Indefinite
3044 res = 0x80000000;
3045 BID_RETURN (res);
3046 }
3047 // else cases that can be rounded to a 32-bit int fall through
3048 // to '1 <= q + exp <= 10'
3049 }
3050 } else { // if n > 0 and q + exp = 10
3051 // if n >= 2^32 - 1/2 then n is too large
3052 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
3053 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
3054 if (q <= 11) {
3055 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
3056 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3057 if (tmp64 >= 0x9fffffffbull) {
3058 // set invalid flag
3059 *pfpsf |= INVALID_EXCEPTION;
3060 // return Integer Indefinite
3061 res = 0x80000000;
3062 BID_RETURN (res);
3063 }
3064 // else cases that can be rounded to a 32-bit int fall through
3065 // to '1 <= q + exp <= 10'
3066 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3067 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
3068 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
3069 // (scale 2^32-1/2 up)
3070 tmp64 = 0x9fffffffbull;
3071 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3072 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
3073 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3074 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
3075 }
3076 if (C1.w[1] > C.w[1]
3077 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
3078 // set invalid flag
3079 *pfpsf |= INVALID_EXCEPTION;
3080 // return Integer Indefinite
3081 res = 0x80000000;
3082 BID_RETURN (res);
3083 }
3084 // else cases that can be rounded to a 32-bit int fall through
3085 // to '1 <= q + exp <= 10'
3086 }
3087 }
3088 }
3089 // n is not too large to be converted to int32: -1/2 < n < 2^32 - 1/2
3090 // Note: some of the cases tested for above fall through to this point
3091 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
3092 // return 0
3093 res = 0x00000000;
3094 BID_RETURN (res);
3095 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
3096 // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
3097 // res = 0
3098 // else if x > 0
3099 // res = +1
3100 // else // if x < 0
3101 // invalid exc
3102 ind = q - 1;
3103 if (ind <= 18) { // 0 <= ind <= 18
3104 if ((C1.w[1] == 0) && (C1.w[0] < midpoint64[ind])) {
3105 res = 0x00000000; // return 0
3106 } else if (!x_sign) { // n > 0
3107 res = 0x00000001; // return +1
3108 } else {
3109 res = 0x80000000;
3110 *pfpsf |= INVALID_EXCEPTION;
3111 }
3112 } else { // 19 <= ind <= 33
3113 if ((C1.w[1] < midpoint128[ind - 19].w[1])
3114 || ((C1.w[1] == midpoint128[ind - 19].w[1])
3115 && (C1.w[0] < midpoint128[ind - 19].w[0]))) {
3116 res = 0x00000000; // return 0
3117 } else if (!x_sign) { // n > 0
3118 res = 0x00000001; // return +1
3119 } else {
3120 res = 0x80000000;
3121 *pfpsf |= INVALID_EXCEPTION;
3122 }
3123 }
3124 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
3125 if (x_sign) { // x <= -1
3126 // set invalid flag
3127 *pfpsf |= INVALID_EXCEPTION;
3128 // return Integer Indefinite
3129 res = 0x80000000;
3130 BID_RETURN (res);
3131 }
3132 // 1 <= x < 2^31-1/2 so x can be rounded
3133 // to nearest-away to a 32-bit signed integer
3134 if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
3135 ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
3136 // chop off ind digits from the lower part of C1
3137 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
3138 tmp64 = C1.w[0];
3139 if (ind <= 19) {
3140 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
3141 } else {
3142 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
3143 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
3144 }
3145 if (C1.w[0] < tmp64)
3146 C1.w[1]++;
3147 // calculate C* and f*
3148 // C* is actually floor(C*) in this case
3149 // C* and f* need shifting and masking, as shown by
3150 // shiftright128[] and maskhigh128[]
3151 // 1 <= x <= 33
3152 // kx = 10^(-x) = ten2mk128[ind - 1]
3153 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
3154 // the approximation of 10^(-x) was rounded up to 118 bits
3155 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
3156 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
3157 Cstar.w[1] = P256.w[3];
3158 Cstar.w[0] = P256.w[2];
3159 } else { // 22 <= ind - 1 <= 33
3160 Cstar.w[1] = 0;
3161 Cstar.w[0] = P256.w[3];
3162 }
3163 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
3164 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
3165 // if (0 < f* < 10^(-x)) then the result is a midpoint
3166 // if floor(C*) is even then C* = floor(C*) - logical right
3167 // shift; C* has p decimal digits, correct by Prop. 1)
3168 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
3169 // shift; C* has p decimal digits, correct by Pr. 1)
3170 // else
3171 // C* = floor(C*) (logical right shift; C has p decimal digits,
3172 // correct by Property 1)
3173 // n = C* * 10^(e+x)
3174
3175 // shift right C* by Ex-128 = shiftright128[ind]
3176 shift = shiftright128[ind - 1]; // 0 <= shift <= 102
3177 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
3178 Cstar.w[0] =
3179 (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
3180 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
3181 } else { // 22 <= ind - 1 <= 33
3182 Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
3183 }
3184 // if the result was a midpoint, it was already rounded away from zero
3185 res = Cstar.w[0]; // always positive
3186 // no need to check for midpoints - already rounded away from zero!
3187 } else if (exp == 0) {
3188 // 1 <= q <= 10
3189 // res = +C (exact)
3190 res = C1.w[0];
3191 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
3192 // res = +C * 10^exp (exact)
3193 res = C1.w[0] * ten2k64[exp];
3194 }
3195 }
3196 }
3197
3198 BID_RETURN (res);
3199 }
3200
3201 /*****************************************************************************
3202 * BID128_to_uint32_xrninta
3203 ****************************************************************************/
3204
3205 BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
3206 bid128_to_uint32_xrninta, x)
3207
3208 unsigned int res;
3209 UINT64 x_sign;
3210 UINT64 x_exp;
3211 int exp; // unbiased exponent
3212 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
3213 UINT64 tmp64, tmp64A;
3214 BID_UI64DOUBLE tmp1;
3215 unsigned int x_nr_bits;
3216 int q, ind, shift;
3217 UINT128 C1, C;
3218 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
3219 UINT256 fstar;
3220 UINT256 P256;
3221 unsigned int tmp_inexact = 0;
3222
3223 // unpack x
3224 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
3225 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
3226 C1.w[1] = x.w[1] & MASK_COEFF;
3227 C1.w[0] = x.w[0];
3228
3229 // check for NaN or Infinity
3230 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
3231 // x is special
3232 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
3233 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
3234 // set invalid flag
3235 *pfpsf |= INVALID_EXCEPTION;
3236 // return Integer Indefinite
3237 res = 0x80000000;
3238 } else { // x is QNaN
3239 // set invalid flag
3240 *pfpsf |= INVALID_EXCEPTION;
3241 // return Integer Indefinite
3242 res = 0x80000000;
3243 }
3244 BID_RETURN (res);
3245 } else { // x is not a NaN, so it must be infinity
3246 if (!x_sign) { // x is +inf
3247 // set invalid flag
3248 *pfpsf |= INVALID_EXCEPTION;
3249 // return Integer Indefinite
3250 res = 0x80000000;
3251 } else { // x is -inf
3252 // set invalid flag
3253 *pfpsf |= INVALID_EXCEPTION;
3254 // return Integer Indefinite
3255 res = 0x80000000;
3256 }
3257 BID_RETURN (res);
3258 }
3259 }
3260 // check for non-canonical values (after the check for special values)
3261 if ((C1.w[1] > 0x0001ed09bead87c0ull)
3262 || (C1.w[1] == 0x0001ed09bead87c0ull
3263 && (C1.w[0] > 0x378d8e63ffffffffull))
3264 || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
3265 res = 0x00000000;
3266 BID_RETURN (res);
3267 } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
3268 // x is 0
3269 res = 0x00000000;
3270 BID_RETURN (res);
3271 } else { // x is not special and is not zero
3272
3273 // q = nr. of decimal digits in x
3274 // determine first the nr. of bits in x
3275 if (C1.w[1] == 0) {
3276 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
3277 // split the 64-bit value in two 32-bit halves to avoid rounding errors
3278 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
3279 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
3280 x_nr_bits =
3281 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3282 } else { // x < 2^32
3283 tmp1.d = (double) (C1.w[0]); // exact conversion
3284 x_nr_bits =
3285 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3286 }
3287 } else { // if x < 2^53
3288 tmp1.d = (double) C1.w[0]; // exact conversion
3289 x_nr_bits =
3290 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3291 }
3292 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
3293 tmp1.d = (double) C1.w[1]; // exact conversion
3294 x_nr_bits =
3295 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3296 }
3297 q = nr_digits[x_nr_bits - 1].digits;
3298 if (q == 0) {
3299 q = nr_digits[x_nr_bits - 1].digits1;
3300 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
3301 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
3302 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
3303 q++;
3304 }
3305 exp = (x_exp >> 49) - 6176;
3306 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
3307 // set invalid flag
3308 *pfpsf |= INVALID_EXCEPTION;
3309 // return Integer Indefinite
3310 res = 0x80000000;
3311 BID_RETURN (res);
3312 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
3313 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
3314 // so x rounded to an integer may or may not fit in a signed 32-bit int
3315 // the cases that do not fit are identified here; the ones that fit
3316 // fall through and will be handled with other cases further,
3317 // under '1 <= q + exp <= 10'
3318 if (x_sign) { // if n < 0 and q + exp = 10
3319 // if n <= -1/2 then n is too large
3320 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1/2
3321 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05, 1<=q<=34
3322 if (q <= 11) {
3323 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
3324 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3325 if (tmp64 >= 0x05ull) {
3326 // set invalid flag
3327 *pfpsf |= INVALID_EXCEPTION;
3328 // return Integer Indefinite
3329 res = 0x80000000;
3330 BID_RETURN (res);
3331 }
3332 // else cases that can be rounded to a 32-bit int fall through
3333 // to '1 <= q + exp <= 10'
3334 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3335 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05 <=>
3336 // C >= 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
3337 // (scale 1/2 up)
3338 tmp64 = 0x05ull;
3339 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3340 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
3341 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3342 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
3343 }
3344 if (C1.w[1] > C.w[1]
3345 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
3346 // set invalid flag
3347 *pfpsf |= INVALID_EXCEPTION;
3348 // return Integer Indefinite
3349 res = 0x80000000;
3350 BID_RETURN (res);
3351 }
3352 // else cases that can be rounded to a 32-bit int fall through
3353 // to '1 <= q + exp <= 10'
3354 }
3355 } else { // if n > 0 and q + exp = 10
3356 // if n >= 2^32 - 1/2 then n is too large
3357 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
3358 // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
3359 if (q <= 11) {
3360 tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
3361 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3362 if (tmp64 >= 0x9fffffffbull) {
3363 // set invalid flag
3364 *pfpsf |= INVALID_EXCEPTION;
3365 // return Integer Indefinite
3366 res = 0x80000000;
3367 BID_RETURN (res);
3368 }
3369 // else cases that can be rounded to a 32-bit int fall through
3370 // to '1 <= q + exp <= 10'
3371 } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3372 // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
3373 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
3374 // (scale 2^32-1/2 up)
3375 tmp64 = 0x9fffffffbull;
3376 if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3377 __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
3378 } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3379 __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
3380 }
3381 if (C1.w[1] > C.w[1]
3382 || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
3383 // set invalid flag
3384 *pfpsf |= INVALID_EXCEPTION;
3385 // return Integer Indefinite
3386 res = 0x80000000;
3387 BID_RETURN (res);
3388 }
3389 // else cases that can be rounded to a 32-bit int fall through
3390 // to '1 <= q + exp <= 10'
3391 }
3392 }
3393 }
3394 // n is not too large to be converted to int32: -1/2 < n < 2^32 - 1/2
3395 // Note: some of the cases tested for above fall through to this point
3396 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
3397 // set inexact flag
3398 *pfpsf |= INEXACT_EXCEPTION;
3399 // return 0
3400 res = 0x00000000;
3401 BID_RETURN (res);
3402 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
3403 // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
3404 // res = 0
3405 // else if x > 0
3406 // res = +1
3407 // else // if x < 0
3408 // invalid exc
3409 ind = q - 1;
3410 if (ind <= 18) { // 0 <= ind <= 18
3411 if ((C1.w[1] == 0) && (C1.w[0] < midpoint64[ind])) {
3412 res = 0x00000000; // return 0
3413 } else if (!x_sign) { // n > 0
3414 res = 0x00000001; // return +1
3415 } else {
3416 res = 0x80000000;
3417 *pfpsf |= INVALID_EXCEPTION;
3418 BID_RETURN (res);
3419 }
3420 } else { // 19 <= ind <= 33
3421 if ((C1.w[1] < midpoint128[ind - 19].w[1])
3422 || ((C1.w[1] == midpoint128[ind - 19].w[1])
3423 && (C1.w[0] < midpoint128[ind - 19].w[0]))) {
3424 res = 0x00000000; // return 0
3425 } else if (!x_sign) { // n > 0
3426 res = 0x00000001; // return +1
3427 } else {
3428 res = 0x80000000;
3429 *pfpsf |= INVALID_EXCEPTION;
3430 BID_RETURN (res);
3431 }
3432 }
3433 // set inexact flag
3434 *pfpsf |= INEXACT_EXCEPTION;
3435 } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
3436 if (x_sign) { // x <= -1
3437 // set invalid flag
3438 *pfpsf |= INVALID_EXCEPTION;
3439 // return Integer Indefinite
3440 res = 0x80000000;
3441 BID_RETURN (res);
3442 }
3443 // 1 <= x < 2^31-1/2 so x can be rounded
3444 // to nearest-away to a 32-bit signed integer
3445 if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
3446 ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
3447 // chop off ind digits from the lower part of C1
3448 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
3449 tmp64 = C1.w[0];
3450 if (ind <= 19) {
3451 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
3452 } else {
3453 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
3454 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
3455 }
3456 if (C1.w[0] < tmp64)
3457 C1.w[1]++;
3458 // calculate C* and f*
3459 // C* is actually floor(C*) in this case
3460 // C* and f* need shifting and masking, as shown by
3461 // shiftright128[] and maskhigh128[]
3462 // 1 <= x <= 33
3463 // kx = 10^(-x) = ten2mk128[ind - 1]
3464 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
3465 // the approximation of 10^(-x) was rounded up to 118 bits
3466 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
3467 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
3468 Cstar.w[1] = P256.w[3];
3469 Cstar.w[0] = P256.w[2];
3470 fstar.w[3] = 0;
3471 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
3472 fstar.w[1] = P256.w[1];
3473 fstar.w[0] = P256.w[0];
3474 } else { // 22 <= ind - 1 <= 33
3475 Cstar.w[1] = 0;
3476 Cstar.w[0] = P256.w[3];
3477 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
3478 fstar.w[2] = P256.w[2];
3479 fstar.w[1] = P256.w[1];
3480 fstar.w[0] = P256.w[0];
3481 }
3482 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
3483 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
3484 // if (0 < f* < 10^(-x)) then the result is a midpoint
3485 // if floor(C*) is even then C* = floor(C*) - logical right
3486 // shift; C* has p decimal digits, correct by Prop. 1)
3487 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
3488 // shift; C* has p decimal digits, correct by Pr. 1)
3489 // else
3490 // C* = floor(C*) (logical right shift; C has p decimal digits,
3491 // correct by Property 1)
3492 // n = C* * 10^(e+x)
3493
3494 // shift right C* by Ex-128 = shiftright128[ind]
3495 shift = shiftright128[ind - 1]; // 0 <= shift <= 102
3496 if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
3497 Cstar.w[0] =
3498 (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
3499 // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
3500 } else { // 22 <= ind - 1 <= 33
3501 Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
3502 }
3503 // if the result was a midpoint, it was already rounded away from zero
3504 // determine inexactness of the rounding of C*
3505 // if (0 < f* - 1/2 < 10^(-x)) then
3506 // the result is exact
3507 // else // if (f* - 1/2 > T*) then
3508 // the result is inexact
3509 if (ind - 1 <= 2) {
3510 if (fstar.w[1] > 0x8000000000000000ull ||
3511 (fstar.w[1] == 0x8000000000000000ull
3512 && fstar.w[0] > 0x0ull)) {
3513 // f* > 1/2 and the result may be exact
3514 tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
3515 if (tmp64 > ten2mk128trunc[ind - 1].w[1]
3516 || (tmp64 == ten2mk128trunc[ind - 1].w[1]
3517 && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
3518 // set the inexact flag
3519 // *pfpsf |= INEXACT_EXCEPTION;
3520 tmp_inexact = 1;
3521 } // else the result is exact
3522 } else { // the result is inexact; f2* <= 1/2
3523 // set the inexact flag
3524 // *pfpsf |= INEXACT_EXCEPTION;
3525 tmp_inexact = 1;
3526 }
3527 } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
3528 if (fstar.w[3] > 0x0 ||
3529 (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
3530 (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
3531 (fstar.w[1] || fstar.w[0]))) {
3532 // f2* > 1/2 and the result may be exact
3533 // Calculate f2* - 1/2
3534 tmp64 = fstar.w[2] - onehalf128[ind - 1];
3535 tmp64A = fstar.w[3];
3536 if (tmp64 > fstar.w[2])
3537 tmp64A--;
3538 if (tmp64A || tmp64
3539 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
3540 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
3541 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
3542 // set the inexact flag
3543 // *pfpsf |= INEXACT_EXCEPTION;
3544 tmp_inexact = 1;
3545 } // else the result is exact
3546 } else { // the result is inexact; f2* <= 1/2
3547 // set the inexact flag
3548 // *pfpsf |= INEXACT_EXCEPTION;
3549 tmp_inexact = 1;
3550 }
3551 } else { // if 22 <= ind <= 33
3552 if (fstar.w[3] > onehalf128[ind - 1] ||
3553 (fstar.w[3] == onehalf128[ind - 1] &&
3554 (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
3555 // f2* > 1/2 and the result may be exact
3556 // Calculate f2* - 1/2
3557 tmp64 = fstar.w[3] - onehalf128[ind - 1];
3558 if (tmp64 || fstar.w[2]
3559 || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
3560 || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
3561 && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
3562 // set the inexact flag
3563 // *pfpsf |= INEXACT_EXCEPTION;
3564 tmp_inexact = 1;
3565 } // else the result is exact
3566 } else { // the result is inexact; f2* <= 1/2
3567 // set the inexact flag
3568 // *pfpsf |= INEXACT_EXCEPTION;
3569 tmp_inexact = 1;
3570 }
3571 }
3572 // no need to check for midpoints - already rounded away from zero!
3573 res = Cstar.w[0]; // the result is positive
3574 if (tmp_inexact)
3575 *pfpsf |= INEXACT_EXCEPTION;
3576 } else if (exp == 0) {
3577 // 1 <= q <= 10
3578 // res = +C (exact)
3579 res = C1.w[0];
3580 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
3581 // res = +C * 10^exp (exact)
3582 res = C1.w[0] * ten2k64[exp];
3583 }
3584 }
3585 }
3586
3587 BID_RETURN (res);
3588 }