1 /* mpfr_exp2 -- power of 2 function 2^y 2 3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4 Contributed by the Arenaire and Caramel projects, INRIA. 5 6 This file is part of the GNU MPFR Library. 7 8 The GNU MPFR Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MPFR Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23 #define MPFR_NEED_LONGLONG_H 24 #include "mpfr-impl.h" 25 26 /* The computation of y = 2^z is done by * 27 * y = exp(z*log(2)). The result is exact iff z is an integer. */ 28 29 int 30 mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 31 { 32 int inexact; 33 long xint; 34 mpfr_t xfrac; 35 MPFR_SAVE_EXPO_DECL (expo); 36 37 MPFR_LOG_FUNC 38 (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, rnd_mode), 39 ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec(y), mpfr_log_prec, y, 40 inexact)); 41 42 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 43 { 44 if (MPFR_IS_NAN (x)) 45 { 46 MPFR_SET_NAN (y); 47 MPFR_RET_NAN; 48 } 49 else if (MPFR_IS_INF (x)) 50 { 51 if (MPFR_IS_POS (x)) 52 MPFR_SET_INF (y); 53 else 54 MPFR_SET_ZERO (y); 55 MPFR_SET_POS (y); 56 MPFR_RET (0); 57 } 58 else /* 2^0 = 1 */ 59 { 60 MPFR_ASSERTD (MPFR_IS_ZERO(x)); 61 return mpfr_set_ui (y, 1, rnd_mode); 62 } 63 } 64 65 /* since the smallest representable non-zero float is 1/2*2^__gmpfr_emin, 66 if x < __gmpfr_emin - 1, the result is either 1/2*2^__gmpfr_emin or 0 */ 67 MPFR_ASSERTN (MPFR_EMIN_MIN >= LONG_MIN + 2); 68 if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emin - 1) < 0)) 69 { 70 mpfr_rnd_t rnd2 = rnd_mode; 71 /* in round to nearest mode, round to zero when x <= __gmpfr_emin-2 */ 72 if (rnd_mode == MPFR_RNDN && 73 mpfr_cmp_si_2exp (x, __gmpfr_emin - 2, 0) <= 0) 74 rnd2 = MPFR_RNDZ; 75 return mpfr_underflow (y, rnd2, 1); 76 } 77 78 MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX); 79 if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emax) >= 0)) 80 return mpfr_overflow (y, rnd_mode, 1); 81 82 /* We now know that emin - 1 <= x < emax. */ 83 84 MPFR_SAVE_EXPO_MARK (expo); 85 86 /* 2^x = 1 + x*log(2) + O(x^2) for x near zero, and for |x| <= 1 we have 87 |2^x - 1| <= x < 2^EXP(x). If x > 0 we must round away from 0 (dir=1); 88 if x < 0 we must round toward 0 (dir=0). */ 89 MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y, __gmpfr_one, - MPFR_GET_EXP (x), 0, 90 MPFR_SIGN(x) > 0, rnd_mode, expo, {}); 91 92 xint = mpfr_get_si (x, MPFR_RNDZ); 93 mpfr_init2 (xfrac, MPFR_PREC (x)); 94 mpfr_sub_si (xfrac, x, xint, MPFR_RNDN); /* exact */ 95 96 if (MPFR_IS_ZERO (xfrac)) 97 { 98 mpfr_set_ui (y, 1, MPFR_RNDN); 99 inexact = 0; 100 } 101 else 102 { 103 /* Declaration of the intermediary variable */ 104 mpfr_t t; 105 106 /* Declaration of the size variable */ 107 mpfr_prec_t Ny = MPFR_PREC(y); /* target precision */ 108 mpfr_prec_t Nt; /* working precision */ 109 mpfr_exp_t err; /* error */ 110 MPFR_ZIV_DECL (loop); 111 112 /* compute the precision of intermediary variable */ 113 /* the optimal number of bits : see algorithms.tex */ 114 Nt = Ny + 5 + MPFR_INT_CEIL_LOG2 (Ny); 115 116 /* initialise of intermediary variable */ 117 mpfr_init2 (t, Nt); 118 119 /* First computation */ 120 MPFR_ZIV_INIT (loop, Nt); 121 for (;;) 122 { 123 /* compute exp(x*ln(2))*/ 124 mpfr_const_log2 (t, MPFR_RNDU); /* ln(2) */ 125 mpfr_mul (t, xfrac, t, MPFR_RNDU); /* xfrac * ln(2) */ 126 err = Nt - (MPFR_GET_EXP (t) + 2); /* Estimate of the error */ 127 mpfr_exp (t, t, MPFR_RNDN); /* exp(xfrac * ln(2)) */ 128 129 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode))) 130 break; 131 132 /* Actualisation of the precision */ 133 MPFR_ZIV_NEXT (loop, Nt); 134 mpfr_set_prec (t, Nt); 135 } 136 MPFR_ZIV_FREE (loop); 137 138 inexact = mpfr_set (y, t, rnd_mode); 139 140 mpfr_clear (t); 141 } 142 143 mpfr_clear (xfrac); 144 mpfr_clear_flags (); 145 mpfr_mul_2si (y, y, xint, MPFR_RNDN); /* exact or overflow */ 146 /* Note: We can have an overflow only when t was rounded up to 2. */ 147 MPFR_ASSERTD (MPFR_IS_PURE_FP (y) || inexact > 0); 148 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); 149 MPFR_SAVE_EXPO_FREE (expo); 150 return mpfr_check_range (y, inexact, rnd_mode); 151 } 152