1 /* Copyright (C) 2007-2019 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 * BID64 divide 26 ***************************************************************************** 27 * 28 * Algorithm description: 29 * 30 * if(coefficient_x<coefficient_y) 31 * p = number_digits(coefficient_y) - number_digits(coefficient_x) 32 * A = coefficient_x*10^p 33 * B = coefficient_y 34 * CA= A*10^(15+j), j=0 for A>=B, 1 otherwise 35 * Q = 0 36 * else 37 * get Q=(int)(coefficient_x/coefficient_y) 38 * (based on double precision divide) 39 * check for exact divide case 40 * Let R = coefficient_x - Q*coefficient_y 41 * Let m=16-number_digits(Q) 42 * CA=R*10^m, Q=Q*10^m 43 * B = coefficient_y 44 * endif 45 * if (CA<2^64) 46 * Q += CA/B (64-bit unsigned divide) 47 * else 48 * get final Q using double precision divide, followed by 3 integer 49 * iterations 50 * if exact result, eliminate trailing zeros 51 * check for underflow 52 * round coefficient to nearest 53 * 54 ****************************************************************************/ 55 56 #include "bid_internal.h" 57 #include "bid_div_macros.h" 58 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 59 #include <fenv.h> 60 61 #define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT 62 #endif 63 64 extern UINT32 convert_table[5][128][2]; 65 extern SINT8 factors[][2]; 66 extern UINT8 packed_10000_zeros[]; 67 68 69 #if DECIMAL_CALL_BY_REFERENCE 70 71 void 72 bid64_div (UINT64 * pres, UINT64 * px, 73 UINT64 * 74 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 75 _EXC_INFO_PARAM) { 76 UINT64 x, y; 77 #else 78 79 UINT64 80 bid64_div (UINT64 x, 81 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM 82 _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 83 #endif 84 UINT128 CA, CT; 85 UINT64 sign_x, sign_y, coefficient_x, coefficient_y, A, B, QX, PD; 86 UINT64 A2, Q, Q2, B2, B4, B5, R, T, DU, res; 87 UINT64 valid_x, valid_y; 88 SINT64 D; 89 int_double t_scale, tempq, temp_b; 90 int_float tempx, tempy; 91 double da, db, dq, da_h, da_l; 92 int exponent_x, exponent_y, bin_expon_cx; 93 int diff_expon, ed1, ed2, bin_index; 94 int rmode, amount; 95 int nzeros, i, j, k, d5; 96 UINT32 QX32, tdigit[3], digit, digit_h, digit_low; 97 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 98 fexcept_t binaryflags = 0; 99 #endif 100 101 #if DECIMAL_CALL_BY_REFERENCE 102 #if !DECIMAL_GLOBAL_ROUNDING 103 _IDEC_round rnd_mode = *prnd_mode; 104 #endif 105 x = *px; 106 y = *py; 107 #endif 108 109 valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x); 110 valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y); 111 112 // unpack arguments, check for NaN or Infinity 113 if (!valid_x) { 114 // x is Inf. or NaN 115 #ifdef SET_STATUS_FLAGS 116 if ((y & SNAN_MASK64) == SNAN_MASK64) // y is sNaN 117 __set_status_flags (pfpsf, INVALID_EXCEPTION); 118 #endif 119 120 // test if x is NaN 121 if ((x & NAN_MASK64) == NAN_MASK64) { 122 #ifdef SET_STATUS_FLAGS 123 if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN 124 __set_status_flags (pfpsf, INVALID_EXCEPTION); 125 #endif 126 BID_RETURN (coefficient_x & QUIET_MASK64); 127 } 128 // x is Infinity? 129 if ((x & INFINITY_MASK64) == INFINITY_MASK64) { 130 // check if y is Inf or NaN 131 if ((y & INFINITY_MASK64) == INFINITY_MASK64) { 132 // y==Inf, return NaN 133 if ((y & NAN_MASK64) == INFINITY_MASK64) { // Inf/Inf 134 #ifdef SET_STATUS_FLAGS 135 __set_status_flags (pfpsf, INVALID_EXCEPTION); 136 #endif 137 BID_RETURN (NAN_MASK64); 138 } 139 } else { 140 // otherwise return +/-Inf 141 BID_RETURN (((x ^ y) & 0x8000000000000000ull) | 142 INFINITY_MASK64); 143 } 144 } 145 // x==0 146 if (((y & INFINITY_MASK64) != INFINITY_MASK64) 147 && !(coefficient_y)) { 148 // y==0 , return NaN 149 #ifdef SET_STATUS_FLAGS 150 __set_status_flags (pfpsf, INVALID_EXCEPTION); 151 #endif 152 BID_RETURN (NAN_MASK64); 153 } 154 if (((y & INFINITY_MASK64) != INFINITY_MASK64)) { 155 if ((y & SPECIAL_ENCODING_MASK64) == SPECIAL_ENCODING_MASK64) 156 exponent_y = ((UINT32) (y >> 51)) & 0x3ff; 157 else 158 exponent_y = ((UINT32) (y >> 53)) & 0x3ff; 159 sign_y = y & 0x8000000000000000ull; 160 161 exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS; 162 if (exponent_x > DECIMAL_MAX_EXPON_64) 163 exponent_x = DECIMAL_MAX_EXPON_64; 164 else if (exponent_x < 0) 165 exponent_x = 0; 166 BID_RETURN ((sign_x ^ sign_y) | (((UINT64) exponent_x) << 53)); 167 } 168 169 } 170 if (!valid_y) { 171 // y is Inf. or NaN 172 173 // test if y is NaN 174 if ((y & NAN_MASK64) == NAN_MASK64) { 175 #ifdef SET_STATUS_FLAGS 176 if ((y & SNAN_MASK64) == SNAN_MASK64) // sNaN 177 __set_status_flags (pfpsf, INVALID_EXCEPTION); 178 #endif 179 BID_RETURN (coefficient_y & QUIET_MASK64); 180 } 181 // y is Infinity? 182 if ((y & INFINITY_MASK64) == INFINITY_MASK64) { 183 // return +/-0 184 BID_RETURN (((x ^ y) & 0x8000000000000000ull)); 185 } 186 // y is 0 187 #ifdef SET_STATUS_FLAGS 188 __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION); 189 #endif 190 BID_RETURN ((sign_x ^ sign_y) | INFINITY_MASK64); 191 } 192 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 193 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS); 194 #endif 195 diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS; 196 197 if (coefficient_x < coefficient_y) { 198 // get number of decimal digits for c_x, c_y 199 200 //--- get number of bits in the coefficients of x and y --- 201 tempx.d = (float) coefficient_x; 202 tempy.d = (float) coefficient_y; 203 bin_index = (tempy.i - tempx.i) >> 23; 204 205 A = coefficient_x * power10_index_binexp[bin_index]; 206 B = coefficient_y; 207 208 temp_b.d = (double) B; 209 210 // compare A, B 211 DU = (A - B) >> 63; 212 ed1 = 15 + (int) DU; 213 ed2 = estimate_decimal_digits[bin_index] + ed1; 214 T = power10_table_128[ed1].w[0]; 215 __mul_64x64_to_128 (CA, A, T); 216 217 Q = 0; 218 diff_expon = diff_expon - ed2; 219 220 // adjust double precision db, to ensure that later A/B - (int)(da/db) > -1 221 if (coefficient_y < 0x0020000000000000ull) { 222 temp_b.i += 1; 223 db = temp_b.d; 224 } else 225 db = (double) (B + 2 + (B & 1)); 226 227 } else { 228 // get c_x/c_y 229 230 // set last bit before conversion to DP 231 A2 = coefficient_x | 1; 232 da = (double) A2; 233 234 db = (double) coefficient_y; 235 236 tempq.d = da / db; 237 Q = (UINT64) tempq.d; 238 239 R = coefficient_x - coefficient_y * Q; 240 241 // will use to get number of dec. digits of Q 242 bin_expon_cx = (tempq.i >> 52) - 0x3ff; 243 244 // R<0 ? 245 D = ((SINT64) R) >> 63; 246 Q += D; 247 R += (coefficient_y & D); 248 249 // exact result ? 250 if (((SINT64) R) <= 0) { 251 // can have R==-1 for coeff_y==1 252 res = 253 get_BID64 (sign_x ^ sign_y, diff_expon, (Q + R), rnd_mode, 254 pfpsf); 255 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 256 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 257 #endif 258 BID_RETURN (res); 259 } 260 // get decimal digits of Q 261 DU = power10_index_binexp[bin_expon_cx] - Q - 1; 262 DU >>= 63; 263 264 ed2 = 16 - estimate_decimal_digits[bin_expon_cx] - (int) DU; 265 266 T = power10_table_128[ed2].w[0]; 267 __mul_64x64_to_128 (CA, R, T); 268 B = coefficient_y; 269 270 Q *= power10_table_128[ed2].w[0]; 271 diff_expon -= ed2; 272 273 } 274 275 if (!CA.w[1]) { 276 Q2 = CA.w[0] / B; 277 B2 = B + B; 278 B4 = B2 + B2; 279 R = CA.w[0] - Q2 * B; 280 Q += Q2; 281 } else { 282 283 // 2^64 284 t_scale.i = 0x43f0000000000000ull; 285 // convert CA to DP 286 da_h = CA.w[1]; 287 da_l = CA.w[0]; 288 da = da_h * t_scale.d + da_l; 289 290 // quotient 291 dq = da / db; 292 Q2 = (UINT64) dq; 293 294 // get w[0] remainder 295 R = CA.w[0] - Q2 * B; 296 297 // R<0 ? 298 D = ((SINT64) R) >> 63; 299 Q2 += D; 300 R += (B & D); 301 302 // now R<6*B 303 304 // quick divide 305 306 // 4*B 307 B2 = B + B; 308 B4 = B2 + B2; 309 310 R = R - B4; 311 // R<0 ? 312 D = ((SINT64) R) >> 63; 313 // restore R if negative 314 R += (B4 & D); 315 Q2 += ((~D) & 4); 316 317 R = R - B2; 318 // R<0 ? 319 D = ((SINT64) R) >> 63; 320 // restore R if negative 321 R += (B2 & D); 322 Q2 += ((~D) & 2); 323 324 R = R - B; 325 // R<0 ? 326 D = ((SINT64) R) >> 63; 327 // restore R if negative 328 R += (B & D); 329 Q2 += ((~D) & 1); 330 331 Q += Q2; 332 } 333 334 #ifdef SET_STATUS_FLAGS 335 if (R) { 336 // set status flags 337 __set_status_flags (pfpsf, INEXACT_EXCEPTION); 338 } 339 #ifndef LEAVE_TRAILING_ZEROS 340 else 341 #endif 342 #else 343 #ifndef LEAVE_TRAILING_ZEROS 344 if (!R) 345 #endif 346 #endif 347 #ifndef LEAVE_TRAILING_ZEROS 348 { 349 // eliminate trailing zeros 350 351 // check whether CX, CY are short 352 if ((coefficient_x <= 1024) && (coefficient_y <= 1024)) { 353 i = (int) coefficient_y - 1; 354 j = (int) coefficient_x - 1; 355 // difference in powers of 2 factors for Y and X 356 nzeros = ed2 - factors[i][0] + factors[j][0]; 357 // difference in powers of 5 factors 358 d5 = ed2 - factors[i][1] + factors[j][1]; 359 if (d5 < nzeros) 360 nzeros = d5; 361 362 __mul_64x64_to_128 (CT, Q, reciprocals10_64[nzeros]); 363 364 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128 365 amount = short_recip_scale[nzeros]; 366 Q = CT.w[1] >> amount; 367 368 diff_expon += nzeros; 369 } else { 370 tdigit[0] = Q & 0x3ffffff; 371 tdigit[1] = 0; 372 QX = Q >> 26; 373 QX32 = QX; 374 nzeros = 0; 375 376 for (j = 0; QX32; j++, QX32 >>= 7) { 377 k = (QX32 & 127); 378 tdigit[0] += convert_table[j][k][0]; 379 tdigit[1] += convert_table[j][k][1]; 380 if (tdigit[0] >= 100000000) { 381 tdigit[0] -= 100000000; 382 tdigit[1]++; 383 } 384 } 385 386 digit = tdigit[0]; 387 if (!digit && !tdigit[1]) 388 nzeros += 16; 389 else { 390 if (!digit) { 391 nzeros += 8; 392 digit = tdigit[1]; 393 } 394 // decompose digit 395 PD = (UINT64) digit *0x068DB8BBull; 396 digit_h = (UINT32) (PD >> 40); 397 digit_low = digit - digit_h * 10000; 398 399 if (!digit_low) 400 nzeros += 4; 401 else 402 digit_h = digit_low; 403 404 if (!(digit_h & 1)) 405 nzeros += 406 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >> 407 (digit_h & 7)); 408 } 409 410 if (nzeros) { 411 __mul_64x64_to_128 (CT, Q, reciprocals10_64[nzeros]); 412 413 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128 414 amount = short_recip_scale[nzeros]; 415 Q = CT.w[1] >> amount; 416 } 417 diff_expon += nzeros; 418 419 } 420 if (diff_expon >= 0) { 421 res = 422 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, Q, 423 rnd_mode, pfpsf); 424 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 425 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 426 #endif 427 BID_RETURN (res); 428 } 429 } 430 #endif 431 432 if (diff_expon >= 0) { 433 #ifdef IEEE_ROUND_NEAREST 434 // round to nearest code 435 // R*10 436 R += R; 437 R = (R << 2) + R; 438 B5 = B4 + B; 439 440 // compare 10*R to 5*B 441 R = B5 - R; 442 // correction for (R==0 && (Q&1)) 443 R -= (Q & 1); 444 // R<0 ? 445 D = ((UINT64) R) >> 63; 446 Q += D; 447 #else 448 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY 449 // round to nearest code 450 // R*10 451 R += R; 452 R = (R << 2) + R; 453 B5 = B4 + B; 454 455 // compare 10*R to 5*B 456 R = B5 - R; 457 // correction for (R==0 && (Q&1)) 458 R -= (Q & 1); 459 // R<0 ? 460 D = ((UINT64) R) >> 63; 461 Q += D; 462 #else 463 rmode = rnd_mode; 464 if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2) 465 rmode = 3 - rmode; 466 switch (rmode) { 467 case 0: // round to nearest code 468 case ROUNDING_TIES_AWAY: 469 // R*10 470 R += R; 471 R = (R << 2) + R; 472 B5 = B4 + B; 473 // compare 10*R to 5*B 474 R = B5 - R; 475 // correction for (R==0 && (Q&1)) 476 R -= ((Q | (rmode >> 2)) & 1); 477 // R<0 ? 478 D = ((UINT64) R) >> 63; 479 Q += D; 480 break; 481 case ROUNDING_DOWN: 482 case ROUNDING_TO_ZERO: 483 break; 484 default: // rounding up 485 Q++; 486 break; 487 } 488 #endif 489 #endif 490 491 res = 492 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, Q, rnd_mode, 493 pfpsf); 494 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 495 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 496 #endif 497 BID_RETURN (res); 498 } else { 499 // UF occurs 500 501 #ifdef SET_STATUS_FLAGS 502 if ((diff_expon + 16 < 0)) { 503 // set status flags 504 __set_status_flags (pfpsf, INEXACT_EXCEPTION); 505 } 506 #endif 507 rmode = rnd_mode; 508 res = 509 get_BID64_UF (sign_x ^ sign_y, diff_expon, Q, R, rmode, pfpsf); 510 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 511 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 512 #endif 513 BID_RETURN (res); 514 515 } 516 } 517 518 519 520 TYPE0_FUNCTION_ARGTYPE1_ARG128 (UINT64, bid64dq_div, UINT64, x, y) 521 UINT256 CA4 = 522 { {0x0ull, 0x0ull, 0x0ull, 0x0ull} }, CA4r, P256, QB256; 523 UINT128 CX, CY, T128, CQ, CQ2, CR, CA, TP128, Qh, Tmp; 524 UINT64 sign_x, sign_y, T, carry64, D, Q_low, QX, valid_y, PD, res; 525 int_float fx, fy, f64; 526 UINT32 QX32, tdigit[3], digit, digit_h, digit_low; 527 int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2, 528 digits_q, amount; 529 int nzeros, i, j, k, d5, done = 0; 530 unsigned rmode; 531 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 532 fexcept_t binaryflags = 0; 533 #endif 534 535 valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y); 536 537 // unpack arguments, check for NaN or Infinity 538 CX.w[1] = 0; 539 if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], (x))) { 540 #ifdef SET_STATUS_FLAGS 541 if (((y.w[1] & SNAN_MASK64) == SNAN_MASK64) || // y is sNaN 542 ((x & SNAN_MASK64) == SNAN_MASK64)) 543 __set_status_flags (pfpsf, INVALID_EXCEPTION); 544 #endif 545 // test if x is NaN 546 if (((x) & 0x7c00000000000000ull) == 0x7c00000000000000ull) { 547 res = CX.w[0]; 548 BID_RETURN (res & QUIET_MASK64); 549 } 550 // x is Infinity? 551 if (((x) & 0x7800000000000000ull) == 0x7800000000000000ull) { 552 // check if y is Inf. 553 if (((y.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull)) 554 // return NaN 555 { 556 #ifdef SET_STATUS_FLAGS 557 __set_status_flags (pfpsf, INVALID_EXCEPTION); 558 #endif 559 res = 0x7c00000000000000ull; 560 BID_RETURN (res); 561 } 562 if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull)) { 563 // otherwise return +/-Inf 564 res = 565 (((x) ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull; 566 BID_RETURN (res); 567 } 568 } 569 // x is 0 570 if ((y.w[1] & INFINITY_MASK64) != INFINITY_MASK64) { 571 if ((!CY.w[0]) && !(CY.w[1] & 0x0001ffffffffffffull)) { 572 #ifdef SET_STATUS_FLAGS 573 __set_status_flags (pfpsf, INVALID_EXCEPTION); 574 #endif 575 // x=y=0, return NaN 576 res = 0x7c00000000000000ull; 577 BID_RETURN (res); 578 } 579 // return 0 580 res = ((x) ^ y.w[1]) & 0x8000000000000000ull; 581 exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128; 582 if (exponent_x > DECIMAL_MAX_EXPON_64) 583 exponent_x = DECIMAL_MAX_EXPON_64; 584 else if (exponent_x < 0) 585 exponent_x = 0; 586 res |= (((UINT64) exponent_x) << 53); 587 BID_RETURN (res); 588 } 589 } 590 exponent_x += (DECIMAL_EXPONENT_BIAS_128 - DECIMAL_EXPONENT_BIAS); 591 if (!valid_y) { 592 // y is Inf. or NaN 593 594 // test if y is NaN 595 if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) { 596 #ifdef SET_STATUS_FLAGS 597 if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN 598 __set_status_flags (pfpsf, INVALID_EXCEPTION); 599 #endif 600 Tmp.w[1] = (CY.w[1] & 0x00003fffffffffffull); 601 Tmp.w[0] = CY.w[0]; 602 TP128 = reciprocals10_128[18]; 603 __mul_128x128_high (Qh, Tmp, TP128); 604 amount = recip_scale[18]; 605 __shr_128 (Tmp, Qh, amount); 606 res = (CY.w[1] & 0xfc00000000000000ull) | Tmp.w[0]; 607 BID_RETURN (res); 608 } 609 // y is Infinity? 610 if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) { 611 // return +/-0 612 res = sign_x ^ sign_y; 613 BID_RETURN (res); 614 } 615 // y is 0, return +/-Inf 616 res = 617 (((x) ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull; 618 #ifdef SET_STATUS_FLAGS 619 __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION); 620 #endif 621 BID_RETURN (res); 622 } 623 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 624 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS); 625 #endif 626 diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS; 627 628 if (__unsigned_compare_gt_128 (CY, CX)) { 629 // CX < CY 630 631 // 2^64 632 f64.i = 0x5f800000; 633 634 // fx ~ CX, fy ~ CY 635 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0]; 636 fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0]; 637 // expon_cy - expon_cx 638 bin_index = (fy.i - fx.i) >> 23; 639 640 if (CX.w[1]) { 641 T = power10_index_binexp_128[bin_index].w[0]; 642 __mul_64x128_short (CA, T, CX); 643 } else { 644 T128 = power10_index_binexp_128[bin_index]; 645 __mul_64x128_short (CA, CX.w[0], T128); 646 } 647 648 ed2 = 15; 649 if (__unsigned_compare_gt_128 (CY, CA)) 650 ed2++; 651 652 T128 = power10_table_128[ed2]; 653 __mul_128x128_to_256 (CA4, CA, T128); 654 655 ed2 += estimate_decimal_digits[bin_index]; 656 CQ.w[0] = CQ.w[1] = 0; 657 diff_expon = diff_expon - ed2; 658 659 } else { 660 // get CQ = CX/CY 661 __div_128_by_128 (&CQ, &CR, CX, CY); 662 663 // get number of decimal digits in CQ 664 // 2^64 665 f64.i = 0x5f800000; 666 fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0]; 667 // binary expon. of CQ 668 bin_expon = (fx.i - 0x3f800000) >> 23; 669 670 digits_q = estimate_decimal_digits[bin_expon]; 671 TP128.w[0] = power10_index_binexp_128[bin_expon].w[0]; 672 TP128.w[1] = power10_index_binexp_128[bin_expon].w[1]; 673 if (__unsigned_compare_ge_128 (CQ, TP128)) 674 digits_q++; 675 676 if (digits_q <= 16) { 677 if (!CR.w[1] && !CR.w[0]) { 678 res = get_BID64 (sign_x ^ sign_y, diff_expon, 679 CQ.w[0], rnd_mode, pfpsf); 680 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 681 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 682 #endif 683 BID_RETURN (res); 684 } 685 686 ed2 = 16 - digits_q; 687 T128.w[0] = power10_table_128[ed2].w[0]; 688 __mul_64x128_to_192 (CA4, (T128.w[0]), CR); 689 diff_expon = diff_expon - ed2; 690 CQ.w[0] *= T128.w[0]; 691 } else { 692 ed2 = digits_q - 16; 693 diff_expon += ed2; 694 T128 = reciprocals10_128[ed2]; 695 __mul_128x128_to_256 (P256, CQ, T128); 696 amount = recip_scale[ed2]; 697 CQ.w[0] = (P256.w[2] >> amount) | (P256.w[3] << (64 - amount)); 698 CQ.w[1] = 0; 699 700 __mul_64x64_to_128 (CQ2, CQ.w[0], (power10_table_128[ed2].w[0])); 701 702 __mul_64x64_to_128 (QB256, CQ2.w[0], CY.w[0]); 703 QB256.w[1] += CQ2.w[0] * CY.w[1] + CQ2.w[1] * CY.w[0]; 704 705 CA4.w[1] = CX.w[1] - QB256.w[1]; 706 CA4.w[0] = CX.w[0] - QB256.w[0]; 707 if (CX.w[0] < QB256.w[0]) 708 CA4.w[1]--; 709 if (CR.w[0] || CR.w[1]) 710 CA4.w[0] |= 1; 711 done = 1; 712 713 } 714 715 } 716 if (!done) { 717 __div_256_by_128 (&CQ, &CA4, CY); 718 } 719 720 721 722 #ifdef SET_STATUS_FLAGS 723 if (CA4.w[0] || CA4.w[1]) { 724 // set status flags 725 __set_status_flags (pfpsf, INEXACT_EXCEPTION); 726 } 727 #ifndef LEAVE_TRAILING_ZEROS 728 else 729 #endif 730 #else 731 #ifndef LEAVE_TRAILING_ZEROS 732 if (!CA4.w[0] && !CA4.w[1]) 733 #endif 734 #endif 735 #ifndef LEAVE_TRAILING_ZEROS 736 // check whether result is exact 737 { 738 // check whether CX, CY are short 739 if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) { 740 i = (int) CY.w[0] - 1; 741 j = (int) CX.w[0] - 1; 742 // difference in powers of 2 factors for Y and X 743 nzeros = ed2 - factors[i][0] + factors[j][0]; 744 // difference in powers of 5 factors 745 d5 = ed2 - factors[i][1] + factors[j][1]; 746 if (d5 < nzeros) 747 nzeros = d5; 748 // get P*(2^M[extra_digits])/10^extra_digits 749 __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]); 750 751 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 752 amount = recip_scale[nzeros]; 753 __shr_128_long (CQ, Qh, amount); 754 755 diff_expon += nzeros; 756 } else { 757 // decompose Q as Qh*10^17 + Ql 758 Q_low = CQ.w[0]; 759 760 { 761 tdigit[0] = Q_low & 0x3ffffff; 762 tdigit[1] = 0; 763 QX = Q_low >> 26; 764 QX32 = QX; 765 nzeros = 0; 766 767 for (j = 0; QX32; j++, QX32 >>= 7) { 768 k = (QX32 & 127); 769 tdigit[0] += convert_table[j][k][0]; 770 tdigit[1] += convert_table[j][k][1]; 771 if (tdigit[0] >= 100000000) { 772 tdigit[0] -= 100000000; 773 tdigit[1]++; 774 } 775 } 776 777 if (tdigit[1] >= 100000000) { 778 tdigit[1] -= 100000000; 779 if (tdigit[1] >= 100000000) 780 tdigit[1] -= 100000000; 781 } 782 783 digit = tdigit[0]; 784 if (!digit && !tdigit[1]) 785 nzeros += 16; 786 else { 787 if (!digit) { 788 nzeros += 8; 789 digit = tdigit[1]; 790 } 791 // decompose digit 792 PD = (UINT64) digit *0x068DB8BBull; 793 digit_h = (UINT32) (PD >> 40); 794 digit_low = digit - digit_h * 10000; 795 796 if (!digit_low) 797 nzeros += 4; 798 else 799 digit_h = digit_low; 800 801 if (!(digit_h & 1)) 802 nzeros += 803 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >> 804 (digit_h & 7)); 805 } 806 807 if (nzeros) { 808 // get P*(2^M[extra_digits])/10^extra_digits 809 __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]); 810 811 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 812 amount = recip_scale[nzeros]; 813 __shr_128 (CQ, Qh, amount); 814 } 815 diff_expon += nzeros; 816 817 } 818 } 819 if(diff_expon>=0){ 820 res = 821 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0], 822 rnd_mode, pfpsf); 823 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 824 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 825 #endif 826 BID_RETURN (res); 827 } 828 } 829 #endif 830 831 if (diff_expon >= 0) { 832 #ifdef IEEE_ROUND_NEAREST 833 // rounding 834 // 2*CA4 - CY 835 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63); 836 CA4r.w[0] = CA4.w[0] + CA4.w[0]; 837 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]); 838 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64; 839 840 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0; 841 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D); 842 843 CQ.w[0] += carry64; 844 #else 845 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY 846 // rounding 847 // 2*CA4 - CY 848 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63); 849 CA4r.w[0] = CA4.w[0] + CA4.w[0]; 850 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]); 851 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64; 852 853 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1; 854 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D; 855 856 CQ.w[0] += carry64; 857 if (CQ.w[0] < carry64) 858 CQ.w[1]++; 859 #else 860 rmode = rnd_mode; 861 if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2) 862 rmode = 3 - rmode; 863 switch (rmode) { 864 case ROUNDING_TO_NEAREST: // round to nearest code 865 // rounding 866 // 2*CA4 - CY 867 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63); 868 CA4r.w[0] = CA4.w[0] + CA4.w[0]; 869 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]); 870 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64; 871 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0; 872 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D); 873 CQ.w[0] += carry64; 874 if (CQ.w[0] < carry64) 875 CQ.w[1]++; 876 break; 877 case ROUNDING_TIES_AWAY: 878 // rounding 879 // 2*CA4 - CY 880 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63); 881 CA4r.w[0] = CA4.w[0] + CA4.w[0]; 882 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]); 883 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64; 884 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1; 885 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D; 886 CQ.w[0] += carry64; 887 if (CQ.w[0] < carry64) 888 CQ.w[1]++; 889 break; 890 case ROUNDING_DOWN: 891 case ROUNDING_TO_ZERO: 892 break; 893 default: // rounding up 894 CQ.w[0]++; 895 if (!CQ.w[0]) 896 CQ.w[1]++; 897 break; 898 } 899 #endif 900 #endif 901 902 res = 903 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0], rnd_mode, 904 pfpsf); 905 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 906 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 907 #endif 908 BID_RETURN (res); 909 } else { 910 // UF occurs 911 912 #ifdef SET_STATUS_FLAGS 913 if ((diff_expon + 16 < 0)) { 914 // set status flags 915 __set_status_flags (pfpsf, INEXACT_EXCEPTION); 916 } 917 #endif 918 rmode = rnd_mode; 919 res = 920 get_BID64_UF (sign_x ^ sign_y, diff_expon, CQ.w[0], CA4.w[1] | CA4.w[0], rmode, pfpsf); 921 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 922 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 923 #endif 924 BID_RETURN (res); 925 926 } 927 928 } 929 930 931 //#define LEAVE_TRAILING_ZEROS 932 933 TYPE0_FUNCTION_ARG128_ARGTYPE2 (UINT64, bid64qd_div, x, UINT64, y) 934 935 UINT256 CA4 = 936 { {0x0ull, 0x0ull, 0x0ull, 0x0ull} }, CA4r, P256, QB256; 937 UINT128 CX, CY, T128, CQ, CQ2, CR, CA, TP128, Qh, Tmp; 938 UINT64 sign_x, sign_y, T, carry64, D, Q_low, QX, PD, res, valid_y; 939 int_float fx, fy, f64; 940 UINT32 QX32, tdigit[3], digit, digit_h, digit_low; 941 int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2, 942 digits_q, amount; 943 int nzeros, i, j, k, d5, done = 0; 944 unsigned rmode; 945 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 946 fexcept_t binaryflags = 0; 947 #endif 948 949 valid_y = unpack_BID64 (&sign_y, &exponent_y, &CY.w[0], (y)); 950 951 // unpack arguments, check for NaN or Infinity 952 if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) { 953 // test if x is NaN 954 if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) { 955 #ifdef SET_STATUS_FLAGS 956 if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull || // sNaN 957 (y & 0x7e00000000000000ull) == 0x7e00000000000000ull) 958 __set_status_flags (pfpsf, INVALID_EXCEPTION); 959 #endif 960 Tmp.w[1] = (CX.w[1] & 0x00003fffffffffffull); 961 Tmp.w[0] = CX.w[0]; 962 TP128 = reciprocals10_128[18]; 963 __mul_128x128_high (Qh, Tmp, TP128); 964 amount = recip_scale[18]; 965 __shr_128 (Tmp, Qh, amount); 966 res = (CX.w[1] & 0xfc00000000000000ull) | Tmp.w[0]; 967 BID_RETURN (res); 968 } 969 // x is Infinity? 970 if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) { 971 // check if y is Inf. 972 if (((y & 0x7c00000000000000ull) == 0x7800000000000000ull)) 973 // return NaN 974 { 975 #ifdef SET_STATUS_FLAGS 976 __set_status_flags (pfpsf, INVALID_EXCEPTION); 977 #endif 978 res = 0x7c00000000000000ull; 979 BID_RETURN (res); 980 } 981 if (((y & 0x7c00000000000000ull) != 0x7c00000000000000ull)) { 982 // otherwise return +/-Inf 983 res = 984 ((x.w[1] ^ (y)) & 0x8000000000000000ull) | 0x7800000000000000ull; 985 BID_RETURN (res); 986 } 987 } 988 // x is 0 989 if (((y & INFINITY_MASK64) != INFINITY_MASK64) && 990 !(CY.w[0])) { 991 #ifdef SET_STATUS_FLAGS 992 __set_status_flags (pfpsf, INVALID_EXCEPTION); 993 #endif 994 // x=y=0, return NaN 995 res = 0x7c00000000000000ull; 996 BID_RETURN (res); 997 } 998 // return 0 999 if (((y & 0x7800000000000000ull) != 0x7800000000000000ull)) { 1000 if (!CY.w[0]) { 1001 #ifdef SET_STATUS_FLAGS 1002 __set_status_flags (pfpsf, INVALID_EXCEPTION); 1003 #endif 1004 res = 0x7c00000000000000ull; 1005 BID_RETURN (res); 1006 } 1007 exponent_x = 1008 exponent_x - exponent_y - DECIMAL_EXPONENT_BIAS_128 + 1009 (DECIMAL_EXPONENT_BIAS << 1); 1010 if (exponent_x > DECIMAL_MAX_EXPON_64) 1011 exponent_x = DECIMAL_MAX_EXPON_64; 1012 else if (exponent_x < 0) 1013 exponent_x = 0; 1014 res = (sign_x ^ sign_y) | (((UINT64) exponent_x) << 53); 1015 BID_RETURN (res); 1016 } 1017 } 1018 CY.w[1] = 0; 1019 if (!valid_y) { 1020 // y is Inf. or NaN 1021 1022 // test if y is NaN 1023 if ((y & NAN_MASK64) == NAN_MASK64) { 1024 #ifdef SET_STATUS_FLAGS 1025 if ((y & SNAN_MASK64) == SNAN_MASK64) // sNaN 1026 __set_status_flags (pfpsf, INVALID_EXCEPTION); 1027 #endif 1028 BID_RETURN (CY.w[0] & QUIET_MASK64); 1029 } 1030 // y is Infinity? 1031 if (((y) & 0x7800000000000000ull) == 0x7800000000000000ull) { 1032 // return +/-0 1033 res = sign_x ^ sign_y; 1034 BID_RETURN (res); 1035 } 1036 // y is 0, return +/-Inf 1037 res = 1038 ((x.w[1] ^ (y)) & 0x8000000000000000ull) | 0x7800000000000000ull; 1039 #ifdef SET_STATUS_FLAGS 1040 __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION); 1041 #endif 1042 BID_RETURN (res); 1043 } 1044 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 1045 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS); 1046 #endif 1047 diff_expon = 1048 exponent_x - exponent_y - DECIMAL_EXPONENT_BIAS_128 + 1049 (DECIMAL_EXPONENT_BIAS << 1); 1050 1051 if (__unsigned_compare_gt_128 (CY, CX)) { 1052 // CX < CY 1053 1054 // 2^64 1055 f64.i = 0x5f800000; 1056 1057 // fx ~ CX, fy ~ CY 1058 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0]; 1059 fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0]; 1060 // expon_cy - expon_cx 1061 bin_index = (fy.i - fx.i) >> 23; 1062 1063 if (CX.w[1]) { 1064 T = power10_index_binexp_128[bin_index].w[0]; 1065 __mul_64x128_short (CA, T, CX); 1066 } else { 1067 T128 = power10_index_binexp_128[bin_index]; 1068 __mul_64x128_short (CA, CX.w[0], T128); 1069 } 1070 1071 ed2 = 15; 1072 if (__unsigned_compare_gt_128 (CY, CA)) 1073 ed2++; 1074 1075 T128 = power10_table_128[ed2]; 1076 __mul_128x128_to_256 (CA4, CA, T128); 1077 1078 ed2 += estimate_decimal_digits[bin_index]; 1079 CQ.w[0] = CQ.w[1] = 0; 1080 diff_expon = diff_expon - ed2; 1081 1082 } else { 1083 // get CQ = CX/CY 1084 __div_128_by_128 (&CQ, &CR, CX, CY); 1085 1086 // get number of decimal digits in CQ 1087 // 2^64 1088 f64.i = 0x5f800000; 1089 fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0]; 1090 // binary expon. of CQ 1091 bin_expon = (fx.i - 0x3f800000) >> 23; 1092 1093 digits_q = estimate_decimal_digits[bin_expon]; 1094 TP128.w[0] = power10_index_binexp_128[bin_expon].w[0]; 1095 TP128.w[1] = power10_index_binexp_128[bin_expon].w[1]; 1096 if (__unsigned_compare_ge_128 (CQ, TP128)) 1097 digits_q++; 1098 1099 if (digits_q <= 16) { 1100 if (!CR.w[1] && !CR.w[0]) { 1101 res = get_BID64 (sign_x ^ sign_y, diff_expon, 1102 CQ.w[0], rnd_mode, pfpsf); 1103 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 1104 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 1105 #endif 1106 BID_RETURN (res); 1107 } 1108 1109 ed2 = 16 - digits_q; 1110 T128.w[0] = power10_table_128[ed2].w[0]; 1111 __mul_64x128_to_192 (CA4, (T128.w[0]), CR); 1112 diff_expon = diff_expon - ed2; 1113 CQ.w[0] *= T128.w[0]; 1114 } else { 1115 ed2 = digits_q - 16; 1116 diff_expon += ed2; 1117 T128 = reciprocals10_128[ed2]; 1118 __mul_128x128_to_256 (P256, CQ, T128); 1119 amount = recip_scale[ed2]; 1120 CQ.w[0] = (P256.w[2] >> amount) | (P256.w[3] << (64 - amount)); 1121 CQ.w[1] = 0; 1122 1123 __mul_64x64_to_128 (CQ2, CQ.w[0], (power10_table_128[ed2].w[0])); 1124 1125 __mul_64x64_to_128 (QB256, CQ2.w[0], CY.w[0]); 1126 QB256.w[1] += CQ2.w[0] * CY.w[1] + CQ2.w[1] * CY.w[0]; 1127 1128 CA4.w[1] = CX.w[1] - QB256.w[1]; 1129 CA4.w[0] = CX.w[0] - QB256.w[0]; 1130 if (CX.w[0] < QB256.w[0]) 1131 CA4.w[1]--; 1132 if (CR.w[0] || CR.w[1]) 1133 CA4.w[0] |= 1; 1134 done = 1; 1135 if(CA4.w[1]|CA4.w[0]) { 1136 __mul_64x128_low(CY, (power10_table_128[ed2].w[0]),CY); 1137 } 1138 1139 } 1140 1141 } 1142 1143 if (!done) { 1144 __div_256_by_128 (&CQ, &CA4, CY); 1145 } 1146 1147 #ifdef SET_STATUS_FLAGS 1148 if (CA4.w[0] || CA4.w[1]) { 1149 // set status flags 1150 __set_status_flags (pfpsf, INEXACT_EXCEPTION); 1151 } 1152 #ifndef LEAVE_TRAILING_ZEROS 1153 else 1154 #endif 1155 #else 1156 #ifndef LEAVE_TRAILING_ZEROS 1157 if (!CA4.w[0] && !CA4.w[1]) 1158 #endif 1159 #endif 1160 #ifndef LEAVE_TRAILING_ZEROS 1161 // check whether result is exact 1162 { 1163 if(!done) { 1164 // check whether CX, CY are short 1165 if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) { 1166 i = (int) CY.w[0] - 1; 1167 j = (int) CX.w[0] - 1; 1168 // difference in powers of 2 factors for Y and X 1169 nzeros = ed2 - factors[i][0] + factors[j][0]; 1170 // difference in powers of 5 factors 1171 d5 = ed2 - factors[i][1] + factors[j][1]; 1172 if (d5 < nzeros) 1173 nzeros = d5; 1174 // get P*(2^M[extra_digits])/10^extra_digits 1175 __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]); 1176 //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2]; 1177 1178 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 1179 amount = recip_scale[nzeros]; 1180 __shr_128_long (CQ, Qh, amount); 1181 1182 diff_expon += nzeros; 1183 } else { 1184 // decompose Q as Qh*10^17 + Ql 1185 //T128 = reciprocals10_128[17]; 1186 Q_low = CQ.w[0]; 1187 1188 { 1189 tdigit[0] = Q_low & 0x3ffffff; 1190 tdigit[1] = 0; 1191 QX = Q_low >> 26; 1192 QX32 = QX; 1193 nzeros = 0; 1194 1195 for (j = 0; QX32; j++, QX32 >>= 7) { 1196 k = (QX32 & 127); 1197 tdigit[0] += convert_table[j][k][0]; 1198 tdigit[1] += convert_table[j][k][1]; 1199 if (tdigit[0] >= 100000000) { 1200 tdigit[0] -= 100000000; 1201 tdigit[1]++; 1202 } 1203 } 1204 1205 if (tdigit[1] >= 100000000) { 1206 tdigit[1] -= 100000000; 1207 if (tdigit[1] >= 100000000) 1208 tdigit[1] -= 100000000; 1209 } 1210 1211 digit = tdigit[0]; 1212 if (!digit && !tdigit[1]) 1213 nzeros += 16; 1214 else { 1215 if (!digit) { 1216 nzeros += 8; 1217 digit = tdigit[1]; 1218 } 1219 // decompose digit 1220 PD = (UINT64) digit *0x068DB8BBull; 1221 digit_h = (UINT32) (PD >> 40); 1222 digit_low = digit - digit_h * 10000; 1223 1224 if (!digit_low) 1225 nzeros += 4; 1226 else 1227 digit_h = digit_low; 1228 1229 if (!(digit_h & 1)) 1230 nzeros += 1231 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >> 1232 (digit_h & 7)); 1233 } 1234 1235 if (nzeros) { 1236 // get P*(2^M[extra_digits])/10^extra_digits 1237 __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]); 1238 1239 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 1240 amount = recip_scale[nzeros]; 1241 __shr_128 (CQ, Qh, amount); 1242 } 1243 diff_expon += nzeros; 1244 1245 } 1246 } 1247 } 1248 if(diff_expon>=0){ 1249 res = 1250 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0], 1251 rnd_mode, pfpsf); 1252 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 1253 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 1254 #endif 1255 BID_RETURN (res); 1256 } 1257 } 1258 #endif 1259 1260 if (diff_expon >= 0) { 1261 #ifdef IEEE_ROUND_NEAREST 1262 // rounding 1263 // 2*CA4 - CY 1264 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63); 1265 CA4r.w[0] = CA4.w[0] + CA4.w[0]; 1266 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]); 1267 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64; 1268 1269 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0; 1270 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D); 1271 1272 CQ.w[0] += carry64; 1273 //if(CQ.w[0]<carry64) 1274 //CQ.w[1] ++; 1275 #else 1276 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY 1277 // rounding 1278 // 2*CA4 - CY 1279 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63); 1280 CA4r.w[0] = CA4.w[0] + CA4.w[0]; 1281 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]); 1282 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64; 1283 1284 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1; 1285 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D; 1286 1287 CQ.w[0] += carry64; 1288 if (CQ.w[0] < carry64) 1289 CQ.w[1]++; 1290 #else 1291 rmode = rnd_mode; 1292 if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2) 1293 rmode = 3 - rmode; 1294 switch (rmode) { 1295 case ROUNDING_TO_NEAREST: // round to nearest code 1296 // rounding 1297 // 2*CA4 - CY 1298 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63); 1299 CA4r.w[0] = CA4.w[0] + CA4.w[0]; 1300 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]); 1301 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64; 1302 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0; 1303 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D); 1304 CQ.w[0] += carry64; 1305 if (CQ.w[0] < carry64) 1306 CQ.w[1]++; 1307 break; 1308 case ROUNDING_TIES_AWAY: 1309 // rounding 1310 // 2*CA4 - CY 1311 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63); 1312 CA4r.w[0] = CA4.w[0] + CA4.w[0]; 1313 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]); 1314 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64; 1315 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1; 1316 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D; 1317 CQ.w[0] += carry64; 1318 if (CQ.w[0] < carry64) 1319 CQ.w[1]++; 1320 break; 1321 case ROUNDING_DOWN: 1322 case ROUNDING_TO_ZERO: 1323 break; 1324 default: // rounding up 1325 CQ.w[0]++; 1326 if (!CQ.w[0]) 1327 CQ.w[1]++; 1328 break; 1329 } 1330 #endif 1331 #endif 1332 1333 1334 res = 1335 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0], rnd_mode, 1336 pfpsf); 1337 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 1338 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 1339 #endif 1340 BID_RETURN (res); 1341 } else { 1342 // UF occurs 1343 1344 #ifdef SET_STATUS_FLAGS 1345 if ((diff_expon + 16 < 0)) { 1346 // set status flags 1347 __set_status_flags (pfpsf, INEXACT_EXCEPTION); 1348 } 1349 #endif 1350 rmode = rnd_mode; 1351 res = 1352 get_BID64_UF (sign_x ^ sign_y, diff_expon, CQ.w[0], CA4.w[1] | CA4.w[0], rmode, pfpsf); 1353 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 1354 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 1355 #endif 1356 BID_RETURN (res); 1357 1358 } 1359 1360 } 1361 1362 //#define LEAVE_TRAILING_ZEROS 1363 1364 extern UINT32 convert_table[5][128][2]; 1365 extern SINT8 factors[][2]; 1366 extern UINT8 packed_10000_zeros[]; 1367 1368 1369 //UINT64* bid64_div128x128(UINT64 res, UINT128 *px, UINT128 *py, unsigned rnd_mode, unsigned *pfpsf) 1370 1371 TYPE0_FUNCTION_ARG128_ARG128 (UINT64, bid64qq_div, x, y) 1372 UINT256 CA4 = 1373 { {0x0ull, 0x0ull, 0x0ull, 0x0ull} }, CA4r, P256, QB256; 1374 UINT128 CX, CY, T128, CQ, CQ2, CR, CA, TP128, Qh, Tmp; 1375 UINT64 sign_x, sign_y, T, carry64, D, Q_low, QX, valid_y, PD, res; 1376 int_float fx, fy, f64; 1377 UINT32 QX32, tdigit[3], digit, digit_h, digit_low; 1378 int exponent_x, exponent_y, bin_index, bin_expon, diff_expon, ed2, 1379 digits_q, amount; 1380 int nzeros, i, j, k, d5, done = 0; 1381 unsigned rmode; 1382 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 1383 fexcept_t binaryflags = 0; 1384 #endif 1385 1386 valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y); 1387 1388 // unpack arguments, check for NaN or Infinity 1389 if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) { 1390 // test if x is NaN 1391 if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) { 1392 #ifdef SET_STATUS_FLAGS 1393 if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull || // sNaN 1394 (y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) 1395 __set_status_flags (pfpsf, INVALID_EXCEPTION); 1396 #endif 1397 Tmp.w[1] = (CX.w[1] & 0x00003fffffffffffull); 1398 Tmp.w[0] = CX.w[0]; 1399 TP128 = reciprocals10_128[18]; 1400 __mul_128x128_high (Qh, Tmp, TP128); 1401 amount = recip_scale[18]; 1402 __shr_128 (Tmp, Qh, amount); 1403 res = (CX.w[1] & 0xfc00000000000000ull) | Tmp.w[0]; 1404 BID_RETURN (res); 1405 } 1406 // x is Infinity? 1407 if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) { 1408 // check if y is Inf. 1409 if (((y.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull)) 1410 // return NaN 1411 { 1412 #ifdef SET_STATUS_FLAGS 1413 __set_status_flags (pfpsf, INVALID_EXCEPTION); 1414 #endif 1415 res = 0x7c00000000000000ull; 1416 BID_RETURN (res); 1417 } 1418 if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull)) { 1419 // otherwise return +/-Inf 1420 res = 1421 ((x.w[1] ^ y. 1422 w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull; 1423 BID_RETURN (res); 1424 } 1425 } 1426 // x is 0 1427 if (((y.w[1] & 0x7800000000000000ull) != 0x7800000000000000ull)) { 1428 if ((!CY.w[0]) && !(CY.w[1] & 0x0001ffffffffffffull)) { 1429 #ifdef SET_STATUS_FLAGS 1430 __set_status_flags (pfpsf, INVALID_EXCEPTION); 1431 #endif 1432 // x=y=0, return NaN 1433 res = 0x7c00000000000000ull; 1434 BID_RETURN (res); 1435 } 1436 // return 0 1437 res = (x.w[1] ^ y.w[1]) & 0x8000000000000000ull; 1438 exponent_x = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS; 1439 if (exponent_x > DECIMAL_MAX_EXPON_64) 1440 exponent_x = DECIMAL_MAX_EXPON_64; 1441 else if (exponent_x < 0) 1442 exponent_x = 0; 1443 res |= (((UINT64) exponent_x) << 53); 1444 BID_RETURN (res); 1445 } 1446 } 1447 if (!valid_y) { 1448 // y is Inf. or NaN 1449 1450 // test if y is NaN 1451 if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) { 1452 #ifdef SET_STATUS_FLAGS 1453 if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN 1454 __set_status_flags (pfpsf, INVALID_EXCEPTION); 1455 #endif 1456 Tmp.w[1] = (CY.w[1] & 0x00003fffffffffffull); 1457 Tmp.w[0] = CY.w[0]; 1458 TP128 = reciprocals10_128[18]; 1459 __mul_128x128_high (Qh, Tmp, TP128); 1460 amount = recip_scale[18]; 1461 __shr_128 (Tmp, Qh, amount); 1462 res = (CY.w[1] & 0xfc00000000000000ull) | Tmp.w[0]; 1463 BID_RETURN (res); 1464 } 1465 // y is Infinity? 1466 if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) { 1467 // return +/-0 1468 res = sign_x ^ sign_y; 1469 BID_RETURN (res); 1470 } 1471 // y is 0, return +/-Inf 1472 res = 1473 ((x.w[1] ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull; 1474 #ifdef SET_STATUS_FLAGS 1475 __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION); 1476 #endif 1477 BID_RETURN (res); 1478 } 1479 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 1480 (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS); 1481 #endif 1482 diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS; 1483 1484 if (__unsigned_compare_gt_128 (CY, CX)) { 1485 // CX < CY 1486 1487 // 2^64 1488 f64.i = 0x5f800000; 1489 1490 // fx ~ CX, fy ~ CY 1491 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0]; 1492 fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0]; 1493 // expon_cy - expon_cx 1494 bin_index = (fy.i - fx.i) >> 23; 1495 1496 if (CX.w[1]) { 1497 T = power10_index_binexp_128[bin_index].w[0]; 1498 __mul_64x128_short (CA, T, CX); 1499 } else { 1500 T128 = power10_index_binexp_128[bin_index]; 1501 __mul_64x128_short (CA, CX.w[0], T128); 1502 } 1503 1504 ed2 = 15; 1505 if (__unsigned_compare_gt_128 (CY, CA)) 1506 ed2++; 1507 1508 T128 = power10_table_128[ed2]; 1509 __mul_128x128_to_256 (CA4, CA, T128); 1510 1511 ed2 += estimate_decimal_digits[bin_index]; 1512 CQ.w[0] = CQ.w[1] = 0; 1513 diff_expon = diff_expon - ed2; 1514 1515 } else { 1516 // get CQ = CX/CY 1517 __div_128_by_128 (&CQ, &CR, CX, CY); 1518 1519 // get number of decimal digits in CQ 1520 // 2^64 1521 f64.i = 0x5f800000; 1522 fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0]; 1523 // binary expon. of CQ 1524 bin_expon = (fx.i - 0x3f800000) >> 23; 1525 1526 digits_q = estimate_decimal_digits[bin_expon]; 1527 TP128.w[0] = power10_index_binexp_128[bin_expon].w[0]; 1528 TP128.w[1] = power10_index_binexp_128[bin_expon].w[1]; 1529 if (__unsigned_compare_ge_128 (CQ, TP128)) 1530 digits_q++; 1531 1532 if (digits_q <= 16) { 1533 if (!CR.w[1] && !CR.w[0]) { 1534 res = get_BID64 (sign_x ^ sign_y, diff_expon, 1535 CQ.w[0], rnd_mode, pfpsf); 1536 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 1537 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 1538 #endif 1539 BID_RETURN (res); 1540 } 1541 1542 ed2 = 16 - digits_q; 1543 T128.w[0] = power10_table_128[ed2].w[0]; 1544 __mul_64x128_to_192 (CA4, (T128.w[0]), CR); 1545 diff_expon = diff_expon - ed2; 1546 CQ.w[0] *= T128.w[0]; 1547 } else { 1548 ed2 = digits_q - 16; 1549 diff_expon += ed2; 1550 T128 = reciprocals10_128[ed2]; 1551 __mul_128x128_to_256 (P256, CQ, T128); 1552 amount = recip_scale[ed2]; 1553 CQ.w[0] = (P256.w[2] >> amount) | (P256.w[3] << (64 - amount)); 1554 CQ.w[1] = 0; 1555 1556 __mul_64x64_to_128 (CQ2, CQ.w[0], (power10_table_128[ed2].w[0])); 1557 1558 __mul_64x64_to_128 (QB256, CQ2.w[0], CY.w[0]); 1559 QB256.w[1] += CQ2.w[0] * CY.w[1] + CQ2.w[1] * CY.w[0]; 1560 1561 CA4.w[1] = CX.w[1] - QB256.w[1]; 1562 CA4.w[0] = CX.w[0] - QB256.w[0]; 1563 if (CX.w[0] < QB256.w[0]) 1564 CA4.w[1]--; 1565 if (CR.w[0] || CR.w[1]) 1566 CA4.w[0] |= 1; 1567 done = 1; 1568 if(CA4.w[1]|CA4.w[0]) { 1569 __mul_64x128_low(CY, (power10_table_128[ed2].w[0]),CY); 1570 } 1571 } 1572 1573 } 1574 1575 if (!done) { 1576 __div_256_by_128 (&CQ, &CA4, CY); 1577 } 1578 1579 1580 1581 #ifdef SET_STATUS_FLAGS 1582 if (CA4.w[0] || CA4.w[1]) { 1583 // set status flags 1584 __set_status_flags (pfpsf, INEXACT_EXCEPTION); 1585 } 1586 #ifndef LEAVE_TRAILING_ZEROS 1587 else 1588 #endif 1589 #else 1590 #ifndef LEAVE_TRAILING_ZEROS 1591 if (!CA4.w[0] && !CA4.w[1]) 1592 #endif 1593 #endif 1594 #ifndef LEAVE_TRAILING_ZEROS 1595 // check whether result is exact 1596 { 1597 if(!done) { 1598 // check whether CX, CY are short 1599 if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) { 1600 i = (int) CY.w[0] - 1; 1601 j = (int) CX.w[0] - 1; 1602 // difference in powers of 2 factors for Y and X 1603 nzeros = ed2 - factors[i][0] + factors[j][0]; 1604 // difference in powers of 5 factors 1605 d5 = ed2 - factors[i][1] + factors[j][1]; 1606 if (d5 < nzeros) 1607 nzeros = d5; 1608 // get P*(2^M[extra_digits])/10^extra_digits 1609 __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]); 1610 //__mul_128x128_to_256(P256, CQ, reciprocals10_128[nzeros]);Qh.w[1]=P256.w[3];Qh.w[0]=P256.w[2]; 1611 1612 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 1613 amount = recip_scale[nzeros]; 1614 __shr_128_long (CQ, Qh, amount); 1615 1616 diff_expon += nzeros; 1617 } else { 1618 // decompose Q as Qh*10^17 + Ql 1619 //T128 = reciprocals10_128[17]; 1620 Q_low = CQ.w[0]; 1621 1622 { 1623 tdigit[0] = Q_low & 0x3ffffff; 1624 tdigit[1] = 0; 1625 QX = Q_low >> 26; 1626 QX32 = QX; 1627 nzeros = 0; 1628 1629 for (j = 0; QX32; j++, QX32 >>= 7) { 1630 k = (QX32 & 127); 1631 tdigit[0] += convert_table[j][k][0]; 1632 tdigit[1] += convert_table[j][k][1]; 1633 if (tdigit[0] >= 100000000) { 1634 tdigit[0] -= 100000000; 1635 tdigit[1]++; 1636 } 1637 } 1638 1639 if (tdigit[1] >= 100000000) { 1640 tdigit[1] -= 100000000; 1641 if (tdigit[1] >= 100000000) 1642 tdigit[1] -= 100000000; 1643 } 1644 1645 digit = tdigit[0]; 1646 if (!digit && !tdigit[1]) 1647 nzeros += 16; 1648 else { 1649 if (!digit) { 1650 nzeros += 8; 1651 digit = tdigit[1]; 1652 } 1653 // decompose digit 1654 PD = (UINT64) digit *0x068DB8BBull; 1655 digit_h = (UINT32) (PD >> 40); 1656 digit_low = digit - digit_h * 10000; 1657 1658 if (!digit_low) 1659 nzeros += 4; 1660 else 1661 digit_h = digit_low; 1662 1663 if (!(digit_h & 1)) 1664 nzeros += 1665 3 & (UINT32) (packed_10000_zeros[digit_h >> 3] >> 1666 (digit_h & 7)); 1667 } 1668 1669 if (nzeros) { 1670 // get P*(2^M[extra_digits])/10^extra_digits 1671 __mul_128x128_high (Qh, CQ, reciprocals10_128[nzeros]); 1672 1673 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 1674 amount = recip_scale[nzeros]; 1675 __shr_128 (CQ, Qh, amount); 1676 } 1677 diff_expon += nzeros; 1678 1679 } 1680 } 1681 } 1682 if(diff_expon>=0){ 1683 res = 1684 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0], 1685 rnd_mode, pfpsf); 1686 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 1687 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 1688 #endif 1689 BID_RETURN (res); 1690 } 1691 } 1692 #endif 1693 1694 if(diff_expon>=0) { 1695 1696 #ifdef IEEE_ROUND_NEAREST 1697 // rounding 1698 // 2*CA4 - CY 1699 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63); 1700 CA4r.w[0] = CA4.w[0] + CA4.w[0]; 1701 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]); 1702 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64; 1703 1704 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0; 1705 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D); 1706 1707 CQ.w[0] += carry64; 1708 //if(CQ.w[0]<carry64) 1709 //CQ.w[1] ++; 1710 #else 1711 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY 1712 // rounding 1713 // 2*CA4 - CY 1714 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63); 1715 CA4r.w[0] = CA4.w[0] + CA4.w[0]; 1716 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]); 1717 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64; 1718 1719 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1; 1720 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D; 1721 1722 CQ.w[0] += carry64; 1723 if (CQ.w[0] < carry64) 1724 CQ.w[1]++; 1725 #else 1726 rmode = rnd_mode; 1727 if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2) 1728 rmode = 3 - rmode; 1729 switch (rmode) { 1730 case ROUNDING_TO_NEAREST: // round to nearest code 1731 // rounding 1732 // 2*CA4 - CY 1733 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63); 1734 CA4r.w[0] = CA4.w[0] + CA4.w[0]; 1735 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]); 1736 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64; 1737 D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0; 1738 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D); 1739 CQ.w[0] += carry64; 1740 if (CQ.w[0] < carry64) 1741 CQ.w[1]++; 1742 break; 1743 case ROUNDING_TIES_AWAY: 1744 // rounding 1745 // 2*CA4 - CY 1746 CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63); 1747 CA4r.w[0] = CA4.w[0] + CA4.w[0]; 1748 __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]); 1749 CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64; 1750 D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1; 1751 carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D; 1752 CQ.w[0] += carry64; 1753 if (CQ.w[0] < carry64) 1754 CQ.w[1]++; 1755 break; 1756 case ROUNDING_DOWN: 1757 case ROUNDING_TO_ZERO: 1758 break; 1759 default: // rounding up 1760 CQ.w[0]++; 1761 if (!CQ.w[0]) 1762 CQ.w[1]++; 1763 break; 1764 } 1765 #endif 1766 #endif 1767 1768 1769 res = 1770 fast_get_BID64_check_OF (sign_x ^ sign_y, diff_expon, CQ.w[0], rnd_mode, 1771 pfpsf); 1772 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 1773 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 1774 #endif 1775 BID_RETURN (res); 1776 } else { 1777 // UF occurs 1778 1779 #ifdef SET_STATUS_FLAGS 1780 if ((diff_expon + 16 < 0)) { 1781 // set status flags 1782 __set_status_flags (pfpsf, INEXACT_EXCEPTION); 1783 } 1784 #endif 1785 rmode = rnd_mode; 1786 res = 1787 get_BID64_UF (sign_x ^ sign_y, diff_expon, CQ.w[0], CA4.w[1] | CA4.w[0], rmode, pfpsf); 1788 #ifdef UNCHANGED_BINARY_STATUS_FLAGS 1789 (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS); 1790 #endif 1791 BID_RETURN (res); 1792 1793 } 1794 1795 } 1796