1 /*
2 Copyright (C) 2011 Fredrik Johansson
3
4 This file is part of FLINT.
5
6 FLINT 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 <math.h>
13 #include "arith.h"
14
15
16 static slong
_bell_series_cutoff(slong n)17 _bell_series_cutoff(slong n)
18 {
19 double N, log_N, log_pow, log_fac;
20
21 N = n;
22 log_N = (N==0 ? 0 : log(N));
23 log_pow = n * log_N;
24 log_fac = N*log_N - N;
25 while (log_pow - log_fac >= -2)
26 {
27 N++;
28 log_N = log(N);
29 log_pow = n * log_N;
30 log_fac += log_N;
31 }
32 return N;
33 }
34
35 static void
_mpz_bell_bsplit(mpz_t P,mpz_t Q,slong a,slong b,slong n,slong bmax)36 _mpz_bell_bsplit(mpz_t P, mpz_t Q, slong a, slong b, slong n, slong bmax)
37 {
38 if (b - a < 20)
39 {
40 mpz_t u;
41 slong k;
42 mpz_init(u);
43 flint_mpz_set_ui(P, UWORD(0));
44 flint_mpz_set_ui(Q, UWORD(0));
45 flint_mpz_set_ui(Q, (b - 1 == bmax) ? UWORD(1) : b);
46 for (k = b - 1; k >= a; k--)
47 {
48 flint_mpz_set_ui(u, k);
49 flint_mpz_pow_ui(u, u, n);
50 mpz_addmul(P, Q, u);
51 if (k != a)
52 flint_mpz_mul_ui(Q, Q, k);
53 }
54 mpz_clear(u);
55 }
56 else
57 {
58 slong m;
59 mpz_t P1, Q2;
60 m = (a + b) / 2;
61 mpz_init(P1);
62 mpz_init(Q2);
63 _mpz_bell_bsplit(P1, Q, a, m, n, bmax);
64 _mpz_bell_bsplit(P, Q2, m, b, n, bmax);
65 mpz_mul(Q, Q, Q2);
66 mpz_addmul(P, P1, Q2);
67 mpz_clear(P1);
68 mpz_clear(Q2);
69 }
70 }
71
72 void
arith_bell_number_bsplit(fmpz_t b,ulong n)73 arith_bell_number_bsplit(fmpz_t b, ulong n)
74 {
75 slong N;
76 flint_bitcnt_t prec;
77 mpz_t P, Q;
78 mpfr_t Pf, Qf, E, one;
79
80 N = _bell_series_cutoff(n);
81
82 mpz_init(P);
83 mpz_init(Q);
84
85 _mpz_bell_bsplit(P, Q, 1, N + 1, n, N);
86
87 prec = mpz_sizeinbase(P, 2) - mpz_sizeinbase(Q, 2) + 10;
88
89 mpfr_init2(Pf, prec);
90 mpfr_init2(Qf, prec);
91 mpfr_init2(E, prec);
92 mpfr_init2(one, 2);
93
94 mpfr_set_z(Pf, P, GMP_RNDN);
95 mpfr_set_z(Qf, Q, GMP_RNDN);
96 mpfr_set_ui(one, 1, GMP_RNDN);
97 mpfr_exp(E, one, GMP_RNDN);
98 mpfr_mul(Qf, Qf, E, GMP_RNDN);
99 mpfr_div(Pf, Pf, Qf, GMP_RNDN);
100 mpfr_get_z(P, Pf, GMP_RNDN);
101
102 fmpz_set_mpz(b, P);
103
104 mpfr_clear(one);
105 mpfr_clear(Pf);
106 mpfr_clear(Qf);
107 mpfr_clear(E);
108 mpz_clear(P);
109 mpz_clear(Q);
110 }
111