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