1 #include "cado.h" // IWYU pragma: keep
2 #include <cstddef> // for NULL, size_t
3 #include <algorithm> // for is_sorted, min
4 #include <array> // for array
5 #include <iterator> // for back_insert_iterator, back_ins...
6 #include <map> // for map, operator!=, _Rb_tree_iter...
7 #include <mutex> // for lock_guard, mutex
8 #include <utility> // for swap, pair
9 #include <vector> // for vector, vector<>::iterator
10 #include <gmp.h> // for mpz_divisible_ui_p, mpz_get_ui
11
12 #include "fb-types.h" // for fbprime_t
13 #include "fb.hpp" // for fb_factorbase::key_type, fb_en...
14 #include "cxx_mpz.hpp" // cxx_mpz
15 #include "getprime.h" // prime_info getprime_mt
16 #include "las-sieve-shared-data.hpp" // for sieve_shared_data::side_data
17 #include "lock_guarded_container.hpp" // for lock_guarded_container
18 #include "macros.h" // for MIN, ASSERT, ASSERT_ALWAYS
19 #include "mpz_poly.h"
20 #include "rootfinder.h" // mpz_poly_roots_ulong
21 #include "trialdiv.hpp" // for trialdiv_data, trialdiv_data::...
22
23
24 template<typename T>
25 static unsigned long
append_prime_list(T inserter,prime_info pi,unsigned long pmax,cxx_mpz_poly const & f,gmp_randstate_ptr rstate,int minroots=1)26 append_prime_list (T inserter, prime_info pi, unsigned long pmax, cxx_mpz_poly const & f, gmp_randstate_ptr rstate, int minroots = 1)
27 {
28 unsigned long p;
29 ASSERT_ALWAYS(minroots <= f->deg);
30 if (f->deg == 1) {
31 for (; (p = getprime_mt (pi)) < pmax; )
32 *inserter++ = p;
33 } else {
34 for (; (p = getprime_mt (pi)) < pmax; )
35 if (mpz_divisible_ui_p (f->coeff[f->deg], p) ||
36 mpz_poly_roots_ulong (NULL, f, p, rstate) >= minroots)
37 *inserter++ = p;
38 }
39 return p;
40 }
41
42
get_trialdiv_data(fb_factorbase::key_type fbK,fb_factorbase::slicing const * fbs)43 trialdiv_data const * sieve_shared_data::side_data::get_trialdiv_data(fb_factorbase::key_type fbK, fb_factorbase::slicing const * fbs)
44 {
45 std::lock_guard<std::mutex> foo(trialdiv_data_cache.mutex());
46 auto it = trialdiv_data_cache.find(fbK);
47 if (it != trialdiv_data_cache.end()) {
48 return &it->second;
49 }
50
51 ASSERT_ALWAYS(fbs != NULL);
52
53 /* Note that since we have trialdiv_data_cache.mutex() unlock, we may
54 * safely access the random state in this->rstate
55 */
56
57 /* Now compute the trialdiv data for these thresholds. */
58
59 /* Our trial division needs odd divisors, 2 is handled by mpz_even_p().
60 If the FB primes to trial divide contain 2, we skip over it.
61 We assume that if 2 is in the list, it is the first list entry,
62 and that it appears at most once. */
63
64 unsigned long pmax = std::min((unsigned long) fbK.thresholds[0],
65 trialdiv_data::max_p);
66
67 std::vector<unsigned long> trialdiv_primes = fbs->small_sieve_entries.skipped;
68
69 /* Maybe we can use the factor base. If we have one, of course ! */
70 unsigned long pmax_sofar = 0;
71
72 for(auto const & pp : fbs->small_sieve_entries.rest) {
73 if (pp.k > 1) continue;
74 trialdiv_primes.push_back(pp.p);
75 }
76 if (!trialdiv_primes.empty()) {
77 cxx_mpz zz(trialdiv_primes.back());
78 mpz_nextprime(zz, zz);
79 pmax_sofar = MIN(pmax, mpz_get_ui(zz));
80 }
81
82 if (pmax_sofar < pmax) {
83 /* we need some more. */
84 prime_info pi;
85 prime_info_init(pi);
86 unsigned long p;
87 /* first seek to the end of the fb. */
88 for ( ; (p = getprime_mt (pi)) < pmax_sofar ; );
89
90 for(int minroots = 1 ; minroots <= f->deg ; minroots++) {
91 p = append_prime_list(std::back_inserter(trialdiv_primes),
92 pi, MIN(pmax, minroots * fbK.td_thresh), f, rstate, minroots);
93 }
94 prime_info_clear (pi);
95 }
96
97 ASSERT(std::is_sorted(trialdiv_primes.begin(), trialdiv_primes.end()));
98 // std::sort(trialdiv_primes.begin(), trialdiv_primes.end());
99
100 /* note that we might have several "2"'s in the factor base because
101 * of powers: when bucket-sieving powers, we separate factor base
102 * entries that lead to distinct exponents */
103 size_t skip2 = 0;
104 for( ; skip2 < trialdiv_primes.size() && trialdiv_primes[skip2] == 2 ; skip2++);
105
106 trialdiv_data td(trialdiv_primes, skip2);
107 trialdiv_data_cache[fbK];
108 std::swap(trialdiv_data_cache[fbK], td);
109
110 return &trialdiv_data_cache[fbK];
111 }
112