contrib/mpfr/src/set_z_exp.c

   1 /* mpfr_set_z_2exp -- set a floating-point number from a multiple-precision
   2                       integer and an exponent
   3
   4 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
   5 Contributed by the AriC 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 /* set f to the integer z multiplied by 2^e */
  28 int
  29 mpfr_set_z_2exp (mpfr_ptr f, mpz_srcptr z, mpfr_exp_t e, mpfr_rnd_t rnd_mode)
  30 {
  31   mp_size_t fn, zn, dif, en;
  32   int k, sign_z, inex;
  33   mp_limb_t *fp, *zp;
  34   mpfr_exp_t exp;
  35
  36   sign_z = mpz_sgn (z);
  37   if (MPFR_UNLIKELY (sign_z == 0)) /* ignore the exponent for 0 */
  38     {
  39       MPFR_SET_ZERO(f);
  40       MPFR_SET_POS(f);
  41       MPFR_RET(0);
  42     }
  43   MPFR_ASSERTD (sign_z == MPFR_SIGN_POS || sign_z == MPFR_SIGN_NEG);
  44
  45   zn = ABS(SIZ(z)); /* limb size of z */
  46   /* compute en = floor(e/GMP_NUMB_BITS) */
  47   en = (e >= 0) ? e / GMP_NUMB_BITS : (e + 1) / GMP_NUMB_BITS - 1;
  48   MPFR_ASSERTD (zn >= 1);
  49   if (MPFR_UNLIKELY (zn + en > MPFR_EMAX_MAX / GMP_NUMB_BITS + 1))
  50     return mpfr_overflow (f, rnd_mode, sign_z);
  51   /* because zn + en >= MPFR_EMAX_MAX / GMP_NUMB_BITS + 2
  52      implies (zn + en) * GMP_NUMB_BITS >= MPFR_EMAX_MAX + GMP_NUMB_BITS + 1
  53      and exp = zn * GMP_NUMB_BITS + e - k
  54              >= (zn + en) * GMP_NUMB_BITS - k > MPFR_EMAX_MAX */
  55
  56   fp = MPFR_MANT (f);
  57   fn = MPFR_LIMB_SIZE (f);
  58   dif = zn - fn;
  59   zp = PTR(z);
  60   count_leading_zeros (k, zp[zn-1]);
  61
  62   /* now zn + en <= MPFR_EMAX_MAX / GMP_NUMB_BITS + 1
  63      thus (zn + en) * GMP_NUMB_BITS <= MPFR_EMAX_MAX + GMP_NUMB_BITS
  64      and exp = zn * GMP_NUMB_BITS + e - k
  65              <= (zn + en) * GMP_NUMB_BITS - k + GMP_NUMB_BITS - 1
  66              <= MPFR_EMAX_MAX + 2 * GMP_NUMB_BITS - 1 */
  67   exp = (mpfr_prec_t) zn * GMP_NUMB_BITS + e - k;
  68   /* The exponent will be exp or exp + 1 (due to rounding) */
  69   if (MPFR_UNLIKELY (exp > __gmpfr_emax))
  70     return mpfr_overflow (f, rnd_mode, sign_z);
  71   if (MPFR_UNLIKELY (exp + 1 < __gmpfr_emin))
  72     return mpfr_underflow (f, rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode,
  73                            sign_z);
  74
  75   if (MPFR_LIKELY (dif >= 0))
  76     {
  77       mp_limb_t rb, sb, ulp;
  78       int sh;
  79
  80       /* number has to be truncated */
  81       if (MPFR_LIKELY (k != 0))
  82         {
  83           mpn_lshift (fp, &zp[dif], fn, k);
  84           if (MPFR_LIKELY (dif > 0))
  85             fp[0] |= zp[dif - 1] >> (GMP_NUMB_BITS - k);
  86         }
  87       else
  88         MPN_COPY (fp, zp + dif, fn);
  89
  90       /* Compute Rounding Bit and Sticky Bit */
  91       MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (f) );
  92       if (MPFR_LIKELY (sh != 0))
  93         {
  94           mp_limb_t mask = MPFR_LIMB_ONE << (sh-1);
  95           mp_limb_t limb = fp[0];
  96           rb = limb & mask;
  97           sb = limb & (mask-1);
  98           ulp = 2*mask;
  99           fp[0] = limb & ~(ulp-1);
 100         }
 101       else /* sh == 0 */
 102         {
 103           mp_limb_t mask = MPFR_LIMB_ONE << (GMP_NUMB_BITS - 1 - k);
 104           if (MPFR_LIKELY (dif > 0))
 105             {
 106               rb = zp[--dif] & mask;
 107               sb = zp[dif] & (mask-1);
 108             }
 109           else
 110             rb = sb = 0;
 111           k = 0;
 112           ulp = MPFR_LIMB_ONE;
 113         }
 114       if (MPFR_UNLIKELY (sb == 0) && MPFR_LIKELY (dif > 0))
 115         {
 116           sb = zp[--dif];
 117           if (MPFR_LIKELY (k != 0))
 118             sb &= MPFR_LIMB_MASK (GMP_NUMB_BITS - k);
 119           if (MPFR_UNLIKELY (sb == 0) && MPFR_LIKELY (dif > 0))
 120             do {
 121               sb = zp[--dif];
 122             } while (dif > 0 && sb == 0);
 123         }
 124
 125       /* Rounding */
 126       if (MPFR_LIKELY (rnd_mode == MPFR_RNDN))
 127         {
 128           if (rb == 0 || MPFR_UNLIKELY (sb == 0 && (fp[0] & ulp) == 0))
 129             goto trunc;
 130           else
 131             goto addoneulp;
 132         }
 133       else /* Not Nearest */
 134         {
 135           if (MPFR_LIKELY (MPFR_IS_LIKE_RNDZ (rnd_mode, sign_z < 0))
 136               || MPFR_UNLIKELY ( (sb | rb) == 0 ))
 137             goto trunc;
 138           else
 139             goto addoneulp;
 140         }
 141
 142     trunc:
 143       inex = MPFR_LIKELY ((sb | rb) != 0) ? -1 : 0;
 144       goto end;
 145
 146     addoneulp:
 147       inex = 1;
 148       if (MPFR_UNLIKELY (mpn_add_1 (fp, fp, fn, ulp)))
 149         {
 150           /* Pow 2 case */
 151           if (MPFR_UNLIKELY (exp == __gmpfr_emax))
 152             return mpfr_overflow (f, rnd_mode, sign_z);
 153           exp ++;
 154           fp[fn-1] = MPFR_LIMB_HIGHBIT;
 155         }
 156     end:
 157       (void) 0;
 158     }
 159   else   /* dif < 0: Mantissa F is strictly bigger than z's one */
 160     {
 161       if (MPFR_LIKELY (k != 0))
 162         mpn_lshift (fp - dif, zp, zn, k);
 163       else
 164         MPN_COPY (fp - dif, zp, zn);
 165       /* fill with zeroes */
 166       MPN_ZERO (fp, -dif);
 167       inex = 0; /* result is exact */
 168     }
 169
 170   if (MPFR_UNLIKELY (exp < __gmpfr_emin))
 171     {
 172       if (rnd_mode == MPFR_RNDN && inex == 0 && mpfr_powerof2_raw (f))
 173         rnd_mode = MPFR_RNDZ;
 174       return mpfr_underflow (f, rnd_mode, sign_z);
 175     }
 176
 177   MPFR_SET_EXP (f, exp);
 178   MPFR_SET_SIGN (f, sign_z);
 179   MPFR_RET (inex*sign_z);
 180 }