1/* Copyright (C) 2007-2022 Free Software Foundation, Inc. 2 3This file is part of GCC. 4 5GCC is free software; you can redistribute it and/or modify it under 6the terms of the GNU General Public License as published by the Free 7Software Foundation; either version 3, or (at your option) any later 8version. 9 10GCC is distributed in the hope that it will be useful, but WITHOUT ANY 11WARRANTY; without even the implied warranty of MERCHANTABILITY or 12FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 13for more details. 14 15Under Section 7 of GPL version 3, you are granted additional 16permissions described in the GCC Runtime Library Exception, version 173.1, as published by the Free Software Foundation. 18 19You should have received a copy of the GNU General Public License and 20a copy of the GCC Runtime Library Exception along with this program; 21see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 22<http://www.gnu.org/licenses/>. */ 23 24/***************************************************************************** 25 * 26 * Helper add functions (for fma) 27 * 28 * __BID_INLINE__ UINT64 get_add64( 29 * UINT64 sign_x, int exponent_x, UINT64 coefficient_x, 30 * UINT64 sign_y, int exponent_y, UINT64 coefficient_y, 31 * int rounding_mode) 32 * 33 * __BID_INLINE__ UINT64 get_add128( 34 * UINT64 sign_x, int exponent_x, UINT64 coefficient_x, 35 * UINT64 sign_y, int final_exponent_y, UINT128 CY, 36 * int extra_digits, int rounding_mode) 37 * 38 ***************************************************************************** 39 * 40 * Algorithm description: 41 * 42 * get_add64: same as BID64 add, but arguments are unpacked and there 43 * are no special case checks 44 * 45 * get_add128: add 64-bit coefficient to 128-bit product (which contains 46 * 16+extra_digits decimal digits), 47 * return BID64 result 48 * - the exponents are compared and the two coefficients are 49 * properly aligned for addition/subtraction 50 * - multiple paths are needed 51 * - final result exponent is calculated and the lower term is 52 * rounded first if necessary, to avoid manipulating 53 * coefficients longer than 128 bits 54 * 55 ****************************************************************************/ 56 57#ifndef _INLINE_BID_ADD_H_ 58#define _INLINE_BID_ADD_H_ 59 60#include "bid_internal.h" 61 62#define MAX_FORMAT_DIGITS 16 63#define DECIMAL_EXPONENT_BIAS 398 64#define MASK_BINARY_EXPONENT 0x7ff0000000000000ull 65#define BINARY_EXPONENT_BIAS 0x3ff 66#define UPPER_EXPON_LIMIT 51 67 68/////////////////////////////////////////////////////////////////////// 69// 70// get_add64() is essentially the same as bid_add(), except that 71// the arguments are unpacked 72// 73////////////////////////////////////////////////////////////////////// 74__BID_INLINE__ UINT64 75get_add64 (UINT64 sign_x, int exponent_x, UINT64 coefficient_x, 76 UINT64 sign_y, int exponent_y, UINT64 coefficient_y, 77 int rounding_mode, unsigned *fpsc) { 78 UINT128 CA, CT, CT_new; 79 UINT64 sign_a, sign_b, coefficient_a, coefficient_b, sign_s, sign_ab, 80 rem_a; 81 UINT64 saved_ca, saved_cb, C0_64, C64, remainder_h, T1, carry, tmp, 82 C64_new; 83 int_double tempx; 84 int exponent_a, exponent_b, diff_dec_expon; 85 int bin_expon_ca, extra_digits, amount, scale_k, scale_ca; 86 unsigned rmode, status; 87 88 // sort arguments by exponent 89 if (exponent_x <= exponent_y) { 90 sign_a = sign_y; 91 exponent_a = exponent_y; 92 coefficient_a = coefficient_y; 93 sign_b = sign_x; 94 exponent_b = exponent_x; 95 coefficient_b = coefficient_x; 96 } else { 97 sign_a = sign_x; 98 exponent_a = exponent_x; 99 coefficient_a = coefficient_x; 100 sign_b = sign_y; 101 exponent_b = exponent_y; 102 coefficient_b = coefficient_y; 103 } 104 105 // exponent difference 106 diff_dec_expon = exponent_a - exponent_b; 107 108 /* get binary coefficients of x and y */ 109 110 //--- get number of bits in the coefficients of x and y --- 111 112 tempx.d = (double) coefficient_a; 113 bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; 114 115 if (!coefficient_a) { 116 return get_BID64 (sign_b, exponent_b, coefficient_b, rounding_mode, 117 fpsc); 118 } 119 if (diff_dec_expon > MAX_FORMAT_DIGITS) { 120 // normalize a to a 16-digit coefficient 121 122 scale_ca = estimate_decimal_digits[bin_expon_ca]; 123 if (coefficient_a >= power10_table_128[scale_ca].w[0]) 124 scale_ca++; 125 126 scale_k = 16 - scale_ca; 127 128 coefficient_a *= power10_table_128[scale_k].w[0]; 129 130 diff_dec_expon -= scale_k; 131 exponent_a -= scale_k; 132 133 /* get binary coefficients of x and y */ 134 135 //--- get number of bits in the coefficients of x and y --- 136 tempx.d = (double) coefficient_a; 137 bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; 138 139 if (diff_dec_expon > MAX_FORMAT_DIGITS) { 140#ifdef SET_STATUS_FLAGS 141 if (coefficient_b) { 142 __set_status_flags (fpsc, INEXACT_EXCEPTION); 143 } 144#endif 145 146#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 147#ifndef IEEE_ROUND_NEAREST 148 if (((rounding_mode) & 3) && coefficient_b) // not ROUNDING_TO_NEAREST 149 { 150 switch (rounding_mode) { 151 case ROUNDING_DOWN: 152 if (sign_b) { 153 coefficient_a -= ((((SINT64) sign_a) >> 63) | 1); 154 if (coefficient_a < 1000000000000000ull) { 155 exponent_a--; 156 coefficient_a = 9999999999999999ull; 157 } else if (coefficient_a >= 10000000000000000ull) { 158 exponent_a++; 159 coefficient_a = 1000000000000000ull; 160 } 161 } 162 break; 163 case ROUNDING_UP: 164 if (!sign_b) { 165 coefficient_a += ((((SINT64) sign_a) >> 63) | 1); 166 if (coefficient_a < 1000000000000000ull) { 167 exponent_a--; 168 coefficient_a = 9999999999999999ull; 169 } else if (coefficient_a >= 10000000000000000ull) { 170 exponent_a++; 171 coefficient_a = 1000000000000000ull; 172 } 173 } 174 break; 175 default: // RZ 176 if (sign_a != sign_b) { 177 coefficient_a--; 178 if (coefficient_a < 1000000000000000ull) { 179 exponent_a--; 180 coefficient_a = 9999999999999999ull; 181 } 182 } 183 break; 184 } 185 } else 186#endif 187#endif 188 // check special case here 189 if ((coefficient_a == 1000000000000000ull) 190 && (diff_dec_expon == MAX_FORMAT_DIGITS + 1) 191 && (sign_a ^ sign_b) 192 && (coefficient_b > 5000000000000000ull)) { 193 coefficient_a = 9999999999999999ull; 194 exponent_a--; 195 } 196 197 return get_BID64 (sign_a, exponent_a, coefficient_a, 198 rounding_mode, fpsc); 199 } 200 } 201 // test whether coefficient_a*10^(exponent_a-exponent_b) may exceed 2^62 202 if (bin_expon_ca + estimate_bin_expon[diff_dec_expon] < 60) { 203 // coefficient_a*10^(exponent_a-exponent_b)<2^63 204 205 // multiply by 10^(exponent_a-exponent_b) 206 coefficient_a *= power10_table_128[diff_dec_expon].w[0]; 207 208 // sign mask 209 sign_b = ((SINT64) sign_b) >> 63; 210 // apply sign to coeff. of b 211 coefficient_b = (coefficient_b + sign_b) ^ sign_b; 212 213 // apply sign to coefficient a 214 sign_a = ((SINT64) sign_a) >> 63; 215 coefficient_a = (coefficient_a + sign_a) ^ sign_a; 216 217 coefficient_a += coefficient_b; 218 // get sign 219 sign_s = ((SINT64) coefficient_a) >> 63; 220 coefficient_a = (coefficient_a + sign_s) ^ sign_s; 221 sign_s &= 0x8000000000000000ull; 222 223 // coefficient_a < 10^16 ? 224 if (coefficient_a < power10_table_128[MAX_FORMAT_DIGITS].w[0]) { 225#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 226#ifndef IEEE_ROUND_NEAREST 227 if (rounding_mode == ROUNDING_DOWN && (!coefficient_a) 228 && sign_a != sign_b) 229 sign_s = 0x8000000000000000ull; 230#endif 231#endif 232 return get_BID64 (sign_s, exponent_b, coefficient_a, 233 rounding_mode, fpsc); 234 } 235 // otherwise rounding is necessary 236 237 // already know coefficient_a<10^19 238 // coefficient_a < 10^17 ? 239 if (coefficient_a < power10_table_128[17].w[0]) 240 extra_digits = 1; 241 else if (coefficient_a < power10_table_128[18].w[0]) 242 extra_digits = 2; 243 else 244 extra_digits = 3; 245 246#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 247#ifndef IEEE_ROUND_NEAREST 248 rmode = rounding_mode; 249 if (sign_s && (unsigned) (rmode - 1) < 2) 250 rmode = 3 - rmode; 251#else 252 rmode = 0; 253#endif 254#else 255 rmode = 0; 256#endif 257 coefficient_a += round_const_table[rmode][extra_digits]; 258 259 // get P*(2^M[extra_digits])/10^extra_digits 260 __mul_64x64_to_128 (CT, coefficient_a, 261 reciprocals10_64[extra_digits]); 262 263 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128 264 amount = short_recip_scale[extra_digits]; 265 C64 = CT.w[1] >> amount; 266 267 } else { 268 // coefficient_a*10^(exponent_a-exponent_b) is large 269 sign_s = sign_a; 270 271#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 272#ifndef IEEE_ROUND_NEAREST 273 rmode = rounding_mode; 274 if (sign_s && (unsigned) (rmode - 1) < 2) 275 rmode = 3 - rmode; 276#else 277 rmode = 0; 278#endif 279#else 280 rmode = 0; 281#endif 282 283 // check whether we can take faster path 284 scale_ca = estimate_decimal_digits[bin_expon_ca]; 285 286 sign_ab = sign_a ^ sign_b; 287 sign_ab = ((SINT64) sign_ab) >> 63; 288 289 // T1 = 10^(16-diff_dec_expon) 290 T1 = power10_table_128[16 - diff_dec_expon].w[0]; 291 292 // get number of digits in coefficient_a 293 //P_ca = power10_table_128[scale_ca].w[0]; 294 //P_ca_m1 = power10_table_128[scale_ca-1].w[0]; 295 if (coefficient_a >= power10_table_128[scale_ca].w[0]) { 296 scale_ca++; 297 //P_ca_m1 = P_ca; 298 //P_ca = power10_table_128[scale_ca].w[0]; 299 } 300 301 scale_k = 16 - scale_ca; 302 303 // apply sign 304 //Ts = (T1 + sign_ab) ^ sign_ab; 305 306 // test range of ca 307 //X = coefficient_a + Ts - P_ca_m1; 308 309 // addition 310 saved_ca = coefficient_a - T1; 311 coefficient_a = 312 (SINT64) saved_ca *(SINT64) power10_table_128[scale_k].w[0]; 313 extra_digits = diff_dec_expon - scale_k; 314 315 // apply sign 316 saved_cb = (coefficient_b + sign_ab) ^ sign_ab; 317 // add 10^16 and rounding constant 318 coefficient_b = 319 saved_cb + 10000000000000000ull + 320 round_const_table[rmode][extra_digits]; 321 322 // get P*(2^M[extra_digits])/10^extra_digits 323 __mul_64x64_to_128 (CT, coefficient_b, 324 reciprocals10_64[extra_digits]); 325 326 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128 327 amount = short_recip_scale[extra_digits]; 328 C0_64 = CT.w[1] >> amount; 329 330 // result coefficient 331 C64 = C0_64 + coefficient_a; 332 // filter out difficult (corner) cases 333 // the following test is equivalent to 334 // ( (initial_coefficient_a + Ts) < P_ca && 335 // (initial_coefficient_a + Ts) > P_ca_m1 ), 336 // which ensures the number of digits in coefficient_a does not change 337 // after adding (the appropriately scaled and rounded) coefficient_b 338 if ((UINT64) (C64 - 1000000000000000ull - 1) > 339 9000000000000000ull - 2) { 340 if (C64 >= 10000000000000000ull) { 341 // result has more than 16 digits 342 if (!scale_k) { 343 // must divide coeff_a by 10 344 saved_ca = saved_ca + T1; 345 __mul_64x64_to_128 (CA, saved_ca, 0x3333333333333334ull); 346 //reciprocals10_64[1]); 347 coefficient_a = CA.w[1] >> 1; 348 rem_a = 349 saved_ca - (coefficient_a << 3) - (coefficient_a << 1); 350 coefficient_a = coefficient_a - T1; 351 352 saved_cb += 353 /*90000000000000000 */ +rem_a * 354 power10_table_128[diff_dec_expon].w[0]; 355 } else 356 coefficient_a = 357 (SINT64) (saved_ca - T1 - 358 (T1 << 3)) * (SINT64) power10_table_128[scale_k - 359 1].w[0]; 360 361 extra_digits++; 362 coefficient_b = 363 saved_cb + 100000000000000000ull + 364 round_const_table[rmode][extra_digits]; 365 366 // get P*(2^M[extra_digits])/10^extra_digits 367 __mul_64x64_to_128 (CT, coefficient_b, 368 reciprocals10_64[extra_digits]); 369 370 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128 371 amount = short_recip_scale[extra_digits]; 372 C0_64 = CT.w[1] >> amount; 373 374 // result coefficient 375 C64 = C0_64 + coefficient_a; 376 } else if (C64 <= 1000000000000000ull) { 377 // less than 16 digits in result 378 coefficient_a = 379 (SINT64) saved_ca *(SINT64) power10_table_128[scale_k + 380 1].w[0]; 381 //extra_digits --; 382 exponent_b--; 383 coefficient_b = 384 (saved_cb << 3) + (saved_cb << 1) + 100000000000000000ull + 385 round_const_table[rmode][extra_digits]; 386 387 // get P*(2^M[extra_digits])/10^extra_digits 388 __mul_64x64_to_128 (CT_new, coefficient_b, 389 reciprocals10_64[extra_digits]); 390 391 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128 392 amount = short_recip_scale[extra_digits]; 393 C0_64 = CT_new.w[1] >> amount; 394 395 // result coefficient 396 C64_new = C0_64 + coefficient_a; 397 if (C64_new < 10000000000000000ull) { 398 C64 = C64_new; 399#ifdef SET_STATUS_FLAGS 400 CT = CT_new; 401#endif 402 } else 403 exponent_b++; 404 } 405 406 } 407 408 } 409 410#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 411#ifndef IEEE_ROUND_NEAREST 412 if (rmode == 0) //ROUNDING_TO_NEAREST 413#endif 414 if (C64 & 1) { 415 // check whether fractional part of initial_P/10^extra_digits 416 // is exactly .5 417 // this is the same as fractional part of 418 // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero 419 420 // get remainder 421 remainder_h = CT.w[1] << (64 - amount); 422 423 // test whether fractional part is 0 424 if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) { 425 C64--; 426 } 427 } 428#endif 429 430#ifdef SET_STATUS_FLAGS 431 status = INEXACT_EXCEPTION; 432 433 // get remainder 434 remainder_h = CT.w[1] << (64 - amount); 435 436 switch (rmode) { 437 case ROUNDING_TO_NEAREST: 438 case ROUNDING_TIES_AWAY: 439 // test whether fractional part is 0 440 if ((remainder_h == 0x8000000000000000ull) 441 && (CT.w[0] < reciprocals10_64[extra_digits])) 442 status = EXACT_STATUS; 443 break; 444 case ROUNDING_DOWN: 445 case ROUNDING_TO_ZERO: 446 if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) 447 status = EXACT_STATUS; 448 break; 449 default: 450 // round up 451 __add_carry_out (tmp, carry, CT.w[0], 452 reciprocals10_64[extra_digits]); 453 if ((remainder_h >> (64 - amount)) + carry >= 454 (((UINT64) 1) << amount)) 455 status = EXACT_STATUS; 456 break; 457 } 458 __set_status_flags (fpsc, status); 459 460#endif 461 462 return get_BID64 (sign_s, exponent_b + extra_digits, C64, 463 rounding_mode, fpsc); 464} 465 466 467/////////////////////////////////////////////////////////////////// 468// round 128-bit coefficient and return result in BID64 format 469// do not worry about midpoint cases 470////////////////////////////////////////////////////////////////// 471static UINT64 472__bid_simple_round64_sticky (UINT64 sign, int exponent, UINT128 P, 473 int extra_digits, int rounding_mode, 474 unsigned *fpsc) { 475 UINT128 Q_high, Q_low, C128; 476 UINT64 C64; 477 int amount, rmode; 478 479 rmode = rounding_mode; 480#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 481#ifndef IEEE_ROUND_NEAREST 482 if (sign && (unsigned) (rmode - 1) < 2) 483 rmode = 3 - rmode; 484#endif 485#endif 486 __add_128_64 (P, P, round_const_table[rmode][extra_digits]); 487 488 // get P*(2^M[extra_digits])/10^extra_digits 489 __mul_128x128_full (Q_high, Q_low, P, 490 reciprocals10_128[extra_digits]); 491 492 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 493 amount = recip_scale[extra_digits]; 494 __shr_128 (C128, Q_high, amount); 495 496 C64 = __low_64 (C128); 497 498#ifdef SET_STATUS_FLAGS 499 500 __set_status_flags (fpsc, INEXACT_EXCEPTION); 501 502#endif 503 504 return get_BID64 (sign, exponent, C64, rounding_mode, fpsc); 505} 506 507/////////////////////////////////////////////////////////////////// 508// round 128-bit coefficient and return result in BID64 format 509/////////////////////////////////////////////////////////////////// 510static UINT64 511__bid_full_round64 (UINT64 sign, int exponent, UINT128 P, 512 int extra_digits, int rounding_mode, 513 unsigned *fpsc) { 514 UINT128 Q_high, Q_low, C128, Stemp, PU; 515 UINT64 remainder_h, C64, carry, CY; 516 int amount, amount2, rmode, status = 0; 517 518 if (exponent < 0) { 519 if (exponent >= -16 && (extra_digits + exponent < 0)) { 520 extra_digits = -exponent; 521#ifdef SET_STATUS_FLAGS 522 if (extra_digits > 0) { 523 rmode = rounding_mode; 524 if (sign && (unsigned) (rmode - 1) < 2) 525 rmode = 3 - rmode; 526 __add_128_128 (PU, P, 527 round_const_table_128[rmode][extra_digits]); 528 if (__unsigned_compare_gt_128 529 (power10_table_128[extra_digits + 15], PU)) 530 status = UNDERFLOW_EXCEPTION; 531 } 532#endif 533 } 534 } 535 536 if (extra_digits > 0) { 537 exponent += extra_digits; 538 rmode = rounding_mode; 539 if (sign && (unsigned) (rmode - 1) < 2) 540 rmode = 3 - rmode; 541 __add_128_128 (P, P, round_const_table_128[rmode][extra_digits]); 542 543 // get P*(2^M[extra_digits])/10^extra_digits 544 __mul_128x128_full (Q_high, Q_low, P, 545 reciprocals10_128[extra_digits]); 546 547 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 548 amount = recip_scale[extra_digits]; 549 __shr_128_long (C128, Q_high, amount); 550 551 C64 = __low_64 (C128); 552 553#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 554#ifndef IEEE_ROUND_NEAREST 555 if (rmode == 0) //ROUNDING_TO_NEAREST 556#endif 557 if (C64 & 1) { 558 // check whether fractional part of initial_P/10^extra_digits 559 // is exactly .5 560 561 // get remainder 562 amount2 = 64 - amount; 563 remainder_h = 0; 564 remainder_h--; 565 remainder_h >>= amount2; 566 remainder_h = remainder_h & Q_high.w[0]; 567 568 if (!remainder_h 569 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] 570 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] 571 && Q_low.w[0] < 572 reciprocals10_128[extra_digits].w[0]))) { 573 C64--; 574 } 575 } 576#endif 577 578#ifdef SET_STATUS_FLAGS 579 status |= INEXACT_EXCEPTION; 580 581 // get remainder 582 remainder_h = Q_high.w[0] << (64 - amount); 583 584 switch (rmode) { 585 case ROUNDING_TO_NEAREST: 586 case ROUNDING_TIES_AWAY: 587 // test whether fractional part is 0 588 if (remainder_h == 0x8000000000000000ull 589 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] 590 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] 591 && Q_low.w[0] < 592 reciprocals10_128[extra_digits].w[0]))) 593 status = EXACT_STATUS; 594 break; 595 case ROUNDING_DOWN: 596 case ROUNDING_TO_ZERO: 597 if (!remainder_h 598 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] 599 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] 600 && Q_low.w[0] < 601 reciprocals10_128[extra_digits].w[0]))) 602 status = EXACT_STATUS; 603 break; 604 default: 605 // round up 606 __add_carry_out (Stemp.w[0], CY, Q_low.w[0], 607 reciprocals10_128[extra_digits].w[0]); 608 __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1], 609 reciprocals10_128[extra_digits].w[1], CY); 610 if ((remainder_h >> (64 - amount)) + carry >= 611 (((UINT64) 1) << amount)) 612 status = EXACT_STATUS; 613 } 614 615 __set_status_flags (fpsc, status); 616 617#endif 618 } else { 619 C64 = P.w[0]; 620 if (!C64) { 621 sign = 0; 622 if (rounding_mode == ROUNDING_DOWN) 623 sign = 0x8000000000000000ull; 624 } 625 } 626 return get_BID64 (sign, exponent, C64, rounding_mode, fpsc); 627} 628 629///////////////////////////////////////////////////////////////////////////////// 630// round 192-bit coefficient (P, remainder_P) and return result in BID64 format 631// the lowest 64 bits (remainder_P) are used for midpoint checking only 632//////////////////////////////////////////////////////////////////////////////// 633static UINT64 634__bid_full_round64_remainder (UINT64 sign, int exponent, UINT128 P, 635 int extra_digits, UINT64 remainder_P, 636 int rounding_mode, unsigned *fpsc, 637 unsigned uf_status) { 638 UINT128 Q_high, Q_low, C128, Stemp; 639 UINT64 remainder_h, C64, carry, CY; 640 int amount, amount2, rmode, status = uf_status; 641 642 rmode = rounding_mode; 643 if (sign && (unsigned) (rmode - 1) < 2) 644 rmode = 3 - rmode; 645 if (rmode == ROUNDING_UP && remainder_P) { 646 P.w[0]++; 647 if (!P.w[0]) 648 P.w[1]++; 649 } 650 651 if (extra_digits) { 652 __add_128_64 (P, P, round_const_table[rmode][extra_digits]); 653 654 // get P*(2^M[extra_digits])/10^extra_digits 655 __mul_128x128_full (Q_high, Q_low, P, 656 reciprocals10_128[extra_digits]); 657 658 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 659 amount = recip_scale[extra_digits]; 660 __shr_128 (C128, Q_high, amount); 661 662 C64 = __low_64 (C128); 663 664#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 665#ifndef IEEE_ROUND_NEAREST 666 if (rmode == 0) //ROUNDING_TO_NEAREST 667#endif 668 if (!remainder_P && (C64 & 1)) { 669 // check whether fractional part of initial_P/10^extra_digits 670 // is exactly .5 671 672 // get remainder 673 amount2 = 64 - amount; 674 remainder_h = 0; 675 remainder_h--; 676 remainder_h >>= amount2; 677 remainder_h = remainder_h & Q_high.w[0]; 678 679 if (!remainder_h 680 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] 681 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] 682 && Q_low.w[0] < 683 reciprocals10_128[extra_digits].w[0]))) { 684 C64--; 685 } 686 } 687#endif 688 689#ifdef SET_STATUS_FLAGS 690 status |= INEXACT_EXCEPTION; 691 692 if (!remainder_P) { 693 // get remainder 694 remainder_h = Q_high.w[0] << (64 - amount); 695 696 switch (rmode) { 697 case ROUNDING_TO_NEAREST: 698 case ROUNDING_TIES_AWAY: 699 // test whether fractional part is 0 700 if (remainder_h == 0x8000000000000000ull 701 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] 702 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] 703 && Q_low.w[0] < 704 reciprocals10_128[extra_digits].w[0]))) 705 status = EXACT_STATUS; 706 break; 707 case ROUNDING_DOWN: 708 case ROUNDING_TO_ZERO: 709 if (!remainder_h 710 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] 711 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] 712 && Q_low.w[0] < 713 reciprocals10_128[extra_digits].w[0]))) 714 status = EXACT_STATUS; 715 break; 716 default: 717 // round up 718 __add_carry_out (Stemp.w[0], CY, Q_low.w[0], 719 reciprocals10_128[extra_digits].w[0]); 720 __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1], 721 reciprocals10_128[extra_digits].w[1], CY); 722 if ((remainder_h >> (64 - amount)) + carry >= 723 (((UINT64) 1) << amount)) 724 status = EXACT_STATUS; 725 } 726 } 727 __set_status_flags (fpsc, status); 728 729#endif 730 } else { 731 C64 = P.w[0]; 732#ifdef SET_STATUS_FLAGS 733 if (remainder_P) { 734 __set_status_flags (fpsc, uf_status | INEXACT_EXCEPTION); 735 } 736#endif 737 } 738 739 return get_BID64 (sign, exponent + extra_digits, C64, rounding_mode, 740 fpsc); 741} 742 743 744/////////////////////////////////////////////////////////////////// 745// get P/10^extra_digits 746// result fits in 64 bits 747/////////////////////////////////////////////////////////////////// 748__BID_INLINE__ UINT64 749__truncate (UINT128 P, int extra_digits) 750// extra_digits <= 16 751{ 752 UINT128 Q_high, Q_low, C128; 753 UINT64 C64; 754 int amount; 755 756 // get P*(2^M[extra_digits])/10^extra_digits 757 __mul_128x128_full (Q_high, Q_low, P, 758 reciprocals10_128[extra_digits]); 759 760 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 761 amount = recip_scale[extra_digits]; 762 __shr_128 (C128, Q_high, amount); 763 764 C64 = __low_64 (C128); 765 766 return C64; 767} 768 769 770/////////////////////////////////////////////////////////////////// 771// return number of decimal digits in 128-bit value X 772/////////////////////////////////////////////////////////////////// 773__BID_INLINE__ int 774__get_dec_digits64 (UINT128 X) { 775 int_double tempx; 776 int digits_x, bin_expon_cx; 777 778 if (!X.w[1]) { 779 //--- get number of bits in the coefficients of x and y --- 780 tempx.d = (double) X.w[0]; 781 bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; 782 // get number of decimal digits in the coeff_x 783 digits_x = estimate_decimal_digits[bin_expon_cx]; 784 if (X.w[0] >= power10_table_128[digits_x].w[0]) 785 digits_x++; 786 return digits_x; 787 } 788 tempx.d = (double) X.w[1]; 789 bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; 790 // get number of decimal digits in the coeff_x 791 digits_x = estimate_decimal_digits[bin_expon_cx + 64]; 792 if (__unsigned_compare_ge_128 (X, power10_table_128[digits_x])) 793 digits_x++; 794 795 return digits_x; 796} 797 798 799//////////////////////////////////////////////////////////////////////////////// 800// 801// add 64-bit coefficient to 128-bit coefficient, return result in BID64 format 802// 803//////////////////////////////////////////////////////////////////////////////// 804__BID_INLINE__ UINT64 805get_add128 (UINT64 sign_x, int exponent_x, UINT64 coefficient_x, 806 UINT64 sign_y, int final_exponent_y, UINT128 CY, 807 int extra_digits, int rounding_mode, unsigned *fpsc) { 808 UINT128 CY_L, CX, FS, F, CT, ST, T2; 809 UINT64 CYh, CY0L, T, S, coefficient_y, remainder_y; 810 SINT64 D = 0; 811 int_double tempx; 812 int diff_dec_expon, extra_digits2, exponent_y, status; 813 int extra_dx, diff_dec2, bin_expon_cx, digits_x, rmode; 814 815 // CY has more than 16 decimal digits 816 817 exponent_y = final_exponent_y - extra_digits; 818 819#ifdef IEEE_ROUND_NEAREST_TIES_AWAY 820 rounding_mode = 0; 821#endif 822#ifdef IEEE_ROUND_NEAREST 823 rounding_mode = 0; 824#endif 825 826 if (exponent_x > exponent_y) { 827 // normalize x 828 //--- get number of bits in the coefficients of x and y --- 829 tempx.d = (double) coefficient_x; 830 bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; 831 // get number of decimal digits in the coeff_x 832 digits_x = estimate_decimal_digits[bin_expon_cx]; 833 if (coefficient_x >= power10_table_128[digits_x].w[0]) 834 digits_x++; 835 836 extra_dx = 16 - digits_x; 837 coefficient_x *= power10_table_128[extra_dx].w[0]; 838 if ((sign_x ^ sign_y) && (coefficient_x == 1000000000000000ull)) { 839 extra_dx++; 840 coefficient_x = 10000000000000000ull; 841 } 842 exponent_x -= extra_dx; 843 844 if (exponent_x > exponent_y) { 845 846 // exponent_x > exponent_y 847 diff_dec_expon = exponent_x - exponent_y; 848 849 if (exponent_x <= final_exponent_y + 1) { 850 __mul_64x64_to_128 (CX, coefficient_x, 851 power10_table_128[diff_dec_expon].w[0]); 852 853 if (sign_x == sign_y) { 854 __add_128_128 (CT, CY, CX); 855 if ((exponent_x > 856 final_exponent_y) /*&& (final_exponent_y>0) */ ) 857 extra_digits++; 858 if (__unsigned_compare_ge_128 859 (CT, power10_table_128[16 + extra_digits])) 860 extra_digits++; 861 } else { 862 __sub_128_128 (CT, CY, CX); 863 if (((SINT64) CT.w[1]) < 0) { 864 CT.w[0] = 0 - CT.w[0]; 865 CT.w[1] = 0 - CT.w[1]; 866 if (CT.w[0]) 867 CT.w[1]--; 868 sign_y = sign_x; 869 } else if (!(CT.w[1] | CT.w[0])) { 870 sign_y = 871 (rounding_mode != 872 ROUNDING_DOWN) ? 0 : 0x8000000000000000ull; 873 } 874 if ((exponent_x + 1 >= 875 final_exponent_y) /*&& (final_exponent_y>=0) */ ) { 876 extra_digits = __get_dec_digits64 (CT) - 16; 877 if (extra_digits <= 0) { 878 if (!CT.w[0] && rounding_mode == ROUNDING_DOWN) 879 sign_y = 0x8000000000000000ull; 880 return get_BID64 (sign_y, exponent_y, CT.w[0], 881 rounding_mode, fpsc); 882 } 883 } else 884 if (__unsigned_compare_gt_128 885 (power10_table_128[15 + extra_digits], CT)) 886 extra_digits--; 887 } 888 889 return __bid_full_round64 (sign_y, exponent_y, CT, extra_digits, 890 rounding_mode, fpsc); 891 } 892 // diff_dec2+extra_digits is the number of digits to eliminate from 893 // argument CY 894 diff_dec2 = exponent_x - final_exponent_y; 895 896 if (diff_dec2 >= 17) { 897#ifndef IEEE_ROUND_NEAREST 898#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 899 if ((rounding_mode) & 3) { 900 switch (rounding_mode) { 901 case ROUNDING_UP: 902 if (!sign_y) { 903 D = ((SINT64) (sign_x ^ sign_y)) >> 63; 904 D = D + D + 1; 905 coefficient_x += D; 906 } 907 break; 908 case ROUNDING_DOWN: 909 if (sign_y) { 910 D = ((SINT64) (sign_x ^ sign_y)) >> 63; 911 D = D + D + 1; 912 coefficient_x += D; 913 } 914 break; 915 case ROUNDING_TO_ZERO: 916 if (sign_y != sign_x) { 917 D = 0 - 1; 918 coefficient_x += D; 919 } 920 break; 921 } 922 if (coefficient_x < 1000000000000000ull) { 923 coefficient_x -= D; 924 coefficient_x = 925 D + (coefficient_x << 1) + (coefficient_x << 3); 926 exponent_x--; 927 } 928 } 929#endif 930#endif 931#ifdef SET_STATUS_FLAGS 932 if (CY.w[1] | CY.w[0]) 933 __set_status_flags (fpsc, INEXACT_EXCEPTION); 934#endif 935 return get_BID64 (sign_x, exponent_x, coefficient_x, 936 rounding_mode, fpsc); 937 } 938 // here exponent_x <= 16+final_exponent_y 939 940 // truncate CY to 16 dec. digits 941 CYh = __truncate (CY, extra_digits); 942 943 // get remainder 944 T = power10_table_128[extra_digits].w[0]; 945 __mul_64x64_to_64 (CY0L, CYh, T); 946 947 remainder_y = CY.w[0] - CY0L; 948 949 // align coeff_x, CYh 950 __mul_64x64_to_128 (CX, coefficient_x, 951 power10_table_128[diff_dec2].w[0]); 952 953 if (sign_x == sign_y) { 954 __add_128_64 (CT, CX, CYh); 955 if (__unsigned_compare_ge_128 956 (CT, power10_table_128[16 + diff_dec2])) 957 diff_dec2++; 958 } else { 959 if (remainder_y) 960 CYh++; 961 __sub_128_64 (CT, CX, CYh); 962 if (__unsigned_compare_gt_128 963 (power10_table_128[15 + diff_dec2], CT)) 964 diff_dec2--; 965 } 966 967 return __bid_full_round64_remainder (sign_x, final_exponent_y, CT, 968 diff_dec2, remainder_y, 969 rounding_mode, fpsc, 0); 970 } 971 } 972 // Here (exponent_x <= exponent_y) 973 { 974 diff_dec_expon = exponent_y - exponent_x; 975 976 if (diff_dec_expon > MAX_FORMAT_DIGITS) { 977 rmode = rounding_mode; 978 979 if ((sign_x ^ sign_y)) { 980 if (!CY.w[0]) 981 CY.w[1]--; 982 CY.w[0]--; 983 if (__unsigned_compare_gt_128 984 (power10_table_128[15 + extra_digits], CY)) { 985 if (rmode & 3) { 986 extra_digits--; 987 final_exponent_y--; 988 } else { 989 CY.w[0] = 1000000000000000ull; 990 CY.w[1] = 0; 991 extra_digits = 0; 992 } 993 } 994 } 995 __scale128_10 (CY, CY); 996 extra_digits++; 997 CY.w[0] |= 1; 998 999 return __bid_simple_round64_sticky (sign_y, final_exponent_y, CY, 1000 extra_digits, rmode, fpsc); 1001 } 1002 // apply sign to coeff_x 1003 sign_x ^= sign_y; 1004 sign_x = ((SINT64) sign_x) >> 63; 1005 CX.w[0] = (coefficient_x + sign_x) ^ sign_x; 1006 CX.w[1] = sign_x; 1007 1008 // check whether CY (rounded to 16 digits) and CX have 1009 // any digits in the same position 1010 diff_dec2 = final_exponent_y - exponent_x; 1011 1012 if (diff_dec2 <= 17) { 1013 // align CY to 10^ex 1014 S = power10_table_128[diff_dec_expon].w[0]; 1015 __mul_64x128_short (CY_L, S, CY); 1016 1017 __add_128_128 (ST, CY_L, CX); 1018 extra_digits2 = __get_dec_digits64 (ST) - 16; 1019 return __bid_full_round64 (sign_y, exponent_x, ST, extra_digits2, 1020 rounding_mode, fpsc); 1021 } 1022 // truncate CY to 16 dec. digits 1023 CYh = __truncate (CY, extra_digits); 1024 1025 // get remainder 1026 T = power10_table_128[extra_digits].w[0]; 1027 __mul_64x64_to_64 (CY0L, CYh, T); 1028 1029 coefficient_y = CY.w[0] - CY0L; 1030 // add rounding constant 1031 rmode = rounding_mode; 1032 if (sign_y && (unsigned) (rmode - 1) < 2) 1033 rmode = 3 - rmode; 1034#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 1035#ifndef IEEE_ROUND_NEAREST 1036 if (!(rmode & 3)) //ROUNDING_TO_NEAREST 1037#endif 1038#endif 1039 { 1040 coefficient_y += round_const_table[rmode][extra_digits]; 1041 } 1042 // align coefficient_y, coefficient_x 1043 S = power10_table_128[diff_dec_expon].w[0]; 1044 __mul_64x64_to_128 (F, coefficient_y, S); 1045 1046 // fraction 1047 __add_128_128 (FS, F, CX); 1048 1049#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 1050#ifndef IEEE_ROUND_NEAREST 1051 if (rmode == 0) //ROUNDING_TO_NEAREST 1052#endif 1053 { 1054 // rounding code, here RN_EVEN 1055 // 10^(extra_digits+diff_dec_expon) 1056 T2 = power10_table_128[diff_dec_expon + extra_digits]; 1057 if (__unsigned_compare_gt_128 (FS, T2) 1058 || ((CYh & 1) && __test_equal_128 (FS, T2))) { 1059 CYh++; 1060 __sub_128_128 (FS, FS, T2); 1061 } 1062 } 1063#endif 1064#ifndef IEEE_ROUND_NEAREST 1065#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 1066 if (rmode == 4) //ROUNDING_TO_NEAREST 1067#endif 1068 { 1069 // rounding code, here RN_AWAY 1070 // 10^(extra_digits+diff_dec_expon) 1071 T2 = power10_table_128[diff_dec_expon + extra_digits]; 1072 if (__unsigned_compare_ge_128 (FS, T2)) { 1073 CYh++; 1074 __sub_128_128 (FS, FS, T2); 1075 } 1076 } 1077#endif 1078#ifndef IEEE_ROUND_NEAREST 1079#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 1080 switch (rmode) { 1081 case ROUNDING_DOWN: 1082 case ROUNDING_TO_ZERO: 1083 if ((SINT64) FS.w[1] < 0) { 1084 CYh--; 1085 if (CYh < 1000000000000000ull) { 1086 CYh = 9999999999999999ull; 1087 final_exponent_y--; 1088 } 1089 } else { 1090 T2 = power10_table_128[diff_dec_expon + extra_digits]; 1091 if (__unsigned_compare_ge_128 (FS, T2)) { 1092 CYh++; 1093 __sub_128_128 (FS, FS, T2); 1094 } 1095 } 1096 break; 1097 case ROUNDING_UP: 1098 if ((SINT64) FS.w[1] < 0) 1099 break; 1100 T2 = power10_table_128[diff_dec_expon + extra_digits]; 1101 if (__unsigned_compare_gt_128 (FS, T2)) { 1102 CYh += 2; 1103 __sub_128_128 (FS, FS, T2); 1104 } else if ((FS.w[1] == T2.w[1]) && (FS.w[0] == T2.w[0])) { 1105 CYh++; 1106 FS.w[1] = FS.w[0] = 0; 1107 } else if (FS.w[1] | FS.w[0]) 1108 CYh++; 1109 break; 1110 } 1111#endif 1112#endif 1113 1114#ifdef SET_STATUS_FLAGS 1115 status = INEXACT_EXCEPTION; 1116#ifndef IEEE_ROUND_NEAREST 1117#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 1118 if (!(rmode & 3)) 1119#endif 1120#endif 1121 { 1122 // RN modes 1123 if ((FS.w[1] == 1124 round_const_table_128[0][diff_dec_expon + extra_digits].w[1]) 1125 && (FS.w[0] == 1126 round_const_table_128[0][diff_dec_expon + 1127 extra_digits].w[0])) 1128 status = EXACT_STATUS; 1129 } 1130#ifndef IEEE_ROUND_NEAREST 1131#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 1132 else if (!FS.w[1] && !FS.w[0]) 1133 status = EXACT_STATUS; 1134#endif 1135#endif 1136 1137 __set_status_flags (fpsc, status); 1138#endif 1139 1140 return get_BID64 (sign_y, final_exponent_y, CYh, rounding_mode, 1141 fpsc); 1142 } 1143 1144} 1145 1146////////////////////////////////////////////////////////////////////////// 1147// 1148// If coefficient_z is less than 16 digits long, normalize to 16 digits 1149// 1150///////////////////////////////////////////////////////////////////////// 1151static UINT64 1152BID_normalize (UINT64 sign_z, int exponent_z, 1153 UINT64 coefficient_z, UINT64 round_dir, int round_flag, 1154 int rounding_mode, unsigned *fpsc) { 1155 SINT64 D; 1156 int_double tempx; 1157 int digits_z, bin_expon, scale, rmode; 1158 1159#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 1160#ifndef IEEE_ROUND_NEAREST 1161 rmode = rounding_mode; 1162 if (sign_z && (unsigned) (rmode - 1) < 2) 1163 rmode = 3 - rmode; 1164#else 1165 if (coefficient_z >= power10_table_128[15].w[0]) 1166 return z; 1167#endif 1168#endif 1169 1170 //--- get number of bits in the coefficients of x and y --- 1171 tempx.d = (double) coefficient_z; 1172 bin_expon = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; 1173 // get number of decimal digits in the coeff_x 1174 digits_z = estimate_decimal_digits[bin_expon]; 1175 if (coefficient_z >= power10_table_128[digits_z].w[0]) 1176 digits_z++; 1177 1178 scale = 16 - digits_z; 1179 exponent_z -= scale; 1180 if (exponent_z < 0) { 1181 scale += exponent_z; 1182 exponent_z = 0; 1183 } 1184 coefficient_z *= power10_table_128[scale].w[0]; 1185 1186#ifdef SET_STATUS_FLAGS 1187 if (round_flag) { 1188 __set_status_flags (fpsc, INEXACT_EXCEPTION); 1189 if (coefficient_z < 1000000000000000ull) 1190 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION); 1191 else if ((coefficient_z == 1000000000000000ull) && !exponent_z 1192 && ((SINT64) (round_dir ^ sign_z) < 0) && round_flag 1193 && (rmode == ROUNDING_DOWN || rmode == ROUNDING_TO_ZERO)) 1194 __set_status_flags (fpsc, UNDERFLOW_EXCEPTION); 1195 } 1196#endif 1197 1198#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 1199#ifndef IEEE_ROUND_NEAREST 1200 if (round_flag && (rmode & 3)) { 1201 D = round_dir ^ sign_z; 1202 1203 if (rmode == ROUNDING_UP) { 1204 if (D >= 0) 1205 coefficient_z++; 1206 } else { 1207 if (D < 0) 1208 coefficient_z--; 1209 if (coefficient_z < 1000000000000000ull && exponent_z) { 1210 coefficient_z = 9999999999999999ull; 1211 exponent_z--; 1212 } 1213 } 1214 } 1215#endif 1216#endif 1217 1218 return get_BID64 (sign_z, exponent_z, coefficient_z, rounding_mode, 1219 fpsc); 1220} 1221 1222 1223////////////////////////////////////////////////////////////////////////// 1224// 1225// 0*10^ey + cz*10^ez, ey<ez 1226// 1227////////////////////////////////////////////////////////////////////////// 1228 1229__BID_INLINE__ UINT64 1230add_zero64 (int exponent_y, UINT64 sign_z, int exponent_z, 1231 UINT64 coefficient_z, unsigned *prounding_mode, 1232 unsigned *fpsc) { 1233 int_double tempx; 1234 int bin_expon, scale_k, scale_cz; 1235 int diff_expon; 1236 1237 diff_expon = exponent_z - exponent_y; 1238 1239 tempx.d = (double) coefficient_z; 1240 bin_expon = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff; 1241 scale_cz = estimate_decimal_digits[bin_expon]; 1242 if (coefficient_z >= power10_table_128[scale_cz].w[0]) 1243 scale_cz++; 1244 1245 scale_k = 16 - scale_cz; 1246 if (diff_expon < scale_k) 1247 scale_k = diff_expon; 1248 coefficient_z *= power10_table_128[scale_k].w[0]; 1249 1250 return get_BID64 (sign_z, exponent_z - scale_k, coefficient_z, 1251 *prounding_mode, fpsc); 1252} 1253#endif 1254