1 /* mpfr_log -- natural logarithm of a floating-point number
2
3 Copyright 1999-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
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.
12
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.
17
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 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 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 log(x) is done using the formula :
27 if we want p bits of the result,
28
29 pi
30 log(x) ~ ------------ - m log 2
31 2 AG(1,4/s)
32
33 where s = x 2^m > 2^(p/2)
34
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.
39 */
40
41 int
mpfr_log(mpfr_ptr r,mpfr_srcptr a,mpfr_rnd_t rnd_mode)42 mpfr_log (mpfr_ptr r, mpfr_srcptr a, mpfr_rnd_t rnd_mode)
43 {
44 int inexact;
45 mpfr_prec_t p, q;
46 mpfr_t tmp1, tmp2;
47 mpfr_exp_t exp_a;
48 MPFR_SAVE_EXPO_DECL (expo);
49 MPFR_ZIV_DECL (loop);
50 MPFR_GROUP_DECL(group);
51
52 MPFR_LOG_FUNC
53 (("a[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (a), mpfr_log_prec, a, rnd_mode),
54 ("r[%Pd]=%.*Rg inexact=%d", mpfr_get_prec (r), mpfr_log_prec, r,
55 inexact));
56
57 /* Special cases */
58 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a)))
59 {
60 /* If a is NaN, the result is NaN */
61 if (MPFR_IS_NAN (a))
62 {
63 MPFR_SET_NAN (r);
64 MPFR_RET_NAN;
65 }
66 /* check for infinity before zero */
67 else if (MPFR_IS_INF (a))
68 {
69 if (MPFR_IS_NEG (a))
70 /* log(-Inf) = NaN */
71 {
72 MPFR_SET_NAN (r);
73 MPFR_RET_NAN;
74 }
75 else /* log(+Inf) = +Inf */
76 {
77 MPFR_SET_INF (r);
78 MPFR_SET_POS (r);
79 MPFR_RET (0);
80 }
81 }
82 else /* a is zero */
83 {
84 MPFR_ASSERTD (MPFR_IS_ZERO (a));
85 MPFR_SET_INF (r);
86 MPFR_SET_NEG (r);
87 MPFR_SET_DIVBY0 ();
88 MPFR_RET (0); /* log(0) is an exact -infinity */
89 }
90 }
91
92 /* If a is negative, the result is NaN */
93 if (MPFR_UNLIKELY (MPFR_IS_NEG (a)))
94 {
95 MPFR_SET_NAN (r);
96 MPFR_RET_NAN;
97 }
98
99 exp_a = MPFR_GET_EXP (a);
100
101 /* If a is 1, the result is +0 */
102 if (MPFR_UNLIKELY (exp_a == 1 && mpfr_cmp_ui (a, 1) == 0))
103 {
104 MPFR_SET_ZERO (r);
105 MPFR_SET_POS (r);
106 MPFR_RET (0); /* only "normal" case where the result is exact */
107 }
108
109 q = MPFR_PREC (r);
110
111 /* use initial precision about q+2*lg(q)+cte */
112 p = q + 2 * MPFR_INT_CEIL_LOG2 (q) + 10;
113 /* % ~(mpfr_prec_t)GMP_NUMB_BITS ;
114 m=q; while (m) { p++; m >>= 1; } */
115 /* if (MPFR_LIKELY(p % GMP_NUMB_BITS != 0))
116 p += GMP_NUMB_BITS - (p%GMP_NUMB_BITS); */
117
118 MPFR_SAVE_EXPO_MARK (expo);
119 MPFR_GROUP_INIT_2 (group, p, tmp1, tmp2);
120
121 MPFR_ZIV_INIT (loop, p);
122 for (;;)
123 {
124 mpfr_t scaled_a;
125 mpfr_exp_t m;
126 mpfr_exp_t cancel;
127
128 /* Calculus of m (depends on p)
129 If mpfr_exp_t has N bits, then both (p + 3) / 2 and |exp_a| fit
130 on N-2 bits, so that there cannot be an overflow. */
131 m = (p + 3) / 2 - exp_a;
132
133 /* In standard configuration (_MPFR_EXP_FORMAT <= 3), one has
134 mpfr_exp_t <= long, so that the following assertion is always
135 true. This assertion is needed for the mpfr_mul_si below. */
136 MPFR_ASSERTN (m >= LONG_MIN && m <= LONG_MAX);
137
138 /* FIXME: Redo the error analysis. The error concerning the AGM
139 should be explained since 4/s is inexact (one needs a bound
140 on its derivative). */
141 MPFR_ALIAS (scaled_a, a, MPFR_SIGN_POS, (p + 3) / 2); /* s=a*2^m */
142 /* [FIXME] and one can have the equality, even if p is even.
143 This means that if a is a power of 2 and p is even, then
144 s = (1/2) * 2^((p+2)/2) = 2^(p/2), so that the condition
145 s > 2^(p/2) from algorithms.tex is not satisfied. */
146 mpfr_div (tmp1, __gmpfr_four, scaled_a, MPFR_RNDF); /* 4/s, err<=2 ulps */
147 mpfr_agm (tmp2, __gmpfr_one, tmp1, MPFR_RNDN); /* AG(1,4/s),err<=3 ulps */
148 mpfr_mul_2ui (tmp2, tmp2, 1, MPFR_RNDN); /* 2*AG(1,4/s), err<=3 ulps */
149 mpfr_const_pi (tmp1, MPFR_RNDN); /* compute pi, err<=1ulp */
150 mpfr_div (tmp2, tmp1, tmp2, MPFR_RNDN); /* pi/2*AG(1,4/s), err<=5ulps */
151 mpfr_const_log2 (tmp1, MPFR_RNDN); /* compute log(2), err<=1ulp */
152 mpfr_mul_si (tmp1, tmp1, m, MPFR_RNDN); /* compute m*log(2),err<=2ulps */
153 mpfr_sub (tmp1, tmp2, tmp1, MPFR_RNDN); /* log(a), err<=7ulps+cancel */
154
155 if (MPFR_LIKELY (MPFR_IS_PURE_FP (tmp1) && MPFR_IS_PURE_FP (tmp2)))
156 {
157 cancel = MPFR_GET_EXP (tmp2) - MPFR_GET_EXP (tmp1);
158 MPFR_LOG_MSG (("canceled bits=%ld\n", (long) cancel));
159 MPFR_LOG_VAR (tmp1);
160 if (MPFR_UNLIKELY (cancel < 0))
161 cancel = 0;
162
163 /* we have 7 ulps of error from the above roundings,
164 4 ulps from the 4/s^2 second order term,
165 plus the canceled bits */
166 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp1, p - cancel - 4, q, rnd_mode)))
167 break;
168
169 /* VL: I think it is better to have an increment that it isn't
170 too low; in particular, the increment must be positive even
171 if cancel = 0 (can this occur?). */
172 p += cancel + MPFR_INT_CEIL_LOG2 (p);
173 }
174 else
175 {
176 /* TODO: find why this case can occur and what is best to do
177 with it. */
178 p += MPFR_INT_CEIL_LOG2 (p);
179 }
180
181 MPFR_ZIV_NEXT (loop, p);
182 MPFR_GROUP_REPREC_2 (group, p, tmp1, tmp2);
183 }
184 MPFR_ZIV_FREE (loop);
185 inexact = mpfr_set (r, tmp1, rnd_mode);
186 /* We clean */
187 MPFR_GROUP_CLEAR (group);
188
189 MPFR_SAVE_EXPO_FREE (expo);
190 return mpfr_check_range (r, inexact, rnd_mode);
191 }
192