1/* mpfr_expm1 -- Compute exp(x)-1 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 expm1 is done by 27 expm1(x)=exp(x)-1 28 */ 29 30int 31mpfr_expm1 (mpfr_ptr y, mpfr_srcptr x , mpfr_rnd_t rnd_mode) 32{ 33 int inexact; 34 mpfr_exp_t ex; 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 /* check for inf or -inf (expm1(-inf)=-1) */ 50 else if (MPFR_IS_INF (x)) 51 { 52 if (MPFR_IS_POS (x)) 53 { 54 MPFR_SET_INF (y); 55 MPFR_SET_POS (y); 56 MPFR_RET (0); 57 } 58 else 59 return mpfr_set_si (y, -1, rnd_mode); 60 } 61 else 62 { 63 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 64 MPFR_SET_ZERO (y); /* expm1(+/- 0) = +/- 0 */ 65 MPFR_SET_SAME_SIGN (y, x); 66 MPFR_RET (0); 67 } 68 } 69 70 ex = MPFR_GET_EXP (x); 71 if (ex < 0) 72 { 73 /* For -1 < x < 0, abs(expm1(x)-x) < x^2/2. 74 For 0 < x < 1, abs(expm1(x)-x) < x^2. */ 75 if (MPFR_IS_POS (x)) 76 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 0, 1, rnd_mode, {}); 77 else 78 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 1, 0, rnd_mode, {}); 79 } 80 81 MPFR_SAVE_EXPO_MARK (expo); 82 83 if (MPFR_IS_NEG (x) && ex > 5) /* x <= -32 */ 84 { 85 mpfr_t minus_one, t; 86 mpfr_exp_t err; 87 88 mpfr_init2 (minus_one, 2); 89 mpfr_init2 (t, 64); 90 mpfr_set_si (minus_one, -1, MPFR_RNDN); 91 mpfr_const_log2 (t, MPFR_RNDU); /* round upward since x is negative */ 92 mpfr_div (t, x, t, MPFR_RNDU); /* > x / ln(2) */ 93 err = mpfr_cmp_si (t, MPFR_EMIN_MIN >= -LONG_MAX ? 94 MPFR_EMIN_MIN : -LONG_MAX) <= 0 ? 95 - (MPFR_EMIN_MIN >= -LONG_MAX ? MPFR_EMIN_MIN : -LONG_MAX) : 96 - mpfr_get_si (t, MPFR_RNDU); 97 /* exp(x) = 2^(x/ln(2)) 98 <= 2^max(MPFR_EMIN_MIN,-LONG_MAX,ceil(x/ln(2)+epsilon)) 99 with epsilon > 0 */ 100 mpfr_clear (t); 101 MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (y, minus_one, err, 0, 0, rnd_mode, 102 expo, { mpfr_clear (minus_one); }); 103 mpfr_clear (minus_one); 104 } 105 106 /* General case */ 107 { 108 /* Declaration of the intermediary variable */ 109 mpfr_t t; 110 /* Declaration of the size variable */ 111 mpfr_prec_t Ny = MPFR_PREC(y); /* target precision */ 112 mpfr_prec_t Nt; /* working precision */ 113 mpfr_exp_t err, exp_te; /* error */ 114 MPFR_ZIV_DECL (loop); 115 116 /* compute the precision of intermediary variable */ 117 /* the optimal number of bits : see algorithms.tex */ 118 Nt = Ny + MPFR_INT_CEIL_LOG2 (Ny) + 6; 119 120 /* if |x| is smaller than 2^(-e), we will loose about e bits in the 121 subtraction exp(x) - 1 */ 122 if (ex < 0) 123 Nt += - ex; 124 125 /* initialize auxiliary variable */ 126 mpfr_init2 (t, Nt); 127 128 /* First computation of expm1 */ 129 MPFR_ZIV_INIT (loop, Nt); 130 for (;;) 131 { 132 MPFR_BLOCK_DECL (flags); 133 134 /* exp(x) may overflow and underflow */ 135 MPFR_BLOCK (flags, mpfr_exp (t, x, MPFR_RNDN)); 136 if (MPFR_OVERFLOW (flags)) 137 { 138 inexact = mpfr_overflow (y, rnd_mode, MPFR_SIGN_POS); 139 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW); 140 break; 141 } 142 else if (MPFR_UNDERFLOW (flags)) 143 { 144 inexact = mpfr_set_si (y, -1, rnd_mode); 145 MPFR_ASSERTD (inexact == 0); 146 inexact = -1; 147 if (MPFR_IS_LIKE_RNDZ (rnd_mode, 1)) 148 { 149 inexact = 1; 150 mpfr_nexttozero (y); 151 } 152 break; 153 } 154 155 exp_te = MPFR_GET_EXP (t); /* FIXME: exp(x) may overflow! */ 156 mpfr_sub_ui (t, t, 1, MPFR_RNDN); /* exp(x)-1 */ 157 158 /* error estimate */ 159 /*err=Nt-(__gmpfr_ceil_log2(1+pow(2,MPFR_EXP(te)-MPFR_EXP(t))));*/ 160 err = Nt - (MAX (exp_te - MPFR_GET_EXP (t), 0) + 1); 161 162 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode))) 163 { 164 inexact = mpfr_set (y, t, rnd_mode); 165 break; 166 } 167 168 /* increase the precision */ 169 MPFR_ZIV_NEXT (loop, Nt); 170 mpfr_set_prec (t, Nt); 171 } 172 MPFR_ZIV_FREE (loop); 173 174 mpfr_clear (t); 175 } 176 177 MPFR_SAVE_EXPO_FREE (expo); 178 return mpfr_check_range (y, inexact, rnd_mode); 179} 180