1 #include "cado.h" // IWYU pragma: keep
2
3 #include <inttypes.h> // for PRIu32
4 #include <stdlib.h> // for exit, EXIT_FAILURE
5 #include <algorithm> // for max
6 #include <cstdint> // for uint32_t
7 #include <memory> // for allocator_traits<>::value_type
8 #include <vector> // for vector
9
10 #include "las-cofactor.hpp"
11
12 #include "macros.h" // for ASSERT_ALWAYS, ASSERT
13
14 #include "ecm/facul.hpp" // for FACUL_SMOOTH, facul_both, FACUL_MAYBE, FAC...
15 #include "las-siever-config.hpp" // for siever_config::side_config, siever_...
16 #include "modredc_15ul.h" // for modredc15ul_clearmod, modredc15ul_initmod_int
17 #include "modredc_2ul2.h" // for modredc2ul2_clearmod, modredc2ul2_initmod_int
18 #include "modredc_ul.h" // for modredcul_clearmod, modredcul_initmod_ul
19 #include "params.h"
20
declare_usage(cxx_param_list & pl)21 void cofactorization_statistics::declare_usage(cxx_param_list & pl)
22 {
23 param_list_decl_usage(pl, "stats-cofact", "write statistics about the cofactorization step in file xxx");
24 }
25
26 // las_info::{init,clear,print}_cof_stats
cofactorization_statistics(param_list_ptr pl)27 cofactorization_statistics::cofactorization_statistics(param_list_ptr pl)
28 {
29 const char * statsfilename = param_list_lookup_string (pl, "stats-cofact");
30 if (!statsfilename) {
31 file = NULL;
32 return;
33 }
34 file = fopen (statsfilename, "w");
35 if (file == NULL) {
36 fprintf (stderr, "Error, cannot create file %s\n", statsfilename);
37 exit (EXIT_FAILURE);
38 }
39 }
40
call(int bits0,int bits1)41 void cofactorization_statistics::call(int bits0, int bits1)
42 {
43 if (!file) return;
44 std::lock_guard<std::mutex> dummy(lock);
45 size_t s0 = cof_call.size();
46 if ((size_t) bits0 >= s0) {
47 size_t news0 = std::max((size_t) bits0+1, s0 + s0/2);
48 cof_call.insert(cof_call.end(), news0-s0, std::vector<uint32_t>());
49 cof_success.insert(cof_success.end(), news0-s0, std::vector<uint32_t>());
50 s0 = news0;
51 }
52 size_t s1 = cof_call[bits0].size();
53 if ((size_t) bits1 >= s1) {
54 size_t news1 = std::max((size_t) bits1+1, s1 + s1/2);
55 cof_call[bits0].insert(cof_call[bits0].end(), news1-s1, 0);
56 cof_success[bits0].insert(cof_success[bits0].end(), news1-s1, 0);
57 s1 = news1;
58 }
59 /* no need to use a mutex here: either we use one thread only
60 to compute the cofactorization data and if several threads
61 the order is irrelevant. The only problem that can happen
62 is when two threads increase the value at the same time,
63 and it is increased by 1 instead of 2, but this should
64 happen rarely. */
65 cof_call[bits0][bits1]++;
66 }
67
print()68 void cofactorization_statistics::print()
69 {
70 if (!file) return;
71 for(size_t bits0 = 0 ; bits0 < cof_call.size() ; ++bits0) {
72 for(size_t bits1 = 0 ; bits1 < cof_call[bits0].size() ; ++bits1) {
73 fprintf (file, "%zu %zu %" PRIu32 " %" PRIu32 "\n",
74 bits0, bits1,
75 cof_call[bits0][bits1],
76 cof_success[bits0][bits1]);
77 }
78 }
79 }
80
~cofactorization_statistics()81 cofactorization_statistics::~cofactorization_statistics()
82 {
83 if (!file) return;
84 fclose (file);
85 }
86 //
87
88 /* {{{ factor_leftover_norm */
89
90 #define NMILLER_RABIN 1 /* in the worst case, what can happen is that a
91 composite number is declared as prime, thus
92 a relation might be missed, but this will not
93 affect correctness */
94 #define IS_PROBAB_PRIME(X) (0 != mpz_probab_prime_p((X), NMILLER_RABIN))
95 #define BITSIZE(X) (mpz_sizeinbase((X), 2))
96 #define NFACTORS 8 /* maximal number of large primes */
97
98 /************************ cofactorization ********************************/
99
100 /* {{{ cofactoring area */
101
102 /* Return 0 if the leftover norm n cannot yield a relation.
103
104 Possible cases, where qj represents a prime in [B,L], and rj a prime > L
105 (assuming L < B^2, which might be false for the DLP descent):
106 (0) n >= 2^mfb
107 (a) n < L: 1 or q1
108 (b) L < n < B^2: r1 -> cannot yield a relation
109 (c) B^2 < n < B*L: r1 or q1*q2
110 (d) B*L < n < L^2: r1 or q1*q2 or q1*r2
111 (e) L^2 < n < B^3: r1 or q1*r2 or r1*r2 -> cannot yield a relation
112 (f) B^3 < n < B^2*L: r1 or q1*r2 or r1*r2 or q1*q2*q3
113 (g) B^2*L < n < L^3: r1 or q1*r2 or r1*r2
114 (h) L^3 < n < B^4: r1 or q1*r2, r1*r2 or q1*q2*r3 or q1*r2*r3 or r1*r2*r3
115 -> cannot yield a relation
116 */
117 int
check_leftover_norm(cxx_mpz const & n,siever_config::side_config const & scs)118 check_leftover_norm (cxx_mpz const & n, siever_config::side_config const & scs)
119 {
120 size_t s = mpz_sizeinbase (n, 2);
121 unsigned int lpb = scs.lpb;
122 unsigned int mfb = scs.mfb;
123 unsigned int klpb;
124 double nd, kB, B;
125
126 if (s > mfb)
127 return 0; /* n has more than mfb bits, which is the given limit */
128
129 if (scs.lim == 0) {
130 /* special case when not sieving */
131 return 1;
132 }
133
134 if (s <= lpb)
135 return 1; /* case (a) */
136 /* Note that in the case where L > B^2, if we're below L it's still fine of
137 course, but we have no guarantee that our cofactor is prime... */
138
139 nd = mpz_get_d (n);
140 B = (double) scs.lim;
141 kB = B * B;
142 for (klpb = lpb; klpb < s; klpb += lpb, kB *= B)
143 {
144 /* invariant: klpb = k * lpb, kB = B^(k+1) */
145 if (nd < kB) /* L^k < n < B^(k+1) */
146 return 0;
147 }
148
149 /* Here we have L < n < 2^mfb. If n is composite and we wrongly consider
150 it prime, we'll return 0, thus we'll potentially miss a relation, but
151 we won't output a relation with a composite ideal, thus a base-2 strong
152 prime test is enough. */
153
154 // TODO: maybe we should pass the modulus to the facul machinery
155 // instead of reconstructing it.
156 int prime=0;
157 if (s <= MODREDCUL_MAXBITS) {
158 modulusredcul_t m;
159 ASSERT(mpz_fits_ulong_p(n));
160 modredcul_initmod_ul (m, mpz_get_ui(n));
161 prime = modredcul_sprp2(m);
162 modredcul_clearmod (m);
163 } else if (s <= MODREDC15UL_MAXBITS) {
164 modulusredc15ul_t m;
165 unsigned long t[2];
166 modintredc15ul_t nn;
167 size_t written;
168 mpz_export (t, &written, -1, sizeof(unsigned long), 0, 0, n);
169 ASSERT_ALWAYS(written <= 2);
170 modredc15ul_intset_uls (nn, t, written);
171 modredc15ul_initmod_int (m, nn);
172 prime = modredc15ul_sprp2(m);
173 modredc15ul_clearmod (m);
174 } else if (s <= MODREDC2UL2_MAXBITS) {
175 modulusredc2ul2_t m;
176 unsigned long t[2];
177 modintredc2ul2_t nn;
178 size_t written;
179 mpz_export (t, &written, -1, sizeof(unsigned long), 0, 0, n);
180 ASSERT_ALWAYS(written <= 2);
181 modredc2ul2_intset_uls (nn, t, written);
182 modredc2ul2_initmod_int (m, nn);
183 prime = modredc2ul2_sprp2(m);
184 modredc2ul2_clearmod (m);
185 } else {
186 prime = mpz_probab_prime_p (n, 1);
187 }
188 if (prime)
189 return 0; /* n is a pseudo-prime larger than L */
190 return 1;
191 }
192
193 /* This is the header-comment for the old factor_leftover_norm()
194 * function, that is now deleted */
195 /* This function was contributed by Jerome Milan (and bugs were introduced
196 by Paul Zimmermann :-).
197 Input: n - the number to be factored (leftover norm). Must be composite!
198 Assumed to have no factor < B (factor base bound).
199 L - large prime bound is L=2^l
200 Assumes n > 0.
201 Return value:
202 -1 if n has a prime factor larger than L
203 1 if all prime factors of n are < L
204 0 if n could not be completely factored
205 Output:
206 the prime factors of n are factors->data[0..factors->length-1],
207 with corresponding multiplicities multis[0..factors->length-1].
208 */
209
210 /* This is the same function as factor_leftover_norm() but it works
211 with both norms! It is used when we want to factor these norms
212 simultaneously and not one after the other.
213 Return values:
214 -1 one of the cofactors is not smooth
215 0 unable to fully factor one of the cofactors
216 1 both cofactors are smooth
217 */
218
factor_both_leftover_norms(std::array<cxx_mpz,2> & n,std::array<std::vector<cxx_mpz>,2> & factors,std::array<unsigned long,2> const & Bs,facul_strategies_t const * strat)219 int factor_both_leftover_norms(
220 std::array<cxx_mpz, 2> & n,
221 std::array<std::vector<cxx_mpz>, 2> & factors,
222 std::array<unsigned long, 2> const & Bs,
223 facul_strategies_t const * strat)
224 {
225 int is_smooth[2] = {FACUL_MAYBE, FACUL_MAYBE};
226 /* To remember if a cofactor is already factored.*/
227
228 for (int side = 0; side < 2; side++) {
229 factors[side].clear();
230
231 double B = (double) Bs[side];
232 /* If n < B^2, then n is prime, since all primes < B have been removed */
233 if (mpz_get_d (n[side]) < B * B)
234 is_smooth[side] = FACUL_SMOOTH;
235 }
236
237 /* call the facul library */
238 std::array<int, 2> facul_code = facul_both (factors, n, strat, is_smooth);
239
240 if (is_smooth[0] != FACUL_SMOOTH || is_smooth[1] != FACUL_SMOOTH) {
241 if (is_smooth[0] == FACUL_NOT_SMOOTH || is_smooth[1] == FACUL_NOT_SMOOTH)
242 return -1;
243 else
244 return 0;
245 }
246
247 /* now we know both cofactors are smooth */
248 for (int side = 0; side < 2; side++) {
249 /* facul_code[side] is the number of found (smooth) factors */
250 for (int i = 0; i < facul_code[side]; i++) {
251 /* we know that factors found by facul_both() are primes < L */
252 mpz_divexact (n[side], n[side], factors[side][i]);
253 /* repeated factors should not be a problem, since they will
254 be dealt correctly in the filtering */
255 }
256
257 /* since the cofactor is smooth, n[side] is a prime < L here */
258 if (mpz_cmp_ui (n[side], 1) > 0) {
259 /* 1 is special */
260 factors[side].push_back(n[side]);
261 }
262 }
263 return 1; /* both cofactors are smooth */
264 }
265
266
267 /*}}}*/
268 /*}}}*/
269
270