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