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