contrib/mpfr/src/add1.c

   1 /* mpfr_add1 -- internal function to perform a "real" addition
   2
   3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Caramel 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 3 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.LESSER.  If not, see
  20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  22
  23 #include "mpfr-impl.h"
  24
  25 /* compute sign(b) * (|b| + |c|), assuming b and c have same sign,
  26    and are not NaN, Inf, nor zero. */
  27 int
  28 mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
  29 {
  30   mp_limb_t *ap, *bp, *cp;
  31   mpfr_prec_t aq, bq, cq, aq2;
  32   mp_size_t an, bn, cn;
  33   mpfr_exp_t difw, exp;
  34   int sh, rb, fb, inex;
  35   mpfr_uexp_t diff_exp;
  36   MPFR_TMP_DECL(marker);
  37
  38   MPFR_ASSERTD(MPFR_IS_PURE_FP(b));
  39   MPFR_ASSERTD(MPFR_IS_PURE_FP(c));
  40
  41   MPFR_TMP_MARK(marker);
  42
  43   aq = MPFR_PREC(a);
  44   bq = MPFR_PREC(b);
  45   cq = MPFR_PREC(c);
  46
  47   an = (aq-1)/GMP_NUMB_BITS+1; /* number of limbs of a */
  48   aq2 = (mpfr_prec_t) an * GMP_NUMB_BITS;
  49   sh = aq2 - aq;                  /* non-significant bits in low limb */
  50
  51   bn = (bq-1)/GMP_NUMB_BITS+1; /* number of limbs of b */
  52   cn = (cq-1)/GMP_NUMB_BITS+1; /* number of limbs of c */
  53
  54   ap = MPFR_MANT(a);
  55   bp = MPFR_MANT(b);
  56   cp = MPFR_MANT(c);
  57
  58   if (MPFR_UNLIKELY(ap == bp))
  59     {
  60       bp = MPFR_TMP_LIMBS_ALLOC (bn);
  61       MPN_COPY (bp, ap, bn);
  62       if (ap == cp)
  63         { cp = bp; }
  64     }
  65   else if (MPFR_UNLIKELY(ap == cp))
  66     {
  67       cp = MPFR_TMP_LIMBS_ALLOC (cn);
  68       MPN_COPY(cp, ap, cn);
  69     }
  70
  71   exp = MPFR_GET_EXP (b);
  72   MPFR_SET_SAME_SIGN(a, b);
  73   MPFR_UPDATE2_RND_MODE(rnd_mode, MPFR_SIGN(b));
  74   /* now rnd_mode is either MPFR_RNDN, MPFR_RNDZ or MPFR_RNDA */
  75   /* Note: exponents can be negative, but the unsigned subtraction is
  76      a modular subtraction, so that one gets the correct result. */
  77   diff_exp = (mpfr_uexp_t) exp - MPFR_GET_EXP(c);
  78
  79   /*
  80    * 1. Compute the significant part A', the non-significant bits of A
  81    * are taken into account.
  82    *
  83    * 2. Perform the rounding. At each iteration, we remember:
  84    *     _ r = rounding bit
  85    *     _ f = following bits (same value)
  86    * where the result has the form: [number A]rfff...fff + a remaining
  87    * value in the interval [0,2) ulp. We consider the most significant
  88    * bits of the remaining value to update the result; a possible carry
  89    * is immediately taken into account and A is updated accordingly. As
  90    * soon as the bits f don't have the same value, A can be rounded.
  91    * Variables:
  92    *     _ rb = rounding bit (0 or 1).
  93    *     _ fb = following bits (0 or 1), then sticky bit.
  94    * If fb == 0, the only thing that can change is the sticky bit.
  95    */
  96
  97   rb = fb = -1; /* means: not initialized */
  98
  99   if (MPFR_UNLIKELY (MPFR_UEXP (aq2) <= diff_exp))
 100     { /* c does not overlap with a' */
 101       if (MPFR_UNLIKELY(an > bn))
 102         { /* a has more limbs than b */
 103           /* copy b to the most significant limbs of a */
 104           MPN_COPY(ap + (an - bn), bp, bn);
 105           /* zero the least significant limbs of a */
 106           MPN_ZERO(ap, an - bn);
 107         }
 108       else /* an <= bn */
 109         {
 110           /* copy the most significant limbs of b to a */
 111           MPN_COPY(ap, bp + (bn - an), an);
 112         }
 113     }
 114   else /* aq2 > diff_exp */
 115     { /* c overlaps with a' */
 116       mp_limb_t *a2p;
 117       mp_limb_t cc;
 118       mpfr_prec_t dif;
 119       mp_size_t difn, k;
 120       int shift;
 121
 122       /* copy c (shifted) into a */
 123
 124       dif = aq2 - diff_exp;
 125       /* dif is the number of bits of c which overlap with a' */
 126
 127       difn = (dif-1)/GMP_NUMB_BITS + 1;
 128       /* only the highest difn limbs from c have to be considered */
 129       if (MPFR_UNLIKELY(difn > cn))
 130         {
 131           /* c doesn't have enough limbs; take into account the virtual
 132              zero limbs now by zeroing the least significant limbs of a' */
 133           MPFR_ASSERTD(difn - cn <= an);
 134           MPN_ZERO(ap, difn - cn);
 135           difn = cn;
 136         }
 137       k = diff_exp / GMP_NUMB_BITS;
 138
 139       /* zero the most significant k limbs of a */
 140       a2p = ap + (an - k);
 141       MPN_ZERO(a2p, k);
 142
 143       shift = diff_exp % GMP_NUMB_BITS;
 144
 145       if (MPFR_LIKELY(shift))
 146         {
 147           MPFR_ASSERTD(a2p - difn >= ap);
 148           cc = mpn_rshift(a2p - difn, cp + (cn - difn), difn, shift);
 149           if (MPFR_UNLIKELY(a2p - difn > ap))
 150             *(a2p - difn - 1) = cc;
 151         }
 152       else
 153         MPN_COPY(a2p - difn, cp + (cn - difn), difn);
 154
 155       /* add b to a */
 156       cc = MPFR_UNLIKELY(an > bn)
 157         ? mpn_add_n(ap + (an - bn), ap + (an - bn), bp, bn)
 158         : mpn_add_n(ap, ap, bp + (bn - an), an);
 159
 160       if (MPFR_UNLIKELY(cc)) /* carry */
 161         {
 162           if (MPFR_UNLIKELY(exp == __gmpfr_emax))
 163             {
 164               inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
 165               goto end_of_add;
 166             }
 167           exp++;
 168           rb = (ap[0] >> sh) & 1; /* LSB(a) --> rounding bit after the shift */
 169           if (MPFR_LIKELY(sh))
 170             {
 171               mp_limb_t mask, bb;
 172
 173               mask = MPFR_LIMB_MASK (sh);
 174               bb = ap[0] & mask;
 175               ap[0] &= (~mask) << 1;
 176               if (bb == 0)
 177                 fb = 0;
 178               else if (bb == mask)
 179                 fb = 1;
 180             }
 181           mpn_rshift(ap, ap, an, 1);
 182           ap[an-1] += MPFR_LIMB_HIGHBIT;
 183           if (sh && fb < 0)
 184             goto rounding;
 185         } /* cc */
 186     } /* aq2 > diff_exp */
 187
 188   /* non-significant bits of a */
 189   if (MPFR_LIKELY(rb < 0 && sh))
 190     {
 191       mp_limb_t mask, bb;
 192
 193       mask = MPFR_LIMB_MASK (sh);
 194       bb = ap[0] & mask;
 195       ap[0] &= ~mask;
 196       rb = bb >> (sh - 1);
 197       if (MPFR_LIKELY(sh > 1))
 198         {
 199           mask >>= 1;
 200           bb &= mask;
 201           if (bb == 0)
 202             fb = 0;
 203           else if (bb == mask)
 204             fb = 1;
 205           else
 206             goto rounding;
 207         }
 208     }
 209
 210   /* determine rounding and sticky bits (and possible carry) */
 211
 212   difw = (mpfr_exp_t) an - (mpfr_exp_t) (diff_exp / GMP_NUMB_BITS);
 213   /* difw is the number of limbs from b (regarded as having an infinite
 214      precision) that have already been combined with c; -n if the next
 215      n limbs from b won't be combined with c. */
 216
 217   if (MPFR_UNLIKELY(bn > an))
 218     { /* there are still limbs from b that haven't been taken into account */
 219       mp_size_t bk;
 220
 221       if (fb == 0 && difw <= 0)
 222         {
 223           fb = 1; /* c hasn't been taken into account ==> sticky bit != 0 */
 224           goto rounding;
 225         }
 226
 227       bk = bn - an; /* index of lowest considered limb from b, > 0 */
 228       while (difw < 0)
 229         { /* ulp(next limb from b) > msb(c) */
 230           mp_limb_t bb;
 231
 232           bb = bp[--bk];
 233
 234           MPFR_ASSERTD(fb != 0);
 235           if (fb > 0)
 236             {
 237               if (bb != MP_LIMB_T_MAX)
 238                 {
 239                   fb = 1; /* c hasn't been taken into account
 240                              ==> sticky bit != 0 */
 241                   goto rounding;
 242                 }
 243             }
 244           else /* fb not initialized yet */
 245             {
 246               if (rb < 0) /* rb not initialized yet */
 247                 {
 248                   rb = bb >> (GMP_NUMB_BITS - 1);
 249                   bb |= MPFR_LIMB_HIGHBIT;
 250                 }
 251               fb = 1;
 252               if (bb != MP_LIMB_T_MAX)
 253                 goto rounding;
 254             }
 255
 256           if (bk == 0)
 257             { /* b has entirely been read */
 258               fb = 1; /* c hasn't been taken into account
 259                          ==> sticky bit != 0 */
 260               goto rounding;
 261             }
 262
 263           difw++;
 264         } /* while */
 265       MPFR_ASSERTD(bk > 0 && difw >= 0);
 266
 267       if (difw <= cn)
 268         {
 269           mp_size_t ck;
 270           mp_limb_t cprev;
 271           int difs;
 272
 273           ck = cn - difw;
 274           difs = diff_exp % GMP_NUMB_BITS;
 275
 276           if (difs == 0 && ck == 0)
 277             goto c_read;
 278
 279           cprev = ck == cn ? 0 : cp[ck];
 280
 281           if (fb < 0)
 282             {
 283               mp_limb_t bb, cc;
 284
 285               if (difs)
 286                 {
 287                   cc = cprev << (GMP_NUMB_BITS - difs);
 288                   if (--ck >= 0)
 289                     {
 290                       cprev = cp[ck];
 291                       cc += cprev >> difs;
 292                     }
 293                 }
 294               else
 295                 cc = cp[--ck];
 296
 297               bb = bp[--bk] + cc;
 298
 299               if (bb < cc /* carry */
 300                   && (rb < 0 || (rb ^= 1) == 0)
 301                   && mpn_add_1(ap, ap, an, MPFR_LIMB_ONE << sh))
 302                 {
 303                   if (exp == __gmpfr_emax)
 304                     {
 305                       inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
 306                       goto end_of_add;
 307                     }
 308                   exp++;
 309                   ap[an-1] = MPFR_LIMB_HIGHBIT;
 310                   rb = 0;
 311                 }
 312
 313               if (rb < 0) /* rb not initialized yet */
 314                 {
 315                   rb = bb >> (GMP_NUMB_BITS - 1);
 316                   bb <<= 1;
 317                   bb |= bb >> (GMP_NUMB_BITS - 1);
 318                 }
 319
 320               fb = bb != 0;
 321               if (fb && bb != MP_LIMB_T_MAX)
 322                 goto rounding;
 323             } /* fb < 0 */
 324
 325           while (bk > 0)
 326             {
 327               mp_limb_t bb, cc;
 328
 329               if (difs)
 330                 {
 331                   if (ck < 0)
 332                     goto c_read;
 333                   cc = cprev << (GMP_NUMB_BITS - difs);
 334                   if (--ck >= 0)
 335                     {
 336                       cprev = cp[ck];
 337                       cc += cprev >> difs;
 338                     }
 339                 }
 340               else
 341                 {
 342                   if (ck == 0)
 343                     goto c_read;
 344                   cc = cp[--ck];
 345                 }
 346
 347               bb = bp[--bk] + cc;
 348               if (bb < cc) /* carry */
 349                 {
 350                   fb ^= 1;
 351                   if (fb)
 352                     goto rounding;
 353                   rb ^= 1;
 354                   if (rb == 0 && mpn_add_1(ap, ap, an, MPFR_LIMB_ONE << sh))
 355                     {
 356                       if (MPFR_UNLIKELY(exp == __gmpfr_emax))
 357                         {
 358                           inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
 359                           goto end_of_add;
 360                         }
 361                       exp++;
 362                       ap[an-1] = MPFR_LIMB_HIGHBIT;
 363                     }
 364                 } /* bb < cc */
 365
 366               if (!fb && bb != 0)
 367                 {
 368                   fb = 1;
 369                   goto rounding;
 370                 }
 371               if (fb && bb != MP_LIMB_T_MAX)
 372                 goto rounding;
 373             } /* while */
 374
 375           /* b has entirely been read */
 376
 377           if (fb || ck < 0)
 378             goto rounding;
 379           if (difs && cprev << (GMP_NUMB_BITS - difs))
 380             {
 381               fb = 1;
 382               goto rounding;
 383             }
 384           while (ck)
 385             {
 386               if (cp[--ck])
 387                 {
 388                   fb = 1;
 389                   goto rounding;
 390                 }
 391             } /* while */
 392         } /* difw <= cn */
 393       else
 394         { /* c has entirely been read */
 395         c_read:
 396           if (fb < 0) /* fb not initialized yet */
 397             {
 398               mp_limb_t bb;
 399
 400               MPFR_ASSERTD(bk > 0);
 401               bb = bp[--bk];
 402               if (rb < 0) /* rb not initialized yet */
 403                 {
 404                   rb = bb >> (GMP_NUMB_BITS - 1);
 405                   bb &= ~MPFR_LIMB_HIGHBIT;
 406                 }
 407               fb = bb != 0;
 408             } /* fb < 0 */
 409           if (fb)
 410             goto rounding;
 411           while (bk)
 412             {
 413               if (bp[--bk])
 414                 {
 415                   fb = 1;
 416                   goto rounding;
 417                 }
 418             } /* while */
 419         } /* difw > cn */
 420     } /* bn > an */
 421   else if (fb != 1) /* if fb == 1, the sticky bit is 1 (no possible carry) */
 422     { /* b has entirely been read */
 423       if (difw > cn)
 424         { /* c has entirely been read */
 425           if (rb < 0)
 426             rb = 0;
 427           fb = 0;
 428         }
 429       else if (diff_exp > MPFR_UEXP (aq2))
 430         { /* b is followed by at least a zero bit, then by c */
 431           if (rb < 0)
 432             rb = 0;
 433           fb = 1;
 434         }
 435       else
 436         {
 437           mp_size_t ck;
 438           int difs;
 439
 440           MPFR_ASSERTD(difw >= 0 && cn >= difw);
 441           ck = cn - difw;
 442           difs = diff_exp % GMP_NUMB_BITS;
 443
 444           if (difs == 0 && ck == 0)
 445             { /* c has entirely been read */
 446               if (rb < 0)
 447                 rb = 0;
 448               fb = 0;
 449             }
 450           else
 451             {
 452               mp_limb_t cc;
 453
 454               cc = difs ? (MPFR_ASSERTD(ck < cn),
 455                            cp[ck] << (GMP_NUMB_BITS - difs)) : cp[--ck];
 456               if (rb < 0)
 457                 {
 458                   rb = cc >> (GMP_NUMB_BITS - 1);
 459                   cc &= ~MPFR_LIMB_HIGHBIT;
 460                 }
 461               while (cc == 0)
 462                 {
 463                   if (ck == 0)
 464                     {
 465                       fb = 0;
 466                       goto rounding;
 467                     }
 468                   cc = cp[--ck];
 469                 } /* while */
 470               fb = 1;
 471             }
 472         }
 473     } /* fb != 1 */
 474
 475  rounding:
 476   /* rnd_mode should be one of MPFR_RNDN, MPFR_RNDZ or MPFR_RNDA */
 477   if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
 478     {
 479       if (fb == 0)
 480         {
 481           if (rb == 0)
 482             {
 483               inex = 0;
 484               goto set_exponent;
 485             }
 486           /* round to even */
 487           if (ap[0] & (MPFR_LIMB_ONE << sh))
 488             goto rndn_away;
 489           else
 490             goto rndn_zero;
 491         }
 492       if (rb == 0)
 493         {
 494         rndn_zero:
 495           inex = MPFR_IS_NEG(a) ? 1 : -1;
 496           goto set_exponent;
 497         }
 498       else
 499         {
 500         rndn_away:
 501           inex = MPFR_IS_POS(a) ? 1 : -1;
 502           goto add_one_ulp;
 503         }
 504     }
 505   else if (rnd_mode == MPFR_RNDZ)
 506     {
 507       inex = rb || fb ? (MPFR_IS_NEG(a) ? 1 : -1) : 0;
 508       goto set_exponent;
 509     }
 510   else
 511     {
 512       MPFR_ASSERTN (rnd_mode == MPFR_RNDA);
 513       inex = rb || fb ? (MPFR_IS_POS(a) ? 1 : -1) : 0;
 514       if (inex)
 515         goto add_one_ulp;
 516       else
 517         goto set_exponent;
 518     }
 519
 520  add_one_ulp: /* add one unit in last place to a */
 521   if (MPFR_UNLIKELY(mpn_add_1 (ap, ap, an, MPFR_LIMB_ONE << sh)))
 522     {
 523       if (MPFR_UNLIKELY(exp == __gmpfr_emax))
 524         {
 525           inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
 526           goto end_of_add;
 527         }
 528       exp++;
 529       ap[an-1] = MPFR_LIMB_HIGHBIT;
 530     }
 531
 532  set_exponent:
 533   MPFR_SET_EXP (a, exp);
 534
 535  end_of_add:
 536   MPFR_TMP_FREE(marker);
 537   MPFR_RET (inex);
 538 }