1 /* Copyright (C) 2007, 2009 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_round_integral_exact 28 ****************************************************************************/ 29 30 #if DECIMAL_CALL_BY_REFERENCE 31 void 32 bid64_round_integral_exact (UINT64 * pres, 33 UINT64 * 34 px _RND_MODE_PARAM _EXC_FLAGS_PARAM 35 _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 36 UINT64 x = *px; 37 #if !DECIMAL_GLOBAL_ROUNDING 38 unsigned int rnd_mode = *prnd_mode; 39 #endif 40 #else 41 UINT64 42 bid64_round_integral_exact (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM 43 _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 44 #endif 45 46 UINT64 res = 0xbaddbaddbaddbaddull; 47 UINT64 x_sign; 48 int exp; // unbiased exponent 49 // Note: C1 represents the significand (UINT64) 50 BID_UI64DOUBLE tmp1; 51 int x_nr_bits; 52 int q, ind, shift; 53 UINT64 C1; 54 // UINT64 res is C* at first - represents up to 16 decimal digits <= 54 bits 55 UINT128 fstar = { {0x0ull, 0x0ull} }; 56 UINT128 P128; 57 58 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 59 60 // check for NaNs and infinities 61 if ((x & MASK_NAN) == MASK_NAN) { // check for NaN 62 if ((x & 0x0003ffffffffffffull) > 999999999999999ull) 63 x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits 64 else 65 x = x & 0xfe03ffffffffffffull; // clear G6-G12 66 if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN 67 // set invalid flag 68 *pfpsf |= INVALID_EXCEPTION; 69 // return quiet (SNaN) 70 res = x & 0xfdffffffffffffffull; 71 } else { // QNaN 72 res = x; 73 } 74 BID_RETURN (res); 75 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity 76 res = x_sign | 0x7800000000000000ull; 77 BID_RETURN (res); 78 } 79 // unpack x 80 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 81 // if the steering bits are 11 (condition will be 0), then 82 // the exponent is G[0:w+1] 83 exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398; 84 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 85 if (C1 > 9999999999999999ull) { // non-canonical 86 C1 = 0; 87 } 88 } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) 89 exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398; 90 C1 = (x & MASK_BINARY_SIG1); 91 } 92 93 // if x is 0 or non-canonical return 0 preserving the sign bit and 94 // the preferred exponent of MAX(Q(x), 0) 95 if (C1 == 0) { 96 if (exp < 0) 97 exp = 0; 98 res = x_sign | (((UINT64) exp + 398) << 53); 99 BID_RETURN (res); 100 } 101 // x is a finite non-zero number (not 0, non-canonical, or special) 102 103 switch (rnd_mode) { 104 case ROUNDING_TO_NEAREST: 105 case ROUNDING_TIES_AWAY: 106 // return 0 if (exp <= -(p+1)) 107 if (exp <= -17) { 108 res = x_sign | 0x31c0000000000000ull; 109 *pfpsf |= INEXACT_EXCEPTION; 110 BID_RETURN (res); 111 } 112 break; 113 case ROUNDING_DOWN: 114 // return 0 if (exp <= -p) 115 if (exp <= -16) { 116 if (x_sign) { 117 res = 0xb1c0000000000001ull; 118 } else { 119 res = 0x31c0000000000000ull; 120 } 121 *pfpsf |= INEXACT_EXCEPTION; 122 BID_RETURN (res); 123 } 124 break; 125 case ROUNDING_UP: 126 // return 0 if (exp <= -p) 127 if (exp <= -16) { 128 if (x_sign) { 129 res = 0xb1c0000000000000ull; 130 } else { 131 res = 0x31c0000000000001ull; 132 } 133 *pfpsf |= INEXACT_EXCEPTION; 134 BID_RETURN (res); 135 } 136 break; 137 case ROUNDING_TO_ZERO: 138 // return 0 if (exp <= -p) 139 if (exp <= -16) { 140 res = x_sign | 0x31c0000000000000ull; 141 *pfpsf |= INEXACT_EXCEPTION; 142 BID_RETURN (res); 143 } 144 break; 145 } // end switch () 146 147 // q = nr. of decimal digits in x (1 <= q <= 54) 148 // determine first the nr. of bits in x 149 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 150 q = 16; 151 } else { // if x < 2^53 152 tmp1.d = (double) C1; // exact conversion 153 x_nr_bits = 154 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 155 q = nr_digits[x_nr_bits - 1].digits; 156 if (q == 0) { 157 q = nr_digits[x_nr_bits - 1].digits1; 158 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 159 q++; 160 } 161 } 162 163 if (exp >= 0) { // -exp <= 0 164 // the argument is an integer already 165 res = x; 166 BID_RETURN (res); 167 } 168 169 switch (rnd_mode) { 170 case ROUNDING_TO_NEAREST: 171 if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q 172 // need to shift right -exp digits from the coefficient; exp will be 0 173 ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' 174 // chop off ind digits from the lower part of C1 175 // C1 = C1 + 1/2 * 10^x where the result C1 fits in 64 bits 176 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate 177 C1 = C1 + midpoint64[ind - 1]; 178 // calculate C* and f* 179 // C* is actually floor(C*) in this case 180 // C* and f* need shifting and masking, as shown by 181 // shiftright128[] and maskhigh128[] 182 // 1 <= x <= 16 183 // kx = 10^(-x) = ten2mk64[ind - 1] 184 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 185 // the approximation of 10^(-x) was rounded up to 64 bits 186 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); 187 188 // if (0 < f* < 10^(-x)) then the result is a midpoint 189 // if floor(C*) is even then C* = floor(C*) - logical right 190 // shift; C* has p decimal digits, correct by Prop. 1) 191 // else if floor(C*) is odd C* = floor(C*)-1 (logical right 192 // shift; C* has p decimal digits, correct by Pr. 1) 193 // else 194 // C* = floor(C*) (logical right shift; C has p decimal digits, 195 // correct by Property 1) 196 // n = C* * 10^(e+x) 197 198 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 199 res = P128.w[1]; 200 fstar.w[1] = 0; 201 fstar.w[0] = P128.w[0]; 202 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 203 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 204 res = (P128.w[1] >> shift); 205 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 206 fstar.w[0] = P128.w[0]; 207 } 208 // if (0 < f* < 10^(-x)) then the result is a midpoint 209 // since round_to_even, subtract 1 if current result is odd 210 if ((res & 0x0000000000000001ull) && (fstar.w[1] == 0) 211 && (fstar.w[0] < ten2mk64[ind - 1])) { 212 res--; 213 } 214 // determine inexactness of the rounding of C* 215 // if (0 < f* - 1/2 < 10^(-x)) then 216 // the result is exact 217 // else // if (f* - 1/2 > T*) then 218 // the result is inexact 219 if (ind - 1 <= 2) { 220 if (fstar.w[0] > 0x8000000000000000ull) { 221 // f* > 1/2 and the result may be exact 222 // fstar.w[0] - 0x8000000000000000ull is f* - 1/2 223 if ((fstar.w[0] - 0x8000000000000000ull) > ten2mk64[ind - 1]) { 224 // set the inexact flag 225 *pfpsf |= INEXACT_EXCEPTION; 226 } // else the result is exact 227 } else { // the result is inexact; f2* <= 1/2 228 // set the inexact flag 229 *pfpsf |= INEXACT_EXCEPTION; 230 } 231 } else { // if 3 <= ind - 1 <= 21 232 if (fstar.w[1] > onehalf128[ind - 1] || 233 (fstar.w[1] == onehalf128[ind - 1] && fstar.w[0])) { 234 // f2* > 1/2 and the result may be exact 235 // Calculate f2* - 1/2 236 if (fstar.w[1] > onehalf128[ind - 1] 237 || fstar.w[0] > ten2mk64[ind - 1]) { 238 // set the inexact flag 239 *pfpsf |= INEXACT_EXCEPTION; 240 } // else the result is exact 241 } else { // the result is inexact; f2* <= 1/2 242 // set the inexact flag 243 *pfpsf |= INEXACT_EXCEPTION; 244 } 245 } 246 // set exponent to zero as it was negative before. 247 res = x_sign | 0x31c0000000000000ull | res; 248 BID_RETURN (res); 249 } else { // if exp < 0 and q + exp < 0 250 // the result is +0 or -0 251 res = x_sign | 0x31c0000000000000ull; 252 *pfpsf |= INEXACT_EXCEPTION; 253 BID_RETURN (res); 254 } 255 break; 256 case ROUNDING_TIES_AWAY: 257 if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q 258 // need to shift right -exp digits from the coefficient; exp will be 0 259 ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' 260 // chop off ind digits from the lower part of C1 261 // C1 = C1 + 1/2 * 10^x where the result C1 fits in 64 bits 262 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate 263 C1 = C1 + midpoint64[ind - 1]; 264 // calculate C* and f* 265 // C* is actually floor(C*) in this case 266 // C* and f* need shifting and masking, as shown by 267 // shiftright128[] and maskhigh128[] 268 // 1 <= x <= 16 269 // kx = 10^(-x) = ten2mk64[ind - 1] 270 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 271 // the approximation of 10^(-x) was rounded up to 64 bits 272 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); 273 274 // if (0 < f* < 10^(-x)) then the result is a midpoint 275 // C* = floor(C*) - logical right shift; C* has p decimal digits, 276 // correct by Prop. 1) 277 // else 278 // C* = floor(C*) (logical right shift; C has p decimal digits, 279 // correct by Property 1) 280 // n = C* * 10^(e+x) 281 282 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 283 res = P128.w[1]; 284 fstar.w[1] = 0; 285 fstar.w[0] = P128.w[0]; 286 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 287 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 288 res = (P128.w[1] >> shift); 289 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 290 fstar.w[0] = P128.w[0]; 291 } 292 // midpoints are already rounded correctly 293 // determine inexactness of the rounding of C* 294 // if (0 < f* - 1/2 < 10^(-x)) then 295 // the result is exact 296 // else // if (f* - 1/2 > T*) then 297 // the result is inexact 298 if (ind - 1 <= 2) { 299 if (fstar.w[0] > 0x8000000000000000ull) { 300 // f* > 1/2 and the result may be exact 301 // fstar.w[0] - 0x8000000000000000ull is f* - 1/2 302 if ((fstar.w[0] - 0x8000000000000000ull) > ten2mk64[ind - 1]) { 303 // set the inexact flag 304 *pfpsf |= INEXACT_EXCEPTION; 305 } // else the result is exact 306 } else { // the result is inexact; f2* <= 1/2 307 // set the inexact flag 308 *pfpsf |= INEXACT_EXCEPTION; 309 } 310 } else { // if 3 <= ind - 1 <= 21 311 if (fstar.w[1] > onehalf128[ind - 1] || 312 (fstar.w[1] == onehalf128[ind - 1] && fstar.w[0])) { 313 // f2* > 1/2 and the result may be exact 314 // Calculate f2* - 1/2 315 if (fstar.w[1] > onehalf128[ind - 1] 316 || fstar.w[0] > ten2mk64[ind - 1]) { 317 // set the inexact flag 318 *pfpsf |= INEXACT_EXCEPTION; 319 } // else the result is exact 320 } else { // the result is inexact; f2* <= 1/2 321 // set the inexact flag 322 *pfpsf |= INEXACT_EXCEPTION; 323 } 324 } 325 // set exponent to zero as it was negative before. 326 res = x_sign | 0x31c0000000000000ull | res; 327 BID_RETURN (res); 328 } else { // if exp < 0 and q + exp < 0 329 // the result is +0 or -0 330 res = x_sign | 0x31c0000000000000ull; 331 *pfpsf |= INEXACT_EXCEPTION; 332 BID_RETURN (res); 333 } 334 break; 335 case ROUNDING_DOWN: 336 if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q 337 // need to shift right -exp digits from the coefficient; exp will be 0 338 ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' 339 // chop off ind digits from the lower part of C1 340 // C1 fits in 64 bits 341 // calculate C* and f* 342 // C* is actually floor(C*) in this case 343 // C* and f* need shifting and masking, as shown by 344 // shiftright128[] and maskhigh128[] 345 // 1 <= x <= 16 346 // kx = 10^(-x) = ten2mk64[ind - 1] 347 // C* = C1 * 10^(-x) 348 // the approximation of 10^(-x) was rounded up to 64 bits 349 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); 350 351 // C* = floor(C*) (logical right shift; C has p decimal digits, 352 // correct by Property 1) 353 // if (0 < f* < 10^(-x)) then the result is exact 354 // n = C* * 10^(e+x) 355 356 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 357 res = P128.w[1]; 358 fstar.w[1] = 0; 359 fstar.w[0] = P128.w[0]; 360 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 361 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 362 res = (P128.w[1] >> shift); 363 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 364 fstar.w[0] = P128.w[0]; 365 } 366 // if (f* > 10^(-x)) then the result is inexact 367 if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1])) { 368 if (x_sign) { 369 // if negative and not exact, increment magnitude 370 res++; 371 } 372 *pfpsf |= INEXACT_EXCEPTION; 373 } 374 // set exponent to zero as it was negative before. 375 res = x_sign | 0x31c0000000000000ull | res; 376 BID_RETURN (res); 377 } else { // if exp < 0 and q + exp <= 0 378 // the result is +0 or -1 379 if (x_sign) { 380 res = 0xb1c0000000000001ull; 381 } else { 382 res = 0x31c0000000000000ull; 383 } 384 *pfpsf |= INEXACT_EXCEPTION; 385 BID_RETURN (res); 386 } 387 break; 388 case ROUNDING_UP: 389 if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q 390 // need to shift right -exp digits from the coefficient; exp will be 0 391 ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' 392 // chop off ind digits from the lower part of C1 393 // C1 fits in 64 bits 394 // calculate C* and f* 395 // C* is actually floor(C*) in this case 396 // C* and f* need shifting and masking, as shown by 397 // shiftright128[] and maskhigh128[] 398 // 1 <= x <= 16 399 // kx = 10^(-x) = ten2mk64[ind - 1] 400 // C* = C1 * 10^(-x) 401 // the approximation of 10^(-x) was rounded up to 64 bits 402 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); 403 404 // C* = floor(C*) (logical right shift; C has p decimal digits, 405 // correct by Property 1) 406 // if (0 < f* < 10^(-x)) then the result is exact 407 // n = C* * 10^(e+x) 408 409 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 410 res = P128.w[1]; 411 fstar.w[1] = 0; 412 fstar.w[0] = P128.w[0]; 413 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 414 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 415 res = (P128.w[1] >> shift); 416 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 417 fstar.w[0] = P128.w[0]; 418 } 419 // if (f* > 10^(-x)) then the result is inexact 420 if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1])) { 421 if (!x_sign) { 422 // if positive and not exact, increment magnitude 423 res++; 424 } 425 *pfpsf |= INEXACT_EXCEPTION; 426 } 427 // set exponent to zero as it was negative before. 428 res = x_sign | 0x31c0000000000000ull | res; 429 BID_RETURN (res); 430 } else { // if exp < 0 and q + exp <= 0 431 // the result is -0 or +1 432 if (x_sign) { 433 res = 0xb1c0000000000000ull; 434 } else { 435 res = 0x31c0000000000001ull; 436 } 437 *pfpsf |= INEXACT_EXCEPTION; 438 BID_RETURN (res); 439 } 440 break; 441 case ROUNDING_TO_ZERO: 442 if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q 443 // need to shift right -exp digits from the coefficient; exp will be 0 444 ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' 445 // chop off ind digits from the lower part of C1 446 // C1 fits in 127 bits 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 <= 16 452 // kx = 10^(-x) = ten2mk64[ind - 1] 453 // C* = C1 * 10^(-x) 454 // the approximation of 10^(-x) was rounded up to 64 bits 455 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); 456 457 // C* = floor(C*) (logical right shift; C has p decimal digits, 458 // correct by Property 1) 459 // if (0 < f* < 10^(-x)) then the result is exact 460 // n = C* * 10^(e+x) 461 462 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 463 res = P128.w[1]; 464 fstar.w[1] = 0; 465 fstar.w[0] = P128.w[0]; 466 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 467 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 468 res = (P128.w[1] >> shift); 469 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 470 fstar.w[0] = P128.w[0]; 471 } 472 // if (f* > 10^(-x)) then the result is inexact 473 if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1])) { 474 *pfpsf |= INEXACT_EXCEPTION; 475 } 476 // set exponent to zero as it was negative before. 477 res = x_sign | 0x31c0000000000000ull | res; 478 BID_RETURN (res); 479 } else { // if exp < 0 and q + exp < 0 480 // the result is +0 or -0 481 res = x_sign | 0x31c0000000000000ull; 482 *pfpsf |= INEXACT_EXCEPTION; 483 BID_RETURN (res); 484 } 485 break; 486 } // end switch () 487 BID_RETURN (res); 488 } 489 490 /***************************************************************************** 491 * BID64_round_integral_nearest_even 492 ****************************************************************************/ 493 494 #if DECIMAL_CALL_BY_REFERENCE 495 void 496 bid64_round_integral_nearest_even (UINT64 * pres, 497 UINT64 * 498 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 499 _EXC_INFO_PARAM) { 500 UINT64 x = *px; 501 #else 502 UINT64 503 bid64_round_integral_nearest_even (UINT64 x _EXC_FLAGS_PARAM 504 _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 505 #endif 506 507 UINT64 res = 0xbaddbaddbaddbaddull; 508 UINT64 x_sign; 509 int exp; // unbiased exponent 510 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) 511 BID_UI64DOUBLE tmp1; 512 int x_nr_bits; 513 int q, ind, shift; 514 UINT64 C1; 515 UINT128 fstar; 516 UINT128 P128; 517 518 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 519 520 // check for NaNs and infinities 521 if ((x & MASK_NAN) == MASK_NAN) { // check for NaN 522 if ((x & 0x0003ffffffffffffull) > 999999999999999ull) 523 x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits 524 else 525 x = x & 0xfe03ffffffffffffull; // clear G6-G12 526 if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN 527 // set invalid flag 528 *pfpsf |= INVALID_EXCEPTION; 529 // return quiet (SNaN) 530 res = x & 0xfdffffffffffffffull; 531 } else { // QNaN 532 res = x; 533 } 534 BID_RETURN (res); 535 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity 536 res = x_sign | 0x7800000000000000ull; 537 BID_RETURN (res); 538 } 539 // unpack x 540 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 541 // if the steering bits are 11 (condition will be 0), then 542 // the exponent is G[0:w+1] 543 exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398; 544 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 545 if (C1 > 9999999999999999ull) { // non-canonical 546 C1 = 0; 547 } 548 } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) 549 exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398; 550 C1 = (x & MASK_BINARY_SIG1); 551 } 552 553 // if x is 0 or non-canonical 554 if (C1 == 0) { 555 if (exp < 0) 556 exp = 0; 557 res = x_sign | (((UINT64) exp + 398) << 53); 558 BID_RETURN (res); 559 } 560 // x is a finite non-zero number (not 0, non-canonical, or special) 561 562 // return 0 if (exp <= -(p+1)) 563 if (exp <= -17) { 564 res = x_sign | 0x31c0000000000000ull; 565 BID_RETURN (res); 566 } 567 // q = nr. of decimal digits in x (1 <= q <= 54) 568 // determine first the nr. of bits in x 569 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 570 q = 16; 571 } else { // if x < 2^53 572 tmp1.d = (double) C1; // exact conversion 573 x_nr_bits = 574 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 575 q = nr_digits[x_nr_bits - 1].digits; 576 if (q == 0) { 577 q = nr_digits[x_nr_bits - 1].digits1; 578 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 579 q++; 580 } 581 } 582 583 if (exp >= 0) { // -exp <= 0 584 // the argument is an integer already 585 res = x; 586 BID_RETURN (res); 587 } else if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q 588 // need to shift right -exp digits from the coefficient; the exp will be 0 589 ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' 590 // chop off ind digits from the lower part of C1 591 // C1 = C1 + 1/2 * 10^x where the result C1 fits in 64 bits 592 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate 593 C1 = C1 + midpoint64[ind - 1]; 594 // calculate C* and f* 595 // C* is actually floor(C*) in this case 596 // C* and f* need shifting and masking, as shown by 597 // shiftright128[] and maskhigh128[] 598 // 1 <= x <= 16 599 // kx = 10^(-x) = ten2mk64[ind - 1] 600 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 601 // the approximation of 10^(-x) was rounded up to 64 bits 602 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); 603 604 // if (0 < f* < 10^(-x)) then the result is a midpoint 605 // if floor(C*) is even then C* = floor(C*) - logical right 606 // shift; C* has p decimal digits, correct by Prop. 1) 607 // else if floor(C*) is odd C* = floor(C*)-1 (logical right 608 // shift; C* has p decimal digits, correct by Pr. 1) 609 // else 610 // C* = floor(C*) (logical right shift; C has p decimal digits, 611 // correct by Property 1) 612 // n = C* * 10^(e+x) 613 614 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 615 res = P128.w[1]; 616 fstar.w[1] = 0; 617 fstar.w[0] = P128.w[0]; 618 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 619 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 620 res = (P128.w[1] >> shift); 621 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 622 fstar.w[0] = P128.w[0]; 623 } 624 // if (0 < f* < 10^(-x)) then the result is a midpoint 625 // since round_to_even, subtract 1 if current result is odd 626 if ((res & 0x0000000000000001ull) && (fstar.w[1] == 0) 627 && (fstar.w[0] < ten2mk64[ind - 1])) { 628 res--; 629 } 630 // set exponent to zero as it was negative before. 631 res = x_sign | 0x31c0000000000000ull | res; 632 BID_RETURN (res); 633 } else { // if exp < 0 and q + exp < 0 634 // the result is +0 or -0 635 res = x_sign | 0x31c0000000000000ull; 636 BID_RETURN (res); 637 } 638 } 639 640 /***************************************************************************** 641 * BID64_round_integral_negative 642 *****************************************************************************/ 643 644 #if DECIMAL_CALL_BY_REFERENCE 645 void 646 bid64_round_integral_negative (UINT64 * pres, 647 UINT64 * 648 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 649 _EXC_INFO_PARAM) { 650 UINT64 x = *px; 651 #else 652 UINT64 653 bid64_round_integral_negative (UINT64 x _EXC_FLAGS_PARAM 654 _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 655 #endif 656 657 UINT64 res = 0xbaddbaddbaddbaddull; 658 UINT64 x_sign; 659 int exp; // unbiased exponent 660 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) 661 BID_UI64DOUBLE tmp1; 662 int x_nr_bits; 663 int q, ind, shift; 664 UINT64 C1; 665 // UINT64 res is C* at first - represents up to 34 decimal digits ~ 113 bits 666 UINT128 fstar; 667 UINT128 P128; 668 669 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 670 671 // check for NaNs and infinities 672 if ((x & MASK_NAN) == MASK_NAN) { // check for NaN 673 if ((x & 0x0003ffffffffffffull) > 999999999999999ull) 674 x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits 675 else 676 x = x & 0xfe03ffffffffffffull; // clear G6-G12 677 if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN 678 // set invalid flag 679 *pfpsf |= INVALID_EXCEPTION; 680 // return quiet (SNaN) 681 res = x & 0xfdffffffffffffffull; 682 } else { // QNaN 683 res = x; 684 } 685 BID_RETURN (res); 686 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity 687 res = x_sign | 0x7800000000000000ull; 688 BID_RETURN (res); 689 } 690 // unpack x 691 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 692 // if the steering bits are 11 (condition will be 0), then 693 // the exponent is G[0:w+1] 694 exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398; 695 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 696 if (C1 > 9999999999999999ull) { // non-canonical 697 C1 = 0; 698 } 699 } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) 700 exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398; 701 C1 = (x & MASK_BINARY_SIG1); 702 } 703 704 // if x is 0 or non-canonical 705 if (C1 == 0) { 706 if (exp < 0) 707 exp = 0; 708 res = x_sign | (((UINT64) exp + 398) << 53); 709 BID_RETURN (res); 710 } 711 // x is a finite non-zero number (not 0, non-canonical, or special) 712 713 // return 0 if (exp <= -p) 714 if (exp <= -16) { 715 if (x_sign) { 716 res = 0xb1c0000000000001ull; 717 } else { 718 res = 0x31c0000000000000ull; 719 } 720 BID_RETURN (res); 721 } 722 // q = nr. of decimal digits in x (1 <= q <= 54) 723 // determine first the nr. of bits in x 724 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 725 q = 16; 726 } else { // if x < 2^53 727 tmp1.d = (double) C1; // exact conversion 728 x_nr_bits = 729 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 730 q = nr_digits[x_nr_bits - 1].digits; 731 if (q == 0) { 732 q = nr_digits[x_nr_bits - 1].digits1; 733 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 734 q++; 735 } 736 } 737 738 if (exp >= 0) { // -exp <= 0 739 // the argument is an integer already 740 res = x; 741 BID_RETURN (res); 742 } else if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q 743 // need to shift right -exp digits from the coefficient; the exp will be 0 744 ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' 745 // chop off ind digits from the lower part of C1 746 // C1 fits in 64 bits 747 // calculate C* and f* 748 // C* is actually floor(C*) in this case 749 // C* and f* need shifting and masking, as shown by 750 // shiftright128[] and maskhigh128[] 751 // 1 <= x <= 16 752 // kx = 10^(-x) = ten2mk64[ind - 1] 753 // C* = C1 * 10^(-x) 754 // the approximation of 10^(-x) was rounded up to 64 bits 755 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); 756 757 // C* = floor(C*) (logical right shift; C has p decimal digits, 758 // correct by Property 1) 759 // if (0 < f* < 10^(-x)) then the result is exact 760 // n = C* * 10^(e+x) 761 762 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 763 res = P128.w[1]; 764 fstar.w[1] = 0; 765 fstar.w[0] = P128.w[0]; 766 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 767 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 768 res = (P128.w[1] >> shift); 769 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 770 fstar.w[0] = P128.w[0]; 771 } 772 // if (f* > 10^(-x)) then the result is inexact 773 if (x_sign 774 && ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1]))) { 775 // if negative and not exact, increment magnitude 776 res++; 777 } 778 // set exponent to zero as it was negative before. 779 res = x_sign | 0x31c0000000000000ull | res; 780 BID_RETURN (res); 781 } else { // if exp < 0 and q + exp <= 0 782 // the result is +0 or -1 783 if (x_sign) { 784 res = 0xb1c0000000000001ull; 785 } else { 786 res = 0x31c0000000000000ull; 787 } 788 BID_RETURN (res); 789 } 790 } 791 792 /***************************************************************************** 793 * BID64_round_integral_positive 794 ****************************************************************************/ 795 796 #if DECIMAL_CALL_BY_REFERENCE 797 void 798 bid64_round_integral_positive (UINT64 * pres, 799 UINT64 * 800 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 801 _EXC_INFO_PARAM) { 802 UINT64 x = *px; 803 #else 804 UINT64 805 bid64_round_integral_positive (UINT64 x _EXC_FLAGS_PARAM 806 _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 807 #endif 808 809 UINT64 res = 0xbaddbaddbaddbaddull; 810 UINT64 x_sign; 811 int exp; // unbiased exponent 812 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) 813 BID_UI64DOUBLE tmp1; 814 int x_nr_bits; 815 int q, ind, shift; 816 UINT64 C1; 817 // UINT64 res is C* at first - represents up to 34 decimal digits ~ 113 bits 818 UINT128 fstar; 819 UINT128 P128; 820 821 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 822 823 // check for NaNs and infinities 824 if ((x & MASK_NAN) == MASK_NAN) { // check for NaN 825 if ((x & 0x0003ffffffffffffull) > 999999999999999ull) 826 x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits 827 else 828 x = x & 0xfe03ffffffffffffull; // clear G6-G12 829 if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN 830 // set invalid flag 831 *pfpsf |= INVALID_EXCEPTION; 832 // return quiet (SNaN) 833 res = x & 0xfdffffffffffffffull; 834 } else { // QNaN 835 res = x; 836 } 837 BID_RETURN (res); 838 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity 839 res = x_sign | 0x7800000000000000ull; 840 BID_RETURN (res); 841 } 842 // unpack x 843 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 844 // if the steering bits are 11 (condition will be 0), then 845 // the exponent is G[0:w+1] 846 exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398; 847 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 848 if (C1 > 9999999999999999ull) { // non-canonical 849 C1 = 0; 850 } 851 } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) 852 exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398; 853 C1 = (x & MASK_BINARY_SIG1); 854 } 855 856 // if x is 0 or non-canonical 857 if (C1 == 0) { 858 if (exp < 0) 859 exp = 0; 860 res = x_sign | (((UINT64) exp + 398) << 53); 861 BID_RETURN (res); 862 } 863 // x is a finite non-zero number (not 0, non-canonical, or special) 864 865 // return 0 if (exp <= -p) 866 if (exp <= -16) { 867 if (x_sign) { 868 res = 0xb1c0000000000000ull; 869 } else { 870 res = 0x31c0000000000001ull; 871 } 872 BID_RETURN (res); 873 } 874 // q = nr. of decimal digits in x (1 <= q <= 54) 875 // determine first the nr. of bits in x 876 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 877 q = 16; 878 } else { // if x < 2^53 879 tmp1.d = (double) C1; // exact conversion 880 x_nr_bits = 881 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 882 q = nr_digits[x_nr_bits - 1].digits; 883 if (q == 0) { 884 q = nr_digits[x_nr_bits - 1].digits1; 885 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 886 q++; 887 } 888 } 889 890 if (exp >= 0) { // -exp <= 0 891 // the argument is an integer already 892 res = x; 893 BID_RETURN (res); 894 } else if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q 895 // need to shift right -exp digits from the coefficient; the exp will be 0 896 ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' 897 // chop off ind digits from the lower part of C1 898 // C1 fits in 64 bits 899 // calculate C* and f* 900 // C* is actually floor(C*) in this case 901 // C* and f* need shifting and masking, as shown by 902 // shiftright128[] and maskhigh128[] 903 // 1 <= x <= 16 904 // kx = 10^(-x) = ten2mk64[ind - 1] 905 // C* = C1 * 10^(-x) 906 // the approximation of 10^(-x) was rounded up to 64 bits 907 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); 908 909 // C* = floor(C*) (logical right shift; C has p decimal digits, 910 // correct by Property 1) 911 // if (0 < f* < 10^(-x)) then the result is exact 912 // n = C* * 10^(e+x) 913 914 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 915 res = P128.w[1]; 916 fstar.w[1] = 0; 917 fstar.w[0] = P128.w[0]; 918 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 919 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 920 res = (P128.w[1] >> shift); 921 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 922 fstar.w[0] = P128.w[0]; 923 } 924 // if (f* > 10^(-x)) then the result is inexact 925 if (!x_sign 926 && ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind - 1]))) { 927 // if positive and not exact, increment magnitude 928 res++; 929 } 930 // set exponent to zero as it was negative before. 931 res = x_sign | 0x31c0000000000000ull | res; 932 BID_RETURN (res); 933 } else { // if exp < 0 and q + exp <= 0 934 // the result is -0 or +1 935 if (x_sign) { 936 res = 0xb1c0000000000000ull; 937 } else { 938 res = 0x31c0000000000001ull; 939 } 940 BID_RETURN (res); 941 } 942 } 943 944 /***************************************************************************** 945 * BID64_round_integral_zero 946 ****************************************************************************/ 947 948 #if DECIMAL_CALL_BY_REFERENCE 949 void 950 bid64_round_integral_zero (UINT64 * pres, 951 UINT64 * 952 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 953 _EXC_INFO_PARAM) { 954 UINT64 x = *px; 955 #else 956 UINT64 957 bid64_round_integral_zero (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 958 _EXC_INFO_PARAM) { 959 #endif 960 961 UINT64 res = 0xbaddbaddbaddbaddull; 962 UINT64 x_sign; 963 int exp; // unbiased exponent 964 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) 965 BID_UI64DOUBLE tmp1; 966 int x_nr_bits; 967 int q, ind, shift; 968 UINT64 C1; 969 // UINT64 res is C* at first - represents up to 34 decimal digits ~ 113 bits 970 UINT128 P128; 971 972 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 973 974 // check for NaNs and infinities 975 if ((x & MASK_NAN) == MASK_NAN) { // check for NaN 976 if ((x & 0x0003ffffffffffffull) > 999999999999999ull) 977 x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits 978 else 979 x = x & 0xfe03ffffffffffffull; // clear G6-G12 980 if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN 981 // set invalid flag 982 *pfpsf |= INVALID_EXCEPTION; 983 // return quiet (SNaN) 984 res = x & 0xfdffffffffffffffull; 985 } else { // QNaN 986 res = x; 987 } 988 BID_RETURN (res); 989 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity 990 res = x_sign | 0x7800000000000000ull; 991 BID_RETURN (res); 992 } 993 // unpack x 994 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 995 // if the steering bits are 11 (condition will be 0), then 996 // the exponent is G[0:w+1] 997 exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398; 998 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 999 if (C1 > 9999999999999999ull) { // non-canonical 1000 C1 = 0; 1001 } 1002 } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) 1003 exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398; 1004 C1 = (x & MASK_BINARY_SIG1); 1005 } 1006 1007 // if x is 0 or non-canonical 1008 if (C1 == 0) { 1009 if (exp < 0) 1010 exp = 0; 1011 res = x_sign | (((UINT64) exp + 398) << 53); 1012 BID_RETURN (res); 1013 } 1014 // x is a finite non-zero number (not 0, non-canonical, or special) 1015 1016 // return 0 if (exp <= -p) 1017 if (exp <= -16) { 1018 res = x_sign | 0x31c0000000000000ull; 1019 BID_RETURN (res); 1020 } 1021 // q = nr. of decimal digits in x (1 <= q <= 54) 1022 // determine first the nr. of bits in x 1023 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 1024 q = 16; 1025 } else { // if x < 2^53 1026 tmp1.d = (double) C1; // exact conversion 1027 x_nr_bits = 1028 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1029 q = nr_digits[x_nr_bits - 1].digits; 1030 if (q == 0) { 1031 q = nr_digits[x_nr_bits - 1].digits1; 1032 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 1033 q++; 1034 } 1035 } 1036 1037 if (exp >= 0) { // -exp <= 0 1038 // the argument is an integer already 1039 res = x; 1040 BID_RETURN (res); 1041 } else if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q 1042 // need to shift right -exp digits from the coefficient; the exp will be 0 1043 ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' 1044 // chop off ind digits from the lower part of C1 1045 // C1 fits in 127 bits 1046 // calculate C* and f* 1047 // C* is actually floor(C*) in this case 1048 // C* and f* need shifting and masking, as shown by 1049 // shiftright128[] and maskhigh128[] 1050 // 1 <= x <= 16 1051 // kx = 10^(-x) = ten2mk64[ind - 1] 1052 // C* = C1 * 10^(-x) 1053 // the approximation of 10^(-x) was rounded up to 64 bits 1054 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); 1055 1056 // C* = floor(C*) (logical right shift; C has p decimal digits, 1057 // correct by Property 1) 1058 // if (0 < f* < 10^(-x)) then the result is exact 1059 // n = C* * 10^(e+x) 1060 1061 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 1062 res = P128.w[1]; 1063 // redundant fstar.w[1] = 0; 1064 // redundant fstar.w[0] = P128.w[0]; 1065 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 1066 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 1067 res = (P128.w[1] >> shift); 1068 // redundant fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; 1069 // redundant fstar.w[0] = P128.w[0]; 1070 } 1071 // if (f* > 10^(-x)) then the result is inexact 1072 // if ((fstar.w[1] != 0) || (fstar.w[0] >= ten2mk64[ind-1])){ 1073 // // redundant 1074 // } 1075 // set exponent to zero as it was negative before. 1076 res = x_sign | 0x31c0000000000000ull | res; 1077 BID_RETURN (res); 1078 } else { // if exp < 0 and q + exp < 0 1079 // the result is +0 or -0 1080 res = x_sign | 0x31c0000000000000ull; 1081 BID_RETURN (res); 1082 } 1083 } 1084 1085 /***************************************************************************** 1086 * BID64_round_integral_nearest_away 1087 ****************************************************************************/ 1088 1089 #if DECIMAL_CALL_BY_REFERENCE 1090 void 1091 bid64_round_integral_nearest_away (UINT64 * pres, 1092 UINT64 * 1093 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 1094 _EXC_INFO_PARAM) { 1095 UINT64 x = *px; 1096 #else 1097 UINT64 1098 bid64_round_integral_nearest_away (UINT64 x _EXC_FLAGS_PARAM 1099 _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 1100 #endif 1101 1102 UINT64 res = 0xbaddbaddbaddbaddull; 1103 UINT64 x_sign; 1104 int exp; // unbiased exponent 1105 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) 1106 BID_UI64DOUBLE tmp1; 1107 int x_nr_bits; 1108 int q, ind, shift; 1109 UINT64 C1; 1110 UINT128 P128; 1111 1112 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 1113 1114 // check for NaNs and infinities 1115 if ((x & MASK_NAN) == MASK_NAN) { // check for NaN 1116 if ((x & 0x0003ffffffffffffull) > 999999999999999ull) 1117 x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits 1118 else 1119 x = x & 0xfe03ffffffffffffull; // clear G6-G12 1120 if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN 1121 // set invalid flag 1122 *pfpsf |= INVALID_EXCEPTION; 1123 // return quiet (SNaN) 1124 res = x & 0xfdffffffffffffffull; 1125 } else { // QNaN 1126 res = x; 1127 } 1128 BID_RETURN (res); 1129 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity 1130 res = x_sign | 0x7800000000000000ull; 1131 BID_RETURN (res); 1132 } 1133 // unpack x 1134 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 1135 // if the steering bits are 11 (condition will be 0), then 1136 // the exponent is G[0:w+1] 1137 exp = ((x & MASK_BINARY_EXPONENT2) >> 51) - 398; 1138 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 1139 if (C1 > 9999999999999999ull) { // non-canonical 1140 C1 = 0; 1141 } 1142 } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) 1143 exp = ((x & MASK_BINARY_EXPONENT1) >> 53) - 398; 1144 C1 = (x & MASK_BINARY_SIG1); 1145 } 1146 1147 // if x is 0 or non-canonical 1148 if (C1 == 0) { 1149 if (exp < 0) 1150 exp = 0; 1151 res = x_sign | (((UINT64) exp + 398) << 53); 1152 BID_RETURN (res); 1153 } 1154 // x is a finite non-zero number (not 0, non-canonical, or special) 1155 1156 // return 0 if (exp <= -(p+1)) 1157 if (exp <= -17) { 1158 res = x_sign | 0x31c0000000000000ull; 1159 BID_RETURN (res); 1160 } 1161 // q = nr. of decimal digits in x (1 <= q <= 54) 1162 // determine first the nr. of bits in x 1163 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 1164 q = 16; 1165 } else { // if x < 2^53 1166 tmp1.d = (double) C1; // exact conversion 1167 x_nr_bits = 1168 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1169 q = nr_digits[x_nr_bits - 1].digits; 1170 if (q == 0) { 1171 q = nr_digits[x_nr_bits - 1].digits1; 1172 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 1173 q++; 1174 } 1175 } 1176 1177 if (exp >= 0) { // -exp <= 0 1178 // the argument is an integer already 1179 res = x; 1180 BID_RETURN (res); 1181 } else if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q 1182 // need to shift right -exp digits from the coefficient; the exp will be 0 1183 ind = -exp; // 1 <= ind <= 16; ind is a synonym for 'x' 1184 // chop off ind digits from the lower part of C1 1185 // C1 = C1 + 1/2 * 10^x where the result C1 fits in 64 bits 1186 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate 1187 C1 = C1 + midpoint64[ind - 1]; 1188 // calculate C* and f* 1189 // C* is actually floor(C*) in this case 1190 // C* and f* need shifting and masking, as shown by 1191 // shiftright128[] and maskhigh128[] 1192 // 1 <= x <= 16 1193 // kx = 10^(-x) = ten2mk64[ind - 1] 1194 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 1195 // the approximation of 10^(-x) was rounded up to 64 bits 1196 __mul_64x64_to_128 (P128, C1, ten2mk64[ind - 1]); 1197 1198 // if (0 < f* < 10^(-x)) then the result is a midpoint 1199 // C* = floor(C*) - logical right shift; C* has p decimal digits, 1200 // correct by Prop. 1) 1201 // else 1202 // C* = floor(C*) (logical right shift; C has p decimal digits, 1203 // correct by Property 1) 1204 // n = C* * 10^(e+x) 1205 1206 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 1207 res = P128.w[1]; 1208 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 1209 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 1210 res = (P128.w[1] >> shift); 1211 } 1212 // midpoints are already rounded correctly 1213 // set exponent to zero as it was negative before. 1214 res = x_sign | 0x31c0000000000000ull | res; 1215 BID_RETURN (res); 1216 } else { // if exp < 0 and q + exp < 0 1217 // the result is +0 or -0 1218 res = x_sign | 0x31c0000000000000ull; 1219 BID_RETURN (res); 1220 } 1221 } 1222