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