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_round_integral_exact 28 ****************************************************************************/ 29 30#if DECIMAL_CALL_BY_REFERENCE 31void 32bid64_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 41UINT64 42bid64_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 495void 496bid64_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 502UINT64 503bid64_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 645void 646bid64_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 652UINT64 653bid64_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 797void 798bid64_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 804UINT64 805bid64_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 949void 950bid64_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 956UINT64 957bid64_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 1090void 1091bid64_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 1097UINT64 1098bid64_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