1 /* mpfr_log10 -- logarithm in base 10.
3 Copyright 2001-2015 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 /* The computation of r=log10(a)
28 r=log10(a)=log(a)/log(10)
32 mpfr_log10 (mpfr_ptr r
, mpfr_srcptr a
, mpfr_rnd_t rnd_mode
)
35 MPFR_SAVE_EXPO_DECL (expo
);
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
));
42 /* If a is NaN, the result is NaN */
43 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a
)))
50 /* check for infinity before zero */
51 else if (MPFR_IS_INF (a
))
54 /* log10(-Inf) = NaN */
59 else /* log10(+Inf) = +Inf */
63 MPFR_RET (0); /* exact */
68 MPFR_ASSERTD (MPFR_IS_ZERO (a
));
72 MPFR_RET (0); /* log10(0) is an exact -infinity */
76 /* If a is negative, the result is NaN */
77 if (MPFR_UNLIKELY (MPFR_IS_NEG (a
)))
83 /* If a is 1, the result is 0 */
84 if (mpfr_cmp_ui (a
, 1) == 0)
88 MPFR_RET (0); /* result is exact */
91 MPFR_SAVE_EXPO_MARK (expo
);
95 /* Declaration of the intermediary variable */
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 */
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
);
107 /* initialise of intermediary variables */
111 /* First computation of log10 */
112 MPFR_ZIV_INIT (loop
, Nt
);
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) */
121 /* estimation of the error */
123 if (MPFR_LIKELY (MPFR_CAN_ROUND (t
, err
, Ny
, rnd_mode
)))
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? */
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)
135 /* actualisation of the precision */
136 MPFR_ZIV_NEXT (loop
, Nt
);
137 mpfr_set_prec (t
, Nt
);
138 mpfr_set_prec (tt
, Nt
);
140 MPFR_ZIV_FREE (loop
);
142 inexact
= mpfr_set (r
, t
, rnd_mode
);
148 MPFR_SAVE_EXPO_FREE (expo
);
149 return mpfr_check_range (r
, inexact
, rnd_mode
);