1 /* mpfr_const_euler -- Euler's constant
2
3 Copyright 2001-2023 Free Software Foundation, Inc.
4 Contributed by Fredrik Johansson.
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 /* The approximation error bound uses Theorem 1 and Remark 2 in
24 https://arxiv.org/pdf/1312.0039v1.pdf */
25
26 #define MPFR_NEED_LONGLONG_H
27 #include "mpfr-impl.h"
28
29 /* Declare the cache */
30 MPFR_DECL_INIT_CACHE (__gmpfr_cache_const_euler, mpfr_const_euler_internal)
31
32 /* Set User Interface */
33 #undef mpfr_const_euler
34 int
35 mpfr_const_euler (mpfr_ptr x, mpfr_rnd_t rnd_mode) {
36 return mpfr_cache (x, __gmpfr_cache_const_euler, rnd_mode);
37 }
38
39
40 typedef struct
41 {
42 mpz_t P;
43 mpz_t Q;
44 mpz_t T;
45 mpz_t C;
46 mpz_t D;
47 mpz_t V;
48 } mpfr_const_euler_bs_struct;
49
50 typedef mpfr_const_euler_bs_struct mpfr_const_euler_bs_t[1];
51
52 static void
53 mpfr_const_euler_bs_init (mpfr_const_euler_bs_t s)
54 {
55 mpz_init (s->P);
56 mpz_init (s->Q);
57 mpz_init (s->T);
58 mpz_init (s->C);
59 mpz_init (s->D);
60 mpz_init (s->V);
61 }
62
63 static void
64 mpfr_const_euler_bs_clear (mpfr_const_euler_bs_t s)
65 {
66 mpz_clear (s->P);
67 mpz_clear (s->Q);
68 mpz_clear (s->T);
69 mpz_clear (s->C);
70 mpz_clear (s->D);
71 mpz_clear (s->V);
72 }
73
74 static void
75 mpfr_const_euler_bs_1 (mpfr_const_euler_bs_t s,
76 unsigned long n1, unsigned long n2, unsigned long N,
77 int cont)
78 {
79 if (n2 - n1 == 1)
80 {
81 mpz_set_ui (s->P, N);
82 mpz_mul (s->P, s->P, s->P);
83 mpz_set_ui (s->Q, n1 + 1);
84 mpz_mul (s->Q, s->Q, s->Q);
85 mpz_set_ui (s->C, 1);
86 mpz_set_ui (s->D, n1 + 1);
87 mpz_set (s->T, s->P);
88 mpz_set (s->V, s->P);
89 }
90 else
91 {
92 mpfr_const_euler_bs_t L, R;
93 mpz_t t, u, v;
94 unsigned long m = (n1 + n2) / 2;
95
96 mpfr_const_euler_bs_init (L);
97 mpfr_const_euler_bs_init (R);
98 mpfr_const_euler_bs_1 (L, n1, m, N, 1);
99 mpfr_const_euler_bs_1 (R, m, n2, N, 1);
100
101 mpz_init (t);
102 mpz_init (u);
103 mpz_init (v);
104
105 if (cont)
106 mpz_mul (s->P, L->P, R->P);
107
108 mpz_mul (s->Q, L->Q, R->Q);
109 mpz_mul (s->D, L->D, R->D);
110
111 /* T = LP RT + RQ LT*/
112 mpz_mul (t, L->P, R->T);
113 mpz_mul (v, R->Q, L->T);
114 mpz_add (s->T, t, v);
115
116 /* C = LC RD + RC LD */
117 if (cont)
118 {
119 mpz_mul (s->C, L->C, R->D);
120 mpz_addmul (s->C, R->C, L->D);
121 }
122
123 /* V = RD (RQ LV + LC LP RT) + LD LP RV */
124 mpz_mul (u, L->P, R->V);
125 mpz_mul (u, u, L->D);
126 mpz_mul (v, R->Q, L->V);
127 mpz_addmul (v, t, L->C);
128 mpz_mul (v, v, R->D);
129 mpz_add (s->V, u, v);
130
131 mpfr_const_euler_bs_clear (L);
132 mpfr_const_euler_bs_clear (R);
133 mpz_clear (t);
134 mpz_clear (u);
135 mpz_clear (v);
136 }
137 }
138
139 static void
140 mpfr_const_euler_bs_2 (mpz_t P, mpz_t Q, mpz_t T,
141 unsigned long n1, unsigned long n2, unsigned long N,
142 int cont)
143 {
144 if (n2 - n1 == 1)
145 {
146 if (n1 == 0)
147 {
148 mpz_set_ui (P, 1);
149 mpz_set_ui (Q, 4 * N);
150 }
151 else
152 {
153 mpz_set_ui (P, 2 * n1 - 1);
154 mpz_pow_ui (P, P, 3);
155 mpz_set_ui (Q, 32 * n1);
156 mpz_mul_ui (Q, Q, N);
157 mpz_mul_ui (Q, Q, N);
158 }
159 mpz_set (T, P);
160 }
161 else
162 {
163 mpz_t P2, Q2, T2;
164 unsigned long m = (n1 + n2) / 2;
165
166 mpz_init (P2);
167 mpz_init (Q2);
168 mpz_init (T2);
169 mpfr_const_euler_bs_2 (P, Q, T, n1, m, N, 1);
170 mpfr_const_euler_bs_2 (P2, Q2, T2, m, n2, N, 1);
171 mpz_mul (T, T, Q2);
172 mpz_mul (T2, T2, P);
173 mpz_add (T, T, T2);
174 if (cont)
175 mpz_mul (P, P, P2);
176 mpz_mul (Q, Q, Q2);
177 mpz_clear (P2);
178 mpz_clear (Q2);
179 mpz_clear (T2);
180 }
181 }
182
183 int
184 mpfr_const_euler_internal (mpfr_ptr x, mpfr_rnd_t rnd)
185 {
186 mpfr_const_euler_bs_t sum;
187 mpz_t t, u, v;
188 unsigned long n, N;
189 mpfr_prec_t prec, wp, magn;
190 mpfr_t y;
191 int inexact;
192 MPFR_ZIV_DECL (loop);
193
194 prec = mpfr_get_prec (x);
195 wp = prec + MPFR_INT_CEIL_LOG2 (prec) + 5;
196
197 mpfr_init2 (y, wp);
198 mpfr_const_euler_bs_init (sum);
199 mpz_init (t);
200 mpz_init (u);
201 mpz_init (v);
202
203 MPFR_ZIV_INIT (loop, wp);
204 for (;;)
205 {
206 /* The approximation error is bounded by 24 exp(-8n) when
207 n > 1, which is smaller than 2^-wp if
208 n > (wp + log_2(24)) * (log(2)/8).
209 Note log2(24) < 5 and log(2)/8 < 866434 / 10000000. */
210 mpz_set_ui (t, wp + 5);
211 mpz_mul_ui (t, t, 866434);
212 mpz_cdiv_q_ui (t, t, 10000000);
213 n = mpz_get_ui (t);
214
215 /* It is sufficient to take N >= alpha*n + 1
216 where alpha = 3/LambertW(3/e) = 4.970625759544... */
217 mpz_set_ui (t, n);
218 mpz_mul_ui (t, t, 4970626);
219 mpz_cdiv_q_ui (t, t, 1000000);
220 mpz_add_ui (t, t, 1);
221 N = mpz_get_ui (t);
222
223 /* V / ((T + Q) * D) = S / I
224 where S = sum_{k=0}^{N-1} H_k n^(2k) / (k!)^2,
225 I = sum_{k=0}^{N-1} n^(2k) / (k!)^2 */
226 mpfr_const_euler_bs_1 (sum, 0, N, n, 0);
227 mpz_add (sum->T, sum->T, sum->Q);
228 mpz_mul (t, sum->T, sum->D);
229 mpz_mul_2exp (u, sum->V, wp);
230 mpz_tdiv_q (v, u, t);
231 /* v * 2^-wp = S/I with error < 1 */
232
233 /* C / (D * V) = U where
234 U = (1/(4n)) sum_{k=0}^{2n-1} [(2k)!]^3 / ((k!)^4 8^(2k) (2n)^(2k)) */
235 mpfr_const_euler_bs_2 (sum->C, sum->D, sum->V, 0, 2*n, n, 0);
236 mpz_mul (t, sum->Q, sum->Q);
237 mpz_mul (t, t, sum->V);
238 mpz_mul (u, sum->T, sum->T);
239 mpz_mul (u, u, sum->D);
240 mpz_mul_2exp (t, t, wp);
241 mpz_tdiv_q (t, t, u);
242 /* t * 2^-wp = U/I^2 with error < 1 */
243
244 /* gamma = S/I - U/I^2 - log(n) with error at most 2^-wp */
245 mpz_sub (v, v, t);
246 /* v * 2^-wp now equals gamma + log(n) with error at most 3*2^-wp */
247
248 /* log(n) < 2^ceil(log2(n)) */
249 magn = MPFR_INT_CEIL_LOG2(n);
250 mpfr_set_prec (y, wp + magn);
251 mpfr_set_ui (y, n, MPFR_RNDZ); /* exact */
252 mpfr_log (y, y, MPFR_RNDZ); /* error < 2^-wp */
253
254 mpfr_mul_2ui (y, y, wp, MPFR_RNDZ);
255 mpfr_z_sub (y, v, y, MPFR_RNDZ);
256 mpfr_div_2ui (y, y, wp, MPFR_RNDZ);
257 /* rounding error from the last subtraction < 2^-wp */
258 /* so y = gamma with error < 5*2^-wp */
259
260 if (MPFR_LIKELY (MPFR_CAN_ROUND (y, wp - 3, prec, rnd)))
261 break;
262
263 MPFR_ZIV_NEXT (loop, wp);
264 }
265
266 MPFR_ZIV_FREE (loop);
267 inexact = mpfr_set (x, y, rnd);
268
269 mpfr_clear (y);
270 mpz_clear (t);
271 mpz_clear (u);
272 mpz_clear (v);
273 mpfr_const_euler_bs_clear (sum);
274
275 return inexact; /* always inexact */
276 }