contrib/mpfr/src/add1sp.c

   1 /* mpfr_add1sp -- internal function to perform a "real" addition
   2    All the op must have the same precision
   3
   4 Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   5 Contributed by the Arenaire and Caramel projects, INRIA.
   6
   7 This file is part of the GNU MPFR Library.
   8
   9 The GNU MPFR Library is free software; you can redistribute it and/or modify
  10 it under the terms of the GNU Lesser General Public License as published by
  11 the Free Software Foundation; either version 3 of the License, or (at your
  12 option) any later version.
  13
  14 The GNU MPFR Library is distributed in the hope that it will be useful, but
  15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  16 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  17 License for more details.
  18
  19 You should have received a copy of the GNU Lesser General Public License
  20 along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
  21 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  23
  24 #define MPFR_NEED_LONGLONG_H
  25 #include "mpfr-impl.h"
  26
  27 /* Check if we have to check the result of mpfr_add1sp with mpfr_add1 */
  28 #ifdef WANT_ASSERT
  29 # if WANT_ASSERT >= 2
  30
  31 int mpfr_add1sp2 (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t);
  32 int mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
  33 {
  34   mpfr_t tmpa, tmpb, tmpc;
  35   int inexb, inexc, inexact, inexact2;
  36
  37   mpfr_init2 (tmpa, MPFR_PREC (a));
  38   mpfr_init2 (tmpb, MPFR_PREC (b));
  39   mpfr_init2 (tmpc, MPFR_PREC (c));
  40
  41   inexb = mpfr_set (tmpb, b, MPFR_RNDN);
  42   MPFR_ASSERTN (inexb == 0);
  43
  44   inexc = mpfr_set (tmpc, c, MPFR_RNDN);
  45   MPFR_ASSERTN (inexc == 0);
  46
  47   inexact2 = mpfr_add1 (tmpa, tmpb, tmpc, rnd_mode);
  48   inexact  = mpfr_add1sp2 (a, b, c, rnd_mode);
  49
  50   if (mpfr_cmp (tmpa, a) || inexact != inexact2)
  51     {
  52       fprintf (stderr, "add1 & add1sp return different values for %s\n"
  53                "Prec_a = %lu, Prec_b = %lu, Prec_c = %lu\nB = ",
  54                mpfr_print_rnd_mode (rnd_mode),
  55                (unsigned long) MPFR_PREC (a),
  56                (unsigned long) MPFR_PREC (b),
  57                (unsigned long) MPFR_PREC (c));
  58       mpfr_fprint_binary (stderr, tmpb);
  59       fprintf (stderr, "\nC = ");
  60       mpfr_fprint_binary (stderr, tmpc);
  61       fprintf (stderr, "\n\nadd1  : ");
  62       mpfr_fprint_binary (stderr, tmpa);
  63       fprintf (stderr, "\nadd1sp: ");
  64       mpfr_fprint_binary (stderr, a);
  65       fprintf (stderr, "\nInexact sp = %d | Inexact = %d\n",
  66                inexact, inexact2);
  67       MPFR_ASSERTN (0);
  68     }
  69   mpfr_clears (tmpa, tmpb, tmpc, (mpfr_ptr) 0);
  70   return inexact;
  71 }
  72 #  define mpfr_add1sp mpfr_add1sp2
  73 # endif
  74 #endif
  75
  76 /* Debugging support */
  77 #ifdef DEBUG
  78 # undef DEBUG
  79 # define DEBUG(x) (x)
  80 #else
  81 # define DEBUG(x) /**/
  82 #endif
  83
  84 /* compute sign(b) * (|b| + |c|)
  85    Returns 0 iff result is exact,
  86    a negative value when the result is less than the exact value,
  87    a positive value otherwise. */
  88 int
  89 mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
  90 {
  91   mpfr_uexp_t d;
  92   mpfr_prec_t p;
  93   unsigned int sh;
  94   mp_size_t n;
  95   mp_limb_t *ap, *cp;
  96   mpfr_exp_t bx;
  97   mp_limb_t limb;
  98   int inexact;
  99   MPFR_TMP_DECL(marker);
 100
 101   MPFR_TMP_MARK(marker);
 102
 103   MPFR_ASSERTD(MPFR_PREC(a) == MPFR_PREC(b) && MPFR_PREC(b) == MPFR_PREC(c));
 104   MPFR_ASSERTD(MPFR_IS_PURE_FP(b));
 105   MPFR_ASSERTD(MPFR_IS_PURE_FP(c));
 106   MPFR_ASSERTD(MPFR_GET_EXP(b) >= MPFR_GET_EXP(c));
 107
 108   /* Read prec and num of limbs */
 109   p = MPFR_PREC(b);
 110   n = (p+GMP_NUMB_BITS-1)/GMP_NUMB_BITS;
 111   MPFR_UNSIGNED_MINUS_MODULO(sh, p);
 112   bx = MPFR_GET_EXP(b);
 113   d = (mpfr_uexp_t) (bx - MPFR_GET_EXP(c));
 114
 115   DEBUG (printf ("New add1sp with diff=%lu\n", (unsigned long) d));
 116
 117   if (MPFR_UNLIKELY(d == 0))
 118     {
 119       /* d==0 */
 120       DEBUG( mpfr_print_mant_binary("C= ", MPFR_MANT(c), p) );
 121       DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b), p) );
 122       bx++;                                /* exp + 1 */
 123       ap = MPFR_MANT(a);
 124       limb = mpn_add_n(ap, MPFR_MANT(b), MPFR_MANT(c), n);
 125       DEBUG( mpfr_print_mant_binary("A= ", ap, p) );
 126       MPFR_ASSERTD(limb != 0);             /* There must be a carry */
 127       limb = ap[0];                        /* Get LSB (In fact, LSW) */
 128       mpn_rshift(ap, ap, n, 1);            /* Shift mantissa A */
 129       ap[n-1] |= MPFR_LIMB_HIGHBIT;        /* Set MSB */
 130       ap[0]   &= ~MPFR_LIMB_MASK(sh);      /* Clear LSB bit */
 131       if (MPFR_LIKELY((limb&(MPFR_LIMB_ONE<<sh)) == 0)) /* Check exact case */
 132         { inexact = 0; goto set_exponent; }
 133       /* Zero: Truncate
 134          Nearest: Even Rule => truncate or add 1
 135          Away: Add 1 */
 136       if (MPFR_LIKELY(rnd_mode==MPFR_RNDN))
 137         {
 138           if (MPFR_LIKELY((ap[0]&(MPFR_LIMB_ONE<<sh))==0))
 139             { inexact = -1; goto set_exponent; }
 140           else
 141             goto add_one_ulp;
 142         }
 143       MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b));
 144       if (rnd_mode==MPFR_RNDZ)
 145         { inexact = -1; goto set_exponent; }
 146       else
 147         goto add_one_ulp;
 148     }
 149   else if (MPFR_UNLIKELY (d >= p))
 150     {
 151       if (MPFR_LIKELY (d > p))
 152         {
 153           /* d > p : Copy B in A */
 154           /* Away:    Add 1
 155              Nearest: Trunc
 156              Zero:    Trunc */
 157           if (MPFR_LIKELY (rnd_mode==MPFR_RNDN
 158                            || MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (b))))
 159             {
 160             copy_set_exponent:
 161               ap = MPFR_MANT (a);
 162               MPN_COPY (ap, MPFR_MANT(b), n);
 163               inexact = -1;
 164               goto set_exponent;
 165             }
 166           else
 167             {
 168             copy_add_one_ulp:
 169               ap = MPFR_MANT(a);
 170               MPN_COPY (ap, MPFR_MANT(b), n);
 171               goto add_one_ulp;
 172             }
 173         }
 174       else
 175         {
 176           /* d==p : Copy B in A */
 177           /* Away:    Add 1
 178              Nearest: Even Rule if C is a power of 2, else Add 1
 179              Zero:    Trunc */
 180           if (MPFR_LIKELY(rnd_mode==MPFR_RNDN))
 181             {
 182               /* Check if C was a power of 2 */
 183               cp = MPFR_MANT(c);
 184               if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
 185                 {
 186                   mp_size_t k = n-1;
 187                   do {
 188                     k--;
 189                   } while (k>=0 && cp[k]==0);
 190                   if (MPFR_UNLIKELY(k<0))
 191                     /* Power of 2: Even rule */
 192                     if ((MPFR_MANT (b)[0]&(MPFR_LIMB_ONE<<sh))==0)
 193                       goto copy_set_exponent;
 194                 }
 195               /* Not a Power of 2 */
 196               goto copy_add_one_ulp;
 197             }
 198           else if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (b)))
 199             goto copy_set_exponent;
 200           else
 201             goto copy_add_one_ulp;
 202         }
 203     }
 204   else
 205     {
 206       mp_limb_t mask;
 207       mp_limb_t bcp, bcp1; /* Cp and C'p+1 */
 208
 209       /* General case: 1 <= d < p */
 210       cp = MPFR_TMP_LIMBS_ALLOC (n);
 211
 212       /* Shift c in temporary allocated place */
 213       {
 214         mpfr_uexp_t dm;
 215         mp_size_t m;
 216
 217         dm = d % GMP_NUMB_BITS;
 218         m = d / GMP_NUMB_BITS;
 219         if (MPFR_UNLIKELY(dm == 0))
 220           {
 221             /* dm = 0 and m > 0: Just copy */
 222             MPFR_ASSERTD(m!=0);
 223             MPN_COPY(cp, MPFR_MANT(c)+m, n-m);
 224             MPN_ZERO(cp+n-m, m);
 225           }
 226         else if (MPFR_LIKELY(m == 0))
 227           {
 228             /* dm >=1 and m == 0: just shift */
 229             MPFR_ASSERTD(dm >= 1);
 230             mpn_rshift(cp, MPFR_MANT(c), n, dm);
 231           }
 232         else
 233           {
 234             /* dm > 0 and m > 0: shift and zero  */
 235             mpn_rshift(cp, MPFR_MANT(c)+m, n-m, dm);
 236             MPN_ZERO(cp+n-m, m);
 237           }
 238       }
 239
 240       DEBUG( mpfr_print_mant_binary("Before", MPFR_MANT(c), p) );
 241       DEBUG( mpfr_print_mant_binary("B=    ", MPFR_MANT(b), p) );
 242       DEBUG( mpfr_print_mant_binary("After ", cp, p) );
 243
 244       /* Compute bcp=Cp and bcp1=C'p+1 */
 245       if (MPFR_LIKELY (sh > 0))
 246         {
 247           /* Try to compute them from C' rather than C */
 248           bcp = (cp[0] & (MPFR_LIMB_ONE<<(sh-1))) ;
 249           if (MPFR_LIKELY(cp[0]&MPFR_LIMB_MASK(sh-1)))
 250             bcp1 = 1;
 251           else
 252             {
 253               /* We can't compute C'p+1 from C'. Compute it from C */
 254               /* Start from bit x=p-d+sh in mantissa C
 255                  (+sh since we have already looked sh bits in C'!) */
 256               mpfr_prec_t x = p-d+sh-1;
 257               if (MPFR_LIKELY(x>p))
 258                 /* We are already looked at all the bits of c, so C'p+1 = 0*/
 259                 bcp1 = 0;
 260               else
 261                 {
 262                   mp_limb_t *tp = MPFR_MANT(c);
 263                   mp_size_t kx = n-1 - (x / GMP_NUMB_BITS);
 264                   mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
 265                   DEBUG (printf ("(First) x=%lu Kx=%ld Sx=%lu\n",
 266                                  (unsigned long) x, (long) kx,
 267                                  (unsigned long) sx));
 268                   /* Looks at the last bits of limb kx (if sx=0 does nothing)*/
 269                   if (tp[kx] & MPFR_LIMB_MASK(sx))
 270                     bcp1 = 1;
 271                   else
 272                     {
 273                       /*kx += (sx==0);*/
 274                       /*If sx==0, tp[kx] hasn't been checked*/
 275                       do {
 276                         kx--;
 277                       } while (kx>=0 && tp[kx]==0);
 278                       bcp1 = (kx >= 0);
 279                     }
 280                 }
 281             }
 282         }
 283       else /* sh == 0 */
 284         {
 285           /* Compute Cp and C'p+1 from C with sh=0 */
 286           mp_limb_t *tp = MPFR_MANT(c);
 287           /* Start from bit x=p-d in mantissa C */
 288           mpfr_prec_t  x = p-d;
 289           mp_size_t   kx = n-1 - (x / GMP_NUMB_BITS);
 290           mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
 291           MPFR_ASSERTD(p >= d);
 292           bcp = tp[kx] & (MPFR_LIMB_ONE<<sx);
 293           /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
 294           if (tp[kx]&MPFR_LIMB_MASK(sx))
 295             bcp1 = 1;
 296           else
 297             {
 298               do {
 299                 kx--;
 300               } while (kx>=0 && tp[kx]==0);
 301               bcp1 = (kx>=0);
 302             }
 303         }
 304       DEBUG (printf("sh=%u Cp=%lu C'p+1=%lu\n", sh,
 305                     (unsigned long) bcp, (unsigned long) bcp1));
 306
 307       /* Clean shifted C' */
 308       mask = ~MPFR_LIMB_MASK(sh);
 309       cp[0] &= mask;
 310
 311       /* Add the mantissa c from b in a */
 312       ap = MPFR_MANT(a);
 313       limb = mpn_add_n (ap, MPFR_MANT(b), cp, n);
 314       DEBUG( mpfr_print_mant_binary("Add=  ", ap, p) );
 315
 316       /* Check for overflow */
 317       if (MPFR_UNLIKELY (limb))
 318         {
 319           limb = ap[0] & (MPFR_LIMB_ONE<<sh); /* Get LSB */
 320           mpn_rshift (ap, ap, n, 1);          /* Shift mantissa*/
 321           bx++;                               /* Fix exponent */
 322           ap[n-1] |= MPFR_LIMB_HIGHBIT;       /* Set MSB */
 323           ap[0]   &= mask;                    /* Clear LSB bit */
 324           bcp1    |= bcp;                     /* Recompute C'p+1 */
 325           bcp      = limb;                    /* Recompute Cp */
 326           DEBUG (printf ("(Overflow) Cp=%lu C'p+1=%lu\n",
 327                          (unsigned long) bcp, (unsigned long) bcp1));
 328           DEBUG (mpfr_print_mant_binary ("Add=  ", ap, p));
 329         }
 330
 331       /* Round:
 332           Zero: Truncate but could be exact.
 333           Away: Add 1 if Cp or C'p+1 !=0
 334           Nearest: Truncate but could be exact if Cp==0
 335                    Add 1 if C'p+1 !=0,
 336                    Even rule else */
 337       if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
 338         {
 339           if (MPFR_LIKELY(bcp == 0))
 340             { inexact = MPFR_LIKELY(bcp1) ? -1 : 0; goto set_exponent; }
 341           else if (MPFR_UNLIKELY(bcp1==0) && (ap[0]&(MPFR_LIMB_ONE<<sh))==0)
 342             { inexact = -1; goto set_exponent; }
 343           else
 344             goto add_one_ulp;
 345         }
 346       MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b));
 347       if (rnd_mode == MPFR_RNDZ)
 348         {
 349           inexact = MPFR_LIKELY(bcp || bcp1) ? -1 : 0;
 350           goto set_exponent;
 351         }
 352       else
 353         {
 354           if (MPFR_UNLIKELY(bcp==0 && bcp1==0))
 355             { inexact = 0; goto set_exponent; }
 356           else
 357             goto add_one_ulp;
 358         }
 359     }
 360   MPFR_ASSERTN(0);
 361
 362  add_one_ulp:
 363   /* add one unit in last place to a */
 364   DEBUG( printf("AddOneUlp\n") );
 365   if (MPFR_UNLIKELY( mpn_add_1(ap, ap, n, MPFR_LIMB_ONE<<sh) ))
 366     {
 367       /* Case 100000x0 = 0x1111x1 + 1*/
 368       DEBUG( printf("Pow of 2\n") );
 369       bx++;
 370       ap[n-1] = MPFR_LIMB_HIGHBIT;
 371     }
 372   inexact = 1;
 373
 374  set_exponent:
 375   if (MPFR_UNLIKELY(bx > __gmpfr_emax)) /* Check for overflow */
 376     {
 377       DEBUG( printf("Overflow\n") );
 378       MPFR_TMP_FREE(marker);
 379       MPFR_SET_SAME_SIGN(a,b);
 380       return mpfr_overflow(a, rnd_mode, MPFR_SIGN(a));
 381     }
 382   MPFR_SET_EXP (a, bx);
 383   MPFR_SET_SAME_SIGN(a,b);
 384
 385   MPFR_TMP_FREE(marker);
 386   MPFR_RET (inexact * MPFR_INT_SIGN (a));
 387 }