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