1/* mpfr_set_z_2exp -- set a floating-point number from a multiple-precision
2                      integer and an exponent
3
4Copyright 1999-2023 Free Software Foundation, Inc.
5Contributed by the AriC and Caramba 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
21https://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;
32  int k, sign_z, inex;
33  mp_limb_t *fp, *zp;
34  mpfr_exp_t exp, nmax;
35  mpfr_uexp_t uexp;
36
37  sign_z = mpz_sgn (z);
38  if (MPFR_UNLIKELY (sign_z == 0)) /* ignore the exponent for 0 */
39    {
40      MPFR_SET_ZERO(f);
41      MPFR_SET_POS(f);
42      MPFR_RET(0);
43    }
44  MPFR_ASSERTD (sign_z == MPFR_SIGN_POS || sign_z == MPFR_SIGN_NEG);
45
46  zn = ABSIZ(z); /* limb size of z */
47  MPFR_ASSERTD (zn >= 1);
48  nmax = MPFR_EMAX_MAX / GMP_NUMB_BITS + 1;
49  /* Detect early overflow with zn + en > nmax,
50     where en = floor(e / GMP_NUMB_BITS).
51     This is checked without an integer overflow (even assuming some
52     future version of GMP, where limitations may be removed). */
53  if (MPFR_UNLIKELY (e >= 0 ?
54                     zn > nmax - e / GMP_NUMB_BITS :
55                     zn + (e + 1) / GMP_NUMB_BITS - 1 > nmax))
56    return mpfr_overflow (f, rnd_mode, sign_z);
57  /* because zn + en >= MPFR_EMAX_MAX / GMP_NUMB_BITS + 2
58     implies (zn + en) * GMP_NUMB_BITS >= MPFR_EMAX_MAX + GMP_NUMB_BITS + 1
59     and exp = zn * GMP_NUMB_BITS + e - k
60             >= (zn + en) * GMP_NUMB_BITS - k > MPFR_EMAX_MAX */
61
62  fp = MPFR_MANT (f);
63  fn = MPFR_LIMB_SIZE (f);
64  dif = zn - fn;
65  zp = PTR(z);
66  count_leading_zeros (k, zp[zn-1]);
67
68  /* now zn + en <= MPFR_EMAX_MAX / GMP_NUMB_BITS + 1
69     thus (zn + en) * GMP_NUMB_BITS <= MPFR_EMAX_MAX + GMP_NUMB_BITS
70     and exp = zn * GMP_NUMB_BITS + e - k
71             <= (zn + en) * GMP_NUMB_BITS - k + GMP_NUMB_BITS - 1
72             <= MPFR_EMAX_MAX + 2 * GMP_NUMB_BITS - 1 */
73  /* We need to compute exp = zn * GMP_NUMB_BITS + e - k with well-defined
74     operations (no integer overflows / no implementation-defined results).
75     The mathematical result of zn * GMP_NUMB_BITS may be larger than
76     the largest value of mpfr_exp_t while exp could still be less than
77     __gmpfr_emax. Thanks to early overflow detection, we can compute the
78     result in modular arithmetic, using mpfr_uexp_t, and convert it to
79     mpfr_exp_t. */
80  uexp = (mpfr_uexp_t) zn * GMP_NUMB_BITS + (mpfr_uexp_t) e - k;
81
82  /* Convert to signed in a portable way (see doc/README.dev).
83     On most platforms, this can be optimized to identity (no-op). */
84  exp = uexp > MPFR_EXP_MAX ? -1 - (mpfr_exp_t) ~uexp : (mpfr_exp_t) uexp;
85
86  /* The exponent will be exp or exp + 1 (due to rounding) */
87
88  if (MPFR_UNLIKELY (exp > __gmpfr_emax))
89    return mpfr_overflow (f, rnd_mode, sign_z);
90  if (MPFR_UNLIKELY (exp + 1 < __gmpfr_emin))
91    return mpfr_underflow (f, rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode,
92                           sign_z);
93
94  if (MPFR_LIKELY (dif >= 0))
95    {
96      mp_limb_t rb, sb, ulp;
97      int sh;
98
99      /* number has to be truncated */
100      if (MPFR_LIKELY (k != 0))
101        {
102          mpn_lshift (fp, &zp[dif], fn, k);
103          if (MPFR_UNLIKELY (dif > 0))
104            fp[0] |= zp[dif - 1] >> (GMP_NUMB_BITS - k);
105        }
106      else
107        MPN_COPY (fp, zp + dif, fn);
108
109      /* Compute Rounding Bit and Sticky Bit */
110      MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (f) );
111      if (MPFR_LIKELY (sh != 0))
112        {
113          mp_limb_t mask = MPFR_LIMB_ONE << (sh-1);
114          mp_limb_t limb = fp[0];
115          rb = limb & mask;
116          sb = limb & (mask-1);
117          ulp = 2*mask;
118          fp[0] = limb & ~(ulp-1);
119        }
120      else /* sh == 0 */
121        {
122          mp_limb_t mask = MPFR_LIMB_ONE << (GMP_NUMB_BITS - 1 - k);
123          if (MPFR_UNLIKELY (dif > 0))
124            {
125              rb = zp[--dif] & mask;
126              sb = zp[dif] & (mask-1);
127            }
128          else
129            rb = sb = 0;
130          k = 0;
131          ulp = MPFR_LIMB_ONE;
132        }
133      if (MPFR_UNLIKELY (sb == 0 && dif > 0))
134        {
135          sb = zp[--dif];
136          if (MPFR_LIKELY (k != 0))
137            sb &= MPFR_LIMB_MASK (GMP_NUMB_BITS - k);
138          if (MPFR_UNLIKELY (sb == 0) && MPFR_LIKELY (dif > 0))
139            do {
140              sb = zp[--dif];
141            } while (dif > 0 && sb == 0);
142        }
143
144      /* Rounding */
145      if (MPFR_LIKELY (rnd_mode == MPFR_RNDN))
146        {
147          if (rb == 0 || MPFR_UNLIKELY (sb == 0 && (fp[0] & ulp) == 0))
148            goto trunc;
149          else
150            goto addoneulp;
151        }
152      else /* Not Nearest */
153        {
154          if (MPFR_LIKELY (MPFR_IS_LIKE_RNDZ (rnd_mode, sign_z < 0))
155              || MPFR_UNLIKELY ( (sb | rb) == 0 ))
156            goto trunc;
157          else
158            goto addoneulp;
159        }
160
161    trunc:
162      inex = - ((sb | rb) != 0);
163      goto end;
164
165    addoneulp:
166      inex = 1;
167      if (MPFR_UNLIKELY (mpn_add_1 (fp, fp, fn, ulp)))
168        {
169          /* Pow 2 case */
170          if (MPFR_UNLIKELY (exp == __gmpfr_emax))
171            return mpfr_overflow (f, rnd_mode, sign_z);
172          exp ++;
173          fp[fn-1] = MPFR_LIMB_HIGHBIT;
174        }
175    end:
176      (void) 0;
177    }
178  else   /* dif < 0: Mantissa F is strictly bigger than z's one */
179    {
180      if (MPFR_LIKELY (k != 0))
181        mpn_lshift (fp - dif, zp, zn, k);
182      else
183        MPN_COPY (fp - dif, zp, zn);
184      /* fill with zeroes */
185      MPN_ZERO (fp, -dif);
186      inex = 0; /* result is exact */
187    }
188
189  if (MPFR_UNLIKELY (exp < __gmpfr_emin))
190    {
191      if (rnd_mode == MPFR_RNDN && inex == 0 && mpfr_powerof2_raw (f))
192        rnd_mode = MPFR_RNDZ;
193      return mpfr_underflow (f, rnd_mode, sign_z);
194    }
195
196  MPFR_SET_EXP (f, exp);
197  MPFR_SET_SIGN (f, sign_z);
198  MPFR_RET (inex*sign_z);
199}
200