1/* mpfr_div -- divide two floating-point numbers 2 3Copyright 1999, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23/* References: 24 [1] Short Division of Long Integers, David Harvey and Paul Zimmermann, 25 Proceedings of the 20th Symposium on Computer Arithmetic (ARITH-20), 26 July 25-27, 2011, pages 7-14. 27*/ 28 29#define MPFR_NEED_LONGLONG_H 30#include "mpfr-impl.h" 31 32#ifdef DEBUG2 33#define mpfr_mpn_print(ap,n) mpfr_mpn_print3 (ap,n,MPFR_LIMB_ZERO) 34static void 35mpfr_mpn_print3 (mpfr_limb_ptr ap, mp_size_t n, mp_limb_t cy) 36{ 37 mp_size_t i; 38 for (i = 0; i < n; i++) 39 printf ("+%lu*2^%lu", (unsigned long) ap[i], (unsigned long) 40 (GMP_NUMB_BITS * i)); 41 if (cy) 42 printf ("+2^%lu", (unsigned long) (GMP_NUMB_BITS * n)); 43 printf ("\n"); 44} 45#endif 46 47/* check if {ap, an} is zero */ 48static int 49mpfr_mpn_cmpzero (mpfr_limb_ptr ap, mp_size_t an) 50{ 51 while (an > 0) 52 if (MPFR_LIKELY(ap[--an] != MPFR_LIMB_ZERO)) 53 return 1; 54 return 0; 55} 56 57/* compare {ap, an} and {bp, bn} >> extra, 58 aligned by the more significant limbs. 59 Takes into account bp[0] for extra=1. 60*/ 61static int 62mpfr_mpn_cmp_aux (mpfr_limb_ptr ap, mp_size_t an, 63 mpfr_limb_ptr bp, mp_size_t bn, int extra) 64{ 65 int cmp = 0; 66 mp_size_t k; 67 mp_limb_t bb; 68 69 if (an >= bn) 70 { 71 k = an - bn; 72 while (cmp == 0 && bn > 0) 73 { 74 bn --; 75 bb = (extra) ? ((bp[bn+1] << (GMP_NUMB_BITS - 1)) | (bp[bn] >> 1)) 76 : bp[bn]; 77 cmp = (ap[k + bn] > bb) ? 1 : ((ap[k + bn] < bb) ? -1 : 0); 78 } 79 bb = (extra) ? bp[0] << (GMP_NUMB_BITS - 1) : MPFR_LIMB_ZERO; 80 while (cmp == 0 && k > 0) 81 { 82 k--; 83 cmp = (ap[k] > bb) ? 1 : ((ap[k] < bb) ? -1 : 0); 84 bb = MPFR_LIMB_ZERO; /* ensure we consider only once bp[0] & 1 */ 85 } 86 if (cmp == 0 && bb != MPFR_LIMB_ZERO) 87 cmp = -1; 88 } 89 else /* an < bn */ 90 { 91 k = bn - an; 92 while (cmp == 0 && an > 0) 93 { 94 an --; 95 bb = (extra) ? ((bp[k+an+1] << (GMP_NUMB_BITS - 1)) | (bp[k+an] >> 1)) 96 : bp[k+an]; 97 if (ap[an] > bb) 98 cmp = 1; 99 else if (ap[an] < bb) 100 cmp = -1; 101 } 102 while (cmp == 0 && k > 0) 103 { 104 k--; 105 bb = (extra) ? ((bp[k+1] << (GMP_NUMB_BITS - 1)) | (bp[k] >> 1)) 106 : bp[k]; 107 cmp = (bb != MPFR_LIMB_ZERO) ? -1 : 0; 108 } 109 if (cmp == 0 && extra && (bp[0] & MPFR_LIMB_ONE)) 110 cmp = -1; 111 } 112 return cmp; 113} 114 115/* {ap, n} <- {ap, n} - {bp, n} >> extra - cy, with cy = 0 or 1. 116 Return borrow out. 117*/ 118static mp_limb_t 119mpfr_mpn_sub_aux (mpfr_limb_ptr ap, mpfr_limb_ptr bp, mp_size_t n, 120 mp_limb_t cy, int extra) 121{ 122 mp_limb_t bb, rp; 123 124 MPFR_ASSERTD (cy <= 1); 125 while (n--) 126 { 127 bb = (extra) ? ((bp[1] << (GMP_NUMB_BITS-1)) | (bp[0] >> 1)) : bp[0]; 128 rp = ap[0] - bb - cy; 129 cy = (ap[0] < bb) || (cy && ~rp == MPFR_LIMB_ZERO) ? 130 MPFR_LIMB_ONE : MPFR_LIMB_ZERO; 131 ap[0] = rp; 132 ap ++; 133 bp ++; 134 } 135 MPFR_ASSERTD (cy <= 1); 136 return cy; 137} 138 139int 140mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode) 141{ 142 mp_size_t q0size = MPFR_LIMB_SIZE(q); /* number of limbs of destination */ 143 mp_size_t usize = MPFR_LIMB_SIZE(u); 144 mp_size_t vsize = MPFR_LIMB_SIZE(v); 145 mp_size_t qsize; /* number of limbs wanted for the computed quotient */ 146 mp_size_t qqsize; 147 mp_size_t k; 148 mpfr_limb_ptr q0p = MPFR_MANT(q), qp; 149 mpfr_limb_ptr up = MPFR_MANT(u); 150 mpfr_limb_ptr vp = MPFR_MANT(v); 151 mpfr_limb_ptr ap; 152 mpfr_limb_ptr bp; 153 mp_limb_t qh; 154 mp_limb_t sticky_u = MPFR_LIMB_ZERO; 155 mp_limb_t low_u; 156 mp_limb_t sticky_v = MPFR_LIMB_ZERO; 157 mp_limb_t sticky; 158 mp_limb_t sticky3; 159 mp_limb_t round_bit = MPFR_LIMB_ZERO; 160 mpfr_exp_t qexp; 161 int sign_quotient; 162 int extra_bit; 163 int sh, sh2; 164 int inex; 165 int like_rndz; 166 MPFR_TMP_DECL(marker); 167 168 MPFR_LOG_FUNC ( 169 ("u[%Pu]=%.*Rg v[%Pu]=%.*Rg rnd=%d", 170 mpfr_get_prec(u), mpfr_log_prec, u, 171 mpfr_get_prec (v),mpfr_log_prec, v, rnd_mode), 172 ("q[%Pu]=%.*Rg inexact=%d", mpfr_get_prec(q), mpfr_log_prec, q, inex)); 173 174 /************************************************************************** 175 * * 176 * This part of the code deals with special cases * 177 * * 178 **************************************************************************/ 179 180 if (MPFR_UNLIKELY(MPFR_ARE_SINGULAR(u,v))) 181 { 182 if (MPFR_IS_NAN(u) || MPFR_IS_NAN(v)) 183 { 184 MPFR_SET_NAN(q); 185 MPFR_RET_NAN; 186 } 187 sign_quotient = MPFR_MULT_SIGN( MPFR_SIGN(u) , MPFR_SIGN(v) ); 188 MPFR_SET_SIGN(q, sign_quotient); 189 if (MPFR_IS_INF(u)) 190 { 191 if (MPFR_IS_INF(v)) 192 { 193 MPFR_SET_NAN(q); 194 MPFR_RET_NAN; 195 } 196 else 197 { 198 MPFR_SET_INF(q); 199 MPFR_RET(0); 200 } 201 } 202 else if (MPFR_IS_INF(v)) 203 { 204 MPFR_SET_ZERO (q); 205 MPFR_RET (0); 206 } 207 else if (MPFR_IS_ZERO (v)) 208 { 209 if (MPFR_IS_ZERO (u)) 210 { 211 MPFR_SET_NAN(q); 212 MPFR_RET_NAN; 213 } 214 else 215 { 216 MPFR_ASSERTD (! MPFR_IS_INF (u)); 217 MPFR_SET_INF(q); 218 mpfr_set_divby0 (); 219 MPFR_RET(0); 220 } 221 } 222 else 223 { 224 MPFR_ASSERTD (MPFR_IS_ZERO (u)); 225 MPFR_SET_ZERO (q); 226 MPFR_RET (0); 227 } 228 } 229 230 /************************************************************************** 231 * * 232 * End of the part concerning special values. * 233 * * 234 **************************************************************************/ 235 236 MPFR_TMP_MARK(marker); 237 238 /* set sign */ 239 sign_quotient = MPFR_MULT_SIGN( MPFR_SIGN(u) , MPFR_SIGN(v) ); 240 MPFR_SET_SIGN(q, sign_quotient); 241 242 /* determine if an extra bit comes from the division, i.e. if the 243 significand of u (as a fraction in [1/2, 1[) is larger than that 244 of v */ 245 if (MPFR_LIKELY(up[usize - 1] != vp[vsize - 1])) 246 extra_bit = (up[usize - 1] > vp[vsize - 1]) ? 1 : 0; 247 else /* most significant limbs are equal, must look at further limbs */ 248 { 249 mp_size_t l; 250 251 k = usize - 1; 252 l = vsize - 1; 253 while (k != 0 && l != 0 && up[--k] == vp[--l]); 254 /* now k=0 or l=0 or up[k] != vp[l] */ 255 if (up[k] > vp[l]) 256 extra_bit = 1; 257 else if (up[k] < vp[l]) 258 extra_bit = 0; 259 /* now up[k] = vp[l], thus either k=0 or l=0 */ 260 else if (l == 0) /* no more divisor limb */ 261 extra_bit = 1; 262 else /* k=0: no more dividend limb */ 263 extra_bit = mpfr_mpn_cmpzero (vp, l) == 0; 264 } 265#ifdef DEBUG 266 printf ("extra_bit=%d\n", extra_bit); 267#endif 268 269 /* set exponent */ 270 qexp = MPFR_GET_EXP (u) - MPFR_GET_EXP (v) + extra_bit; 271 272 /* sh is the number of zero bits in the low limb of the quotient */ 273 MPFR_UNSIGNED_MINUS_MODULO(sh, MPFR_PREC(q)); 274 275 like_rndz = rnd_mode == MPFR_RNDZ || 276 rnd_mode == (sign_quotient < 0 ? MPFR_RNDU : MPFR_RNDD); 277 278 /************************************************************************** 279 * * 280 * We first try Mulders' short division (for large operands) * 281 * * 282 **************************************************************************/ 283 284 if (MPFR_UNLIKELY(q0size >= MPFR_DIV_THRESHOLD && 285 vsize >= MPFR_DIV_THRESHOLD)) 286 { 287 mp_size_t n = q0size + 1; /* we will perform a short (2n)/n division */ 288 mpfr_limb_ptr ap, bp, qp; 289 mpfr_prec_t p; 290 291 /* since Mulders' short division clobbers the dividend, we have to 292 copy it */ 293 ap = MPFR_TMP_LIMBS_ALLOC (n + n); 294 if (usize >= n + n) /* truncate the dividend */ 295 MPN_COPY(ap, up + usize - (n + n), n + n); 296 else /* zero-pad the dividend */ 297 { 298 MPN_COPY(ap + (n + n) - usize, up, usize); 299 MPN_ZERO(ap, (n + n) - usize); 300 } 301 302 if (vsize >= n) /* truncate the divisor */ 303 bp = vp + vsize - n; 304 else /* zero-pad the divisor */ 305 { 306 bp = MPFR_TMP_LIMBS_ALLOC (n); 307 MPN_COPY(bp + n - vsize, vp, vsize); 308 MPN_ZERO(bp, n - vsize); 309 } 310 311 qp = MPFR_TMP_LIMBS_ALLOC (n); 312 qh = mpfr_divhigh_n (qp, ap, bp, n); 313 /* in all cases, the error is at most (2n+2) ulps on qh*B^n+{qp,n}, 314 cf algorithms.tex */ 315 316 p = n * GMP_NUMB_BITS - MPFR_INT_CEIL_LOG2 (2 * n + 2); 317 /* if qh is 1, then we need only PREC(q)-1 bits of {qp,n}, 318 if rnd=RNDN, we need to be able to round with a directed rounding 319 and one more bit */ 320 if (MPFR_LIKELY (mpfr_round_p (qp, n, p, 321 MPFR_PREC(q) + (rnd_mode == MPFR_RNDN) - qh))) 322 { 323 /* we can round correctly whatever the rounding mode */ 324 if (qh == 0) 325 MPN_COPY (q0p, qp + 1, q0size); 326 else 327 { 328 mpn_rshift (q0p, qp + 1, q0size, 1); 329 q0p[q0size - 1] ^= MPFR_LIMB_HIGHBIT; 330 } 331 q0p[0] &= ~MPFR_LIMB_MASK(sh); /* put to zero low sh bits */ 332 333 if (rnd_mode == MPFR_RNDN) /* round to nearest */ 334 { 335 /* we know we can round, thus we are never in the even rule case: 336 if the round bit is 0, we truncate 337 if the round bit is 1, we add 1 */ 338 if (qh == 0) 339 { 340 if (sh > 0) 341 round_bit = (qp[1] >> (sh - 1)) & 1; 342 else 343 round_bit = qp[0] >> (GMP_NUMB_BITS - 1); 344 } 345 else /* qh = 1 */ 346 round_bit = (qp[1] >> sh) & 1; 347 if (round_bit == 0) 348 { 349 inex = -1; 350 goto truncate; 351 } 352 else /* round_bit = 1 */ 353 goto add_one_ulp; 354 } 355 else if (like_rndz == 0) /* round away */ 356 goto add_one_ulp; 357 /* else round to zero: nothing to do */ 358 else 359 { 360 inex = -1; 361 goto truncate; 362 } 363 } 364 } 365 366 /************************************************************************** 367 * * 368 * Mulders' short division failed: we revert to integer division * 369 * * 370 **************************************************************************/ 371 372 if (MPFR_UNLIKELY(rnd_mode == MPFR_RNDN && sh == 0)) 373 { /* we compute the quotient with one more limb, in order to get 374 the round bit in the quotient, and the remainder only contains 375 sticky bits */ 376 qsize = q0size + 1; 377 /* need to allocate memory for the quotient */ 378 qp = MPFR_TMP_LIMBS_ALLOC (qsize); 379 } 380 else 381 { 382 qsize = q0size; 383 qp = q0p; /* directly put the quotient in the destination */ 384 } 385 qqsize = qsize + qsize; 386 387 /* prepare the dividend */ 388 ap = MPFR_TMP_LIMBS_ALLOC (qqsize); 389 if (MPFR_LIKELY(qqsize > usize)) /* use the full dividend */ 390 { 391 k = qqsize - usize; /* k > 0 */ 392 MPN_ZERO(ap, k); 393 if (extra_bit) 394 ap[k - 1] = mpn_rshift (ap + k, up, usize, 1); 395 else 396 MPN_COPY(ap + k, up, usize); 397 } 398 else /* truncate the dividend */ 399 { 400 k = usize - qqsize; 401 if (extra_bit) 402 sticky_u = mpn_rshift (ap, up + k, qqsize, 1); 403 else 404 MPN_COPY(ap, up + k, qqsize); 405 sticky_u = sticky_u || mpfr_mpn_cmpzero (up, k); 406 } 407 low_u = sticky_u; 408 409 /* now sticky_u is non-zero iff the truncated part of u is non-zero */ 410 411 /* prepare the divisor */ 412 if (MPFR_LIKELY(vsize >= qsize)) 413 { 414 k = vsize - qsize; 415 if (qp != vp) 416 bp = vp + k; /* avoid copying the divisor */ 417 else /* need to copy, since mpn_divrem doesn't allow overlap 418 between quotient and divisor, necessarily k = 0 419 since quotient and divisor are the same mpfr variable */ 420 { 421 bp = MPFR_TMP_LIMBS_ALLOC (qsize); 422 MPN_COPY(bp, vp, vsize); 423 } 424 sticky_v = sticky_v || mpfr_mpn_cmpzero (vp, k); 425 k = 0; 426 } 427 else /* vsize < qsize: small divisor case */ 428 { 429 bp = vp; 430 k = qsize - vsize; 431 } 432 433 /************************************************************************** 434 * * 435 * Here we perform the real division of {ap+k,qqsize-k} by {bp,qsize-k} * 436 * * 437 **************************************************************************/ 438 439 /* if Mulders' short division failed, we revert to division with remainder */ 440 qh = mpn_divrem (qp, 0, ap + k, qqsize - k, bp, qsize - k); 441 /* warning: qh may be 1 if u1 == v1, but u < v */ 442#ifdef DEBUG2 443 printf ("q="); mpfr_mpn_print (qp, qsize); 444 printf ("r="); mpfr_mpn_print (ap, qsize); 445#endif 446 447 k = qsize; 448 sticky_u = sticky_u || mpfr_mpn_cmpzero (ap, k); 449 450 sticky = sticky_u | sticky_v; 451 452 /* now sticky is non-zero iff one of the following holds: 453 (a) the truncated part of u is non-zero 454 (b) the truncated part of v is non-zero 455 (c) the remainder from division is non-zero */ 456 457 if (MPFR_LIKELY(qsize == q0size)) 458 { 459 sticky3 = qp[0] & MPFR_LIMB_MASK(sh); /* does nothing when sh=0 */ 460 sh2 = sh; 461 } 462 else /* qsize = q0size + 1: only happens when rnd_mode=MPFR_RNDN and sh=0 */ 463 { 464 MPN_COPY (q0p, qp + 1, q0size); 465 sticky3 = qp[0]; 466 sh2 = GMP_NUMB_BITS; 467 } 468 qp[0] ^= sticky3; 469 /* sticky3 contains the truncated bits from the quotient, 470 including the round bit, and 1 <= sh2 <= GMP_NUMB_BITS 471 is the number of bits in sticky3 */ 472 inex = (sticky != MPFR_LIMB_ZERO) || (sticky3 != MPFR_LIMB_ZERO); 473#ifdef DEBUG 474 printf ("sticky=%lu sticky3=%lu inex=%d\n", 475 (unsigned long) sticky, (unsigned long) sticky3, inex); 476#endif 477 478 /* to round, we distinguish two cases: 479 (a) vsize <= qsize: we used the full divisor 480 (b) vsize > qsize: the divisor was truncated 481 */ 482 483#ifdef DEBUG 484 printf ("vsize=%lu qsize=%lu\n", 485 (unsigned long) vsize, (unsigned long) qsize); 486#endif 487 if (MPFR_LIKELY(vsize <= qsize)) /* use the full divisor */ 488 { 489 if (MPFR_LIKELY(rnd_mode == MPFR_RNDN)) 490 { 491 round_bit = sticky3 & (MPFR_LIMB_ONE << (sh2 - 1)); 492 sticky = (sticky3 ^ round_bit) | sticky_u; 493 } 494 else if (like_rndz || inex == 0) 495 sticky = (inex == 0) ? MPFR_LIMB_ZERO : MPFR_LIMB_ONE; 496 else /* round away from zero */ 497 sticky = MPFR_LIMB_ONE; 498 goto case_1; 499 } 500 else /* vsize > qsize: need to truncate the divisor */ 501 { 502 if (inex == 0) 503 goto truncate; 504 else 505 { 506 /* We know the estimated quotient is an upper bound of the exact 507 quotient (with rounding toward zero), with a difference of at 508 most 2 in qp[0]. 509 Thus we can round except when sticky3 is 000...000 or 000...001 510 for directed rounding, and 100...000 or 100...001 for rounding 511 to nearest. (For rounding to nearest, we cannot determine the 512 inexact flag for 000...000 or 000...001.) 513 */ 514 mp_limb_t sticky3orig = sticky3; 515 if (rnd_mode == MPFR_RNDN) 516 { 517 round_bit = sticky3 & (MPFR_LIMB_ONE << (sh2 - 1)); 518 sticky3 = sticky3 ^ round_bit; 519#ifdef DEBUG 520 printf ("rb=%lu sb=%lu\n", 521 (unsigned long) round_bit, (unsigned long) sticky3); 522#endif 523 } 524 if (sticky3 != MPFR_LIMB_ZERO && sticky3 != MPFR_LIMB_ONE) 525 { 526 sticky = sticky3; 527 goto case_1; 528 } 529 else /* hard case: we have to compare q1 * v0 and r + low(u), 530 where q1 * v0 has qsize + (vsize-qsize) = vsize limbs, and 531 r + low(u) has qsize + (usize-2*qsize) = usize-qsize limbs */ 532 { 533 mp_size_t l; 534 mpfr_limb_ptr sp; 535 int cmp_s_r; 536 mp_limb_t qh2; 537 538 sp = MPFR_TMP_LIMBS_ALLOC (vsize); 539 k = vsize - qsize; 540 /* sp <- {qp, qsize} * {vp, vsize-qsize} */ 541 qp[0] ^= sticky3orig; /* restore original quotient */ 542 if (qsize >= k) 543 mpn_mul (sp, qp, qsize, vp, k); 544 else 545 mpn_mul (sp, vp, k, qp, qsize); 546 if (qh) 547 qh2 = mpn_add_n (sp + qsize, sp + qsize, vp, k); 548 else 549 qh2 = (mp_limb_t) 0; 550 qp[0] ^= sticky3orig; /* restore truncated quotient */ 551 552 /* compare qh2 + {sp, k + qsize} to {ap, qsize} + low(u) */ 553 cmp_s_r = (qh2 != 0) ? 1 : mpn_cmp (sp + k, ap, qsize); 554 if (cmp_s_r == 0) /* compare {sp, k} and low(u) */ 555 { 556 cmp_s_r = (usize >= qqsize) ? 557 mpfr_mpn_cmp_aux (sp, k, up, usize - qqsize, extra_bit) : 558 mpfr_mpn_cmpzero (sp, k); 559 } 560#ifdef DEBUG 561 printf ("cmp(q*v0,r+u0)=%d\n", cmp_s_r); 562#endif 563 /* now cmp_s_r > 0 if {sp, vsize} > {ap, qsize} + low(u) 564 cmp_s_r = 0 if {sp, vsize} = {ap, qsize} + low(u) 565 cmp_s_r < 0 if {sp, vsize} < {ap, qsize} + low(u) */ 566 if (cmp_s_r <= 0) /* quotient is in [q1, q1+1) */ 567 { 568 sticky = (cmp_s_r == 0) ? sticky3 : MPFR_LIMB_ONE; 569 goto case_1; 570 } 571 else /* cmp_s_r > 0, quotient is < q1: to determine if it is 572 in [q1-2,q1-1] or in [q1-1,q1], we need to subtract 573 the low part u0 of the dividend u0 from q*v0 */ 574 { 575 mp_limb_t cy = MPFR_LIMB_ZERO; 576 577 /* subtract low(u)>>extra_bit if non-zero */ 578 if (qh2 != 0) /* whatever the value of {up, m + k}, it 579 will be smaller than qh2 + {sp, k} */ 580 cmp_s_r = 1; 581 else 582 { 583 if (low_u != MPFR_LIMB_ZERO) 584 { 585 mp_size_t m; 586 l = usize - qqsize; /* number of low limbs in u */ 587 m = (l > k) ? l - k : 0; 588 cy = (extra_bit) ? 589 (up[m] & MPFR_LIMB_ONE) : MPFR_LIMB_ZERO; 590 if (l >= k) /* u0 has more limbs than s: 591 first look if {up, m} is not zero, 592 and compare {sp, k} and {up + m, k} */ 593 { 594 cy = cy || mpfr_mpn_cmpzero (up, m); 595 low_u = cy; 596 cy = mpfr_mpn_sub_aux (sp, up + m, k, 597 cy, extra_bit); 598 } 599 else /* l < k: s has more limbs than u0 */ 600 { 601 low_u = MPFR_LIMB_ZERO; 602 if (cy != MPFR_LIMB_ZERO) 603 cy = mpn_sub_1 (sp + k - l - 1, sp + k - l - 1, 604 1, MPFR_LIMB_HIGHBIT); 605 cy = mpfr_mpn_sub_aux (sp + k - l, up, l, 606 cy, extra_bit); 607 } 608 } 609 MPFR_ASSERTD (cy <= 1); 610 cy = mpn_sub_1 (sp + k, sp + k, qsize, cy); 611 /* subtract r */ 612 cy += mpn_sub_n (sp + k, sp + k, ap, qsize); 613 MPFR_ASSERTD (cy <= 1); 614 /* now compare {sp, ssize} to v */ 615 cmp_s_r = mpn_cmp (sp, vp, vsize); 616 if (cmp_s_r == 0 && low_u != MPFR_LIMB_ZERO) 617 cmp_s_r = 1; /* since in fact we subtracted 618 less than 1 */ 619 } 620#ifdef DEBUG 621 printf ("cmp(q*v0-(r+u0),v)=%d\n", cmp_s_r); 622#endif 623 if (cmp_s_r <= 0) /* q1-1 <= u/v < q1 */ 624 { 625 if (sticky3 == MPFR_LIMB_ONE) 626 { /* q1-1 is either representable (directed rounding), 627 or the middle of two numbers (nearest) */ 628 sticky = (cmp_s_r) ? MPFR_LIMB_ONE : MPFR_LIMB_ZERO; 629 goto case_1; 630 } 631 /* now necessarily sticky3=0 */ 632 else if (round_bit == MPFR_LIMB_ZERO) 633 { /* round_bit=0, sticky3=0: q1-1 is exact only 634 when sh=0 */ 635 inex = (cmp_s_r || sh) ? -1 : 0; 636 if (rnd_mode == MPFR_RNDN || 637 (! like_rndz && inex != 0)) 638 { 639 inex = 1; 640 goto truncate_check_qh; 641 } 642 else /* round down */ 643 goto sub_one_ulp; 644 } 645 else /* sticky3=0, round_bit=1 ==> rounding to nearest */ 646 { 647 inex = cmp_s_r; 648 goto truncate; 649 } 650 } 651 else /* q1-2 < u/v < q1-1 */ 652 { 653 /* if rnd=MPFR_RNDN, the result is q1 when 654 q1-2 >= q1-2^(sh-1), i.e. sh >= 2, 655 otherwise (sh=1) it is q1-2 */ 656 if (rnd_mode == MPFR_RNDN) /* sh > 0 */ 657 { 658 /* Case sh=1: sb=0 always, and q1-rb is exactly 659 representable, like q1-rb-2. 660 rb action 661 0 subtract two ulps, inex=-1 662 1 truncate, inex=1 663 664 Case sh>1: one ulp is 2^(sh-1) >= 2 665 rb sb action 666 0 0 truncate, inex=1 667 0 1 truncate, inex=1 668 1 x truncate, inex=-1 669 */ 670 if (sh == 1) 671 { 672 if (round_bit == MPFR_LIMB_ZERO) 673 { 674 inex = -1; 675 sh = 0; 676 goto sub_two_ulp; 677 } 678 else 679 { 680 inex = 1; 681 goto truncate_check_qh; 682 } 683 } 684 else /* sh > 1 */ 685 { 686 inex = (round_bit == MPFR_LIMB_ZERO) ? 1 : -1; 687 goto truncate_check_qh; 688 } 689 } 690 else if (like_rndz) 691 { 692 /* the result is down(q1-2), i.e. subtract one 693 ulp if sh > 0, and two ulps if sh=0 */ 694 inex = -1; 695 if (sh > 0) 696 goto sub_one_ulp; 697 else 698 goto sub_two_ulp; 699 } 700 /* if round away from zero, the result is up(q1-1), 701 which is q1 unless sh = 0, where it is q1-1 */ 702 else 703 { 704 inex = 1; 705 if (sh > 0) 706 goto truncate_check_qh; 707 else /* sh = 0 */ 708 goto sub_one_ulp; 709 } 710 } 711 } 712 } 713 } 714 } 715 716 case_1: /* quotient is in [q1, q1+1), 717 round_bit is the round_bit (0 for directed rounding), 718 sticky the sticky bit */ 719 if (like_rndz || (round_bit == MPFR_LIMB_ZERO && sticky == MPFR_LIMB_ZERO)) 720 { 721 inex = round_bit == MPFR_LIMB_ZERO && sticky == MPFR_LIMB_ZERO ? 0 : -1; 722 goto truncate; 723 } 724 else if (rnd_mode == MPFR_RNDN) /* sticky <> 0 or round <> 0 */ 725 { 726 if (round_bit == MPFR_LIMB_ZERO) /* necessarily sticky <> 0 */ 727 { 728 inex = -1; 729 goto truncate; 730 } 731 /* round_bit = 1 */ 732 else if (sticky != MPFR_LIMB_ZERO) 733 goto add_one_ulp; /* inex=1 */ 734 else /* round_bit=1, sticky=0 */ 735 goto even_rule; 736 } 737 else /* round away from zero, sticky <> 0 */ 738 goto add_one_ulp; /* with inex=1 */ 739 740 sub_two_ulp: 741 /* we cannot subtract MPFR_LIMB_MPFR_LIMB_ONE << (sh+1) since this is 742 undefined for sh = GMP_NUMB_BITS */ 743 qh -= mpn_sub_1 (q0p, q0p, q0size, MPFR_LIMB_ONE << sh); 744 /* go through */ 745 746 sub_one_ulp: 747 qh -= mpn_sub_1 (q0p, q0p, q0size, MPFR_LIMB_ONE << sh); 748 /* go through truncate_check_qh */ 749 750 truncate_check_qh: 751 if (qh) 752 { 753 qexp ++; 754 q0p[q0size - 1] = MPFR_LIMB_HIGHBIT; 755 } 756 goto truncate; 757 758 even_rule: /* has to set inex */ 759 inex = (q0p[0] & (MPFR_LIMB_ONE << sh)) ? 1 : -1; 760 if (inex < 0) 761 goto truncate; 762 /* else go through add_one_ulp */ 763 764 add_one_ulp: 765 inex = 1; /* always here */ 766 if (mpn_add_1 (q0p, q0p, q0size, MPFR_LIMB_ONE << sh)) 767 { 768 qexp ++; 769 q0p[q0size - 1] = MPFR_LIMB_HIGHBIT; 770 } 771 772 truncate: /* inex already set */ 773 774 MPFR_TMP_FREE(marker); 775 776 /* check for underflow/overflow */ 777 if (MPFR_UNLIKELY(qexp > __gmpfr_emax)) 778 return mpfr_overflow (q, rnd_mode, sign_quotient); 779 else if (MPFR_UNLIKELY(qexp < __gmpfr_emin)) 780 { 781 if (rnd_mode == MPFR_RNDN && ((qexp < __gmpfr_emin - 1) || 782 (inex >= 0 && mpfr_powerof2_raw (q)))) 783 rnd_mode = MPFR_RNDZ; 784 return mpfr_underflow (q, rnd_mode, sign_quotient); 785 } 786 MPFR_SET_EXP(q, qexp); 787 788 inex *= sign_quotient; 789 MPFR_RET (inex); 790} 791