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