1/* mpfr_log1p -- Compute log(1+x) 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 log1p is done by 27 log1p(x)=log(1+x) */ 28 29int 30mpfr_log1p (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 31{ 32 int comp, inexact; 33 mpfr_exp_t ex; 34 MPFR_SAVE_EXPO_DECL (expo); 35 36 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 37 { 38 if (MPFR_IS_NAN (x)) 39 { 40 MPFR_SET_NAN (y); 41 MPFR_RET_NAN; 42 } 43 /* check for inf or -inf (result is not defined) */ 44 else if (MPFR_IS_INF (x)) 45 { 46 if (MPFR_IS_POS (x)) 47 { 48 MPFR_SET_INF (y); 49 MPFR_SET_POS (y); 50 MPFR_RET (0); 51 } 52 else 53 { 54 MPFR_SET_NAN (y); 55 MPFR_RET_NAN; 56 } 57 } 58 else /* x is zero */ 59 { 60 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 61 MPFR_SET_ZERO (y); /* log1p(+/- 0) = +/- 0 */ 62 MPFR_SET_SAME_SIGN (y, x); 63 MPFR_RET (0); 64 } 65 } 66 67 ex = MPFR_GET_EXP (x); 68 if (ex < 0) /* -0.5 < x < 0.5 */ 69 { 70 /* For x > 0, abs(log(1+x)-x) < x^2/2. 71 For x > -0.5, abs(log(1+x)-x) < x^2. */ 72 if (MPFR_IS_POS (x)) 73 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex - 1, 0, 0, rnd_mode, {}); 74 else 75 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, - ex, 0, 1, rnd_mode, {}); 76 } 77 78 comp = mpfr_cmp_si (x, -1); 79 /* log1p(x) is undefined for x < -1 */ 80 if (MPFR_UNLIKELY(comp <= 0)) 81 { 82 if (comp == 0) 83 /* x=0: log1p(-1)=-inf (division by zero) */ 84 { 85 MPFR_SET_INF (y); 86 MPFR_SET_NEG (y); 87 MPFR_RET (0); 88 } 89 MPFR_SET_NAN (y); 90 MPFR_RET_NAN; 91 } 92 93 MPFR_SAVE_EXPO_MARK (expo); 94 95 /* General case */ 96 { 97 /* Declaration of the intermediary variable */ 98 mpfr_t t; 99 /* Declaration of the size variable */ 100 mpfr_prec_t Ny = MPFR_PREC(y); /* target precision */ 101 mpfr_prec_t Nt; /* working precision */ 102 mpfr_exp_t err; /* error */ 103 MPFR_ZIV_DECL (loop); 104 105 /* compute the precision of intermediary variable */ 106 /* the optimal number of bits : see algorithms.tex */ 107 Nt = Ny + MPFR_INT_CEIL_LOG2 (Ny) + 6; 108 109 /* if |x| is smaller than 2^(-e), we will loose about e bits 110 in log(1+x) */ 111 if (MPFR_EXP(x) < 0) 112 Nt += -MPFR_EXP(x); 113 114 /* initialise of intermediary variable */ 115 mpfr_init2 (t, Nt); 116 117 /* First computation of log1p */ 118 MPFR_ZIV_INIT (loop, Nt); 119 for (;;) 120 { 121 /* compute log1p */ 122 inexact = mpfr_add_ui (t, x, 1, MPFR_RNDN); /* 1+x */ 123 /* if inexact = 0, then t = x+1, and the result is simply log(t) */ 124 if (inexact == 0) 125 { 126 inexact = mpfr_log (y, t, rnd_mode); 127 goto end; 128 } 129 mpfr_log (t, t, MPFR_RNDN); /* log(1+x) */ 130 131 /* the error is bounded by (1/2+2^(1-EXP(t))*ulp(t) (cf algorithms.tex) 132 if EXP(t)>=2, then error <= ulp(t) 133 if EXP(t)<=1, then error <= 2^(2-EXP(t))*ulp(t) */ 134 err = Nt - MAX (0, 2 - MPFR_GET_EXP (t)); 135 136 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode))) 137 break; 138 139 /* increase the precision */ 140 MPFR_ZIV_NEXT (loop, Nt); 141 mpfr_set_prec (t, Nt); 142 } 143 inexact = mpfr_set (y, t, rnd_mode); 144 145 end: 146 MPFR_ZIV_FREE (loop); 147 mpfr_clear (t); 148 } 149 150 MPFR_SAVE_EXPO_FREE (expo); 151 return mpfr_check_range (y, inexact, rnd_mode); 152} 153