xref: /dflybsd-src/contrib/mpfr/src/sub1.c (revision 2786097444a0124b5d33763854de247e230c6629)

/* mpfr_sub1 -- internal function to perform a "real" subtraction

Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
Contributed by the AriC and Caramel projects, INRIA.

This file is part of the GNU MPFR Library.

The GNU MPFR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

The GNU MPFR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */

#include "mpfr-impl.h"

/* compute sign(b) * (|b| - |c|), with |b| > |c|, diff_exp = EXP(b) - EXP(c)
   Returns 0 iff result is exact,
   a negative value when the result is less than the exact value,
   a positive value otherwise.
*/

int
mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
  int sign;
  mpfr_uexp_t diff_exp;
  mpfr_prec_t cancel, cancel1;
  mp_size_t cancel2, an, bn, cn, cn0;
  mp_limb_t *ap, *bp, *cp;
  mp_limb_t carry, bb, cc;
  int inexact, shift_b, shift_c, add_exp = 0;
  int cmp_low = 0; /* used for rounding to nearest: 0 if low(b) = low(c),
                      negative if low(b) < low(c), positive if low(b)>low(c) */
  int sh, k;
  MPFR_TMP_DECL(marker);

  MPFR_TMP_MARK(marker);
  ap = MPFR_MANT(a);
  an = MPFR_LIMB_SIZE(a);

  sign = mpfr_cmp2 (b, c, &cancel);
  if (MPFR_UNLIKELY(sign == 0))
    {
      if (rnd_mode == MPFR_RNDD)
        MPFR_SET_NEG (a);
      else
        MPFR_SET_POS (a);
      MPFR_SET_ZERO (a);
      MPFR_RET (0);
    }

  /*
   * If subtraction: sign(a) = sign * sign(b)
   * If addition: sign(a) = sign of the larger argument in absolute value.
   *
   * Both cases can be simplidied in:
   * if (sign>0)
   *    if addition: sign(a) = sign * sign(b) = sign(b)
   *    if subtraction, b is greater, so sign(a) = sign(b)
   * else
   *    if subtraction, sign(a) = - sign(b)
   *    if addition, sign(a) = sign(c) (since c is greater)
   *      But if it is an addition, sign(b) and sign(c) are opposed!
   *      So sign(a) = - sign(b)
   */

  if (sign < 0) /* swap b and c so that |b| > |c| */
    {
      mpfr_srcptr t;
      MPFR_SET_OPPOSITE_SIGN (a,b);
      t = b; b = c; c = t;
    }
  else
    MPFR_SET_SAME_SIGN (a,b);

  /* Check if c is too small.
     A more precise test is to replace 2 by
      (rnd == MPFR_RNDN) + mpfr_power2_raw (b)
      but it is more expensive and not very useful */
  if (MPFR_UNLIKELY (MPFR_GET_EXP (c) <= MPFR_GET_EXP (b)
                     - (mpfr_exp_t) MAX (MPFR_PREC (a), MPFR_PREC (b)) - 2))
    {
      /* Remember, we can't have an exact result! */
      /*   A.AAAAAAAAAAAAAAAAA
         = B.BBBBBBBBBBBBBBB
          -                     C.CCCCCCCCCCCCC */
      /* A = S*ABS(B) +/- ulp(a) */
      MPFR_SET_EXP (a, MPFR_GET_EXP (b));
      MPFR_RNDRAW_EVEN (inexact, a, MPFR_MANT (b), MPFR_PREC (b),
                        rnd_mode, MPFR_SIGN (a),
                        if (MPFR_UNLIKELY ( ++MPFR_EXP (a) > __gmpfr_emax))
                        inexact = mpfr_overflow (a, rnd_mode, MPFR_SIGN (a)));
      /* inexact = mpfr_set4 (a, b, rnd_mode, MPFR_SIGN (a));  */
      if (inexact == 0)
        {
          /* a = b (Exact)
             But we know it isn't (Since we have to remove `c')
             So if we round to Zero, we have to remove one ulp.
             Otherwise the result is correctly rounded. */
          if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a)))
            {
              mpfr_nexttozero (a);
              MPFR_RET (- MPFR_INT_SIGN (a));
            }
          MPFR_RET (MPFR_INT_SIGN (a));
        }
      else
        {
          /*   A.AAAAAAAAAAAAAA
             = B.BBBBBBBBBBBBBBB
              -                   C.CCCCCCCCCCCCC */
          /* It isn't exact so Prec(b) > Prec(a) and the last
             Prec(b)-Prec(a) bits of `b' are not zeros.
             Which means that removing c from b can't generate a carry
             execpt in case of even rounding.
             In all other case the result and the inexact flag should be
             correct (We can't have an exact result).
             In case of EVEN rounding:
               1.BBBBBBBBBBBBBx10
             -                     1.CCCCCCCCCCCC
             = 1.BBBBBBBBBBBBBx01  Rounded to Prec(b)
             = 1.BBBBBBBBBBBBBx    Nearest / Rounded to Prec(a)
             Set gives:
               1.BBBBBBBBBBBBB0   if inexact == EVEN_INEX  (x == 0)
               1.BBBBBBBBBBBBB1+1 if inexact == -EVEN_INEX (x == 1)
             which means we get a wrong rounded result if x==1,
             i.e. inexact= MPFR_EVEN_INEX */
          if (MPFR_UNLIKELY (inexact == MPFR_EVEN_INEX*MPFR_INT_SIGN (a)))
            {
              mpfr_nexttozero (a);
              inexact = -MPFR_INT_SIGN (a);
            }
          MPFR_RET (inexact);
        }
    }

  diff_exp = (mpfr_uexp_t) MPFR_GET_EXP (b) - MPFR_GET_EXP (c);

  /* reserve a space to store b aligned with the result, i.e. shifted by
     (-cancel) % GMP_NUMB_BITS to the right */
  bn      = MPFR_LIMB_SIZE (b);
  MPFR_UNSIGNED_MINUS_MODULO (shift_b, cancel);
  cancel1 = (cancel + shift_b) / GMP_NUMB_BITS;

  /* the high cancel1 limbs from b should not be taken into account */
  if (MPFR_UNLIKELY (shift_b == 0))
    {
      bp = MPFR_MANT(b); /* no need of an extra space */
      /* Ensure ap != bp */
      if (MPFR_UNLIKELY (ap == bp))
        {
          bp = MPFR_TMP_LIMBS_ALLOC (bn);
          MPN_COPY (bp, ap, bn);
        }
    }
  else
    {
      bp = MPFR_TMP_LIMBS_ALLOC (bn + 1);
      bp[0] = mpn_rshift (bp + 1, MPFR_MANT(b), bn++, shift_b);
    }

  /* reserve a space to store c aligned with the result, i.e. shifted by
      (diff_exp-cancel) % GMP_NUMB_BITS to the right */
  cn      = MPFR_LIMB_SIZE(c);
  if ((UINT_MAX % GMP_NUMB_BITS) == (GMP_NUMB_BITS-1)
      && ((-(unsigned) 1)%GMP_NUMB_BITS > 0))
    shift_c = ((mpfr_uexp_t) diff_exp - cancel) % GMP_NUMB_BITS;
  else
    {
      shift_c = diff_exp - (cancel % GMP_NUMB_BITS);
      shift_c = (shift_c + GMP_NUMB_BITS) % GMP_NUMB_BITS;
    }
  MPFR_ASSERTD( shift_c >= 0 && shift_c < GMP_NUMB_BITS);

  if (MPFR_UNLIKELY(shift_c == 0))
    {
       cp = MPFR_MANT(c);
      /* Ensure ap != cp */
      if (ap == cp)
        {
          cp = MPFR_TMP_LIMBS_ALLOC (cn);
          MPN_COPY(cp, ap, cn);
        }
    }
 else
    {
      cp = MPFR_TMP_LIMBS_ALLOC (cn + 1);
      cp[0] = mpn_rshift (cp + 1, MPFR_MANT(c), cn++, shift_c);
    }

#ifdef DEBUG
  printf ("rnd=%s shift_b=%d shift_c=%d diffexp=%lu\n",
          mpfr_print_rnd_mode (rnd_mode), shift_b, shift_c,
          (unsigned long) diff_exp);
#endif

  MPFR_ASSERTD (ap != cp);
  MPFR_ASSERTD (bp != cp);

  /* here we have shift_c = (diff_exp - cancel) % GMP_NUMB_BITS,
        0 <= shift_c < GMP_NUMB_BITS
     thus we want cancel2 = ceil((cancel - diff_exp) / GMP_NUMB_BITS) */

  /* Possible optimization with a C99 compiler (i.e. well-defined
     integer division): if MPFR_PREC_MAX is reduced to
     ((mpfr_prec_t)((mpfr_uprec_t)(~(mpfr_uprec_t)0)>>1) - GMP_NUMB_BITS + 1)
     and diff_exp is of type mpfr_exp_t (no need for mpfr_uexp_t, since
     the sum or difference of 2 exponents must be representable, as used
     by the multiplication code), then the computation of cancel2 could
     be simplified to
       cancel2 = (cancel - (diff_exp - shift_c)) / GMP_NUMB_BITS;
     because cancel, diff_exp and shift_c are all non-negative and
     these variables are signed. */

  MPFR_ASSERTD (cancel >= 0);
  if (cancel >= diff_exp)
    /* Note that cancel is signed and will be converted to mpfr_uexp_t
       (type of diff_exp) in the expression below, so that this will
       work even if cancel is very large and diff_exp = 0. */
    cancel2 = (cancel - diff_exp + (GMP_NUMB_BITS - 1)) / GMP_NUMB_BITS;
  else
    cancel2 = - (mp_size_t) ((diff_exp - cancel) / GMP_NUMB_BITS);
  /* the high cancel2 limbs from b should not be taken into account */
#ifdef DEBUG
  printf ("cancel=%lu cancel1=%lu cancel2=%ld\n",
          (unsigned long) cancel, (unsigned long) cancel1, (long) cancel2);
#endif

  /*               ap[an-1]        ap[0]
             <----------------+-----------|---->
             <----------PREC(a)----------><-sh->
 cancel1
 limbs        bp[bn-cancel1-1]
 <--...-----><----------------+-----------+----------->
  cancel2
  limbs       cp[cn-cancel2-1]                                    cancel2 >= 0
    <--...--><----------------+----------------+---------------->
                (-cancel2)                                        cancel2 < 0
                   limbs      <----------------+---------------->
  */

  /* first part: put in ap[0..an-1] the value of high(b) - high(c),
     where high(b) consists of the high an+cancel1 limbs of b,
     and high(c) consists of the high an+cancel2 limbs of c.
   */

  /* copy high(b) into a */
  if (MPFR_LIKELY(an + (mp_size_t) cancel1 <= bn))
    /* a: <----------------+-----------|---->
       b: <-----------------------------------------> */
      MPN_COPY (ap, bp + bn - (an + cancel1), an);
  else
    /* a: <----------------+-----------|---->
       b: <-------------------------> */
    if ((mp_size_t) cancel1 < bn) /* otherwise b does not overlap with a */
      {
        MPN_ZERO (ap, an + cancel1 - bn);
        MPN_COPY (ap + (an + cancel1 - bn), bp, bn - cancel1);
      }
    else
      MPN_ZERO (ap, an);

#ifdef DEBUG
  printf("after copying high(b), a="); mpfr_print_binary(a); putchar('\n');
#endif

  /* subtract high(c) */
  if (MPFR_LIKELY(an + cancel2 > 0)) /* otherwise c does not overlap with a */
    {
      mp_limb_t *ap2;

      if (cancel2 >= 0)
        {
          if (an + cancel2 <= cn)
            /* a: <----------------------------->
               c: <-----------------------------------------> */
            mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an);
          else
            /* a: <---------------------------->
               c: <-------------------------> */
            {
              ap2 = ap + an + (cancel2 - cn);
              if (cn > cancel2)
                mpn_sub_n (ap2, ap2, cp, cn - cancel2);
            }
        }
      else /* cancel2 < 0 */
        {
          mp_limb_t borrow;

          if (an + cancel2 <= cn)
            /* a: <----------------------------->
               c: <-----------------------------> */
            borrow = mpn_sub_n (ap, ap, cp + cn - (an + cancel2),
                                an + cancel2);
          else
            /* a: <---------------------------->
               c: <----------------> */
            {
              ap2 = ap + an + cancel2 - cn;
              borrow = mpn_sub_n (ap2, ap2, cp, cn);
            }
          ap2 = ap + an + cancel2;
          mpn_sub_1 (ap2, ap2, -cancel2, borrow);
        }
    }

#ifdef DEBUG
  printf("after subtracting high(c), a=");
  mpfr_print_binary(a);
  putchar('\n');
#endif

  /* now perform rounding */
  sh = (mpfr_prec_t) an * GMP_NUMB_BITS - MPFR_PREC(a);
  /* last unused bits from a */
  carry = ap[0] & MPFR_LIMB_MASK (sh);
  ap[0] -= carry;

  if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
    {
      if (MPFR_LIKELY(sh))
        {
          /* can decide except when carry = 2^(sh-1) [middle]
             or carry = 0 [truncate, but cannot decide inexact flag] */
          if (carry > (MPFR_LIMB_ONE << (sh - 1)))
            goto add_one_ulp;
          else if ((0 < carry) && (carry < (MPFR_LIMB_ONE << (sh - 1))))
            {
              inexact = -1; /* result if smaller than exact value */
              goto truncate;
            }
          /* now carry = 2^(sh-1), in which case cmp_low=2,
             or carry = 0, in which case cmp_low=0 */
          cmp_low = (carry == 0) ? 0 : 2;
        }
    }
  else /* directed rounding: set rnd_mode to RNDZ iff toward zero */
    {
      if (MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd_mode, MPFR_IS_NEG(a)))
        rnd_mode = MPFR_RNDZ;

      if (carry)
        {
          if (rnd_mode == MPFR_RNDZ)
            {
              inexact = -1;
              goto truncate;
            }
          else /* round away */
            goto add_one_ulp;
        }
    }

  /* we have to consider the low (bn - (an+cancel1)) limbs from b,
     and the (cn - (an+cancel2)) limbs from c. */
  bn -= an + cancel1;
  cn0 = cn;
  cn -= an + cancel2;

#ifdef DEBUG
  printf ("last sh=%d bits from a are %lu, bn=%ld, cn=%ld\n",
          sh, (unsigned long) carry, (long) bn, (long) cn);
#endif

  /* for rounding to nearest, we couldn't conclude up to here in the following
     cases:
     1. sh = 0, then cmp_low=0: we can either truncate, subtract one ulp
        or add one ulp: -1 ulp < low(b)-low(c) < 1 ulp
     2. sh > 0 but the low sh bits from high(b)-high(c) equal 2^(sh-1):
        -0.5 ulp <= -1/2^sh < low(b)-low(c)-0.5 < 1/2^sh <= 0.5 ulp
        we can't decide the rounding, in that case cmp_low=2:
        either we truncate and flag=-1, or we add one ulp and flag=1
     3. the low sh>0 bits from high(b)-high(c) equal 0: we know we have to
        truncate but we can't decide the ternary value, here cmp_low=0:
        -0.5 ulp <= -1/2^sh < low(b)-low(c) < 1/2^sh <= 0.5 ulp
        we always truncate and inexact can be any of -1,0,1
  */

  /* note: here cn might exceed cn0, in which case we consider a zero limb */
  for (k = 0; (bn > 0) || (cn > 0); k = 1)
    {
      /* if cmp_low < 0, we know low(b) - low(c) < 0
         if cmp_low > 0, we know low(b) - low(c) > 0
            (more precisely if cmp_low = 2, low(b) - low(c) = 0.5 ulp so far)
         if cmp_low = 0, so far low(b) - low(c) = 0 */

      /* get next limbs */
      bb = (bn > 0) ? bp[--bn] : 0;
      if ((cn > 0) && (cn-- <= cn0))
        cc = cp[cn];
      else
        cc = 0;

      /* cmp_low compares low(b) and low(c) */
      if (cmp_low == 0) /* case 1 or 3 */
        cmp_low = (bb < cc) ? -2+k : (bb > cc) ? 1 : 0;

      /* Case 1 for k=0 splits into 7 subcases:
         1a: bb > cc + half
         1b: bb = cc + half
         1c: 0 < bb - cc < half
         1d: bb = cc
         1e: -half < bb - cc < 0
         1f: bb - cc = -half
         1g: bb - cc < -half

         Case 2 splits into 3 subcases:
         2a: bb > cc
         2b: bb = cc
         2c: bb < cc

         Case 3 splits into 3 subcases:
         3a: bb > cc
         3b: bb = cc
         3c: bb < cc
      */

      /* the case rounding to nearest with sh=0 is special since one couldn't
         subtract above 1/2 ulp in the trailing limb of the result */
      if (rnd_mode == MPFR_RNDN && sh == 0 && k == 0) /* case 1 for k=0 */
        {
          mp_limb_t half = MPFR_LIMB_HIGHBIT;

          /* add one ulp if bb > cc + half
             truncate if cc - half < bb < cc + half
             sub one ulp if bb < cc - half
          */

          if (cmp_low < 0) /* bb < cc: -1 ulp < low(b) - low(c) < 0,
                              cases 1e, 1f and 1g */
            {
              if (cc >= half)
                cc -= half;
              else /* since bb < cc < half, bb+half < 2*half */
                bb += half;
              /* now we have bb < cc + half:
                 we have to subtract one ulp if bb < cc,
                 and truncate if bb > cc */
            }
          else if (cmp_low >= 0) /* bb >= cc, cases 1a to 1d */
            {
              if (cc < half)
                cc += half;
              else /* since bb >= cc >= half, bb - half >= 0 */
                bb -= half;
              /* now we have bb > cc - half: we have to add one ulp if bb > cc,
                 and truncate if bb < cc */
              if (cmp_low > 0)
                cmp_low = 2;
            }
        }

#ifdef DEBUG
      printf ("k=%u bb=%lu cc=%lu cmp_low=%d\n", k,
              (unsigned long) bb, (unsigned long) cc, cmp_low);
#endif
      if (cmp_low < 0) /* low(b) - low(c) < 0: either truncate or subtract
                          one ulp */
        {
          if (rnd_mode == MPFR_RNDZ)
            goto sub_one_ulp; /* set inexact=-1 */
          else if (rnd_mode != MPFR_RNDN) /* round away */
            {
              inexact = 1;
              goto truncate;
            }
          else /* round to nearest */
            {
              /* If cmp_low < 0 and bb > cc, then -0.5 ulp < low(b)-low(c) < 0,
                 whatever the value of sh.
                 If sh>0, then cmp_low < 0 implies that the initial neglected
                 sh bits were 0 (otherwise cmp_low=2 initially), thus the
                 weight of the new bits is less than 0.5 ulp too.
                 If k > 0 (and sh=0) this means that either the first neglected
                 limbs bb and cc were equal (thus cmp_low was 0 for k=0),
                 or we had bb - cc = -0.5 ulp or 0.5 ulp.
                 The last case is not possible here since we would have
                 cmp_low > 0 which is sticky.
                 In the first case (where we have cmp_low = -1), we truncate,
                 whereas in the 2nd case we have cmp_low = -2 and we subtract
                 one ulp.
              */
              if (bb > cc || sh > 0 || cmp_low == -1)
                {  /* -0.5 ulp < low(b)-low(c) < 0,
                      bb > cc corresponds to cases 1e and 1f1
                      sh > 0 corresponds to cases 3c and 3b3
                      cmp_low = -1 corresponds to case 1d3 (also 3b3) */
                  inexact = 1;
                  goto truncate;
                }
              else if (bb < cc) /* here sh = 0 and low(b)-low(c) < -0.5 ulp,
                                   this corresponds to cases 1g and 1f3 */
                goto sub_one_ulp;
              /* the only case where we can't conclude is sh=0 and bb=cc,
                 i.e., we have low(b) - low(c) = -0.5 ulp (up to now), thus
                 we don't know if we must truncate or subtract one ulp.
                 Note: for sh=0 we can't have low(b) - low(c) = -0.5 ulp up to
                 now, since low(b) - low(c) > 1/2^sh */
            }
        }
      else if (cmp_low > 0) /* 0 < low(b) - low(c): either truncate or
                               add one ulp */
        {
          if (rnd_mode == MPFR_RNDZ)
            {
              inexact = -1;
              goto truncate;
            }
          else if (rnd_mode != MPFR_RNDN) /* round away */
            goto add_one_ulp;
          else /* round to nearest */
            {
              if (bb > cc)
                {
                  /* if sh=0, then bb>cc means that low(b)-low(c) > 0.5 ulp,
                     and similarly when cmp_low=2 */
                  if (cmp_low == 2) /* cases 1a, 1b1, 2a and 2b1 */
                    goto add_one_ulp;
                  /* sh > 0 and cmp_low > 0: this implies that the sh initial
                     neglected bits were 0, and the remaining low(b)-low(c)>0,
                     but its weight is less than 0.5 ulp */
                  else /* 0 < low(b) - low(c) < 0.5 ulp, this corresponds to
                          cases 3a, 1d1 and 3b1 */
                    {
                      inexact = -1;
                      goto truncate;
                    }
                }
              else if (bb < cc) /* 0 < low(b) - low(c) < 0.5 ulp, cases 1c,
                                   1b3, 2b3 and 2c */
                {
                  inexact = -1;
                  goto truncate;
                }
              /* the only case where we can't conclude is bb=cc, i.e.,
                 low(b) - low(c) = 0.5 ulp (up to now), thus we don't know
                 if we must truncate or add one ulp. */
            }
        }
      /* after k=0, we cannot conclude in the following cases, we split them
         according to the values of bb and cc for k=1:
         1b. sh=0 and cmp_low = 1 and bb-cc = half [around 0.5 ulp]
             1b1. bb > cc: add one ulp, inex = 1
             1b2: bb = cc: cannot conclude
             1b3: bb < cc: truncate, inex = -1
         1d. sh=0 and cmp_low = 0 and bb-cc = 0 [around 0]
             1d1: bb > cc: truncate, inex = -1
             1d2: bb = cc: cannot conclude
             1d3: bb < cc: truncate, inex = +1
         1f. sh=0 and cmp_low = -1 and bb-cc = -half [around -0.5 ulp]
             1f1: bb > cc: truncate, inex = +1
             1f2: bb = cc: cannot conclude
             1f3: bb < cc: sub one ulp, inex = -1
         2b. sh > 0 and cmp_low = 2 and bb=cc [around 0.5 ulp]
             2b1. bb > cc: add one ulp, inex = 1
             2b2: bb = cc: cannot conclude
             2b3: bb < cc: truncate, inex = -1
         3b. sh > 0 and cmp_low = 0 [around 0]
             3b1. bb > cc: truncate, inex = -1
             3b2: bb = cc: cannot conclude
             3b3: bb < cc: truncate, inex = +1
      */
    }

  if ((rnd_mode == MPFR_RNDN) && cmp_low != 0)
    {
      /* even rounding rule */
      if ((ap[0] >> sh) & 1)
        {
          if (cmp_low < 0)
            goto sub_one_ulp;
          else
            goto add_one_ulp;
        }
      else
        inexact = (cmp_low > 0) ? -1 : 1;
    }
  else
    inexact = 0;
  goto truncate;

 sub_one_ulp: /* sub one unit in last place to a */
  mpn_sub_1 (ap, ap, an, MPFR_LIMB_ONE << sh);
  inexact = -1;
  goto end_of_sub;

 add_one_ulp: /* add one unit in last place to a */
  if (MPFR_UNLIKELY(mpn_add_1 (ap, ap, an, MPFR_LIMB_ONE << sh)))
    /* result is a power of 2: 11111111111111 + 1 = 1000000000000000 */
    {
      ap[an-1] = MPFR_LIMB_HIGHBIT;
      add_exp = 1;
    }
  inexact = 1; /* result larger than exact value */

 truncate:
  if (MPFR_UNLIKELY((ap[an-1] >> (GMP_NUMB_BITS - 1)) == 0))
    /* case 1 - epsilon */
    {
      ap[an-1] = MPFR_LIMB_HIGHBIT;
      add_exp = 1;
    }

 end_of_sub:
  /* we have to set MPFR_EXP(a) to MPFR_EXP(b) - cancel + add_exp, taking
     care of underflows/overflows in that computation, and of the allowed
     exponent range */
  if (MPFR_LIKELY(cancel))
    {
      mpfr_exp_t exp_a;

      cancel -= add_exp; /* OK: add_exp is an int equal to 0 or 1 */
      exp_a = MPFR_GET_EXP (b) - cancel;
      if (MPFR_UNLIKELY(exp_a < __gmpfr_emin))
        {
          MPFR_TMP_FREE(marker);
          if (rnd_mode == MPFR_RNDN &&
              (exp_a < __gmpfr_emin - 1 ||
               (inexact >= 0 && mpfr_powerof2_raw (a))))
            rnd_mode = MPFR_RNDZ;
          return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
        }
      MPFR_SET_EXP (a, exp_a);
    }
  else /* cancel = 0: MPFR_EXP(a) <- MPFR_EXP(b) + add_exp */
    {
      /* in case cancel = 0, add_exp can still be 1, in case b is just
         below a power of two, c is very small, prec(a) < prec(b),
         and rnd=away or nearest */
      mpfr_exp_t exp_b;

      exp_b = MPFR_GET_EXP (b);
      if (MPFR_UNLIKELY(add_exp && exp_b == __gmpfr_emax))
        {
          MPFR_TMP_FREE(marker);
          return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
        }
      MPFR_SET_EXP (a, exp_b + add_exp);
    }
  MPFR_TMP_FREE(marker);
#ifdef DEBUG
  printf ("result is a="); mpfr_print_binary(a); putchar('\n');
#endif
  /* check that result is msb-normalized */
  MPFR_ASSERTD(ap[an-1] > ~ap[an-1]);
  MPFR_RET (inexact * MPFR_INT_SIGN (a));
}