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