1/* Copyright (C) 2007-2022 Free Software Foundation, Inc. 2 3This file is part of GCC. 4 5GCC is free software; you can redistribute it and/or modify it under 6the terms of the GNU General Public License as published by the Free 7Software Foundation; either version 3, or (at your option) any later 8version. 9 10GCC is distributed in the hope that it will be useful, but WITHOUT ANY 11WARRANTY; without even the implied warranty of MERCHANTABILITY or 12FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 13for more details. 14 15Under Section 7 of GPL version 3, you are granted additional 16permissions described in the GCC Runtime Library Exception, version 173.1, as published by the Free Software Foundation. 18 19You should have received a copy of the GNU General Public License and 20a copy of the GCC Runtime Library Exception along with this program; 21see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 22<http://www.gnu.org/licenses/>. */ 23 24#include "bid_internal.h" 25 26/***************************************************************************** 27 * BID64_to_uint32_rnint 28 ****************************************************************************/ 29 30#if DECIMAL_CALL_BY_REFERENCE 31void 32bid64_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 37unsigned int 38bid64_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 268void 269bid64_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 274unsigned int 275bid64_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 545void 546bid64_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 551unsigned int 552bid64_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 731void 732bid64_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 737unsigned int 738bid64_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 942void 943bid64_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 948unsigned int 949bid64_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 1161void 1162bid64_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 1167unsigned int 1168bid64_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 1386void 1387bid64_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 1392unsigned int 1393bid64_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 1581void 1582bid64_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 1587unsigned int 1588bid64_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 1801void 1802bid64_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 1807unsigned int 1808bid64_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 2017void 2018bid64_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 2023unsigned int 2024bid64_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