1 /* mpfr_log -- natural logarithm of a floating-point number
3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 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 log(x) is done using the formula :
27 if we want p bits of the result,
30 log(x) ~ ------------ - m log 2
33 where s = x 2^m > 2^(p/2)
35 More precisely, if F(x) = int(1/sqrt(1-(1-x^2)*sin(t)^2), t=0..PI/2),
36 then for s>=1.26 we have log(s) < F(4/s) < log(s)*(1+4/s^2)
37 from which we deduce pi/2/AG(1,4/s)*(1-4/s^2) < log(s) < pi/2/AG(1,4/s)
38 so the relative error 4/s^2 is < 4/2^p i.e. 4 ulps.
42 mpfr_log (mpfr_ptr r, mpfr_srcptr a, mpfr_rnd_t rnd_mode)
47 MPFR_SAVE_EXPO_DECL (expo);
49 MPFR_GROUP_DECL(group);
52 (("a[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (a), mpfr_log_prec, a, rnd_mode),
53 ("r[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (r), mpfr_log_prec, r,
57 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a)))
59 /* If a is NaN, the result is NaN */
65 /* check for infinity before zero */
66 else if (MPFR_IS_INF (a))
74 else /* log(+Inf) = +Inf */
83 MPFR_ASSERTD (MPFR_IS_ZERO (a));
87 MPFR_RET (0); /* log(0) is an exact -infinity */
90 /* If a is negative, the result is NaN */
91 else if (MPFR_UNLIKELY (MPFR_IS_NEG (a)))
96 /* If a is 1, the result is 0 */
97 else if (MPFR_UNLIKELY (MPFR_GET_EXP (a) == 1 && mpfr_cmp_ui (a, 1) == 0))
101 MPFR_RET (0); /* only "normal" case where the result is exact */
106 /* use initial precision about q+lg(q)+5 */
107 p = q + 5 + 2 * MPFR_INT_CEIL_LOG2 (q);
108 /* % ~(mpfr_prec_t)GMP_NUMB_BITS ;
109 m=q; while (m) { p++; m >>= 1; } */
110 /* if (MPFR_LIKELY(p % GMP_NUMB_BITS != 0))
111 p += GMP_NUMB_BITS - (p%GMP_NUMB_BITS); */
113 MPFR_SAVE_EXPO_MARK (expo);
114 MPFR_GROUP_INIT_2 (group, p, tmp1, tmp2);
116 MPFR_ZIV_INIT (loop, p);
122 /* Calculus of m (depends on p) */
123 m = (p + 1) / 2 - MPFR_GET_EXP (a) + 1;
125 mpfr_mul_2si (tmp2, a, m, MPFR_RNDN); /* s=a*2^m, err<=1 ulp */
126 mpfr_div (tmp1, __gmpfr_four, tmp2, MPFR_RNDN);/* 4/s, err<=2 ulps */
127 mpfr_agm (tmp2, __gmpfr_one, tmp1, MPFR_RNDN); /* AG(1,4/s),err<=3 ulps */
128 mpfr_mul_2ui (tmp2, tmp2, 1, MPFR_RNDN); /* 2*AG(1,4/s), err<=3 ulps */
129 mpfr_const_pi (tmp1, MPFR_RNDN); /* compute pi, err<=1ulp */
130 mpfr_div (tmp2, tmp1, tmp2, MPFR_RNDN); /* pi/2*AG(1,4/s), err<=5ulps */
131 mpfr_const_log2 (tmp1, MPFR_RNDN); /* compute log(2), err<=1ulp */
132 mpfr_mul_si (tmp1, tmp1, m, MPFR_RNDN); /* compute m*log(2),err<=2ulps */
133 mpfr_sub (tmp1, tmp2, tmp1, MPFR_RNDN); /* log(a), err<=7ulps+cancel */
135 if (MPFR_LIKELY (MPFR_IS_PURE_FP (tmp1) && MPFR_IS_PURE_FP (tmp2)))
137 cancel = MPFR_GET_EXP (tmp2) - MPFR_GET_EXP (tmp1);
138 MPFR_LOG_MSG (("canceled bits=%ld\n", (long) cancel));
140 if (MPFR_UNLIKELY (cancel < 0))
143 /* we have 7 ulps of error from the above roundings,
144 4 ulps from the 4/s^2 second order term,
145 plus the canceled bits */
146 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp1, p-cancel-4, q, rnd_mode)))
149 /* VL: I think it is better to have an increment that it isn't
150 too low; in particular, the increment must be positive even
151 if cancel = 0 (can this occur?). */
152 p += cancel >= 8 ? cancel : 8;
156 /* TODO: find why this case can occur and what is best to do
161 MPFR_ZIV_NEXT (loop, p);
162 MPFR_GROUP_REPREC_2 (group, p, tmp1, tmp2);
164 MPFR_ZIV_FREE (loop);
165 inexact = mpfr_set (r, tmp1, rnd_mode);
167 MPFR_GROUP_CLEAR (group);
169 MPFR_SAVE_EXPO_FREE (expo);
170 return mpfr_check_range (r, inexact, rnd_mode);