1 /* mpfr_log -- natural logarithm of a floating-point number
3 Copyright 1999-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 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
);