Import mpfr-2.4.1

[dragonfly.git] / contrib / mpfr / add1.c

/* mpfr_add1 -- internal function to perform a "real" addition

Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
Contributed by the Arenaire and Cacao 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 2.1 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.LIB.  If not, 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|) */
int
mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
{
  mp_limb_t *ap, *bp, *cp;
  mp_prec_t aq, bq, cq, aq2;
  mp_size_t an, bn, cn;
  mp_exp_t difw, exp;
  int sh, rb, fb, inex;
  mp_exp_unsigned_t diff_exp;
  MPFR_TMP_DECL(marker);

  MPFR_ASSERTD(MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c));

  MPFR_TMP_MARK(marker);

  aq = MPFR_PREC(a);
  bq = MPFR_PREC(b);
  cq = MPFR_PREC(c);

  an = (aq-1)/BITS_PER_MP_LIMB+1; /* number of limbs of a */
  aq2 = (mp_prec_t) an * BITS_PER_MP_LIMB;
  sh = aq2 - aq;                  /* non-significant bits in low limb */

  bn = (bq-1)/BITS_PER_MP_LIMB+1; /* number of limbs of b */
  cn = (cq-1)/BITS_PER_MP_LIMB+1; /* number of limbs of c */

  ap = MPFR_MANT(a);
  bp = MPFR_MANT(b);
  cp = MPFR_MANT(c);

  if (MPFR_UNLIKELY(ap == bp))
    {
      bp = (mp_ptr) MPFR_TMP_ALLOC (bn * BYTES_PER_MP_LIMB);
      MPN_COPY (bp, ap, bn);
      if (ap == cp)
        { cp = bp; }
    }
  else if (MPFR_UNLIKELY(ap == cp))
    {
      cp = (mp_ptr) MPFR_TMP_ALLOC (cn * BYTES_PER_MP_LIMB);
      MPN_COPY(cp, ap, cn);
    }

  exp = MPFR_GET_EXP (b);
  MPFR_SET_SAME_SIGN(a, b);
  diff_exp = (mp_exp_unsigned_t) exp - MPFR_GET_EXP(c);

  /*
   * 1. Compute the significant part A', the non-significant bits of A
   * are taken into account.
   *
   * 2. Perform the rounding. At each iteration, we remember:
   *     _ r = rounding bit
   *     _ f = following bits (same value)
   * where the result has the form: [number A]rfff...fff + a remaining
   * value in the interval [0,2) ulp. We consider the most significant
   * bits of the remaining value to update the result; a possible carry
   * is immediately taken into account and A is updated accordingly. As
   * soon as the bits f don't have the same value, A can be rounded.
   * Variables:
   *     _ rb = rounding bit (0 or 1).
   *     _ fb = following bits (0 or 1), then sticky bit.
   * If fb == 0, the only thing that can change is the sticky bit.
   */

  rb = fb = -1; /* means: not initialized */

  if (MPFR_UNLIKELY(aq2 <= diff_exp))
    { /* c does not overlap with a' */
      if (MPFR_UNLIKELY(an > bn))
        { /* a has more limbs than b */
          /* copy b to the most significant limbs of a */
          MPN_COPY(ap + (an - bn), bp, bn);
          /* zero the least significant limbs of a */
          MPN_ZERO(ap, an - bn);
        }
      else /* an <= bn */
        {
          /* copy the most significant limbs of b to a */
          MPN_COPY(ap, bp + (bn - an), an);
        }
    }
  else /* aq2 > diff_exp */
    { /* c overlaps with a' */
      mp_limb_t *a2p;
      mp_limb_t cc;
      mp_prec_t dif;
      mp_size_t difn, k;
      int shift;

      /* copy c (shifted) into a */

      dif = aq2 - diff_exp;
      /* dif is the number of bits of c which overlap with a' */

      difn = (dif-1)/BITS_PER_MP_LIMB + 1;
      /* only the highest difn limbs from c have to be considered */
      if (MPFR_UNLIKELY(difn > cn))
        {
          /* c doesn't have enough limbs; take into account the virtual
             zero limbs now by zeroing the least significant limbs of a' */
          MPFR_ASSERTD(difn - cn <= an);
          MPN_ZERO(ap, difn - cn);
          difn = cn;
        }
      k = diff_exp / BITS_PER_MP_LIMB;

      /* zero the most significant k limbs of a */
      a2p = ap + (an - k);
      MPN_ZERO(a2p, k);

      shift = diff_exp % BITS_PER_MP_LIMB;

      if (MPFR_LIKELY(shift))
        {
          MPFR_ASSERTD(a2p - difn >= ap);
          cc = mpn_rshift(a2p - difn, cp + (cn - difn), difn, shift);
          if (MPFR_UNLIKELY(a2p - difn > ap))
            *(a2p - difn - 1) = cc;
        }
      else
        MPN_COPY(a2p - difn, cp + (cn - difn), difn);

      /* add b to a */
      cc = MPFR_UNLIKELY(an > bn)
        ? mpn_add_n(ap + (an - bn), ap + (an - bn), bp, bn)
        : mpn_add_n(ap, ap, bp + (bn - an), an);

      if (MPFR_UNLIKELY(cc)) /* carry */
        {
          if (MPFR_UNLIKELY(exp == __gmpfr_emax))
            {
              inex = mpfr_overflow(a, rnd_mode, MPFR_SIGN(a));
              goto end_of_add;
            }
          exp++;
          rb = (ap[0] >> sh) & 1; /* LSB(a) --> rounding bit after the shift */
          if (MPFR_LIKELY(sh))
            {
              mp_limb_t mask, bb;

              mask = MPFR_LIMB_MASK (sh);
              bb = ap[0] & mask;
              ap[0] &= (~mask) << 1;
              if (bb == 0)
                fb = 0;
              else if (bb == mask)
                fb = 1;
            }
          mpn_rshift(ap, ap, an, 1);
          ap[an-1] += MPFR_LIMB_HIGHBIT;
          if (sh && fb < 0)
            goto rounding;
        } /* cc */
    } /* aq2 > diff_exp */

  /* non-significant bits of a */
  if (MPFR_LIKELY(rb < 0 && sh))
    {
      mp_limb_t mask, bb;

      mask = MPFR_LIMB_MASK (sh);
      bb = ap[0] & mask;
      ap[0] &= ~mask;
      rb = bb >> (sh - 1);
      if (MPFR_LIKELY(sh > 1))
        {
          mask >>= 1;
          bb &= mask;
          if (bb == 0)
            fb = 0;
          else if (bb == mask)
            fb = 1;
          else
            goto rounding;
        }
    }

  /* determine rounding and sticky bits (and possible carry) */

  difw = (mp_exp_t) an - (mp_exp_t) (diff_exp / BITS_PER_MP_LIMB);
  /* difw is the number of limbs from b (regarded as having an infinite
     precision) that have already been combined with c; -n if the next
     n limbs from b won't be combined with c. */

  if (MPFR_UNLIKELY(bn > an))
    { /* there are still limbs from b that haven't been taken into account */
      mp_size_t bk;

      if (fb == 0 && difw <= 0)
        {
          fb = 1; /* c hasn't been taken into account ==> sticky bit != 0 */
          goto rounding;
        }

      bk = bn - an; /* index of lowest considered limb from b, > 0 */
      while (difw < 0)
        { /* ulp(next limb from b) > msb(c) */
          mp_limb_t bb;

          bb = bp[--bk];

          MPFR_ASSERTD(fb != 0);
          if (fb > 0)
            {
              if (bb != MP_LIMB_T_MAX)
                {
                  fb = 1; /* c hasn't been taken into account
                             ==> sticky bit != 0 */
                  goto rounding;
                }
            }
          else /* fb not initialized yet */
            {
              if (rb < 0) /* rb not initialized yet */
                {
                  rb = bb >> (BITS_PER_MP_LIMB - 1);
                  bb |= MPFR_LIMB_HIGHBIT;
                }
              fb = 1;
              if (bb != MP_LIMB_T_MAX)
                goto rounding;
            }

          if (bk == 0)
            { /* b has entirely been read */
              fb = 1; /* c hasn't been taken into account
                         ==> sticky bit != 0 */
              goto rounding;
            }

          difw++;
        } /* while */
      MPFR_ASSERTD(bk > 0 && difw >= 0);

      if (difw <= cn)
        {
          mp_size_t ck;
          mp_limb_t cprev;
          int difs;

          ck = cn - difw;
          difs = diff_exp % BITS_PER_MP_LIMB;

          if (difs == 0 && ck == 0)
            goto c_read;

          cprev = ck == cn ? 0 : cp[ck];

          if (fb < 0)
            {
              mp_limb_t bb, cc;

              if (difs)
                {
                  cc = cprev << (BITS_PER_MP_LIMB - difs);
                  if (--ck >= 0)
                    {
                      cprev = cp[ck];
                      cc += cprev >> difs;
                    }
                }
              else
                cc = cp[--ck];

              bb = bp[--bk] + cc;

              if (bb < cc /* carry */
                  && (rb < 0 || (rb ^= 1) == 0)
                  && mpn_add_1(ap, ap, an, MPFR_LIMB_ONE << sh))
                {
                  if (exp == __gmpfr_emax)
                    {
                      inex = mpfr_overflow(a, rnd_mode, MPFR_SIGN(a));
                      goto end_of_add;
                    }
                  exp++;
                  ap[an-1] = MPFR_LIMB_HIGHBIT;
                  rb = 0;
                }

              if (rb < 0) /* rb not initialized yet */
                {
                  rb = bb >> (BITS_PER_MP_LIMB - 1);
                  bb <<= 1;
                  bb |= bb >> (BITS_PER_MP_LIMB - 1);
                }

              fb = bb != 0;
              if (fb && bb != MP_LIMB_T_MAX)
                goto rounding;
            } /* fb < 0 */

          while (bk > 0)
            {
              mp_limb_t bb, cc;

              if (difs)
                {
                  if (ck < 0)
                    goto c_read;
                  cc = cprev << (BITS_PER_MP_LIMB - difs);
                  if (--ck >= 0)
                    {
                      cprev = cp[ck];
                      cc += cprev >> difs;
                    }
                }
              else
                {
                  if (ck == 0)
                    goto c_read;
                  cc = cp[--ck];
                }

              bb = bp[--bk] + cc;
              if (bb < cc) /* carry */
                {
                  fb ^= 1;
                  if (fb)
                    goto rounding;
                  rb ^= 1;
                  if (rb == 0 && mpn_add_1(ap, ap, an, MPFR_LIMB_ONE << sh))
                    {
                      if (MPFR_UNLIKELY(exp == __gmpfr_emax))
                        {
                          inex = mpfr_overflow(a, rnd_mode, MPFR_SIGN(a));
                          goto end_of_add;
                        }
                      exp++;
                      ap[an-1] = MPFR_LIMB_HIGHBIT;
                    }
                } /* bb < cc */

              if (!fb && bb != 0)
                {
                  fb = 1;
                  goto rounding;
                }
              if (fb && bb != MP_LIMB_T_MAX)
                goto rounding;
            } /* while */

          /* b has entirely been read */

          if (fb || ck < 0)
            goto rounding;
          if (difs && cprev << (BITS_PER_MP_LIMB - difs))
            {
              fb = 1;
              goto rounding;
            }
          while (ck)
            {
              if (cp[--ck])
                {
                  fb = 1;
                  goto rounding;
                }
            } /* while */
        } /* difw <= cn */
      else
        { /* c has entirely been read */
        c_read:
          if (fb < 0) /* fb not initialized yet */
            {
              mp_limb_t bb;

              MPFR_ASSERTD(bk > 0);
              bb = bp[--bk];
              if (rb < 0) /* rb not initialized yet */
                {
                  rb = bb >> (BITS_PER_MP_LIMB - 1);
                  bb &= ~MPFR_LIMB_HIGHBIT;
                }
              fb = bb != 0;
            } /* fb < 0 */
          if (fb)
            goto rounding;
          while (bk)
            {
              if (bp[--bk])
                {
                  fb = 1;
                  goto rounding;
                }
            } /* while */
        } /* difw > cn */
    } /* bn > an */
  else if (fb != 1) /* if fb == 1, the sticky bit is 1 (no possible carry) */
    { /* b has entirely been read */
      if (difw > cn)
        { /* c has entirely been read */
          if (rb < 0)
            rb = 0;
          fb = 0;
        }
      else if (diff_exp > aq2)
        { /* b is followed by at least a zero bit, then by c */
          if (rb < 0)
            rb = 0;
          fb = 1;
        }
      else
        {
          mp_size_t ck;
          int difs;

          MPFR_ASSERTD(difw >= 0 && cn >= difw);
          ck = cn - difw;
          difs = diff_exp % BITS_PER_MP_LIMB;

          if (difs == 0 && ck == 0)
            { /* c has entirely been read */
              if (rb < 0)
                rb = 0;
              fb = 0;
            }
          else
            {
              mp_limb_t cc;

              cc = difs ? (MPFR_ASSERTD(ck < cn),
                           cp[ck] << (BITS_PER_MP_LIMB - difs)) : cp[--ck];
              if (rb < 0)
                {
                  rb = cc >> (BITS_PER_MP_LIMB - 1);
                  cc &= ~MPFR_LIMB_HIGHBIT;
                }
              while (cc == 0)
                {
                  if (ck == 0)
                    {
                      fb = 0;
                      goto rounding;
                    }
                  cc = cp[--ck];
                } /* while */
              fb = 1;
            }
        }
    } /* fb != 1 */

 rounding:
  switch (rnd_mode)
    {
    case GMP_RNDN:
      if (fb == 0)
        {
          if (rb == 0)
            {
              inex = 0;
              goto set_exponent;
            }
          /* round to even */
          if (ap[0] & (MPFR_LIMB_ONE << sh))
            goto rndn_away;
          else
            goto rndn_zero;
        }
      if (rb == 0)
        {
        rndn_zero:
          inex = MPFR_IS_NEG(a) ? 1 : -1;
          goto set_exponent;
        }
      else
        {
        rndn_away:
          inex = MPFR_IS_POS(a) ? 1 : -1;
          goto add_one_ulp;
        }

    case GMP_RNDZ:
      inex = rb || fb ? (MPFR_IS_NEG(a) ? 1 : -1) : 0;
      goto set_exponent;

    case GMP_RNDU:
      inex = rb || fb;
      if (inex && MPFR_IS_POS(a))
        goto add_one_ulp;
      else
        goto set_exponent;

    case GMP_RNDD:
      inex = - (rb || fb);
      if (inex && MPFR_IS_NEG(a))
        goto add_one_ulp;
      else
        goto set_exponent;

    default:
      MPFR_ASSERTN(0);
      inex = 0;
      goto set_exponent;
    }

 add_one_ulp: /* add one unit in last place to a */
  if (MPFR_UNLIKELY(mpn_add_1(ap, ap, an, MPFR_LIMB_ONE << sh)))
    {
      if (MPFR_UNLIKELY(exp == __gmpfr_emax))
        {
          inex = mpfr_overflow(a, rnd_mode, MPFR_SIGN(a));
          goto end_of_add;
        }
      exp++;
      ap[an-1] = MPFR_LIMB_HIGHBIT;
    }

 set_exponent:
  MPFR_SET_EXP (a, exp);

 end_of_add:
  MPFR_TMP_FREE(marker);
  MPFR_RET (inex);
}