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 "arith.h"
13
14 void
arith_bell_number_multi_mod(fmpz_t res,ulong n)15 arith_bell_number_multi_mod(fmpz_t res, ulong n)
16 {
17 fmpz_comb_temp_t temp;
18 fmpz_comb_t comb;
19 nmod_t mod;
20 mp_ptr primes, residues;
21 slong k, num_primes;
22 flint_bitcnt_t size, prime_bits;
23
24 size = arith_bell_number_size(n);
25 prime_bits = FLINT_BITS - 1;
26 num_primes = (size + prime_bits - 1) / prime_bits;
27
28 primes = flint_malloc(num_primes * sizeof(mp_limb_t));
29 residues = flint_malloc(num_primes * sizeof(mp_limb_t));
30
31 primes[0] = n_nextprime(UWORD(1) << prime_bits, 0);
32 for (k = 1; k < num_primes; k++)
33 primes[k] = n_nextprime(primes[k-1], 0);
34
35 for (k = 0; k < num_primes; k++)
36 {
37 nmod_init(&mod, primes[k]);
38 residues[k] = arith_bell_number_nmod(n, mod);
39 }
40
41 fmpz_comb_init(comb, primes, num_primes);
42 fmpz_comb_temp_init(temp, comb);
43
44 fmpz_multi_CRT_ui(res, residues, comb, temp, 0);
45
46 fmpz_comb_clear(comb);
47 fmpz_comb_temp_clear(temp);
48
49 flint_free(primes);
50 flint_free(residues);
51 }
52