]> git.ipfire.org Git - thirdparty/gcc.git/blame - libgcc/config/libbid/bid128_to_uint64.c
projects / thirdparty / gcc.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update copyright years.
[thirdparty/gcc.git] / libgcc / config / libbid / bid128_to_uint64.c
CommitLineData
8d9254fc 1/* Copyright (C) 2007-2020 Free Software Foundation, Inc.
200359e8
L		 2
3This file is part of GCC.
4
5GCC is free software; you can redistribute it and/or modify it under
6the terms of the GNU General Public License as published by the Free
748086b7 7Software Foundation; either version 3, or (at your option) any later
200359e8
L		 8version.
9
200359e8
L		 10GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11WARRANTY; without even the implied warranty of MERCHANTABILITY or
12FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13for more details.
14
748086b7
JJ		 15Under Section 7 of GPL version 3, you are granted additional
16permissions described in the GCC Runtime Library Exception, version
173.1, as published by the Free Software Foundation.
18
19You should have received a copy of the GNU General Public License and
20a copy of the GCC Runtime Library Exception along with this program;
21see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
22<http://www.gnu.org/licenses/>. */
200359e8
L		 23
24#include "bid_internal.h"
25
26/*****************************************************************************
27 * BID128_to_uint64_rnint
28 ****************************************************************************/
29
b2a00c89
L		 30BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
31 bid128_to_uint64_rnint, x)
200359e8 32
b2a00c89
L		 33 UINT64 res;
34 UINT64 x_sign;
35 UINT64 x_exp;
36 int exp; // unbiased exponent
200359e8 37 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 38 UINT64 tmp64;
39 BID_UI64DOUBLE tmp1;
40 unsigned int x_nr_bits;
41 int q, ind, shift;
42 UINT128 C1, C;
43 UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
44 UINT256 fstar;
45 UINT256 P256;
200359e8
L		 46
47 // unpack x
b2a00c89
L		 48x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
49x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
50C1.w[1] = x.w[1] & MASK_COEFF;
51C1.w[0] = x.w[0];
200359e8
L		 52
53 // check for NaN or Infinity
b2a00c89 54if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 55 // x is special
b2a00c89
L		 56if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
57 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
58 // set invalid flag
59 *pfpsf |= INVALID_EXCEPTION;
60 // return Integer Indefinite
61 res = 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;
200359e8 80 }
b2a00c89
L		 81 BID_RETURN (res);
82}
83}
200359e8 84 // check for non-canonical values (after the check for special values)
b2a00c89
L		 85if ((C1.w[1] > 0x0001ed09bead87c0ull)
86 || (C1.w[1] == 0x0001ed09bead87c0ull
87 && (C1.w[0] > 0x378d8e63ffffffffull))
88 || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
89 res = 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
200359e8
L		 108 x_nr_bits =
109 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
110 }
b2a00c89
L		 111 } else { // if x < 2^53
112 tmp1.d = (double) C1.w[0]; // exact conversion
200359e8 113 x_nr_bits =
b2a00c89 114 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200359e8 115 }
b2a00c89
L		 116 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
117 tmp1.d = (double) C1.w[1]; // exact conversion
118 x_nr_bits =
119 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
120 }
121 q = nr_digits[x_nr_bits - 1].digits;
122 if (q == 0) {
123 q = nr_digits[x_nr_bits - 1].digits1;
124 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
125 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
126 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
127 q++;
128 }
129 exp = (x_exp >> 49) - 6176;
200359e8 130
b2a00c89
L		 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) {
200359e8
L		 151 // set invalid flag
152 *pfpsf |= INVALID_EXCEPTION;
153 // return Integer Indefinite
154 res = 0x8000000000000000ull;
155 BID_RETURN (res);
156 }
b2a00c89
L		 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);
200359e8 169 }
b2a00c89
L		 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
200359e8 184 res = 0x8000000000000000ull;
b2a00c89
L		 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
200359e8 211 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 212 // return Integer Indefinite
213 res = 0x8000000000000000ull;
214 BID_RETURN (res);
200359e8 215 }
b2a00c89
L		 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
200359e8 225 res = 0x8000000000000000ull;
b2a00c89
L		 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
200359e8 238 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 239 // return Integer Indefinite
240 res = 0x8000000000000000ull;
241 BID_RETURN (res);
200359e8 242 }
b2a00c89
L		 243 // else cases that can be rounded to a 64-bit int fall through
244 // to '1 <= q + exp <= 20'
200359e8 245 }
b2a00c89
L		 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;
200359e8 269 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 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 {
200359e8 279 res = 0x8000000000000000ull;
b2a00c89 280 *pfpsf |= INVALID_EXCEPTION;
200359e8 281 }
b2a00c89
L		 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];
200359e8 302 } else {
b2a00c89
L		 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]
200359e8 366 }
b2a00c89
L		 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];
200359e8
L		 376 }
377 }
b2a00c89
L		 378}
379
380BID_RETURN (res);
200359e8
L		 381}
382
383/*****************************************************************************
384 * BID128_to_uint64_xrnint
385 ****************************************************************************/
386
b2a00c89
L		 387BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
388 bid128_to_uint64_xrnint, x)
200359e8 389
b2a00c89
L		 390 UINT64 res;
391 UINT64 x_sign;
392 UINT64 x_exp;
393 int exp; // unbiased exponent
200359e8 394 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 403
404 // unpack x
b2a00c89
L		 405x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
406x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
407C1.w[1] = x.w[1] & MASK_COEFF;
408C1.w[0] = x.w[0];
200359e8
L		 409
410 // check for NaN or Infinity
b2a00c89 411if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 412 // x is special
b2a00c89
L		 413if ((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;
200359e8 424 }
b2a00c89
L		 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}
200359e8 441 // check for non-canonical values (after the check for special values)
b2a00c89
L		 442if ((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
200359e8
L		 465 x_nr_bits =
466 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
467 }
b2a00c89
L		 468 } else { // if x < 2^53
469 tmp1.d = (double) C1.w[0]; // exact conversion
200359e8 470 x_nr_bits =
b2a00c89 471 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200359e8 472 }
b2a00c89
L		 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;
200359e8 487
b2a00c89
L		 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) {
200359e8
L		 508 // set invalid flag
509 *pfpsf |= INVALID_EXCEPTION;
510 // return Integer Indefinite
511 res = 0x8000000000000000ull;
512 BID_RETURN (res);
513 }
b2a00c89
L		 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);
200359e8 526 }
b2a00c89
L		 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
200359e8 554 res = 0x8000000000000000ull;
b2a00c89
L		 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
200359e8 568 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 569 // return Integer Indefinite
570 res = 0x8000000000000000ull;
200359e8
L		 571 BID_RETURN (res);
572 }
b2a00c89
L		 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
200359e8 582 res = 0x8000000000000000ull;
b2a00c89
L		 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
200359e8 595 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 596 // return Integer Indefinite
597 res = 0x8000000000000000ull;
200359e8
L		 598 BID_RETURN (res);
599 }
b2a00c89
L		 600 // else cases that can be rounded to a 64-bit int fall through
601 // to '1 <= q + exp <= 20'
200359e8 602 }
b2a00c89
L		 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;
200359e8 628 *pfpsf |= INVALID_EXCEPTION;
b2a00c89
L		 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 {
200359e8 639 res = 0x8000000000000000ull;
b2a00c89 640 *pfpsf |= INVALID_EXCEPTION;
200359e8
L		 641 BID_RETURN (res);
642 }
b2a00c89
L		 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])) {
200359e8
L		 730 // set the inexact flag
731 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L		 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])) {
200359e8
L		 752 // set the inexact flag
753 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L		 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])) {
200359e8
L		 770 // set the inexact flag
771 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L		 772 } // else the result is exact
773 } else { // the result is inexact; f2* <= 1/2
774 // set the inexact flag
775 *pfpsf |= INEXACT_EXCEPTION;
200359e8 776 }
b2a00c89 777 }
200359e8 778
b2a00c89
L		 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]
200359e8 792 }
b2a00c89
L		 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];
200359e8
L		 802 }
803 }
b2a00c89
L		 804}
805
806BID_RETURN (res);
200359e8
L		 807}
808
809/*****************************************************************************
810 * BID128_to_uint64_floor
811 ****************************************************************************/
812
b2a00c89
L		 813BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
814 bid128_to_uint64_floor, x)
200359e8 815
b2a00c89
L		 816 UINT64 res;
817 UINT64 x_sign;
818 UINT64 x_exp;
819 int exp; // unbiased exponent
200359e8 820 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 827
828 // unpack x
b2a00c89
L		 829x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
830x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
831C1.w[1] = x.w[1] & MASK_COEFF;
832C1.w[0] = x.w[0];
200359e8
L		 833
834 // check for NaN or Infinity
b2a00c89 835if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 836 // x is special
b2a00c89
L		 837if ((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)
866if ((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) {
200359e8
L		 943 // set invalid flag
944 *pfpsf |= INVALID_EXCEPTION;
945 // return Integer Indefinite
946 res = 0x8000000000000000ull;
b2a00c89
L		 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) {
200359e8
L		 956 // set invalid flag
957 *pfpsf |= INVALID_EXCEPTION;
958 // return Integer Indefinite
959 res = 0x8000000000000000ull;
b2a00c89 960 BID_RETURN (res);
200359e8 961 }
b2a00c89
L		 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) {
200359e8
L		 968 // set invalid flag
969 *pfpsf |= INVALID_EXCEPTION;
970 // return Integer Indefinite
971 res = 0x8000000000000000ull;
b2a00c89
L		 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) {
200359e8
L		 980 // set invalid flag
981 *pfpsf |= INVALID_EXCEPTION;
982 // return Integer Indefinite
983 res = 0x8000000000000000ull;
b2a00c89 984 BID_RETURN (res);
200359e8 985 }
b2a00c89
L		 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'
200359e8
L		 1002 }
1003 }
b2a00c89
L		 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
200359e8
L		 1008 res = 0x0000000000000000ull;
1009 BID_RETURN (res);
b2a00c89
L		 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];
200359e8 1032 }
b2a00c89
L		 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
200359e8 1047 }
b2a00c89
L		 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];
200359e8
L		 1057 }
1058 }
b2a00c89
L		 1059}
1060
1061BID_RETURN (res);
200359e8
L		 1062}
1063
1064/*****************************************************************************
1065 * BID128_to_uint64_xfloor
1066 ****************************************************************************/
1067
b2a00c89
L		 1068BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
1069 bid128_to_uint64_xfloor, x)
200359e8 1070
b2a00c89
L		 1071 UINT64 res;
1072 UINT64 x_sign;
1073 UINT64 x_exp;
1074 int exp; // unbiased exponent
200359e8 1075 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 1083
1084 // unpack x
b2a00c89
L		 1085x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1086x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
1087C1.w[1] = x.w[1] & MASK_COEFF;
1088C1.w[0] = x.w[0];
200359e8
L		 1089
1090 // check for NaN or Infinity
b2a00c89 1091if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 1092 // x is special
b2a00c89
L		 1093if ((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)
1122if ((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) {
200359e8
L		 1199 // set invalid flag
1200 *pfpsf |= INVALID_EXCEPTION;
1201 // return Integer Indefinite
1202 res = 0x8000000000000000ull;
b2a00c89
L		 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) {
200359e8
L		 1212 // set invalid flag
1213 *pfpsf |= INVALID_EXCEPTION;
1214 // return Integer Indefinite
1215 res = 0x8000000000000000ull;
b2a00c89 1216 BID_RETURN (res);
200359e8 1217 }
b2a00c89
L		 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) {
200359e8
L		 1224 // set invalid flag
1225 *pfpsf |= INVALID_EXCEPTION;
1226 // return Integer Indefinite
1227 res = 0x8000000000000000ull;
b2a00c89
L		 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) {
200359e8
L		 1236 // set invalid flag
1237 *pfpsf |= INVALID_EXCEPTION;
1238 // return Integer Indefinite
1239 res = 0x8000000000000000ull;
b2a00c89 1240 BID_RETURN (res);
200359e8 1241 }
b2a00c89
L		 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'
200359e8
L		 1258 }
1259 }
b2a00c89
L		 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
200359e8
L		 1266 res = 0x0000000000000000ull;
1267 BID_RETURN (res);
b2a00c89
L		 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 }
200359e8 1342
b2a00c89
L		 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];
200359e8 1352 }
b2a00c89
L		 1353 }
1354}
1355
1356BID_RETURN (res);
1357}
1358
1359/*****************************************************************************
1360 * BID128_to_uint64_ceil
1361 ****************************************************************************/
1362
1363BID128_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
1380x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1381x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
1382C1.w[1] = x.w[1] & MASK_COEFF;
1383C1.w[0] = x.w[0];
1384
1385 // check for NaN or Infinity
1386if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1387 // x is special
1388if ((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)
1417if ((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
200359e8
L		 1440 x_nr_bits =
1441 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1442 }
b2a00c89
L		 1443 } else { // if x < 2^53
1444 tmp1.d = (double) C1.w[0]; // exact conversion
200359e8 1445 x_nr_bits =
b2a00c89 1446 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200359e8 1447 }
b2a00c89
L		 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
200359e8 1507 if (q == 1) {
b2a00c89
L		 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)) {
200359e8
L		 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) {
b2a00c89
L		 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)) {
200359e8
L		 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) {
b2a00c89
L		 1534 // C > 0xffffffffffffffff
1535 if (C1.w[1]) {
200359e8
L		 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) {
b2a00c89
L		 1545 // C > 0x9fffffffffffffff6
1546 if (C1.w[1] > 0x09 || (C1.w[1] == 0x09
1547 && C1.w[0] > 0xfffffffffffffff6ull)) {
200359e8
L		 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'
b2a00c89
L		 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])) {
200359e8
L		 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 }
b2a00c89
L		 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)
200359e8 1578 res = 0x0000000000000000ull;
b2a00c89
L		 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;
200359e8 1590 BID_RETURN (res);
b2a00c89
L		 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];
200359e8 1621 }
b2a00c89
L		 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];
200359e8
L		 1686 }
1687 }
b2a00c89
L		 1688}
1689
1690BID_RETURN (res);
200359e8
L		 1691}
1692
1693/*****************************************************************************
b2a00c89 1694 * BID128_to_uint64_xceil
200359e8
L		 1695 ****************************************************************************/
1696
b2a00c89
L		 1697BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
1698 bid128_to_uint64_xceil, x)
200359e8 1699
b2a00c89
L		 1700 UINT64 res;
1701 UINT64 x_sign;
1702 UINT64 x_exp;
1703 int exp; // unbiased exponent
200359e8 1704 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 1712
1713 // unpack x
b2a00c89
L		 1714x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1715x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
1716C1.w[1] = x.w[1] & MASK_COEFF;
1717C1.w[0] = x.w[0];
200359e8
L		 1718
1719 // check for NaN or Infinity
b2a00c89 1720if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 1721 // x is special
b2a00c89
L		 1722if ((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;
200359e8 1746 }
b2a00c89
L		 1747 BID_RETURN (res);
1748}
1749}
200359e8 1750 // check for non-canonical values (after the check for special values)
b2a00c89
L		 1751if ((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
200359e8
L		 1774 x_nr_bits =
1775 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1776 }
b2a00c89
L		 1777 } else { // if x < 2^53
1778 tmp1.d = (double) C1.w[0]; // exact conversion
200359e8 1779 x_nr_bits =
b2a00c89 1780 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200359e8 1781 }
b2a00c89
L		 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) {
200359e8
L		 1816 // set invalid flag
1817 *pfpsf |= INVALID_EXCEPTION;
1818 // return Integer Indefinite
1819 res = 0x8000000000000000ull;
1820 BID_RETURN (res);
1821 }
b2a00c89
L		 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
200359e8
L		 1829 // set invalid flag
1830 *pfpsf |= INVALID_EXCEPTION;
1831 // return Integer Indefinite
1832 res = 0x8000000000000000ull;
1833 BID_RETURN (res);
1834 }
b2a00c89
L		 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);
200359e8 1851 }
b2a00c89
L		 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);
200359e8 1864 }
b2a00c89
L		 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);
200359e8 1875 }
b2a00c89
L		 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'
200359e8 1904 }
200359e8
L		 1905 }
1906 }
b2a00c89
L		 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;
200359e8 1917 BID_RETURN (res);
b2a00c89
L		 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
200359e8
L		 1922 // set invalid flag
1923 *pfpsf |= INVALID_EXCEPTION;
1924 // return Integer Indefinite
1925 res = 0x8000000000000000ull;
1926 BID_RETURN (res);
b2a00c89
L		 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]++;
200359e8 1988 }
b2a00c89
L		 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]++;
200359e8 2000 }
b2a00c89
L		 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]++;
200359e8 2013 }
b2a00c89
L		 2014 // set the inexact flag
2015 *pfpsf |= INEXACT_EXCEPTION;
2016 } // else the result is exact
200359e8 2017 }
200359e8 2018
b2a00c89
L		 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];
200359e8
L		 2028 }
2029 }
b2a00c89
L		 2030}
2031
2032BID_RETURN (res);
200359e8
L		 2033}
2034
2035/*****************************************************************************
2036 * BID128_to_uint64_int
2037 ****************************************************************************/
2038
b2a00c89
L		 2039BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64, bid128_to_uint64_int,
2040 x)
200359e8 2041
b2a00c89
L		 2042 UINT64 res;
2043 UINT64 x_sign;
2044 UINT64 x_exp;
2045 int exp; // unbiased exponent
200359e8 2046 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 2053
2054 // unpack x
b2a00c89
L		 2055x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
2056x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
2057C1.w[1] = x.w[1] & MASK_COEFF;
2058C1.w[0] = x.w[0];
200359e8
L		 2059
2060 // check for NaN or Infinity
b2a00c89 2061if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 2062 // x is special
b2a00c89
L		 2063if ((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;
200359e8 2087 }
b2a00c89
L		 2088 BID_RETURN (res);
2089}
2090}
200359e8 2091 // check for non-canonical values (after the check for special values)
b2a00c89
L		 2092if ((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
200359e8
L		 2115 x_nr_bits =
2116 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2117 }
b2a00c89
L		 2118 } else { // if x < 2^53
2119 tmp1.d = (double) C1.w[0]; // exact conversion
200359e8 2120 x_nr_bits =
b2a00c89 2121 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200359e8 2122 }
b2a00c89
L		 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;
200359e8 2137
b2a00c89
L		 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) {
200359e8
L		 2158 // set invalid flag
2159 *pfpsf |= INVALID_EXCEPTION;
2160 // return Integer Indefinite
2161 res = 0x8000000000000000ull;
2162 BID_RETURN (res);
2163 }
b2a00c89
L		 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
200359e8
L		 2171 // set invalid flag
2172 *pfpsf |= INVALID_EXCEPTION;
2173 // return Integer Indefinite
2174 res = 0x8000000000000000ull;
2175 BID_RETURN (res);
2176 }
b2a00c89
L		 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);
200359e8 2193 }
b2a00c89
L		 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);
200359e8 2206 }
b2a00c89
L		 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'
200359e8
L		 2248 }
2249 }
2250 }
b2a00c89
L		 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
2317BID_RETURN (res);
200359e8
L		 2318}
2319
2320/*****************************************************************************
2321 * BID128_to_uint64_xint
2322 ****************************************************************************/
2323
b2a00c89
L		 2324BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64, bid128_to_uint64_xint,
2325 x)
200359e8 2326
b2a00c89
L		 2327 UINT64 res;
2328 UINT64 x_sign;
2329 UINT64 x_exp;
2330 int exp; // unbiased exponent
200359e8 2331 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 2339
2340 // unpack x
b2a00c89
L		 2341x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
2342x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
2343C1.w[1] = x.w[1] & MASK_COEFF;
2344C1.w[0] = x.w[0];
200359e8
L		 2345
2346 // check for NaN or Infinity
b2a00c89 2347if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 2348 // x is special
b2a00c89
L		 2349if ((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;
200359e8 2360 }
b2a00c89
L		 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}
200359e8 2377 // check for non-canonical values (after the check for special values)
b2a00c89
L		 2378if ((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
200359e8
L		 2401 x_nr_bits =
2402 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2403 }
b2a00c89
L		 2404 } else { // if x < 2^53
2405 tmp1.d = (double) C1.w[0]; // exact conversion
200359e8 2406 x_nr_bits =
b2a00c89 2407 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200359e8 2408 }
b2a00c89
L		 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) {
200359e8
L		 2443 // set invalid flag
2444 *pfpsf |= INVALID_EXCEPTION;
2445 // return Integer Indefinite
2446 res = 0x8000000000000000ull;
2447 BID_RETURN (res);
2448 }
b2a00c89
L		 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
200359e8
L		 2456 // set invalid flag
2457 *pfpsf |= INVALID_EXCEPTION;
2458 // return Integer Indefinite
2459 res = 0x8000000000000000ull;
2460 BID_RETURN (res);
2461 }
b2a00c89
L		 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);
200359e8 2491 }
b2a00c89
L		 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);
200359e8 2503 }
b2a00c89
L		 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);
200359e8 2515 }
b2a00c89
L		 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
200359e8 2598 }
b2a00c89
L		 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];
200359e8
L		 2637 }
2638 }
b2a00c89
L		 2639}
2640
2641BID_RETURN (res);
200359e8
L		 2642}
2643
2644/*****************************************************************************
2645 * BID128_to_uint64_rninta
2646 ****************************************************************************/
2647
b2a00c89
L		 2648BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
2649 bid128_to_uint64_rninta, x)
200359e8 2650
b2a00c89
L		 2651 UINT64 res;
2652 UINT64 x_sign;
2653 UINT64 x_exp;
2654 int exp; // unbiased exponent
200359e8 2655 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 2663
2664 // unpack x
b2a00c89
L		 2665x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
2666x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
2667C1.w[1] = x.w[1] & MASK_COEFF;
2668C1.w[0] = x.w[0];
200359e8
L		 2669
2670 // check for NaN or Infinity
b2a00c89 2671if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 2672 // x is special
b2a00c89
L		 2673if ((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)
2702if ((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);
200359e8 2727 }
b2a00c89
L		 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
200359e8
L		 2780 // set invalid flag
2781 *pfpsf |= INVALID_EXCEPTION;
2782 // return Integer Indefinite
2783 res = 0x8000000000000000ull;
b2a00c89 2784 BID_RETURN (res);
200359e8 2785 }
b2a00c89
L		 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);
200359e8 2802 }
b2a00c89
L		 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)) {
200359e8
L		 2810 // set invalid flag
2811 *pfpsf |= INVALID_EXCEPTION;
2812 // return Integer Indefinite
2813 res = 0x8000000000000000ull;
2814 BID_RETURN (res);
2815 }
b2a00c89
L		 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);
200359e8 2831 }
b2a00c89
L		 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)) {
200359e8
L		 2838 // set invalid flag
2839 *pfpsf |= INVALID_EXCEPTION;
2840 // return Integer Indefinite
2841 res = 0x8000000000000000ull;
2842 BID_RETURN (res);
2843 }
b2a00c89
L		 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])) {
200359e8
L		 2853 // set invalid flag
2854 *pfpsf |= INVALID_EXCEPTION;
2855 // return Integer Indefinite
2856 res = 0x8000000000000000ull;
2857 BID_RETURN (res);
2858 }
b2a00c89
L		 2859 // else cases that can be rounded to a 64-bit int fall through
2860 // to '1 <= q + exp <= 20'
200359e8 2861 }
b2a00c89
L		 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 {
200359e8
L		 2884 // set invalid flag
2885 *pfpsf |= INVALID_EXCEPTION;
2886 // return Integer Indefinite
2887 res = 0x8000000000000000ull;
2888 BID_RETURN (res);
2889 }
b2a00c89
L		 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];
200359e8 2924 } else {
b2a00c89
L		 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];
200359e8 2945 }
b2a00c89
L		 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];
200359e8
L		 2978 }
2979 }
b2a00c89
L		 2980}
2981
2982BID_RETURN (res);
200359e8
L		 2983}
2984
2985/*****************************************************************************
2986 * BID128_to_uint64_xrninta
2987 ****************************************************************************/
2988
b2a00c89
L		 2989BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (UINT64,
2990 bid128_to_uint64_xrninta, x)
200359e8 2991
b2a00c89
L		 2992 UINT64 res;
2993 UINT64 x_sign;
2994 UINT64 x_exp;
2995 int exp; // unbiased exponent
200359e8 2996 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
b2a00c89
L		 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;
200359e8
L		 3005
3006 // unpack x
b2a00c89
L		 3007x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
3008x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
3009C1.w[1] = x.w[1] & MASK_COEFF;
3010C1.w[0] = x.w[0];
200359e8
L		 3011
3012 // check for NaN or Infinity
b2a00c89 3013if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 3014 // x is special
b2a00c89
L		 3015if ((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;
200359e8 3039 }
b2a00c89
L		 3040 BID_RETURN (res);
3041}
3042}
200359e8 3043 // check for non-canonical values (after the check for special values)
b2a00c89
L		 3044if ((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
200359e8
L		 3067 x_nr_bits =
3068 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3069 }
b2a00c89
L		 3070 } else { // if x < 2^53
3071 tmp1.d = (double) C1.w[0]; // exact conversion
200359e8 3072 x_nr_bits =
b2a00c89 3073 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200359e8 3074 }
b2a00c89
L		 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;
200359e8 3089
b2a00c89
L		 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) {
200359e8
L		 3110 // set invalid flag
3111 *pfpsf |= INVALID_EXCEPTION;
3112 // return Integer Indefinite
3113 res = 0x8000000000000000ull;
3114 BID_RETURN (res);
3115 }
b2a00c89
L		 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);
200359e8 3128 }
b2a00c89
L		 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
200359e8 3156 res = 0x8000000000000000ull;
b2a00c89
L		 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)) {
200359e8
L		 3169 // set invalid flag
3170 *pfpsf |= INVALID_EXCEPTION;
3171 // return Integer Indefinite
3172 res = 0x8000000000000000ull;
3173 BID_RETURN (res);
3174 }
b2a00c89
L		 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
200359e8 3184 res = 0x8000000000000000ull;
b2a00c89
L		 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
200359e8
L		 3197 *pfpsf |= INVALID_EXCEPTION;
3198 // return Integer Indefinite
3199 res = 0x8000000000000000ull;
3200 BID_RETURN (res);
3201 }
b2a00c89
L		 3202 // else cases that can be rounded to a 64-bit int fall through
3203 // to '1 <= q + exp <= 20'
200359e8 3204 }
b2a00c89
L		 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;
200359e8
L		 3230 // set invalid flag
3231 *pfpsf |= INVALID_EXCEPTION;
3232 // return Integer Indefinite
3233 res = 0x8000000000000000ull;
3234 BID_RETURN (res);
3235 }
b2a00c89
L		 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])) {
200359e8
L		 3337 // set the inexact flag
3338 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L		 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])) {
200359e8
L		 3359 // set the inexact flag
3360 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L		 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])) {
200359e8
L		 3377 // set the inexact flag
3378 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L		 3379 } // else the result is exact
3380 } else { // the result is inexact; f2* <= 1/2
3381 // set the inexact flag
3382 *pfpsf |= INEXACT_EXCEPTION;
200359e8 3383 }
200359e8 3384 }
b2a00c89
L		 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];
200359e8
L		 3396 }
3397 }
b2a00c89
L		 3398}
3399
3400BID_RETURN (res);
200359e8 3401}
Mirror of https://gcc.gnu.org/git/gcc.git