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