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_uint64_rnint 28 ****************************************************************************/ 29 30#if DECIMAL_CALL_BY_REFERENCE 31void 32bid64_to_uint64_rnint (UINT64 * pres, UINT64 * px 33 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 34 _EXC_INFO_PARAM) { 35 UINT64 x = *px; 36#else 37UINT64 38bid64_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 266void 267bid64_to_uint64_xrnint (UINT64 * pres, UINT64 * px 268 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 269 _EXC_INFO_PARAM) { 270 UINT64 x = *px; 271#else 272UINT64 273bid64_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 542void 543bid64_to_uint64_floor (UINT64 * pres, UINT64 * px 544 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 545 _EXC_INFO_PARAM) { 546 UINT64 x = *px; 547#else 548UINT64 549bid64_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 727void 728bid64_to_uint64_xfloor (UINT64 * pres, UINT64 * px 729 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 730 _EXC_INFO_PARAM) { 731 UINT64 x = *px; 732#else 733UINT64 734bid64_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 938void 939bid64_to_uint64_ceil (UINT64 * pres, UINT64 * px 940 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 941 _EXC_INFO_PARAM) { 942 UINT64 x = *px; 943#else 944UINT64 945bid64_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 1157void 1158bid64_to_uint64_xceil (UINT64 * pres, UINT64 * px 1159 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 1160 _EXC_INFO_PARAM) { 1161 UINT64 x = *px; 1162#else 1163UINT64 1164bid64_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 1382void 1383bid64_to_uint64_int (UINT64 * pres, UINT64 * px 1384 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) 1385{ 1386 UINT64 x = *px; 1387#else 1388UINT64 1389bid64_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 1576void 1577bid64_to_uint64_xint (UINT64 * pres, UINT64 * px 1578 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 1579 _EXC_INFO_PARAM) { 1580 UINT64 x = *px; 1581#else 1582UINT64 1583bid64_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 1796void 1797bid64_to_uint64_rninta (UINT64 * pres, UINT64 * px 1798 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 1799 _EXC_INFO_PARAM) { 1800 UINT64 x = *px; 1801#else 1802UINT64 1803bid64_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 2016void 2017bid64_to_uint64_xrninta (UINT64 * pres, UINT64 * px 2018 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 2019 _EXC_INFO_PARAM) { 2020 UINT64 x = *px; 2021#else 2022UINT64 2023bid64_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