1/* mpfr_exp2 -- power of 2 function 2^y 2 3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 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 29int 30mpfr_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