1 /* mpfr_set_z_2exp -- set a floating-point number from a multiple-precision
2 integer and an exponent
3
4 Copyright 1999-2023 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramba 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 https://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;
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 }