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