1/* mpfr_set_z_2exp -- set a floating-point number from a multiple-precision
2                      integer and an exponent
3
4Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
5Contributed by the AriC and Caramel projects, INRIA.
6
7This file is part of the GNU MPFR Library.
8
9The GNU MPFR Library is free software; you can redistribute it and/or modify
10it under the terms of the GNU Lesser General Public License as published by
11the Free Software Foundation; either version 3 of the License, or (at your
12option) any later version.
13
14The GNU MPFR Library is distributed in the hope that it will be useful, but
15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
17License for more details.
18
19You should have received a copy of the GNU Lesser General Public License
20along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
21http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2251 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 */
28int
29mpfr_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}
181