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