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