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