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