1/* mpfr_log10 -- logarithm in base 10. 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 r=log10(a) 27 28 r=log10(a)=log(a)/log(10) 29 */ 30 31int 32mpfr_log10 (mpfr_ptr r, mpfr_srcptr a, mpfr_rnd_t rnd_mode) 33{ 34 int inexact; 35 MPFR_SAVE_EXPO_DECL (expo); 36 37 MPFR_LOG_FUNC 38 (("a[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (a), mpfr_log_prec, a, rnd_mode), 39 ("r[%Pu]=%.*Rg inexact=%d", 40 mpfr_get_prec (r), mpfr_log_prec, r, inexact)); 41 42 /* If a is NaN, the result is NaN */ 43 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a))) 44 { 45 if (MPFR_IS_NAN (a)) 46 { 47 MPFR_SET_NAN (r); 48 MPFR_RET_NAN; 49 } 50 /* check for infinity before zero */ 51 else if (MPFR_IS_INF (a)) 52 { 53 if (MPFR_IS_NEG (a)) 54 /* log10(-Inf) = NaN */ 55 { 56 MPFR_SET_NAN (r); 57 MPFR_RET_NAN; 58 } 59 else /* log10(+Inf) = +Inf */ 60 { 61 MPFR_SET_INF (r); 62 MPFR_SET_POS (r); 63 MPFR_RET (0); /* exact */ 64 } 65 } 66 else /* a = 0 */ 67 { 68 MPFR_ASSERTD (MPFR_IS_ZERO (a)); 69 MPFR_SET_INF (r); 70 MPFR_SET_NEG (r); 71 mpfr_set_divby0 (); 72 MPFR_RET (0); /* log10(0) is an exact -infinity */ 73 } 74 } 75 76 /* If a is negative, the result is NaN */ 77 if (MPFR_UNLIKELY (MPFR_IS_NEG (a))) 78 { 79 MPFR_SET_NAN (r); 80 MPFR_RET_NAN; 81 } 82 83 /* If a is 1, the result is 0 */ 84 if (mpfr_cmp_ui (a, 1) == 0) 85 { 86 MPFR_SET_ZERO (r); 87 MPFR_SET_POS (r); 88 MPFR_RET (0); /* result is exact */ 89 } 90 91 MPFR_SAVE_EXPO_MARK (expo); 92 93 /* General case */ 94 { 95 /* Declaration of the intermediary variable */ 96 mpfr_t t, tt; 97 MPFR_ZIV_DECL (loop); 98 /* Declaration of the size variable */ 99 mpfr_prec_t Ny = MPFR_PREC(r); /* Precision of output variable */ 100 mpfr_prec_t Nt; /* Precision of the intermediary variable */ 101 mpfr_exp_t err; /* Precision of error */ 102 103 /* compute the precision of intermediary variable */ 104 /* the optimal number of bits : see algorithms.tex */ 105 Nt = Ny + 4 + MPFR_INT_CEIL_LOG2 (Ny); 106 107 /* initialise of intermediary variables */ 108 mpfr_init2 (t, Nt); 109 mpfr_init2 (tt, Nt); 110 111 /* First computation of log10 */ 112 MPFR_ZIV_INIT (loop, Nt); 113 for (;;) 114 { 115 /* compute log10 */ 116 mpfr_set_ui (t, 10, MPFR_RNDN); /* 10 */ 117 mpfr_log (t, t, MPFR_RNDD); /* log(10) */ 118 mpfr_log (tt, a, MPFR_RNDN); /* log(a) */ 119 mpfr_div (t, tt, t, MPFR_RNDN); /* log(a)/log(10) */ 120 121 /* estimation of the error */ 122 err = Nt - 4; 123 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode))) 124 break; 125 126 /* log10(10^n) is exact: 127 FIXME: Can we have 10^n exactly representable as a mpfr_t 128 but n can't fit an unsigned long? */ 129 if (MPFR_IS_POS (t) 130 && mpfr_integer_p (t) && mpfr_fits_ulong_p (t, MPFR_RNDN) 131 && !mpfr_ui_pow_ui (tt, 10, mpfr_get_ui (t, MPFR_RNDN), MPFR_RNDN) 132 && mpfr_cmp (a, tt) == 0) 133 break; 134 135 /* actualisation of the precision */ 136 MPFR_ZIV_NEXT (loop, Nt); 137 mpfr_set_prec (t, Nt); 138 mpfr_set_prec (tt, Nt); 139 } 140 MPFR_ZIV_FREE (loop); 141 142 inexact = mpfr_set (r, t, rnd_mode); 143 144 mpfr_clear (t); 145 mpfr_clear (tt); 146 } 147 148 MPFR_SAVE_EXPO_FREE (expo); 149 return mpfr_check_range (r, inexact, rnd_mode); 150} 151