contrib/mpfr/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 Free Software Foundation, Inc.
   5 Contributed by the Arenaire and Cacao 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 2.1 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.LIB.  If not, write to
  21 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  22 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, mp_rnd_t);
  32 int mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_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, GMP_RNDN);
  42   MPFR_ASSERTN (inexb == 0);
  43
  44   inexc = mpfr_set (tmpc, c, GMP_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                MPFR_PREC (a), MPFR_PREC (b), MPFR_PREC (c));
  56       mpfr_fprint_binary (stderr, tmpb);
  57       fprintf (stderr, "\nC = ");
  58       mpfr_fprint_binary (stderr, tmpc);
  59       fprintf (stderr, "\n\nadd1  : ");
  60       mpfr_fprint_binary (stderr, tmpa);
  61       fprintf (stderr, "\nadd1sp: ");
  62       mpfr_fprint_binary (stderr, a);
  63       fprintf (stderr, "\nInexact sp = %d | Inexact = %d\n",
  64                inexact, inexact2);
  65       MPFR_ASSERTN (0);
  66     }
  67   mpfr_clears (tmpa, tmpb, tmpc, (mpfr_ptr) 0);
  68   return inexact;
  69 }
  70 #  define mpfr_add1sp mpfr_add1sp2
  71 # endif
  72 #endif
  73
  74 /* Debugging support */
  75 #ifdef DEBUG
  76 # undef DEBUG
  77 # define DEBUG(x) (x)
  78 #else
  79 # define DEBUG(x) /**/
  80 #endif
  81
  82 /* compute sign(b) * (|b| + |c|)
  83    Returns 0 iff result is exact,
  84    a negative value when the result is less than the exact value,
  85    a positive value otherwise. */
  86 int
  87 mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode)
  88 {
  89   mp_exp_unsigned_t d;
  90   mp_prec_t p;
  91   unsigned int sh;
  92   mp_size_t n;
  93   mp_limb_t *ap, *cp;
  94   mp_exp_t  bx;
  95   mp_limb_t limb;
  96   int inexact;
  97   MPFR_TMP_DECL(marker);
  98
  99   MPFR_TMP_MARK(marker);
 100
 101   MPFR_ASSERTD(MPFR_PREC(a) == MPFR_PREC(b) && MPFR_PREC(b) == MPFR_PREC(c));
 102   MPFR_ASSERTD(MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c));
 103   MPFR_ASSERTD(MPFR_GET_EXP(b) >= MPFR_GET_EXP(c));
 104
 105   /* Read prec and num of limbs */
 106   p = MPFR_PREC(b);
 107   n = (p+BITS_PER_MP_LIMB-1)/BITS_PER_MP_LIMB;
 108   MPFR_UNSIGNED_MINUS_MODULO(sh, p);
 109   bx = MPFR_GET_EXP(b);
 110   d = (mpfr_uexp_t) (bx - MPFR_GET_EXP(c));
 111
 112   DEBUG (printf ("New add1sp with diff=%lu\n", (unsigned long) d));
 113
 114   if (MPFR_UNLIKELY(d == 0))
 115     {
 116       /* d==0 */
 117       DEBUG( mpfr_print_mant_binary("C= ", MPFR_MANT(c), p) );
 118       DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b), p) );
 119       bx++;                                /* exp + 1 */
 120       ap = MPFR_MANT(a);
 121       limb = mpn_add_n(ap, MPFR_MANT(b), MPFR_MANT(c), n);
 122       DEBUG( mpfr_print_mant_binary("A= ", ap, p) );
 123       MPFR_ASSERTD(limb != 0);             /* There must be a carry */
 124       limb = ap[0];                        /* Get LSB (In fact, LSW) */
 125       mpn_rshift(ap, ap, n, 1);            /* Shift mantissa A */
 126       ap[n-1] |= MPFR_LIMB_HIGHBIT;        /* Set MSB */
 127       ap[0]   &= ~MPFR_LIMB_MASK(sh);      /* Clear LSB bit */
 128       if (MPFR_LIKELY((limb&(MPFR_LIMB_ONE<<sh)) == 0)) /* Check exact case */
 129         { inexact = 0; goto set_exponent; }
 130       /* Zero: Truncate
 131          Nearest: Even Rule => truncate or add 1
 132          Away: Add 1 */
 133       if (MPFR_LIKELY(rnd_mode==GMP_RNDN))
 134         {
 135           if (MPFR_LIKELY((ap[0]&(MPFR_LIMB_ONE<<sh))==0))
 136             { inexact = -1; goto set_exponent; }
 137           else
 138             goto add_one_ulp;
 139         }
 140       MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b));
 141       if (rnd_mode==GMP_RNDZ)
 142         { inexact = -1; goto set_exponent; }
 143       else
 144         goto add_one_ulp;
 145     }
 146   else if (MPFR_UNLIKELY (d >= p))
 147     {
 148       if (MPFR_LIKELY (d > p))
 149         {
 150           /* d > p : Copy B in A */
 151           /* Away:    Add 1
 152              Nearest: Trunc
 153              Zero:    Trunc */
 154           if (MPFR_LIKELY (rnd_mode==GMP_RNDN
 155                            || MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (b))))
 156             {
 157             copy_set_exponent:
 158               ap = MPFR_MANT (a);
 159               MPN_COPY (ap, MPFR_MANT(b), n);
 160               inexact = -1;
 161               goto set_exponent;
 162             }
 163           else
 164             {
 165             copy_add_one_ulp:
 166               ap = MPFR_MANT(a);
 167               MPN_COPY (ap, MPFR_MANT(b), n);
 168               goto add_one_ulp;
 169             }
 170         }
 171       else
 172         {
 173           /* d==p : Copy B in A */
 174           /* Away:    Add 1
 175              Nearest: Even Rule if C is a power of 2, else Add 1
 176              Zero:    Trunc */
 177           if (MPFR_LIKELY(rnd_mode==GMP_RNDN))
 178             {
 179               /* Check if C was a power of 2 */
 180               cp = MPFR_MANT(c);
 181               if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
 182                 {
 183                   mp_size_t k = n-1;
 184                   do {
 185                     k--;
 186                   } while (k>=0 && cp[k]==0);
 187                   if (MPFR_UNLIKELY(k<0))
 188                     /* Power of 2: Even rule */
 189                     if ((MPFR_MANT (b)[0]&(MPFR_LIMB_ONE<<sh))==0)
 190                       goto copy_set_exponent;
 191                 }
 192               /* Not a Power of 2 */
 193               goto copy_add_one_ulp;
 194             }
 195           else if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (b)))
 196             goto copy_set_exponent;
 197           else
 198             goto copy_add_one_ulp;
 199         }
 200     }
 201   else
 202     {
 203       mp_limb_t mask;
 204       mp_limb_t bcp, bcp1; /* Cp and C'p+1 */
 205
 206       /* General case: 1 <= d < p */
 207       cp = (mp_limb_t*) MPFR_TMP_ALLOC(n * BYTES_PER_MP_LIMB);
 208
 209       /* Shift c in temporary allocated place */
 210       {
 211         mp_exp_unsigned_t dm;
 212         mp_size_t m;
 213
 214         dm = d % BITS_PER_MP_LIMB;
 215         m = d / BITS_PER_MP_LIMB;
 216         if (MPFR_UNLIKELY(dm == 0))
 217           {
 218             /* dm = 0 and m > 0: Just copy */
 219             MPFR_ASSERTD(m!=0);
 220             MPN_COPY(cp, MPFR_MANT(c)+m, n-m);
 221             MPN_ZERO(cp+n-m, m);
 222           }
 223         else if (MPFR_LIKELY(m == 0))
 224           {
 225             /* dm >=1 and m == 0: just shift */
 226             MPFR_ASSERTD(dm >= 1);
 227             mpn_rshift(cp, MPFR_MANT(c), n, dm);
 228           }
 229         else
 230           {
 231             /* dm > 0 and m > 0: shift and zero  */
 232             mpn_rshift(cp, MPFR_MANT(c)+m, n-m, dm);
 233             MPN_ZERO(cp+n-m, m);
 234           }
 235       }
 236
 237       DEBUG( mpfr_print_mant_binary("Before", MPFR_MANT(c), p) );
 238       DEBUG( mpfr_print_mant_binary("B=    ", MPFR_MANT(b), p) );
 239       DEBUG( mpfr_print_mant_binary("After ", cp, p) );
 240
 241       /* Compute bcp=Cp and bcp1=C'p+1 */
 242       if (MPFR_LIKELY (sh > 0))
 243         {
 244           /* Try to compute them from C' rather than C */
 245           bcp = (cp[0] & (MPFR_LIMB_ONE<<(sh-1))) ;
 246           if (MPFR_LIKELY(cp[0]&MPFR_LIMB_MASK(sh-1)))
 247             bcp1 = 1;
 248           else
 249             {
 250               /* We can't compute C'p+1 from C'. Compute it from C */
 251               /* Start from bit x=p-d+sh in mantissa C
 252                  (+sh since we have already looked sh bits in C'!) */
 253               mpfr_prec_t x = p-d+sh-1;
 254               if (MPFR_LIKELY(x>p))
 255                 /* We are already looked at all the bits of c, so C'p+1 = 0*/
 256                 bcp1 = 0;
 257               else
 258                 {
 259                   mp_limb_t *tp = MPFR_MANT(c);
 260                   mp_size_t kx = n-1 - (x / BITS_PER_MP_LIMB);
 261                   mpfr_prec_t sx = BITS_PER_MP_LIMB-1-(x%BITS_PER_MP_LIMB);
 262                   DEBUG (printf ("(First) x=%lu Kx=%ld Sx=%lu\n",
 263                                  (unsigned long) x, (long) kx,
 264                                  (unsigned long) sx));
 265                   /* Looks at the last bits of limb kx (if sx=0 does nothing)*/
 266                   if (tp[kx] & MPFR_LIMB_MASK(sx))
 267                     bcp1 = 1;
 268                   else
 269                     {
 270                       /*kx += (sx==0);*/
 271                       /*If sx==0, tp[kx] hasn't been checked*/
 272                       do {
 273                         kx--;
 274                       } while (kx>=0 && tp[kx]==0);
 275                       bcp1 = (kx >= 0);
 276                     }
 277                 }
 278             }
 279         }
 280       else /* sh == 0 */
 281         {
 282           /* Compute Cp and C'p+1 from C with sh=0 */
 283           mp_limb_t *tp = MPFR_MANT(c);
 284           /* Start from bit x=p-d in mantissa C */
 285           mpfr_prec_t  x = p-d;
 286           mp_size_t   kx = n-1 - (x / BITS_PER_MP_LIMB);
 287           mpfr_prec_t sx = BITS_PER_MP_LIMB-1-(x%BITS_PER_MP_LIMB);
 288           MPFR_ASSERTD(p >= d);
 289           bcp = tp[kx] & (MPFR_LIMB_ONE<<sx);
 290           /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
 291           if (tp[kx]&MPFR_LIMB_MASK(sx))
 292             bcp1 = 1;
 293           else
 294             {
 295               do {
 296                 kx--;
 297               } while (kx>=0 && tp[kx]==0);
 298               bcp1 = (kx>=0);
 299             }
 300         }
 301       DEBUG (printf("sh=%u Cp=%lu C'p+1=%lu\n", sh,
 302                     (unsigned long) bcp, (unsigned long) bcp1));
 303
 304       /* Clean shifted C' */
 305       mask = ~MPFR_LIMB_MASK(sh);
 306       cp[0] &= mask;
 307
 308       /* Add the mantissa c from b in a */
 309       ap = MPFR_MANT(a);
 310       limb = mpn_add_n (ap, MPFR_MANT(b), cp, n);
 311       DEBUG( mpfr_print_mant_binary("Add=  ", ap, p) );
 312
 313       /* Check for overflow */
 314       if (MPFR_UNLIKELY (limb))
 315         {
 316           limb = ap[0] & (MPFR_LIMB_ONE<<sh); /* Get LSB */
 317           mpn_rshift (ap, ap, n, 1);          /* Shift mantissa*/
 318           bx++;                               /* Fix exponent */
 319           ap[n-1] |= MPFR_LIMB_HIGHBIT;       /* Set MSB */
 320           ap[0]   &= mask;                    /* Clear LSB bit */
 321           bcp1    |= bcp;                     /* Recompute C'p+1 */
 322           bcp      = limb;                    /* Recompute Cp */
 323           DEBUG (printf ("(Overflow) Cp=%lu C'p+1=%lu\n",
 324                          (unsigned long) bcp, (unsigned long) bcp1));
 325           DEBUG (mpfr_print_mant_binary ("Add=  ", ap, p));
 326         }
 327
 328       /* Round:
 329           Zero: Truncate but could be exact.
 330           Away: Add 1 if Cp or C'p+1 !=0
 331           Nearest: Truncate but could be exact if Cp==0
 332                    Add 1 if C'p+1 !=0,
 333                    Even rule else */
 334       if (MPFR_LIKELY(rnd_mode == GMP_RNDN))
 335         {
 336           if (MPFR_LIKELY(bcp == 0))
 337             { inexact = MPFR_LIKELY(bcp1) ? -1 : 0; goto set_exponent; }
 338           else if (MPFR_UNLIKELY(bcp1==0) && (ap[0]&(MPFR_LIMB_ONE<<sh))==0)
 339             { inexact = -1; goto set_exponent; }
 340           else
 341             goto add_one_ulp;
 342         }
 343       MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b));
 344       if (rnd_mode == GMP_RNDZ)
 345         {
 346           inexact = MPFR_LIKELY(bcp || bcp1) ? -1 : 0;
 347           goto set_exponent;
 348         }
 349       else
 350         {
 351           if (MPFR_UNLIKELY(bcp==0 && bcp1==0))
 352             { inexact = 0; goto set_exponent; }
 353           else
 354             goto add_one_ulp;
 355         }
 356     }
 357   MPFR_ASSERTN(0);
 358
 359  add_one_ulp:
 360   /* add one unit in last place to a */
 361   DEBUG( printf("AddOneUlp\n") );
 362   if (MPFR_UNLIKELY( mpn_add_1(ap, ap, n, MPFR_LIMB_ONE<<sh) ))
 363     {
 364       /* Case 100000x0 = 0x1111x1 + 1*/
 365       DEBUG( printf("Pow of 2\n") );
 366       bx++;
 367       ap[n-1] = MPFR_LIMB_HIGHBIT;
 368     }
 369   inexact = 1;
 370
 371  set_exponent:
 372   if (MPFR_UNLIKELY(bx > __gmpfr_emax)) /* Check for overflow */
 373     {
 374       DEBUG( printf("Overflow\n") );
 375       MPFR_TMP_FREE(marker);
 376       MPFR_SET_SAME_SIGN(a,b);
 377       return mpfr_overflow(a, rnd_mode, MPFR_SIGN(a));
 378     }
 379   MPFR_SET_EXP (a, bx);
 380   MPFR_SET_SAME_SIGN(a,b);
 381
 382   MPFR_TMP_FREE(marker);
 383   MPFR_RET (inexact * MPFR_INT_SIGN (a));
 384 }