1 /* mpfr_const_log2 -- compute natural logarithm of 2
2
3 Copyright 1999, 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
24 #include "mpfr-impl.h"
25
26 /* Declare the cache */
27 #ifndef MPFR_USE_LOGGING
MPFR_DECL_INIT_CACHE(__gmpfr_cache_const_log2,mpfr_const_log2_internal)28 MPFR_DECL_INIT_CACHE (__gmpfr_cache_const_log2, mpfr_const_log2_internal)
29 #else
30 MPFR_DECL_INIT_CACHE (__gmpfr_normal_log2, mpfr_const_log2_internal)
31 MPFR_DECL_INIT_CACHE (__gmpfr_logging_log2, mpfr_const_log2_internal)
32 MPFR_THREAD_VAR (mpfr_cache_ptr, __gmpfr_cache_const_log2, __gmpfr_normal_log2)
33 #endif
34
35 /* Set User interface */
36 #undef mpfr_const_log2
37 int
38 mpfr_const_log2 (mpfr_ptr x, mpfr_rnd_t rnd_mode) {
39 return mpfr_cache (x, __gmpfr_cache_const_log2, rnd_mode);
40 }
41
42 /* Auxiliary function: Compute the terms from n1 to n2 (excluded)
43 3/4*sum((-1)^n*n!^2/2^n/(2*n+1)!, n = n1..n2-1).
44 Numerator is T[0], denominator is Q[0],
45 Compute P[0] only when need_P is non-zero.
46 Need 1+ceil(log(n2-n1)/log(2)) cells in T[],P[],Q[].
47 */
48 static void
S(mpz_t * T,mpz_t * P,mpz_t * Q,unsigned long n1,unsigned long n2,int need_P)49 S (mpz_t *T, mpz_t *P, mpz_t *Q, unsigned long n1, unsigned long n2, int need_P)
50 {
51 if (n2 == n1 + 1)
52 {
53 if (n1 == 0)
54 mpz_set_ui (P[0], 3);
55 else
56 {
57 mpz_set_ui (P[0], n1);
58 mpz_neg (P[0], P[0]);
59 }
60 /* since n1 <= N, where N is the value from mpfr_const_log2_internal(),
61 and N = w / 3 + 1, where w <= PREC_MAX <= ULONG_MAX, then
62 N <= floor(ULONG_MAX/3) + 1, thus 2*N+1 <= ULONG_MAX */
63 MPFR_STAT_STATIC_ASSERT (MPFR_PREC_MAX <= ULONG_MAX);
64 mpz_set_ui (Q[0], 2 * n1 + 1);
65 mpz_mul_2exp (Q[0], Q[0], 2);
66 mpz_set (T[0], P[0]);
67 }
68 else
69 {
70 unsigned long m = (n1 / 2) + (n2 / 2) + (n1 & 1UL & n2);
71 mp_bitcnt_t v, w;
72
73 S (T, P, Q, n1, m, 1);
74 S (T + 1, P + 1, Q + 1, m, n2, need_P);
75 mpz_mul (T[0], T[0], Q[1]);
76 mpz_mul (T[1], T[1], P[0]);
77 mpz_add (T[0], T[0], T[1]);
78 if (need_P)
79 mpz_mul (P[0], P[0], P[1]);
80 mpz_mul (Q[0], Q[0], Q[1]);
81
82 /* remove common trailing zeroes if any */
83 v = mpz_scan1 (T[0], 0);
84 if (v > 0)
85 {
86 w = mpz_scan1 (Q[0], 0);
87 if (w < v)
88 v = w;
89 if (need_P)
90 {
91 w = mpz_scan1 (P[0], 0);
92 if (w < v)
93 v = w;
94 }
95 /* now v = min(val(T), val(Q), val(P)) */
96 if (v > 0)
97 {
98 mpz_fdiv_q_2exp (T[0], T[0], v);
99 mpz_fdiv_q_2exp (Q[0], Q[0], v);
100 if (need_P)
101 mpz_fdiv_q_2exp (P[0], P[0], v);
102 }
103 }
104 }
105 }
106
107 /* Don't need to save / restore exponent range: the cache does it */
108 int
mpfr_const_log2_internal(mpfr_ptr x,mpfr_rnd_t rnd_mode)109 mpfr_const_log2_internal (mpfr_ptr x, mpfr_rnd_t rnd_mode)
110 {
111 unsigned long n = MPFR_PREC (x);
112 mpfr_prec_t w; /* working precision */
113 unsigned long N;
114 mpz_t *T, *P, *Q;
115 mpfr_t t, q;
116 int inexact;
117 unsigned long lgN, i;
118 MPFR_GROUP_DECL(group);
119 MPFR_TMP_DECL(marker);
120 MPFR_ZIV_DECL(loop);
121
122 MPFR_LOG_FUNC (
123 ("rnd_mode=%d", rnd_mode),
124 ("x[%Pd]=%.*Rg inex=%d", mpfr_get_prec(x), mpfr_log_prec, x, inexact));
125
126 w = n + MPFR_INT_CEIL_LOG2 (n) + 3;
127
128 MPFR_TMP_MARK(marker);
129 MPFR_GROUP_INIT_2(group, w, t, q);
130
131 MPFR_ZIV_INIT (loop, w);
132 for (;;)
133 {
134 N = w / 3 + 1;
135
136 /* the following are needed for error analysis (see algorithms.tex) */
137 MPFR_ASSERTD(w >= 3 && N >= 2);
138
139 lgN = MPFR_INT_CEIL_LOG2 (N) + 1;
140 T = (mpz_t *) MPFR_TMP_ALLOC (3 * lgN * sizeof (mpz_t));
141 P = T + lgN;
142 Q = T + 2*lgN;
143 for (i = 0; i < lgN; i++)
144 {
145 mpz_init (T[i]);
146 mpz_init (P[i]);
147 mpz_init (Q[i]);
148 }
149
150 S (T, P, Q, 0, N, 0);
151
152 mpfr_set_z (t, T[0], MPFR_RNDN);
153 mpfr_set_z (q, Q[0], MPFR_RNDN);
154 mpfr_div (t, t, q, MPFR_RNDN);
155
156 for (i = 0; i < lgN; i++)
157 {
158 mpz_clear (T[i]);
159 mpz_clear (P[i]);
160 mpz_clear (Q[i]);
161 }
162
163 if (MPFR_CAN_ROUND (t, w - 2, n, rnd_mode))
164 break;
165
166 MPFR_ZIV_NEXT (loop, w);
167 MPFR_GROUP_REPREC_2(group, w, t, q);
168 }
169 MPFR_ZIV_FREE (loop);
170
171 inexact = mpfr_set (x, t, rnd_mode);
172
173 MPFR_GROUP_CLEAR(group);
174 MPFR_TMP_FREE(marker);
175
176 return inexact;
177 }
178