1 /* mpfr_log10p1 -- Compute log10(1+x)
2
3 Copyright 2001-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 /* needed for MPFR_INT_CEIL_LOG2 */
24 #include "mpfr-impl.h"
25
26 #define ULSIZE (sizeof (unsigned long) * CHAR_BIT)
27
28 /* Return non-zero if log10(1+x) is exactly representable in infinite
29 precision, and in such case the returned value is k such that 1+x = 10^k
30 (the case k=0 cannot happen since we assume x<>0). */
31 static mpfr_exp_t
32 mpfr_log10p1_exact_p (mpfr_srcptr x)
33 {
34 /* log10(1+x) is exactly representable when 1+x is a power of 10,
35 we thus simply compute 1+x with enough precision and check whether
36 the addition is exact. This routine is called with extended exponent
37 range, thus no need to extend it. */
38 mpfr_t t;
39 int inex, ret = 0;
40
41 MPFR_ASSERTD(!MPFR_IS_SINGULAR(x));
42 if (MPFR_IS_NEG(x) || MPFR_EXP(x) <= 3) /* x < 8 */
43 return 0;
44 mpfr_init2 (t, MPFR_PREC(x));
45 inex = mpfr_add_ui (t, x, 1, MPFR_RNDZ);
46 if (inex == 0) /* otherwise 1+x = 2^k, and cannot be a power of 10 */
47 {
48 mpfr_prec_t trailing_x = mpfr_min_prec (x);
49 mpfr_prec_t trailing_t = mpfr_min_prec (t);
50 if (trailing_x > trailing_t)
51 {
52 mpfr_prec_t k = trailing_x - trailing_t;
53 /* if 1+x = 10^k, then t has k more trailing zeros than x */
54 mpz_t z;
55 mpfr_t y;
56 mpz_init (z);
57 mpz_ui_pow_ui (z, 5, k);
58 mpfr_init2 (y, mpz_sizeinbase (z, 2));
59 mpfr_set_z_2exp (y, z, k, MPFR_RNDZ);
60 if (mpfr_equal_p (t, y))
61 ret = k;
62 mpfr_clear (y);
63 mpz_clear (z);
64 }
65 }
66 mpfr_clear (t);
67 return ret;
68 }
69
70 /* Deal with the case where x is small, so that log10(1+x) ~ x/log(10).
71 In case we can round correctly, put in y the correctly-rounded value,
72 and return the corresponding ternary value (which cannot be zero).
73 Otherwise return 0.
74 This routine cannot be called only once after the first failure of Ziv's
75 strategy, since it might be that it fails the first time, thus we need
76 to pass the (increasing) working precision 'prec'.
77 In case of underflow, we set y to 0, and let the caller call
78 mpfr_underflow. */
79 static int
80 mpfr_log10p1_small (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode,
81 mpfr_prec_t prec)
82 {
83 mpfr_t t;
84 mpfr_exp_t e = MPFR_GET_EXP(x);
85 int inex;
86
87 /* for |x| < 1/2, |log10(x+1) - x/log(10)| < x^2/log(10) */
88 if (e > - (mpfr_exp_t) MPFR_PREC(y))
89 return 0; /* the term in x^2 will contribute */
90 /* now e = EXP(x) <= -PREC(y) <= -1 which ensures |x| < 1/2 */
91 mpfr_init2 (t, prec);
92 mpfr_log_ui (t, 10, MPFR_RNDN);
93 MPFR_SET_EXP (t, MPFR_GET_EXP (t) - 2);
94 /* we divide x by log(10)/4 which is smaller than 1 to avoid any underflow */
95 mpfr_div (t, x, t, MPFR_RNDN);
96 if (MPFR_GET_EXP (t) < __gmpfr_emin + 2) /* underflow case */
97 {
98 MPFR_SET_ZERO(y); /* the sign does not matter */
99 inex = 1;
100 }
101 else
102 {
103 MPFR_SET_EXP (t, MPFR_GET_EXP (t) - 2);
104 /* t = x/log(10) * (1 + theta)^2 where |theta| < 2^-prec.
105 For prec>=2, |(1 + theta)^2 - 1| < 3*theta thus the error is
106 bounded by 3 ulps. The error term in x^2 is bounded by |t*x|,
107 which is less than |t|*2^e < 2^(EXP(t)+e). */
108 e += prec;
109 /* now the error is bounded by 2^e+3 ulps */
110 e = (e >= 2) ? e + 1 : 3;
111 /* now the error is bounded by 2^e ulps */
112 if (MPFR_CAN_ROUND (t, prec - e, MPFR_PREC(y), rnd_mode))
113 inex = mpfr_set (y, t, rnd_mode);
114 else
115 inex = 0;
116 }
117 mpfr_clear (t);
118 return inex;
119 }
120
121 /* The computation of log10p1 is done by log10p1(x) = log1p(x)/log(2) */
122 int
123 mpfr_log10p1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
124 {
125 int comp, inexact, nloop;
126 mpfr_t t, lg10;
127 mpfr_prec_t Ny = MPFR_PREC(y), prec;
128 MPFR_ZIV_DECL (loop);
129 MPFR_SAVE_EXPO_DECL (expo);
130
131 MPFR_LOG_FUNC
132 (("x[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
133 ("y[%Pd]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y,
134 inexact));
135
136 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
137 return mpfr_log1p (y, x, rnd_mode); /* same result for singular cases */
138
139 comp = mpfr_cmp_si (x, -1);
140 /* log10p1(x) is undefined for x < -1 */
141 if (MPFR_UNLIKELY(comp <= 0))
142 {
143 if (comp == 0)
144 /* x=0: log10p1(-1)=-inf (divide-by-zero exception) */
145 {
146 MPFR_SET_INF (y);
147 MPFR_SET_NEG (y);
148 MPFR_SET_DIVBY0 ();
149 MPFR_RET (0);
150 }
151 MPFR_SET_NAN (y);
152 MPFR_RET_NAN;
153 }
154
155 MPFR_SAVE_EXPO_MARK (expo);
156
157 prec = Ny + MPFR_INT_CEIL_LOG2 (Ny) + 6;
158
159 mpfr_init2 (t, prec);
160 mpfr_init2 (lg10, prec);
161
162 MPFR_ZIV_INIT (loop, prec);
163 for (nloop = 0; ; nloop++)
164 {
165 mpfr_log1p (t, x, MPFR_RNDN);
166 mpfr_log_ui (lg10, 10, MPFR_RNDN);
167 mpfr_div (t, t, lg10, MPFR_RNDN);
168 /* t = log10(1+x) * (1 + theta)^3 where |theta| < 2^-prec,
169 for prec >= 2 we have |(1 + theta)^3 - 1| < 4*theta. */
170 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, prec - 2, Ny, rnd_mode)))
171 break;
172
173 if (nloop == 0)
174 {
175 /* check for exact cases */
176 mpfr_exp_t k;
177
178 MPFR_LOG_MSG (("check for exact cases\n", 0));
179 k = mpfr_log10p1_exact_p (x);
180 if (k != 0) /* 1+x = 10^k */
181 {
182 inexact = mpfr_set_si (y, k, rnd_mode);
183 goto end;
184 }
185 }
186
187 /* inexact will be the non-zero ternary value if rounding could be
188 done, otherwise it is set to 0. */
189 inexact = mpfr_log10p1_small (y, x, rnd_mode, prec);
190 if (inexact)
191 goto end;
192
193 MPFR_ZIV_NEXT (loop, prec);
194 mpfr_set_prec (t, prec);
195 mpfr_set_prec (lg10, prec);
196 }
197 inexact = mpfr_set (y, t, rnd_mode);
198
199 end:
200 MPFR_ZIV_FREE (loop);
201 mpfr_clear (t);
202 mpfr_clear (lg10);
203
204 MPFR_SAVE_EXPO_FREE (expo);
205 if (MPFR_IS_ZERO(y)) /* underflow from mpfr_log10p1_small */
206 return mpfr_underflow (y, (rnd_mode == MPFR_RNDN) ? MPFR_RNDZ : rnd_mode,
207 1);
208 else
209 return mpfr_check_range (y, inexact, rnd_mode);
210 }