1 /*
2 Copyright (C) 2015 Fredrik Johansson
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "flint/arith.h"
13 #include "flint/double_extras.h"
14 #include "arb.h"
15
16 /* \sum_{k=0}^{a-1} \frac{k^n}{k!} */
17 /* b = a * \frac{(a-1)^n}{(a-1)!} */
18 /* assumes n > 0, a >= 0 */
19 static void
lower_bound(mag_t bound,const fmpz_t a,const fmpz_t n)20 lower_bound(mag_t bound, const fmpz_t a, const fmpz_t n)
21 {
22 arb_t t, u;
23 slong wp;
24
25 if (fmpz_is_zero(a))
26 {
27 mag_zero(bound);
28 return;
29 }
30
31 wp = 10 + fmpz_bits(n);
32
33 arb_init(t);
34 arb_init(u);
35
36 /* decreasing condition: a * (a-1)^n < a^n */
37 arb_set_fmpz(t, a);
38 arb_pow_fmpz(t, t, n, wp);
39
40 arb_set_fmpz(u, a);
41 arb_sub_ui(u, u, 1, wp);
42 arb_pow_fmpz(u, u, n, wp);
43 arb_mul_fmpz(u, u, a, wp);
44
45 if (arb_lt(u, t))
46 {
47 arb_gamma_fmpz(t, a, wp);
48 arb_div(t, u, t, wp);
49 arb_get_mag(bound, t);
50 }
51 else
52 {
53 mag_inf(bound);
54 }
55
56 arb_clear(t);
57 arb_clear(u);
58 }
59
60 /*
61 b^n [ ((b+1)/b)^n ((b+2)/b)^n ]
62 --- [ 1 + ----------- + ----------- + .... ]
63 b! [ (b+1) (b+1)(b+2) ]
64 */
65 static void
upper_bound(mag_t bound,const fmpz_t b,const fmpz_t n)66 upper_bound(mag_t bound, const fmpz_t b, const fmpz_t n)
67 {
68 arb_t t, u;
69 slong wp;
70
71 wp = 10 + 2 * fmpz_bits(n);
72
73 arb_init(t);
74 arb_init(u);
75
76 /* decreasing condition: (1+1/b)^n / (b+1) < 1 */
77 /* geometric series factor: 1 + t^2 + t^3 + ... = 1/(1-t) */
78 arb_one(t);
79 arb_div_fmpz(t, t, b, wp);
80 arb_add_ui(t, t, 1, wp);
81 arb_pow_fmpz(t, t, n, wp);
82 arb_set_fmpz(u, b);
83 arb_add_ui(u, u, 1, wp);
84 arb_div(t, t, u, wp);
85
86 arb_one(u);
87 arb_sub(u, u, t, wp);
88
89 if (arb_is_positive(u))
90 {
91 arb_set_fmpz(t, b);
92 arb_pow_fmpz(t, t, n, wp);
93 arb_div(t, t, u, wp);
94
95 arb_set_fmpz(u, b);
96 arb_add_ui(u, u, 1, wp);
97 arb_gamma(u, u, wp);
98
99 arb_div(t, t, u, wp);
100 arb_get_mag(bound, t);
101 }
102 else
103 {
104 mag_inf(bound);
105 }
106
107 arb_clear(t);
108 arb_clear(u);
109 }
110
111 /* approximate; need not give a correct bound, but
112 should be accurate so that we find near-optimal cutoffs
113 (we compute correct bounds elsewhere) */
114 static void
_arb_bell_mag(fmpz_t mmag,const fmpz_t n,const fmpz_t k)115 _arb_bell_mag(fmpz_t mmag, const fmpz_t n, const fmpz_t k)
116 {
117 if (fmpz_cmp_ui(k, 1) <= 0)
118 {
119 fmpz_set(mmag, k);
120 }
121 else if (fmpz_bits(n) < 50)
122 {
123 double dn, dk, z, u, lg;
124
125 dn = fmpz_get_d(n);
126 dk = fmpz_get_d(k);
127
128 z = dk + 1.0;
129 u = 1.0 / z;
130 lg = 0.91893853320467274178 + (z-0.5)*log(z) - z;
131 lg = lg + u * (0.08333333333333333 - 0.00277777777777777778 * (u * u)
132 + 0.00079365079365079365079 * ((u * u) * (u * u)));
133 u = (dn * log(dk) - lg) * 1.4426950408889634074 + 1.0;
134 fmpz_set_d(mmag, u);
135 }
136 else
137 {
138 arb_t t, u;
139 arf_t bound;
140 slong prec;
141
142 arb_init(t);
143 arb_init(u);
144 arf_init(bound);
145
146 prec = 10 + 1.1 * fmpz_bits(n);
147
148 arb_log_fmpz(t, k, prec);
149 arb_mul_fmpz(t, t, n, prec);
150
151 arb_set_fmpz(u, k);
152 arb_add_ui(u, u, 1, prec);
153 arb_lgamma(u, u, prec);
154
155 arb_sub(t, t, u, prec);
156
157 arb_const_log2(u, prec);
158 arb_div(t, t, u, prec);
159
160 arf_set_mag(bound, arb_radref(t));
161 arf_add(bound, arb_midref(t), bound, prec, ARF_RND_CEIL);
162 arf_get_fmpz(mmag, bound, ARF_RND_CEIL);
163
164 arb_clear(t);
165 arb_clear(u);
166 arf_clear(bound);
167 }
168 }
169
170 void
arb_bell_find_cutoffs(fmpz_t A,fmpz_t B,fmpz_t M,fmpz_t Mmag,const fmpz_t n,slong prec)171 arb_bell_find_cutoffs(fmpz_t A, fmpz_t B, fmpz_t M, fmpz_t Mmag, const fmpz_t n, slong prec)
172 {
173 fmpz_t a, amag, b, bmag, m, mmag, w, wmag, Amag, Bmag;
174
175 fmpz_init(a); fmpz_init(amag);
176 fmpz_init(b); fmpz_init(bmag);
177 fmpz_init(m); fmpz_init(mmag);
178 fmpz_init(w); fmpz_init(wmag);
179 fmpz_init(Amag); fmpz_init(Bmag);
180
181 if (fmpz_bits(n) < 53 && 0)
182 {
183 double dn = fmpz_get_d(n);
184
185 fmpz_set_d(M, dn / d_lambertw(dn));
186 _arb_bell_mag(Mmag, n, M);
187 }
188 else
189 {
190 /* do ternary search for M */
191 fmpz_zero(a);
192 fmpz_mul_ui(b, n, 2);
193 fmpz_zero(amag);
194 fmpz_zero(bmag);
195
196 while (_fmpz_sub_small(b, a) > 4)
197 {
198 fmpz_sub(m, b, a);
199 fmpz_tdiv_q_ui(m, m, 3);
200 fmpz_mul_2exp(w, m, 1);
201 fmpz_add(m, a, m);
202 fmpz_add(w, a, w);
203
204 _arb_bell_mag(mmag, n, m);
205 _arb_bell_mag(wmag, n, w);
206
207 if (fmpz_cmp(mmag, wmag) < 0)
208 {
209 fmpz_swap(a, m);
210 fmpz_swap(amag, mmag);
211 }
212 else
213 {
214 fmpz_swap(b, w);
215 fmpz_swap(bmag, wmag);
216 }
217 }
218
219 fmpz_set(M, a);
220 fmpz_set(Mmag, amag);
221 }
222
223 /* bisect for A */
224 fmpz_zero(a);
225 fmpz_zero(amag);
226 fmpz_set(b, M);
227 fmpz_set(bmag, Mmag);
228
229 while (_fmpz_sub_small(b, a) > 4)
230 {
231 fmpz_sub(m, b, a);
232 fmpz_tdiv_q_2exp(m, m, 1);
233 fmpz_add(m, a, m);
234
235 _arb_bell_mag(mmag, n, m);
236
237 /* mmag < Mmag - p */
238 if (_fmpz_sub_small(mmag, Mmag) < -prec)
239 {
240 fmpz_swap(a, m);
241 fmpz_swap(amag, mmag);
242 }
243 else
244 {
245 fmpz_swap(b, m);
246 fmpz_swap(bmag, mmag);
247 }
248 }
249
250 fmpz_set(A, a);
251 fmpz_set(Amag, amag);
252
253 /* bisect for B */
254 fmpz_set(a, M);
255 fmpz_set(amag, Mmag);
256 fmpz_mul_ui(b, n, 2);
257 fmpz_zero(bmag);
258
259 while (_fmpz_sub_small(b, a) > 4)
260 {
261 fmpz_sub(m, b, a);
262 fmpz_tdiv_q_2exp(m, m, 1);
263 fmpz_add(m, a, m);
264
265 _arb_bell_mag(mmag, n, m);
266
267 /* mmag < Mmag - p */
268 if (_fmpz_sub_small(mmag, Mmag) < -prec || fmpz_sgn(mmag) <= 0)
269 {
270 fmpz_swap(b, m);
271 fmpz_swap(bmag, mmag);
272 }
273 else
274 {
275 fmpz_swap(a, m);
276 fmpz_swap(amag, mmag);
277 }
278 }
279
280 fmpz_set(B, a);
281 fmpz_set(Bmag, amag);
282
283 fmpz_clear(a); fmpz_clear(amag);
284 fmpz_clear(b); fmpz_clear(bmag);
285 fmpz_clear(m); fmpz_clear(mmag);
286 fmpz_clear(w); fmpz_clear(wmag);
287 fmpz_clear(Amag); fmpz_clear(Bmag);
288 }
289
290 void
arb_bell_fmpz(arb_t res,const fmpz_t n,slong prec)291 arb_bell_fmpz(arb_t res, const fmpz_t n, slong prec)
292 {
293 fmpz_t a, b, m, mmag, c;
294 arb_t t;
295 mag_t bound;
296
297 if (fmpz_sgn(n) < 0)
298 {
299 arb_zero(res);
300 return;
301 }
302
303 if (fmpz_fits_si(n))
304 {
305 slong nn = fmpz_get_si(n);
306
307 /* compute exactly if we would be computing at least half the bits
308 of the exact number */
309 if (nn < 50 || prec > 0.5 * nn * log(0.7*nn / log(nn)) * 1.442695041)
310 {
311 fmpz_init(a);
312 arith_bell_number(a, nn);
313 arb_set_round_fmpz(res, a, prec);
314 fmpz_clear(a);
315 return;
316 }
317 }
318
319 fmpz_init(a);
320 fmpz_init(b);
321 fmpz_init(m);
322 fmpz_init(mmag);
323 fmpz_init(c);
324 arb_init(t);
325 mag_init(bound);
326
327 arb_bell_find_cutoffs(a, b, m, mmag, n, 1.03 * prec + fmpz_bits(n) + 2);
328
329 /* cutoff: n > 2^12 * prec^2 */
330 fmpz_set_ui(c, prec);
331 fmpz_mul_ui(c, c, prec);
332 fmpz_mul_2exp(c, c, 12);
333
334 if (fmpz_cmp(n, c) > 0)
335 arb_bell_sum_taylor(res, n, a, b, mmag, prec + 2);
336 else
337 arb_bell_sum_bsplit(res, n, a, b, mmag, prec + 2);
338
339 lower_bound(bound, a, n);
340 arb_add_error_mag(res, bound);
341
342 upper_bound(bound, b, n);
343 arb_add_error_mag(res, bound);
344
345 arb_const_e(t, prec + 2);
346 arb_div(res, res, t, prec);
347
348 fmpz_clear(a);
349 fmpz_clear(b);
350 fmpz_clear(m);
351 fmpz_clear(mmag);
352 fmpz_clear(c);
353 arb_clear(t);
354 mag_clear(bound);
355 }
356
357 void
arb_bell_ui(arb_t res,ulong n,slong prec)358 arb_bell_ui(arb_t res, ulong n, slong prec)
359 {
360 fmpz_t t;
361 fmpz_init(t);
362 fmpz_set_ui(t, n);
363 arb_bell_fmpz(res, t, prec);
364 fmpz_clear(t);
365 }
366
367