1 /* mpfr_log2p1 -- Compute log2(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 log2(1+x) is exactly representable in infinite precision,
29 and in such case the returned value is k such that 1+x = 2^k (the case k=0
30 cannot happen since we assume x<>0) */
31 static mpfr_exp_t
32 mpfr_log2p1_isexact (mpfr_srcptr x)
33 {
34 /* log2(1+x) is exactly representable when 1+x is a power of two,
35 we thus simply compute 1+x with 1-bit 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;
40 mpfr_exp_t e;
41
42 mpfr_init2 (t, 1);
43 inex = mpfr_add_ui (t, x, 1, MPFR_RNDZ);
44 e = MPFR_GET_EXP (t);
45 mpfr_clear (t);
46 return inex == 0 ? e - 1 : 0;
47 }
48
49 /* in case x=2^k and we can decide of the correct rounding,
50 put the correctly-rounded value in y and return the corresponding
51 ternary value (which is necessarily non-zero),
52 otherwise return 0 */
53 static int
54 mpfr_log2p1_special (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
55 {
56 mpfr_exp_t expx = MPFR_GET_EXP(x);
57 mpfr_exp_t k = expx - 1, expk;
58 mpfr_prec_t prec;
59 mpfr_t t;
60 int inex;
61
62 if (k <= 0 || mpfr_cmp_si_2exp (x, 1, k) != 0)
63 return 0;
64 /* k < log2(1+x) < k + 1/x/log(2) < k + 2/x */
65 expk = MPFR_INT_CEIL_LOG2(k); /* exponent of k */
66 /* 2/x < 2^(2-EXP(x)) thus if 2-EXP(x) < expk - PREC(y) - 1,
67 we have 2/x < 1/4*ulp(k) and we can decide the correct rounding */
68 if (2 - expx >= expk - MPFR_PREC(y) - 1)
69 return 0;
70 prec = (MPFR_PREC(y) + 2 <= ULSIZE) ? ULSIZE : MPFR_PREC(y) + 2;
71 mpfr_init2 (t, prec);
72 mpfr_set_ui (t, k, MPFR_RNDZ); /* exact since prec >= ULSIZE */
73 mpfr_nextabove (t);
74 /* now k < t < k + 2/x and round(t) = round(log2(1+x)) */
75 inex = mpfr_set (y, t, rnd_mode);
76 mpfr_clear (t);
77 /* Warning: for RNDF, the mpfr_set calls above might return 0 */
78 return (rnd_mode == MPFR_RNDF) ? 1 : inex;
79 }
80
81 /* The computation of log2p1 is done by log2p1(x) = log1p(x)/log(2) */
82 int
83 mpfr_log2p1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
84 {
85 int comp, inexact, nloop;
86 mpfr_t t, lg2;
87 mpfr_prec_t Ny = MPFR_PREC(y), prec;
88 MPFR_ZIV_DECL (loop);
89 MPFR_SAVE_EXPO_DECL (expo);
90
91 MPFR_LOG_FUNC
92 (("x[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
93 ("y[%Pd]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y,
94 inexact));
95
96 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
97 return mpfr_log1p (y, x, rnd_mode); /* same result for singular cases */
98
99 comp = mpfr_cmp_si (x, -1);
100 /* log2p1(x) is undefined for x < -1 */
101 if (MPFR_UNLIKELY(comp <= 0))
102 {
103 if (comp == 0)
104 /* x=0: log2p1(-1)=-inf (divide-by-zero exception) */
105 {
106 MPFR_SET_INF (y);
107 MPFR_SET_NEG (y);
108 MPFR_SET_DIVBY0 ();
109 MPFR_RET (0);
110 }
111 MPFR_SET_NAN (y);
112 MPFR_RET_NAN;
113 }
114
115 MPFR_SAVE_EXPO_MARK (expo);
116
117 prec = Ny + MPFR_INT_CEIL_LOG2 (Ny) + 6;
118
119 mpfr_init2 (t, prec);
120 mpfr_init2 (lg2, prec);
121
122 MPFR_ZIV_INIT (loop, prec);
123 for (nloop = 0; ; nloop++)
124 {
125 mpfr_log1p (t, x, MPFR_RNDN);
126 mpfr_const_log2 (lg2, MPFR_RNDN);
127 mpfr_div (t, t, lg2, MPFR_RNDN);
128 /* t = log2(1+x) * (1 + theta)^3 where |theta| < 2^-prec,
129 for prec >= 2 we have |(1 + theta)^3 - 1| < 4*theta.
130 Note: contrary to log10p1, no underflow is possible in extended
131 exponent range, since for tiny x, |log2(1+x)| ~ |x|/log(2) >= |x|,
132 and x is representable, thus x/log(2) too. */
133 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, prec - 2, Ny, rnd_mode)))
134 break;
135
136 if (nloop == 0)
137 {
138 /* check for exact cases */
139 mpfr_exp_t k;
140
141 MPFR_LOG_MSG (("check for exact cases\n", 0));
142 k = mpfr_log2p1_isexact (x);
143 if (k != 0) /* 1+x = 2^k */
144 {
145 inexact = mpfr_set_si (y, k, rnd_mode);
146 goto end;
147 }
148
149 /* if x = 2^k with huge k, Ziv's loop will fail */
150 inexact = mpfr_log2p1_special (y, x, rnd_mode);
151 if (inexact != 0)
152 goto end;
153 }
154
155 MPFR_ZIV_NEXT (loop, prec);
156 mpfr_set_prec (t, prec);
157 mpfr_set_prec (lg2, prec);
158 }
159 inexact = mpfr_set (y, t, rnd_mode);
160
161 end:
162 MPFR_ZIV_FREE (loop);
163 mpfr_clear (t);
164 mpfr_clear (lg2);
165
166 MPFR_SAVE_EXPO_FREE (expo);
167 return mpfr_check_range (y, inexact, rnd_mode);
168 }