1 /* Copyright (C) 2007-2024 Free Software Foundation, Inc. 2 3 This file is part of GCC. 4 5 GCC is free software; you can redistribute it and/or modify it under 6 the terms of the GNU General Public License as published by the Free 7 Software Foundation; either version 3, or (at your option) any later 8 version. 9 10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY 11 WARRANTY; without even the implied warranty of MERCHANTABILITY or 12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 13 for more details. 14 15 Under Section 7 of GPL version 3, you are granted additional 16 permissions described in the GCC Runtime Library Exception, version 17 3.1, as published by the Free Software Foundation. 18 19 You should have received a copy of the GNU General Public License and 20 a copy of the GCC Runtime Library Exception along with this program; 21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 22 <http://www.gnu.org/licenses/>. */ 23 24 #include "bid_internal.h" 25 26 /***************************************************************************** 27 * BID64_to_uint32_rnint 28 ****************************************************************************/ 29 30 #if DECIMAL_CALL_BY_REFERENCE 31 void 32 bid64_to_uint32_rnint (unsigned int *pres, UINT64 * px 33 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 34 _EXC_INFO_PARAM) { 35 UINT64 x = *px; 36 #else 37 unsigned int 38 bid64_to_uint32_rnint (UINT64 x 39 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 40 _EXC_INFO_PARAM) { 41 #endif 42 unsigned int res; 43 UINT64 x_sign; 44 UINT64 x_exp; 45 int exp; // unbiased exponent 46 // Note: C1 represents x_significand (UINT64) 47 UINT64 tmp64; 48 BID_UI64DOUBLE tmp1; 49 unsigned int x_nr_bits; 50 int q, ind, shift; 51 UINT64 C1; 52 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits 53 UINT128 fstar; 54 UINT128 P128; 55 56 // check for NaN or Infinity 57 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) { 58 // set invalid flag 59 *pfpsf |= INVALID_EXCEPTION; 60 // return Integer Indefinite 61 res = 0x80000000; 62 BID_RETURN (res); 63 } 64 // unpack x 65 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 66 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] => 67 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 68 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased 69 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 70 if (C1 > 9999999999999999ull) { // non-canonical 71 x_exp = 0; 72 C1 = 0; 73 } 74 } else { 75 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased 76 C1 = x & MASK_BINARY_SIG1; 77 } 78 79 // check for zeros (possibly from non-canonical values) 80 if (C1 == 0x0ull) { 81 // x is 0 82 res = 0x00000000; 83 BID_RETURN (res); 84 } 85 // x is not special and is not zero 86 87 // q = nr. of decimal digits in x (1 <= q <= 54) 88 // determine first the nr. of bits in x 89 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 90 // split the 64-bit value in two 32-bit halves to avoid rounding errors 91 if (C1 >= 0x0000000100000000ull) { // x >= 2^32 92 tmp1.d = (double) (C1 >> 32); // exact conversion 93 x_nr_bits = 94 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 95 } else { // x < 2^32 96 tmp1.d = (double) C1; // exact conversion 97 x_nr_bits = 98 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 99 } 100 } else { // if x < 2^53 101 tmp1.d = (double) C1; // exact conversion 102 x_nr_bits = 103 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 104 } 105 q = nr_digits[x_nr_bits - 1].digits; 106 if (q == 0) { 107 q = nr_digits[x_nr_bits - 1].digits1; 108 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 109 q++; 110 } 111 exp = x_exp - 398; // unbiased exponent 112 113 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) 114 // set invalid flag 115 *pfpsf |= INVALID_EXCEPTION; 116 // return Integer Indefinite 117 res = 0x80000000; 118 BID_RETURN (res); 119 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) 120 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... 121 // so x rounded to an integer may or may not fit in an unsigned 32-bit int 122 // the cases that do not fit are identified here; the ones that fit 123 // fall through and will be handled with other cases further, 124 // under '1 <= q + exp <= 10' 125 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1/2 126 // => set invalid flag 127 *pfpsf |= INVALID_EXCEPTION; 128 // return Integer Indefinite 129 res = 0x80000000; 130 BID_RETURN (res); 131 } else { // if n > 0 and q + exp = 10 132 // if n >= 2^32 - 1/2 then n is too large 133 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2 134 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=16 135 // <=> C * 10^(11-q) >= 0x9fffffffb, 1<=q<=16 136 if (q <= 11) { 137 // Note: C * 10^(11-q) has 10 or 11 digits; 0x9fffffffb has 11 digits 138 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int 139 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) 140 if (tmp64 >= 0x9fffffffbull) { 141 // set invalid flag 142 *pfpsf |= INVALID_EXCEPTION; 143 // return Integer Indefinite 144 res = 0x80000000; 145 BID_RETURN (res); 146 } 147 // else cases that can be rounded to a 32-bit unsigned int fall through 148 // to '1 <= q + exp <= 10' 149 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2 150 // C * 10^(11-q) >= 0x9fffffffb <=> 151 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 5 152 // (scale 2^32-1/2 up) 153 // Note: 0x9fffffffb*10^(q-11) has q-1 or q digits, where q <= 16 154 tmp64 = 0x9fffffffbull * ten2k64[q - 11]; 155 if (C1 >= tmp64) { 156 // set invalid flag 157 *pfpsf |= INVALID_EXCEPTION; 158 // return Integer Indefinite 159 res = 0x80000000; 160 BID_RETURN (res); 161 } 162 // else cases that can be rounded to a 32-bit int fall through 163 // to '1 <= q + exp <= 10' 164 } 165 } 166 } 167 // n is not too large to be converted to int32 if -1/2 <= n < 2^32 - 1/2 168 // Note: some of the cases tested for above fall through to this point 169 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1) 170 // return 0 171 res = 0x00000000; 172 BID_RETURN (res); 173 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1) 174 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1) 175 // res = 0 176 // else if x > 0 177 // res = +1 178 // else // if x < 0 179 // invalid exc 180 ind = q - 1; 181 if (C1 <= midpoint64[ind]) { 182 res = 0x00000000; // return 0 183 } else if (x_sign) { // n < 0 184 // set invalid flag 185 *pfpsf |= INVALID_EXCEPTION; 186 // return Integer Indefinite 187 res = 0x80000000; 188 BID_RETURN (res); 189 } else { // n > 0 190 res = 0x00000001; // return +1 191 } 192 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9) 193 // -2^32-1/2 <= x <= -1 or 1 <= x < 2^32-1/2 so if positive, x can be 194 // rounded to nearest to a 32-bit unsigned integer 195 if (x_sign) { // x <= -1 196 // set invalid flag 197 *pfpsf |= INVALID_EXCEPTION; 198 // return Integer Indefinite 199 res = 0x80000000; 200 BID_RETURN (res); 201 } 202 // 1 <= x < 2^32-1/2 so x can be rounded 203 // to nearest to a 32-bit unsigned integer 204 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10 205 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x' 206 // chop off ind digits from the lower part of C1 207 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits 208 C1 = C1 + midpoint64[ind - 1]; 209 // calculate C* and f* 210 // C* is actually floor(C*) in this case 211 // C* and f* need shifting and masking, as shown by 212 // shiftright128[] and maskhigh128[] 213 // 1 <= x <= 15 214 // kx = 10^(-x) = ten2mk64[ind - 1] 215 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 216 // the approximation of 10^(-x) was rounded up to 54 bits 217 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); 218 Cstar = P128.w[1]; 219 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 220 fstar.w[0] = P128.w[0]; 221 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g. 222 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999 223 // if (0 < f* < 10^(-x)) then the result is a midpoint 224 // if floor(C*) is even then C* = floor(C*) - logical right 225 // shift; C* has p decimal digits, correct by Prop. 1) 226 // else if floor(C*) is odd C* = floor(C*)-1 (logical right 227 // shift; C* has p decimal digits, correct by Pr. 1) 228 // else 229 // C* = floor(C*) (logical right shift; C has p decimal digits, 230 // correct by Property 1) 231 // n = C* * 10^(e+x) 232 233 // shift right C* by Ex-64 = shiftright128[ind] 234 shift = shiftright128[ind - 1]; // 0 <= shift <= 39 235 Cstar = Cstar >> shift; 236 237 // if the result was a midpoint it was rounded away from zero, so 238 // it will need a correction 239 // check for midpoints 240 if ((fstar.w[1] == 0) && fstar.w[0] && 241 (fstar.w[0] <= ten2mk128trunc[ind - 1].w[1])) { 242 // ten2mk128trunc[ind -1].w[1] is identical to 243 // ten2mk128[ind -1].w[1] 244 // the result is a midpoint; round to nearest 245 if (Cstar & 0x01) { // Cstar is odd; MP in [EVEN, ODD] 246 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1 247 Cstar--; // Cstar is now even 248 } // else MP in [ODD, EVEN] 249 } 250 res = Cstar; // the result is positive 251 } else if (exp == 0) { 252 // 1 <= q <= 10 253 // res = +C (exact) 254 res = C1; // the result is positive 255 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 256 // res = +C * 10^exp (exact) 257 res = C1 * ten2k64[exp]; // the result is positive 258 } 259 } 260 BID_RETURN (res); 261 } 262 263 /***************************************************************************** 264 * BID64_to_uint32_xrnint 265 ****************************************************************************/ 266 267 #if DECIMAL_CALL_BY_REFERENCE 268 void 269 bid64_to_uint32_xrnint (unsigned int *pres, UINT64 * px 270 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 271 _EXC_INFO_PARAM) { 272 UINT64 x = *px; 273 #else 274 unsigned int 275 bid64_to_uint32_xrnint (UINT64 x 276 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 277 _EXC_INFO_PARAM) { 278 #endif 279 unsigned int res; 280 UINT64 x_sign; 281 UINT64 x_exp; 282 int exp; // unbiased exponent 283 // Note: C1 represents x_significand (UINT64) 284 UINT64 tmp64; 285 BID_UI64DOUBLE tmp1; 286 unsigned int x_nr_bits; 287 int q, ind, shift; 288 UINT64 C1; 289 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits 290 UINT128 fstar; 291 UINT128 P128; 292 293 // check for NaN or Infinity 294 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) { 295 // set invalid flag 296 *pfpsf |= INVALID_EXCEPTION; 297 // return Integer Indefinite 298 res = 0x80000000; 299 BID_RETURN (res); 300 } 301 // unpack x 302 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 303 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] => 304 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 305 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased 306 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 307 if (C1 > 9999999999999999ull) { // non-canonical 308 x_exp = 0; 309 C1 = 0; 310 } 311 } else { 312 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased 313 C1 = x & MASK_BINARY_SIG1; 314 } 315 316 // check for zeros (possibly from non-canonical values) 317 if (C1 == 0x0ull) { 318 // x is 0 319 res = 0x00000000; 320 BID_RETURN (res); 321 } 322 // x is not special and is not zero 323 324 // q = nr. of decimal digits in x (1 <= q <= 54) 325 // determine first the nr. of bits in x 326 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 327 // split the 64-bit value in two 32-bit halves to avoid rounding errors 328 if (C1 >= 0x0000000100000000ull) { // x >= 2^32 329 tmp1.d = (double) (C1 >> 32); // exact conversion 330 x_nr_bits = 331 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 332 } else { // x < 2^32 333 tmp1.d = (double) C1; // exact conversion 334 x_nr_bits = 335 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 336 } 337 } else { // if x < 2^53 338 tmp1.d = (double) C1; // exact conversion 339 x_nr_bits = 340 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 341 } 342 q = nr_digits[x_nr_bits - 1].digits; 343 if (q == 0) { 344 q = nr_digits[x_nr_bits - 1].digits1; 345 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 346 q++; 347 } 348 exp = x_exp - 398; // unbiased exponent 349 350 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) 351 // set invalid flag 352 *pfpsf |= INVALID_EXCEPTION; 353 // return Integer Indefinite 354 res = 0x80000000; 355 BID_RETURN (res); 356 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) 357 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... 358 // so x rounded to an integer may or may not fit in an unsigned 32-bit int 359 // the cases that do not fit are identified here; the ones that fit 360 // fall through and will be handled with other cases further, 361 // under '1 <= q + exp <= 10' 362 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1/2 363 // => set invalid flag 364 *pfpsf |= INVALID_EXCEPTION; 365 // return Integer Indefinite 366 res = 0x80000000; 367 BID_RETURN (res); 368 } else { // if n > 0 and q + exp = 10 369 // if n >= 2^32 - 1/2 then n is too large 370 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2 371 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=16 372 // <=> C * 10^(11-q) >= 0x9fffffffb, 1<=q<=16 373 if (q <= 11) { 374 // Note: C * 10^(11-q) has 10 or 11 digits; 0x9fffffffb has 11 digits 375 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int 376 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) 377 if (tmp64 >= 0x9fffffffbull) { 378 // set invalid flag 379 *pfpsf |= INVALID_EXCEPTION; 380 // return Integer Indefinite 381 res = 0x80000000; 382 BID_RETURN (res); 383 } 384 // else cases that can be rounded to a 32-bit unsigned int fall through 385 // to '1 <= q + exp <= 10' 386 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2 387 // C * 10^(11-q) >= 0x9fffffffb <=> 388 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 5 389 // (scale 2^32-1/2 up) 390 // Note: 0x9fffffffb*10^(q-11) has q-1 or q digits, where q <= 16 391 tmp64 = 0x9fffffffbull * ten2k64[q - 11]; 392 if (C1 >= tmp64) { 393 // set invalid flag 394 *pfpsf |= INVALID_EXCEPTION; 395 // return Integer Indefinite 396 res = 0x80000000; 397 BID_RETURN (res); 398 } 399 // else cases that can be rounded to a 32-bit int fall through 400 // to '1 <= q + exp <= 10' 401 } 402 } 403 } 404 // n is not too large to be converted to int32 if -1/2 <= n < 2^32 - 1/2 405 // Note: some of the cases tested for above fall through to this point 406 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1) 407 // set inexact flag 408 *pfpsf |= INEXACT_EXCEPTION; 409 // return 0 410 res = 0x00000000; 411 BID_RETURN (res); 412 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1) 413 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1) 414 // res = 0 415 // else if x > 0 416 // res = +1 417 // else // if x < 0 418 // invalid exc 419 ind = q - 1; 420 if (C1 <= midpoint64[ind]) { 421 res = 0x00000000; // return 0 422 } else if (x_sign) { // n < 0 423 // set invalid flag 424 *pfpsf |= INVALID_EXCEPTION; 425 // return Integer Indefinite 426 res = 0x80000000; 427 BID_RETURN (res); 428 } else { // n > 0 429 res = 0x00000001; // return +1 430 } 431 // set inexact flag 432 *pfpsf |= INEXACT_EXCEPTION; 433 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9) 434 // -2^32-1/2 <= x <= -1 or 1 <= x < 2^32-1/2 so if positive, x can be 435 // rounded to nearest to a 32-bit unsigned integer 436 if (x_sign) { // x <= -1 437 // set invalid flag 438 *pfpsf |= INVALID_EXCEPTION; 439 // return Integer Indefinite 440 res = 0x80000000; 441 BID_RETURN (res); 442 } 443 // 1 <= x < 2^32-1/2 so x can be rounded 444 // to nearest to a 32-bit unsigned integer 445 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10 446 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x' 447 // chop off ind digits from the lower part of C1 448 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits 449 C1 = C1 + midpoint64[ind - 1]; 450 // calculate C* and f* 451 // C* is actually floor(C*) in this case 452 // C* and f* need shifting and masking, as shown by 453 // shiftright128[] and maskhigh128[] 454 // 1 <= x <= 15 455 // kx = 10^(-x) = ten2mk64[ind - 1] 456 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 457 // the approximation of 10^(-x) was rounded up to 54 bits 458 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); 459 Cstar = P128.w[1]; 460 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 461 fstar.w[0] = P128.w[0]; 462 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g. 463 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999 464 // if (0 < f* < 10^(-x)) then the result is a midpoint 465 // if floor(C*) is even then C* = floor(C*) - logical right 466 // shift; C* has p decimal digits, correct by Prop. 1) 467 // else if floor(C*) is odd C* = floor(C*)-1 (logical right 468 // shift; C* has p decimal digits, correct by Pr. 1) 469 // else 470 // C* = floor(C*) (logical right shift; C has p decimal digits, 471 // correct by Property 1) 472 // n = C* * 10^(e+x) 473 474 // shift right C* by Ex-64 = shiftright128[ind] 475 shift = shiftright128[ind - 1]; // 0 <= shift <= 39 476 Cstar = Cstar >> shift; 477 // determine inexactness of the rounding of C* 478 // if (0 < f* - 1/2 < 10^(-x)) then 479 // the result is exact 480 // else // if (f* - 1/2 > T*) then 481 // the result is inexact 482 if (ind - 1 <= 2) { // fstar.w[1] is 0 483 if (fstar.w[0] > 0x8000000000000000ull) { 484 // f* > 1/2 and the result may be exact 485 tmp64 = fstar.w[0] - 0x8000000000000000ull; // f* - 1/2 486 if ((tmp64 > ten2mk128trunc[ind - 1].w[1])) { 487 // ten2mk128trunc[ind -1].w[1] is identical to 488 // ten2mk128[ind -1].w[1] 489 // set the inexact flag 490 *pfpsf |= INEXACT_EXCEPTION; 491 } // else the result is exact 492 } else { // the result is inexact; f2* <= 1/2 493 // set the inexact flag 494 *pfpsf |= INEXACT_EXCEPTION; 495 } 496 } else { // if 3 <= ind - 1 <= 14 497 if (fstar.w[1] > onehalf128[ind - 1] || 498 (fstar.w[1] == onehalf128[ind - 1] && fstar.w[0])) { 499 // f2* > 1/2 and the result may be exact 500 // Calculate f2* - 1/2 501 tmp64 = fstar.w[1] - onehalf128[ind - 1]; 502 if (tmp64 || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) { 503 // ten2mk128trunc[ind -1].w[1] is identical to 504 // ten2mk128[ind -1].w[1] 505 // set the inexact flag 506 *pfpsf |= INEXACT_EXCEPTION; 507 } // else the result is exact 508 } else { // the result is inexact; f2* <= 1/2 509 // set the inexact flag 510 *pfpsf |= INEXACT_EXCEPTION; 511 } 512 } 513 514 // if the result was a midpoint it was rounded away from zero, so 515 // it will need a correction 516 // check for midpoints 517 if ((fstar.w[1] == 0) && fstar.w[0] && 518 (fstar.w[0] <= ten2mk128trunc[ind - 1].w[1])) { 519 // ten2mk128trunc[ind -1].w[1] is identical to 520 // ten2mk128[ind -1].w[1] 521 // the result is a midpoint; round to nearest 522 if (Cstar & 0x01) { // Cstar is odd; MP in [EVEN, ODD] 523 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1 524 Cstar--; // Cstar is now even 525 } // else MP in [ODD, EVEN] 526 } 527 res = Cstar; // the result is positive 528 } else if (exp == 0) { 529 // 1 <= q <= 10 530 // res = +C (exact) 531 res = C1; // the result is positive 532 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 533 // res = +C * 10^exp (exact) 534 res = C1 * ten2k64[exp]; // the result is positive 535 } 536 } 537 BID_RETURN (res); 538 } 539 540 /***************************************************************************** 541 * BID64_to_uint32_floor 542 ****************************************************************************/ 543 544 #if DECIMAL_CALL_BY_REFERENCE 545 void 546 bid64_to_uint32_floor (unsigned int *pres, UINT64 * px 547 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 548 _EXC_INFO_PARAM) { 549 UINT64 x = *px; 550 #else 551 unsigned int 552 bid64_to_uint32_floor (UINT64 x 553 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 554 _EXC_INFO_PARAM) { 555 #endif 556 unsigned int res; 557 UINT64 x_sign; 558 UINT64 x_exp; 559 int exp; // unbiased exponent 560 // Note: C1 represents x_significand (UINT64) 561 UINT64 tmp64; 562 BID_UI64DOUBLE tmp1; 563 unsigned int x_nr_bits; 564 int q, ind, shift; 565 UINT64 C1; 566 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits 567 UINT128 P128; 568 569 // check for NaN or Infinity 570 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) { 571 // set invalid flag 572 *pfpsf |= INVALID_EXCEPTION; 573 // return Integer Indefinite 574 res = 0x80000000; 575 BID_RETURN (res); 576 } 577 // unpack x 578 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 579 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] => 580 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 581 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased 582 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 583 if (C1 > 9999999999999999ull) { // non-canonical 584 x_exp = 0; 585 C1 = 0; 586 } 587 } else { 588 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased 589 C1 = x & MASK_BINARY_SIG1; 590 } 591 592 // check for zeros (possibly from non-canonical values) 593 if (C1 == 0x0ull) { 594 // x is 0 595 res = 0x00000000; 596 BID_RETURN (res); 597 } 598 // x is not special and is not zero 599 600 if (x_sign) { // if n < 0 the conversion is invalid 601 // set invalid flag 602 *pfpsf |= INVALID_EXCEPTION; 603 // return Integer Indefinite 604 res = 0x80000000; 605 BID_RETURN (res); 606 } 607 // q = nr. of decimal digits in x (1 <= q <= 54) 608 // determine first the nr. of bits in x 609 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 610 // split the 64-bit value in two 32-bit halves to avoid rounding errors 611 if (C1 >= 0x0000000100000000ull) { // x >= 2^32 612 tmp1.d = (double) (C1 >> 32); // exact conversion 613 x_nr_bits = 614 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 615 } else { // x < 2^32 616 tmp1.d = (double) C1; // exact conversion 617 x_nr_bits = 618 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 619 } 620 } else { // if x < 2^53 621 tmp1.d = (double) C1; // exact conversion 622 x_nr_bits = 623 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 624 } 625 q = nr_digits[x_nr_bits - 1].digits; 626 if (q == 0) { 627 q = nr_digits[x_nr_bits - 1].digits1; 628 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 629 q++; 630 } 631 exp = x_exp - 398; // unbiased exponent 632 633 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) 634 // set invalid flag 635 *pfpsf |= INVALID_EXCEPTION; 636 // return Integer Indefinite 637 res = 0x80000000; 638 BID_RETURN (res); 639 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) 640 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... 641 // so x rounded to an integer may or may not fit in an unsigned 32-bit int 642 // the cases that do not fit are identified here; the ones that fit 643 // fall through and will be handled with other cases further, 644 // under '1 <= q + exp <= 10' 645 // n > 0 and q + exp = 10 646 // if n >= 2^32 then n is too large 647 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32 648 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=16 649 // <=> C * 10^(11-q) >= 0xa00000000, 1<=q<=16 650 if (q <= 11) { 651 // Note: C * 10^(11-q) has 10 or 11 digits; 0xa00000000 has 11 digits 652 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int 653 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) 654 if (tmp64 >= 0xa00000000ull) { 655 // set invalid flag 656 *pfpsf |= INVALID_EXCEPTION; 657 // return Integer Indefinite 658 res = 0x80000000; 659 BID_RETURN (res); 660 } 661 // else cases that can be rounded to a 32-bit unsigned int fall through 662 // to '1 <= q + exp <= 10' 663 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2 664 // C * 10^(11-q) >= 0xa00000000 <=> 665 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 5 666 // (scale 2^32-1/2 up) 667 // Note: 0xa00000000*10^(q-11) has q-1 or q digits, where q <= 16 668 tmp64 = 0xa00000000ull * ten2k64[q - 11]; 669 if (C1 >= tmp64) { 670 // set invalid flag 671 *pfpsf |= INVALID_EXCEPTION; 672 // return Integer Indefinite 673 res = 0x80000000; 674 BID_RETURN (res); 675 } 676 // else cases that can be rounded to a 32-bit int fall through 677 // to '1 <= q + exp <= 10' 678 } 679 } 680 // n is not too large to be converted to int32 if -1 < n < 2^32 681 // Note: some of the cases tested for above fall through to this point 682 if ((q + exp) <= 0) { // n = +0.[0...0]c(0)c(1)...c(q-1) 683 // return 0 684 res = 0x00000000; 685 BID_RETURN (res); 686 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9) 687 // 1 <= x < 2^32 so x can be rounded 688 // to nearest to a 32-bit unsigned integer 689 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10 690 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x' 691 // chop off ind digits from the lower part of C1 692 // C1 fits in 64 bits 693 // calculate C* and f* 694 // C* is actually floor(C*) in this case 695 // C* and f* need shifting and masking, as shown by 696 // shiftright128[] and maskhigh128[] 697 // 1 <= x <= 15 698 // kx = 10^(-x) = ten2mk64[ind - 1] 699 // C* = C1 * 10^(-x) 700 // the approximation of 10^(-x) was rounded up to 54 bits 701 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); 702 Cstar = P128.w[1]; 703 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g. 704 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999 705 // C* = floor(C*) (logical right shift; C has p decimal digits, 706 // correct by Property 1) 707 // n = C* * 10^(e+x) 708 709 // shift right C* by Ex-64 = shiftright128[ind] 710 shift = shiftright128[ind - 1]; // 0 <= shift <= 39 711 Cstar = Cstar >> shift; 712 713 res = Cstar; // the result is positive 714 } else if (exp == 0) { 715 // 1 <= q <= 10 716 // res = +C (exact) 717 res = C1; // the result is positive 718 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 719 // res = +C * 10^exp (exact) 720 res = C1 * ten2k64[exp]; // the result is positive 721 } 722 } 723 BID_RETURN (res); 724 } 725 726 /***************************************************************************** 727 * BID64_to_uint32_xfloor 728 ****************************************************************************/ 729 730 #if DECIMAL_CALL_BY_REFERENCE 731 void 732 bid64_to_uint32_xfloor (unsigned int *pres, UINT64 * px 733 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 734 _EXC_INFO_PARAM) { 735 UINT64 x = *px; 736 #else 737 unsigned int 738 bid64_to_uint32_xfloor (UINT64 x 739 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 740 _EXC_INFO_PARAM) { 741 #endif 742 unsigned int res; 743 UINT64 x_sign; 744 UINT64 x_exp; 745 int exp; // unbiased exponent 746 // Note: C1 represents x_significand (UINT64) 747 UINT64 tmp64; 748 BID_UI64DOUBLE tmp1; 749 unsigned int x_nr_bits; 750 int q, ind, shift; 751 UINT64 C1; 752 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits 753 UINT128 fstar; 754 UINT128 P128; 755 756 // check for NaN or Infinity 757 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) { 758 // set invalid flag 759 *pfpsf |= INVALID_EXCEPTION; 760 // return Integer Indefinite 761 res = 0x80000000; 762 BID_RETURN (res); 763 } 764 // unpack x 765 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 766 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] => 767 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 768 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased 769 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 770 if (C1 > 9999999999999999ull) { // non-canonical 771 x_exp = 0; 772 C1 = 0; 773 } 774 } else { 775 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased 776 C1 = x & MASK_BINARY_SIG1; 777 } 778 779 // check for zeros (possibly from non-canonical values) 780 if (C1 == 0x0ull) { 781 // x is 0 782 res = 0x00000000; 783 BID_RETURN (res); 784 } 785 // x is not special and is not zero 786 787 if (x_sign) { // if n < 0 the conversion is invalid 788 // set invalid flag 789 *pfpsf |= INVALID_EXCEPTION; 790 // return Integer Indefinite 791 res = 0x80000000; 792 BID_RETURN (res); 793 } 794 // q = nr. of decimal digits in x (1 <= q <= 54) 795 // determine first the nr. of bits in x 796 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 797 // split the 64-bit value in two 32-bit halves to avoid rounding errors 798 if (C1 >= 0x0000000100000000ull) { // x >= 2^32 799 tmp1.d = (double) (C1 >> 32); // exact conversion 800 x_nr_bits = 801 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 802 } else { // x < 2^32 803 tmp1.d = (double) C1; // exact conversion 804 x_nr_bits = 805 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 806 } 807 } else { // if x < 2^53 808 tmp1.d = (double) C1; // exact conversion 809 x_nr_bits = 810 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 811 } 812 q = nr_digits[x_nr_bits - 1].digits; 813 if (q == 0) { 814 q = nr_digits[x_nr_bits - 1].digits1; 815 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 816 q++; 817 } 818 exp = x_exp - 398; // unbiased exponent 819 820 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) 821 // set invalid flag 822 *pfpsf |= INVALID_EXCEPTION; 823 // return Integer Indefinite 824 res = 0x80000000; 825 BID_RETURN (res); 826 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) 827 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... 828 // so x rounded to an integer may or may not fit in an unsigned 32-bit int 829 // the cases that do not fit are identified here; the ones that fit 830 // fall through and will be handled with other cases further, 831 // under '1 <= q + exp <= 10' 832 // if n > 0 and q + exp = 10 833 // if n >= 2^32 then n is too large 834 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32 835 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=16 836 // <=> C * 10^(11-q) >= 0xa00000000, 1<=q<=16 837 if (q <= 11) { 838 // Note: C * 10^(11-q) has 10 or 11 digits; 0xa00000000 has 11 digits 839 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int 840 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) 841 if (tmp64 >= 0xa00000000ull) { 842 // set invalid flag 843 *pfpsf |= INVALID_EXCEPTION; 844 // return Integer Indefinite 845 res = 0x80000000; 846 BID_RETURN (res); 847 } 848 // else cases that can be rounded to a 32-bit unsigned int fall through 849 // to '1 <= q + exp <= 10' 850 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2 851 // C * 10^(11-q) >= 0xa00000000 <=> 852 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 5 853 // (scale 2^32-1/2 up) 854 // Note: 0xa00000000*10^(q-11) has q-1 or q digits, where q <= 16 855 tmp64 = 0xa00000000ull * ten2k64[q - 11]; 856 if (C1 >= tmp64) { 857 // set invalid flag 858 *pfpsf |= INVALID_EXCEPTION; 859 // return Integer Indefinite 860 res = 0x80000000; 861 BID_RETURN (res); 862 } 863 // else cases that can be rounded to a 32-bit int fall through 864 // to '1 <= q + exp <= 10' 865 } 866 } 867 // n is not too large to be converted to int32 if -1 < n < 2^32 868 // Note: some of the cases tested for above fall through to this point 869 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1) 870 // set inexact flag 871 *pfpsf |= INEXACT_EXCEPTION; 872 // return 0 873 res = 0x00000000; 874 BID_RETURN (res); 875 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9) 876 // 1 <= x < 2^32 so x can be rounded 877 // to nearest to a 32-bit unsigned integer 878 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10 879 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x' 880 // chop off ind digits from the lower part of C1 881 // C1 fits in 64 bits 882 // calculate C* and f* 883 // C* is actually floor(C*) in this case 884 // C* and f* need shifting and masking, as shown by 885 // shiftright128[] and maskhigh128[] 886 // 1 <= x <= 15 887 // kx = 10^(-x) = ten2mk64[ind - 1] 888 // C* = C1 * 10^(-x) 889 // the approximation of 10^(-x) was rounded up to 54 bits 890 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); 891 Cstar = P128.w[1]; 892 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 893 fstar.w[0] = P128.w[0]; 894 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g. 895 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999 896 // C* = floor(C*) (logical right shift; C has p decimal digits, 897 // correct by Property 1) 898 // n = C* * 10^(e+x) 899 900 // shift right C* by Ex-64 = shiftright128[ind] 901 shift = shiftright128[ind - 1]; // 0 <= shift <= 39 902 Cstar = Cstar >> shift; 903 // determine inexactness of the rounding of C* 904 // if (0 < f* < 10^(-x)) then 905 // the result is exact 906 // else // if (f* > T*) then 907 // the result is inexact 908 if (ind - 1 <= 2) { 909 if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) { 910 // ten2mk128trunc[ind -1].w[1] is identical to 911 // ten2mk128[ind -1].w[1] 912 // set the inexact flag 913 *pfpsf |= INEXACT_EXCEPTION; 914 } // else the result is exact 915 } else { // if 3 <= ind - 1 <= 14 916 if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) { 917 // ten2mk128trunc[ind -1].w[1] is identical to 918 // ten2mk128[ind -1].w[1] 919 // set the inexact flag 920 *pfpsf |= INEXACT_EXCEPTION; 921 } // else the result is exact 922 } 923 924 res = Cstar; // the result is positive 925 } else if (exp == 0) { 926 // 1 <= q <= 10 927 // res = +C (exact) 928 res = C1; // the result is positive 929 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 930 // res = +C * 10^exp (exact) 931 res = C1 * ten2k64[exp]; // the result is positive 932 } 933 } 934 BID_RETURN (res); 935 } 936 937 /***************************************************************************** 938 * BID64_to_uint32_ceil 939 ****************************************************************************/ 940 941 #if DECIMAL_CALL_BY_REFERENCE 942 void 943 bid64_to_uint32_ceil (unsigned int *pres, UINT64 * px 944 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 945 _EXC_INFO_PARAM) { 946 UINT64 x = *px; 947 #else 948 unsigned int 949 bid64_to_uint32_ceil (UINT64 x 950 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 951 _EXC_INFO_PARAM) { 952 #endif 953 unsigned int res; 954 UINT64 x_sign; 955 UINT64 x_exp; 956 int exp; // unbiased exponent 957 // Note: C1 represents x_significand (UINT64) 958 UINT64 tmp64; 959 BID_UI64DOUBLE tmp1; 960 unsigned int x_nr_bits; 961 int q, ind, shift; 962 UINT64 C1; 963 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits 964 UINT128 fstar; 965 UINT128 P128; 966 967 // check for NaN or Infinity 968 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) { 969 // set invalid flag 970 *pfpsf |= INVALID_EXCEPTION; 971 // return Integer Indefinite 972 res = 0x80000000; 973 BID_RETURN (res); 974 } 975 // unpack x 976 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 977 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] => 978 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 979 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased 980 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 981 if (C1 > 9999999999999999ull) { // non-canonical 982 x_exp = 0; 983 C1 = 0; 984 } 985 } else { 986 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased 987 C1 = x & MASK_BINARY_SIG1; 988 } 989 990 // check for zeros (possibly from non-canonical values) 991 if (C1 == 0x0ull) { 992 // x is 0 993 res = 0x00000000; 994 BID_RETURN (res); 995 } 996 // x is not special and is not zero 997 998 // q = nr. of decimal digits in x (1 <= q <= 54) 999 // determine first the nr. of bits in x 1000 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 1001 // split the 64-bit value in two 32-bit halves to avoid rounding errors 1002 if (C1 >= 0x0000000100000000ull) { // x >= 2^32 1003 tmp1.d = (double) (C1 >> 32); // exact conversion 1004 x_nr_bits = 1005 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1006 } else { // x < 2^32 1007 tmp1.d = (double) C1; // exact conversion 1008 x_nr_bits = 1009 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1010 } 1011 } else { // if x < 2^53 1012 tmp1.d = (double) C1; // exact conversion 1013 x_nr_bits = 1014 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1015 } 1016 q = nr_digits[x_nr_bits - 1].digits; 1017 if (q == 0) { 1018 q = nr_digits[x_nr_bits - 1].digits1; 1019 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 1020 q++; 1021 } 1022 exp = x_exp - 398; // unbiased exponent 1023 1024 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) 1025 // set invalid flag 1026 *pfpsf |= INVALID_EXCEPTION; 1027 // return Integer Indefinite 1028 res = 0x80000000; 1029 BID_RETURN (res); 1030 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) 1031 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... 1032 // so x rounded to an integer may or may not fit in an unsigned 32-bit int 1033 // the cases that do not fit are identified here; the ones that fit 1034 // fall through and will be handled with other cases further, 1035 // under '1 <= q + exp <= 10' 1036 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1 1037 // => set invalid flag 1038 *pfpsf |= INVALID_EXCEPTION; 1039 // return Integer Indefinite 1040 res = 0x80000000; 1041 BID_RETURN (res); 1042 } else { // if n > 0 and q + exp = 10 1043 // if n > 2^32 - 1 then n is too large 1044 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^32 - 1 1045 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6, 1<=q<=16 1046 // <=> C * 10^(11-q) > 0x9fffffff6, 1<=q<=16 1047 if (q <= 11) { 1048 // Note: C * 10^(11-q) has 10 or 11 digits; 0x9fffffff6 has 11 digits 1049 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int 1050 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) 1051 if (tmp64 > 0x9fffffff6ull) { 1052 // set invalid flag 1053 *pfpsf |= INVALID_EXCEPTION; 1054 // return Integer Indefinite 1055 res = 0x80000000; 1056 BID_RETURN (res); 1057 } 1058 // else cases that can be rounded to a 32-bit unsigned int fall through 1059 // to '1 <= q + exp <= 10' 1060 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2 1061 // C * 10^(11-q) > 0x9fffffff6 <=> 1062 // C > 0x9fffffff6 * 10^(q-11) where 1 <= q - 11 <= 5 1063 // (scale 2^32-1 up) 1064 // Note: 0x9fffffff6*10^(q-11) has q-1 or q digits, where q <= 16 1065 tmp64 = 0x9fffffff6ull * ten2k64[q - 11]; 1066 if (C1 > tmp64) { 1067 // set invalid flag 1068 *pfpsf |= INVALID_EXCEPTION; 1069 // return Integer Indefinite 1070 res = 0x80000000; 1071 BID_RETURN (res); 1072 } 1073 // else cases that can be rounded to a 32-bit int fall through 1074 // to '1 <= q + exp <= 10' 1075 } 1076 } 1077 } 1078 // n is not too large to be converted to int32 if -1 < n < 2^32 1079 // Note: some of the cases tested for above fall through to this point 1080 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1) 1081 // return 0 or 1 1082 if (x_sign) 1083 res = 0x00000000; 1084 else 1085 res = 0x00000001; 1086 BID_RETURN (res); 1087 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9) 1088 // x <= -1 or 1 <= x <= 2^32 - 1 so if positive, x can be 1089 // rounded to nearest to a 32-bit unsigned integer 1090 if (x_sign) { // x <= -1 1091 // set invalid flag 1092 *pfpsf |= INVALID_EXCEPTION; 1093 // return Integer Indefinite 1094 res = 0x80000000; 1095 BID_RETURN (res); 1096 } 1097 // 1 <= x <= 2^32 - 1 so x can be rounded 1098 // to nearest to a 32-bit unsigned integer 1099 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10 1100 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x' 1101 // chop off ind digits from the lower part of C1 1102 // C1 fits in 64 bits 1103 // calculate C* and f* 1104 // C* is actually floor(C*) in this case 1105 // C* and f* need shifting and masking, as shown by 1106 // shiftright128[] and maskhigh128[] 1107 // 1 <= x <= 15 1108 // kx = 10^(-x) = ten2mk64[ind - 1] 1109 // C* = C1 * 10^(-x) 1110 // the approximation of 10^(-x) was rounded up to 54 bits 1111 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); 1112 Cstar = P128.w[1]; 1113 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 1114 fstar.w[0] = P128.w[0]; 1115 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g. 1116 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999 1117 // C* = floor(C*) (logical right shift; C has p decimal digits, 1118 // correct by Property 1) 1119 // n = C* * 10^(e+x) 1120 1121 // shift right C* by Ex-64 = shiftright128[ind] 1122 shift = shiftright128[ind - 1]; // 0 <= shift <= 39 1123 Cstar = Cstar >> shift; 1124 // determine inexactness of the rounding of C* 1125 // if (0 < f* < 10^(-x)) then 1126 // the result is exact 1127 // else // if (f* > T*) then 1128 // the result is inexact 1129 if (ind - 1 <= 2) { // fstar.w[1] is 0 1130 if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) { 1131 // ten2mk128trunc[ind -1].w[1] is identical to 1132 // ten2mk128[ind -1].w[1] 1133 Cstar++; 1134 } // else the result is exact 1135 } else { // if 3 <= ind - 1 <= 14 1136 if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) { 1137 // ten2mk128trunc[ind -1].w[1] is identical to 1138 // ten2mk128[ind -1].w[1] 1139 Cstar++; 1140 } // else the result is exact 1141 } 1142 1143 res = Cstar; // the result is positive 1144 } else if (exp == 0) { 1145 // 1 <= q <= 10 1146 // res = +C (exact) 1147 res = C1; // the result is positive 1148 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 1149 // res = +C * 10^exp (exact) 1150 res = C1 * ten2k64[exp]; // the result is positive 1151 } 1152 } 1153 BID_RETURN (res); 1154 } 1155 1156 /***************************************************************************** 1157 * BID64_to_uint32_xceil 1158 ****************************************************************************/ 1159 1160 #if DECIMAL_CALL_BY_REFERENCE 1161 void 1162 bid64_to_uint32_xceil (unsigned int *pres, UINT64 * px 1163 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 1164 _EXC_INFO_PARAM) { 1165 UINT64 x = *px; 1166 #else 1167 unsigned int 1168 bid64_to_uint32_xceil (UINT64 x 1169 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 1170 _EXC_INFO_PARAM) { 1171 #endif 1172 unsigned int res; 1173 UINT64 x_sign; 1174 UINT64 x_exp; 1175 int exp; // unbiased exponent 1176 // Note: C1 represents x_significand (UINT64) 1177 UINT64 tmp64; 1178 BID_UI64DOUBLE tmp1; 1179 unsigned int x_nr_bits; 1180 int q, ind, shift; 1181 UINT64 C1; 1182 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits 1183 UINT128 fstar; 1184 UINT128 P128; 1185 1186 // check for NaN or Infinity 1187 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) { 1188 // set invalid flag 1189 *pfpsf |= INVALID_EXCEPTION; 1190 // return Integer Indefinite 1191 res = 0x80000000; 1192 BID_RETURN (res); 1193 } 1194 // unpack x 1195 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 1196 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] => 1197 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 1198 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased 1199 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 1200 if (C1 > 9999999999999999ull) { // non-canonical 1201 x_exp = 0; 1202 C1 = 0; 1203 } 1204 } else { 1205 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased 1206 C1 = x & MASK_BINARY_SIG1; 1207 } 1208 1209 // check for zeros (possibly from non-canonical values) 1210 if (C1 == 0x0ull) { 1211 // x is 0 1212 res = 0x00000000; 1213 BID_RETURN (res); 1214 } 1215 // x is not special and is not zero 1216 1217 // q = nr. of decimal digits in x (1 <= q <= 54) 1218 // determine first the nr. of bits in x 1219 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 1220 // split the 64-bit value in two 32-bit halves to avoid rounding errors 1221 if (C1 >= 0x0000000100000000ull) { // x >= 2^32 1222 tmp1.d = (double) (C1 >> 32); // exact conversion 1223 x_nr_bits = 1224 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1225 } else { // x < 2^32 1226 tmp1.d = (double) C1; // exact conversion 1227 x_nr_bits = 1228 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1229 } 1230 } else { // if x < 2^53 1231 tmp1.d = (double) C1; // exact conversion 1232 x_nr_bits = 1233 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1234 } 1235 q = nr_digits[x_nr_bits - 1].digits; 1236 if (q == 0) { 1237 q = nr_digits[x_nr_bits - 1].digits1; 1238 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 1239 q++; 1240 } 1241 exp = x_exp - 398; // unbiased exponent 1242 1243 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) 1244 // set invalid flag 1245 *pfpsf |= INVALID_EXCEPTION; 1246 // return Integer Indefinite 1247 res = 0x80000000; 1248 BID_RETURN (res); 1249 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) 1250 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... 1251 // so x rounded to an integer may or may not fit in an unsigned 32-bit int 1252 // the cases that do not fit are identified here; the ones that fit 1253 // fall through and will be handled with other cases further, 1254 // under '1 <= q + exp <= 10' 1255 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1 1256 // => set invalid flag 1257 *pfpsf |= INVALID_EXCEPTION; 1258 // return Integer Indefinite 1259 res = 0x80000000; 1260 BID_RETURN (res); 1261 } else { // if n > 0 and q + exp = 10 1262 // if n > 2^32 - 1 then n is too large 1263 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^32 - 1 1264 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6, 1<=q<=16 1265 // <=> C * 10^(11-q) > 0x9fffffff6, 1<=q<=16 1266 if (q <= 11) { 1267 // Note: C * 10^(11-q) has 10 or 11 digits; 0x9fffffff6 has 11 digits 1268 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int 1269 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) 1270 if (tmp64 > 0x9fffffff6ull) { 1271 // set invalid flag 1272 *pfpsf |= INVALID_EXCEPTION; 1273 // return Integer Indefinite 1274 res = 0x80000000; 1275 BID_RETURN (res); 1276 } 1277 // else cases that can be rounded to a 32-bit unsigned int fall through 1278 // to '1 <= q + exp <= 10' 1279 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2 1280 // C * 10^(11-q) > 0x9fffffff6 <=> 1281 // C > 0x9fffffff6 * 10^(q-11) where 1 <= q - 11 <= 5 1282 // (scale 2^32-1 up) 1283 // Note: 0x9fffffff6*10^(q-11) has q-1 or q digits, where q <= 16 1284 tmp64 = 0x9fffffff6ull * ten2k64[q - 11]; 1285 if (C1 > tmp64) { 1286 // set invalid flag 1287 *pfpsf |= INVALID_EXCEPTION; 1288 // return Integer Indefinite 1289 res = 0x80000000; 1290 BID_RETURN (res); 1291 } 1292 // else cases that can be rounded to a 32-bit int fall through 1293 // to '1 <= q + exp <= 10' 1294 } 1295 } 1296 } 1297 // n is not too large to be converted to int32 if -1 < n < 2^32 1298 // Note: some of the cases tested for above fall through to this point 1299 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1) 1300 // set inexact flag 1301 *pfpsf |= INEXACT_EXCEPTION; 1302 // return 0 or 1 1303 if (x_sign) 1304 res = 0x00000000; 1305 else 1306 res = 0x00000001; 1307 BID_RETURN (res); 1308 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9) 1309 // x <= -1 or 1 <= x < 2^32 so if positive, x can be 1310 // rounded to nearest to a 32-bit unsigned integer 1311 if (x_sign) { // x <= -1 1312 // set invalid flag 1313 *pfpsf |= INVALID_EXCEPTION; 1314 // return Integer Indefinite 1315 res = 0x80000000; 1316 BID_RETURN (res); 1317 } 1318 // 1 <= x < 2^32 so x can be rounded 1319 // to nearest to a 32-bit unsigned integer 1320 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10 1321 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x' 1322 // chop off ind digits from the lower part of C1 1323 // C1 fits in 64 bits 1324 // calculate C* and f* 1325 // C* is actually floor(C*) in this case 1326 // C* and f* need shifting and masking, as shown by 1327 // shiftright128[] and maskhigh128[] 1328 // 1 <= x <= 15 1329 // kx = 10^(-x) = ten2mk64[ind - 1] 1330 // C* = C1 * 10^(-x) 1331 // the approximation of 10^(-x) was rounded up to 54 bits 1332 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); 1333 Cstar = P128.w[1]; 1334 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 1335 fstar.w[0] = P128.w[0]; 1336 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g. 1337 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999 1338 // C* = floor(C*) (logical right shift; C has p decimal digits, 1339 // correct by Property 1) 1340 // n = C* * 10^(e+x) 1341 1342 // shift right C* by Ex-64 = shiftright128[ind] 1343 shift = shiftright128[ind - 1]; // 0 <= shift <= 39 1344 Cstar = Cstar >> shift; 1345 // determine inexactness of the rounding of C* 1346 // if (0 < f* < 10^(-x)) then 1347 // the result is exact 1348 // else // if (f* > T*) then 1349 // the result is inexact 1350 if (ind - 1 <= 2) { // fstar.w[1] is 0 1351 if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) { 1352 // ten2mk128trunc[ind -1].w[1] is identical to 1353 // ten2mk128[ind -1].w[1] 1354 Cstar++; 1355 // set the inexact flag 1356 *pfpsf |= INEXACT_EXCEPTION; 1357 } // else the result is exact 1358 } else { // if 3 <= ind - 1 <= 14 1359 if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) { 1360 // ten2mk128trunc[ind -1].w[1] is identical to 1361 // ten2mk128[ind -1].w[1] 1362 Cstar++; 1363 // set the inexact flag 1364 *pfpsf |= INEXACT_EXCEPTION; 1365 } // else the result is exact 1366 } 1367 1368 res = Cstar; // the result is positive 1369 } else if (exp == 0) { 1370 // 1 <= q <= 10 1371 // res = +C (exact) 1372 res = C1; // the result is positive 1373 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 1374 // res = +C * 10^exp (exact) 1375 res = C1 * ten2k64[exp]; // the result is positive 1376 } 1377 } 1378 BID_RETURN (res); 1379 } 1380 1381 /***************************************************************************** 1382 * BID64_to_uint32_int 1383 ****************************************************************************/ 1384 1385 #if DECIMAL_CALL_BY_REFERENCE 1386 void 1387 bid64_to_uint32_int (unsigned int *pres, UINT64 * px 1388 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) 1389 { 1390 UINT64 x = *px; 1391 #else 1392 unsigned int 1393 bid64_to_uint32_int (UINT64 x 1394 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) 1395 { 1396 #endif 1397 unsigned int res; 1398 UINT64 x_sign; 1399 UINT64 x_exp; 1400 int exp; // unbiased exponent 1401 // Note: C1 represents x_significand (UINT64) 1402 UINT64 tmp64; 1403 BID_UI64DOUBLE tmp1; 1404 unsigned int x_nr_bits; 1405 int q, ind, shift; 1406 UINT64 C1; 1407 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits 1408 UINT128 P128; 1409 1410 // check for NaN or Infinity 1411 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) { 1412 // set invalid flag 1413 *pfpsf |= INVALID_EXCEPTION; 1414 // return Integer Indefinite 1415 res = 0x80000000; 1416 BID_RETURN (res); 1417 } 1418 // unpack x 1419 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 1420 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] => 1421 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 1422 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased 1423 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 1424 if (C1 > 9999999999999999ull) { // non-canonical 1425 x_exp = 0; 1426 C1 = 0; 1427 } 1428 } else { 1429 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased 1430 C1 = x & MASK_BINARY_SIG1; 1431 } 1432 1433 // check for zeros (possibly from non-canonical values) 1434 if (C1 == 0x0ull) { 1435 // x is 0 1436 res = 0x00000000; 1437 BID_RETURN (res); 1438 } 1439 // x is not special and is not zero 1440 1441 // q = nr. of decimal digits in x (1 <= q <= 54) 1442 // determine first the nr. of bits in x 1443 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 1444 // split the 64-bit value in two 32-bit halves to avoid rounding errors 1445 if (C1 >= 0x0000000100000000ull) { // x >= 2^32 1446 tmp1.d = (double) (C1 >> 32); // exact conversion 1447 x_nr_bits = 1448 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1449 } else { // x < 2^32 1450 tmp1.d = (double) C1; // exact conversion 1451 x_nr_bits = 1452 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1453 } 1454 } else { // if x < 2^53 1455 tmp1.d = (double) C1; // exact conversion 1456 x_nr_bits = 1457 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1458 } 1459 q = nr_digits[x_nr_bits - 1].digits; 1460 if (q == 0) { 1461 q = nr_digits[x_nr_bits - 1].digits1; 1462 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 1463 q++; 1464 } 1465 exp = x_exp - 398; // unbiased exponent 1466 1467 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) 1468 // set invalid flag 1469 *pfpsf |= INVALID_EXCEPTION; 1470 // return Integer Indefinite 1471 res = 0x80000000; 1472 BID_RETURN (res); 1473 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) 1474 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... 1475 // so x rounded to an integer may or may not fit in an unsigned 32-bit int 1476 // the cases that do not fit are identified here; the ones that fit 1477 // fall through and will be handled with other cases further, 1478 // under '1 <= q + exp <= 10' 1479 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1 1480 // => set invalid flag 1481 *pfpsf |= INVALID_EXCEPTION; 1482 // return Integer Indefinite 1483 res = 0x80000000; 1484 BID_RETURN (res); 1485 } else { // if n > 0 and q + exp = 10 1486 // if n >= 2^32 then n is too large 1487 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32 1488 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=16 1489 // <=> C * 10^(11-q) >= 0xa00000000, 1<=q<=16 1490 if (q <= 11) { 1491 // Note: C * 10^(11-q) has 10 or 11 digits; 0xa00000000 has 11 digits 1492 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int 1493 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) 1494 if (tmp64 >= 0xa00000000ull) { 1495 // set invalid flag 1496 *pfpsf |= INVALID_EXCEPTION; 1497 // return Integer Indefinite 1498 res = 0x80000000; 1499 BID_RETURN (res); 1500 } 1501 // else cases that can be rounded to a 32-bit unsigned int fall through 1502 // to '1 <= q + exp <= 10' 1503 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2 1504 // C * 10^(11-q) >= 0xa00000000 <=> 1505 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 5 1506 // (scale 2^32-1/2 up) 1507 // Note: 0xa00000000*10^(q-11) has q-1 or q digits, where q <= 16 1508 tmp64 = 0xa00000000ull * ten2k64[q - 11]; 1509 if (C1 >= tmp64) { 1510 // set invalid flag 1511 *pfpsf |= INVALID_EXCEPTION; 1512 // return Integer Indefinite 1513 res = 0x80000000; 1514 BID_RETURN (res); 1515 } 1516 // else cases that can be rounded to a 32-bit int fall through 1517 // to '1 <= q + exp <= 10' 1518 } 1519 } 1520 } 1521 // n is not too large to be converted to int32 if -1 < n < 2^32 1522 // Note: some of the cases tested for above fall through to this point 1523 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1) 1524 // return 0 1525 res = 0x00000000; 1526 BID_RETURN (res); 1527 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9) 1528 // x <= -1 or 1 <= x < 2^32 so if positive, x can be 1529 // rounded to nearest to a 32-bit unsigned integer 1530 if (x_sign) { // x <= -1 1531 // set invalid flag 1532 *pfpsf |= INVALID_EXCEPTION; 1533 // return Integer Indefinite 1534 res = 0x80000000; 1535 BID_RETURN (res); 1536 } 1537 // 1 <= x < 2^32 so x can be rounded 1538 // to nearest to a 32-bit unsigned integer 1539 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10 1540 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x' 1541 // chop off ind digits from the lower part of C1 1542 // C1 fits in 64 bits 1543 // calculate C* and f* 1544 // C* is actually floor(C*) in this case 1545 // C* and f* need shifting and masking, as shown by 1546 // shiftright128[] and maskhigh128[] 1547 // 1 <= x <= 15 1548 // kx = 10^(-x) = ten2mk64[ind - 1] 1549 // C* = C1 * 10^(-x) 1550 // the approximation of 10^(-x) was rounded up to 54 bits 1551 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); 1552 Cstar = P128.w[1]; 1553 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g. 1554 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999 1555 // C* = floor(C*) (logical right shift; C has p decimal digits, 1556 // correct by Property 1) 1557 // n = C* * 10^(e+x) 1558 1559 // shift right C* by Ex-64 = shiftright128[ind] 1560 shift = shiftright128[ind - 1]; // 0 <= shift <= 39 1561 Cstar = Cstar >> shift; 1562 1563 res = Cstar; // the result is positive 1564 } else if (exp == 0) { 1565 // 1 <= q <= 10 1566 // res = +C (exact) 1567 res = C1; // the result is positive 1568 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 1569 // res = +C * 10^exp (exact) 1570 res = C1 * ten2k64[exp]; // the result is positive 1571 } 1572 } 1573 BID_RETURN (res); 1574 } 1575 1576 /***************************************************************************** 1577 * BID64_to_uint32_xint 1578 ****************************************************************************/ 1579 1580 #if DECIMAL_CALL_BY_REFERENCE 1581 void 1582 bid64_to_uint32_xint (unsigned int *pres, UINT64 * px 1583 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 1584 _EXC_INFO_PARAM) { 1585 UINT64 x = *px; 1586 #else 1587 unsigned int 1588 bid64_to_uint32_xint (UINT64 x 1589 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 1590 _EXC_INFO_PARAM) { 1591 #endif 1592 unsigned int res; 1593 UINT64 x_sign; 1594 UINT64 x_exp; 1595 int exp; // unbiased exponent 1596 // Note: C1 represents x_significand (UINT64) 1597 UINT64 tmp64; 1598 BID_UI64DOUBLE tmp1; 1599 unsigned int x_nr_bits; 1600 int q, ind, shift; 1601 UINT64 C1; 1602 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits 1603 UINT128 fstar; 1604 UINT128 P128; 1605 1606 // check for NaN or Infinity 1607 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) { 1608 // set invalid flag 1609 *pfpsf |= INVALID_EXCEPTION; 1610 // return Integer Indefinite 1611 res = 0x80000000; 1612 BID_RETURN (res); 1613 } 1614 // unpack x 1615 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 1616 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] => 1617 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 1618 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased 1619 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 1620 if (C1 > 9999999999999999ull) { // non-canonical 1621 x_exp = 0; 1622 C1 = 0; 1623 } 1624 } else { 1625 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased 1626 C1 = x & MASK_BINARY_SIG1; 1627 } 1628 1629 // check for zeros (possibly from non-canonical values) 1630 if (C1 == 0x0ull) { 1631 // x is 0 1632 res = 0x00000000; 1633 BID_RETURN (res); 1634 } 1635 // x is not special and is not zero 1636 1637 // q = nr. of decimal digits in x (1 <= q <= 54) 1638 // determine first the nr. of bits in x 1639 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 1640 // split the 64-bit value in two 32-bit halves to avoid rounding errors 1641 if (C1 >= 0x0000000100000000ull) { // x >= 2^32 1642 tmp1.d = (double) (C1 >> 32); // exact conversion 1643 x_nr_bits = 1644 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1645 } else { // x < 2^32 1646 tmp1.d = (double) C1; // exact conversion 1647 x_nr_bits = 1648 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1649 } 1650 } else { // if x < 2^53 1651 tmp1.d = (double) C1; // exact conversion 1652 x_nr_bits = 1653 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1654 } 1655 q = nr_digits[x_nr_bits - 1].digits; 1656 if (q == 0) { 1657 q = nr_digits[x_nr_bits - 1].digits1; 1658 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 1659 q++; 1660 } 1661 exp = x_exp - 398; // unbiased exponent 1662 1663 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) 1664 // set invalid flag 1665 *pfpsf |= INVALID_EXCEPTION; 1666 // return Integer Indefinite 1667 res = 0x80000000; 1668 BID_RETURN (res); 1669 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) 1670 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... 1671 // so x rounded to an integer may or may not fit in an unsigned 32-bit int 1672 // the cases that do not fit are identified here; the ones that fit 1673 // fall through and will be handled with other cases further, 1674 // under '1 <= q + exp <= 10' 1675 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1 1676 // => set invalid flag 1677 *pfpsf |= INVALID_EXCEPTION; 1678 // return Integer Indefinite 1679 res = 0x80000000; 1680 BID_RETURN (res); 1681 } else { // if n > 0 and q + exp = 10 1682 // if n >= 2^32 then n is too large 1683 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32 1684 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=16 1685 // <=> C * 10^(11-q) >= 0xa00000000, 1<=q<=16 1686 if (q <= 11) { 1687 // Note: C * 10^(11-q) has 10 or 11 digits; 0xa00000000 has 11 digits 1688 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int 1689 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) 1690 if (tmp64 >= 0xa00000000ull) { 1691 // set invalid flag 1692 *pfpsf |= INVALID_EXCEPTION; 1693 // return Integer Indefinite 1694 res = 0x80000000; 1695 BID_RETURN (res); 1696 } 1697 // else cases that can be rounded to a 32-bit unsigned int fall through 1698 // to '1 <= q + exp <= 10' 1699 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2 1700 // C * 10^(11-q) >= 0xa00000000 <=> 1701 // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 5 1702 // (scale 2^32-1/2 up) 1703 // Note: 0xa00000000*10^(q-11) has q-1 or q digits, where q <= 16 1704 tmp64 = 0xa00000000ull * ten2k64[q - 11]; 1705 if (C1 >= tmp64) { 1706 // set invalid flag 1707 *pfpsf |= INVALID_EXCEPTION; 1708 // return Integer Indefinite 1709 res = 0x80000000; 1710 BID_RETURN (res); 1711 } 1712 // else cases that can be rounded to a 32-bit int fall through 1713 // to '1 <= q + exp <= 10' 1714 } 1715 } 1716 } 1717 // n is not too large to be converted to int32 if -1 < n < 2^32 1718 // Note: some of the cases tested for above fall through to this point 1719 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1) 1720 // set inexact flag 1721 *pfpsf |= INEXACT_EXCEPTION; 1722 // return 0 1723 res = 0x00000000; 1724 BID_RETURN (res); 1725 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9) 1726 // x <= -1 or 1 <= x < 2^32 so if positive, x can be 1727 // rounded to nearest to a 32-bit unsigned integer 1728 if (x_sign) { // x <= -1 1729 // set invalid flag 1730 *pfpsf |= INVALID_EXCEPTION; 1731 // return Integer Indefinite 1732 res = 0x80000000; 1733 BID_RETURN (res); 1734 } 1735 // 1 <= x < 2^32 so x can be rounded 1736 // to nearest to a 32-bit unsigned integer 1737 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10 1738 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x' 1739 // chop off ind digits from the lower part of C1 1740 // C1 fits in 64 bits 1741 // calculate C* and f* 1742 // C* is actually floor(C*) in this case 1743 // C* and f* need shifting and masking, as shown by 1744 // shiftright128[] and maskhigh128[] 1745 // 1 <= x <= 15 1746 // kx = 10^(-x) = ten2mk64[ind - 1] 1747 // C* = C1 * 10^(-x) 1748 // the approximation of 10^(-x) was rounded up to 54 bits 1749 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); 1750 Cstar = P128.w[1]; 1751 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 1752 fstar.w[0] = P128.w[0]; 1753 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g. 1754 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999 1755 // C* = floor(C*) (logical right shift; C has p decimal digits, 1756 // correct by Property 1) 1757 // n = C* * 10^(e+x) 1758 1759 // shift right C* by Ex-64 = shiftright128[ind] 1760 shift = shiftright128[ind - 1]; // 0 <= shift <= 39 1761 Cstar = Cstar >> shift; 1762 // determine inexactness of the rounding of C* 1763 // if (0 < f* < 10^(-x)) then 1764 // the result is exact 1765 // else // if (f* > T*) then 1766 // the result is inexact 1767 if (ind - 1 <= 2) { // fstar.w[1] is 0 1768 if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) { 1769 // ten2mk128trunc[ind -1].w[1] is identical to 1770 // ten2mk128[ind -1].w[1] 1771 // set the inexact flag 1772 *pfpsf |= INEXACT_EXCEPTION; 1773 } // else the result is exact 1774 } else { // if 3 <= ind - 1 <= 14 1775 if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) { 1776 // ten2mk128trunc[ind -1].w[1] is identical to 1777 // ten2mk128[ind -1].w[1] 1778 // set the inexact flag 1779 *pfpsf |= INEXACT_EXCEPTION; 1780 } // else the result is exact 1781 } 1782 1783 res = Cstar; // the result is positive 1784 } else if (exp == 0) { 1785 // 1 <= q <= 10 1786 // res = +C (exact) 1787 res = C1; // the result is positive 1788 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 1789 // res = +C * 10^exp (exact) 1790 res = C1 * ten2k64[exp]; // the result is positive 1791 } 1792 } 1793 BID_RETURN (res); 1794 } 1795 1796 /***************************************************************************** 1797 * BID64_to_uint32_rninta 1798 ****************************************************************************/ 1799 1800 #if DECIMAL_CALL_BY_REFERENCE 1801 void 1802 bid64_to_uint32_rninta (unsigned int *pres, UINT64 * px 1803 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 1804 _EXC_INFO_PARAM) { 1805 UINT64 x = *px; 1806 #else 1807 unsigned int 1808 bid64_to_uint32_rninta (UINT64 x 1809 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 1810 _EXC_INFO_PARAM) { 1811 #endif 1812 unsigned int res; 1813 UINT64 x_sign; 1814 UINT64 x_exp; 1815 int exp; // unbiased exponent 1816 // Note: C1 represents x_significand (UINT64) 1817 UINT64 tmp64; 1818 BID_UI64DOUBLE tmp1; 1819 unsigned int x_nr_bits; 1820 int q, ind, shift; 1821 UINT64 C1; 1822 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits 1823 UINT128 P128; 1824 1825 // check for NaN or Infinity 1826 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) { 1827 // set invalid flag 1828 *pfpsf |= INVALID_EXCEPTION; 1829 // return Integer Indefinite 1830 res = 0x80000000; 1831 BID_RETURN (res); 1832 } 1833 // unpack x 1834 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 1835 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] => 1836 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 1837 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased 1838 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 1839 if (C1 > 9999999999999999ull) { // non-canonical 1840 x_exp = 0; 1841 C1 = 0; 1842 } 1843 } else { 1844 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased 1845 C1 = x & MASK_BINARY_SIG1; 1846 } 1847 1848 // check for zeros (possibly from non-canonical values) 1849 if (C1 == 0x0ull) { 1850 // x is 0 1851 res = 0x00000000; 1852 BID_RETURN (res); 1853 } 1854 // x is not special and is not zero 1855 1856 // q = nr. of decimal digits in x (1 <= q <= 54) 1857 // determine first the nr. of bits in x 1858 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 1859 // split the 64-bit value in two 32-bit halves to avoid rounding errors 1860 if (C1 >= 0x0000000100000000ull) { // x >= 2^32 1861 tmp1.d = (double) (C1 >> 32); // exact conversion 1862 x_nr_bits = 1863 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1864 } else { // x < 2^32 1865 tmp1.d = (double) C1; // exact conversion 1866 x_nr_bits = 1867 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1868 } 1869 } else { // if x < 2^53 1870 tmp1.d = (double) C1; // exact conversion 1871 x_nr_bits = 1872 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1873 } 1874 q = nr_digits[x_nr_bits - 1].digits; 1875 if (q == 0) { 1876 q = nr_digits[x_nr_bits - 1].digits1; 1877 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 1878 q++; 1879 } 1880 exp = x_exp - 398; // unbiased exponent 1881 1882 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) 1883 // set invalid flag 1884 *pfpsf |= INVALID_EXCEPTION; 1885 // return Integer Indefinite 1886 res = 0x80000000; 1887 BID_RETURN (res); 1888 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) 1889 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... 1890 // so x rounded to an integer may or may not fit in an unsigned 32-bit int 1891 // the cases that do not fit are identified here; the ones that fit 1892 // fall through and will be handled with other cases further, 1893 // under '1 <= q + exp <= 10' 1894 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1/2 1895 // => set invalid flag 1896 *pfpsf |= INVALID_EXCEPTION; 1897 // return Integer Indefinite 1898 res = 0x80000000; 1899 BID_RETURN (res); 1900 } else { // if n > 0 and q + exp = 10 1901 // if n >= 2^32 - 1/2 then n is too large 1902 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2 1903 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=16 1904 // <=> C * 10^(11-q) >= 0x9fffffffb, 1<=q<=16 1905 if (q <= 11) { 1906 // Note: C * 10^(11-q) has 10 or 11 digits; 0x9fffffffb has 11 digits 1907 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int 1908 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) 1909 if (tmp64 >= 0x9fffffffbull) { 1910 // set invalid flag 1911 *pfpsf |= INVALID_EXCEPTION; 1912 // return Integer Indefinite 1913 res = 0x80000000; 1914 BID_RETURN (res); 1915 } 1916 // else cases that can be rounded to a 32-bit unsigned int fall through 1917 // to '1 <= q + exp <= 10' 1918 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2 1919 // C * 10^(11-q) >= 0x9fffffffb <=> 1920 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 5 1921 // (scale 2^32-1/2 up) 1922 // Note: 0x9fffffffb*10^(q-11) has q-1 or q digits, where q <= 16 1923 tmp64 = 0x9fffffffbull * ten2k64[q - 11]; 1924 if (C1 >= tmp64) { 1925 // set invalid flag 1926 *pfpsf |= INVALID_EXCEPTION; 1927 // return Integer Indefinite 1928 res = 0x80000000; 1929 BID_RETURN (res); 1930 } 1931 // else cases that can be rounded to a 32-bit int fall through 1932 // to '1 <= q + exp <= 10' 1933 } 1934 } 1935 } 1936 // n is not too large to be converted to int32 if -1/2 < n < 2^32 - 1/2 1937 // Note: some of the cases tested for above fall through to this point 1938 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1) 1939 // return 0 1940 res = 0x00000000; 1941 BID_RETURN (res); 1942 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1) 1943 // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1) 1944 // res = 0 1945 // else if x > 0 1946 // res = +1 1947 // else // if x < 0 1948 // invalid exc 1949 ind = q - 1; 1950 if (C1 < midpoint64[ind]) { 1951 res = 0x00000000; // return 0 1952 } else if (x_sign) { // n < 0 1953 // set invalid flag 1954 *pfpsf |= INVALID_EXCEPTION; 1955 // return Integer Indefinite 1956 res = 0x80000000; 1957 BID_RETURN (res); 1958 } else { // n > 0 1959 res = 0x00000001; // return +1 1960 } 1961 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9) 1962 // -2^32-1/2 <= x <= -1 or 1 <= x < 2^32-1/2 so if positive, x can be 1963 // rounded to nearest to a 32-bit unsigned integer 1964 if (x_sign) { // x <= -1 1965 // set invalid flag 1966 *pfpsf |= INVALID_EXCEPTION; 1967 // return Integer Indefinite 1968 res = 0x80000000; 1969 BID_RETURN (res); 1970 } 1971 // 1 <= x < 2^32-1/2 so x can be rounded 1972 // to nearest to a 32-bit unsigned integer 1973 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10 1974 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x' 1975 // chop off ind digits from the lower part of C1 1976 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits 1977 C1 = C1 + midpoint64[ind - 1]; 1978 // calculate C* and f* 1979 // C* is actually floor(C*) in this case 1980 // C* and f* need shifting and masking, as shown by 1981 // shiftright128[] and maskhigh128[] 1982 // 1 <= x <= 15 1983 // kx = 10^(-x) = ten2mk64[ind - 1] 1984 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 1985 // the approximation of 10^(-x) was rounded up to 54 bits 1986 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); 1987 Cstar = P128.w[1]; 1988 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g. 1989 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999 1990 // C* = floor(C*) (logical right shift; C has p decimal digits, 1991 // correct by Property 1) 1992 // n = C* * 10^(e+x) 1993 1994 // shift right C* by Ex-64 = shiftright128[ind] 1995 shift = shiftright128[ind - 1]; // 0 <= shift <= 39 1996 Cstar = Cstar >> shift; 1997 1998 // if the result was a midpoint it was rounded away from zero 1999 res = Cstar; // the result is positive 2000 } else if (exp == 0) { 2001 // 1 <= q <= 10 2002 // res = +C (exact) 2003 res = C1; // the result is positive 2004 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 2005 // res = +C * 10^exp (exact) 2006 res = C1 * ten2k64[exp]; // the result is positive 2007 } 2008 } 2009 BID_RETURN (res); 2010 } 2011 2012 /***************************************************************************** 2013 * BID64_to_uint32_xrninta 2014 ****************************************************************************/ 2015 2016 #if DECIMAL_CALL_BY_REFERENCE 2017 void 2018 bid64_to_uint32_xrninta (unsigned int *pres, UINT64 * px 2019 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 2020 _EXC_INFO_PARAM) { 2021 UINT64 x = *px; 2022 #else 2023 unsigned int 2024 bid64_to_uint32_xrninta (UINT64 x 2025 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 2026 _EXC_INFO_PARAM) { 2027 #endif 2028 unsigned int res; 2029 UINT64 x_sign; 2030 UINT64 x_exp; 2031 int exp; // unbiased exponent 2032 // Note: C1 represents x_significand (UINT64) 2033 UINT64 tmp64; 2034 BID_UI64DOUBLE tmp1; 2035 unsigned int x_nr_bits; 2036 int q, ind, shift; 2037 UINT64 C1; 2038 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits 2039 UINT128 fstar; 2040 UINT128 P128; 2041 2042 // check for NaN or Infinity 2043 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) { 2044 // set invalid flag 2045 *pfpsf |= INVALID_EXCEPTION; 2046 // return Integer Indefinite 2047 res = 0x80000000; 2048 BID_RETURN (res); 2049 } 2050 // unpack x 2051 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 2052 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] => 2053 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 2054 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased 2055 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 2056 if (C1 > 9999999999999999ull) { // non-canonical 2057 x_exp = 0; 2058 C1 = 0; 2059 } 2060 } else { 2061 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased 2062 C1 = x & MASK_BINARY_SIG1; 2063 } 2064 2065 // check for zeros (possibly from non-canonical values) 2066 if (C1 == 0x0ull) { 2067 // x is 0 2068 res = 0x00000000; 2069 BID_RETURN (res); 2070 } 2071 // x is not special and is not zero 2072 2073 // q = nr. of decimal digits in x (1 <= q <= 54) 2074 // determine first the nr. of bits in x 2075 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 2076 // split the 64-bit value in two 32-bit halves to avoid rounding errors 2077 if (C1 >= 0x0000000100000000ull) { // x >= 2^32 2078 tmp1.d = (double) (C1 >> 32); // exact conversion 2079 x_nr_bits = 2080 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 2081 } else { // x < 2^32 2082 tmp1.d = (double) C1; // exact conversion 2083 x_nr_bits = 2084 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 2085 } 2086 } else { // if x < 2^53 2087 tmp1.d = (double) C1; // exact conversion 2088 x_nr_bits = 2089 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 2090 } 2091 q = nr_digits[x_nr_bits - 1].digits; 2092 if (q == 0) { 2093 q = nr_digits[x_nr_bits - 1].digits1; 2094 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 2095 q++; 2096 } 2097 exp = x_exp - 398; // unbiased exponent 2098 2099 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) 2100 // set invalid flag 2101 *pfpsf |= INVALID_EXCEPTION; 2102 // return Integer Indefinite 2103 res = 0x80000000; 2104 BID_RETURN (res); 2105 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) 2106 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... 2107 // so x rounded to an integer may or may not fit in an unsigned 32-bit int 2108 // the cases that do not fit are identified here; the ones that fit 2109 // fall through and will be handled with other cases further, 2110 // under '1 <= q + exp <= 10' 2111 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1/2 2112 // => set invalid flag 2113 *pfpsf |= INVALID_EXCEPTION; 2114 // return Integer Indefinite 2115 res = 0x80000000; 2116 BID_RETURN (res); 2117 } else { // if n > 0 and q + exp = 10 2118 // if n >= 2^32 - 1/2 then n is too large 2119 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2 2120 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=16 2121 // <=> C * 10^(11-q) >= 0x9fffffffb, 1<=q<=16 2122 if (q <= 11) { 2123 // Note: C * 10^(11-q) has 10 or 11 digits; 0x9fffffffb has 11 digits 2124 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int 2125 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) 2126 if (tmp64 >= 0x9fffffffbull) { 2127 // set invalid flag 2128 *pfpsf |= INVALID_EXCEPTION; 2129 // return Integer Indefinite 2130 res = 0x80000000; 2131 BID_RETURN (res); 2132 } 2133 // else cases that can be rounded to a 32-bit unsigned int fall through 2134 // to '1 <= q + exp <= 10' 2135 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2 2136 // C * 10^(11-q) >= 0x9fffffffb <=> 2137 // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 5 2138 // (scale 2^32-1/2 up) 2139 // Note: 0x9fffffffb*10^(q-11) has q-1 or q digits, where q <= 16 2140 tmp64 = 0x9fffffffbull * ten2k64[q - 11]; 2141 if (C1 >= tmp64) { 2142 // set invalid flag 2143 *pfpsf |= INVALID_EXCEPTION; 2144 // return Integer Indefinite 2145 res = 0x80000000; 2146 BID_RETURN (res); 2147 } 2148 // else cases that can be rounded to a 32-bit int fall through 2149 // to '1 <= q + exp <= 10' 2150 } 2151 } 2152 } 2153 // n is not too large to be converted to int32 if -1/2 < n < 2^32 - 1/2 2154 // Note: some of the cases tested for above fall through to this point 2155 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1) 2156 // set inexact flag 2157 *pfpsf |= INEXACT_EXCEPTION; 2158 // return 0 2159 res = 0x00000000; 2160 BID_RETURN (res); 2161 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1) 2162 // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1) 2163 // res = 0 2164 // else if x > 0 2165 // res = +1 2166 // else // if x < 0 2167 // invalid exc 2168 ind = q - 1; 2169 if (C1 < midpoint64[ind]) { 2170 res = 0x00000000; // return 0 2171 } else if (x_sign) { // n < 0 2172 // set invalid flag 2173 *pfpsf |= INVALID_EXCEPTION; 2174 // return Integer Indefinite 2175 res = 0x80000000; 2176 BID_RETURN (res); 2177 } else { // n > 0 2178 res = 0x00000001; // return +1 2179 } 2180 // set inexact flag 2181 *pfpsf |= INEXACT_EXCEPTION; 2182 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9) 2183 // -2^32-1/2 <= x <= -1 or 1 <= x < 2^32-1/2 so if positive, x can be 2184 // rounded to nearest to a 32-bit unsigned integer 2185 if (x_sign) { // x <= -1 2186 // set invalid flag 2187 *pfpsf |= INVALID_EXCEPTION; 2188 // return Integer Indefinite 2189 res = 0x80000000; 2190 BID_RETURN (res); 2191 } 2192 // 1 <= x < 2^32-1/2 so x can be rounded 2193 // to nearest to a 32-bit unsigned integer 2194 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10 2195 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x' 2196 // chop off ind digits from the lower part of C1 2197 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits 2198 C1 = C1 + midpoint64[ind - 1]; 2199 // calculate C* and f* 2200 // C* is actually floor(C*) in this case 2201 // C* and f* need shifting and masking, as shown by 2202 // shiftright128[] and maskhigh128[] 2203 // 1 <= x <= 15 2204 // kx = 10^(-x) = ten2mk64[ind - 1] 2205 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 2206 // the approximation of 10^(-x) was rounded up to 54 bits 2207 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); 2208 Cstar = P128.w[1]; 2209 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 2210 fstar.w[0] = P128.w[0]; 2211 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g. 2212 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999 2213 // C* = floor(C*) (logical right shift; C has p decimal digits, 2214 // correct by Property 1) 2215 // n = C* * 10^(e+x) 2216 2217 // shift right C* by Ex-64 = shiftright128[ind] 2218 shift = shiftright128[ind - 1]; // 0 <= shift <= 39 2219 Cstar = Cstar >> shift; 2220 2221 // determine inexactness of the rounding of C* 2222 // if (0 < f* - 1/2 < 10^(-x)) then 2223 // the result is exact 2224 // else // if (f* - 1/2 > T*) then 2225 // the result is inexact 2226 if (ind - 1 <= 2) { // fstar.w[1] is 0 2227 if (fstar.w[0] > 0x8000000000000000ull) { 2228 // f* > 1/2 and the result may be exact 2229 tmp64 = fstar.w[0] - 0x8000000000000000ull; // f* - 1/2 2230 if ((tmp64 > ten2mk128trunc[ind - 1].w[1])) { 2231 // ten2mk128trunc[ind -1].w[1] is identical to 2232 // ten2mk128[ind -1].w[1] 2233 // set the inexact flag 2234 *pfpsf |= INEXACT_EXCEPTION; 2235 } // else the result is exact 2236 } else { // the result is inexact; f2* <= 1/2 2237 // set the inexact flag 2238 *pfpsf |= INEXACT_EXCEPTION; 2239 } 2240 } else { // if 3 <= ind - 1 <= 14 2241 if (fstar.w[1] > onehalf128[ind - 1] || 2242 (fstar.w[1] == onehalf128[ind - 1] && fstar.w[0])) { 2243 // f2* > 1/2 and the result may be exact 2244 // Calculate f2* - 1/2 2245 tmp64 = fstar.w[1] - onehalf128[ind - 1]; 2246 if (tmp64 || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) { 2247 // ten2mk128trunc[ind -1].w[1] is identical to 2248 // ten2mk128[ind -1].w[1] 2249 // set the inexact flag 2250 *pfpsf |= INEXACT_EXCEPTION; 2251 } // else the result is exact 2252 } else { // the result is inexact; f2* <= 1/2 2253 // set the inexact flag 2254 *pfpsf |= INEXACT_EXCEPTION; 2255 } 2256 } 2257 2258 // if the result was a midpoint it was rounded away from zero 2259 res = Cstar; // the result is positive 2260 } else if (exp == 0) { 2261 // 1 <= q <= 10 2262 // res = +C (exact) 2263 res = C1; // the result is positive 2264 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 2265 // res = +C * 10^exp (exact) 2266 res = C1 * ten2k64[exp]; // the result is positive 2267 } 2268 } 2269 BID_RETURN (res); 2270 } 2271
Indexes created Sat Feb 28 05:31:39 UTC 2026