contrib/mpfr/div.c

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