1/* Copyright (C) 2007-2022 Free Software Foundation, Inc. 2 3This file is part of GCC. 4 5GCC is free software; you can redistribute it and/or modify it under 6the terms of the GNU General Public License as published by the Free 7Software Foundation; either version 3, or (at your option) any later 8version. 9 10GCC is distributed in the hope that it will be useful, but WITHOUT ANY 11WARRANTY; without even the implied warranty of MERCHANTABILITY or 12FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 13for more details. 14 15Under Section 7 of GPL version 3, you are granted additional 16permissions described in the GCC Runtime Library Exception, version 173.1, as published by the Free Software Foundation. 18 19You should have received a copy of the GNU General Public License and 20a copy of the GCC Runtime Library Exception along with this program; 21see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 22<http://www.gnu.org/licenses/>. */ 23 24#include "bid_internal.h" 25 26/***************************************************************************** 27 * BID128_to_int32_rnint 28 ****************************************************************************/ 29 30BID128_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 47x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 48x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions 49C1.w[1] = x.w[1] & MASK_COEFF; 50C1.w[0] = x.w[0]; 51 52 // check for NaN or Infinity 53if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 54 // x is special 55if ((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) 84if ((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 343BID_RETURN (res); 344} 345 346/***************************************************************************** 347 * BID128_to_int32_xrnint 348 ****************************************************************************/ 349 350BID128_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 368x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 369x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions 370C1.w[1] = x.w[1] & MASK_COEFF; 371C1.w[0] = x.w[0]; 372 373 // check for NaN or Infinity 374if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 375 // x is special 376if ((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) 405if ((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 727BID_RETURN (res); 728} 729 730/***************************************************************************** 731 * BID128_to_int32_floor 732 ****************************************************************************/ 733 734BID128_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 755x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 756x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions 757C1.w[1] = x.w[1] & MASK_COEFF; 758C1.w[0] = x.w[0]; 759 760 // check for NaN or Infinity 761if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 762 // x is special 763if ((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) 792if ((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 1100BID_RETURN (res); 1101} 1102 1103 1104/***************************************************************************** 1105 * BID128_to_int32_xfloor 1106 ****************************************************************************/ 1107 1108BID128_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 1130x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 1131x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions 1132C1.w[1] = x.w[1] & MASK_COEFF; 1133C1.w[0] = x.w[0]; 1134 1135 // check for NaN or Infinity 1136if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 1137 // x is special 1138if ((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) 1167if ((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 1489BID_RETURN (res); 1490} 1491 1492/***************************************************************************** 1493 * BID128_to_int32_ceil 1494 ****************************************************************************/ 1495 1496BID128_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 1517x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 1518x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions 1519C1.w[1] = x.w[1] & MASK_COEFF; 1520C1.w[0] = x.w[0]; 1521 1522 // check for NaN or Infinity 1523if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 1524 // x is special 1525if ((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) 1554if ((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 1862BID_RETURN (res); 1863} 1864 1865/***************************************************************************** 1866 * BID128_to_int32_xceil 1867 ****************************************************************************/ 1868 1869BID128_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 1890x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 1891x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions 1892C1.w[1] = x.w[1] & MASK_COEFF; 1893C1.w[0] = x.w[0]; 1894 1895 // check for NaN or Infinity 1896if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 1897 // x is special 1898if ((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) 1927if ((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 2249BID_RETURN (res); 2250} 2251 2252/***************************************************************************** 2253 * BID128_to_int32_int 2254 ****************************************************************************/ 2255 2256BID128_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 2275x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 2276x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions 2277C1.w[1] = x.w[1] & MASK_COEFF; 2278C1.w[0] = x.w[0]; 2279 2280 // check for NaN or Infinity 2281if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 2282 // x is special 2283if ((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) 2312if ((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 2609BID_RETURN (res); 2610} 2611 2612/***************************************************************************** 2613 * BID128_to_int32_xint 2614 ****************************************************************************/ 2615 2616BID128_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 2635x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 2636x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions 2637C1.w[1] = x.w[1] & MASK_COEFF; 2638C1.w[0] = x.w[0]; 2639 2640 // check for NaN or Infinity 2641if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 2642 // x is special 2643if ((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) 2672if ((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 2984BID_RETURN (res); 2985} 2986 2987/***************************************************************************** 2988 * BID128_to_int32_rninta 2989 ****************************************************************************/ 2990 2991BID128_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 3008x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 3009x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions 3010C1.w[1] = x.w[1] & MASK_COEFF; 3011C1.w[0] = x.w[0]; 3012 3013 // check for NaN or Infinity 3014if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 3015 // x is special 3016if ((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) 3045if ((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 3285BID_RETURN (res); 3286} 3287 3288/***************************************************************************** 3289 * BID128_to_int32_xrninta 3290 ****************************************************************************/ 3291 3292BID128_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 3310x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 3311x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions 3312C1.w[1] = x.w[1] & MASK_COEFF; 3313C1.w[0] = x.w[0]; 3314 3315 // check for NaN or Infinity 3316if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 3317 // x is special 3318if ((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) 3347if ((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 3658BID_RETURN (res); 3659} 3660