1 /* Copyright (C) 2007-2020 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 #define BID_128RES 25 26 #include "bid_internal.h" 27 28 /***************************************************************************** 29 * BID128_round_integral_exact 30 ****************************************************************************/ 31 32 BID128_FUNCTION_ARG1 (bid128_round_integral_exact, x) 33 34 UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} 35 }; 36 UINT64 x_sign; 37 UINT64 x_exp; 38 int exp; // unbiased exponent 39 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) 40 UINT64 tmp64; 41 BID_UI64DOUBLE tmp1; 42 unsigned int x_nr_bits; 43 int q, ind, shift; 44 UINT128 C1; 45 UINT256 fstar; 46 UINT256 P256; 47 48 // check for NaN or Infinity 49 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 50 // x is special 51 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN 52 // if x = NaN, then res = Q (x) 53 // check first for non-canonical NaN payload 54 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) || 55 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) && 56 (x.w[0] > 0x38c15b09ffffffffull))) { 57 x.w[1] = x.w[1] & 0xffffc00000000000ull; 58 x.w[0] = 0x0ull; 59 } 60 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN 61 // set invalid flag 62 *pfpsf |= INVALID_EXCEPTION; 63 // return quiet (x) 64 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16] 65 res.w[0] = x.w[0]; 66 } else { // x is QNaN 67 // return x 68 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16] 69 res.w[0] = x.w[0]; 70 } 71 BID_RETURN (res) 72 } else { // x is not a NaN, so it must be infinity 73 if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf 74 // return +inf 75 res.w[1] = 0x7800000000000000ull; 76 res.w[0] = 0x0000000000000000ull; 77 } else { // x is -inf 78 // return -inf 79 res.w[1] = 0xf800000000000000ull; 80 res.w[0] = 0x0000000000000000ull; 81 } 82 BID_RETURN (res); 83 } 84 } 85 // unpack x 86 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 87 C1.w[1] = x.w[1] & MASK_COEFF; 88 C1.w[0] = x.w[0]; 89 90 // check for non-canonical values (treated as zero) 91 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11 92 // non-canonical 93 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits 94 C1.w[1] = 0; // significand high 95 C1.w[0] = 0; // significand low 96 } else { // G0_G1 != 11 97 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits 98 if (C1.w[1] > 0x0001ed09bead87c0ull || 99 (C1.w[1] == 0x0001ed09bead87c0ull 100 && C1.w[0] > 0x378d8e63ffffffffull)) { 101 // x is non-canonical if coefficient is larger than 10^34 -1 102 C1.w[1] = 0; 103 C1.w[0] = 0; 104 } else { // canonical 105 ; 106 } 107 } 108 109 // test for input equal to zero 110 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { 111 // x is 0 112 // return 0 preserving the sign bit and the preferred exponent 113 // of MAX(Q(x), 0) 114 if (x_exp <= (0x1820ull << 49)) { 115 res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull; 116 } else { 117 res.w[1] = x_sign | x_exp; 118 } 119 res.w[0] = 0x0000000000000000ull; 120 BID_RETURN (res); 121 } 122 // x is not special and is not zero 123 124 switch (rnd_mode) { 125 case ROUNDING_TO_NEAREST: 126 case ROUNDING_TIES_AWAY: 127 // if (exp <= -(p+1)) return 0.0 128 if (x_exp <= 0x2ffa000000000000ull) { // 0x2ffa000000000000ull == -35 129 res.w[1] = x_sign | 0x3040000000000000ull; 130 res.w[0] = 0x0000000000000000ull; 131 *pfpsf |= INEXACT_EXCEPTION; 132 BID_RETURN (res); 133 } 134 break; 135 case ROUNDING_DOWN: 136 // if (exp <= -p) return -1.0 or +0.0 137 if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffa000000000000ull == -34 138 if (x_sign) { 139 // if negative, return negative 1, because we know coefficient 140 // is non-zero (would have been caught above) 141 res.w[1] = 0xb040000000000000ull; 142 res.w[0] = 0x0000000000000001ull; 143 } else { 144 // if positive, return positive 0, because we know coefficient is 145 // non-zero (would have been caught above) 146 res.w[1] = 0x3040000000000000ull; 147 res.w[0] = 0x0000000000000000ull; 148 } 149 *pfpsf |= INEXACT_EXCEPTION; 150 BID_RETURN (res); 151 } 152 break; 153 case ROUNDING_UP: 154 // if (exp <= -p) return -0.0 or +1.0 155 if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34 156 if (x_sign) { 157 // if negative, return negative 0, because we know the coefficient 158 // is non-zero (would have been caught above) 159 res.w[1] = 0xb040000000000000ull; 160 res.w[0] = 0x0000000000000000ull; 161 } else { 162 // if positive, return positive 1, because we know coefficient is 163 // non-zero (would have been caught above) 164 res.w[1] = 0x3040000000000000ull; 165 res.w[0] = 0x0000000000000001ull; 166 } 167 *pfpsf |= INEXACT_EXCEPTION; 168 BID_RETURN (res); 169 } 170 break; 171 case ROUNDING_TO_ZERO: 172 // if (exp <= -p) return -0.0 or +0.0 173 if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34 174 res.w[1] = x_sign | 0x3040000000000000ull; 175 res.w[0] = 0x0000000000000000ull; 176 *pfpsf |= INEXACT_EXCEPTION; 177 BID_RETURN (res); 178 } 179 break; 180 } 181 182 // q = nr. of decimal digits in x 183 // determine first the nr. of bits in x 184 if (C1.w[1] == 0) { 185 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 186 // split the 64-bit value in two 32-bit halves to avoid rounding errors 187 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 188 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion 189 x_nr_bits = 190 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 191 } else { // x < 2^32 192 tmp1.d = (double) (C1.w[0]); // exact conversion 193 x_nr_bits = 194 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 195 } 196 } else { // if x < 2^53 197 tmp1.d = (double) C1.w[0]; // exact conversion 198 x_nr_bits = 199 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 200 } 201 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) 202 tmp1.d = (double) C1.w[1]; // exact conversion 203 x_nr_bits = 204 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 205 } 206 207 q = nr_digits[x_nr_bits - 1].digits; 208 if (q == 0) { 209 q = nr_digits[x_nr_bits - 1].digits1; 210 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi || 211 (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi && 212 C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) 213 q++; 214 } 215 exp = (x_exp >> 49) - 6176; 216 if (exp >= 0) { // -exp <= 0 217 // the argument is an integer already 218 res.w[1] = x.w[1]; 219 res.w[0] = x.w[0]; 220 BID_RETURN (res); 221 } 222 // exp < 0 223 switch (rnd_mode) { 224 case ROUNDING_TO_NEAREST: 225 if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q 226 // need to shift right -exp digits from the coefficient; exp will be 0 227 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x' 228 // chop off ind digits from the lower part of C1 229 // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits 230 tmp64 = C1.w[0]; 231 if (ind <= 19) { 232 C1.w[0] = C1.w[0] + midpoint64[ind - 1]; 233 } else { 234 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; 235 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; 236 } 237 if (C1.w[0] < tmp64) 238 C1.w[1]++; 239 // calculate C* and f* 240 // C* is actually floor(C*) in this case 241 // C* and f* need shifting and masking, as shown by 242 // shiftright128[] and maskhigh128[] 243 // 1 <= x <= 34 244 // kx = 10^(-x) = ten2mk128[ind - 1] 245 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 246 // the approximation of 10^(-x) was rounded up to 118 bits 247 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); 248 // determine the value of res and fstar 249 250 // determine inexactness of the rounding of C* 251 // if (0 < f* - 1/2 < 10^(-x)) then 252 // the result is exact 253 // else // if (f* - 1/2 > T*) then 254 // the result is inexact 255 // Note: we are going to use ten2mk128[] instead of ten2mk128trunc[] 256 257 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 258 // redundant shift = shiftright128[ind - 1]; // shift = 0 259 res.w[1] = P256.w[3]; 260 res.w[0] = P256.w[2]; 261 // redundant fstar.w[3] = 0; 262 // redundant fstar.w[2] = 0; 263 fstar.w[1] = P256.w[1]; 264 fstar.w[0] = P256.w[0]; 265 // fraction f* < 10^(-x) <=> midpoint 266 // f* is in the right position to be compared with 267 // 10^(-x) from ten2mk128[] 268 // if 0 < fstar < 10^(-x), subtract 1 if odd (for rounding to even) 269 if ((res.w[0] & 0x0000000000000001ull) && // is result odd? 270 ((fstar.w[1] < (ten2mk128[ind - 1].w[1])) 271 || ((fstar.w[1] == ten2mk128[ind - 1].w[1]) 272 && (fstar.w[0] < ten2mk128[ind - 1].w[0])))) { 273 // subract 1 to make even 274 if (res.w[0]-- == 0) { 275 res.w[1]--; 276 } 277 } 278 if (fstar.w[1] > 0x8000000000000000ull || 279 (fstar.w[1] == 0x8000000000000000ull 280 && fstar.w[0] > 0x0ull)) { 281 // f* > 1/2 and the result may be exact 282 tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2 283 if (tmp64 > ten2mk128[ind - 1].w[1] || 284 (tmp64 == ten2mk128[ind - 1].w[1] && 285 fstar.w[0] >= ten2mk128[ind - 1].w[0])) { 286 // set the inexact flag 287 *pfpsf |= INEXACT_EXCEPTION; 288 } // else the result is exact 289 } else { // the result is inexact; f2* <= 1/2 290 // set the inexact flag 291 *pfpsf |= INEXACT_EXCEPTION; 292 } 293 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 294 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 295 res.w[1] = (P256.w[3] >> shift); 296 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift); 297 // redundant fstar.w[3] = 0; 298 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; 299 fstar.w[1] = P256.w[1]; 300 fstar.w[0] = P256.w[0]; 301 // fraction f* < 10^(-x) <=> midpoint 302 // f* is in the right position to be compared with 303 // 10^(-x) from ten2mk128[] 304 if ((res.w[0] & 0x0000000000000001ull) && // is result odd? 305 fstar.w[2] == 0 && (fstar.w[1] < ten2mk128[ind - 1].w[1] || 306 (fstar.w[1] == ten2mk128[ind - 1].w[1] && 307 fstar.w[0] < ten2mk128[ind - 1].w[0]))) { 308 // subract 1 to make even 309 if (res.w[0]-- == 0) { 310 res.w[1]--; 311 } 312 } 313 if (fstar.w[2] > onehalf128[ind - 1] || 314 (fstar.w[2] == onehalf128[ind - 1] 315 && (fstar.w[1] || fstar.w[0]))) { 316 // f2* > 1/2 and the result may be exact 317 // Calculate f2* - 1/2 318 tmp64 = fstar.w[2] - onehalf128[ind - 1]; 319 if (tmp64 || fstar.w[1] > ten2mk128[ind - 1].w[1] || 320 (fstar.w[1] == ten2mk128[ind - 1].w[1] && 321 fstar.w[0] >= ten2mk128[ind - 1].w[0])) { 322 // set the inexact flag 323 *pfpsf |= INEXACT_EXCEPTION; 324 } // else the result is exact 325 } else { // the result is inexact; f2* <= 1/2 326 // set the inexact flag 327 *pfpsf |= INEXACT_EXCEPTION; 328 } 329 } else { // 22 <= ind - 1 <= 33 330 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38 331 res.w[1] = 0; 332 res.w[0] = P256.w[3] >> shift; 333 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; 334 fstar.w[2] = P256.w[2]; 335 fstar.w[1] = P256.w[1]; 336 fstar.w[0] = P256.w[0]; 337 // fraction f* < 10^(-x) <=> midpoint 338 // f* is in the right position to be compared with 339 // 10^(-x) from ten2mk128[] 340 if ((res.w[0] & 0x0000000000000001ull) && // is result odd? 341 fstar.w[3] == 0 && fstar.w[2] == 0 && 342 (fstar.w[1] < ten2mk128[ind - 1].w[1] || 343 (fstar.w[1] == ten2mk128[ind - 1].w[1] && 344 fstar.w[0] < ten2mk128[ind - 1].w[0]))) { 345 // subract 1 to make even 346 if (res.w[0]-- == 0) { 347 res.w[1]--; 348 } 349 } 350 if (fstar.w[3] > onehalf128[ind - 1] || 351 (fstar.w[3] == onehalf128[ind - 1] && 352 (fstar.w[2] || fstar.w[1] || fstar.w[0]))) { 353 // f2* > 1/2 and the result may be exact 354 // Calculate f2* - 1/2 355 tmp64 = fstar.w[3] - onehalf128[ind - 1]; 356 if (tmp64 || fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] 357 || (fstar.w[1] == ten2mk128[ind - 1].w[1] 358 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) { 359 // set the inexact flag 360 *pfpsf |= INEXACT_EXCEPTION; 361 } // else the result is exact 362 } else { // the result is inexact; f2* <= 1/2 363 // set the inexact flag 364 *pfpsf |= INEXACT_EXCEPTION; 365 } 366 } 367 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1]; 368 BID_RETURN (res); 369 } else { // if ((q + exp) < 0) <=> q < -exp 370 // the result is +0 or -0 371 res.w[1] = x_sign | 0x3040000000000000ull; 372 res.w[0] = 0x0000000000000000ull; 373 *pfpsf |= INEXACT_EXCEPTION; 374 BID_RETURN (res); 375 } 376 break; 377 case ROUNDING_TIES_AWAY: 378 if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q 379 // need to shift right -exp digits from the coefficient; exp will be 0 380 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x' 381 // chop off ind digits from the lower part of C1 382 // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits 383 tmp64 = C1.w[0]; 384 if (ind <= 19) { 385 C1.w[0] = C1.w[0] + midpoint64[ind - 1]; 386 } else { 387 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; 388 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; 389 } 390 if (C1.w[0] < tmp64) 391 C1.w[1]++; 392 // calculate C* and f* 393 // C* is actually floor(C*) in this case 394 // C* and f* need shifting and masking, as shown by 395 // shiftright128[] and maskhigh128[] 396 // 1 <= x <= 34 397 // kx = 10^(-x) = ten2mk128[ind - 1] 398 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 399 // the approximation of 10^(-x) was rounded up to 118 bits 400 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); 401 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g. 402 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999 403 // if (0 < f* < 10^(-x)) then the result is a midpoint 404 // if floor(C*) is even then C* = floor(C*) - logical right 405 // shift; C* has p decimal digits, correct by Prop. 1) 406 // else if floor(C*) is odd C* = floor(C*)-1 (logical right 407 // shift; C* has p decimal digits, correct by Pr. 1) 408 // else 409 // C* = floor(C*) (logical right shift; C has p decimal digits, 410 // correct by Property 1) 411 // n = C* * 10^(e+x) 412 413 // determine also the inexactness of the rounding of C* 414 // if (0 < f* - 1/2 < 10^(-x)) then 415 // the result is exact 416 // else // if (f* - 1/2 > T*) then 417 // the result is inexact 418 // Note: we are going to use ten2mk128[] instead of ten2mk128trunc[] 419 // shift right C* by Ex-128 = shiftright128[ind] 420 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 421 // redundant shift = shiftright128[ind - 1]; // shift = 0 422 res.w[1] = P256.w[3]; 423 res.w[0] = P256.w[2]; 424 // redundant fstar.w[3] = 0; 425 // redundant fstar.w[2] = 0; 426 fstar.w[1] = P256.w[1]; 427 fstar.w[0] = P256.w[0]; 428 if (fstar.w[1] > 0x8000000000000000ull || 429 (fstar.w[1] == 0x8000000000000000ull 430 && fstar.w[0] > 0x0ull)) { 431 // f* > 1/2 and the result may be exact 432 tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2 433 if ((tmp64 > ten2mk128[ind - 1].w[1] || 434 (tmp64 == ten2mk128[ind - 1].w[1] && 435 fstar.w[0] >= ten2mk128[ind - 1].w[0]))) { 436 // set the inexact flag 437 *pfpsf |= INEXACT_EXCEPTION; 438 } // else the result is exact 439 } else { // the result is inexact; f2* <= 1/2 440 // set the inexact flag 441 *pfpsf |= INEXACT_EXCEPTION; 442 } 443 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 444 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 445 res.w[1] = (P256.w[3] >> shift); 446 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift); 447 // redundant fstar.w[3] = 0; 448 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; 449 fstar.w[1] = P256.w[1]; 450 fstar.w[0] = P256.w[0]; 451 if (fstar.w[2] > onehalf128[ind - 1] || 452 (fstar.w[2] == onehalf128[ind - 1] 453 && (fstar.w[1] || fstar.w[0]))) { 454 // f2* > 1/2 and the result may be exact 455 // Calculate f2* - 1/2 456 tmp64 = fstar.w[2] - onehalf128[ind - 1]; 457 if (tmp64 || fstar.w[1] > ten2mk128[ind - 1].w[1] || 458 (fstar.w[1] == ten2mk128[ind - 1].w[1] && 459 fstar.w[0] >= ten2mk128[ind - 1].w[0])) { 460 // set the inexact flag 461 *pfpsf |= INEXACT_EXCEPTION; 462 } // else the result is exact 463 } else { // the result is inexact; f2* <= 1/2 464 // set the inexact flag 465 *pfpsf |= INEXACT_EXCEPTION; 466 } 467 } else { // 22 <= ind - 1 <= 33 468 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38 469 res.w[1] = 0; 470 res.w[0] = P256.w[3] >> shift; 471 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; 472 fstar.w[2] = P256.w[2]; 473 fstar.w[1] = P256.w[1]; 474 fstar.w[0] = P256.w[0]; 475 if (fstar.w[3] > onehalf128[ind - 1] || 476 (fstar.w[3] == onehalf128[ind - 1] && 477 (fstar.w[2] || fstar.w[1] || fstar.w[0]))) { 478 // f2* > 1/2 and the result may be exact 479 // Calculate f2* - 1/2 480 tmp64 = fstar.w[3] - onehalf128[ind - 1]; 481 if (tmp64 || fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] 482 || (fstar.w[1] == ten2mk128[ind - 1].w[1] 483 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) { 484 // set the inexact flag 485 *pfpsf |= INEXACT_EXCEPTION; 486 } // else the result is exact 487 } else { // the result is inexact; f2* <= 1/2 488 // set the inexact flag 489 *pfpsf |= INEXACT_EXCEPTION; 490 } 491 } 492 // if the result was a midpoint, it was already rounded away from zero 493 res.w[1] |= x_sign | 0x3040000000000000ull; 494 BID_RETURN (res); 495 } else { // if ((q + exp) < 0) <=> q < -exp 496 // the result is +0 or -0 497 res.w[1] = x_sign | 0x3040000000000000ull; 498 res.w[0] = 0x0000000000000000ull; 499 *pfpsf |= INEXACT_EXCEPTION; 500 BID_RETURN (res); 501 } 502 break; 503 case ROUNDING_DOWN: 504 if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q 505 // need to shift right -exp digits from the coefficient; exp will be 0 506 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x' 507 // (number of digits to be chopped off) 508 // chop off ind digits from the lower part of C1 509 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate 510 // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP 511 // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE 512 // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE 513 // tmp64 = C1.w[0]; 514 // if (ind <= 19) { 515 // C1.w[0] = C1.w[0] + midpoint64[ind - 1]; 516 // } else { 517 // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; 518 // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; 519 // } 520 // if (C1.w[0] < tmp64) C1.w[1]++; 521 // if carry-out from C1.w[0], increment C1.w[1] 522 // calculate C* and f* 523 // C* is actually floor(C*) in this case 524 // C* and f* need shifting and masking, as shown by 525 // shiftright128[] and maskhigh128[] 526 // 1 <= x <= 34 527 // kx = 10^(-x) = ten2mk128[ind - 1] 528 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 529 // the approximation of 10^(-x) was rounded up to 118 bits 530 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); 531 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 532 res.w[1] = P256.w[3]; 533 res.w[0] = P256.w[2]; 534 // redundant fstar.w[3] = 0; 535 // redundant fstar.w[2] = 0; 536 // redundant fstar.w[1] = P256.w[1]; 537 // redundant fstar.w[0] = P256.w[0]; 538 // fraction f* > 10^(-x) <=> inexact 539 // f* is in the right position to be compared with 540 // 10^(-x) from ten2mk128[] 541 if ((P256.w[1] > ten2mk128[ind - 1].w[1]) 542 || (P256.w[1] == ten2mk128[ind - 1].w[1] 543 && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) { 544 *pfpsf |= INEXACT_EXCEPTION; 545 // if positive, the truncated value is already the correct result 546 if (x_sign) { // if negative 547 if (++res.w[0] == 0) { 548 res.w[1]++; 549 } 550 } 551 } 552 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 553 shift = shiftright128[ind - 1]; // 0 <= shift <= 102 554 res.w[1] = (P256.w[3] >> shift); 555 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift); 556 // redundant fstar.w[3] = 0; 557 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; 558 fstar.w[1] = P256.w[1]; 559 fstar.w[0] = P256.w[0]; 560 // fraction f* > 10^(-x) <=> inexact 561 // f* is in the right position to be compared with 562 // 10^(-x) from ten2mk128[] 563 if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] || 564 (fstar.w[1] == ten2mk128[ind - 1].w[1] && 565 fstar.w[0] >= ten2mk128[ind - 1].w[0])) { 566 *pfpsf |= INEXACT_EXCEPTION; 567 // if positive, the truncated value is already the correct result 568 if (x_sign) { // if negative 569 if (++res.w[0] == 0) { 570 res.w[1]++; 571 } 572 } 573 } 574 } else { // 22 <= ind - 1 <= 33 575 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38 576 res.w[1] = 0; 577 res.w[0] = P256.w[3] >> shift; 578 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; 579 fstar.w[2] = P256.w[2]; 580 fstar.w[1] = P256.w[1]; 581 fstar.w[0] = P256.w[0]; 582 // fraction f* > 10^(-x) <=> inexact 583 // f* is in the right position to be compared with 584 // 10^(-x) from ten2mk128[] 585 if (fstar.w[3] || fstar.w[2] 586 || fstar.w[1] > ten2mk128[ind - 1].w[1] 587 || (fstar.w[1] == ten2mk128[ind - 1].w[1] 588 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) { 589 *pfpsf |= INEXACT_EXCEPTION; 590 // if positive, the truncated value is already the correct result 591 if (x_sign) { // if negative 592 if (++res.w[0] == 0) { 593 res.w[1]++; 594 } 595 } 596 } 597 } 598 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1]; 599 BID_RETURN (res); 600 } else { // if exp < 0 and q + exp <= 0 601 if (x_sign) { // negative rounds down to -1.0 602 res.w[1] = 0xb040000000000000ull; 603 res.w[0] = 0x0000000000000001ull; 604 } else { // positive rpunds down to +0.0 605 res.w[1] = 0x3040000000000000ull; 606 res.w[0] = 0x0000000000000000ull; 607 } 608 *pfpsf |= INEXACT_EXCEPTION; 609 BID_RETURN (res); 610 } 611 break; 612 case ROUNDING_UP: 613 if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q 614 // need to shift right -exp digits from the coefficient; exp will be 0 615 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x' 616 // (number of digits to be chopped off) 617 // chop off ind digits from the lower part of C1 618 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate 619 // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP 620 // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE 621 // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE 622 // tmp64 = C1.w[0]; 623 // if (ind <= 19) { 624 // C1.w[0] = C1.w[0] + midpoint64[ind - 1]; 625 // } else { 626 // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; 627 // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; 628 // } 629 // if (C1.w[0] < tmp64) C1.w[1]++; 630 // if carry-out from C1.w[0], increment C1.w[1] 631 // calculate C* and f* 632 // C* is actually floor(C*) in this case 633 // C* and f* need shifting and masking, as shown by 634 // shiftright128[] and maskhigh128[] 635 // 1 <= x <= 34 636 // kx = 10^(-x) = ten2mk128[ind - 1] 637 // C* = C1 * 10^(-x) 638 // the approximation of 10^(-x) was rounded up to 118 bits 639 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); 640 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 641 res.w[1] = P256.w[3]; 642 res.w[0] = P256.w[2]; 643 // redundant fstar.w[3] = 0; 644 // redundant fstar.w[2] = 0; 645 // redundant fstar.w[1] = P256.w[1]; 646 // redundant fstar.w[0] = P256.w[0]; 647 // fraction f* > 10^(-x) <=> inexact 648 // f* is in the right position to be compared with 649 // 10^(-x) from ten2mk128[] 650 if ((P256.w[1] > ten2mk128[ind - 1].w[1]) 651 || (P256.w[1] == ten2mk128[ind - 1].w[1] 652 && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) { 653 *pfpsf |= INEXACT_EXCEPTION; 654 // if negative, the truncated value is already the correct result 655 if (!x_sign) { // if positive 656 if (++res.w[0] == 0) { 657 res.w[1]++; 658 } 659 } 660 } 661 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 662 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 663 res.w[1] = (P256.w[3] >> shift); 664 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift); 665 // redundant fstar.w[3] = 0; 666 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; 667 fstar.w[1] = P256.w[1]; 668 fstar.w[0] = P256.w[0]; 669 // fraction f* > 10^(-x) <=> inexact 670 // f* is in the right position to be compared with 671 // 10^(-x) from ten2mk128[] 672 if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] || 673 (fstar.w[1] == ten2mk128[ind - 1].w[1] && 674 fstar.w[0] >= ten2mk128[ind - 1].w[0])) { 675 *pfpsf |= INEXACT_EXCEPTION; 676 // if negative, the truncated value is already the correct result 677 if (!x_sign) { // if positive 678 if (++res.w[0] == 0) { 679 res.w[1]++; 680 } 681 } 682 } 683 } else { // 22 <= ind - 1 <= 33 684 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38 685 res.w[1] = 0; 686 res.w[0] = P256.w[3] >> shift; 687 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; 688 fstar.w[2] = P256.w[2]; 689 fstar.w[1] = P256.w[1]; 690 fstar.w[0] = P256.w[0]; 691 // fraction f* > 10^(-x) <=> inexact 692 // f* is in the right position to be compared with 693 // 10^(-x) from ten2mk128[] 694 if (fstar.w[3] || fstar.w[2] 695 || fstar.w[1] > ten2mk128[ind - 1].w[1] 696 || (fstar.w[1] == ten2mk128[ind - 1].w[1] 697 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) { 698 *pfpsf |= INEXACT_EXCEPTION; 699 // if negative, the truncated value is already the correct result 700 if (!x_sign) { // if positive 701 if (++res.w[0] == 0) { 702 res.w[1]++; 703 } 704 } 705 } 706 } 707 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1]; 708 BID_RETURN (res); 709 } else { // if exp < 0 and q + exp <= 0 710 if (x_sign) { // negative rounds up to -0.0 711 res.w[1] = 0xb040000000000000ull; 712 res.w[0] = 0x0000000000000000ull; 713 } else { // positive rpunds up to +1.0 714 res.w[1] = 0x3040000000000000ull; 715 res.w[0] = 0x0000000000000001ull; 716 } 717 *pfpsf |= INEXACT_EXCEPTION; 718 BID_RETURN (res); 719 } 720 break; 721 case ROUNDING_TO_ZERO: 722 if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q 723 // need to shift right -exp digits from the coefficient; exp will be 0 724 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x' 725 // (number of digits to be chopped off) 726 // chop off ind digits from the lower part of C1 727 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate 728 // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP 729 // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE 730 // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE 731 //tmp64 = C1.w[0]; 732 // if (ind <= 19) { 733 // C1.w[0] = C1.w[0] + midpoint64[ind - 1]; 734 // } else { 735 // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; 736 // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; 737 // } 738 // if (C1.w[0] < tmp64) C1.w[1]++; 739 // if carry-out from C1.w[0], increment C1.w[1] 740 // calculate C* and f* 741 // C* is actually floor(C*) in this case 742 // C* and f* need shifting and masking, as shown by 743 // shiftright128[] and maskhigh128[] 744 // 1 <= x <= 34 745 // kx = 10^(-x) = ten2mk128[ind - 1] 746 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 747 // the approximation of 10^(-x) was rounded up to 118 bits 748 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); 749 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 750 res.w[1] = P256.w[3]; 751 res.w[0] = P256.w[2]; 752 // redundant fstar.w[3] = 0; 753 // redundant fstar.w[2] = 0; 754 // redundant fstar.w[1] = P256.w[1]; 755 // redundant fstar.w[0] = P256.w[0]; 756 // fraction f* > 10^(-x) <=> inexact 757 // f* is in the right position to be compared with 758 // 10^(-x) from ten2mk128[] 759 if ((P256.w[1] > ten2mk128[ind - 1].w[1]) 760 || (P256.w[1] == ten2mk128[ind - 1].w[1] 761 && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) { 762 *pfpsf |= INEXACT_EXCEPTION; 763 } 764 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 765 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 766 res.w[1] = (P256.w[3] >> shift); 767 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift); 768 // redundant fstar.w[3] = 0; 769 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; 770 fstar.w[1] = P256.w[1]; 771 fstar.w[0] = P256.w[0]; 772 // fraction f* > 10^(-x) <=> inexact 773 // f* is in the right position to be compared with 774 // 10^(-x) from ten2mk128[] 775 if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] || 776 (fstar.w[1] == ten2mk128[ind - 1].w[1] && 777 fstar.w[0] >= ten2mk128[ind - 1].w[0])) { 778 *pfpsf |= INEXACT_EXCEPTION; 779 } 780 } else { // 22 <= ind - 1 <= 33 781 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38 782 res.w[1] = 0; 783 res.w[0] = P256.w[3] >> shift; 784 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; 785 fstar.w[2] = P256.w[2]; 786 fstar.w[1] = P256.w[1]; 787 fstar.w[0] = P256.w[0]; 788 // fraction f* > 10^(-x) <=> inexact 789 // f* is in the right position to be compared with 790 // 10^(-x) from ten2mk128[] 791 if (fstar.w[3] || fstar.w[2] 792 || fstar.w[1] > ten2mk128[ind - 1].w[1] 793 || (fstar.w[1] == ten2mk128[ind - 1].w[1] 794 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) { 795 *pfpsf |= INEXACT_EXCEPTION; 796 } 797 } 798 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1]; 799 BID_RETURN (res); 800 } else { // if exp < 0 and q + exp <= 0 the result is +0 or -0 801 res.w[1] = x_sign | 0x3040000000000000ull; 802 res.w[0] = 0x0000000000000000ull; 803 *pfpsf |= INEXACT_EXCEPTION; 804 BID_RETURN (res); 805 } 806 break; 807 } 808 809 BID_RETURN (res); 810 } 811 812 /***************************************************************************** 813 * BID128_round_integral_nearest_even 814 ****************************************************************************/ 815 816 BID128_FUNCTION_ARG1_NORND (bid128_round_integral_nearest_even, x) 817 818 UINT128 res; 819 UINT64 x_sign; 820 UINT64 x_exp; 821 int exp; // unbiased exponent 822 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) 823 UINT64 tmp64; 824 BID_UI64DOUBLE tmp1; 825 unsigned int x_nr_bits; 826 int q, ind, shift; 827 UINT128 C1; 828 // UINT128 res is C* at first - represents up to 34 decimal digits ~ 113 bits 829 UINT256 fstar; 830 UINT256 P256; 831 832 // check for NaN or Infinity 833 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 834 // x is special 835 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN 836 // if x = NaN, then res = Q (x) 837 // check first for non-canonical NaN payload 838 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) || 839 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) && 840 (x.w[0] > 0x38c15b09ffffffffull))) { 841 x.w[1] = x.w[1] & 0xffffc00000000000ull; 842 x.w[0] = 0x0ull; 843 } 844 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN 845 // set invalid flag 846 *pfpsf |= INVALID_EXCEPTION; 847 // return quiet (x) 848 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16] 849 res.w[0] = x.w[0]; 850 } else { // x is QNaN 851 // return x 852 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16] 853 res.w[0] = x.w[0]; 854 } 855 BID_RETURN (res) 856 } else { // x is not a NaN, so it must be infinity 857 if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf 858 // return +inf 859 res.w[1] = 0x7800000000000000ull; 860 res.w[0] = 0x0000000000000000ull; 861 } else { // x is -inf 862 // return -inf 863 res.w[1] = 0xf800000000000000ull; 864 res.w[0] = 0x0000000000000000ull; 865 } 866 BID_RETURN (res); 867 } 868 } 869 // unpack x 870 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 871 C1.w[1] = x.w[1] & MASK_COEFF; 872 C1.w[0] = x.w[0]; 873 874 // check for non-canonical values (treated as zero) 875 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11 876 // non-canonical 877 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits 878 C1.w[1] = 0; // significand high 879 C1.w[0] = 0; // significand low 880 } else { // G0_G1 != 11 881 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits 882 if (C1.w[1] > 0x0001ed09bead87c0ull || 883 (C1.w[1] == 0x0001ed09bead87c0ull 884 && C1.w[0] > 0x378d8e63ffffffffull)) { 885 // x is non-canonical if coefficient is larger than 10^34 -1 886 C1.w[1] = 0; 887 C1.w[0] = 0; 888 } else { // canonical 889 ; 890 } 891 } 892 893 // test for input equal to zero 894 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { 895 // x is 0 896 // return 0 preserving the sign bit and the preferred exponent 897 // of MAX(Q(x), 0) 898 if (x_exp <= (0x1820ull << 49)) { 899 res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull; 900 } else { 901 res.w[1] = x_sign | x_exp; 902 } 903 res.w[0] = 0x0000000000000000ull; 904 BID_RETURN (res); 905 } 906 // x is not special and is not zero 907 908 // if (exp <= -(p+1)) return 0 909 if (x_exp <= 0x2ffa000000000000ull) { // 0x2ffa000000000000ull == -35 910 res.w[1] = x_sign | 0x3040000000000000ull; 911 res.w[0] = 0x0000000000000000ull; 912 BID_RETURN (res); 913 } 914 // q = nr. of decimal digits in x 915 // determine first the nr. of bits in x 916 if (C1.w[1] == 0) { 917 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 918 // split the 64-bit value in two 32-bit halves to avoid rounding errors 919 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 920 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion 921 x_nr_bits = 922 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 923 } else { // x < 2^32 924 tmp1.d = (double) (C1.w[0]); // exact conversion 925 x_nr_bits = 926 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 927 } 928 } else { // if x < 2^53 929 tmp1.d = (double) C1.w[0]; // exact conversion 930 x_nr_bits = 931 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 932 } 933 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) 934 tmp1.d = (double) C1.w[1]; // exact conversion 935 x_nr_bits = 936 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 937 } 938 939 q = nr_digits[x_nr_bits - 1].digits; 940 if (q == 0) { 941 q = nr_digits[x_nr_bits - 1].digits1; 942 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi 943 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi && 944 C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) 945 q++; 946 } 947 exp = (x_exp >> 49) - 6176; 948 if (exp >= 0) { // -exp <= 0 949 // the argument is an integer already 950 res.w[1] = x.w[1]; 951 res.w[0] = x.w[0]; 952 BID_RETURN (res); 953 } else if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q 954 // need to shift right -exp digits from the coefficient; the exp will be 0 955 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x' 956 // chop off ind digits from the lower part of C1 957 // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits 958 tmp64 = C1.w[0]; 959 if (ind <= 19) { 960 C1.w[0] = C1.w[0] + midpoint64[ind - 1]; 961 } else { 962 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; 963 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; 964 } 965 if (C1.w[0] < tmp64) 966 C1.w[1]++; 967 // calculate C* and f* 968 // C* is actually floor(C*) in this case 969 // C* and f* need shifting and masking, as shown by 970 // shiftright128[] and maskhigh128[] 971 // 1 <= x <= 34 972 // kx = 10^(-x) = ten2mk128[ind - 1] 973 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 974 // the approximation of 10^(-x) was rounded up to 118 bits 975 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); 976 // determine the value of res and fstar 977 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 978 // redundant shift = shiftright128[ind - 1]; // shift = 0 979 res.w[1] = P256.w[3]; 980 res.w[0] = P256.w[2]; 981 // redundant fstar.w[3] = 0; 982 // redundant fstar.w[2] = 0; 983 // redundant fstar.w[1] = P256.w[1]; 984 // redundant fstar.w[0] = P256.w[0]; 985 // fraction f* < 10^(-x) <=> midpoint 986 // f* is in the right position to be compared with 987 // 10^(-x) from ten2mk128[] 988 // if 0 < fstar < 10^(-x), subtract 1 if odd (for rounding to even) 989 if ((res.w[0] & 0x0000000000000001ull) && // is result odd? 990 ((P256.w[1] < (ten2mk128[ind - 1].w[1])) 991 || ((P256.w[1] == ten2mk128[ind - 1].w[1]) 992 && (P256.w[0] < ten2mk128[ind - 1].w[0])))) { 993 // subract 1 to make even 994 if (res.w[0]-- == 0) { 995 res.w[1]--; 996 } 997 } 998 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 999 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 1000 res.w[1] = (P256.w[3] >> shift); 1001 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift); 1002 // redundant fstar.w[3] = 0; 1003 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; 1004 fstar.w[1] = P256.w[1]; 1005 fstar.w[0] = P256.w[0]; 1006 // fraction f* < 10^(-x) <=> midpoint 1007 // f* is in the right position to be compared with 1008 // 10^(-x) from ten2mk128[] 1009 if ((res.w[0] & 0x0000000000000001ull) && // is result odd? 1010 fstar.w[2] == 0 && (fstar.w[1] < ten2mk128[ind - 1].w[1] || 1011 (fstar.w[1] == ten2mk128[ind - 1].w[1] && 1012 fstar.w[0] < ten2mk128[ind - 1].w[0]))) { 1013 // subract 1 to make even 1014 if (res.w[0]-- == 0) { 1015 res.w[1]--; 1016 } 1017 } 1018 } else { // 22 <= ind - 1 <= 33 1019 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38 1020 res.w[1] = 0; 1021 res.w[0] = P256.w[3] >> shift; 1022 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; 1023 fstar.w[2] = P256.w[2]; 1024 fstar.w[1] = P256.w[1]; 1025 fstar.w[0] = P256.w[0]; 1026 // fraction f* < 10^(-x) <=> midpoint 1027 // f* is in the right position to be compared with 1028 // 10^(-x) from ten2mk128[] 1029 if ((res.w[0] & 0x0000000000000001ull) && // is result odd? 1030 fstar.w[3] == 0 && fstar.w[2] == 0 1031 && (fstar.w[1] < ten2mk128[ind - 1].w[1] 1032 || (fstar.w[1] == ten2mk128[ind - 1].w[1] 1033 && fstar.w[0] < ten2mk128[ind - 1].w[0]))) { 1034 // subract 1 to make even 1035 if (res.w[0]-- == 0) { 1036 res.w[1]--; 1037 } 1038 } 1039 } 1040 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1]; 1041 BID_RETURN (res); 1042 } else { // if ((q + exp) < 0) <=> q < -exp 1043 // the result is +0 or -0 1044 res.w[1] = x_sign | 0x3040000000000000ull; 1045 res.w[0] = 0x0000000000000000ull; 1046 BID_RETURN (res); 1047 } 1048 } 1049 1050 /***************************************************************************** 1051 * BID128_round_integral_negative 1052 ****************************************************************************/ 1053 1054 BID128_FUNCTION_ARG1_NORND (bid128_round_integral_negative, x) 1055 1056 UINT128 res; 1057 UINT64 x_sign; 1058 UINT64 x_exp; 1059 int exp; // unbiased exponent 1060 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo 1061 // (all are UINT64) 1062 BID_UI64DOUBLE tmp1; 1063 unsigned int x_nr_bits; 1064 int q, ind, shift; 1065 UINT128 C1; 1066 // UINT128 res is C* at first - represents up to 34 decimal digits ~ 1067 // 113 bits 1068 UINT256 fstar; 1069 UINT256 P256; 1070 1071 // check for NaN or Infinity 1072 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 1073 // x is special 1074 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN 1075 // if x = NaN, then res = Q (x) 1076 // check first for non-canonical NaN payload 1077 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) || 1078 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) && 1079 (x.w[0] > 0x38c15b09ffffffffull))) { 1080 x.w[1] = x.w[1] & 0xffffc00000000000ull; 1081 x.w[0] = 0x0ull; 1082 } 1083 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN 1084 // set invalid flag 1085 *pfpsf |= INVALID_EXCEPTION; 1086 // return quiet (x) 1087 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16] 1088 res.w[0] = x.w[0]; 1089 } else { // x is QNaN 1090 // return x 1091 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16] 1092 res.w[0] = x.w[0]; 1093 } 1094 BID_RETURN (res) 1095 } else { // x is not a NaN, so it must be infinity 1096 if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf 1097 // return +inf 1098 res.w[1] = 0x7800000000000000ull; 1099 res.w[0] = 0x0000000000000000ull; 1100 } else { // x is -inf 1101 // return -inf 1102 res.w[1] = 0xf800000000000000ull; 1103 res.w[0] = 0x0000000000000000ull; 1104 } 1105 BID_RETURN (res); 1106 } 1107 } 1108 // unpack x 1109 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 1110 C1.w[1] = x.w[1] & MASK_COEFF; 1111 C1.w[0] = x.w[0]; 1112 1113 // check for non-canonical values (treated as zero) 1114 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11 1115 // non-canonical 1116 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits 1117 C1.w[1] = 0; // significand high 1118 C1.w[0] = 0; // significand low 1119 } else { // G0_G1 != 11 1120 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits 1121 if (C1.w[1] > 0x0001ed09bead87c0ull || 1122 (C1.w[1] == 0x0001ed09bead87c0ull 1123 && C1.w[0] > 0x378d8e63ffffffffull)) { 1124 // x is non-canonical if coefficient is larger than 10^34 -1 1125 C1.w[1] = 0; 1126 C1.w[0] = 0; 1127 } else { // canonical 1128 ; 1129 } 1130 } 1131 1132 // test for input equal to zero 1133 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { 1134 // x is 0 1135 // return 0 preserving the sign bit and the preferred exponent 1136 // of MAX(Q(x), 0) 1137 if (x_exp <= (0x1820ull << 49)) { 1138 res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull; 1139 } else { 1140 res.w[1] = x_sign | x_exp; 1141 } 1142 res.w[0] = 0x0000000000000000ull; 1143 BID_RETURN (res); 1144 } 1145 // x is not special and is not zero 1146 1147 // if (exp <= -p) return -1.0 or +0.0 1148 if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34 1149 if (x_sign) { 1150 // if negative, return negative 1, because we know the coefficient 1151 // is non-zero (would have been caught above) 1152 res.w[1] = 0xb040000000000000ull; 1153 res.w[0] = 0x0000000000000001ull; 1154 } else { 1155 // if positive, return positive 0, because we know coefficient is 1156 // non-zero (would have been caught above) 1157 res.w[1] = 0x3040000000000000ull; 1158 res.w[0] = 0x0000000000000000ull; 1159 } 1160 BID_RETURN (res); 1161 } 1162 // q = nr. of decimal digits in x 1163 // determine first the nr. of bits in x 1164 if (C1.w[1] == 0) { 1165 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 1166 // split the 64-bit value in two 32-bit halves to avoid rounding errors 1167 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 1168 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion 1169 x_nr_bits = 1170 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1171 } else { // x < 2^32 1172 tmp1.d = (double) (C1.w[0]); // exact conversion 1173 x_nr_bits = 1174 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1175 } 1176 } else { // if x < 2^53 1177 tmp1.d = (double) C1.w[0]; // exact conversion 1178 x_nr_bits = 1179 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1180 } 1181 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) 1182 tmp1.d = (double) C1.w[1]; // exact conversion 1183 x_nr_bits = 1184 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1185 } 1186 1187 q = nr_digits[x_nr_bits - 1].digits; 1188 if (q == 0) { 1189 q = nr_digits[x_nr_bits - 1].digits1; 1190 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi || 1191 (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi && 1192 C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) 1193 q++; 1194 } 1195 exp = (x_exp >> 49) - 6176; 1196 if (exp >= 0) { // -exp <= 0 1197 // the argument is an integer already 1198 res.w[1] = x.w[1]; 1199 res.w[0] = x.w[0]; 1200 BID_RETURN (res); 1201 } else if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q 1202 // need to shift right -exp digits from the coefficient; the exp will be 0 1203 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x' 1204 // (number of digits to be chopped off) 1205 // chop off ind digits from the lower part of C1 1206 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate 1207 // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP 1208 // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE 1209 // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE 1210 //tmp64 = C1.w[0]; 1211 // if (ind <= 19) { 1212 // C1.w[0] = C1.w[0] + midpoint64[ind - 1]; 1213 // } else { 1214 // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; 1215 // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; 1216 // } 1217 // if (C1.w[0] < tmp64) C1.w[1]++; 1218 // if carry-out from C1.w[0], increment C1.w[1] 1219 // calculate C* and f* 1220 // C* is actually floor(C*) in this case 1221 // C* and f* need shifting and masking, as shown by 1222 // shiftright128[] and maskhigh128[] 1223 // 1 <= x <= 34 1224 // kx = 10^(-x) = ten2mk128[ind - 1] 1225 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 1226 // the approximation of 10^(-x) was rounded up to 118 bits 1227 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); 1228 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 1229 res.w[1] = P256.w[3]; 1230 res.w[0] = P256.w[2]; 1231 // if positive, the truncated value is already the correct result 1232 if (x_sign) { // if negative 1233 // redundant fstar.w[3] = 0; 1234 // redundant fstar.w[2] = 0; 1235 // redundant fstar.w[1] = P256.w[1]; 1236 // redundant fstar.w[0] = P256.w[0]; 1237 // fraction f* > 10^(-x) <=> inexact 1238 // f* is in the right position to be compared with 1239 // 10^(-x) from ten2mk128[] 1240 if ((P256.w[1] > ten2mk128[ind - 1].w[1]) 1241 || (P256.w[1] == ten2mk128[ind - 1].w[1] 1242 && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) { 1243 if (++res.w[0] == 0) { 1244 res.w[1]++; 1245 } 1246 } 1247 } 1248 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 1249 shift = shiftright128[ind - 1]; // 0 <= shift <= 102 1250 res.w[1] = (P256.w[3] >> shift); 1251 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift); 1252 // if positive, the truncated value is already the correct result 1253 if (x_sign) { // if negative 1254 // redundant fstar.w[3] = 0; 1255 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; 1256 fstar.w[1] = P256.w[1]; 1257 fstar.w[0] = P256.w[0]; 1258 // fraction f* > 10^(-x) <=> inexact 1259 // f* is in the right position to be compared with 1260 // 10^(-x) from ten2mk128[] 1261 if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] || 1262 (fstar.w[1] == ten2mk128[ind - 1].w[1] && 1263 fstar.w[0] >= ten2mk128[ind - 1].w[0])) { 1264 if (++res.w[0] == 0) { 1265 res.w[1]++; 1266 } 1267 } 1268 } 1269 } else { // 22 <= ind - 1 <= 33 1270 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38 1271 res.w[1] = 0; 1272 res.w[0] = P256.w[3] >> shift; 1273 // if positive, the truncated value is already the correct result 1274 if (x_sign) { // if negative 1275 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; 1276 fstar.w[2] = P256.w[2]; 1277 fstar.w[1] = P256.w[1]; 1278 fstar.w[0] = P256.w[0]; 1279 // fraction f* > 10^(-x) <=> inexact 1280 // f* is in the right position to be compared with 1281 // 10^(-x) from ten2mk128[] 1282 if (fstar.w[3] || fstar.w[2] 1283 || fstar.w[1] > ten2mk128[ind - 1].w[1] 1284 || (fstar.w[1] == ten2mk128[ind - 1].w[1] 1285 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) { 1286 if (++res.w[0] == 0) { 1287 res.w[1]++; 1288 } 1289 } 1290 } 1291 } 1292 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1]; 1293 BID_RETURN (res); 1294 } else { // if exp < 0 and q + exp <= 0 1295 if (x_sign) { // negative rounds down to -1.0 1296 res.w[1] = 0xb040000000000000ull; 1297 res.w[0] = 0x0000000000000001ull; 1298 } else { // positive rpunds down to +0.0 1299 res.w[1] = 0x3040000000000000ull; 1300 res.w[0] = 0x0000000000000000ull; 1301 } 1302 BID_RETURN (res); 1303 } 1304 } 1305 1306 /***************************************************************************** 1307 * BID128_round_integral_positive 1308 ****************************************************************************/ 1309 1310 BID128_FUNCTION_ARG1_NORND (bid128_round_integral_positive, x) 1311 1312 UINT128 res; 1313 UINT64 x_sign; 1314 UINT64 x_exp; 1315 int exp; // unbiased exponent 1316 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo 1317 // (all are UINT64) 1318 BID_UI64DOUBLE tmp1; 1319 unsigned int x_nr_bits; 1320 int q, ind, shift; 1321 UINT128 C1; 1322 // UINT128 res is C* at first - represents up to 34 decimal digits ~ 1323 // 113 bits 1324 UINT256 fstar; 1325 UINT256 P256; 1326 1327 // check for NaN or Infinity 1328 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 1329 // x is special 1330 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN 1331 // if x = NaN, then res = Q (x) 1332 // check first for non-canonical NaN payload 1333 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) || 1334 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) && 1335 (x.w[0] > 0x38c15b09ffffffffull))) { 1336 x.w[1] = x.w[1] & 0xffffc00000000000ull; 1337 x.w[0] = 0x0ull; 1338 } 1339 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN 1340 // set invalid flag 1341 *pfpsf |= INVALID_EXCEPTION; 1342 // return quiet (x) 1343 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16] 1344 res.w[0] = x.w[0]; 1345 } else { // x is QNaN 1346 // return x 1347 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16] 1348 res.w[0] = x.w[0]; 1349 } 1350 BID_RETURN (res) 1351 } else { // x is not a NaN, so it must be infinity 1352 if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf 1353 // return +inf 1354 res.w[1] = 0x7800000000000000ull; 1355 res.w[0] = 0x0000000000000000ull; 1356 } else { // x is -inf 1357 // return -inf 1358 res.w[1] = 0xf800000000000000ull; 1359 res.w[0] = 0x0000000000000000ull; 1360 } 1361 BID_RETURN (res); 1362 } 1363 } 1364 // unpack x 1365 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 1366 C1.w[1] = x.w[1] & MASK_COEFF; 1367 C1.w[0] = x.w[0]; 1368 1369 // check for non-canonical values (treated as zero) 1370 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11 1371 // non-canonical 1372 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits 1373 C1.w[1] = 0; // significand high 1374 C1.w[0] = 0; // significand low 1375 } else { // G0_G1 != 11 1376 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits 1377 if (C1.w[1] > 0x0001ed09bead87c0ull || 1378 (C1.w[1] == 0x0001ed09bead87c0ull 1379 && C1.w[0] > 0x378d8e63ffffffffull)) { 1380 // x is non-canonical if coefficient is larger than 10^34 -1 1381 C1.w[1] = 0; 1382 C1.w[0] = 0; 1383 } else { // canonical 1384 ; 1385 } 1386 } 1387 1388 // test for input equal to zero 1389 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { 1390 // x is 0 1391 // return 0 preserving the sign bit and the preferred exponent 1392 // of MAX(Q(x), 0) 1393 if (x_exp <= (0x1820ull << 49)) { 1394 res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull; 1395 } else { 1396 res.w[1] = x_sign | x_exp; 1397 } 1398 res.w[0] = 0x0000000000000000ull; 1399 BID_RETURN (res); 1400 } 1401 // x is not special and is not zero 1402 1403 // if (exp <= -p) return -0.0 or +1.0 1404 if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34 1405 if (x_sign) { 1406 // if negative, return negative 0, because we know the coefficient 1407 // is non-zero (would have been caught above) 1408 res.w[1] = 0xb040000000000000ull; 1409 res.w[0] = 0x0000000000000000ull; 1410 } else { 1411 // if positive, return positive 1, because we know coefficient is 1412 // non-zero (would have been caught above) 1413 res.w[1] = 0x3040000000000000ull; 1414 res.w[0] = 0x0000000000000001ull; 1415 } 1416 BID_RETURN (res); 1417 } 1418 // q = nr. of decimal digits in x 1419 // determine first the nr. of bits in x 1420 if (C1.w[1] == 0) { 1421 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 1422 // split 64-bit value in two 32-bit halves to avoid rounding errors 1423 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 1424 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion 1425 x_nr_bits = 1426 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1427 } else { // x < 2^32 1428 tmp1.d = (double) (C1.w[0]); // exact conversion 1429 x_nr_bits = 1430 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1431 } 1432 } else { // if x < 2^53 1433 tmp1.d = (double) C1.w[0]; // exact conversion 1434 x_nr_bits = 1435 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1436 } 1437 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) 1438 tmp1.d = (double) C1.w[1]; // exact conversion 1439 x_nr_bits = 1440 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1441 } 1442 1443 q = nr_digits[x_nr_bits - 1].digits; 1444 if (q == 0) { 1445 q = nr_digits[x_nr_bits - 1].digits1; 1446 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi || 1447 (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi && 1448 C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) 1449 q++; 1450 } 1451 exp = (x_exp >> 49) - 6176; 1452 if (exp >= 0) { // -exp <= 0 1453 // the argument is an integer already 1454 res.w[1] = x.w[1]; 1455 res.w[0] = x.w[0]; 1456 BID_RETURN (res); 1457 } else if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q 1458 // need to shift right -exp digits from the coefficient; exp will be 0 1459 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x' 1460 // (number of digits to be chopped off) 1461 // chop off ind digits from the lower part of C1 1462 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate 1463 // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP 1464 // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE 1465 // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE 1466 // tmp64 = C1.w[0]; 1467 // if (ind <= 19) { 1468 // C1.w[0] = C1.w[0] + midpoint64[ind - 1]; 1469 // } else { 1470 // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; 1471 // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; 1472 // } 1473 // if (C1.w[0] < tmp64) C1.w[1]++; 1474 // if carry-out from C1.w[0], increment C1.w[1] 1475 // calculate C* and f* 1476 // C* is actually floor(C*) in this case 1477 // C* and f* need shifting and masking, as shown by 1478 // shiftright128[] and maskhigh128[] 1479 // 1 <= x <= 34 1480 // kx = 10^(-x) = ten2mk128[ind - 1] 1481 // C* = C1 * 10^(-x) 1482 // the approximation of 10^(-x) was rounded up to 118 bits 1483 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); 1484 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 1485 res.w[1] = P256.w[3]; 1486 res.w[0] = P256.w[2]; 1487 // if negative, the truncated value is already the correct result 1488 if (!x_sign) { // if positive 1489 // redundant fstar.w[3] = 0; 1490 // redundant fstar.w[2] = 0; 1491 // redundant fstar.w[1] = P256.w[1]; 1492 // redundant fstar.w[0] = P256.w[0]; 1493 // fraction f* > 10^(-x) <=> inexact 1494 // f* is in the right position to be compared with 1495 // 10^(-x) from ten2mk128[] 1496 if ((P256.w[1] > ten2mk128[ind - 1].w[1]) 1497 || (P256.w[1] == ten2mk128[ind - 1].w[1] 1498 && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) { 1499 if (++res.w[0] == 0) { 1500 res.w[1]++; 1501 } 1502 } 1503 } 1504 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 1505 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 1506 res.w[1] = (P256.w[3] >> shift); 1507 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift); 1508 // if negative, the truncated value is already the correct result 1509 if (!x_sign) { // if positive 1510 // redundant fstar.w[3] = 0; 1511 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; 1512 fstar.w[1] = P256.w[1]; 1513 fstar.w[0] = P256.w[0]; 1514 // fraction f* > 10^(-x) <=> inexact 1515 // f* is in the right position to be compared with 1516 // 10^(-x) from ten2mk128[] 1517 if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] || 1518 (fstar.w[1] == ten2mk128[ind - 1].w[1] && 1519 fstar.w[0] >= ten2mk128[ind - 1].w[0])) { 1520 if (++res.w[0] == 0) { 1521 res.w[1]++; 1522 } 1523 } 1524 } 1525 } else { // 22 <= ind - 1 <= 33 1526 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38 1527 res.w[1] = 0; 1528 res.w[0] = P256.w[3] >> shift; 1529 // if negative, the truncated value is already the correct result 1530 if (!x_sign) { // if positive 1531 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; 1532 fstar.w[2] = P256.w[2]; 1533 fstar.w[1] = P256.w[1]; 1534 fstar.w[0] = P256.w[0]; 1535 // fraction f* > 10^(-x) <=> inexact 1536 // f* is in the right position to be compared with 1537 // 10^(-x) from ten2mk128[] 1538 if (fstar.w[3] || fstar.w[2] 1539 || fstar.w[1] > ten2mk128[ind - 1].w[1] 1540 || (fstar.w[1] == ten2mk128[ind - 1].w[1] 1541 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) { 1542 if (++res.w[0] == 0) { 1543 res.w[1]++; 1544 } 1545 } 1546 } 1547 } 1548 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1]; 1549 BID_RETURN (res); 1550 } else { // if exp < 0 and q + exp <= 0 1551 if (x_sign) { // negative rounds up to -0.0 1552 res.w[1] = 0xb040000000000000ull; 1553 res.w[0] = 0x0000000000000000ull; 1554 } else { // positive rpunds up to +1.0 1555 res.w[1] = 0x3040000000000000ull; 1556 res.w[0] = 0x0000000000000001ull; 1557 } 1558 BID_RETURN (res); 1559 } 1560 } 1561 1562 /***************************************************************************** 1563 * BID128_round_integral_zero 1564 ****************************************************************************/ 1565 1566 BID128_FUNCTION_ARG1_NORND (bid128_round_integral_zero, x) 1567 1568 UINT128 res; 1569 UINT64 x_sign; 1570 UINT64 x_exp; 1571 int exp; // unbiased exponent 1572 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo 1573 // (all are UINT64) 1574 BID_UI64DOUBLE tmp1; 1575 unsigned int x_nr_bits; 1576 int q, ind, shift; 1577 UINT128 C1; 1578 // UINT128 res is C* at first - represents up to 34 decimal digits ~ 1579 // 113 bits 1580 UINT256 P256; 1581 1582 // check for NaN or Infinity 1583 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 1584 // x is special 1585 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN 1586 // if x = NaN, then res = Q (x) 1587 // check first for non-canonical NaN payload 1588 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) || 1589 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) && 1590 (x.w[0] > 0x38c15b09ffffffffull))) { 1591 x.w[1] = x.w[1] & 0xffffc00000000000ull; 1592 x.w[0] = 0x0ull; 1593 } 1594 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN 1595 // set invalid flag 1596 *pfpsf |= INVALID_EXCEPTION; 1597 // return quiet (x) 1598 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16] 1599 res.w[0] = x.w[0]; 1600 } else { // x is QNaN 1601 // return x 1602 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16] 1603 res.w[0] = x.w[0]; 1604 } 1605 BID_RETURN (res) 1606 } else { // x is not a NaN, so it must be infinity 1607 if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf 1608 // return +inf 1609 res.w[1] = 0x7800000000000000ull; 1610 res.w[0] = 0x0000000000000000ull; 1611 } else { // x is -inf 1612 // return -inf 1613 res.w[1] = 0xf800000000000000ull; 1614 res.w[0] = 0x0000000000000000ull; 1615 } 1616 BID_RETURN (res); 1617 } 1618 } 1619 // unpack x 1620 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 1621 C1.w[1] = x.w[1] & MASK_COEFF; 1622 C1.w[0] = x.w[0]; 1623 1624 // check for non-canonical values (treated as zero) 1625 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11 1626 // non-canonical 1627 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits 1628 C1.w[1] = 0; // significand high 1629 C1.w[0] = 0; // significand low 1630 } else { // G0_G1 != 11 1631 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits 1632 if (C1.w[1] > 0x0001ed09bead87c0ull || 1633 (C1.w[1] == 0x0001ed09bead87c0ull 1634 && C1.w[0] > 0x378d8e63ffffffffull)) { 1635 // x is non-canonical if coefficient is larger than 10^34 -1 1636 C1.w[1] = 0; 1637 C1.w[0] = 0; 1638 } else { // canonical 1639 ; 1640 } 1641 } 1642 1643 // test for input equal to zero 1644 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { 1645 // x is 0 1646 // return 0 preserving the sign bit and the preferred exponent 1647 // of MAX(Q(x), 0) 1648 if (x_exp <= (0x1820ull << 49)) { 1649 res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull; 1650 } else { 1651 res.w[1] = x_sign | x_exp; 1652 } 1653 res.w[0] = 0x0000000000000000ull; 1654 BID_RETURN (res); 1655 } 1656 // x is not special and is not zero 1657 1658 // if (exp <= -p) return -0.0 or +0.0 1659 if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34 1660 res.w[1] = x_sign | 0x3040000000000000ull; 1661 res.w[0] = 0x0000000000000000ull; 1662 BID_RETURN (res); 1663 } 1664 // q = nr. of decimal digits in x 1665 // determine first the nr. of bits in x 1666 if (C1.w[1] == 0) { 1667 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 1668 // split the 64-bit value in two 32-bit halves to avoid rounding errors 1669 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 1670 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion 1671 x_nr_bits = 1672 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1673 } else { // x < 2^32 1674 tmp1.d = (double) (C1.w[0]); // exact conversion 1675 x_nr_bits = 1676 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1677 } 1678 } else { // if x < 2^53 1679 tmp1.d = (double) C1.w[0]; // exact conversion 1680 x_nr_bits = 1681 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1682 } 1683 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) 1684 tmp1.d = (double) C1.w[1]; // exact conversion 1685 x_nr_bits = 1686 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1687 } 1688 1689 q = nr_digits[x_nr_bits - 1].digits; 1690 if (q == 0) { 1691 q = nr_digits[x_nr_bits - 1].digits1; 1692 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi || 1693 (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi && 1694 C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) 1695 q++; 1696 } 1697 exp = (x_exp >> 49) - 6176; 1698 if (exp >= 0) { // -exp <= 0 1699 // the argument is an integer already 1700 res.w[1] = x.w[1]; 1701 res.w[0] = x.w[0]; 1702 BID_RETURN (res); 1703 } else if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q 1704 // need to shift right -exp digits from the coefficient; the exp will be 0 1705 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x' 1706 // (number of digits to be chopped off) 1707 // chop off ind digits from the lower part of C1 1708 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate 1709 // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP 1710 // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE 1711 // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE 1712 //tmp64 = C1.w[0]; 1713 // if (ind <= 19) { 1714 // C1.w[0] = C1.w[0] + midpoint64[ind - 1]; 1715 // } else { 1716 // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; 1717 // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; 1718 // } 1719 // if (C1.w[0] < tmp64) C1.w[1]++; 1720 // if carry-out from C1.w[0], increment C1.w[1] 1721 // calculate C* and f* 1722 // C* is actually floor(C*) in this case 1723 // C* and f* need shifting and masking, as shown by 1724 // shiftright128[] and maskhigh128[] 1725 // 1 <= x <= 34 1726 // kx = 10^(-x) = ten2mk128[ind - 1] 1727 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 1728 // the approximation of 10^(-x) was rounded up to 118 bits 1729 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); 1730 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 1731 res.w[1] = P256.w[3]; 1732 res.w[0] = P256.w[2]; 1733 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 1734 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 1735 res.w[1] = (P256.w[3] >> shift); 1736 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift); 1737 } else { // 22 <= ind - 1 <= 33 1738 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38 1739 res.w[1] = 0; 1740 res.w[0] = P256.w[3] >> shift; 1741 } 1742 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1]; 1743 BID_RETURN (res); 1744 } else { // if exp < 0 and q + exp <= 0 the result is +0 or -0 1745 res.w[1] = x_sign | 0x3040000000000000ull; 1746 res.w[0] = 0x0000000000000000ull; 1747 BID_RETURN (res); 1748 } 1749 } 1750 1751 /***************************************************************************** 1752 * BID128_round_integral_nearest_away 1753 ****************************************************************************/ 1754 1755 BID128_FUNCTION_ARG1_NORND (bid128_round_integral_nearest_away, x) 1756 1757 UINT128 res; 1758 UINT64 x_sign; 1759 UINT64 x_exp; 1760 int exp; // unbiased exponent 1761 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo 1762 // (all are UINT64) 1763 UINT64 tmp64; 1764 BID_UI64DOUBLE tmp1; 1765 unsigned int x_nr_bits; 1766 int q, ind, shift; 1767 UINT128 C1; 1768 // UINT128 res is C* at first - represents up to 34 decimal digits ~ 1769 // 113 bits 1770 // UINT256 fstar; 1771 UINT256 P256; 1772 1773 // check for NaN or Infinity 1774 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 1775 // x is special 1776 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN 1777 // if x = NaN, then res = Q (x) 1778 // check first for non-canonical NaN payload 1779 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) || 1780 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) && 1781 (x.w[0] > 0x38c15b09ffffffffull))) { 1782 x.w[1] = x.w[1] & 0xffffc00000000000ull; 1783 x.w[0] = 0x0ull; 1784 } 1785 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN 1786 // set invalid flag 1787 *pfpsf |= INVALID_EXCEPTION; 1788 // return quiet (x) 1789 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16] 1790 res.w[0] = x.w[0]; 1791 } else { // x is QNaN 1792 // return x 1793 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16] 1794 res.w[0] = x.w[0]; 1795 } 1796 BID_RETURN (res) 1797 } else { // x is not a NaN, so it must be infinity 1798 if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf 1799 // return +inf 1800 res.w[1] = 0x7800000000000000ull; 1801 res.w[0] = 0x0000000000000000ull; 1802 } else { // x is -inf 1803 // return -inf 1804 res.w[1] = 0xf800000000000000ull; 1805 res.w[0] = 0x0000000000000000ull; 1806 } 1807 BID_RETURN (res); 1808 } 1809 } 1810 // unpack x 1811 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 1812 C1.w[1] = x.w[1] & MASK_COEFF; 1813 C1.w[0] = x.w[0]; 1814 1815 // check for non-canonical values (treated as zero) 1816 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11 1817 // non-canonical 1818 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits 1819 C1.w[1] = 0; // significand high 1820 C1.w[0] = 0; // significand low 1821 } else { // G0_G1 != 11 1822 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits 1823 if (C1.w[1] > 0x0001ed09bead87c0ull || 1824 (C1.w[1] == 0x0001ed09bead87c0ull 1825 && C1.w[0] > 0x378d8e63ffffffffull)) { 1826 // x is non-canonical if coefficient is larger than 10^34 -1 1827 C1.w[1] = 0; 1828 C1.w[0] = 0; 1829 } else { // canonical 1830 ; 1831 } 1832 } 1833 1834 // test for input equal to zero 1835 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { 1836 // x is 0 1837 // return 0 preserving the sign bit and the preferred exponent 1838 // of MAX(Q(x), 0) 1839 if (x_exp <= (0x1820ull << 49)) { 1840 res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull; 1841 } else { 1842 res.w[1] = x_sign | x_exp; 1843 } 1844 res.w[0] = 0x0000000000000000ull; 1845 BID_RETURN (res); 1846 } 1847 // x is not special and is not zero 1848 1849 // if (exp <= -(p+1)) return 0.0 1850 if (x_exp <= 0x2ffa000000000000ull) { // 0x2ffa000000000000ull == -35 1851 res.w[1] = x_sign | 0x3040000000000000ull; 1852 res.w[0] = 0x0000000000000000ull; 1853 BID_RETURN (res); 1854 } 1855 // q = nr. of decimal digits in x 1856 // determine first the nr. of bits in x 1857 if (C1.w[1] == 0) { 1858 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 1859 // split the 64-bit value in two 32-bit halves to avoid rounding errors 1860 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 1861 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion 1862 x_nr_bits = 1863 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1864 } else { // x < 2^32 1865 tmp1.d = (double) (C1.w[0]); // exact conversion 1866 x_nr_bits = 1867 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1868 } 1869 } else { // if x < 2^53 1870 tmp1.d = (double) C1.w[0]; // exact conversion 1871 x_nr_bits = 1872 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1873 } 1874 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) 1875 tmp1.d = (double) C1.w[1]; // exact conversion 1876 x_nr_bits = 1877 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 1878 } 1879 1880 q = nr_digits[x_nr_bits - 1].digits; 1881 if (q == 0) { 1882 q = nr_digits[x_nr_bits - 1].digits1; 1883 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi || 1884 (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi && 1885 C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) 1886 q++; 1887 } 1888 exp = (x_exp >> 49) - 6176; 1889 if (exp >= 0) { // -exp <= 0 1890 // the argument is an integer already 1891 res.w[1] = x.w[1]; 1892 res.w[0] = x.w[0]; 1893 BID_RETURN (res); 1894 } else if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q 1895 // need to shift right -exp digits from the coefficient; the exp will be 0 1896 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x' 1897 // chop off ind digits from the lower part of C1 1898 // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits 1899 tmp64 = C1.w[0]; 1900 if (ind <= 19) { 1901 C1.w[0] = C1.w[0] + midpoint64[ind - 1]; 1902 } else { 1903 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; 1904 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; 1905 } 1906 if (C1.w[0] < tmp64) 1907 C1.w[1]++; 1908 // calculate C* and f* 1909 // C* is actually floor(C*) in this case 1910 // C* and f* need shifting and masking, as shown by 1911 // shiftright128[] and maskhigh128[] 1912 // 1 <= x <= 34 1913 // kx = 10^(-x) = ten2mk128[ind - 1] 1914 // C* = (C1 + 1/2 * 10^x) * 10^(-x) 1915 // the approximation of 10^(-x) was rounded up to 118 bits 1916 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); 1917 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g. 1918 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999 1919 // if (0 < f* < 10^(-x)) then the result is a midpoint 1920 // if floor(C*) is even then C* = floor(C*) - logical right 1921 // shift; C* has p decimal digits, correct by Prop. 1) 1922 // else if floor(C*) is odd C* = floor(C*)-1 (logical right 1923 // shift; C* has p decimal digits, correct by Pr. 1) 1924 // else 1925 // C* = floor(C*) (logical right shift; C has p decimal digits, 1926 // correct by Property 1) 1927 // n = C* * 10^(e+x) 1928 1929 // shift right C* by Ex-128 = shiftright128[ind] 1930 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0 1931 res.w[1] = P256.w[3]; 1932 res.w[0] = P256.w[2]; 1933 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63 1934 shift = shiftright128[ind - 1]; // 3 <= shift <= 63 1935 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift); 1936 res.w[1] = (P256.w[3] >> shift); 1937 } else { // 22 <= ind - 1 <= 33 1938 shift = shiftright128[ind - 1]; // 2 <= shift <= 38 1939 res.w[1] = 0; 1940 res.w[0] = (P256.w[3] >> (shift - 64)); // 2 <= shift - 64 <= 38 1941 } 1942 // if the result was a midpoint, it was already rounded away from zero 1943 res.w[1] |= x_sign | 0x3040000000000000ull; 1944 BID_RETURN (res); 1945 } else { // if ((q + exp) < 0) <=> q < -exp 1946 // the result is +0 or -0 1947 res.w[1] = x_sign | 0x3040000000000000ull; 1948 res.w[0] = 0x0000000000000000ull; 1949 BID_RETURN (res); 1950 } 1951 } 1952