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