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