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