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