]> git.ipfire.org Git - thirdparty/gcc.git/blame - libgcc/config/libbid/bid128_to_uint32.c
projects / thirdparty / gcc.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update copyright years.
[thirdparty/gcc.git] / libgcc / config / libbid / bid128_to_uint32.c
CommitLineData
85ec4feb 1/* Copyright (C) 2007-2018 Free Software Foundation, Inc.
200359e8
L		 2
3This file is part of GCC.
4
5GCC is free software; you can redistribute it and/or modify it under
6the terms of the GNU General Public License as published by the Free
748086b7 7Software Foundation; either version 3, or (at your option) any later
200359e8
L		 8version.
9
200359e8
L		 10GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11WARRANTY; without even the implied warranty of MERCHANTABILITY or
12FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13for more details.
14
748086b7
JJ		 15Under Section 7 of GPL version 3, you are granted additional
16permissions described in the GCC Runtime Library Exception, version
173.1, as published by the Free Software Foundation.
18
19You should have received a copy of the GNU General Public License and
20a copy of the GCC Runtime Library Exception along with this program;
21see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
22<http://www.gnu.org/licenses/>. */
200359e8
L		 23
24#include "bid_internal.h"
25
26/*****************************************************************************
27 * BID128_to_uint32_rnint
28 ****************************************************************************/
29
b2a00c89
L		 30BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
31 bid128_to_uint32_rnint, x)
200359e8 32
b2a00c89
L		 33 unsigned int res;
34 UINT64 x_sign;
35 UINT64 x_exp;
36 int exp; // unbiased exponent
200359e8 37 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 46
47 // unpack x
b2a00c89
L		 48x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
49x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
50C1.w[1] = x.w[1] & MASK_COEFF;
51C1.w[0] = x.w[0];
200359e8
L		 52
53 // check for NaN or Infinity
b2a00c89 54if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 55 // x is special
b2a00c89
L		 56if ((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;
200359e8 80 }
b2a00c89
L		 81 BID_RETURN (res);
82}
83}
200359e8 84 // check for non-canonical values (after the check for special values)
b2a00c89
L		 85if ((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
200359e8
L		 108 x_nr_bits =
109 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
110 }
b2a00c89
L		 111 } else { // if x < 2^53
112 tmp1.d = (double) C1.w[0]; // exact conversion
200359e8 113 x_nr_bits =
b2a00c89 114 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200359e8 115 }
b2a00c89
L		 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
200359e8 153 res = 0x80000000;
b2a00c89
L		 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
200359e8 170 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 171 // return Integer Indefinite
172 res = 0x80000000;
173 BID_RETURN (res);
200359e8 174 }
b2a00c89
L		 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
200359e8 189 res = 0x80000000;
b2a00c89
L		 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
200359e8 207 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 208 // return Integer Indefinite
209 res = 0x80000000;
210 BID_RETURN (res);
200359e8 211 }
b2a00c89
L		 212 // else cases that can be rounded to a 32-bit int fall through
213 // to '1 <= q + exp <= 10'
200359e8 214 }
b2a00c89
L		 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;
200359e8 238 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 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 {
200359e8 248 res = 0x80000000;
b2a00c89
L		 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];
200359e8 272 }
b2a00c89
L		 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];
200359e8 298 }
b2a00c89
L		 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];
200359e8
L		 342 }
343 }
b2a00c89
L		 344}
345
346BID_RETURN (res);
200359e8
L		 347}
348
349/*****************************************************************************
350 * BID128_to_uint32_xrnint
351 ****************************************************************************/
352
b2a00c89
L		 353BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
354 bid128_to_uint32_xrnint, x)
200359e8 355
b2a00c89
L		 356 unsigned int res;
357 UINT64 x_sign;
358 UINT64 x_exp;
359 int exp; // unbiased exponent
200359e8 360 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 370
371 // unpack x
b2a00c89
L		 372x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
373x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
374C1.w[1] = x.w[1] & MASK_COEFF;
375C1.w[0] = x.w[0];
200359e8
L		 376
377 // check for NaN or Infinity
b2a00c89 378if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 379 // x is special
b2a00c89
L		 380if ((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;
200359e8 404 }
b2a00c89
L		 405 BID_RETURN (res);
406}
407}
200359e8 408 // check for non-canonical values (after the check for special values)
b2a00c89
L		 409if ((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
200359e8
L		 432 x_nr_bits =
433 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
434 }
b2a00c89
L		 435 } else { // if x < 2^53
436 tmp1.d = (double) C1.w[0]; // exact conversion
200359e8 437 x_nr_bits =
b2a00c89 438 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200359e8 439 }
b2a00c89
L		 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
200359e8 477 res = 0x80000000;
b2a00c89
L		 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
200359e8 494 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 495 // return Integer Indefinite
496 res = 0x80000000;
200359e8
L		 497 BID_RETURN (res);
498 }
b2a00c89
L		 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
200359e8 513 res = 0x80000000;
b2a00c89
L		 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
200359e8 531 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 532 // return Integer Indefinite
533 res = 0x80000000;
200359e8
L		 534 BID_RETURN (res);
535 }
b2a00c89
L		 536 // else cases that can be rounded to a 32-bit int fall through
537 // to '1 <= q + exp <= 10'
200359e8 538 }
b2a00c89
L		 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;
200359e8 564 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 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 {
200359e8 575 res = 0x80000000;
b2a00c89 576 *pfpsf |= INVALID_EXCEPTION;
200359e8
L		 577 BID_RETURN (res);
578 }
b2a00c89
L		 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])) {
200359e8
L		 664 // set the inexact flag
665 // *pfpsf |= INEXACT_EXCEPTION;
666 tmp_inexact = 1;
b2a00c89
L		 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])) {
200359e8
L		 688 // set the inexact flag
689 // *pfpsf |= INEXACT_EXCEPTION;
690 tmp_inexact = 1;
b2a00c89
L		 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])) {
200359e8
L		 708 // set the inexact flag
709 // *pfpsf |= INEXACT_EXCEPTION;
710 tmp_inexact = 1;
b2a00c89
L		 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]
200359e8 732 }
b2a00c89
L		 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];
200359e8
L		 743 }
744 }
b2a00c89
L		 745}
746
747BID_RETURN (res);
200359e8
L		 748}
749
750/*****************************************************************************
751 * BID128_to_uint32_floor
752 ****************************************************************************/
753
b2a00c89
L		 754BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
755 bid128_to_uint32_floor, x)
200359e8 756
b2a00c89
L		 757 unsigned int res;
758 UINT64 x_sign;
759 UINT64 x_exp;
760 int exp; // unbiased exponent
200359e8 761 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 772
773 // unpack x
b2a00c89
L		 774x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
775x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
776C1.w[1] = x.w[1] & MASK_COEFF;
777C1.w[0] = x.w[0];
200359e8
L		 778
779 // check for NaN or Infinity
b2a00c89 780if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 781 // x is special
b2a00c89
L		 782if ((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)
811if ((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);
200359e8 844 }
b2a00c89
L		 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) {
200359e8
L		 885 // set invalid flag
886 *pfpsf |= INVALID_EXCEPTION;
887 // return Integer Indefinite
888 res = 0x80000000;
b2a00c89
L		 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])) {
200359e8
L		 904 // set invalid flag
905 *pfpsf |= INVALID_EXCEPTION;
906 // return Integer Indefinite
907 res = 0x80000000;
b2a00c89 908 BID_RETURN (res);
200359e8 909 }
b2a00c89
L		 910 // else cases that can be rounded to a 32-bit int fall through
911 // to '1 <= q + exp <= 10'
200359e8
L		 912 }
913 }
b2a00c89
L		 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
200359e8
L		 919 res = 0x00000000;
920 BID_RETURN (res);
b2a00c89
L		 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];
200359e8 933 }
b2a00c89
L		 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;
200359e8 1032 }
b2a00c89
L		 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;
200359e8 1051 }
200359e8 1052 }
b2a00c89
L		 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
200359e8 1058 }
b2a00c89
L		 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];
200359e8
L		 1067 }
1068 }
b2a00c89
L		 1069}
1070
1071BID_RETURN (res);
200359e8
L		 1072}
1073
1074/*****************************************************************************
1075 * BID128_to_uint32_xfloor
1076 ****************************************************************************/
1077
b2a00c89
L		 1078BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
1079 bid128_to_uint32_xfloor, x)
200359e8 1080
b2a00c89
L		 1081 unsigned int res;
1082 UINT64 x_sign;
1083 UINT64 x_exp;
1084 int exp; // unbiased exponent
200359e8 1085 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 1096
1097 // unpack x
b2a00c89
L		 1098x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1099x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
1100C1.w[1] = x.w[1] & MASK_COEFF;
1101C1.w[0] = x.w[0];
200359e8
L		 1102
1103 // check for NaN or Infinity
b2a00c89 1104if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 1105 // x is special
b2a00c89
L		 1106if ((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)
1135if ((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);
200359e8 1168 }
b2a00c89
L		 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) {
200359e8
L		 1208 // set invalid flag
1209 *pfpsf |= INVALID_EXCEPTION;
1210 // return Integer Indefinite
1211 res = 0x80000000;
b2a00c89
L		 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])) {
200359e8
L		 1227 // set invalid flag
1228 *pfpsf |= INVALID_EXCEPTION;
1229 // return Integer Indefinite
1230 res = 0x80000000;
b2a00c89 1231 BID_RETURN (res);
200359e8 1232 }
b2a00c89
L		 1233 // else cases that can be rounded to a 32-bit int fall through
1234 // to '1 <= q + exp <= 10'
200359e8
L		 1235 }
1236 }
b2a00c89
L		 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
200359e8
L		 1244 res = 0x00000000;
1245 BID_RETURN (res);
b2a00c89
L		 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];
200359e8 1258 }
b2a00c89
L		 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];
200359e8 1284 }
b2a00c89
L		 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])) {
200359e8
L		 1320 // set the inexact flag
1321 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L		 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])) {
200359e8
L		 1343 // set the inexact flag
1344 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L		 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])) {
200359e8
L		 1362 // set the inexact flag
1363 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L		 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 }
200359e8 1389 }
b2a00c89
L		 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];
200359e8
L		 1404 }
1405 }
b2a00c89
L		 1406}
1407
1408BID_RETURN (res);
200359e8
L		 1409}
1410
1411/*****************************************************************************
1412 * BID128_to_uint32_ceil
1413 ****************************************************************************/
1414
b2a00c89
L		 1415BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
1416 bid128_to_uint32_ceil, x)
200359e8 1417
b2a00c89
L		 1418 unsigned int res;
1419 UINT64 x_sign;
1420 UINT64 x_exp;
1421 int exp; // unbiased exponent
200359e8 1422 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 1433
1434 // unpack x
b2a00c89
L		 1435x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1436x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
1437C1.w[1] = x.w[1] & MASK_COEFF;
1438C1.w[0] = x.w[0];
200359e8
L		 1439
1440 // check for NaN or Infinity
b2a00c89 1441if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 1442 // x is special
b2a00c89
L		 1443if ((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;
200359e8 1454 }
b2a00c89
L		 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}
200359e8 1471 // check for non-canonical values (after the check for special values)
b2a00c89
L		 1472if ((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
200359e8
L		 1495 x_nr_bits =
1496 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1497 }
b2a00c89
L		 1498 } else { // if x < 2^53
1499 tmp1.d = (double) C1.w[0]; // exact conversion
200359e8 1500 x_nr_bits =
b2a00c89 1501 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200359e8 1502 }
b2a00c89
L		 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;
200359e8 1517
b2a00c89
L		 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
200359e8
L		 1619 // set invalid flag
1620 *pfpsf |= INVALID_EXCEPTION;
1621 // return Integer Indefinite
1622 res = 0x80000000;
1623 BID_RETURN (res);
b2a00c89
L		 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 ;
200359e8
L		 1737 }
1738 }
b2a00c89
L		 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 }
200359e8 1757 }
b2a00c89
L		 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
200359e8 1763 }
b2a00c89
L		 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];
200359e8
L		 1772 }
1773 }
b2a00c89
L		 1774}
1775
1776BID_RETURN (res);
200359e8
L		 1777}
1778
1779/*****************************************************************************
1780 * BID128_to_uint32_xceil
1781 ****************************************************************************/
1782
b2a00c89
L		 1783BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
1784 bid128_to_uint32_xceil, x)
200359e8 1785
b2a00c89
L		 1786 unsigned int res;
1787 UINT64 x_sign;
1788 UINT64 x_exp;
1789 int exp; // unbiased exponent
200359e8 1790 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 1801
1802 // unpack x
b2a00c89
L		 1803x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1804x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
1805C1.w[1] = x.w[1] & MASK_COEFF;
1806C1.w[0] = x.w[0];
200359e8
L		 1807
1808 // check for NaN or Infinity
b2a00c89 1809if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 1810 // x is special
b2a00c89
L		 1811if ((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;
200359e8 1835 }
b2a00c89
L		 1836 BID_RETURN (res);
1837}
1838}
200359e8 1839 // check for non-canonical values (after the check for special values)
b2a00c89
L		 1840if ((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
200359e8
L		 1863 x_nr_bits =
1864 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1865 }
b2a00c89
L		 1866 } else { // if x < 2^53
1867 tmp1.d = (double) C1.w[0]; // exact conversion
200359e8 1868 x_nr_bits =
b2a00c89 1869 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200359e8 1870 }
b2a00c89
L		 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 }
200359e8 1970 }
b2a00c89
L		 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
200359e8
L		 1988 // set invalid flag
1989 *pfpsf |= INVALID_EXCEPTION;
1990 // return Integer Indefinite
1991 res = 0x80000000;
1992 BID_RETURN (res);
200359e8 1993 }
b2a00c89
L		 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];
200359e8 2031 }
b2a00c89
L		 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])) {
200359e8
L		 2067 // set the inexact flag
2068 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L		 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])) {
200359e8
L		 2090 // set the inexact flag
2091 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L		 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])) {
200359e8
L		 2109 // set the inexact flag
2110 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L		 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 }
200359e8 2136 }
b2a00c89
L		 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];
200359e8
L		 2151 }
2152 }
b2a00c89
L		 2153}
2154
2155BID_RETURN (res);
200359e8
L		 2156}
2157
2158/*****************************************************************************
2159 * BID128_to_uint32_int
2160 ****************************************************************************/
2161
b2a00c89
L		 2162BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
2163 bid128_to_uint32_int, x)
200359e8 2164
b2a00c89
L		 2165 int res;
2166 UINT64 x_sign;
2167 UINT64 x_exp;
2168 int exp; // unbiased exponent
200359e8 2169 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 2180
2181 // unpack x
b2a00c89
L		 2182x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
2183x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
2184C1.w[1] = x.w[1] & MASK_COEFF;
2185C1.w[0] = x.w[0];
200359e8
L		 2186
2187 // check for NaN or Infinity
b2a00c89 2188if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 2189 // x is special
b2a00c89
L		 2190if ((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;
200359e8 2214 }
b2a00c89
L		 2215 BID_RETURN (res);
2216}
2217}
200359e8 2218 // check for non-canonical values (after the check for special values)
b2a00c89
L		 2219if ((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
200359e8
L		 2242 x_nr_bits =
2243 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2244 }
b2a00c89
L		 2245 } else { // if x < 2^53
2246 tmp1.d = (double) C1.w[0]; // exact conversion
200359e8 2247 x_nr_bits =
b2a00c89 2248 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200359e8 2249 }
b2a00c89
L		 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 }
200359e8 2350 }
b2a00c89
L		 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
200359e8
L		 2362 // set invalid flag
2363 *pfpsf |= INVALID_EXCEPTION;
2364 // return Integer Indefinite
2365 res = 0x80000000;
2366 BID_RETURN (res);
b2a00c89
L		 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;
200359e8
L		 2479 }
2480 }
b2a00c89
L		 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 }
200359e8 2499 }
b2a00c89
L		 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
200359e8 2505 }
b2a00c89
L		 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];
200359e8
L		 2514 }
2515 }
b2a00c89
L		 2516}
2517
2518BID_RETURN (res);
200359e8
L		 2519}
2520
2521/*****************************************************************************
2522 * BID128_to_uint32_xint
2523 ****************************************************************************/
2524
b2a00c89
L		 2525BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
2526 bid128_to_uint32_xint, x)
200359e8 2527
b2a00c89
L		 2528 int res;
2529 UINT64 x_sign;
2530 UINT64 x_exp;
2531 int exp; // unbiased exponent
200359e8 2532 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 2543
2544 // unpack x
b2a00c89
L		 2545x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
2546x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
2547C1.w[1] = x.w[1] & MASK_COEFF;
2548C1.w[0] = x.w[0];
200359e8
L		 2549
2550 // check for NaN or Infinity
b2a00c89 2551if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 2552 // x is special
b2a00c89
L		 2553if ((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)
2582if ((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);
200359e8 2607 }
b2a00c89
L		 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'
200359e8 2712 }
200359e8
L		 2713 }
2714 }
b2a00c89
L		 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
200359e8
L		 2722 res = 0x00000000;
2723 BID_RETURN (res);
b2a00c89
L		 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
200359e8
L		 2727 // set invalid flag
2728 *pfpsf |= INVALID_EXCEPTION;
2729 // return Integer Indefinite
2730 res = 0x80000000;
2731 BID_RETURN (res);
200359e8 2732 }
b2a00c89
L		 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];
200359e8 2771 }
b2a00c89
L		 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])) {
200359e8
L		 2807 // set the inexact flag
2808 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L		 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])) {
200359e8
L		 2830 // set the inexact flag
2831 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L		 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])) {
200359e8
L		 2849 // set the inexact flag
2850 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L		 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 }
200359e8 2876 }
b2a00c89
L		 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];
200359e8
L		 2891 }
2892 }
b2a00c89
L		 2893}
2894
2895BID_RETURN (res);
200359e8
L		 2896}
2897
2898/*****************************************************************************
2899 * BID128_to_uint32_rninta
2900 ****************************************************************************/
2901
b2a00c89
L		 2902BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
2903 bid128_to_uint32_rninta, x)
200359e8 2904
b2a00c89
L		 2905 unsigned int res;
2906 UINT64 x_sign;
2907 UINT64 x_exp;
2908 int exp; // unbiased exponent
200359e8 2909 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 2917
2918 // unpack x
b2a00c89
L		 2919x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
2920x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
2921C1.w[1] = x.w[1] & MASK_COEFF;
2922C1.w[0] = x.w[0];
200359e8
L		 2923
2924 // check for NaN or Infinity
b2a00c89 2925if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 2926 // x is special
b2a00c89
L		 2927if ((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;
200359e8 2951 }
b2a00c89
L		 2952 BID_RETURN (res);
2953}
2954}
200359e8 2955 // check for non-canonical values (after the check for special values)
b2a00c89
L		 2956if ((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
200359e8
L		 2979 x_nr_bits =
2980 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2981 }
b2a00c89
L		 2982 } else { // if x < 2^53
2983 tmp1.d = (double) C1.w[0]; // exact conversion
200359e8 2984 x_nr_bits =
b2a00c89 2985 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200359e8 2986 }
b2a00c89
L		 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
200359e8 3024 res = 0x80000000;
b2a00c89
L		 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
200359e8 3042 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 3043 // return Integer Indefinite
3044 res = 0x80000000;
3045 BID_RETURN (res);
200359e8 3046 }
b2a00c89
L		 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
200359e8 3061 res = 0x80000000;
b2a00c89
L		 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
200359e8 3079 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 3080 // return Integer Indefinite
3081 res = 0x80000000;
3082 BID_RETURN (res);
200359e8 3083 }
b2a00c89
L		 3084 // else cases that can be rounded to a 32-bit int fall through
3085 // to '1 <= q + exp <= 10'
200359e8 3086 }
b2a00c89
L		 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;
200359e8 3110 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 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 {
200359e8 3120 res = 0x80000000;
b2a00c89 3121 *pfpsf |= INVALID_EXCEPTION;
200359e8 3122 }
b2a00c89
L		 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];
200359e8 3162 }
b2a00c89
L		 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];
200359e8
L		 3194 }
3195 }
b2a00c89
L		 3196}
3197
3198BID_RETURN (res);
200359e8
L		 3199}
3200
3201/*****************************************************************************
3202 * BID128_to_uint32_xrninta
3203 ****************************************************************************/
3204
b2a00c89
L		 3205BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
3206 bid128_to_uint32_xrninta, x)
200359e8 3207
b2a00c89
L		 3208 unsigned int res;
3209 UINT64 x_sign;
3210 UINT64 x_exp;
3211 int exp; // unbiased exponent
200359e8 3212 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 3222
3223 // unpack x
b2a00c89
L		 3224x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
3225x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
3226C1.w[1] = x.w[1] & MASK_COEFF;
3227C1.w[0] = x.w[0];
200359e8
L		 3228
3229 // check for NaN or Infinity
b2a00c89 3230if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 3231 // x is special
b2a00c89
L		 3232if ((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;
200359e8 3256 }
b2a00c89
L		 3257 BID_RETURN (res);
3258}
3259}
200359e8 3260 // check for non-canonical values (after the check for special values)
b2a00c89
L		 3261if ((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
200359e8
L		 3284 x_nr_bits =
3285 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3286 }
b2a00c89
L		 3287 } else { // if x < 2^53
3288 tmp1.d = (double) C1.w[0]; // exact conversion
200359e8 3289 x_nr_bits =
b2a00c89 3290 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200359e8 3291 }
b2a00c89
L		 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
200359e8 3329 res = 0x80000000;
b2a00c89
L		 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
200359e8 3347 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 3348 // return Integer Indefinite
3349 res = 0x80000000;
200359e8
L		 3350 BID_RETURN (res);
3351 }
b2a00c89
L		 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
200359e8 3366 res = 0x80000000;
b2a00c89
L		 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
200359e8 3384 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 3385 // return Integer Indefinite
3386 res = 0x80000000;
200359e8
L		 3387 BID_RETURN (res);
3388 }
b2a00c89
L		 3389 // else cases that can be rounded to a 32-bit int fall through
3390 // to '1 <= q + exp <= 10'
200359e8 3391 }
b2a00c89
L		 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;
200359e8 3417 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 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 {
200359e8 3428 res = 0x80000000;
b2a00c89 3429 *pfpsf |= INVALID_EXCEPTION;
200359e8
L		 3430 BID_RETURN (res);
3431 }
b2a00c89
L		 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])) {
200359e8
L		 3518 // set the inexact flag
3519 // *pfpsf |= INEXACT_EXCEPTION;
3520 tmp_inexact = 1;
b2a00c89
L		 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])) {
200359e8
L		 3542 // set the inexact flag
3543 // *pfpsf |= INEXACT_EXCEPTION;
3544 tmp_inexact = 1;
b2a00c89
L		 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])) {
200359e8
L		 3562 // set the inexact flag
3563 // *pfpsf |= INEXACT_EXCEPTION;
3564 tmp_inexact = 1;
b2a00c89
L		 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;
200359e8 3570 }
200359e8 3571 }
b2a00c89
L		 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];
200359e8
L		 3583 }
3584 }
b2a00c89
L		 3585}
3586
3587BID_RETURN (res);
200359e8 3588}
Mirror of https://gcc.gnu.org/git/gcc.git