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