1 /* Copyright (C) 2007-2024 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 /***************************************************************************** 25 * 26 * BID128 fma x * y + z 27 * 28 ****************************************************************************/ 29 30 #include "bid_internal.h" 31 32 static void 33 rounding_correction (unsigned int rnd_mode, 34 unsigned int is_inexact_lt_midpoint, 35 unsigned int is_inexact_gt_midpoint, 36 unsigned int is_midpoint_lt_even, 37 unsigned int is_midpoint_gt_even, 38 int unbexp, 39 UINT128 * ptrres, _IDEC_flags * ptrfpsf) { 40 // unbiased true exponent unbexp may be larger than emax 41 42 UINT128 res = *ptrres; // expected to have the correct sign and coefficient 43 // (the exponent field is ignored, as unbexp is used instead) 44 UINT64 sign, exp; 45 UINT64 C_hi, C_lo; 46 47 // general correction from RN to RA, RM, RP, RZ 48 // Note: if the result is negative, then is_inexact_lt_midpoint, 49 // is_inexact_gt_midpoint, is_midpoint_lt_even, and is_midpoint_gt_even 50 // have to be considered as if determined for the absolute value of the 51 // result (so they seem to be reversed) 52 53 if (is_inexact_lt_midpoint || is_inexact_gt_midpoint || 54 is_midpoint_lt_even || is_midpoint_gt_even) { 55 *ptrfpsf |= INEXACT_EXCEPTION; 56 } 57 // apply correction to result calculated with unbounded exponent 58 sign = res.w[1] & MASK_SIGN; 59 exp = (UINT64) (unbexp + 6176) << 49; // valid only if expmin<=unbexp<=expmax 60 C_hi = res.w[1] & MASK_COEFF; 61 C_lo = res.w[0]; 62 if ((!sign && ((rnd_mode == ROUNDING_UP && is_inexact_lt_midpoint) || 63 ((rnd_mode == ROUNDING_TIES_AWAY || rnd_mode == ROUNDING_UP) && 64 is_midpoint_gt_even))) || 65 (sign && ((rnd_mode == ROUNDING_DOWN && is_inexact_lt_midpoint) || 66 ((rnd_mode == ROUNDING_TIES_AWAY || rnd_mode == ROUNDING_DOWN) && 67 is_midpoint_gt_even)))) { 68 // C = C + 1 69 C_lo = C_lo + 1; 70 if (C_lo == 0) 71 C_hi = C_hi + 1; 72 if (C_hi == 0x0001ed09bead87c0ull && C_lo == 0x378d8e6400000000ull) { 73 // C = 10^34 => rounding overflow 74 C_hi = 0x0000314dc6448d93ull; 75 C_lo = 0x38c15b0a00000000ull; // 10^33 76 // exp = exp + EXP_P1; 77 unbexp = unbexp + 1; 78 exp = (UINT64) (unbexp + 6176) << 49; 79 } 80 } else if ((is_midpoint_lt_even || is_inexact_gt_midpoint) && 81 ((sign && (rnd_mode == ROUNDING_UP || rnd_mode == ROUNDING_TO_ZERO)) || 82 (!sign && (rnd_mode == ROUNDING_DOWN || rnd_mode == ROUNDING_TO_ZERO)))) { 83 // C = C - 1 84 C_lo = C_lo - 1; 85 if (C_lo == 0xffffffffffffffffull) 86 C_hi--; 87 // check if we crossed into the lower decade 88 if (C_hi == 0x0000314dc6448d93ull && C_lo == 0x38c15b09ffffffffull) { 89 // C = 10^33 - 1 90 if (exp > 0) { 91 C_hi = 0x0001ed09bead87c0ull; // 10^34 - 1 92 C_lo = 0x378d8e63ffffffffull; 93 // exp = exp - EXP_P1; 94 unbexp = unbexp - 1; 95 exp = (UINT64) (unbexp + 6176) << 49; 96 } else { // if exp = 0 the result is tiny & inexact 97 *ptrfpsf |= UNDERFLOW_EXCEPTION; 98 } 99 } 100 } else { 101 ; // the result is already correct 102 } 103 if (unbexp > expmax) { // 6111 104 *ptrfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION); 105 exp = 0; 106 if (!sign) { // result is positive 107 if (rnd_mode == ROUNDING_UP || rnd_mode == ROUNDING_TIES_AWAY) { // +inf 108 C_hi = 0x7800000000000000ull; 109 C_lo = 0x0000000000000000ull; 110 } else { // res = +MAXFP = (10^34-1) * 10^emax 111 C_hi = 0x5fffed09bead87c0ull; 112 C_lo = 0x378d8e63ffffffffull; 113 } 114 } else { // result is negative 115 if (rnd_mode == ROUNDING_DOWN || rnd_mode == ROUNDING_TIES_AWAY) { // -inf 116 C_hi = 0xf800000000000000ull; 117 C_lo = 0x0000000000000000ull; 118 } else { // res = -MAXFP = -(10^34-1) * 10^emax 119 C_hi = 0xdfffed09bead87c0ull; 120 C_lo = 0x378d8e63ffffffffull; 121 } 122 } 123 } 124 // assemble the result 125 res.w[1] = sign | exp | C_hi; 126 res.w[0] = C_lo; 127 *ptrres = res; 128 } 129 130 static void 131 add256 (UINT256 x, UINT256 y, UINT256 * pz) { 132 // *z = x + yl assume the sum fits in 256 bits 133 UINT256 z; 134 z.w[0] = x.w[0] + y.w[0]; 135 if (z.w[0] < x.w[0]) { 136 x.w[1]++; 137 if (x.w[1] == 0x0000000000000000ull) { 138 x.w[2]++; 139 if (x.w[2] == 0x0000000000000000ull) { 140 x.w[3]++; 141 } 142 } 143 } 144 z.w[1] = x.w[1] + y.w[1]; 145 if (z.w[1] < x.w[1]) { 146 x.w[2]++; 147 if (x.w[2] == 0x0000000000000000ull) { 148 x.w[3]++; 149 } 150 } 151 z.w[2] = x.w[2] + y.w[2]; 152 if (z.w[2] < x.w[2]) { 153 x.w[3]++; 154 } 155 z.w[3] = x.w[3] + y.w[3]; // it was assumed that no carry is possible 156 *pz = z; 157 } 158 159 static void 160 sub256 (UINT256 x, UINT256 y, UINT256 * pz) { 161 // *z = x - y; assume x >= y 162 UINT256 z; 163 z.w[0] = x.w[0] - y.w[0]; 164 if (z.w[0] > x.w[0]) { 165 x.w[1]--; 166 if (x.w[1] == 0xffffffffffffffffull) { 167 x.w[2]--; 168 if (x.w[2] == 0xffffffffffffffffull) { 169 x.w[3]--; 170 } 171 } 172 } 173 z.w[1] = x.w[1] - y.w[1]; 174 if (z.w[1] > x.w[1]) { 175 x.w[2]--; 176 if (x.w[2] == 0xffffffffffffffffull) { 177 x.w[3]--; 178 } 179 } 180 z.w[2] = x.w[2] - y.w[2]; 181 if (z.w[2] > x.w[2]) { 182 x.w[3]--; 183 } 184 z.w[3] = x.w[3] - y.w[3]; // no borrow possible, because x >= y 185 *pz = z; 186 } 187 188 189 static int 190 nr_digits256 (UINT256 R256) { 191 int ind; 192 // determine the number of decimal digits in R256 193 if (R256.w[3] == 0x0 && R256.w[2] == 0x0 && R256.w[1] == 0x0) { 194 // between 1 and 19 digits 195 for (ind = 1; ind <= 19; ind++) { 196 if (R256.w[0] < ten2k64[ind]) { 197 break; 198 } 199 } 200 // ind digits 201 } else if (R256.w[3] == 0x0 && R256.w[2] == 0x0 && 202 (R256.w[1] < ten2k128[0].w[1] || 203 (R256.w[1] == ten2k128[0].w[1] 204 && R256.w[0] < ten2k128[0].w[0]))) { 205 // 20 digits 206 ind = 20; 207 } else if (R256.w[3] == 0x0 && R256.w[2] == 0x0) { 208 // between 21 and 38 digits 209 for (ind = 1; ind <= 18; ind++) { 210 if (R256.w[1] < ten2k128[ind].w[1] || 211 (R256.w[1] == ten2k128[ind].w[1] && 212 R256.w[0] < ten2k128[ind].w[0])) { 213 break; 214 } 215 } 216 // ind + 20 digits 217 ind = ind + 20; 218 } else if (R256.w[3] == 0x0 && 219 (R256.w[2] < ten2k256[0].w[2] || 220 (R256.w[2] == ten2k256[0].w[2] && 221 R256.w[1] < ten2k256[0].w[1]) || 222 (R256.w[2] == ten2k256[0].w[2] && 223 R256.w[1] == ten2k256[0].w[1] && 224 R256.w[0] < ten2k256[0].w[0]))) { 225 // 39 digits 226 ind = 39; 227 } else { 228 // between 40 and 68 digits 229 for (ind = 1; ind <= 29; ind++) { 230 if (R256.w[3] < ten2k256[ind].w[3] || 231 (R256.w[3] == ten2k256[ind].w[3] && 232 R256.w[2] < ten2k256[ind].w[2]) || 233 (R256.w[3] == ten2k256[ind].w[3] && 234 R256.w[2] == ten2k256[ind].w[2] && 235 R256.w[1] < ten2k256[ind].w[1]) || 236 (R256.w[3] == ten2k256[ind].w[3] && 237 R256.w[2] == ten2k256[ind].w[2] && 238 R256.w[1] == ten2k256[ind].w[1] && 239 R256.w[0] < ten2k256[ind].w[0])) { 240 break; 241 } 242 } 243 // ind + 39 digits 244 ind = ind + 39; 245 } 246 return (ind); 247 } 248 249 // add/subtract C4 and C3 * 10^scale; this may follow a previous rounding, so 250 // use the rounding information from ptr_is_* to avoid a double rounding error 251 static void 252 add_and_round (int q3, 253 int q4, 254 int e4, 255 int delta, 256 int p34, 257 UINT64 z_sign, 258 UINT64 p_sign, 259 UINT128 C3, 260 UINT256 C4, 261 int rnd_mode, 262 int *ptr_is_midpoint_lt_even, 263 int *ptr_is_midpoint_gt_even, 264 int *ptr_is_inexact_lt_midpoint, 265 int *ptr_is_inexact_gt_midpoint, 266 _IDEC_flags * ptrfpsf, UINT128 * ptrres) { 267 268 int scale; 269 int x0; 270 int ind; 271 UINT64 R64; 272 UINT128 P128, R128; 273 UINT192 P192, R192; 274 UINT256 R256; 275 int is_midpoint_lt_even = 0; 276 int is_midpoint_gt_even = 0; 277 int is_inexact_lt_midpoint = 0; 278 int is_inexact_gt_midpoint = 0; 279 int is_midpoint_lt_even0 = 0; 280 int is_midpoint_gt_even0 = 0; 281 int is_inexact_lt_midpoint0 = 0; 282 int is_inexact_gt_midpoint0 = 0; 283 int incr_exp = 0; 284 int is_tiny = 0; 285 int lt_half_ulp = 0; 286 int eq_half_ulp = 0; 287 // int gt_half_ulp = 0; 288 UINT128 res = *ptrres; 289 290 // scale C3 up by 10^(q4-delta-q3), 0 <= q4-delta-q3 <= 2*P34-2 = 66 291 scale = q4 - delta - q3; // 0 <= scale <= 66 (or 0 <= scale <= 68 if this 292 // comes from Cases (2), (3), (4), (5), (6), with 0 <= |delta| <= 1 293 294 // calculate C3 * 10^scale in R256 (it has at most 67 decimal digits for 295 // Cases (15),(16),(17) and at most 69 for Cases (2),(3),(4),(5),(6)) 296 if (scale == 0) { 297 R256.w[3] = 0x0ull; 298 R256.w[2] = 0x0ull; 299 R256.w[1] = C3.w[1]; 300 R256.w[0] = C3.w[0]; 301 } else if (scale <= 19) { // 10^scale fits in 64 bits 302 P128.w[1] = 0; 303 P128.w[0] = ten2k64[scale]; 304 __mul_128x128_to_256 (R256, P128, C3); 305 } else if (scale <= 38) { // 10^scale fits in 128 bits 306 __mul_128x128_to_256 (R256, ten2k128[scale - 20], C3); 307 } else if (scale <= 57) { // 39 <= scale <= 57 308 // 10^scale fits in 192 bits but C3 * 10^scale fits in 223 or 230 bits 309 // (10^67 has 223 bits; 10^69 has 230 bits); 310 // must split the computation: 311 // 10^scale * C3 = 10*38 * 10^(scale-38) * C3 where 10^38 takes 127 312 // bits and so 10^(scale-38) * C3 fits in 128 bits with certainty 313 // Note that 1 <= scale - 38 <= 19 => 10^(scale-38) fits in 64 bits 314 __mul_64x128_to_128 (R128, ten2k64[scale - 38], C3); 315 // now multiply R128 by 10^38 316 __mul_128x128_to_256 (R256, R128, ten2k128[18]); 317 } else { // 58 <= scale <= 66 318 // 10^scale takes between 193 and 220 bits, 319 // and C3 * 10^scale fits in 223 bits (10^67/10^69 has 223/230 bits) 320 // must split the computation: 321 // 10^scale * C3 = 10*38 * 10^(scale-38) * C3 where 10^38 takes 127 322 // bits and so 10^(scale-38) * C3 fits in 128 bits with certainty 323 // Note that 20 <= scale - 38 <= 30 => 10^(scale-38) fits in 128 bits 324 // Calculate first 10^(scale-38) * C3, which fits in 128 bits; because 325 // 10^(scale-38) takes more than 64 bits, C3 will take less than 64 326 __mul_64x128_to_128 (R128, C3.w[0], ten2k128[scale - 58]); 327 // now calculate 10*38 * 10^(scale-38) * C3 328 __mul_128x128_to_256 (R256, R128, ten2k128[18]); 329 } 330 // C3 * 10^scale is now in R256 331 332 // for Cases (15), (16), (17) C4 > C3 * 10^scale because C4 has at least 333 // one extra digit; for Cases (2), (3), (4), (5), or (6) any order is 334 // possible 335 // add/subtract C4 and C3 * 10^scale; the exponent is e4 336 if (p_sign == z_sign) { // R256 = C4 + R256 337 // calculate R256 = C4 + C3 * 10^scale = C4 + R256 which is exact, 338 // but may require rounding 339 add256 (C4, R256, &R256); 340 } else { // if (p_sign != z_sign) { // R256 = C4 - R256 341 // calculate R256 = C4 - C3 * 10^scale = C4 - R256 or 342 // R256 = C3 * 10^scale - C4 = R256 - C4 which is exact, 343 // but may require rounding 344 345 // compare first R256 = C3 * 10^scale and C4 346 if (R256.w[3] > C4.w[3] || (R256.w[3] == C4.w[3] && R256.w[2] > C4.w[2]) || 347 (R256.w[3] == C4.w[3] && R256.w[2] == C4.w[2] && R256.w[1] > C4.w[1]) || 348 (R256.w[3] == C4.w[3] && R256.w[2] == C4.w[2] && R256.w[1] == C4.w[1] && 349 R256.w[0] >= C4.w[0])) { // C3 * 10^scale >= C4 350 // calculate R256 = C3 * 10^scale - C4 = R256 - C4, which is exact, 351 // but may require rounding 352 sub256 (R256, C4, &R256); 353 // flip p_sign too, because the result has the sign of z 354 p_sign = z_sign; 355 } else { // if C4 > C3 * 10^scale 356 // calculate R256 = C4 - C3 * 10^scale = C4 - R256, which is exact, 357 // but may require rounding 358 sub256 (C4, R256, &R256); 359 } 360 // if the result is pure zero, the sign depends on the rounding mode 361 // (x*y and z had opposite signs) 362 if (R256.w[3] == 0x0ull && R256.w[2] == 0x0ull && 363 R256.w[1] == 0x0ull && R256.w[0] == 0x0ull) { 364 if (rnd_mode != ROUNDING_DOWN) 365 p_sign = 0x0000000000000000ull; 366 else 367 p_sign = 0x8000000000000000ull; 368 // the exponent is max (e4, expmin) 369 if (e4 < -6176) 370 e4 = expmin; 371 // assemble result 372 res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49); 373 res.w[0] = 0x0; 374 *ptrres = res; 375 return; 376 } 377 } 378 379 // determine the number of decimal digits in R256 380 ind = nr_digits256 (R256); 381 382 // the exact result is (-1)^p_sign * R256 * 10^e4 where q (R256) = ind; 383 // round to the destination precision, with unbounded exponent 384 385 if (ind <= p34) { 386 // result rounded to the destination precision with unbounded exponent 387 // is exact 388 if (ind + e4 < p34 + expmin) { 389 is_tiny = 1; // applies to all rounding modes 390 } 391 res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | R256.w[1]; 392 res.w[0] = R256.w[0]; 393 // Note: res is correct only if expmin <= e4 <= expmax 394 } else { // if (ind > p34) 395 // if more than P digits, round to nearest to P digits 396 // round R256 to p34 digits 397 x0 = ind - p34; // 1 <= x0 <= 34 as 35 <= ind <= 68 398 if (ind <= 38) { 399 P128.w[1] = R256.w[1]; 400 P128.w[0] = R256.w[0]; 401 round128_19_38 (ind, x0, P128, &R128, &incr_exp, 402 &is_midpoint_lt_even, &is_midpoint_gt_even, 403 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint); 404 } else if (ind <= 57) { 405 P192.w[2] = R256.w[2]; 406 P192.w[1] = R256.w[1]; 407 P192.w[0] = R256.w[0]; 408 round192_39_57 (ind, x0, P192, &R192, &incr_exp, 409 &is_midpoint_lt_even, &is_midpoint_gt_even, 410 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint); 411 R128.w[1] = R192.w[1]; 412 R128.w[0] = R192.w[0]; 413 } else { // if (ind <= 68) 414 round256_58_76 (ind, x0, R256, &R256, &incr_exp, 415 &is_midpoint_lt_even, &is_midpoint_gt_even, 416 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint); 417 R128.w[1] = R256.w[1]; 418 R128.w[0] = R256.w[0]; 419 } 420 // the rounded result has p34 = 34 digits 421 e4 = e4 + x0 + incr_exp; 422 if (rnd_mode == ROUNDING_TO_NEAREST) { 423 if (e4 < expmin) { 424 is_tiny = 1; // for other rounding modes apply correction 425 } 426 } else { 427 // for RM, RP, RZ, RA apply correction in order to determine tininess 428 // but do not save the result; apply the correction to 429 // (-1)^p_sign * significand * 10^0 430 P128.w[1] = p_sign | 0x3040000000000000ull | R128.w[1]; 431 P128.w[0] = R128.w[0]; 432 rounding_correction (rnd_mode, 433 is_inexact_lt_midpoint, 434 is_inexact_gt_midpoint, is_midpoint_lt_even, 435 is_midpoint_gt_even, 0, &P128, ptrfpsf); 436 scale = ((P128.w[1] & MASK_EXP) >> 49) - 6176; // -1, 0, or +1 437 // the number of digits in the significand is p34 = 34 438 if (e4 + scale < expmin) { 439 is_tiny = 1; 440 } 441 } 442 ind = p34; // the number of decimal digits in the signifcand of res 443 res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | R128.w[1]; // RN 444 res.w[0] = R128.w[0]; 445 // Note: res is correct only if expmin <= e4 <= expmax 446 // set the inexact flag after rounding with bounded exponent, if any 447 } 448 // at this point we have the result rounded with unbounded exponent in 449 // res and we know its tininess: 450 // res = (-1)^p_sign * significand * 10^e4, 451 // where q (significand) = ind <= p34 452 // Note: res is correct only if expmin <= e4 <= expmax 453 454 // check for overflow if RN 455 if (rnd_mode == ROUNDING_TO_NEAREST && (ind + e4) > (p34 + expmax)) { 456 res.w[1] = p_sign | 0x7800000000000000ull; 457 res.w[0] = 0x0000000000000000ull; 458 *ptrres = res; 459 *ptrfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION); 460 return; // BID_RETURN (res) 461 } // else not overflow or not RN, so continue 462 463 // if (e4 >= expmin) we have the result rounded with bounded exponent 464 if (e4 < expmin) { 465 x0 = expmin - e4; // x0 >= 1; the number of digits to chop off of res 466 // where the result rounded [at most] once is 467 // (-1)^p_sign * significand_res * 10^e4 468 469 // avoid double rounding error 470 is_inexact_lt_midpoint0 = is_inexact_lt_midpoint; 471 is_inexact_gt_midpoint0 = is_inexact_gt_midpoint; 472 is_midpoint_lt_even0 = is_midpoint_lt_even; 473 is_midpoint_gt_even0 = is_midpoint_gt_even; 474 is_inexact_lt_midpoint = 0; 475 is_inexact_gt_midpoint = 0; 476 is_midpoint_lt_even = 0; 477 is_midpoint_gt_even = 0; 478 479 if (x0 > ind) { 480 // nothing is left of res when moving the decimal point left x0 digits 481 is_inexact_lt_midpoint = 1; 482 res.w[1] = p_sign | 0x0000000000000000ull; 483 res.w[0] = 0x0000000000000000ull; 484 e4 = expmin; 485 } else if (x0 == ind) { // 1 <= x0 = ind <= p34 = 34 486 // this is <, =, or > 1/2 ulp 487 // compare the ind-digit value in the significand of res with 488 // 1/2 ulp = 5*10^(ind-1), i.e. determine whether it is 489 // less than, equal to, or greater than 1/2 ulp (significand of res) 490 R128.w[1] = res.w[1] & MASK_COEFF; 491 R128.w[0] = res.w[0]; 492 if (ind <= 19) { 493 if (R128.w[0] < midpoint64[ind - 1]) { // < 1/2 ulp 494 lt_half_ulp = 1; 495 is_inexact_lt_midpoint = 1; 496 } else if (R128.w[0] == midpoint64[ind - 1]) { // = 1/2 ulp 497 eq_half_ulp = 1; 498 is_midpoint_gt_even = 1; 499 } else { // > 1/2 ulp 500 // gt_half_ulp = 1; 501 is_inexact_gt_midpoint = 1; 502 } 503 } else { // if (ind <= 38) { 504 if (R128.w[1] < midpoint128[ind - 20].w[1] || 505 (R128.w[1] == midpoint128[ind - 20].w[1] && 506 R128.w[0] < midpoint128[ind - 20].w[0])) { // < 1/2 ulp 507 lt_half_ulp = 1; 508 is_inexact_lt_midpoint = 1; 509 } else if (R128.w[1] == midpoint128[ind - 20].w[1] && 510 R128.w[0] == midpoint128[ind - 20].w[0]) { // = 1/2 ulp 511 eq_half_ulp = 1; 512 is_midpoint_gt_even = 1; 513 } else { // > 1/2 ulp 514 // gt_half_ulp = 1; 515 is_inexact_gt_midpoint = 1; 516 } 517 } 518 if (lt_half_ulp || eq_half_ulp) { 519 // res = +0.0 * 10^expmin 520 res.w[1] = 0x0000000000000000ull; 521 res.w[0] = 0x0000000000000000ull; 522 } else { // if (gt_half_ulp) 523 // res = +1 * 10^expmin 524 res.w[1] = 0x0000000000000000ull; 525 res.w[0] = 0x0000000000000001ull; 526 } 527 res.w[1] = p_sign | res.w[1]; 528 e4 = expmin; 529 } else { // if (1 <= x0 <= ind - 1 <= 33) 530 // round the ind-digit result to ind - x0 digits 531 532 if (ind <= 18) { // 2 <= ind <= 18 533 round64_2_18 (ind, x0, res.w[0], &R64, &incr_exp, 534 &is_midpoint_lt_even, &is_midpoint_gt_even, 535 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint); 536 res.w[1] = 0x0; 537 res.w[0] = R64; 538 } else if (ind <= 38) { 539 P128.w[1] = res.w[1] & MASK_COEFF; 540 P128.w[0] = res.w[0]; 541 round128_19_38 (ind, x0, P128, &res, &incr_exp, 542 &is_midpoint_lt_even, &is_midpoint_gt_even, 543 &is_inexact_lt_midpoint, 544 &is_inexact_gt_midpoint); 545 } 546 e4 = e4 + x0; // expmin 547 // we want the exponent to be expmin, so if incr_exp = 1 then 548 // multiply the rounded result by 10 - it will still fit in 113 bits 549 if (incr_exp) { 550 // 64 x 128 -> 128 551 P128.w[1] = res.w[1] & MASK_COEFF; 552 P128.w[0] = res.w[0]; 553 __mul_64x128_to_128 (res, ten2k64[1], P128); 554 } 555 res.w[1] = 556 p_sign | ((UINT64) (e4 + 6176) << 49) | (res.w[1] & MASK_COEFF); 557 // avoid a double rounding error 558 if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 559 is_midpoint_lt_even) { // double rounding error upward 560 // res = res - 1 561 res.w[0]--; 562 if (res.w[0] == 0xffffffffffffffffull) 563 res.w[1]--; 564 // Note: a double rounding error upward is not possible; for this 565 // the result after the first rounding would have to be 99...95 566 // (35 digits in all), possibly followed by a number of zeros; this 567 // is not possible in Cases (2)-(6) or (15)-(17) which may get here 568 is_midpoint_lt_even = 0; 569 is_inexact_lt_midpoint = 1; 570 } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 571 is_midpoint_gt_even) { // double rounding error downward 572 // res = res + 1 573 res.w[0]++; 574 if (res.w[0] == 0) 575 res.w[1]++; 576 is_midpoint_gt_even = 0; 577 is_inexact_gt_midpoint = 1; 578 } else if (!is_midpoint_lt_even && !is_midpoint_gt_even && 579 !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) { 580 // if this second rounding was exact the result may still be 581 // inexact because of the first rounding 582 if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) { 583 is_inexact_gt_midpoint = 1; 584 } 585 if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) { 586 is_inexact_lt_midpoint = 1; 587 } 588 } else if (is_midpoint_gt_even && 589 (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) { 590 // pulled up to a midpoint 591 is_inexact_lt_midpoint = 1; 592 is_inexact_gt_midpoint = 0; 593 is_midpoint_lt_even = 0; 594 is_midpoint_gt_even = 0; 595 } else if (is_midpoint_lt_even && 596 (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) { 597 // pulled down to a midpoint 598 is_inexact_lt_midpoint = 0; 599 is_inexact_gt_midpoint = 1; 600 is_midpoint_lt_even = 0; 601 is_midpoint_gt_even = 0; 602 } else { 603 ; 604 } 605 } 606 } 607 // res contains the correct result 608 // apply correction if not rounding to nearest 609 if (rnd_mode != ROUNDING_TO_NEAREST) { 610 rounding_correction (rnd_mode, 611 is_inexact_lt_midpoint, is_inexact_gt_midpoint, 612 is_midpoint_lt_even, is_midpoint_gt_even, 613 e4, &res, ptrfpsf); 614 } 615 if (is_midpoint_lt_even || is_midpoint_gt_even || 616 is_inexact_lt_midpoint || is_inexact_gt_midpoint) { 617 // set the inexact flag 618 *ptrfpsf |= INEXACT_EXCEPTION; 619 if (is_tiny) 620 *ptrfpsf |= UNDERFLOW_EXCEPTION; 621 } 622 623 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 624 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 625 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 626 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 627 *ptrres = res; 628 return; 629 } 630 631 632 #if DECIMAL_CALL_BY_REFERENCE 633 static void 634 bid128_ext_fma (int *ptr_is_midpoint_lt_even, 635 int *ptr_is_midpoint_gt_even, 636 int *ptr_is_inexact_lt_midpoint, 637 int *ptr_is_inexact_gt_midpoint, UINT128 * pres, 638 UINT128 * px, UINT128 * py, 639 UINT128 * 640 pz _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 641 _EXC_INFO_PARAM) { 642 UINT128 x = *px, y = *py, z = *pz; 643 #if !DECIMAL_GLOBAL_ROUNDING 644 unsigned int rnd_mode = *prnd_mode; 645 #endif 646 #else 647 static UINT128 648 bid128_ext_fma (int *ptr_is_midpoint_lt_even, 649 int *ptr_is_midpoint_gt_even, 650 int *ptr_is_inexact_lt_midpoint, 651 int *ptr_is_inexact_gt_midpoint, UINT128 x, UINT128 y, 652 UINT128 z _RND_MODE_PARAM _EXC_FLAGS_PARAM 653 _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 654 #endif 655 656 UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} }; 657 UINT64 x_sign, y_sign, z_sign, p_sign, tmp_sign; 658 UINT64 x_exp = 0, y_exp = 0, z_exp = 0, p_exp; 659 int true_p_exp; 660 UINT128 C1, C2, C3; 661 UINT256 C4; 662 int q1 = 0, q2 = 0, q3 = 0, q4; 663 int e1, e2, e3, e4; 664 int scale, ind, delta, x0; 665 int p34 = P34; // used to modify the limit on the number of digits 666 BID_UI64DOUBLE tmp; 667 int x_nr_bits, y_nr_bits, z_nr_bits; 668 unsigned int save_fpsf; 669 int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0; 670 int is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0; 671 int is_midpoint_lt_even0 = 0, is_midpoint_gt_even0 = 0; 672 int is_inexact_lt_midpoint0 = 0, is_inexact_gt_midpoint0 = 0; 673 int incr_exp = 0; 674 int lsb; 675 int lt_half_ulp = 0; 676 int eq_half_ulp = 0; 677 int gt_half_ulp = 0; 678 int is_tiny = 0; 679 UINT64 R64, tmp64; 680 UINT128 P128, R128; 681 UINT192 P192, R192; 682 UINT256 R256; 683 684 // the following are based on the table of special cases for fma; the NaN 685 // behavior is similar to that of the IA-64 Architecture fma 686 687 // identify cases where at least one operand is NaN 688 689 BID_SWAP128 (x); 690 BID_SWAP128 (y); 691 BID_SWAP128 (z); 692 if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN 693 // if x = {0, f, inf, NaN}, y = NaN, z = {0, f, inf, NaN} then res = Q (y) 694 // check first for non-canonical NaN payload 695 if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) || 696 (((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) && 697 (y.w[0] > 0x38c15b09ffffffffull))) { 698 y.w[1] = y.w[1] & 0xffffc00000000000ull; 699 y.w[0] = 0x0ull; 700 } 701 if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN 702 // set invalid flag 703 *pfpsf |= INVALID_EXCEPTION; 704 // return quiet (y) 705 res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16] 706 res.w[0] = y.w[0]; 707 } else { // y is QNaN 708 // return y 709 res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16] 710 res.w[0] = y.w[0]; 711 // if z = SNaN or x = SNaN signal invalid exception 712 if ((z.w[1] & MASK_SNAN) == MASK_SNAN || 713 (x.w[1] & MASK_SNAN) == MASK_SNAN) { 714 // set invalid flag 715 *pfpsf |= INVALID_EXCEPTION; 716 } 717 } 718 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 719 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 720 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 721 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 722 BID_SWAP128 (res); 723 BID_RETURN (res) 724 } else if ((z.w[1] & MASK_NAN) == MASK_NAN) { // z is NAN 725 // if x = {0, f, inf, NaN}, y = {0, f, inf}, z = NaN then res = Q (z) 726 // check first for non-canonical NaN payload 727 if (((z.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) || 728 (((z.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) && 729 (z.w[0] > 0x38c15b09ffffffffull))) { 730 z.w[1] = z.w[1] & 0xffffc00000000000ull; 731 z.w[0] = 0x0ull; 732 } 733 if ((z.w[1] & MASK_SNAN) == MASK_SNAN) { // z is SNAN 734 // set invalid flag 735 *pfpsf |= INVALID_EXCEPTION; 736 // return quiet (z) 737 res.w[1] = z.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16] 738 res.w[0] = z.w[0]; 739 } else { // z is QNaN 740 // return z 741 res.w[1] = z.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16] 742 res.w[0] = z.w[0]; 743 // if x = SNaN signal invalid exception 744 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { 745 // set invalid flag 746 *pfpsf |= INVALID_EXCEPTION; 747 } 748 } 749 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 750 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 751 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 752 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 753 BID_SWAP128 (res); 754 BID_RETURN (res) 755 } else if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN 756 // if x = NaN, y = {0, f, inf}, z = {0, f, inf} then res = Q (x) 757 // check first for non-canonical NaN payload 758 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) || 759 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) && 760 (x.w[0] > 0x38c15b09ffffffffull))) { 761 x.w[1] = x.w[1] & 0xffffc00000000000ull; 762 x.w[0] = 0x0ull; 763 } 764 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN 765 // set invalid flag 766 *pfpsf |= INVALID_EXCEPTION; 767 // return quiet (x) 768 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16] 769 res.w[0] = x.w[0]; 770 } else { // x is QNaN 771 // return x 772 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16] 773 res.w[0] = x.w[0]; 774 } 775 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 776 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 777 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 778 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 779 BID_SWAP128 (res); 780 BID_RETURN (res) 781 } 782 // x, y, z are 0, f, or inf but not NaN => unpack the arguments and check 783 // for non-canonical values 784 785 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 786 C1.w[1] = x.w[1] & MASK_COEFF; 787 C1.w[0] = x.w[0]; 788 if ((x.w[1] & MASK_ANY_INF) != MASK_INF) { // x != inf 789 // if x is not infinity check for non-canonical values - treated as zero 790 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11 791 // non-canonical 792 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits 793 C1.w[1] = 0; // significand high 794 C1.w[0] = 0; // significand low 795 } else { // G0_G1 != 11 796 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits 797 if (C1.w[1] > 0x0001ed09bead87c0ull || 798 (C1.w[1] == 0x0001ed09bead87c0ull && 799 C1.w[0] > 0x378d8e63ffffffffull)) { 800 // x is non-canonical if coefficient is larger than 10^34 -1 801 C1.w[1] = 0; 802 C1.w[0] = 0; 803 } else { // canonical 804 ; 805 } 806 } 807 } 808 y_sign = y.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 809 C2.w[1] = y.w[1] & MASK_COEFF; 810 C2.w[0] = y.w[0]; 811 if ((y.w[1] & MASK_ANY_INF) != MASK_INF) { // y != inf 812 // if y is not infinity check for non-canonical values - treated as zero 813 if ((y.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11 814 // non-canonical 815 y_exp = (y.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits 816 C2.w[1] = 0; // significand high 817 C2.w[0] = 0; // significand low 818 } else { // G0_G1 != 11 819 y_exp = y.w[1] & MASK_EXP; // biased and shifted left 49 bits 820 if (C2.w[1] > 0x0001ed09bead87c0ull || 821 (C2.w[1] == 0x0001ed09bead87c0ull && 822 C2.w[0] > 0x378d8e63ffffffffull)) { 823 // y is non-canonical if coefficient is larger than 10^34 -1 824 C2.w[1] = 0; 825 C2.w[0] = 0; 826 } else { // canonical 827 ; 828 } 829 } 830 } 831 z_sign = z.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 832 C3.w[1] = z.w[1] & MASK_COEFF; 833 C3.w[0] = z.w[0]; 834 if ((z.w[1] & MASK_ANY_INF) != MASK_INF) { // z != inf 835 // if z is not infinity check for non-canonical values - treated as zero 836 if ((z.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11 837 // non-canonical 838 z_exp = (z.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits 839 C3.w[1] = 0; // significand high 840 C3.w[0] = 0; // significand low 841 } else { // G0_G1 != 11 842 z_exp = z.w[1] & MASK_EXP; // biased and shifted left 49 bits 843 if (C3.w[1] > 0x0001ed09bead87c0ull || 844 (C3.w[1] == 0x0001ed09bead87c0ull && 845 C3.w[0] > 0x378d8e63ffffffffull)) { 846 // z is non-canonical if coefficient is larger than 10^34 -1 847 C3.w[1] = 0; 848 C3.w[0] = 0; 849 } else { // canonical 850 ; 851 } 852 } 853 } 854 855 p_sign = x_sign ^ y_sign; // sign of the product 856 857 // identify cases where at least one operand is infinity 858 859 if ((x.w[1] & MASK_ANY_INF) == MASK_INF) { // x = inf 860 if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf 861 if ((z.w[1] & MASK_ANY_INF) == MASK_INF) { // z = inf 862 if (p_sign == z_sign) { 863 res.w[1] = z_sign | MASK_INF; 864 res.w[0] = 0x0; 865 } else { 866 // return QNaN Indefinite 867 res.w[1] = 0x7c00000000000000ull; 868 res.w[0] = 0x0000000000000000ull; 869 // set invalid flag 870 *pfpsf |= INVALID_EXCEPTION; 871 } 872 } else { // z = 0 or z = f 873 res.w[1] = p_sign | MASK_INF; 874 res.w[0] = 0x0; 875 } 876 } else if (C2.w[1] != 0 || C2.w[0] != 0) { // y = f 877 if ((z.w[1] & MASK_ANY_INF) == MASK_INF) { // z = inf 878 if (p_sign == z_sign) { 879 res.w[1] = z_sign | MASK_INF; 880 res.w[0] = 0x0; 881 } else { 882 // return QNaN Indefinite 883 res.w[1] = 0x7c00000000000000ull; 884 res.w[0] = 0x0000000000000000ull; 885 // set invalid flag 886 *pfpsf |= INVALID_EXCEPTION; 887 } 888 } else { // z = 0 or z = f 889 res.w[1] = p_sign | MASK_INF; 890 res.w[0] = 0x0; 891 } 892 } else { // y = 0 893 // return QNaN Indefinite 894 res.w[1] = 0x7c00000000000000ull; 895 res.w[0] = 0x0000000000000000ull; 896 // set invalid flag 897 *pfpsf |= INVALID_EXCEPTION; 898 } 899 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 900 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 901 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 902 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 903 BID_SWAP128 (res); 904 BID_RETURN (res) 905 } else if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf 906 if ((z.w[1] & MASK_ANY_INF) == MASK_INF) { // z = inf 907 // x = f, necessarily 908 if ((p_sign != z_sign) 909 || (C1.w[1] == 0x0ull && C1.w[0] == 0x0ull)) { 910 // return QNaN Indefinite 911 res.w[1] = 0x7c00000000000000ull; 912 res.w[0] = 0x0000000000000000ull; 913 // set invalid flag 914 *pfpsf |= INVALID_EXCEPTION; 915 } else { 916 res.w[1] = z_sign | MASK_INF; 917 res.w[0] = 0x0; 918 } 919 } else if (C1.w[1] == 0x0 && C1.w[0] == 0x0) { // x = 0 920 // z = 0, f, inf 921 // return QNaN Indefinite 922 res.w[1] = 0x7c00000000000000ull; 923 res.w[0] = 0x0000000000000000ull; 924 // set invalid flag 925 *pfpsf |= INVALID_EXCEPTION; 926 } else { 927 // x = f and z = 0, f, necessarily 928 res.w[1] = p_sign | MASK_INF; 929 res.w[0] = 0x0; 930 } 931 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 932 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 933 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 934 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 935 BID_SWAP128 (res); 936 BID_RETURN (res) 937 } else if ((z.w[1] & MASK_ANY_INF) == MASK_INF) { // z = inf 938 // x = 0, f and y = 0, f, necessarily 939 res.w[1] = z_sign | MASK_INF; 940 res.w[0] = 0x0; 941 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 942 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 943 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 944 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 945 BID_SWAP128 (res); 946 BID_RETURN (res) 947 } 948 949 true_p_exp = (x_exp >> 49) - 6176 + (y_exp >> 49) - 6176; 950 if (true_p_exp < -6176) 951 p_exp = 0; // cannot be less than EXP_MIN 952 else 953 p_exp = (UINT64) (true_p_exp + 6176) << 49; 954 955 if (((C1.w[1] == 0x0 && C1.w[0] == 0x0) || (C2.w[1] == 0x0 && C2.w[0] == 0x0)) && C3.w[1] == 0x0 && C3.w[0] == 0x0) { // (x = 0 or y = 0) and z = 0 956 // the result is 0 957 if (p_exp < z_exp) 958 res.w[1] = p_exp; // preferred exponent 959 else 960 res.w[1] = z_exp; // preferred exponent 961 if (p_sign == z_sign) { 962 res.w[1] |= z_sign; 963 res.w[0] = 0x0; 964 } else { // x * y and z have opposite signs 965 if (rnd_mode == ROUNDING_DOWN) { 966 // res = -0.0 967 res.w[1] |= MASK_SIGN; 968 res.w[0] = 0x0; 969 } else { 970 // res = +0.0 971 // res.w[1] |= 0x0; 972 res.w[0] = 0x0; 973 } 974 } 975 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 976 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 977 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 978 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 979 BID_SWAP128 (res); 980 BID_RETURN (res) 981 } 982 // from this point on, we may need to know the number of decimal digits 983 // in the significands of x, y, z when x, y, z != 0 984 985 if (C1.w[1] != 0 || C1.w[0] != 0) { // x = f (non-zero finite) 986 // q1 = nr. of decimal digits in x 987 // determine first the nr. of bits in x 988 if (C1.w[1] == 0) { 989 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 990 // split the 64-bit value in two 32-bit halves to avoid rounding errors 991 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 992 tmp.d = (double) (C1.w[0] >> 32); // exact conversion 993 x_nr_bits = 994 33 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 995 } else { // x < 2^32 996 tmp.d = (double) (C1.w[0]); // exact conversion 997 x_nr_bits = 998 1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 999 } 1000 } else { // if x < 2^53 1001 tmp.d = (double) C1.w[0]; // exact conversion 1002 x_nr_bits = 1003 1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 1004 } 1005 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) 1006 tmp.d = (double) C1.w[1]; // exact conversion 1007 x_nr_bits = 1008 65 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 1009 } 1010 q1 = nr_digits[x_nr_bits - 1].digits; 1011 if (q1 == 0) { 1012 q1 = nr_digits[x_nr_bits - 1].digits1; 1013 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi || 1014 (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi && 1015 C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) 1016 q1++; 1017 } 1018 } 1019 1020 if (C2.w[1] != 0 || C2.w[0] != 0) { // y = f (non-zero finite) 1021 if (C2.w[1] == 0) { 1022 if (C2.w[0] >= 0x0020000000000000ull) { // y >= 2^53 1023 // split the 64-bit value in two 32-bit halves to avoid rounding errors 1024 if (C2.w[0] >= 0x0000000100000000ull) { // y >= 2^32 1025 tmp.d = (double) (C2.w[0] >> 32); // exact conversion 1026 y_nr_bits = 1027 32 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 1028 } else { // y < 2^32 1029 tmp.d = (double) C2.w[0]; // exact conversion 1030 y_nr_bits = 1031 ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 1032 } 1033 } else { // if y < 2^53 1034 tmp.d = (double) C2.w[0]; // exact conversion 1035 y_nr_bits = 1036 ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 1037 } 1038 } else { // C2.w[1] != 0 => nr. bits = 64 + nr_bits (C2.w[1]) 1039 tmp.d = (double) C2.w[1]; // exact conversion 1040 y_nr_bits = 1041 64 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 1042 } 1043 1044 q2 = nr_digits[y_nr_bits].digits; 1045 if (q2 == 0) { 1046 q2 = nr_digits[y_nr_bits].digits1; 1047 if (C2.w[1] > nr_digits[y_nr_bits].threshold_hi || 1048 (C2.w[1] == nr_digits[y_nr_bits].threshold_hi && 1049 C2.w[0] >= nr_digits[y_nr_bits].threshold_lo)) 1050 q2++; 1051 } 1052 } 1053 1054 if (C3.w[1] != 0 || C3.w[0] != 0) { // z = f (non-zero finite) 1055 if (C3.w[1] == 0) { 1056 if (C3.w[0] >= 0x0020000000000000ull) { // z >= 2^53 1057 // split the 64-bit value in two 32-bit halves to avoid rounding errors 1058 if (C3.w[0] >= 0x0000000100000000ull) { // z >= 2^32 1059 tmp.d = (double) (C3.w[0] >> 32); // exact conversion 1060 z_nr_bits = 1061 32 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 1062 } else { // z < 2^32 1063 tmp.d = (double) C3.w[0]; // exact conversion 1064 z_nr_bits = 1065 ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 1066 } 1067 } else { // if z < 2^53 1068 tmp.d = (double) C3.w[0]; // exact conversion 1069 z_nr_bits = 1070 ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 1071 } 1072 } else { // C3.w[1] != 0 => nr. bits = 64 + nr_bits (C3.w[1]) 1073 tmp.d = (double) C3.w[1]; // exact conversion 1074 z_nr_bits = 1075 64 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 1076 } 1077 1078 q3 = nr_digits[z_nr_bits].digits; 1079 if (q3 == 0) { 1080 q3 = nr_digits[z_nr_bits].digits1; 1081 if (C3.w[1] > nr_digits[z_nr_bits].threshold_hi || 1082 (C3.w[1] == nr_digits[z_nr_bits].threshold_hi && 1083 C3.w[0] >= nr_digits[z_nr_bits].threshold_lo)) 1084 q3++; 1085 } 1086 } 1087 1088 if ((C1.w[1] == 0x0 && C1.w[0] == 0x0) || 1089 (C2.w[1] == 0x0 && C2.w[0] == 0x0)) { 1090 // x = 0 or y = 0 1091 // z = f, necessarily; for 0 + z return z, with the preferred exponent 1092 // the result is z, but need to get the preferred exponent 1093 if (z_exp <= p_exp) { // the preferred exponent is z_exp 1094 res.w[1] = z_sign | (z_exp & MASK_EXP) | C3.w[1]; 1095 res.w[0] = C3.w[0]; 1096 } else { // if (p_exp < z_exp) the preferred exponent is p_exp 1097 // return (C3 * 10^scale) * 10^(z_exp - scale) 1098 // where scale = min (p34-q3, (z_exp-p_exp) >> 49) 1099 scale = p34 - q3; 1100 ind = (z_exp - p_exp) >> 49; 1101 if (ind < scale) 1102 scale = ind; 1103 if (scale == 0) { 1104 res.w[1] = z.w[1]; // & MASK_COEFF, which is redundant 1105 res.w[0] = z.w[0]; 1106 } else if (q3 <= 19) { // z fits in 64 bits 1107 if (scale <= 19) { // 10^scale fits in 64 bits 1108 // 64 x 64 C3.w[0] * ten2k64[scale] 1109 __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]); 1110 } else { // 10^scale fits in 128 bits 1111 // 64 x 128 C3.w[0] * ten2k128[scale - 20] 1112 __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]); 1113 } 1114 } else { // z fits in 128 bits, but 10^scale must fit in 64 bits 1115 // 64 x 128 ten2k64[scale] * C3 1116 __mul_128x64_to_128 (res, ten2k64[scale], C3); 1117 } 1118 // subtract scale from the exponent 1119 z_exp = z_exp - ((UINT64) scale << 49); 1120 res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1]; 1121 } 1122 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 1123 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 1124 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 1125 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 1126 BID_SWAP128 (res); 1127 BID_RETURN (res) 1128 } else { 1129 ; // continue with x = f, y = f, z = 0 or x = f, y = f, z = f 1130 } 1131 1132 e1 = (x_exp >> 49) - 6176; // unbiased exponent of x 1133 e2 = (y_exp >> 49) - 6176; // unbiased exponent of y 1134 e3 = (z_exp >> 49) - 6176; // unbiased exponent of z 1135 e4 = e1 + e2; // unbiased exponent of the exact x * y 1136 1137 // calculate C1 * C2 and its number of decimal digits, q4 1138 1139 // the exact product has either q1 + q2 - 1 or q1 + q2 decimal digits 1140 // where 2 <= q1 + q2 <= 68 1141 // calculate C4 = C1 * C2 and determine q 1142 C4.w[3] = C4.w[2] = C4.w[1] = C4.w[0] = 0; 1143 if (q1 + q2 <= 19) { // if 2 <= q1 + q2 <= 19, C4 = C1 * C2 fits in 64 bits 1144 C4.w[0] = C1.w[0] * C2.w[0]; 1145 // if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2 1146 if (C4.w[0] < ten2k64[q1 + q2 - 1]) 1147 q4 = q1 + q2 - 1; // q4 in [1, 18] 1148 else 1149 q4 = q1 + q2; // q4 in [2, 19] 1150 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 64; 1151 } else if (q1 + q2 == 20) { // C4 = C1 * C2 fits in 64 or 128 bits 1152 // q1 <= 19 and q2 <= 19 so both C1 and C2 fit in 64 bits 1153 __mul_64x64_to_128MACH (C4, C1.w[0], C2.w[0]); 1154 // if C4 < 10^(q1+q2-1) = 10^19 then q4 = q1+q2-1 = 19 else q4 = q1+q2 = 20 1155 if (C4.w[1] == 0 && C4.w[0] < ten2k64[19]) { // 19 = q1+q2-1 1156 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 64; 1157 q4 = 19; // 19 = q1 + q2 - 1 1158 } else { 1159 // if (C4.w[1] == 0) 1160 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 64; 1161 // else 1162 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128; 1163 q4 = 20; // 20 = q1 + q2 1164 } 1165 } else if (q1 + q2 <= 38) { // 21 <= q1 + q2 <= 38 1166 // C4 = C1 * C2 fits in 64 or 128 bits 1167 // (64 bits possibly, but only when q1 + q2 = 21 and C4 has 20 digits) 1168 // at least one of C1, C2 has at most 19 decimal digits & fits in 64 bits 1169 if (q1 <= 19) { 1170 __mul_128x64_to_128 (C4, C1.w[0], C2); 1171 } else { // q2 <= 19 1172 __mul_128x64_to_128 (C4, C2.w[0], C1); 1173 } 1174 // if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2 1175 if (C4.w[1] < ten2k128[q1 + q2 - 21].w[1] || 1176 (C4.w[1] == ten2k128[q1 + q2 - 21].w[1] && 1177 C4.w[0] < ten2k128[q1 + q2 - 21].w[0])) { 1178 // if (C4.w[1] == 0) // q4 = 20, necessarily 1179 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 64; 1180 // else 1181 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128; 1182 q4 = q1 + q2 - 1; // q4 in [20, 37] 1183 } else { 1184 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128; 1185 q4 = q1 + q2; // q4 in [21, 38] 1186 } 1187 } else if (q1 + q2 == 39) { // C4 = C1 * C2 fits in 128 or 192 bits 1188 // both C1 and C2 fit in 128 bits (actually in 113 bits) 1189 // may replace this by 128x128_to192 1190 __mul_128x128_to_256 (C4, C1, C2); // C4.w[3] is 0 1191 // if C4 < 10^(q1+q2-1) = 10^38 then q4 = q1+q2-1 = 38 else q4 = q1+q2 = 39 1192 if (C4.w[2] == 0 && (C4.w[1] < ten2k128[18].w[1] || 1193 (C4.w[1] == ten2k128[18].w[1] 1194 && C4.w[0] < ten2k128[18].w[0]))) { 1195 // 18 = 38 - 20 = q1+q2-1 - 20 1196 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128; 1197 q4 = 38; // 38 = q1 + q2 - 1 1198 } else { 1199 // if (C4.w[2] == 0) 1200 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128; 1201 // else 1202 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192; 1203 q4 = 39; // 39 = q1 + q2 1204 } 1205 } else if (q1 + q2 <= 57) { // 40 <= q1 + q2 <= 57 1206 // C4 = C1 * C2 fits in 128 or 192 bits 1207 // (128 bits possibly, but only when q1 + q2 = 40 and C4 has 39 digits) 1208 // both C1 and C2 fit in 128 bits (actually in 113 bits); at most one 1209 // may fit in 64 bits 1210 if (C1.w[1] == 0) { // C1 fits in 64 bits 1211 // __mul_64x128_full (REShi64, RESlo128, A64, B128) 1212 __mul_64x128_full (C4.w[2], C4, C1.w[0], C2); 1213 } else if (C2.w[1] == 0) { // C2 fits in 64 bits 1214 // __mul_64x128_full (REShi64, RESlo128, A64, B128) 1215 __mul_64x128_full (C4.w[2], C4, C2.w[0], C1); 1216 } else { // both C1 and C2 require 128 bits 1217 // may use __mul_128x128_to_192 (C4.w[2], C4.w[0], C2.w[0], C1); 1218 __mul_128x128_to_256 (C4, C1, C2); // C4.w[3] = 0 1219 } 1220 // if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2 1221 if (C4.w[2] < ten2k256[q1 + q2 - 40].w[2] || 1222 (C4.w[2] == ten2k256[q1 + q2 - 40].w[2] && 1223 (C4.w[1] < ten2k256[q1 + q2 - 40].w[1] || 1224 (C4.w[1] == ten2k256[q1 + q2 - 40].w[1] && 1225 C4.w[0] < ten2k256[q1 + q2 - 40].w[0])))) { 1226 // if (C4.w[2] == 0) // q4 = 39, necessarily 1227 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128; 1228 // else 1229 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192; 1230 q4 = q1 + q2 - 1; // q4 in [39, 56] 1231 } else { 1232 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192; 1233 q4 = q1 + q2; // q4 in [40, 57] 1234 } 1235 } else if (q1 + q2 == 58) { // C4 = C1 * C2 fits in 192 or 256 bits 1236 // both C1 and C2 fit in 128 bits (actually in 113 bits); at most one 1237 // may fit in 64 bits 1238 if (C1.w[1] == 0) { // C1 * C2 will fit in 192 bits 1239 __mul_64x128_full (C4.w[2], C4, C1.w[0], C2); // may use 64x128_to_192 1240 } else if (C2.w[1] == 0) { // C1 * C2 will fit in 192 bits 1241 __mul_64x128_full (C4.w[2], C4, C2.w[0], C1); // may use 64x128_to_192 1242 } else { // C1 * C2 will fit in 192 bits or in 256 bits 1243 __mul_128x128_to_256 (C4, C1, C2); 1244 } 1245 // if C4 < 10^(q1+q2-1) = 10^57 then q4 = q1+q2-1 = 57 else q4 = q1+q2 = 58 1246 if (C4.w[3] == 0 && (C4.w[2] < ten2k256[18].w[2] || 1247 (C4.w[2] == ten2k256[18].w[2] 1248 && (C4.w[1] < ten2k256[18].w[1] 1249 || (C4.w[1] == ten2k256[18].w[1] 1250 && C4.w[0] < ten2k256[18].w[0]))))) { 1251 // 18 = 57 - 39 = q1+q2-1 - 39 1252 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192; 1253 q4 = 57; // 57 = q1 + q2 - 1 1254 } else { 1255 // if (C4.w[3] == 0) 1256 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192; 1257 // else 1258 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 256; 1259 q4 = 58; // 58 = q1 + q2 1260 } 1261 } else { // if 59 <= q1 + q2 <= 68 1262 // C4 = C1 * C2 fits in 192 or 256 bits 1263 // (192 bits possibly, but only when q1 + q2 = 59 and C4 has 58 digits) 1264 // both C1 and C2 fit in 128 bits (actually in 113 bits); none fits in 1265 // 64 bits 1266 // may use __mul_128x128_to_192 (C4.w[2], C4.w[0], C2.w[0], C1); 1267 __mul_128x128_to_256 (C4, C1, C2); // C4.w[3] = 0 1268 // if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2 1269 if (C4.w[3] < ten2k256[q1 + q2 - 40].w[3] || 1270 (C4.w[3] == ten2k256[q1 + q2 - 40].w[3] && 1271 (C4.w[2] < ten2k256[q1 + q2 - 40].w[2] || 1272 (C4.w[2] == ten2k256[q1 + q2 - 40].w[2] && 1273 (C4.w[1] < ten2k256[q1 + q2 - 40].w[1] || 1274 (C4.w[1] == ten2k256[q1 + q2 - 40].w[1] && 1275 C4.w[0] < ten2k256[q1 + q2 - 40].w[0])))))) { 1276 // if (C4.w[3] == 0) // q4 = 58, necessarily 1277 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192; 1278 // else 1279 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 256; 1280 q4 = q1 + q2 - 1; // q4 in [58, 67] 1281 } else { 1282 // length of C1 * C2 rounded up to a multiple of 64 bits is len = 256; 1283 q4 = q1 + q2; // q4 in [59, 68] 1284 } 1285 } 1286 1287 if (C3.w[1] == 0x0 && C3.w[0] == 0x0) { // x = f, y = f, z = 0 1288 save_fpsf = *pfpsf; // sticky bits - caller value must be preserved 1289 *pfpsf = 0; 1290 1291 if (q4 > p34) { 1292 1293 // truncate C4 to p34 digits into res 1294 // x = q4-p34, 1 <= x <= 34 because 35 <= q4 <= 68 1295 x0 = q4 - p34; 1296 if (q4 <= 38) { 1297 P128.w[1] = C4.w[1]; 1298 P128.w[0] = C4.w[0]; 1299 round128_19_38 (q4, x0, P128, &res, &incr_exp, 1300 &is_midpoint_lt_even, &is_midpoint_gt_even, 1301 &is_inexact_lt_midpoint, 1302 &is_inexact_gt_midpoint); 1303 } else if (q4 <= 57) { // 35 <= q4 <= 57 1304 P192.w[2] = C4.w[2]; 1305 P192.w[1] = C4.w[1]; 1306 P192.w[0] = C4.w[0]; 1307 round192_39_57 (q4, x0, P192, &R192, &incr_exp, 1308 &is_midpoint_lt_even, &is_midpoint_gt_even, 1309 &is_inexact_lt_midpoint, 1310 &is_inexact_gt_midpoint); 1311 res.w[0] = R192.w[0]; 1312 res.w[1] = R192.w[1]; 1313 } else { // if (q4 <= 68) 1314 round256_58_76 (q4, x0, C4, &R256, &incr_exp, 1315 &is_midpoint_lt_even, &is_midpoint_gt_even, 1316 &is_inexact_lt_midpoint, 1317 &is_inexact_gt_midpoint); 1318 res.w[0] = R256.w[0]; 1319 res.w[1] = R256.w[1]; 1320 } 1321 e4 = e4 + x0; 1322 if (incr_exp) { 1323 e4 = e4 + 1; 1324 } 1325 q4 = p34; 1326 // res is now the coefficient of the result rounded to the destination 1327 // precision, with unbounded exponent; the exponent is e4; q4=digits(res) 1328 } else { // if (q4 <= p34) 1329 // C4 * 10^e4 is the result rounded to the destination precision, with 1330 // unbounded exponent (which is exact) 1331 1332 if ((q4 + e4 <= p34 + expmax) && (e4 > expmax)) { 1333 // e4 is too large, but can be brought within range by scaling up C4 1334 scale = e4 - expmax; // 1 <= scale < P-q4 <= P-1 => 1 <= scale <= P-2 1335 // res = (C4 * 10^scale) * 10^expmax 1336 if (q4 <= 19) { // C4 fits in 64 bits 1337 if (scale <= 19) { // 10^scale fits in 64 bits 1338 // 64 x 64 C4.w[0] * ten2k64[scale] 1339 __mul_64x64_to_128MACH (res, C4.w[0], ten2k64[scale]); 1340 } else { // 10^scale fits in 128 bits 1341 // 64 x 128 C4.w[0] * ten2k128[scale - 20] 1342 __mul_128x64_to_128 (res, C4.w[0], ten2k128[scale - 20]); 1343 } 1344 } else { // C4 fits in 128 bits, but 10^scale must fit in 64 bits 1345 // 64 x 128 ten2k64[scale] * CC43 1346 __mul_128x64_to_128 (res, ten2k64[scale], C4); 1347 } 1348 e4 = e4 - scale; // expmax 1349 q4 = q4 + scale; 1350 } else { 1351 res.w[1] = C4.w[1]; 1352 res.w[0] = C4.w[0]; 1353 } 1354 // res is the coefficient of the result rounded to the destination 1355 // precision, with unbounded exponent (it has q4 digits); the exponent 1356 // is e4 (exact result) 1357 } 1358 1359 // check for overflow 1360 if (q4 + e4 > p34 + expmax) { 1361 if (rnd_mode == ROUNDING_TO_NEAREST) { 1362 res.w[1] = p_sign | 0x7800000000000000ull; // +/-inf 1363 res.w[0] = 0x0000000000000000ull; 1364 *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION); 1365 } else { 1366 res.w[1] = p_sign | res.w[1]; 1367 rounding_correction (rnd_mode, 1368 is_inexact_lt_midpoint, 1369 is_inexact_gt_midpoint, 1370 is_midpoint_lt_even, is_midpoint_gt_even, 1371 e4, &res, pfpsf); 1372 } 1373 *pfpsf |= save_fpsf; 1374 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 1375 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 1376 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 1377 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 1378 BID_SWAP128 (res); 1379 BID_RETURN (res) 1380 } 1381 // check for underflow 1382 if (q4 + e4 < expmin + P34) { 1383 is_tiny = 1; // the result is tiny 1384 if (e4 < expmin) { 1385 // if e4 < expmin, we must truncate more of res 1386 x0 = expmin - e4; // x0 >= 1 1387 is_inexact_lt_midpoint0 = is_inexact_lt_midpoint; 1388 is_inexact_gt_midpoint0 = is_inexact_gt_midpoint; 1389 is_midpoint_lt_even0 = is_midpoint_lt_even; 1390 is_midpoint_gt_even0 = is_midpoint_gt_even; 1391 is_inexact_lt_midpoint = 0; 1392 is_inexact_gt_midpoint = 0; 1393 is_midpoint_lt_even = 0; 1394 is_midpoint_gt_even = 0; 1395 // the number of decimal digits in res is q4 1396 if (x0 < q4) { // 1 <= x0 <= q4-1 => round res to q4 - x0 digits 1397 if (q4 <= 18) { // 2 <= q4 <= 18, 1 <= x0 <= 17 1398 round64_2_18 (q4, x0, res.w[0], &R64, &incr_exp, 1399 &is_midpoint_lt_even, &is_midpoint_gt_even, 1400 &is_inexact_lt_midpoint, 1401 &is_inexact_gt_midpoint); 1402 if (incr_exp) { 1403 // R64 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 17 1404 R64 = ten2k64[q4 - x0]; 1405 } 1406 // res.w[1] = 0; (from above) 1407 res.w[0] = R64; 1408 } else { // if (q4 <= 34) 1409 // 19 <= q4 <= 38 1410 P128.w[1] = res.w[1]; 1411 P128.w[0] = res.w[0]; 1412 round128_19_38 (q4, x0, P128, &res, &incr_exp, 1413 &is_midpoint_lt_even, &is_midpoint_gt_even, 1414 &is_inexact_lt_midpoint, 1415 &is_inexact_gt_midpoint); 1416 if (incr_exp) { 1417 // increase coefficient by a factor of 10; this will be <= 10^33 1418 // R128 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 37 1419 if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19 1420 // res.w[1] = 0; 1421 res.w[0] = ten2k64[q4 - x0]; 1422 } else { // 20 <= q4 - x0 <= 37 1423 res.w[0] = ten2k128[q4 - x0 - 20].w[0]; 1424 res.w[1] = ten2k128[q4 - x0 - 20].w[1]; 1425 } 1426 } 1427 } 1428 e4 = e4 + x0; // expmin 1429 } else if (x0 == q4) { 1430 // the second rounding is for 0.d(0)d(1)...d(q4-1) * 10^emin 1431 // determine relationship with 1/2 ulp 1432 if (q4 <= 19) { 1433 if (res.w[0] < midpoint64[q4 - 1]) { // < 1/2 ulp 1434 lt_half_ulp = 1; 1435 is_inexact_lt_midpoint = 1; 1436 } else if (res.w[0] == midpoint64[q4 - 1]) { // = 1/2 ulp 1437 eq_half_ulp = 1; 1438 is_midpoint_gt_even = 1; 1439 } else { // > 1/2 ulp 1440 // gt_half_ulp = 1; 1441 is_inexact_gt_midpoint = 1; 1442 } 1443 } else { // if (q4 <= 34) 1444 if (res.w[1] < midpoint128[q4 - 20].w[1] || 1445 (res.w[1] == midpoint128[q4 - 20].w[1] && 1446 res.w[0] < midpoint128[q4 - 20].w[0])) { // < 1/2 ulp 1447 lt_half_ulp = 1; 1448 is_inexact_lt_midpoint = 1; 1449 } else if (res.w[1] == midpoint128[q4 - 20].w[1] && 1450 res.w[0] == midpoint128[q4 - 20].w[0]) { // = 1/2 ulp 1451 eq_half_ulp = 1; 1452 is_midpoint_gt_even = 1; 1453 } else { // > 1/2 ulp 1454 // gt_half_ulp = 1; 1455 is_inexact_gt_midpoint = 1; 1456 } 1457 } 1458 if (lt_half_ulp || eq_half_ulp) { 1459 // res = +0.0 * 10^expmin 1460 res.w[1] = 0x0000000000000000ull; 1461 res.w[0] = 0x0000000000000000ull; 1462 } else { // if (gt_half_ulp) 1463 // res = +1 * 10^expmin 1464 res.w[1] = 0x0000000000000000ull; 1465 res.w[0] = 0x0000000000000001ull; 1466 } 1467 e4 = expmin; 1468 } else { // if (x0 > q4) 1469 // the second rounding is for 0.0...d(0)d(1)...d(q4-1) * 10^emin 1470 res.w[1] = 0; 1471 res.w[0] = 0; 1472 e4 = expmin; 1473 is_inexact_lt_midpoint = 1; 1474 } 1475 // avoid a double rounding error 1476 if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 1477 is_midpoint_lt_even) { // double rounding error upward 1478 // res = res - 1 1479 res.w[0]--; 1480 if (res.w[0] == 0xffffffffffffffffull) 1481 res.w[1]--; 1482 // Note: a double rounding error upward is not possible; for this 1483 // the result after the first rounding would have to be 99...95 1484 // (35 digits in all), possibly followed by a number of zeros; this 1485 // not possible for f * f + 0 1486 is_midpoint_lt_even = 0; 1487 is_inexact_lt_midpoint = 1; 1488 } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 1489 is_midpoint_gt_even) { // double rounding error downward 1490 // res = res + 1 1491 res.w[0]++; 1492 if (res.w[0] == 0) 1493 res.w[1]++; 1494 is_midpoint_gt_even = 0; 1495 is_inexact_gt_midpoint = 1; 1496 } else if (!is_midpoint_lt_even && !is_midpoint_gt_even && 1497 !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) { 1498 // if this second rounding was exact the result may still be 1499 // inexact because of the first rounding 1500 if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) { 1501 is_inexact_gt_midpoint = 1; 1502 } 1503 if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) { 1504 is_inexact_lt_midpoint = 1; 1505 } 1506 } else if (is_midpoint_gt_even && 1507 (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) { 1508 // pulled up to a midpoint 1509 is_inexact_lt_midpoint = 1; 1510 is_inexact_gt_midpoint = 0; 1511 is_midpoint_lt_even = 0; 1512 is_midpoint_gt_even = 0; 1513 } else if (is_midpoint_lt_even && 1514 (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) { 1515 // pulled down to a midpoint 1516 is_inexact_lt_midpoint = 0; 1517 is_inexact_gt_midpoint = 1; 1518 is_midpoint_lt_even = 0; 1519 is_midpoint_gt_even = 0; 1520 } else { 1521 ; 1522 } 1523 } else { // if e4 >= emin then q4 < P and the result is tiny and exact 1524 if (e3 < e4) { 1525 // if (e3 < e4) the preferred exponent is e3 1526 // return (C4 * 10^scale) * 10^(e4 - scale) 1527 // where scale = min (p34-q4, (e4 - e3)) 1528 scale = p34 - q4; 1529 ind = e4 - e3; 1530 if (ind < scale) 1531 scale = ind; 1532 if (scale == 0) { 1533 ; // res and e4 are unchanged 1534 } else if (q4 <= 19) { // C4 fits in 64 bits 1535 if (scale <= 19) { // 10^scale fits in 64 bits 1536 // 64 x 64 res.w[0] * ten2k64[scale] 1537 __mul_64x64_to_128MACH (res, res.w[0], ten2k64[scale]); 1538 } else { // 10^scale fits in 128 bits 1539 // 64 x 128 res.w[0] * ten2k128[scale - 20] 1540 __mul_128x64_to_128 (res, res.w[0], ten2k128[scale - 20]); 1541 } 1542 } else { // res fits in 128 bits, but 10^scale must fit in 64 bits 1543 // 64 x 128 ten2k64[scale] * C3 1544 __mul_128x64_to_128 (res, ten2k64[scale], res); 1545 } 1546 // subtract scale from the exponent 1547 e4 = e4 - scale; 1548 } 1549 } 1550 1551 // check for inexact result 1552 if (is_inexact_lt_midpoint || is_inexact_gt_midpoint || 1553 is_midpoint_lt_even || is_midpoint_gt_even) { 1554 // set the inexact flag and the underflow flag 1555 *pfpsf |= INEXACT_EXCEPTION; 1556 *pfpsf |= UNDERFLOW_EXCEPTION; 1557 } 1558 res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | res.w[1]; 1559 if (rnd_mode != ROUNDING_TO_NEAREST) { 1560 rounding_correction (rnd_mode, 1561 is_inexact_lt_midpoint, 1562 is_inexact_gt_midpoint, 1563 is_midpoint_lt_even, is_midpoint_gt_even, 1564 e4, &res, pfpsf); 1565 } 1566 *pfpsf |= save_fpsf; 1567 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 1568 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 1569 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 1570 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 1571 BID_SWAP128 (res); 1572 BID_RETURN (res) 1573 } 1574 // no overflow, and no underflow for rounding to nearest 1575 res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | res.w[1]; 1576 1577 if (rnd_mode != ROUNDING_TO_NEAREST) { 1578 rounding_correction (rnd_mode, 1579 is_inexact_lt_midpoint, 1580 is_inexact_gt_midpoint, 1581 is_midpoint_lt_even, is_midpoint_gt_even, 1582 e4, &res, pfpsf); 1583 // if e4 = expmin && significand < 10^33 => result is tiny (for RD, RZ) 1584 if (e4 == expmin) { 1585 if ((res.w[1] & MASK_COEFF) < 0x0000314dc6448d93ull || 1586 ((res.w[1] & MASK_COEFF) == 0x0000314dc6448d93ull && 1587 res.w[0] < 0x38c15b0a00000000ull)) { 1588 is_tiny = 1; 1589 } 1590 } 1591 } 1592 1593 if (is_inexact_lt_midpoint || is_inexact_gt_midpoint || 1594 is_midpoint_lt_even || is_midpoint_gt_even) { 1595 // set the inexact flag 1596 *pfpsf |= INEXACT_EXCEPTION; 1597 if (is_tiny) 1598 *pfpsf |= UNDERFLOW_EXCEPTION; 1599 } 1600 1601 if ((*pfpsf & INEXACT_EXCEPTION) == 0) { // x * y is exact 1602 // need to ensure that the result has the preferred exponent 1603 p_exp = res.w[1] & MASK_EXP; 1604 if (z_exp < p_exp) { // the preferred exponent is z_exp 1605 // signficand of res in C3 1606 C3.w[1] = res.w[1] & MASK_COEFF; 1607 C3.w[0] = res.w[0]; 1608 // the number of decimal digits of x * y is q4 <= 34 1609 // Note: the coefficient fits in 128 bits 1610 1611 // return (C3 * 10^scale) * 10^(p_exp - scale) 1612 // where scale = min (p34-q4, (p_exp-z_exp) >> 49) 1613 scale = p34 - q4; 1614 ind = (p_exp - z_exp) >> 49; 1615 if (ind < scale) 1616 scale = ind; 1617 // subtract scale from the exponent 1618 p_exp = p_exp - ((UINT64) scale << 49); 1619 if (scale == 0) { 1620 ; // leave res unchanged 1621 } else if (q4 <= 19) { // x * y fits in 64 bits 1622 if (scale <= 19) { // 10^scale fits in 64 bits 1623 // 64 x 64 C3.w[0] * ten2k64[scale] 1624 __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]); 1625 } else { // 10^scale fits in 128 bits 1626 // 64 x 128 C3.w[0] * ten2k128[scale - 20] 1627 __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]); 1628 } 1629 res.w[1] = p_sign | (p_exp & MASK_EXP) | res.w[1]; 1630 } else { // x * y fits in 128 bits, but 10^scale must fit in 64 bits 1631 // 64 x 128 ten2k64[scale] * C3 1632 __mul_128x64_to_128 (res, ten2k64[scale], C3); 1633 res.w[1] = p_sign | (p_exp & MASK_EXP) | res.w[1]; 1634 } 1635 } // else leave the result as it is, because p_exp <= z_exp 1636 } 1637 *pfpsf |= save_fpsf; 1638 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 1639 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 1640 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 1641 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 1642 BID_SWAP128 (res); 1643 BID_RETURN (res) 1644 } // else we have f * f + f 1645 1646 // continue with x = f, y = f, z = f 1647 1648 delta = q3 + e3 - q4 - e4; 1649 delta_ge_zero: 1650 if (delta >= 0) { 1651 1652 if (p34 <= delta - 1 || // Case (1') 1653 (p34 == delta && e3 + 6176 < p34 - q3)) { // Case (1''A) 1654 // check for overflow, which can occur only in Case (1') 1655 if ((q3 + e3) > (p34 + expmax) && p34 <= delta - 1) { 1656 // e3 > expmax implies p34 <= delta-1 and e3 > expmax is a necessary 1657 // condition for (q3 + e3) > (p34 + expmax) 1658 if (rnd_mode == ROUNDING_TO_NEAREST) { 1659 res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf 1660 res.w[0] = 0x0000000000000000ull; 1661 *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION); 1662 } else { 1663 if (p_sign == z_sign) { 1664 is_inexact_lt_midpoint = 1; 1665 } else { 1666 is_inexact_gt_midpoint = 1; 1667 } 1668 // q3 <= p34; if (q3 < p34) scale C3 up by 10^(p34-q3) 1669 scale = p34 - q3; 1670 if (scale == 0) { 1671 res.w[1] = z_sign | C3.w[1]; 1672 res.w[0] = C3.w[0]; 1673 } else { 1674 if (q3 <= 19) { // C3 fits in 64 bits 1675 if (scale <= 19) { // 10^scale fits in 64 bits 1676 // 64 x 64 C3.w[0] * ten2k64[scale] 1677 __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]); 1678 } else { // 10^scale fits in 128 bits 1679 // 64 x 128 C3.w[0] * ten2k128[scale - 20] 1680 __mul_128x64_to_128 (res, C3.w[0], 1681 ten2k128[scale - 20]); 1682 } 1683 } else { // C3 fits in 128 bits, but 10^scale must fit in 64 bits 1684 // 64 x 128 ten2k64[scale] * C3 1685 __mul_128x64_to_128 (res, ten2k64[scale], C3); 1686 } 1687 // the coefficient in res has q3 + scale = p34 digits 1688 } 1689 e3 = e3 - scale; 1690 res.w[1] = z_sign | res.w[1]; 1691 rounding_correction (rnd_mode, 1692 is_inexact_lt_midpoint, 1693 is_inexact_gt_midpoint, 1694 is_midpoint_lt_even, is_midpoint_gt_even, 1695 e3, &res, pfpsf); 1696 } 1697 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 1698 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 1699 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 1700 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 1701 BID_SWAP128 (res); 1702 BID_RETURN (res) 1703 } 1704 // res = z 1705 if (q3 < p34) { // the preferred exponent is z_exp - (p34 - q3) 1706 // return (C3 * 10^scale) * 10^(z_exp - scale) 1707 // where scale = min (p34-q3, z_exp-EMIN) 1708 scale = p34 - q3; 1709 ind = e3 + 6176; 1710 if (ind < scale) 1711 scale = ind; 1712 if (scale == 0) { 1713 res.w[1] = C3.w[1]; 1714 res.w[0] = C3.w[0]; 1715 } else if (q3 <= 19) { // z fits in 64 bits 1716 if (scale <= 19) { // 10^scale fits in 64 bits 1717 // 64 x 64 C3.w[0] * ten2k64[scale] 1718 __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]); 1719 } else { // 10^scale fits in 128 bits 1720 // 64 x 128 C3.w[0] * ten2k128[scale - 20] 1721 __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]); 1722 } 1723 } else { // z fits in 128 bits, but 10^scale must fit in 64 bits 1724 // 64 x 128 ten2k64[scale] * C3 1725 __mul_128x64_to_128 (res, ten2k64[scale], C3); 1726 } 1727 // the coefficient in res has q3 + scale digits 1728 // subtract scale from the exponent 1729 z_exp = z_exp - ((UINT64) scale << 49); 1730 e3 = e3 - scale; 1731 res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1]; 1732 if (scale + q3 < p34) 1733 *pfpsf |= UNDERFLOW_EXCEPTION; 1734 } else { 1735 scale = 0; 1736 res.w[1] = z_sign | ((UINT64) (e3 + 6176) << 49) | C3.w[1]; 1737 res.w[0] = C3.w[0]; 1738 } 1739 1740 // use the following to avoid double rounding errors when operating on 1741 // mixed formats in rounding to nearest, and for correcting the result 1742 // if not rounding to nearest 1743 if ((p_sign != z_sign) && (delta == (q3 + scale + 1))) { 1744 // there is a gap of exactly one digit between the scaled C3 and C4 1745 // C3 * 10^ scale = 10^(q3+scale-1) <=> C3 = 10^(q3-1) is special case 1746 if ((q3 <= 19 && C3.w[0] != ten2k64[q3 - 1]) || 1747 (q3 == 20 && (C3.w[1] != 0 || C3.w[0] != ten2k64[19])) || 1748 (q3 >= 21 && (C3.w[1] != ten2k128[q3 - 21].w[1] || 1749 C3.w[0] != ten2k128[q3 - 21].w[0]))) { 1750 // C3 * 10^ scale != 10^(q3-1) 1751 // if ((res.w[1] & MASK_COEFF) != 0x0000314dc6448d93ull || 1752 // res.w[0] != 0x38c15b0a00000000ull) { // C3 * 10^scale != 10^33 1753 is_inexact_gt_midpoint = 1; // if (z_sign), set as if for abs. value 1754 } else { // if C3 * 10^scale = 10^(q3+scale-1) 1755 // ok from above e3 = (z_exp >> 49) - 6176; 1756 // the result is always inexact 1757 if (q4 == 1) { 1758 R64 = C4.w[0]; 1759 } else { 1760 // if q4 > 1 then truncate C4 from q4 digits to 1 digit; 1761 // x = q4-1, 1 <= x <= 67 and check if this operation is exact 1762 if (q4 <= 18) { // 2 <= q4 <= 18 1763 round64_2_18 (q4, q4 - 1, C4.w[0], &R64, &incr_exp, 1764 &is_midpoint_lt_even, &is_midpoint_gt_even, 1765 &is_inexact_lt_midpoint, 1766 &is_inexact_gt_midpoint); 1767 } else if (q4 <= 38) { 1768 P128.w[1] = C4.w[1]; 1769 P128.w[0] = C4.w[0]; 1770 round128_19_38 (q4, q4 - 1, P128, &R128, &incr_exp, 1771 &is_midpoint_lt_even, 1772 &is_midpoint_gt_even, 1773 &is_inexact_lt_midpoint, 1774 &is_inexact_gt_midpoint); 1775 R64 = R128.w[0]; // one decimal digit 1776 } else if (q4 <= 57) { 1777 P192.w[2] = C4.w[2]; 1778 P192.w[1] = C4.w[1]; 1779 P192.w[0] = C4.w[0]; 1780 round192_39_57 (q4, q4 - 1, P192, &R192, &incr_exp, 1781 &is_midpoint_lt_even, 1782 &is_midpoint_gt_even, 1783 &is_inexact_lt_midpoint, 1784 &is_inexact_gt_midpoint); 1785 R64 = R192.w[0]; // one decimal digit 1786 } else { // if (q4 <= 68) 1787 round256_58_76 (q4, q4 - 1, C4, &R256, &incr_exp, 1788 &is_midpoint_lt_even, 1789 &is_midpoint_gt_even, 1790 &is_inexact_lt_midpoint, 1791 &is_inexact_gt_midpoint); 1792 R64 = R256.w[0]; // one decimal digit 1793 } 1794 if (incr_exp) { 1795 R64 = 10; 1796 } 1797 } 1798 if (q4 == 1 && C4.w[0] == 5) { 1799 is_inexact_lt_midpoint = 0; 1800 is_inexact_gt_midpoint = 0; 1801 is_midpoint_lt_even = 1; 1802 is_midpoint_gt_even = 0; 1803 } else if ((e3 == expmin) || 1804 R64 < 5 || (R64 == 5 && is_inexact_gt_midpoint)) { 1805 // result does not change 1806 is_inexact_lt_midpoint = 0; 1807 is_inexact_gt_midpoint = 1; 1808 is_midpoint_lt_even = 0; 1809 is_midpoint_gt_even = 0; 1810 } else { 1811 is_inexact_lt_midpoint = 1; 1812 is_inexact_gt_midpoint = 0; 1813 is_midpoint_lt_even = 0; 1814 is_midpoint_gt_even = 0; 1815 // result decremented is 10^(q3+scale) - 1 1816 if ((q3 + scale) <= 19) { 1817 res.w[1] = 0; 1818 res.w[0] = ten2k64[q3 + scale]; 1819 } else { // if ((q3 + scale + 1) <= 35) 1820 res.w[1] = ten2k128[q3 + scale - 20].w[1]; 1821 res.w[0] = ten2k128[q3 + scale - 20].w[0]; 1822 } 1823 res.w[0] = res.w[0] - 1; // borrow never occurs 1824 z_exp = z_exp - EXP_P1; 1825 e3 = e3 - 1; 1826 res.w[1] = z_sign | ((UINT64) (e3 + 6176) << 49) | res.w[1]; 1827 } 1828 if (e3 == expmin) { 1829 if (R64 < 5 || (R64 == 5 && !is_inexact_lt_midpoint)) { 1830 ; // result not tiny (in round-to-nearest mode) 1831 } else { 1832 *pfpsf |= UNDERFLOW_EXCEPTION; 1833 } 1834 } 1835 } // end 10^(q3+scale-1) 1836 // set the inexact flag 1837 *pfpsf |= INEXACT_EXCEPTION; 1838 } else { 1839 if (p_sign == z_sign) { 1840 // if (z_sign), set as if for absolute value 1841 is_inexact_lt_midpoint = 1; 1842 } else { // if (p_sign != z_sign) 1843 // if (z_sign), set as if for absolute value 1844 is_inexact_gt_midpoint = 1; 1845 } 1846 *pfpsf |= INEXACT_EXCEPTION; 1847 } 1848 // the result is always inexact => set the inexact flag 1849 // Determine tininess: 1850 // if (exp > expmin) 1851 // the result is not tiny 1852 // else // if exp = emin 1853 // if (q3 + scale < p34) 1854 // the result is tiny 1855 // else // if (q3 + scale = p34) 1856 // if (C3 * 10^scale > 10^33) 1857 // the result is not tiny 1858 // else // if C3 * 10^scale = 10^33 1859 // if (xy * z > 0) 1860 // the result is not tiny 1861 // else // if (xy * z < 0) 1862 // if (z > 0) 1863 // if rnd_mode != RP 1864 // the result is tiny 1865 // else // if RP 1866 // the result is not tiny 1867 // else // if (z < 0) 1868 // if rnd_mode != RM 1869 // the result is tiny 1870 // else // if RM 1871 // the result is not tiny 1872 // endif 1873 // endif 1874 // endif 1875 // endif 1876 // endif 1877 // endif 1878 if ((e3 == expmin && (q3 + scale) < p34) || 1879 (e3 == expmin && (q3 + scale) == p34 && 1880 (res.w[1] & MASK_COEFF) == 0x0000314dc6448d93ull && // 10^33_high 1881 res.w[0] == 0x38c15b0a00000000ull && // 10^33_low 1882 z_sign != p_sign && ((!z_sign && rnd_mode != ROUNDING_UP) || 1883 (z_sign && rnd_mode != ROUNDING_DOWN)))) { 1884 *pfpsf |= UNDERFLOW_EXCEPTION; 1885 } 1886 if (rnd_mode != ROUNDING_TO_NEAREST) { 1887 rounding_correction (rnd_mode, 1888 is_inexact_lt_midpoint, 1889 is_inexact_gt_midpoint, 1890 is_midpoint_lt_even, is_midpoint_gt_even, 1891 e3, &res, pfpsf); 1892 } 1893 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 1894 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 1895 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 1896 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 1897 BID_SWAP128 (res); 1898 BID_RETURN (res) 1899 1900 } else if (p34 == delta) { // Case (1''B) 1901 1902 // because Case (1''A) was treated above, e3 + 6176 >= p34 - q3 1903 // and C3 can be scaled up to p34 digits if needed 1904 1905 // scale C3 to p34 digits if needed 1906 scale = p34 - q3; // 0 <= scale <= p34 - 1 1907 if (scale == 0) { 1908 res.w[1] = C3.w[1]; 1909 res.w[0] = C3.w[0]; 1910 } else if (q3 <= 19) { // z fits in 64 bits 1911 if (scale <= 19) { // 10^scale fits in 64 bits 1912 // 64 x 64 C3.w[0] * ten2k64[scale] 1913 __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]); 1914 } else { // 10^scale fits in 128 bits 1915 // 64 x 128 C3.w[0] * ten2k128[scale - 20] 1916 __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]); 1917 } 1918 } else { // z fits in 128 bits, but 10^scale must fit in 64 bits 1919 // 64 x 128 ten2k64[scale] * C3 1920 __mul_128x64_to_128 (res, ten2k64[scale], C3); 1921 } 1922 // subtract scale from the exponent 1923 z_exp = z_exp - ((UINT64) scale << 49); 1924 e3 = e3 - scale; 1925 // now z_sign, z_exp, and res correspond to a z scaled to p34 = 34 digits 1926 1927 // determine whether x * y is less than, equal to, or greater than 1928 // 1/2 ulp (z) 1929 if (q4 <= 19) { 1930 if (C4.w[0] < midpoint64[q4 - 1]) { // < 1/2 ulp 1931 lt_half_ulp = 1; 1932 } else if (C4.w[0] == midpoint64[q4 - 1]) { // = 1/2 ulp 1933 eq_half_ulp = 1; 1934 } else { // > 1/2 ulp 1935 gt_half_ulp = 1; 1936 } 1937 } else if (q4 <= 38) { 1938 if (C4.w[2] == 0 && (C4.w[1] < midpoint128[q4 - 20].w[1] || 1939 (C4.w[1] == midpoint128[q4 - 20].w[1] && 1940 C4.w[0] < midpoint128[q4 - 20].w[0]))) { // < 1/2 ulp 1941 lt_half_ulp = 1; 1942 } else if (C4.w[2] == 0 && C4.w[1] == midpoint128[q4 - 20].w[1] && 1943 C4.w[0] == midpoint128[q4 - 20].w[0]) { // = 1/2 ulp 1944 eq_half_ulp = 1; 1945 } else { // > 1/2 ulp 1946 gt_half_ulp = 1; 1947 } 1948 } else if (q4 <= 58) { 1949 if (C4.w[3] == 0 && (C4.w[2] < midpoint192[q4 - 39].w[2] || 1950 (C4.w[2] == midpoint192[q4 - 39].w[2] && 1951 C4.w[1] < midpoint192[q4 - 39].w[1]) || 1952 (C4.w[2] == midpoint192[q4 - 39].w[2] && 1953 C4.w[1] == midpoint192[q4 - 39].w[1] && 1954 C4.w[0] < midpoint192[q4 - 39].w[0]))) { // < 1/2 ulp 1955 lt_half_ulp = 1; 1956 } else if (C4.w[3] == 0 && C4.w[2] == midpoint192[q4 - 39].w[2] && 1957 C4.w[1] == midpoint192[q4 - 39].w[1] && 1958 C4.w[0] == midpoint192[q4 - 39].w[0]) { // = 1/2 ulp 1959 eq_half_ulp = 1; 1960 } else { // > 1/2 ulp 1961 gt_half_ulp = 1; 1962 } 1963 } else { 1964 if (C4.w[3] < midpoint256[q4 - 59].w[3] || 1965 (C4.w[3] == midpoint256[q4 - 59].w[3] && 1966 C4.w[2] < midpoint256[q4 - 59].w[2]) || 1967 (C4.w[3] == midpoint256[q4 - 59].w[3] && 1968 C4.w[2] == midpoint256[q4 - 59].w[2] && 1969 C4.w[1] < midpoint256[q4 - 59].w[1]) || 1970 (C4.w[3] == midpoint256[q4 - 59].w[3] && 1971 C4.w[2] == midpoint256[q4 - 59].w[2] && 1972 C4.w[1] == midpoint256[q4 - 59].w[1] && 1973 C4.w[0] < midpoint256[q4 - 59].w[0])) { // < 1/2 ulp 1974 lt_half_ulp = 1; 1975 } else if (C4.w[3] == midpoint256[q4 - 59].w[3] && 1976 C4.w[2] == midpoint256[q4 - 59].w[2] && 1977 C4.w[1] == midpoint256[q4 - 59].w[1] && 1978 C4.w[0] == midpoint256[q4 - 59].w[0]) { // = 1/2 ulp 1979 eq_half_ulp = 1; 1980 } else { // > 1/2 ulp 1981 gt_half_ulp = 1; 1982 } 1983 } 1984 1985 if (p_sign == z_sign) { 1986 if (lt_half_ulp) { 1987 res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1]; 1988 // use the following to avoid double rounding errors when operating on 1989 // mixed formats in rounding to nearest 1990 is_inexact_lt_midpoint = 1; // if (z_sign), as if for absolute value 1991 } else if ((eq_half_ulp && (res.w[0] & 0x01)) || gt_half_ulp) { 1992 // add 1 ulp to the significand 1993 res.w[0]++; 1994 if (res.w[0] == 0x0ull) 1995 res.w[1]++; 1996 // check for rounding overflow, when coeff == 10^34 1997 if ((res.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull && 1998 res.w[0] == 0x378d8e6400000000ull) { // coefficient = 10^34 1999 e3 = e3 + 1; 2000 // coeff = 10^33 2001 z_exp = ((UINT64) (e3 + 6176) << 49) & MASK_EXP; 2002 res.w[1] = 0x0000314dc6448d93ull; 2003 res.w[0] = 0x38c15b0a00000000ull; 2004 } 2005 // end add 1 ulp 2006 res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1]; 2007 if (eq_half_ulp) { 2008 is_midpoint_lt_even = 1; // if (z_sign), as if for absolute value 2009 } else { 2010 is_inexact_gt_midpoint = 1; // if (z_sign), as if for absolute value 2011 } 2012 } else { // if (eq_half_ulp && !(res.w[0] & 0x01)) 2013 // leave unchanged 2014 res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1]; 2015 is_midpoint_gt_even = 1; // if (z_sign), as if for absolute value 2016 } 2017 // the result is always inexact, and never tiny 2018 // set the inexact flag 2019 *pfpsf |= INEXACT_EXCEPTION; 2020 // check for overflow 2021 if (e3 > expmax && rnd_mode == ROUNDING_TO_NEAREST) { 2022 res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf 2023 res.w[0] = 0x0000000000000000ull; 2024 *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION); 2025 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 2026 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 2027 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 2028 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 2029 BID_SWAP128 (res); 2030 BID_RETURN (res) 2031 } 2032 if (rnd_mode != ROUNDING_TO_NEAREST) { 2033 rounding_correction (rnd_mode, 2034 is_inexact_lt_midpoint, 2035 is_inexact_gt_midpoint, 2036 is_midpoint_lt_even, is_midpoint_gt_even, 2037 e3, &res, pfpsf); 2038 z_exp = res.w[1] & MASK_EXP; 2039 } 2040 } else { // if (p_sign != z_sign) 2041 // consider two cases, because C3 * 10^scale = 10^33 is a special case 2042 if (res.w[1] != 0x0000314dc6448d93ull || 2043 res.w[0] != 0x38c15b0a00000000ull) { // C3 * 10^scale != 10^33 2044 if (lt_half_ulp) { 2045 res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1]; 2046 // use the following to avoid double rounding errors when operating 2047 // on mixed formats in rounding to nearest 2048 is_inexact_gt_midpoint = 1; // if (z_sign), as if for absolute value 2049 } else if ((eq_half_ulp && (res.w[0] & 0x01)) || gt_half_ulp) { 2050 // subtract 1 ulp from the significand 2051 res.w[0]--; 2052 if (res.w[0] == 0xffffffffffffffffull) 2053 res.w[1]--; 2054 res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1]; 2055 if (eq_half_ulp) { 2056 is_midpoint_gt_even = 1; // if (z_sign), as if for absolute value 2057 } else { 2058 is_inexact_lt_midpoint = 1; //if(z_sign), as if for absolute value 2059 } 2060 } else { // if (eq_half_ulp && !(res.w[0] & 0x01)) 2061 // leave unchanged 2062 res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1]; 2063 is_midpoint_lt_even = 1; // if (z_sign), as if for absolute value 2064 } 2065 // the result is always inexact, and never tiny 2066 // check for overflow for RN 2067 if (e3 > expmax) { 2068 if (rnd_mode == ROUNDING_TO_NEAREST) { 2069 res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf 2070 res.w[0] = 0x0000000000000000ull; 2071 *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION); 2072 } else { 2073 rounding_correction (rnd_mode, 2074 is_inexact_lt_midpoint, 2075 is_inexact_gt_midpoint, 2076 is_midpoint_lt_even, 2077 is_midpoint_gt_even, e3, &res, 2078 pfpsf); 2079 } 2080 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 2081 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 2082 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 2083 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 2084 BID_SWAP128 (res); 2085 BID_RETURN (res) 2086 } 2087 // set the inexact flag 2088 *pfpsf |= INEXACT_EXCEPTION; 2089 if (rnd_mode != ROUNDING_TO_NEAREST) { 2090 rounding_correction (rnd_mode, 2091 is_inexact_lt_midpoint, 2092 is_inexact_gt_midpoint, 2093 is_midpoint_lt_even, 2094 is_midpoint_gt_even, e3, &res, pfpsf); 2095 } 2096 z_exp = res.w[1] & MASK_EXP; 2097 } else { // if C3 * 10^scale = 10^33 2098 e3 = (z_exp >> 49) - 6176; 2099 if (e3 > expmin) { 2100 // the result is exact if exp > expmin and C4 = d*10^(q4-1), 2101 // where d = 1, 2, 3, ..., 9; it could be tiny too, but exact 2102 if (q4 == 1) { 2103 // if q4 = 1 the result is exact 2104 // result coefficient = 10^34 - C4 2105 res.w[1] = 0x0001ed09bead87c0ull; 2106 res.w[0] = 0x378d8e6400000000ull - C4.w[0]; 2107 z_exp = z_exp - EXP_P1; 2108 e3 = e3 - 1; 2109 res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1]; 2110 } else { 2111 // if q4 > 1 then truncate C4 from q4 digits to 1 digit; 2112 // x = q4-1, 1 <= x <= 67 and check if this operation is exact 2113 if (q4 <= 18) { // 2 <= q4 <= 18 2114 round64_2_18 (q4, q4 - 1, C4.w[0], &R64, &incr_exp, 2115 &is_midpoint_lt_even, 2116 &is_midpoint_gt_even, 2117 &is_inexact_lt_midpoint, 2118 &is_inexact_gt_midpoint); 2119 } else if (q4 <= 38) { 2120 P128.w[1] = C4.w[1]; 2121 P128.w[0] = C4.w[0]; 2122 round128_19_38 (q4, q4 - 1, P128, &R128, &incr_exp, 2123 &is_midpoint_lt_even, 2124 &is_midpoint_gt_even, 2125 &is_inexact_lt_midpoint, 2126 &is_inexact_gt_midpoint); 2127 R64 = R128.w[0]; // one decimal digit 2128 } else if (q4 <= 57) { 2129 P192.w[2] = C4.w[2]; 2130 P192.w[1] = C4.w[1]; 2131 P192.w[0] = C4.w[0]; 2132 round192_39_57 (q4, q4 - 1, P192, &R192, &incr_exp, 2133 &is_midpoint_lt_even, 2134 &is_midpoint_gt_even, 2135 &is_inexact_lt_midpoint, 2136 &is_inexact_gt_midpoint); 2137 R64 = R192.w[0]; // one decimal digit 2138 } else { // if (q4 <= 68) 2139 round256_58_76 (q4, q4 - 1, C4, &R256, &incr_exp, 2140 &is_midpoint_lt_even, 2141 &is_midpoint_gt_even, 2142 &is_inexact_lt_midpoint, 2143 &is_inexact_gt_midpoint); 2144 R64 = R256.w[0]; // one decimal digit 2145 } 2146 if (!is_midpoint_lt_even && !is_midpoint_gt_even && 2147 !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) { 2148 // the result is exact: 10^34 - R64 2149 // incr_exp = 0 with certainty 2150 z_exp = z_exp - EXP_P1; 2151 e3 = e3 - 1; 2152 res.w[1] = 2153 z_sign | (z_exp & MASK_EXP) | 0x0001ed09bead87c0ull; 2154 res.w[0] = 0x378d8e6400000000ull - R64; 2155 } else { 2156 // We want R64 to be the top digit of C4, but we actually 2157 // obtained (C4 * 10^(-q4+1))RN; a correction may be needed, 2158 // because the top digit is (C4 * 10^(-q4+1))RZ 2159 // however, if incr_exp = 1 then R64 = 10 with certainty 2160 if (incr_exp) { 2161 R64 = 10; 2162 } 2163 // the result is inexact as C4 has more than 1 significant digit 2164 // and C3 * 10^scale = 10^33 2165 // example of case that is treated here: 2166 // 100...0 * 10^e3 - 0.41 * 10^e3 = 2167 // 0999...9.59 * 10^e3 -> rounds to 99...96*10^(e3-1) 2168 // note that (e3 > expmin} 2169 // in order to round, subtract R64 from 10^34 and then compare 2170 // C4 - R64 * 10^(q4-1) with 1/2 ulp 2171 // calculate 10^34 - R64 2172 res.w[1] = 0x0001ed09bead87c0ull; 2173 res.w[0] = 0x378d8e6400000000ull - R64; 2174 z_exp = z_exp - EXP_P1; // will be OR-ed with sign & significand 2175 // calculate C4 - R64 * 10^(q4-1); this is a rare case and 2176 // R64 is small, 1 <= R64 <= 9 2177 e3 = e3 - 1; 2178 if (is_inexact_lt_midpoint) { 2179 is_inexact_lt_midpoint = 0; 2180 is_inexact_gt_midpoint = 1; 2181 } else if (is_inexact_gt_midpoint) { 2182 is_inexact_gt_midpoint = 0; 2183 is_inexact_lt_midpoint = 1; 2184 } else if (is_midpoint_lt_even) { 2185 is_midpoint_lt_even = 0; 2186 is_midpoint_gt_even = 1; 2187 } else if (is_midpoint_gt_even) { 2188 is_midpoint_gt_even = 0; 2189 is_midpoint_lt_even = 1; 2190 } else { 2191 ; 2192 } 2193 // the result is always inexact, and never tiny 2194 // check for overflow for RN 2195 if (e3 > expmax) { 2196 if (rnd_mode == ROUNDING_TO_NEAREST) { 2197 res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf 2198 res.w[0] = 0x0000000000000000ull; 2199 *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION); 2200 } else { 2201 rounding_correction (rnd_mode, 2202 is_inexact_lt_midpoint, 2203 is_inexact_gt_midpoint, 2204 is_midpoint_lt_even, 2205 is_midpoint_gt_even, e3, &res, 2206 pfpsf); 2207 } 2208 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 2209 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 2210 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 2211 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 2212 BID_SWAP128 (res); 2213 BID_RETURN (res) 2214 } 2215 // set the inexact flag 2216 *pfpsf |= INEXACT_EXCEPTION; 2217 res.w[1] = 2218 z_sign | ((UINT64) (e3 + 6176) << 49) | res.w[1]; 2219 if (rnd_mode != ROUNDING_TO_NEAREST) { 2220 rounding_correction (rnd_mode, 2221 is_inexact_lt_midpoint, 2222 is_inexact_gt_midpoint, 2223 is_midpoint_lt_even, 2224 is_midpoint_gt_even, e3, &res, 2225 pfpsf); 2226 } 2227 z_exp = res.w[1] & MASK_EXP; 2228 } // end result is inexact 2229 } // end q4 > 1 2230 } else { // if (e3 = emin) 2231 // if e3 = expmin the result is also tiny (the condition for 2232 // tininess is C4 > 050...0 [q4 digits] which is met because 2233 // the msd of C4 is not zero) 2234 // the result is tiny and inexact in all rounding modes; 2235 // it is either 100...0 or 0999...9 (use lt_half_ulp, eq_half_ulp, 2236 // gt_half_ulp to calculate) 2237 // if (lt_half_ulp || eq_half_ulp) res = 10^33 stays unchanged 2238 2239 // p_sign != z_sign so swap gt_half_ulp and lt_half_ulp 2240 if (gt_half_ulp) { // res = 10^33 - 1 2241 res.w[1] = 0x0000314dc6448d93ull; 2242 res.w[0] = 0x38c15b09ffffffffull; 2243 } else { 2244 res.w[1] = 0x0000314dc6448d93ull; 2245 res.w[0] = 0x38c15b0a00000000ull; 2246 } 2247 res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1]; 2248 *pfpsf |= UNDERFLOW_EXCEPTION; // inexact is set later 2249 2250 if (eq_half_ulp) { 2251 is_midpoint_lt_even = 1; // if (z_sign), as if for absolute value 2252 } else if (lt_half_ulp) { 2253 is_inexact_gt_midpoint = 1; //if(z_sign), as if for absolute value 2254 } else { // if (gt_half_ulp) 2255 is_inexact_lt_midpoint = 1; //if(z_sign), as if for absolute value 2256 } 2257 2258 if (rnd_mode != ROUNDING_TO_NEAREST) { 2259 rounding_correction (rnd_mode, 2260 is_inexact_lt_midpoint, 2261 is_inexact_gt_midpoint, 2262 is_midpoint_lt_even, 2263 is_midpoint_gt_even, e3, &res, 2264 pfpsf); 2265 z_exp = res.w[1] & MASK_EXP; 2266 } 2267 } // end e3 = emin 2268 // set the inexact flag (if the result was not exact) 2269 if (is_inexact_lt_midpoint || is_inexact_gt_midpoint || 2270 is_midpoint_lt_even || is_midpoint_gt_even) 2271 *pfpsf |= INEXACT_EXCEPTION; 2272 } // end 10^33 2273 } // end if (p_sign != z_sign) 2274 res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1]; 2275 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 2276 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 2277 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 2278 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 2279 BID_SWAP128 (res); 2280 BID_RETURN (res) 2281 2282 } else if (((q3 <= delta && delta < p34 && p34 < delta + q4) || // Case (2) 2283 (q3 <= delta && delta + q4 <= p34) || // Case (3) 2284 (delta < q3 && p34 < delta + q4) || // Case (4) 2285 (delta < q3 && q3 <= delta + q4 && delta + q4 <= p34) || // Case (5) 2286 (delta + q4 < q3)) && // Case (6) 2287 !(delta <= 1 && p_sign != z_sign)) { // Case (2), (3), (4), (5) or (6) 2288 2289 // the result has the sign of z 2290 2291 if ((q3 <= delta && delta < p34 && p34 < delta + q4) || // Case (2) 2292 (delta < q3 && p34 < delta + q4)) { // Case (4) 2293 // round first the sum x * y + z with unbounded exponent 2294 // scale C3 up by scale = p34 - q3, 1 <= scale <= p34-1, 2295 // 1 <= scale <= 33 2296 // calculate res = C3 * 10^scale 2297 scale = p34 - q3; 2298 x0 = delta + q4 - p34; 2299 } else if (delta + q4 < q3) { // Case (6) 2300 // make Case (6) look like Case (3) or Case (5) with scale = 0 2301 // by scaling up C4 by 10^(q3 - delta - q4) 2302 scale = q3 - delta - q4; // 1 <= scale <= 33 2303 if (q4 <= 19) { // 1 <= scale <= 19; C4 fits in 64 bits 2304 if (scale <= 19) { // 10^scale fits in 64 bits 2305 // 64 x 64 C4.w[0] * ten2k64[scale] 2306 __mul_64x64_to_128MACH (P128, C4.w[0], ten2k64[scale]); 2307 } else { // 10^scale fits in 128 bits 2308 // 64 x 128 C4.w[0] * ten2k128[scale - 20] 2309 __mul_128x64_to_128 (P128, C4.w[0], ten2k128[scale - 20]); 2310 } 2311 } else { // C4 fits in 128 bits, but 10^scale must fit in 64 bits 2312 // 64 x 128 ten2k64[scale] * C4 2313 __mul_128x64_to_128 (P128, ten2k64[scale], C4); 2314 } 2315 C4.w[0] = P128.w[0]; 2316 C4.w[1] = P128.w[1]; 2317 // e4 does not need adjustment, as it is not used from this point on 2318 scale = 0; 2319 x0 = 0; 2320 // now Case (6) looks like Case (3) or Case (5) with scale = 0 2321 } else { // if Case (3) or Case (5) 2322 // Note: Case (3) is similar to Case (2), but scale differs and the 2323 // result is exact, unless it is tiny (so x0 = 0 when calculating the 2324 // result with unbounded exponent) 2325 2326 // calculate first the sum x * y + z with unbounded exponent (exact) 2327 // scale C3 up by scale = delta + q4 - q3, 1 <= scale <= p34-1, 2328 // 1 <= scale <= 33 2329 // calculate res = C3 * 10^scale 2330 scale = delta + q4 - q3; 2331 x0 = 0; 2332 // Note: the comments which follow refer [mainly] to Case (2)] 2333 } 2334 2335 case2_repeat: 2336 if (scale == 0) { // this could happen e.g. if we return to case2_repeat 2337 // or in Case (4) 2338 res.w[1] = C3.w[1]; 2339 res.w[0] = C3.w[0]; 2340 } else if (q3 <= 19) { // 1 <= scale <= 19; z fits in 64 bits 2341 if (scale <= 19) { // 10^scale fits in 64 bits 2342 // 64 x 64 C3.w[0] * ten2k64[scale] 2343 __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]); 2344 } else { // 10^scale fits in 128 bits 2345 // 64 x 128 C3.w[0] * ten2k128[scale - 20] 2346 __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]); 2347 } 2348 } else { // z fits in 128 bits, but 10^scale must fit in 64 bits 2349 // 64 x 128 ten2k64[scale] * C3 2350 __mul_128x64_to_128 (res, ten2k64[scale], C3); 2351 } 2352 // e3 is already calculated 2353 e3 = e3 - scale; 2354 // now res = C3 * 10^scale and e3 = e3 - scale 2355 // Note: C3 * 10^scale could be 10^34 if we returned to case2_repeat 2356 // because the result was too small 2357 2358 // round C4 to nearest to q4 - x0 digits, where x0 = delta + q4 - p34, 2359 // 1 <= x0 <= min (q4 - 1, 2 * p34 - 1) <=> 1 <= x0 <= min (q4 - 1, 67) 2360 // Also: 1 <= q4 - x0 <= p34 -1 => 1 <= q4 - x0 <= 33 (so the result of 2361 // the rounding fits in 128 bits!) 2362 // x0 = delta + q4 - p34 (calculated before reaching case2_repeat) 2363 // because q3 + q4 - x0 <= P => x0 >= q3 + q4 - p34 2364 if (x0 == 0) { // this could happen only if we return to case2_repeat, or 2365 // for Case (3) or Case (6) 2366 R128.w[1] = C4.w[1]; 2367 R128.w[0] = C4.w[0]; 2368 } else if (q4 <= 18) { 2369 // 2 <= q4 <= 18, max(1, q3+q4-p34) <= x0 <= q4 - 1, 1 <= x0 <= 17 2370 round64_2_18 (q4, x0, C4.w[0], &R64, &incr_exp, 2371 &is_midpoint_lt_even, &is_midpoint_gt_even, 2372 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint); 2373 if (incr_exp) { 2374 // R64 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 17 2375 R64 = ten2k64[q4 - x0]; 2376 } 2377 R128.w[1] = 0; 2378 R128.w[0] = R64; 2379 } else if (q4 <= 38) { 2380 // 19 <= q4 <= 38, max(1, q3+q4-p34) <= x0 <= q4 - 1, 1 <= x0 <= 37 2381 P128.w[1] = C4.w[1]; 2382 P128.w[0] = C4.w[0]; 2383 round128_19_38 (q4, x0, P128, &R128, &incr_exp, 2384 &is_midpoint_lt_even, &is_midpoint_gt_even, 2385 &is_inexact_lt_midpoint, 2386 &is_inexact_gt_midpoint); 2387 if (incr_exp) { 2388 // R128 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 37 2389 if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19 2390 R128.w[0] = ten2k64[q4 - x0]; 2391 // R128.w[1] stays 0 2392 } else { // 20 <= q4 - x0 <= 37 2393 R128.w[0] = ten2k128[q4 - x0 - 20].w[0]; 2394 R128.w[1] = ten2k128[q4 - x0 - 20].w[1]; 2395 } 2396 } 2397 } else if (q4 <= 57) { 2398 // 38 <= q4 <= 57, max(1, q3+q4-p34) <= x0 <= q4 - 1, 5 <= x0 <= 56 2399 P192.w[2] = C4.w[2]; 2400 P192.w[1] = C4.w[1]; 2401 P192.w[0] = C4.w[0]; 2402 round192_39_57 (q4, x0, P192, &R192, &incr_exp, 2403 &is_midpoint_lt_even, &is_midpoint_gt_even, 2404 &is_inexact_lt_midpoint, 2405 &is_inexact_gt_midpoint); 2406 // R192.w[2] is always 0 2407 if (incr_exp) { 2408 // R192 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 5, 1 <= q4 - x0 <= 52 2409 if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19 2410 R192.w[0] = ten2k64[q4 - x0]; 2411 // R192.w[1] stays 0 2412 // R192.w[2] stays 0 2413 } else { // 20 <= q4 - x0 <= 33 2414 R192.w[0] = ten2k128[q4 - x0 - 20].w[0]; 2415 R192.w[1] = ten2k128[q4 - x0 - 20].w[1]; 2416 // R192.w[2] stays 0 2417 } 2418 } 2419 R128.w[1] = R192.w[1]; 2420 R128.w[0] = R192.w[0]; 2421 } else { 2422 // 58 <= q4 <= 68, max(1, q3+q4-p34) <= x0 <= q4 - 1, 25 <= x0 <= 67 2423 round256_58_76 (q4, x0, C4, &R256, &incr_exp, 2424 &is_midpoint_lt_even, &is_midpoint_gt_even, 2425 &is_inexact_lt_midpoint, 2426 &is_inexact_gt_midpoint); 2427 // R256.w[3] and R256.w[2] are always 0 2428 if (incr_exp) { 2429 // R256 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 25, 1 <= q4 - x0 <= 43 2430 if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19 2431 R256.w[0] = ten2k64[q4 - x0]; 2432 // R256.w[1] stays 0 2433 // R256.w[2] stays 0 2434 // R256.w[3] stays 0 2435 } else { // 20 <= q4 - x0 <= 33 2436 R256.w[0] = ten2k128[q4 - x0 - 20].w[0]; 2437 R256.w[1] = ten2k128[q4 - x0 - 20].w[1]; 2438 // R256.w[2] stays 0 2439 // R256.w[3] stays 0 2440 } 2441 } 2442 R128.w[1] = R256.w[1]; 2443 R128.w[0] = R256.w[0]; 2444 } 2445 // now add C3 * 10^scale in res and the signed top (q4-x0) digits of C4, 2446 // rounded to nearest, which were copied into R128 2447 if (z_sign == p_sign) { 2448 lsb = res.w[0] & 0x01; // lsb of C3 * 10^scale 2449 // the sum can result in [up to] p34 or p34 + 1 digits 2450 res.w[0] = res.w[0] + R128.w[0]; 2451 res.w[1] = res.w[1] + R128.w[1]; 2452 if (res.w[0] < R128.w[0]) 2453 res.w[1]++; // carry 2454 // if res > 10^34 - 1 need to increase x0 and decrease scale by 1 2455 if (res.w[1] > 0x0001ed09bead87c0ull || 2456 (res.w[1] == 0x0001ed09bead87c0ull && 2457 res.w[0] > 0x378d8e63ffffffffull)) { 2458 // avoid double rounding error 2459 is_inexact_lt_midpoint0 = is_inexact_lt_midpoint; 2460 is_inexact_gt_midpoint0 = is_inexact_gt_midpoint; 2461 is_midpoint_lt_even0 = is_midpoint_lt_even; 2462 is_midpoint_gt_even0 = is_midpoint_gt_even; 2463 is_inexact_lt_midpoint = 0; 2464 is_inexact_gt_midpoint = 0; 2465 is_midpoint_lt_even = 0; 2466 is_midpoint_gt_even = 0; 2467 P128.w[1] = res.w[1]; 2468 P128.w[0] = res.w[0]; 2469 round128_19_38 (35, 1, P128, &res, &incr_exp, 2470 &is_midpoint_lt_even, &is_midpoint_gt_even, 2471 &is_inexact_lt_midpoint, 2472 &is_inexact_gt_midpoint); 2473 // incr_exp is 0 with certainty in this case 2474 // avoid a double rounding error 2475 if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 2476 is_midpoint_lt_even) { // double rounding error upward 2477 // res = res - 1 2478 res.w[0]--; 2479 if (res.w[0] == 0xffffffffffffffffull) 2480 res.w[1]--; 2481 // Note: a double rounding error upward is not possible; for this 2482 // the result after the first rounding would have to be 99...95 2483 // (35 digits in all), possibly followed by a number of zeros; this 2484 // not possible in Cases (2)-(6) or (15)-(17) which may get here 2485 is_midpoint_lt_even = 0; 2486 is_inexact_lt_midpoint = 1; 2487 } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 2488 is_midpoint_gt_even) { // double rounding error downward 2489 // res = res + 1 2490 res.w[0]++; 2491 if (res.w[0] == 0) 2492 res.w[1]++; 2493 is_midpoint_gt_even = 0; 2494 is_inexact_gt_midpoint = 1; 2495 } else if (!is_midpoint_lt_even && !is_midpoint_gt_even && 2496 !is_inexact_lt_midpoint 2497 && !is_inexact_gt_midpoint) { 2498 // if this second rounding was exact the result may still be 2499 // inexact because of the first rounding 2500 if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) { 2501 is_inexact_gt_midpoint = 1; 2502 } 2503 if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) { 2504 is_inexact_lt_midpoint = 1; 2505 } 2506 } else if (is_midpoint_gt_even && 2507 (is_inexact_gt_midpoint0 2508 || is_midpoint_lt_even0)) { 2509 // pulled up to a midpoint 2510 is_inexact_lt_midpoint = 1; 2511 is_inexact_gt_midpoint = 0; 2512 is_midpoint_lt_even = 0; 2513 is_midpoint_gt_even = 0; 2514 } else if (is_midpoint_lt_even && 2515 (is_inexact_lt_midpoint0 2516 || is_midpoint_gt_even0)) { 2517 // pulled down to a midpoint 2518 is_inexact_lt_midpoint = 0; 2519 is_inexact_gt_midpoint = 1; 2520 is_midpoint_lt_even = 0; 2521 is_midpoint_gt_even = 0; 2522 } else { 2523 ; 2524 } 2525 // adjust exponent 2526 e3 = e3 + 1; 2527 if (!is_midpoint_lt_even && !is_midpoint_gt_even && 2528 !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) { 2529 if (is_midpoint_lt_even0 || is_midpoint_gt_even0 || 2530 is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0) { 2531 is_inexact_lt_midpoint = 1; 2532 } 2533 } 2534 } else { 2535 // this is the result rounded with unbounded exponent, unless a 2536 // correction is needed 2537 res.w[1] = res.w[1] & MASK_COEFF; 2538 if (lsb == 1) { 2539 if (is_midpoint_gt_even) { 2540 // res = res + 1 2541 is_midpoint_gt_even = 0; 2542 is_midpoint_lt_even = 1; 2543 res.w[0]++; 2544 if (res.w[0] == 0x0) 2545 res.w[1]++; 2546 // check for rounding overflow 2547 if (res.w[1] == 0x0001ed09bead87c0ull && 2548 res.w[0] == 0x378d8e6400000000ull) { 2549 // res = 10^34 => rounding overflow 2550 res.w[1] = 0x0000314dc6448d93ull; 2551 res.w[0] = 0x38c15b0a00000000ull; // 10^33 2552 e3++; 2553 } 2554 } else if (is_midpoint_lt_even) { 2555 // res = res - 1 2556 is_midpoint_lt_even = 0; 2557 is_midpoint_gt_even = 1; 2558 res.w[0]--; 2559 if (res.w[0] == 0xffffffffffffffffull) 2560 res.w[1]--; 2561 // if the result is pure zero, the sign depends on the rounding 2562 // mode (x*y and z had opposite signs) 2563 if (res.w[1] == 0x0ull && res.w[0] == 0x0ull) { 2564 if (rnd_mode != ROUNDING_DOWN) 2565 z_sign = 0x0000000000000000ull; 2566 else 2567 z_sign = 0x8000000000000000ull; 2568 // the exponent is max (e3, expmin) 2569 res.w[1] = 0x0; 2570 res.w[0] = 0x0; 2571 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 2572 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 2573 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 2574 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 2575 BID_SWAP128 (res); 2576 BID_RETURN (res) 2577 } 2578 } else { 2579 ; 2580 } 2581 } 2582 } 2583 } else { // if (z_sign != p_sign) 2584 lsb = res.w[0] & 0x01; // lsb of C3 * 10^scale; R128 contains rounded C4 2585 // used to swap rounding indicators if p_sign != z_sign 2586 // the sum can result in [up to] p34 or p34 - 1 digits 2587 tmp64 = res.w[0]; 2588 res.w[0] = res.w[0] - R128.w[0]; 2589 res.w[1] = res.w[1] - R128.w[1]; 2590 if (res.w[0] > tmp64) 2591 res.w[1]--; // borrow 2592 // if res < 10^33 and exp > expmin need to decrease x0 and 2593 // increase scale by 1 2594 if (e3 > expmin && ((res.w[1] < 0x0000314dc6448d93ull || 2595 (res.w[1] == 0x0000314dc6448d93ull && 2596 res.w[0] < 0x38c15b0a00000000ull)) || 2597 (is_inexact_lt_midpoint 2598 && res.w[1] == 0x0000314dc6448d93ull 2599 && res.w[0] == 0x38c15b0a00000000ull)) 2600 && x0 >= 1) { 2601 x0 = x0 - 1; 2602 // first restore e3, otherwise it will be too small 2603 e3 = e3 + scale; 2604 scale = scale + 1; 2605 is_inexact_lt_midpoint = 0; 2606 is_inexact_gt_midpoint = 0; 2607 is_midpoint_lt_even = 0; 2608 is_midpoint_gt_even = 0; 2609 incr_exp = 0; 2610 goto case2_repeat; 2611 } 2612 // else this is the result rounded with unbounded exponent; 2613 // because the result has opposite sign to that of C4 which was 2614 // rounded, need to change the rounding indicators 2615 if (is_inexact_lt_midpoint) { 2616 is_inexact_lt_midpoint = 0; 2617 is_inexact_gt_midpoint = 1; 2618 } else if (is_inexact_gt_midpoint) { 2619 is_inexact_gt_midpoint = 0; 2620 is_inexact_lt_midpoint = 1; 2621 } else if (lsb == 0) { 2622 if (is_midpoint_lt_even) { 2623 is_midpoint_lt_even = 0; 2624 is_midpoint_gt_even = 1; 2625 } else if (is_midpoint_gt_even) { 2626 is_midpoint_gt_even = 0; 2627 is_midpoint_lt_even = 1; 2628 } else { 2629 ; 2630 } 2631 } else if (lsb == 1) { 2632 if (is_midpoint_lt_even) { 2633 // res = res + 1 2634 res.w[0]++; 2635 if (res.w[0] == 0x0) 2636 res.w[1]++; 2637 // check for rounding overflow 2638 if (res.w[1] == 0x0001ed09bead87c0ull && 2639 res.w[0] == 0x378d8e6400000000ull) { 2640 // res = 10^34 => rounding overflow 2641 res.w[1] = 0x0000314dc6448d93ull; 2642 res.w[0] = 0x38c15b0a00000000ull; // 10^33 2643 e3++; 2644 } 2645 } else if (is_midpoint_gt_even) { 2646 // res = res - 1 2647 res.w[0]--; 2648 if (res.w[0] == 0xffffffffffffffffull) 2649 res.w[1]--; 2650 // if the result is pure zero, the sign depends on the rounding 2651 // mode (x*y and z had opposite signs) 2652 if (res.w[1] == 0x0ull && res.w[0] == 0x0ull) { 2653 if (rnd_mode != ROUNDING_DOWN) 2654 z_sign = 0x0000000000000000ull; 2655 else 2656 z_sign = 0x8000000000000000ull; 2657 // the exponent is max (e3, expmin) 2658 res.w[1] = 0x0; 2659 res.w[0] = 0x0; 2660 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 2661 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 2662 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 2663 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 2664 BID_SWAP128 (res); 2665 BID_RETURN (res) 2666 } 2667 } else { 2668 ; 2669 } 2670 } else { 2671 ; 2672 } 2673 } 2674 // check for underflow 2675 if (e3 == expmin) { // and if significand < 10^33 => result is tiny 2676 if ((res.w[1] & MASK_COEFF) < 0x0000314dc6448d93ull || 2677 ((res.w[1] & MASK_COEFF) == 0x0000314dc6448d93ull && 2678 res.w[0] < 0x38c15b0a00000000ull)) { 2679 is_tiny = 1; 2680 } 2681 } else if (e3 < expmin) { 2682 // the result is tiny, so we must truncate more of res 2683 is_tiny = 1; 2684 x0 = expmin - e3; 2685 is_inexact_lt_midpoint0 = is_inexact_lt_midpoint; 2686 is_inexact_gt_midpoint0 = is_inexact_gt_midpoint; 2687 is_midpoint_lt_even0 = is_midpoint_lt_even; 2688 is_midpoint_gt_even0 = is_midpoint_gt_even; 2689 is_inexact_lt_midpoint = 0; 2690 is_inexact_gt_midpoint = 0; 2691 is_midpoint_lt_even = 0; 2692 is_midpoint_gt_even = 0; 2693 // determine the number of decimal digits in res 2694 if (res.w[1] == 0x0) { 2695 // between 1 and 19 digits 2696 for (ind = 1; ind <= 19; ind++) { 2697 if (res.w[0] < ten2k64[ind]) { 2698 break; 2699 } 2700 } 2701 // ind digits 2702 } else if (res.w[1] < ten2k128[0].w[1] || 2703 (res.w[1] == ten2k128[0].w[1] 2704 && res.w[0] < ten2k128[0].w[0])) { 2705 // 20 digits 2706 ind = 20; 2707 } else { // between 21 and 38 digits 2708 for (ind = 1; ind <= 18; ind++) { 2709 if (res.w[1] < ten2k128[ind].w[1] || 2710 (res.w[1] == ten2k128[ind].w[1] && 2711 res.w[0] < ten2k128[ind].w[0])) { 2712 break; 2713 } 2714 } 2715 // ind + 20 digits 2716 ind = ind + 20; 2717 } 2718 2719 // at this point ind >= x0; because delta >= 2 on this path, the case 2720 // ind = x0 can occur only in Case (2) or case (3), when C3 has one 2721 // digit (q3 = 1) equal to 1 (C3 = 1), e3 is expmin (e3 = expmin), 2722 // the signs of x * y and z are opposite, and through cancellation 2723 // the most significant decimal digit in res has the weight 2724 // 10^(emin-1); however, it is clear that in this case the most 2725 // significant digit is 9, so the result before rounding is 2726 // 0.9... * 10^emin 2727 // Otherwise, ind > x0 because there are non-zero decimal digits in the 2728 // result with weight of at least 10^emin, and correction for underflow 2729 // can be carried out using the round*_*_2_* () routines 2730 if (x0 == ind) { // the result before rounding is 0.9... * 10^emin 2731 res.w[1] = 0x0; 2732 res.w[0] = 0x1; 2733 is_inexact_gt_midpoint = 1; 2734 } else if (ind <= 18) { // check that 2 <= ind 2735 // 2 <= ind <= 18, 1 <= x0 <= 17 2736 round64_2_18 (ind, x0, res.w[0], &R64, &incr_exp, 2737 &is_midpoint_lt_even, &is_midpoint_gt_even, 2738 &is_inexact_lt_midpoint, 2739 &is_inexact_gt_midpoint); 2740 if (incr_exp) { 2741 // R64 = 10^(ind-x0), 1 <= ind - x0 <= ind - 1, 1 <= ind - x0 <= 17 2742 R64 = ten2k64[ind - x0]; 2743 } 2744 res.w[1] = 0; 2745 res.w[0] = R64; 2746 } else if (ind <= 38) { 2747 // 19 <= ind <= 38 2748 P128.w[1] = res.w[1]; 2749 P128.w[0] = res.w[0]; 2750 round128_19_38 (ind, x0, P128, &res, &incr_exp, 2751 &is_midpoint_lt_even, &is_midpoint_gt_even, 2752 &is_inexact_lt_midpoint, 2753 &is_inexact_gt_midpoint); 2754 if (incr_exp) { 2755 // R128 = 10^(ind-x0), 1 <= ind - x0 <= ind - 1, 1 <= ind - x0 <= 37 2756 if (ind - x0 <= 19) { // 1 <= ind - x0 <= 19 2757 res.w[0] = ten2k64[ind - x0]; 2758 // res.w[1] stays 0 2759 } else { // 20 <= ind - x0 <= 37 2760 res.w[0] = ten2k128[ind - x0 - 20].w[0]; 2761 res.w[1] = ten2k128[ind - x0 - 20].w[1]; 2762 } 2763 } 2764 } 2765 // avoid a double rounding error 2766 if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 2767 is_midpoint_lt_even) { // double rounding error upward 2768 // res = res - 1 2769 res.w[0]--; 2770 if (res.w[0] == 0xffffffffffffffffull) 2771 res.w[1]--; 2772 // Note: a double rounding error upward is not possible; for this 2773 // the result after the first rounding would have to be 99...95 2774 // (35 digits in all), possibly followed by a number of zeros; this 2775 // not possible in Cases (2)-(6) which may get here 2776 is_midpoint_lt_even = 0; 2777 is_inexact_lt_midpoint = 1; 2778 } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 2779 is_midpoint_gt_even) { // double rounding error downward 2780 // res = res + 1 2781 res.w[0]++; 2782 if (res.w[0] == 0) 2783 res.w[1]++; 2784 is_midpoint_gt_even = 0; 2785 is_inexact_gt_midpoint = 1; 2786 } else if (!is_midpoint_lt_even && !is_midpoint_gt_even && 2787 !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) { 2788 // if this second rounding was exact the result may still be 2789 // inexact because of the first rounding 2790 if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) { 2791 is_inexact_gt_midpoint = 1; 2792 } 2793 if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) { 2794 is_inexact_lt_midpoint = 1; 2795 } 2796 } else if (is_midpoint_gt_even && 2797 (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) { 2798 // pulled up to a midpoint 2799 is_inexact_lt_midpoint = 1; 2800 is_inexact_gt_midpoint = 0; 2801 is_midpoint_lt_even = 0; 2802 is_midpoint_gt_even = 0; 2803 } else if (is_midpoint_lt_even && 2804 (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) { 2805 // pulled down to a midpoint 2806 is_inexact_lt_midpoint = 0; 2807 is_inexact_gt_midpoint = 1; 2808 is_midpoint_lt_even = 0; 2809 is_midpoint_gt_even = 0; 2810 } else { 2811 ; 2812 } 2813 // adjust exponent 2814 e3 = e3 + x0; 2815 if (!is_midpoint_lt_even && !is_midpoint_gt_even && 2816 !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) { 2817 if (is_midpoint_lt_even0 || is_midpoint_gt_even0 || 2818 is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0) { 2819 is_inexact_lt_midpoint = 1; 2820 } 2821 } 2822 } else { 2823 ; // not underflow 2824 } 2825 // check for inexact result 2826 if (is_inexact_lt_midpoint || is_inexact_gt_midpoint || 2827 is_midpoint_lt_even || is_midpoint_gt_even) { 2828 // set the inexact flag 2829 *pfpsf |= INEXACT_EXCEPTION; 2830 if (is_tiny) 2831 *pfpsf |= UNDERFLOW_EXCEPTION; 2832 } 2833 // now check for significand = 10^34 (may have resulted from going 2834 // back to case2_repeat) 2835 if (res.w[1] == 0x0001ed09bead87c0ull && 2836 res.w[0] == 0x378d8e6400000000ull) { // if res = 10^34 2837 res.w[1] = 0x0000314dc6448d93ull; // res = 10^33 2838 res.w[0] = 0x38c15b0a00000000ull; 2839 e3 = e3 + 1; 2840 } 2841 res.w[1] = z_sign | ((UINT64) (e3 + 6176) << 49) | res.w[1]; 2842 // check for overflow 2843 if (rnd_mode == ROUNDING_TO_NEAREST && e3 > expmax) { 2844 res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf 2845 res.w[0] = 0x0000000000000000ull; 2846 *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION); 2847 } 2848 if (rnd_mode != ROUNDING_TO_NEAREST) { 2849 rounding_correction (rnd_mode, 2850 is_inexact_lt_midpoint, 2851 is_inexact_gt_midpoint, 2852 is_midpoint_lt_even, is_midpoint_gt_even, 2853 e3, &res, pfpsf); 2854 } 2855 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 2856 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 2857 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 2858 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 2859 BID_SWAP128 (res); 2860 BID_RETURN (res) 2861 2862 } else { 2863 2864 // we get here only if delta <= 1 in Cases (2), (3), (4), (5), or (6) and 2865 // the signs of x*y and z are opposite; in these cases massive 2866 // cancellation can occur, so it is better to scale either C3 or C4 and 2867 // to perform the subtraction before rounding; rounding is performed 2868 // next, depending on the number of decimal digits in the result and on 2869 // the exponent value 2870 // Note: overlow is not possible in this case 2871 // this is similar to Cases (15), (16), and (17) 2872 2873 if (delta + q4 < q3) { // from Case (6) 2874 // Case (6) with 0<= delta <= 1 is similar to Cases (15), (16), and 2875 // (17) if we swap (C3, C4), (q3, q4), (e3, e4), (z_sign, p_sign) 2876 // and call add_and_round; delta stays positive 2877 // C4.w[3] = 0 and C4.w[2] = 0, so swap just the low part of C4 with C3 2878 P128.w[1] = C3.w[1]; 2879 P128.w[0] = C3.w[0]; 2880 C3.w[1] = C4.w[1]; 2881 C3.w[0] = C4.w[0]; 2882 C4.w[1] = P128.w[1]; 2883 C4.w[0] = P128.w[0]; 2884 ind = q3; 2885 q3 = q4; 2886 q4 = ind; 2887 ind = e3; 2888 e3 = e4; 2889 e4 = ind; 2890 tmp_sign = z_sign; 2891 z_sign = p_sign; 2892 p_sign = tmp_sign; 2893 } else { // from Cases (2), (3), (4), (5) 2894 // In Cases (2), (3), (4), (5) with 0 <= delta <= 1 C3 has to be 2895 // scaled up by q4 + delta - q3; this is the same as in Cases (15), 2896 // (16), and (17) if we just change the sign of delta 2897 delta = -delta; 2898 } 2899 add_and_round (q3, q4, e4, delta, p34, z_sign, p_sign, C3, C4, 2900 rnd_mode, &is_midpoint_lt_even, 2901 &is_midpoint_gt_even, &is_inexact_lt_midpoint, 2902 &is_inexact_gt_midpoint, pfpsf, &res); 2903 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 2904 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 2905 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 2906 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 2907 BID_SWAP128 (res); 2908 BID_RETURN (res) 2909 2910 } 2911 2912 } else { // if delta < 0 2913 2914 delta = -delta; 2915 2916 if (p34 < q4 && q4 <= delta) { // Case (7) 2917 2918 // truncate C4 to p34 digits into res 2919 // x = q4-p34, 1 <= x <= 34 because 35 <= q4 <= 68 2920 x0 = q4 - p34; 2921 if (q4 <= 38) { 2922 P128.w[1] = C4.w[1]; 2923 P128.w[0] = C4.w[0]; 2924 round128_19_38 (q4, x0, P128, &res, &incr_exp, 2925 &is_midpoint_lt_even, &is_midpoint_gt_even, 2926 &is_inexact_lt_midpoint, 2927 &is_inexact_gt_midpoint); 2928 } else if (q4 <= 57) { // 35 <= q4 <= 57 2929 P192.w[2] = C4.w[2]; 2930 P192.w[1] = C4.w[1]; 2931 P192.w[0] = C4.w[0]; 2932 round192_39_57 (q4, x0, P192, &R192, &incr_exp, 2933 &is_midpoint_lt_even, &is_midpoint_gt_even, 2934 &is_inexact_lt_midpoint, 2935 &is_inexact_gt_midpoint); 2936 res.w[0] = R192.w[0]; 2937 res.w[1] = R192.w[1]; 2938 } else { // if (q4 <= 68) 2939 round256_58_76 (q4, x0, C4, &R256, &incr_exp, 2940 &is_midpoint_lt_even, &is_midpoint_gt_even, 2941 &is_inexact_lt_midpoint, 2942 &is_inexact_gt_midpoint); 2943 res.w[0] = R256.w[0]; 2944 res.w[1] = R256.w[1]; 2945 } 2946 e4 = e4 + x0; 2947 if (incr_exp) { 2948 e4 = e4 + 1; 2949 } 2950 if (!is_midpoint_lt_even && !is_midpoint_gt_even && 2951 !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) { 2952 // if C4 rounded to p34 digits is exact then the result is inexact, 2953 // in a way that depends on the signs of x * y and z 2954 if (p_sign == z_sign) { 2955 is_inexact_lt_midpoint = 1; 2956 } else { // if (p_sign != z_sign) 2957 if (res.w[1] != 0x0000314dc6448d93ull || 2958 res.w[0] != 0x38c15b0a00000000ull) { // res != 10^33 2959 is_inexact_gt_midpoint = 1; 2960 } else { // res = 10^33 and exact is a special case 2961 // if C3 < 1/2 ulp then res = 10^33 and is_inexact_gt_midpoint = 1 2962 // if C3 = 1/2 ulp then res = 10^33 and is_midpoint_lt_even = 1 2963 // if C3 > 1/2 ulp then res = 10^34-1 and is_inexact_lt_midpoint = 1 2964 // Note: ulp is really ulp/10 (after borrow which propagates to msd) 2965 if (delta > p34 + 1) { // C3 < 1/2 2966 // res = 10^33, unchanged 2967 is_inexact_gt_midpoint = 1; 2968 } else { // if (delta == p34 + 1) 2969 if (q3 <= 19) { 2970 if (C3.w[0] < midpoint64[q3 - 1]) { // C3 < 1/2 ulp 2971 // res = 10^33, unchanged 2972 is_inexact_gt_midpoint = 1; 2973 } else if (C3.w[0] == midpoint64[q3 - 1]) { // C3 = 1/2 ulp 2974 // res = 10^33, unchanged 2975 is_midpoint_lt_even = 1; 2976 } else { // if (C3.w[0] > midpoint64[q3-1]), C3 > 1/2 ulp 2977 res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1 2978 res.w[0] = 0x378d8e63ffffffffull; 2979 e4 = e4 - 1; 2980 is_inexact_lt_midpoint = 1; 2981 } 2982 } else { // if (20 <= q3 <=34) 2983 if (C3.w[1] < midpoint128[q3 - 20].w[1] || 2984 (C3.w[1] == midpoint128[q3 - 20].w[1] && 2985 C3.w[0] < midpoint128[q3 - 20].w[0])) { // C3 < 1/2 ulp 2986 // res = 10^33, unchanged 2987 is_inexact_gt_midpoint = 1; 2988 } else if (C3.w[1] == midpoint128[q3 - 20].w[1] && 2989 C3.w[0] == midpoint128[q3 - 20].w[0]) { // C3 = 1/2 ulp 2990 // res = 10^33, unchanged 2991 is_midpoint_lt_even = 1; 2992 } else { // if (C3 > midpoint128[q3-20]), C3 > 1/2 ulp 2993 res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1 2994 res.w[0] = 0x378d8e63ffffffffull; 2995 e4 = e4 - 1; 2996 is_inexact_lt_midpoint = 1; 2997 } 2998 } 2999 } 3000 } 3001 } 3002 } else if (is_midpoint_lt_even) { 3003 if (z_sign != p_sign) { 3004 // needs correction: res = res - 1 3005 res.w[0] = res.w[0] - 1; 3006 if (res.w[0] == 0xffffffffffffffffull) 3007 res.w[1]--; 3008 // if it is (10^33-1)*10^e4 then the corect result is 3009 // (10^34-1)*10(e4-1) 3010 if (res.w[1] == 0x0000314dc6448d93ull && 3011 res.w[0] == 0x38c15b09ffffffffull) { 3012 res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1 3013 res.w[0] = 0x378d8e63ffffffffull; 3014 e4 = e4 - 1; 3015 } 3016 is_midpoint_lt_even = 0; 3017 is_inexact_lt_midpoint = 1; 3018 } else { // if (z_sign == p_sign) 3019 is_midpoint_lt_even = 0; 3020 is_inexact_gt_midpoint = 1; 3021 } 3022 } else if (is_midpoint_gt_even) { 3023 if (z_sign == p_sign) { 3024 // needs correction: res = res + 1 (cannot cross in the next binade) 3025 res.w[0] = res.w[0] + 1; 3026 if (res.w[0] == 0x0000000000000000ull) 3027 res.w[1]++; 3028 is_midpoint_gt_even = 0; 3029 is_inexact_gt_midpoint = 1; 3030 } else { // if (z_sign != p_sign) 3031 is_midpoint_gt_even = 0; 3032 is_inexact_lt_midpoint = 1; 3033 } 3034 } else { 3035 ; // the rounded result is already correct 3036 } 3037 // check for overflow 3038 if (rnd_mode == ROUNDING_TO_NEAREST && e4 > expmax) { 3039 res.w[1] = p_sign | 0x7800000000000000ull; 3040 res.w[0] = 0x0000000000000000ull; 3041 *pfpsf |= (OVERFLOW_EXCEPTION | INEXACT_EXCEPTION); 3042 } else { // no overflow or not RN 3043 p_exp = ((UINT64) (e4 + 6176) << 49); 3044 res.w[1] = p_sign | (p_exp & MASK_EXP) | res.w[1]; 3045 } 3046 if (rnd_mode != ROUNDING_TO_NEAREST) { 3047 rounding_correction (rnd_mode, 3048 is_inexact_lt_midpoint, 3049 is_inexact_gt_midpoint, 3050 is_midpoint_lt_even, is_midpoint_gt_even, 3051 e4, &res, pfpsf); 3052 } 3053 if (is_inexact_lt_midpoint || is_inexact_gt_midpoint || 3054 is_midpoint_lt_even || is_midpoint_gt_even) { 3055 // set the inexact flag 3056 *pfpsf |= INEXACT_EXCEPTION; 3057 } 3058 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 3059 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 3060 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 3061 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 3062 BID_SWAP128 (res); 3063 BID_RETURN (res) 3064 3065 } else if ((q4 <= p34 && p34 <= delta) || // Case (8) 3066 (q4 <= delta && delta < p34 && p34 < delta + q3) || // Case (9) 3067 (q4 <= delta && delta + q3 <= p34) || // Case (10) 3068 (delta < q4 && q4 <= p34 && p34 < delta + q3) || // Case (13) 3069 (delta < q4 && q4 <= delta + q3 && delta + q3 <= p34) || // Case (14) 3070 (delta + q3 < q4 && q4 <= p34)) { // Case (18) 3071 3072 // Case (8) is similar to Case (1), with C3 and C4 swapped 3073 // Case (9) is similar to Case (2), with C3 and C4 swapped 3074 // Case (10) is similar to Case (3), with C3 and C4 swapped 3075 // Case (13) is similar to Case (4), with C3 and C4 swapped 3076 // Case (14) is similar to Case (5), with C3 and C4 swapped 3077 // Case (18) is similar to Case (6), with C3 and C4 swapped 3078 3079 // swap (C3, C4), (q3, q4), (e3, 34), (z_sign, p_sign), (z_exp, p_exp) 3080 // and go back to delta_ge_zero 3081 // C4.w[3] = 0 and C4.w[2] = 0, so swap just the low part of C4 with C3 3082 P128.w[1] = C3.w[1]; 3083 P128.w[0] = C3.w[0]; 3084 C3.w[1] = C4.w[1]; 3085 C3.w[0] = C4.w[0]; 3086 C4.w[1] = P128.w[1]; 3087 C4.w[0] = P128.w[0]; 3088 ind = q3; 3089 q3 = q4; 3090 q4 = ind; 3091 ind = e3; 3092 e3 = e4; 3093 e4 = ind; 3094 tmp_sign = z_sign; 3095 z_sign = p_sign; 3096 p_sign = tmp_sign; 3097 tmp.ui64 = z_exp; 3098 z_exp = p_exp; 3099 p_exp = tmp.ui64; 3100 goto delta_ge_zero; 3101 3102 } else if ((p34 <= delta && delta < q4 && q4 < delta + q3) || // Case (11) 3103 (delta < p34 && p34 < q4 && q4 < delta + q3)) { // Case (12) 3104 3105 // round C3 to nearest to q3 - x0 digits, where x0 = e4 - e3, 3106 // 1 <= x0 <= q3 - 1 <= p34 - 1 3107 x0 = e4 - e3; // or x0 = delta + q3 - q4 3108 if (q3 <= 18) { // 2 <= q3 <= 18 3109 round64_2_18 (q3, x0, C3.w[0], &R64, &incr_exp, 3110 &is_midpoint_lt_even, &is_midpoint_gt_even, 3111 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint); 3112 // C3.w[1] = 0; 3113 C3.w[0] = R64; 3114 } else if (q3 <= 38) { 3115 round128_19_38 (q3, x0, C3, &R128, &incr_exp, 3116 &is_midpoint_lt_even, &is_midpoint_gt_even, 3117 &is_inexact_lt_midpoint, 3118 &is_inexact_gt_midpoint); 3119 C3.w[1] = R128.w[1]; 3120 C3.w[0] = R128.w[0]; 3121 } 3122 // the rounded result has q3 - x0 digits 3123 // we want the exponent to be e4, so if incr_exp = 1 then 3124 // multiply the rounded result by 10 - it will still fit in 113 bits 3125 if (incr_exp) { 3126 // 64 x 128 -> 128 3127 P128.w[1] = C3.w[1]; 3128 P128.w[0] = C3.w[0]; 3129 __mul_64x128_to_128 (C3, ten2k64[1], P128); 3130 } 3131 e3 = e3 + x0; // this is e4 3132 // now add/subtract the 256-bit C4 and the new (and shorter) 128-bit C3; 3133 // the result will have the sign of x * y; the exponent is e4 3134 R256.w[3] = 0; 3135 R256.w[2] = 0; 3136 R256.w[1] = C3.w[1]; 3137 R256.w[0] = C3.w[0]; 3138 if (p_sign == z_sign) { // R256 = C4 + R256 3139 add256 (C4, R256, &R256); 3140 } else { // if (p_sign != z_sign) { // R256 = C4 - R256 3141 sub256 (C4, R256, &R256); // the result cannot be pure zero 3142 // because the result has opposite sign to that of R256 which was 3143 // rounded, need to change the rounding indicators 3144 lsb = C4.w[0] & 0x01; 3145 if (is_inexact_lt_midpoint) { 3146 is_inexact_lt_midpoint = 0; 3147 is_inexact_gt_midpoint = 1; 3148 } else if (is_inexact_gt_midpoint) { 3149 is_inexact_gt_midpoint = 0; 3150 is_inexact_lt_midpoint = 1; 3151 } else if (lsb == 0) { 3152 if (is_midpoint_lt_even) { 3153 is_midpoint_lt_even = 0; 3154 is_midpoint_gt_even = 1; 3155 } else if (is_midpoint_gt_even) { 3156 is_midpoint_gt_even = 0; 3157 is_midpoint_lt_even = 1; 3158 } else { 3159 ; 3160 } 3161 } else if (lsb == 1) { 3162 if (is_midpoint_lt_even) { 3163 // res = res + 1 3164 R256.w[0]++; 3165 if (R256.w[0] == 0x0) { 3166 R256.w[1]++; 3167 if (R256.w[1] == 0x0) { 3168 R256.w[2]++; 3169 if (R256.w[2] == 0x0) { 3170 R256.w[3]++; 3171 } 3172 } 3173 } 3174 // no check for rounding overflow - R256 was a difference 3175 } else if (is_midpoint_gt_even) { 3176 // res = res - 1 3177 R256.w[0]--; 3178 if (R256.w[0] == 0xffffffffffffffffull) { 3179 R256.w[1]--; 3180 if (R256.w[1] == 0xffffffffffffffffull) { 3181 R256.w[2]--; 3182 if (R256.w[2] == 0xffffffffffffffffull) { 3183 R256.w[3]--; 3184 } 3185 } 3186 } 3187 } else { 3188 ; 3189 } 3190 } else { 3191 ; 3192 } 3193 } 3194 // determine the number of decimal digits in R256 3195 ind = nr_digits256 (R256); // ind >= p34 3196 // if R256 is sum, then ind > p34; if R256 is a difference, then 3197 // ind >= p34; this means that we can calculate the result rounded to 3198 // the destination precision, with unbounded exponent, starting from R256 3199 // and using the indicators from the rounding of C3 to avoid a double 3200 // rounding error 3201 3202 if (ind < p34) { 3203 ; 3204 } else if (ind == p34) { 3205 // the result rounded to the destination precision with 3206 // unbounded exponent 3207 // is (-1)^p_sign * R256 * 10^e4 3208 res.w[1] = R256.w[1]; 3209 res.w[0] = R256.w[0]; 3210 } else { // if (ind > p34) 3211 // if more than P digits, round to nearest to P digits 3212 // round R256 to p34 digits 3213 x0 = ind - p34; // 1 <= x0 <= 34 as 35 <= ind <= 68 3214 // save C3 rounding indicators to help avoid double rounding error 3215 is_inexact_lt_midpoint0 = is_inexact_lt_midpoint; 3216 is_inexact_gt_midpoint0 = is_inexact_gt_midpoint; 3217 is_midpoint_lt_even0 = is_midpoint_lt_even; 3218 is_midpoint_gt_even0 = is_midpoint_gt_even; 3219 // initialize rounding indicators 3220 is_inexact_lt_midpoint = 0; 3221 is_inexact_gt_midpoint = 0; 3222 is_midpoint_lt_even = 0; 3223 is_midpoint_gt_even = 0; 3224 // round to p34 digits; the result fits in 113 bits 3225 if (ind <= 38) { 3226 P128.w[1] = R256.w[1]; 3227 P128.w[0] = R256.w[0]; 3228 round128_19_38 (ind, x0, P128, &R128, &incr_exp, 3229 &is_midpoint_lt_even, &is_midpoint_gt_even, 3230 &is_inexact_lt_midpoint, 3231 &is_inexact_gt_midpoint); 3232 } else if (ind <= 57) { 3233 P192.w[2] = R256.w[2]; 3234 P192.w[1] = R256.w[1]; 3235 P192.w[0] = R256.w[0]; 3236 round192_39_57 (ind, x0, P192, &R192, &incr_exp, 3237 &is_midpoint_lt_even, &is_midpoint_gt_even, 3238 &is_inexact_lt_midpoint, 3239 &is_inexact_gt_midpoint); 3240 R128.w[1] = R192.w[1]; 3241 R128.w[0] = R192.w[0]; 3242 } else { // if (ind <= 68) 3243 round256_58_76 (ind, x0, R256, &R256, &incr_exp, 3244 &is_midpoint_lt_even, &is_midpoint_gt_even, 3245 &is_inexact_lt_midpoint, 3246 &is_inexact_gt_midpoint); 3247 R128.w[1] = R256.w[1]; 3248 R128.w[0] = R256.w[0]; 3249 } 3250 // the rounded result has p34 = 34 digits 3251 e4 = e4 + x0 + incr_exp; 3252 3253 res.w[1] = R128.w[1]; 3254 res.w[0] = R128.w[0]; 3255 3256 // avoid a double rounding error 3257 if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 3258 is_midpoint_lt_even) { // double rounding error upward 3259 // res = res - 1 3260 res.w[0]--; 3261 if (res.w[0] == 0xffffffffffffffffull) 3262 res.w[1]--; 3263 is_midpoint_lt_even = 0; 3264 is_inexact_lt_midpoint = 1; 3265 // Note: a double rounding error upward is not possible; for this 3266 // the result after the first rounding would have to be 99...95 3267 // (35 digits in all), possibly followed by a number of zeros; this 3268 // not possible in Cases (2)-(6) or (15)-(17) which may get here 3269 // if this is 10^33 - 1 make it 10^34 - 1 and decrement exponent 3270 if (res.w[1] == 0x0000314dc6448d93ull && 3271 res.w[0] == 0x38c15b09ffffffffull) { // 10^33 - 1 3272 res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1 3273 res.w[0] = 0x378d8e63ffffffffull; 3274 e4--; 3275 } 3276 } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 3277 is_midpoint_gt_even) { // double rounding error downward 3278 // res = res + 1 3279 res.w[0]++; 3280 if (res.w[0] == 0) 3281 res.w[1]++; 3282 is_midpoint_gt_even = 0; 3283 is_inexact_gt_midpoint = 1; 3284 } else if (!is_midpoint_lt_even && !is_midpoint_gt_even && 3285 !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) { 3286 // if this second rounding was exact the result may still be 3287 // inexact because of the first rounding 3288 if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) { 3289 is_inexact_gt_midpoint = 1; 3290 } 3291 if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) { 3292 is_inexact_lt_midpoint = 1; 3293 } 3294 } else if (is_midpoint_gt_even && 3295 (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) { 3296 // pulled up to a midpoint 3297 is_inexact_lt_midpoint = 1; 3298 is_inexact_gt_midpoint = 0; 3299 is_midpoint_lt_even = 0; 3300 is_midpoint_gt_even = 0; 3301 } else if (is_midpoint_lt_even && 3302 (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) { 3303 // pulled down to a midpoint 3304 is_inexact_lt_midpoint = 0; 3305 is_inexact_gt_midpoint = 1; 3306 is_midpoint_lt_even = 0; 3307 is_midpoint_gt_even = 0; 3308 } else { 3309 ; 3310 } 3311 } 3312 3313 // determine tininess 3314 if (rnd_mode == ROUNDING_TO_NEAREST) { 3315 if (e4 < expmin) { 3316 is_tiny = 1; // for other rounding modes apply correction 3317 } 3318 } else { 3319 // for RM, RP, RZ, RA apply correction in order to determine tininess 3320 // but do not save the result; apply the correction to 3321 // (-1)^p_sign * res * 10^0 3322 P128.w[1] = p_sign | 0x3040000000000000ull | res.w[1]; 3323 P128.w[0] = res.w[0]; 3324 rounding_correction (rnd_mode, 3325 is_inexact_lt_midpoint, 3326 is_inexact_gt_midpoint, 3327 is_midpoint_lt_even, is_midpoint_gt_even, 3328 0, &P128, pfpsf); 3329 scale = ((P128.w[1] & MASK_EXP) >> 49) - 6176; // -1, 0, or +1 3330 // the number of digits in the significand is p34 = 34 3331 if (e4 + scale < expmin) { 3332 is_tiny = 1; 3333 } 3334 } 3335 3336 // the result rounded to the destination precision with unbounded exponent 3337 // is (-1)^p_sign * res * 10^e4 3338 res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | res.w[1]; // RN 3339 // res.w[0] unchanged; 3340 // Note: res is correct only if expmin <= e4 <= expmax 3341 ind = p34; // the number of decimal digits in the signifcand of res 3342 3343 // at this point we have the result rounded with unbounded exponent in 3344 // res and we know its tininess: 3345 // res = (-1)^p_sign * significand * 10^e4, 3346 // where q (significand) = ind = p34 3347 // Note: res is correct only if expmin <= e4 <= expmax 3348 3349 // check for overflow if RN 3350 if (rnd_mode == ROUNDING_TO_NEAREST 3351 && (ind + e4) > (p34 + expmax)) { 3352 res.w[1] = p_sign | 0x7800000000000000ull; 3353 res.w[0] = 0x0000000000000000ull; 3354 *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION); 3355 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 3356 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 3357 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 3358 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 3359 BID_SWAP128 (res); 3360 BID_RETURN (res) 3361 } // else not overflow or not RN, so continue 3362 3363 // from this point on this is similar to the last part of the computation 3364 // for Cases (15), (16), (17) 3365 3366 // if (e4 >= expmin) we have the result rounded with bounded exponent 3367 if (e4 < expmin) { 3368 x0 = expmin - e4; // x0 >= 1; the number of digits to chop off of res 3369 // where the result rounded [at most] once is 3370 // (-1)^p_sign * significand_res * 10^e4 3371 3372 // avoid double rounding error 3373 is_inexact_lt_midpoint0 = is_inexact_lt_midpoint; 3374 is_inexact_gt_midpoint0 = is_inexact_gt_midpoint; 3375 is_midpoint_lt_even0 = is_midpoint_lt_even; 3376 is_midpoint_gt_even0 = is_midpoint_gt_even; 3377 is_inexact_lt_midpoint = 0; 3378 is_inexact_gt_midpoint = 0; 3379 is_midpoint_lt_even = 0; 3380 is_midpoint_gt_even = 0; 3381 3382 if (x0 > ind) { 3383 // nothing is left of res when moving the decimal point left x0 digits 3384 is_inexact_lt_midpoint = 1; 3385 res.w[1] = p_sign | 0x0000000000000000ull; 3386 res.w[0] = 0x0000000000000000ull; 3387 e4 = expmin; 3388 } else if (x0 == ind) { // 1 <= x0 = ind <= p34 = 34 3389 // this is <, =, or > 1/2 ulp 3390 // compare the ind-digit value in the significand of res with 3391 // 1/2 ulp = 5*10^(ind-1), i.e. determine whether it is 3392 // less than, equal to, or greater than 1/2 ulp (significand of res) 3393 R128.w[1] = res.w[1] & MASK_COEFF; 3394 R128.w[0] = res.w[0]; 3395 if (ind <= 19) { 3396 if (R128.w[0] < midpoint64[ind - 1]) { // < 1/2 ulp 3397 lt_half_ulp = 1; 3398 is_inexact_lt_midpoint = 1; 3399 } else if (R128.w[0] == midpoint64[ind - 1]) { // = 1/2 ulp 3400 eq_half_ulp = 1; 3401 is_midpoint_gt_even = 1; 3402 } else { // > 1/2 ulp 3403 gt_half_ulp = 1; 3404 is_inexact_gt_midpoint = 1; 3405 } 3406 } else { // if (ind <= 38) 3407 if (R128.w[1] < midpoint128[ind - 20].w[1] || 3408 (R128.w[1] == midpoint128[ind - 20].w[1] && 3409 R128.w[0] < midpoint128[ind - 20].w[0])) { // < 1/2 ulp 3410 lt_half_ulp = 1; 3411 is_inexact_lt_midpoint = 1; 3412 } else if (R128.w[1] == midpoint128[ind - 20].w[1] && 3413 R128.w[0] == midpoint128[ind - 20].w[0]) { // = 1/2 ulp 3414 eq_half_ulp = 1; 3415 is_midpoint_gt_even = 1; 3416 } else { // > 1/2 ulp 3417 gt_half_ulp = 1; 3418 is_inexact_gt_midpoint = 1; 3419 } 3420 } 3421 if (lt_half_ulp || eq_half_ulp) { 3422 // res = +0.0 * 10^expmin 3423 res.w[1] = 0x0000000000000000ull; 3424 res.w[0] = 0x0000000000000000ull; 3425 } else { // if (gt_half_ulp) 3426 // res = +1 * 10^expmin 3427 res.w[1] = 0x0000000000000000ull; 3428 res.w[0] = 0x0000000000000001ull; 3429 } 3430 res.w[1] = p_sign | res.w[1]; 3431 e4 = expmin; 3432 } else { // if (1 <= x0 <= ind - 1 <= 33) 3433 // round the ind-digit result to ind - x0 digits 3434 3435 if (ind <= 18) { // 2 <= ind <= 18 3436 round64_2_18 (ind, x0, res.w[0], &R64, &incr_exp, 3437 &is_midpoint_lt_even, &is_midpoint_gt_even, 3438 &is_inexact_lt_midpoint, 3439 &is_inexact_gt_midpoint); 3440 res.w[1] = 0x0; 3441 res.w[0] = R64; 3442 } else if (ind <= 38) { 3443 P128.w[1] = res.w[1] & MASK_COEFF; 3444 P128.w[0] = res.w[0]; 3445 round128_19_38 (ind, x0, P128, &res, &incr_exp, 3446 &is_midpoint_lt_even, &is_midpoint_gt_even, 3447 &is_inexact_lt_midpoint, 3448 &is_inexact_gt_midpoint); 3449 } 3450 e4 = e4 + x0; // expmin 3451 // we want the exponent to be expmin, so if incr_exp = 1 then 3452 // multiply the rounded result by 10 - it will still fit in 113 bits 3453 if (incr_exp) { 3454 // 64 x 128 -> 128 3455 P128.w[1] = res.w[1] & MASK_COEFF; 3456 P128.w[0] = res.w[0]; 3457 __mul_64x128_to_128 (res, ten2k64[1], P128); 3458 } 3459 res.w[1] = 3460 p_sign | ((UINT64) (e4 + 6176) << 49) | (res. 3461 w[1] & MASK_COEFF); 3462 // avoid a double rounding error 3463 if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 3464 is_midpoint_lt_even) { // double rounding error upward 3465 // res = res - 1 3466 res.w[0]--; 3467 if (res.w[0] == 0xffffffffffffffffull) 3468 res.w[1]--; 3469 // Note: a double rounding error upward is not possible; for this 3470 // the result after the first rounding would have to be 99...95 3471 // (35 digits in all), possibly followed by a number of zeros; this 3472 // not possible in this underflow case 3473 is_midpoint_lt_even = 0; 3474 is_inexact_lt_midpoint = 1; 3475 } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 3476 is_midpoint_gt_even) { // double rounding error downward 3477 // res = res + 1 3478 res.w[0]++; 3479 if (res.w[0] == 0) 3480 res.w[1]++; 3481 is_midpoint_gt_even = 0; 3482 is_inexact_gt_midpoint = 1; 3483 } else if (!is_midpoint_lt_even && !is_midpoint_gt_even && 3484 !is_inexact_lt_midpoint 3485 && !is_inexact_gt_midpoint) { 3486 // if this second rounding was exact the result may still be 3487 // inexact because of the first rounding 3488 if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) { 3489 is_inexact_gt_midpoint = 1; 3490 } 3491 if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) { 3492 is_inexact_lt_midpoint = 1; 3493 } 3494 } else if (is_midpoint_gt_even && 3495 (is_inexact_gt_midpoint0 3496 || is_midpoint_lt_even0)) { 3497 // pulled up to a midpoint 3498 is_inexact_lt_midpoint = 1; 3499 is_inexact_gt_midpoint = 0; 3500 is_midpoint_lt_even = 0; 3501 is_midpoint_gt_even = 0; 3502 } else if (is_midpoint_lt_even && 3503 (is_inexact_lt_midpoint0 3504 || is_midpoint_gt_even0)) { 3505 // pulled down to a midpoint 3506 is_inexact_lt_midpoint = 0; 3507 is_inexact_gt_midpoint = 1; 3508 is_midpoint_lt_even = 0; 3509 is_midpoint_gt_even = 0; 3510 } else { 3511 ; 3512 } 3513 } 3514 } 3515 // res contains the correct result 3516 // apply correction if not rounding to nearest 3517 if (rnd_mode != ROUNDING_TO_NEAREST) { 3518 rounding_correction (rnd_mode, 3519 is_inexact_lt_midpoint, 3520 is_inexact_gt_midpoint, 3521 is_midpoint_lt_even, is_midpoint_gt_even, 3522 e4, &res, pfpsf); 3523 } 3524 if (is_midpoint_lt_even || is_midpoint_gt_even || 3525 is_inexact_lt_midpoint || is_inexact_gt_midpoint) { 3526 // set the inexact flag 3527 *pfpsf |= INEXACT_EXCEPTION; 3528 if (is_tiny) 3529 *pfpsf |= UNDERFLOW_EXCEPTION; 3530 } 3531 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 3532 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 3533 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 3534 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 3535 BID_SWAP128 (res); 3536 BID_RETURN (res) 3537 3538 } else if ((p34 <= delta && delta + q3 <= q4) || // Case (15) 3539 (delta < p34 && p34 < delta + q3 && delta + q3 <= q4) || //Case (16) 3540 (delta + q3 <= p34 && p34 < q4)) { // Case (17) 3541 3542 // calculate first the result rounded to the destination precision, with 3543 // unbounded exponent 3544 3545 add_and_round (q3, q4, e4, delta, p34, z_sign, p_sign, C3, C4, 3546 rnd_mode, &is_midpoint_lt_even, 3547 &is_midpoint_gt_even, &is_inexact_lt_midpoint, 3548 &is_inexact_gt_midpoint, pfpsf, &res); 3549 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 3550 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 3551 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 3552 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 3553 BID_SWAP128 (res); 3554 BID_RETURN (res) 3555 3556 } else { 3557 ; 3558 } 3559 3560 } // end if delta < 0 3561 3562 *ptr_is_midpoint_lt_even = is_midpoint_lt_even; 3563 *ptr_is_midpoint_gt_even = is_midpoint_gt_even; 3564 *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint; 3565 *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint; 3566 BID_SWAP128 (res); 3567 BID_RETURN (res) 3568 3569 } 3570 3571 3572 #if DECIMAL_CALL_BY_REFERENCE 3573 void 3574 bid128_fma (UINT128 * pres, UINT128 * px, UINT128 * py, UINT128 * pz 3575 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3576 _EXC_INFO_PARAM) { 3577 UINT128 x = *px, y = *py, z = *pz; 3578 #if !DECIMAL_GLOBAL_ROUNDING 3579 unsigned int rnd_mode = *prnd_mode; 3580 #endif 3581 #else 3582 UINT128 3583 bid128_fma (UINT128 x, UINT128 y, UINT128 z 3584 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3585 _EXC_INFO_PARAM) { 3586 #endif 3587 int is_midpoint_lt_even, is_midpoint_gt_even, 3588 is_inexact_lt_midpoint, is_inexact_gt_midpoint; 3589 UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} }; 3590 3591 #if DECIMAL_CALL_BY_REFERENCE 3592 bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3593 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint, 3594 &res, &x, &y, &z 3595 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3596 _EXC_INFO_ARG); 3597 #else 3598 res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3599 &is_inexact_lt_midpoint, 3600 &is_inexact_gt_midpoint, x, y, 3601 z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3602 _EXC_INFO_ARG); 3603 #endif 3604 BID_RETURN (res); 3605 } 3606 3607 3608 #if DECIMAL_CALL_BY_REFERENCE 3609 void 3610 bid128ddd_fma (UINT128 * pres, UINT64 * px, UINT64 * py, UINT64 * pz 3611 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3612 _EXC_INFO_PARAM) { 3613 UINT64 x = *px, y = *py, z = *pz; 3614 #if !DECIMAL_GLOBAL_ROUNDING 3615 unsigned int rnd_mode = *prnd_mode; 3616 #endif 3617 #else 3618 UINT128 3619 bid128ddd_fma (UINT64 x, UINT64 y, UINT64 z 3620 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3621 _EXC_INFO_PARAM) { 3622 #endif 3623 int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0, 3624 is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0; 3625 UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} }; 3626 UINT128 x1, y1, z1; 3627 3628 #if DECIMAL_CALL_BY_REFERENCE 3629 bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3630 bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3631 bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3632 bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3633 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint, 3634 &res, &x1, &y1, &z1 3635 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3636 _EXC_INFO_ARG); 3637 #else 3638 x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3639 y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3640 z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3641 res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3642 &is_inexact_lt_midpoint, 3643 &is_inexact_gt_midpoint, x1, y1, 3644 z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3645 _EXC_INFO_ARG); 3646 #endif 3647 BID_RETURN (res); 3648 } 3649 3650 3651 #if DECIMAL_CALL_BY_REFERENCE 3652 void 3653 bid128ddq_fma (UINT128 * pres, UINT64 * px, UINT64 * py, UINT128 * pz 3654 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3655 _EXC_INFO_PARAM) { 3656 UINT64 x = *px, y = *py; 3657 UINT128 z = *pz; 3658 #if !DECIMAL_GLOBAL_ROUNDING 3659 unsigned int rnd_mode = *prnd_mode; 3660 #endif 3661 #else 3662 UINT128 3663 bid128ddq_fma (UINT64 x, UINT64 y, UINT128 z 3664 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3665 _EXC_INFO_PARAM) { 3666 #endif 3667 int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0, 3668 is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0; 3669 UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} }; 3670 UINT128 x1, y1; 3671 3672 #if DECIMAL_CALL_BY_REFERENCE 3673 bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3674 bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3675 bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3676 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint, 3677 &res, &x1, &y1, &z 3678 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3679 _EXC_INFO_ARG); 3680 #else 3681 x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3682 y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3683 res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3684 &is_inexact_lt_midpoint, 3685 &is_inexact_gt_midpoint, x1, y1, 3686 z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3687 _EXC_INFO_ARG); 3688 #endif 3689 BID_RETURN (res); 3690 } 3691 3692 3693 #if DECIMAL_CALL_BY_REFERENCE 3694 void 3695 bid128dqd_fma (UINT128 * pres, UINT64 * px, UINT128 * py, UINT64 * pz 3696 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3697 _EXC_INFO_PARAM) { 3698 UINT64 x = *px, z = *pz; 3699 #if !DECIMAL_GLOBAL_ROUNDING 3700 unsigned int rnd_mode = *prnd_mode; 3701 #endif 3702 #else 3703 UINT128 3704 bid128dqd_fma (UINT64 x, UINT128 y, UINT64 z 3705 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3706 _EXC_INFO_PARAM) { 3707 #endif 3708 int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0, 3709 is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0; 3710 UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} }; 3711 UINT128 x1, z1; 3712 3713 #if DECIMAL_CALL_BY_REFERENCE 3714 bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3715 bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3716 bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3717 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint, 3718 &res, &x1, py, &z1 3719 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3720 _EXC_INFO_ARG); 3721 #else 3722 x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3723 z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3724 res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3725 &is_inexact_lt_midpoint, 3726 &is_inexact_gt_midpoint, x1, y, 3727 z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3728 _EXC_INFO_ARG); 3729 #endif 3730 BID_RETURN (res); 3731 } 3732 3733 3734 #if DECIMAL_CALL_BY_REFERENCE 3735 void 3736 bid128dqq_fma (UINT128 * pres, UINT64 * px, UINT128 * py, UINT128 * pz 3737 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3738 _EXC_INFO_PARAM) { 3739 UINT64 x = *px; 3740 #if !DECIMAL_GLOBAL_ROUNDING 3741 unsigned int rnd_mode = *prnd_mode; 3742 #endif 3743 #else 3744 UINT128 3745 bid128dqq_fma (UINT64 x, UINT128 y, UINT128 z 3746 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3747 _EXC_INFO_PARAM) { 3748 #endif 3749 int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0, 3750 is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0; 3751 UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} }; 3752 UINT128 x1; 3753 3754 #if DECIMAL_CALL_BY_REFERENCE 3755 bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3756 bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3757 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint, 3758 &res, &x1, py, pz 3759 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3760 _EXC_INFO_ARG); 3761 #else 3762 x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3763 res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3764 &is_inexact_lt_midpoint, 3765 &is_inexact_gt_midpoint, x1, y, 3766 z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3767 _EXC_INFO_ARG); 3768 #endif 3769 BID_RETURN (res); 3770 } 3771 3772 3773 #if DECIMAL_CALL_BY_REFERENCE 3774 void 3775 bid128qdd_fma (UINT128 * pres, UINT128 * px, UINT64 * py, UINT64 * pz 3776 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3777 _EXC_INFO_PARAM) { 3778 UINT64 y = *py, z = *pz; 3779 #if !DECIMAL_GLOBAL_ROUNDING 3780 unsigned int rnd_mode = *prnd_mode; 3781 #endif 3782 #else 3783 UINT128 3784 bid128qdd_fma (UINT128 x, UINT64 y, UINT64 z 3785 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3786 _EXC_INFO_PARAM) { 3787 #endif 3788 int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0, 3789 is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0; 3790 UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} }; 3791 UINT128 y1, z1; 3792 3793 #if DECIMAL_CALL_BY_REFERENCE 3794 bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3795 bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3796 bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3797 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint, 3798 &res, px, &y1, &z1 3799 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3800 _EXC_INFO_ARG); 3801 #else 3802 y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3803 z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3804 res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3805 &is_inexact_lt_midpoint, 3806 &is_inexact_gt_midpoint, x, y1, 3807 z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3808 _EXC_INFO_ARG); 3809 #endif 3810 BID_RETURN (res); 3811 } 3812 3813 3814 #if DECIMAL_CALL_BY_REFERENCE 3815 void 3816 bid128qdq_fma (UINT128 * pres, UINT128 * px, UINT64 * py, UINT128 * pz 3817 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3818 _EXC_INFO_PARAM) { 3819 UINT64 y = *py; 3820 #if !DECIMAL_GLOBAL_ROUNDING 3821 unsigned int rnd_mode = *prnd_mode; 3822 #endif 3823 #else 3824 UINT128 3825 bid128qdq_fma (UINT128 x, UINT64 y, UINT128 z 3826 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3827 _EXC_INFO_PARAM) { 3828 #endif 3829 int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0, 3830 is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0; 3831 UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} }; 3832 UINT128 y1; 3833 3834 #if DECIMAL_CALL_BY_REFERENCE 3835 bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3836 bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3837 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint, 3838 &res, px, &y1, pz 3839 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3840 _EXC_INFO_ARG); 3841 #else 3842 y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3843 res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3844 &is_inexact_lt_midpoint, 3845 &is_inexact_gt_midpoint, x, y1, 3846 z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3847 _EXC_INFO_ARG); 3848 #endif 3849 BID_RETURN (res); 3850 } 3851 3852 3853 #if DECIMAL_CALL_BY_REFERENCE 3854 void 3855 bid128qqd_fma (UINT128 * pres, UINT128 * px, UINT128 * py, UINT64 * pz 3856 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3857 _EXC_INFO_PARAM) { 3858 UINT64 z = *pz; 3859 #if !DECIMAL_GLOBAL_ROUNDING 3860 unsigned int rnd_mode = *prnd_mode; 3861 #endif 3862 #else 3863 UINT128 3864 bid128qqd_fma (UINT128 x, UINT128 y, UINT64 z 3865 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3866 _EXC_INFO_PARAM) { 3867 #endif 3868 int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0, 3869 is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0; 3870 UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} }; 3871 UINT128 z1; 3872 3873 #if DECIMAL_CALL_BY_REFERENCE 3874 bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3875 bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3876 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint, 3877 &res, px, py, &z1 3878 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3879 _EXC_INFO_ARG); 3880 #else 3881 z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3882 res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even, 3883 &is_inexact_lt_midpoint, 3884 &is_inexact_gt_midpoint, x, y, 3885 z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3886 _EXC_INFO_ARG); 3887 #endif 3888 BID_RETURN (res); 3889 } 3890 3891 // Note: bid128qqq_fma is represented by bid128_fma 3892 3893 // Note: bid64ddd_fma is represented by bid64_fma 3894 3895 #if DECIMAL_CALL_BY_REFERENCE 3896 void 3897 bid64ddq_fma (UINT64 * pres, UINT64 * px, UINT64 * py, UINT128 * pz 3898 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3899 _EXC_INFO_PARAM) { 3900 UINT64 x = *px, y = *py; 3901 #if !DECIMAL_GLOBAL_ROUNDING 3902 unsigned int rnd_mode = *prnd_mode; 3903 #endif 3904 #else 3905 UINT64 3906 bid64ddq_fma (UINT64 x, UINT64 y, UINT128 z 3907 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3908 _EXC_INFO_PARAM) { 3909 #endif 3910 UINT64 res1 = 0xbaddbaddbaddbaddull; 3911 UINT128 x1, y1; 3912 3913 #if DECIMAL_CALL_BY_REFERENCE 3914 bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3915 bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3916 bid64qqq_fma (&res1, &x1, &y1, pz 3917 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3918 _EXC_INFO_ARG); 3919 #else 3920 x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3921 y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3922 res1 = bid64qqq_fma (x1, y1, z 3923 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3924 _EXC_INFO_ARG); 3925 #endif 3926 BID_RETURN (res1); 3927 } 3928 3929 3930 #if DECIMAL_CALL_BY_REFERENCE 3931 void 3932 bid64dqd_fma (UINT64 * pres, UINT64 * px, UINT128 * py, UINT64 * pz 3933 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3934 _EXC_INFO_PARAM) { 3935 UINT64 x = *px, z = *pz; 3936 #if !DECIMAL_GLOBAL_ROUNDING 3937 unsigned int rnd_mode = *prnd_mode; 3938 #endif 3939 #else 3940 UINT64 3941 bid64dqd_fma (UINT64 x, UINT128 y, UINT64 z 3942 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3943 _EXC_INFO_PARAM) { 3944 #endif 3945 UINT64 res1 = 0xbaddbaddbaddbaddull; 3946 UINT128 x1, z1; 3947 3948 #if DECIMAL_CALL_BY_REFERENCE 3949 bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3950 bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3951 bid64qqq_fma (&res1, &x1, py, &z1 3952 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3953 _EXC_INFO_ARG); 3954 #else 3955 x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3956 z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3957 res1 = bid64qqq_fma (x1, y, z1 3958 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3959 _EXC_INFO_ARG); 3960 #endif 3961 BID_RETURN (res1); 3962 } 3963 3964 3965 #if DECIMAL_CALL_BY_REFERENCE 3966 void 3967 bid64dqq_fma (UINT64 * pres, UINT64 * px, UINT128 * py, UINT128 * pz 3968 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3969 _EXC_INFO_PARAM) { 3970 UINT64 x = *px; 3971 #if !DECIMAL_GLOBAL_ROUNDING 3972 unsigned int rnd_mode = *prnd_mode; 3973 #endif 3974 #else 3975 UINT64 3976 bid64dqq_fma (UINT64 x, UINT128 y, UINT128 z 3977 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 3978 _EXC_INFO_PARAM) { 3979 #endif 3980 UINT64 res1 = 0xbaddbaddbaddbaddull; 3981 UINT128 x1; 3982 3983 #if DECIMAL_CALL_BY_REFERENCE 3984 bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3985 bid64qqq_fma (&res1, &x1, py, pz 3986 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3987 _EXC_INFO_ARG); 3988 #else 3989 x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 3990 res1 = bid64qqq_fma (x1, y, z 3991 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 3992 _EXC_INFO_ARG); 3993 #endif 3994 BID_RETURN (res1); 3995 } 3996 3997 3998 #if DECIMAL_CALL_BY_REFERENCE 3999 void 4000 bid64qdd_fma (UINT64 * pres, UINT128 * px, UINT64 * py, UINT64 * pz 4001 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 4002 _EXC_INFO_PARAM) { 4003 UINT64 y = *py, z = *pz; 4004 #if !DECIMAL_GLOBAL_ROUNDING 4005 unsigned int rnd_mode = *prnd_mode; 4006 #endif 4007 #else 4008 UINT64 4009 bid64qdd_fma (UINT128 x, UINT64 y, UINT64 z 4010 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 4011 _EXC_INFO_PARAM) { 4012 #endif 4013 UINT64 res1 = 0xbaddbaddbaddbaddull; 4014 UINT128 y1, z1; 4015 4016 #if DECIMAL_CALL_BY_REFERENCE 4017 bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 4018 bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 4019 bid64qqq_fma (&res1, px, &y1, &z1 4020 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 4021 _EXC_INFO_ARG); 4022 #else 4023 y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 4024 z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 4025 res1 = bid64qqq_fma (x, y1, z1 4026 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 4027 _EXC_INFO_ARG); 4028 #endif 4029 BID_RETURN (res1); 4030 } 4031 4032 4033 #if DECIMAL_CALL_BY_REFERENCE 4034 void 4035 bid64qdq_fma (UINT64 * pres, UINT128 * px, UINT64 * py, UINT128 * pz 4036 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 4037 _EXC_INFO_PARAM) { 4038 UINT64 y = *py; 4039 #if !DECIMAL_GLOBAL_ROUNDING 4040 unsigned int rnd_mode = *prnd_mode; 4041 #endif 4042 #else 4043 UINT64 4044 bid64qdq_fma (UINT128 x, UINT64 y, UINT128 z 4045 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 4046 _EXC_INFO_PARAM) { 4047 #endif 4048 UINT64 res1 = 0xbaddbaddbaddbaddull; 4049 UINT128 y1; 4050 4051 #if DECIMAL_CALL_BY_REFERENCE 4052 bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 4053 bid64qqq_fma (&res1, px, &y1, pz 4054 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 4055 _EXC_INFO_ARG); 4056 #else 4057 y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 4058 res1 = bid64qqq_fma (x, y1, z 4059 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 4060 _EXC_INFO_ARG); 4061 #endif 4062 BID_RETURN (res1); 4063 } 4064 4065 4066 #if DECIMAL_CALL_BY_REFERENCE 4067 void 4068 bid64qqd_fma (UINT64 * pres, UINT128 * px, UINT128 * py, UINT64 * pz 4069 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 4070 _EXC_INFO_PARAM) { 4071 UINT64 z = *pz; 4072 #if !DECIMAL_GLOBAL_ROUNDING 4073 unsigned int rnd_mode = *prnd_mode; 4074 #endif 4075 #else 4076 UINT64 4077 bid64qqd_fma (UINT128 x, UINT128 y, UINT64 z 4078 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 4079 _EXC_INFO_PARAM) { 4080 #endif 4081 UINT64 res1 = 0xbaddbaddbaddbaddull; 4082 UINT128 z1; 4083 4084 #if DECIMAL_CALL_BY_REFERENCE 4085 bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 4086 bid64qqq_fma (&res1, px, py, &z1 4087 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 4088 _EXC_INFO_ARG); 4089 #else 4090 z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 4091 res1 = bid64qqq_fma (x, y, z1 4092 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 4093 _EXC_INFO_ARG); 4094 #endif 4095 BID_RETURN (res1); 4096 } 4097 4098 4099 #if DECIMAL_CALL_BY_REFERENCE 4100 void 4101 bid64qqq_fma (UINT64 * pres, UINT128 * px, UINT128 * py, UINT128 * pz 4102 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 4103 _EXC_INFO_PARAM) { 4104 UINT128 x = *px, y = *py, z = *pz; 4105 #if !DECIMAL_GLOBAL_ROUNDING 4106 unsigned int rnd_mode = *prnd_mode; 4107 #endif 4108 #else 4109 UINT64 4110 bid64qqq_fma (UINT128 x, UINT128 y, UINT128 z 4111 _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 4112 _EXC_INFO_PARAM) { 4113 #endif 4114 int is_midpoint_lt_even0 = 0, is_midpoint_gt_even0 = 0, 4115 is_inexact_lt_midpoint0 = 0, is_inexact_gt_midpoint0 = 0; 4116 int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0, 4117 is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0; 4118 int incr_exp; 4119 UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} }; 4120 UINT128 res128 = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} }; 4121 UINT64 res1 = 0xbaddbaddbaddbaddull; 4122 unsigned int save_fpsf; // needed because of the call to bid128_ext_fma 4123 UINT64 sign; 4124 UINT64 exp; 4125 int unbexp; 4126 UINT128 C; 4127 BID_UI64DOUBLE tmp; 4128 int nr_bits; 4129 int q, x0; 4130 int scale; 4131 int lt_half_ulp = 0, eq_half_ulp = 0; 4132 4133 // Note: for rounding modes other than RN or RA, the result can be obtained 4134 // by rounding first to BID128 and then to BID64 4135 4136 save_fpsf = *pfpsf; // sticky bits - caller value must be preserved 4137 *pfpsf = 0; 4138 4139 #if DECIMAL_CALL_BY_REFERENCE 4140 bid128_ext_fma (&is_midpoint_lt_even0, &is_midpoint_gt_even0, 4141 &is_inexact_lt_midpoint0, &is_inexact_gt_midpoint0, 4142 &res, &x, &y, &z 4143 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 4144 _EXC_INFO_ARG); 4145 #else 4146 res = bid128_ext_fma (&is_midpoint_lt_even0, &is_midpoint_gt_even0, 4147 &is_inexact_lt_midpoint0, 4148 &is_inexact_gt_midpoint0, x, y, 4149 z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG 4150 _EXC_INFO_ARG); 4151 #endif 4152 4153 if ((rnd_mode == ROUNDING_DOWN) || (rnd_mode == ROUNDING_UP) || 4154 (rnd_mode == ROUNDING_TO_ZERO) || // no double rounding error is possible 4155 ((res.w[HIGH_128W] & MASK_NAN) == MASK_NAN) || //res=QNaN (cannot be SNaN) 4156 ((res.w[HIGH_128W] & MASK_ANY_INF) == MASK_INF)) { // result is infinity 4157 #if DECIMAL_CALL_BY_REFERENCE 4158 bid128_to_bid64 (&res1, &res _RND_MODE_ARG _EXC_FLAGS_ARG); 4159 #else 4160 res1 = bid128_to_bid64 (res _RND_MODE_ARG _EXC_FLAGS_ARG); 4161 #endif 4162 // determine the unbiased exponent of the result 4163 unbexp = ((res1 >> 53) & 0x3ff) - 398; 4164 4165 // if subnormal, res1 must have exp = -398 4166 // if tiny and inexact set underflow and inexact status flags 4167 if (!((res1 & MASK_NAN) == MASK_NAN) && // res1 not NaN 4168 (unbexp == -398) 4169 && ((res1 & MASK_BINARY_SIG1) < 1000000000000000ull) 4170 && (is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0 4171 || is_midpoint_lt_even0 || is_midpoint_gt_even0)) { 4172 // set the inexact flag and the underflow flag 4173 *pfpsf |= (INEXACT_EXCEPTION | UNDERFLOW_EXCEPTION); 4174 } else if (is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0 || 4175 is_midpoint_lt_even0 || is_midpoint_gt_even0) { 4176 // set the inexact flag and the underflow flag 4177 *pfpsf |= INEXACT_EXCEPTION; 4178 } 4179 4180 *pfpsf |= save_fpsf; 4181 BID_RETURN (res1); 4182 } // else continue, and use rounding to nearest to round to 16 digits 4183 4184 // at this point the result is rounded to nearest (even or away) to 34 digits 4185 // (or less if exact), and it is zero or finite non-zero canonical [sub]normal 4186 sign = res.w[HIGH_128W] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 4187 exp = res.w[HIGH_128W] & MASK_EXP; // biased and shifted left 49 bits 4188 unbexp = (exp >> 49) - 6176; 4189 C.w[1] = res.w[HIGH_128W] & MASK_COEFF; 4190 C.w[0] = res.w[LOW_128W]; 4191 4192 if ((C.w[1] == 0x0 && C.w[0] == 0x0) || // result is zero 4193 (unbexp <= (-398 - 35)) || (unbexp >= (369 + 16))) { 4194 // clear under/overflow 4195 #if DECIMAL_CALL_BY_REFERENCE 4196 bid128_to_bid64 (&res1, &res _RND_MODE_ARG _EXC_FLAGS_ARG); 4197 #else 4198 res1 = bid128_to_bid64 (res _RND_MODE_ARG _EXC_FLAGS_ARG); 4199 #endif 4200 *pfpsf |= save_fpsf; 4201 BID_RETURN (res1); 4202 } // else continue 4203 4204 // -398 - 34 <= unbexp <= 369 + 15 4205 if (rnd_mode == ROUNDING_TIES_AWAY) { 4206 // apply correction, if needed, to make the result rounded to nearest-even 4207 if (is_midpoint_gt_even) { 4208 // res = res - 1 4209 res1--; // res1 is now even 4210 } // else the result is already correctly rounded to nearest-even 4211 } 4212 // at this point the result is finite, non-zero canonical normal or subnormal, 4213 // and in most cases overflow or underflow will not occur 4214 4215 // determine the number of digits q in the result 4216 // q = nr. of decimal digits in x 4217 // determine first the nr. of bits in x 4218 if (C.w[1] == 0) { 4219 if (C.w[0] >= 0x0020000000000000ull) { // x >= 2^53 4220 // split the 64-bit value in two 32-bit halves to avoid rounding errors 4221 if (C.w[0] >= 0x0000000100000000ull) { // x >= 2^32 4222 tmp.d = (double) (C.w[0] >> 32); // exact conversion 4223 nr_bits = 4224 33 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 4225 } else { // x < 2^32 4226 tmp.d = (double) (C.w[0]); // exact conversion 4227 nr_bits = 4228 1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 4229 } 4230 } else { // if x < 2^53 4231 tmp.d = (double) C.w[0]; // exact conversion 4232 nr_bits = 4233 1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 4234 } 4235 } else { // C.w[1] != 0 => nr. bits = 64 + nr_bits (C.w[1]) 4236 tmp.d = (double) C.w[1]; // exact conversion 4237 nr_bits = 4238 65 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff); 4239 } 4240 q = nr_digits[nr_bits - 1].digits; 4241 if (q == 0) { 4242 q = nr_digits[nr_bits - 1].digits1; 4243 if (C.w[1] > nr_digits[nr_bits - 1].threshold_hi || 4244 (C.w[1] == nr_digits[nr_bits - 1].threshold_hi && 4245 C.w[0] >= nr_digits[nr_bits - 1].threshold_lo)) 4246 q++; 4247 } 4248 // if q > 16, round to nearest even to 16 digits (but for underflow it may 4249 // have to be truncated even more) 4250 if (q > 16) { 4251 x0 = q - 16; 4252 if (q <= 18) { 4253 round64_2_18 (q, x0, C.w[0], &res1, &incr_exp, 4254 &is_midpoint_lt_even, &is_midpoint_gt_even, 4255 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint); 4256 } else { // 19 <= q <= 34 4257 round128_19_38 (q, x0, C, &res128, &incr_exp, 4258 &is_midpoint_lt_even, &is_midpoint_gt_even, 4259 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint); 4260 res1 = res128.w[0]; // the result fits in 64 bits 4261 } 4262 unbexp = unbexp + x0; 4263 if (incr_exp) 4264 unbexp++; 4265 q = 16; // need to set in case denormalization is necessary 4266 } else { 4267 // the result does not require a second rounding (and it must have 4268 // been exact in the first rounding, since q <= 16) 4269 res1 = C.w[0]; 4270 } 4271 4272 // avoid a double rounding error 4273 if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 4274 is_midpoint_lt_even) { // double rounding error upward 4275 // res = res - 1 4276 res1--; // res1 becomes odd 4277 is_midpoint_lt_even = 0; 4278 is_inexact_lt_midpoint = 1; 4279 if (res1 == 0x00038d7ea4c67fffull) { // 10^15 - 1 4280 res1 = 0x002386f26fc0ffffull; // 10^16 - 1 4281 unbexp--; 4282 } 4283 } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 4284 is_midpoint_gt_even) { // double rounding error downward 4285 // res = res + 1 4286 res1++; // res1 becomes odd (so it cannot be 10^16) 4287 is_midpoint_gt_even = 0; 4288 is_inexact_gt_midpoint = 1; 4289 } else if (!is_midpoint_lt_even && !is_midpoint_gt_even && 4290 !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) { 4291 // if this second rounding was exact the result may still be 4292 // inexact because of the first rounding 4293 if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) { 4294 is_inexact_gt_midpoint = 1; 4295 } 4296 if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) { 4297 is_inexact_lt_midpoint = 1; 4298 } 4299 } else if (is_midpoint_gt_even && 4300 (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) { 4301 // pulled up to a midpoint 4302 is_inexact_lt_midpoint = 1; 4303 is_inexact_gt_midpoint = 0; 4304 is_midpoint_lt_even = 0; 4305 is_midpoint_gt_even = 0; 4306 } else if (is_midpoint_lt_even && 4307 (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) { 4308 // pulled down to a midpoint 4309 is_inexact_lt_midpoint = 0; 4310 is_inexact_gt_midpoint = 1; 4311 is_midpoint_lt_even = 0; 4312 is_midpoint_gt_even = 0; 4313 } else { 4314 ; 4315 } 4316 // this is the result rounded correctly to nearest even, with unbounded exp. 4317 4318 // check for overflow 4319 if (q + unbexp > P16 + expmax16) { 4320 res1 = sign | 0x7800000000000000ull; 4321 *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION); 4322 *pfpsf |= save_fpsf; 4323 BID_RETURN (res1) 4324 } else if (unbexp > expmax16) { // q + unbexp <= P16 + expmax16 4325 // not overflow; the result must be exact, and we can multiply res1 by 4326 // 10^(unbexp - expmax16) and the product will fit in 16 decimal digits 4327 scale = unbexp - expmax16; 4328 res1 = res1 * ten2k64[scale]; // res1 * 10^scale 4329 unbexp = expmax16; // unbexp - scale 4330 } else { 4331 ; // continue 4332 } 4333 4334 // check for underflow 4335 if (q + unbexp < P16 + expmin16) { 4336 if (unbexp < expmin16) { 4337 // we must truncate more of res 4338 x0 = expmin16 - unbexp; // x0 >= 1 4339 is_inexact_lt_midpoint0 = is_inexact_lt_midpoint; 4340 is_inexact_gt_midpoint0 = is_inexact_gt_midpoint; 4341 is_midpoint_lt_even0 = is_midpoint_lt_even; 4342 is_midpoint_gt_even0 = is_midpoint_gt_even; 4343 is_inexact_lt_midpoint = 0; 4344 is_inexact_gt_midpoint = 0; 4345 is_midpoint_lt_even = 0; 4346 is_midpoint_gt_even = 0; 4347 // the number of decimal digits in res1 is q 4348 if (x0 < q) { // 1 <= x0 <= q-1 => round res to q - x0 digits 4349 // 2 <= q <= 16, 1 <= x0 <= 15 4350 round64_2_18 (q, x0, res1, &res1, &incr_exp, 4351 &is_midpoint_lt_even, &is_midpoint_gt_even, 4352 &is_inexact_lt_midpoint, &is_inexact_gt_midpoint); 4353 if (incr_exp) { 4354 // res1 = 10^(q-x0), 1 <= q - x0 <= q - 1, 1 <= q - x0 <= 15 4355 res1 = ten2k64[q - x0]; 4356 } 4357 unbexp = unbexp + x0; // expmin16 4358 } else if (x0 == q) { 4359 // the second rounding is for 0.d(0)d(1)...d(q-1) * 10^emin 4360 // determine relationship with 1/2 ulp 4361 // q <= 16 4362 if (res1 < midpoint64[q - 1]) { // < 1/2 ulp 4363 lt_half_ulp = 1; 4364 is_inexact_lt_midpoint = 1; 4365 } else if (res1 == midpoint64[q - 1]) { // = 1/2 ulp 4366 eq_half_ulp = 1; 4367 is_midpoint_gt_even = 1; 4368 } else { // > 1/2 ulp 4369 // gt_half_ulp = 1; 4370 is_inexact_gt_midpoint = 1; 4371 } 4372 if (lt_half_ulp || eq_half_ulp) { 4373 // res = +0.0 * 10^expmin16 4374 res1 = 0x0000000000000000ull; 4375 } else { // if (gt_half_ulp) 4376 // res = +1 * 10^expmin16 4377 res1 = 0x0000000000000001ull; 4378 } 4379 unbexp = expmin16; 4380 } else { // if (x0 > q) 4381 // the second rounding is for 0.0...d(0)d(1)...d(q-1) * 10^emin 4382 res1 = 0x0000000000000000ull; 4383 unbexp = expmin16; 4384 is_inexact_lt_midpoint = 1; 4385 } 4386 // avoid a double rounding error 4387 if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 4388 is_midpoint_lt_even) { // double rounding error upward 4389 // res = res - 1 4390 res1--; // res1 becomes odd 4391 is_midpoint_lt_even = 0; 4392 is_inexact_lt_midpoint = 1; 4393 } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 4394 is_midpoint_gt_even) { // double rounding error downward 4395 // res = res + 1 4396 res1++; // res1 becomes odd 4397 is_midpoint_gt_even = 0; 4398 is_inexact_gt_midpoint = 1; 4399 } else if (!is_midpoint_lt_even && !is_midpoint_gt_even && 4400 !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) { 4401 // if this rounding was exact the result may still be 4402 // inexact because of the previous roundings 4403 if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) { 4404 is_inexact_gt_midpoint = 1; 4405 } 4406 if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) { 4407 is_inexact_lt_midpoint = 1; 4408 } 4409 } else if (is_midpoint_gt_even && 4410 (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) { 4411 // pulled up to a midpoint 4412 is_inexact_lt_midpoint = 1; 4413 is_inexact_gt_midpoint = 0; 4414 is_midpoint_lt_even = 0; 4415 is_midpoint_gt_even = 0; 4416 } else if (is_midpoint_lt_even && 4417 (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) { 4418 // pulled down to a midpoint 4419 is_inexact_lt_midpoint = 0; 4420 is_inexact_gt_midpoint = 1; 4421 is_midpoint_lt_even = 0; 4422 is_midpoint_gt_even = 0; 4423 } else { 4424 ; 4425 } 4426 } 4427 // else if unbexp >= emin then q < P (because q + unbexp < P16 + expmin16) 4428 // and the result is tiny and exact 4429 4430 // check for inexact result 4431 if (is_inexact_lt_midpoint || is_inexact_gt_midpoint || 4432 is_midpoint_lt_even || is_midpoint_gt_even || 4433 is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0 || 4434 is_midpoint_lt_even0 || is_midpoint_gt_even0) { 4435 // set the inexact flag and the underflow flag 4436 *pfpsf |= (INEXACT_EXCEPTION | UNDERFLOW_EXCEPTION); 4437 } 4438 } else if (is_inexact_lt_midpoint || is_inexact_gt_midpoint || 4439 is_midpoint_lt_even || is_midpoint_gt_even) { 4440 *pfpsf |= INEXACT_EXCEPTION; 4441 } 4442 // this is the result rounded correctly to nearest, with bounded exponent 4443 4444 if (rnd_mode == ROUNDING_TIES_AWAY && is_midpoint_gt_even) { // correction 4445 // res = res + 1 4446 res1++; // res1 is now odd 4447 } // else the result is already correct 4448 4449 // assemble the result 4450 if (res1 < 0x0020000000000000ull) { // res < 2^53 4451 res1 = sign | ((UINT64) (unbexp + 398) << 53) | res1; 4452 } else { // res1 >= 2^53 4453 res1 = sign | MASK_STEERING_BITS | 4454 ((UINT64) (unbexp + 398) << 51) | (res1 & MASK_BINARY_SIG2); 4455 } 4456 *pfpsf |= save_fpsf; 4457 BID_RETURN (res1); 4458 } 4459
Indexes created Sat Feb 28 05:31:39 UTC 2026