1 #include "cado.h" // IWYU pragma: keep
2 #include <cctype> // for isspace, isdigit
3 #include <cerrno> // for errno
4 #include <climits> // for ULONG_MAX
5 #include <cstdio> // for fprintf, stderr, fclose, fgets, fopen
6 #include <cstdlib> // for exit, strtoul, strtod, EXIT_FAILURE
7 #include <gmp.h> // for mpz_sizeinbase
8 #include "cxx_mpz.hpp" // for cxx_mpz
9 #include "fb-types.h" // for fbprime_t
10 #include "las-multiobj-globals.hpp" // for dlp_descent
11 #include "las-siever-config.hpp"
12 #include "las-todo-entry.hpp" // for las_todo_entry
13 #include "macros.h" // for ASSERT_ALWAYS, MAYBE_UNUSED
14 #include "params.h" // param_list_parse_*
15 #include "verbose.h" // verbose_output_print
16
17 /* siever_config stuff */
18
declare_usage(param_list_ptr pl)19 void siever_config::declare_usage(param_list_ptr pl)
20 {
21 param_list_decl_usage(pl, "I", "set sieving region to 2^I times J, with J <= 2^(I-1) ; -I x is equivalent to -I (2*x-1)");
22 param_list_decl_usage(pl, "A", "set sieving region to (at most) 2^A");
23 param_list_decl_usage(pl, "lim0", "factor base bound on side 0");
24 param_list_decl_usage(pl, "lim1", "factor base bound on side 1");
25 param_list_decl_usage(pl, "lpb0", "set large prime bound on side 0 to 2^lpb0");
26 param_list_decl_usage(pl, "lpb1", "set large prime bound on side 1 to 2^lpb1");
27 param_list_decl_usage(pl, "mfb0", "set rational cofactor bound on side 0 2^mfb0");
28 param_list_decl_usage(pl, "mfb1", "set rational cofactor bound on side 1 2^mfb1");
29 param_list_decl_usage(pl, "lambda0", "lambda value on side 0");
30 param_list_decl_usage(pl, "lambda1", "lambda value on side 1");
31 param_list_decl_usage(pl, "powlim0", "limit on powers on side 0");
32 param_list_decl_usage(pl, "powlim1", "limit on powers on side 1");
33 param_list_decl_usage(pl, "ncurves0", "controls number of curves on side 0");
34 param_list_decl_usage(pl, "ncurves1", "controls number of curves on side 1");
35 param_list_decl_usage(pl, "tdthresh", "trial-divide primes p/r <= ththresh (r=number of roots)");
36 param_list_decl_usage(pl, "skipped", "primes below this bound are not sieved at all");
37 param_list_decl_usage(pl, "bkthresh", "bucket-sieve primes p >= bkthresh (default 2^I)");
38 param_list_decl_usage(pl, "bkthresh1", "2-level bucket-sieve primes in [bkthresh1,lim] (default=lim, meaning inactive)");
39 param_list_decl_usage(pl, "bkmult", "multiplier to use for taking margin in the bucket allocation\n");
40 param_list_decl_usage(pl, "unsievethresh", "Unsieve all p > unsievethresh where p|gcd(a,b)");
41 }
42
display(int side,unsigned int bitsize) const43 void siever_config::display(int side, unsigned int bitsize) const /*{{{*/
44 {
45 if (bitsize == 0) return;
46
47 verbose_output_print(0, 2, "# Sieving parameters for q~2^%d on side %d\n",
48 bitsize, side);
49 /* Strive to keep these output lines untouched */
50 verbose_output_print(0, 2,
51 "# Sieving parameters: lim0=%lu lim1=%lu lpb0=%d lpb1=%d\n",
52 sides[0].lim, sides[1].lim,
53 sides[0].lpb, sides[1].lpb);
54 verbose_output_print(0, 2,
55 "# mfb0=%d mfb1=%d\n",
56 sides[0].mfb, sides[1].mfb);
57 if (sides[0].lambda != 0 || sides[1].lambda != 0) {
58 verbose_output_print(0, 2,
59 "# lambda0=%1.1f lambda1=%1.1f\n",
60 sides[0].lambda, sides[1].lambda);
61 }
62 }/*}}}*/
63
64 /* {{{ Parse default siever config (fill all possible fields). Return
65 * true if the parsed siever config is complete and can be used without
66 * per-special-q info. */
parse_default(siever_config & sc,param_list_ptr pl)67 bool siever_config::parse_default(siever_config & sc, param_list_ptr pl)
68 {
69 /* The default config is not necessarily a complete bit of
70 * information.
71 */
72
73 /*
74 * Note that lim0 lim1 powlim0 powlim1 are also parsed from
75 * fb_factorbase::fb_factorbase, using only the command line (and no
76 * hint file, in particular). This is intended to be the "max" pair
77 * of limits.
78 */
79
80 param_list_parse_double(pl, "lambda0", &(sc.sides[0].lambda));
81 param_list_parse_double(pl, "lambda1", &(sc.sides[1].lambda));
82 /*
83 * Note: the stuff about the config being complete or not is mostly
84 * rubbish now...
85 */
86 int complete = 1;
87 complete &= param_list_parse_ulong(pl, "lim0", &(sc.sides[0].lim));
88 complete &= param_list_parse_ulong(pl, "lim1", &(sc.sides[1].lim));
89 if (param_list_lookup_string(pl, "A")) {
90 complete &= param_list_parse_int (pl, "A", &(sc.logA));
91 if (param_list_lookup_string(pl, "I")) {
92 fprintf(stderr, "# -A and -I are incompatible\n");
93 exit(EXIT_FAILURE);
94 }
95 } else if (param_list_lookup_string(pl, "I")) {
96 int I;
97 complete &= param_list_parse_int (pl, "I", &I);
98 sc.logA = 2 * I - 1;
99 verbose_output_print(0, 1, "# Interpreting -I %d as meaning -A %d\n", I, sc.logA);
100 } else {
101 complete = 0;
102 }
103
104 if (sc.sides[0].lim > 4294967295UL || sc.sides[1].lim > 4294967295UL)
105 {
106 fprintf (stderr, "Error, lim0/lim1 must be < 2^32\n");
107 exit (EXIT_FAILURE);
108 }
109
110 complete &= param_list_parse_int(pl, "lpb0", &(sc.sides[0].lpb));
111 complete &= param_list_parse_int(pl, "mfb0", &(sc.sides[0].mfb));
112 complete &= param_list_parse_int(pl, "lpb1", &(sc.sides[1].lpb));
113 complete &= param_list_parse_int(pl, "mfb1", &(sc.sides[1].mfb));
114 if (!complete) {
115 verbose_output_print(0, 1, "# default siever configuration is incomplete ; required parameters are I, lim[01], lpb[01], mfb[01]\n");
116
117 }
118
119 // Sublattices?
120 sc.sublat_bound = 0; // no sublattices by default.
121 param_list_parse_uint(pl, "sublat", &sc.sublat_bound);
122
123 /* Parse optional siever configuration parameters */
124 param_list_parse_uint(pl, "tdthresh", &(sc.td_thresh));
125 param_list_parse_uint(pl, "skipped", &(sc.skipped));
126
127 if (param_list_parse_uint(pl, "unsievethresh", &(sc.unsieve_thresh))) {
128 verbose_output_print(0, 1, "# Un-sieving primes > %u\n",
129 sc.unsieve_thresh);
130 }
131
132 // XXX note that when the sieving range size varies with the
133 // special-q, we need to accept that the bucket threshold varies,
134 // too.
135 //
136 // As a consequence, we should make it possible to specify extra
137 // parameters in the hint file format (not just I,lim,lpb,mfb).
138 //
139 // The logic that sets bucket_thresh to 2^I by default is done late,
140 // namely when we are about to create a new slicing for the factor
141 // base.
142 /* overrides default only if parameter is given */
143 param_list_parse_ulong(pl, "bkthresh", &(sc.bucket_thresh));
144 param_list_parse_ulong(pl, "bkthresh1", &(sc.bucket_thresh1));
145
146 const char *powlim_params[2] = {"powlim0", "powlim1"};
147 for (int side = 0; side < 2; side++) {
148 if (!param_list_parse_ulong(pl, powlim_params[side],
149 &sc.sides[side].powlim)) {
150 if (sc.bucket_thresh) {
151 sc.sides[side].powlim = sc.bucket_thresh - 1;
152 } else {
153 /* include all powers. We'll discard all those that go to
154 * bucket sieving anyway */
155 sc.sides[side].powlim = ULONG_MAX;
156 }
157 /* This message is also printed by
158 * fb_factorbase::fb_factorbase
159 */
160 /*
161 verbose_output_print(0, 2,
162 "# Using default value of %lu for -%s\n",
163 sc.sides[side].powlim, powlim_params[side]);
164 */
165 }
166 }
167
168 const char *ncurves_params[2] = {"ncurves0", "ncurves1"};
169 for (int side = 0; side < 2; side++)
170 if (!param_list_parse_int(pl, ncurves_params[side],
171 &sc.sides[side].ncurves))
172 sc.sides[side].ncurves = -1;
173
174 return complete;
175 }
176 /* }}} */
177
178 /* returns a set of thresholds that is compatible with the command-line
179 * defaults that we see here, as well as the logI that we've just made
180 * our mind on using.
181 *
182 * XXX NOTE XXX : the scale field is set to 0 by this function.
183 *
184 * XXX NOTE XXX : the nr_workspaces field is set to 0 by this function,
185 * because the caller is expected to set that instead.
186 *
187 *
188 */
instantiate_thresholds(int side) const189 fb_factorbase::key_type siever_config::instantiate_thresholds(int side) const
190 {
191 fbprime_t fbb = sides[side].lim;
192 fbprime_t bucket_thresh = this->bucket_thresh;
193 fbprime_t bucket_thresh1 = this->bucket_thresh1;
194
195 if (bucket_thresh == 0)
196 bucket_thresh = 1UL << logI;
197 if (bucket_thresh < (1UL << logI)) {
198 verbose_output_print(0, 1, "# Warning: with logI = %d,"
199 " we can't have %lu as the bucket threshold. Using %lu\n",
200 logI,
201 (unsigned long) bucket_thresh,
202 1UL << logI);
203 bucket_thresh = 1UL << logI;
204 }
205
206 if (bucket_thresh > fbb) bucket_thresh = fbb;
207 if (bucket_thresh1 == 0 || bucket_thresh1 > fbb) bucket_thresh1 = fbb;
208 if (bucket_thresh > bucket_thresh1) bucket_thresh1 = bucket_thresh;
209
210 return fb_factorbase::key_type {
211 {{bucket_thresh, bucket_thresh1, fbb, fbb}},
212 td_thresh,
213 skipped,
214 0,
215 0
216 };
217 }
get_config_for_q(las_todo_entry const & doing) const218 siever_config siever_config_pool::get_config_for_q(las_todo_entry const & doing) const /*{{{*/
219 {
220 siever_config config = base;
221 unsigned int bitsize = mpz_sizeinbase(doing.p, 2);
222 int side = doing.side;
223
224 /* Do we have a hint table with specifically tuned parameters,
225 * well suited to this problem size ? */
226 siever_config const * adapted = get_hint(side, bitsize);
227 if (adapted) {
228 config = *adapted;
229 verbose_output_print(0, 1, "# Using parameters from hint list for q~2^%d on side %d [%d@%d]\n", bitsize, side, bitsize, side);
230 }
231
232 if (doing.iteration) {
233 verbose_output_print(0, 1, "#\n# NOTE:"
234 " we are re-playing this special-q because of"
235 " %d previous failed attempt(s)\n", doing.iteration);
236 /* update sieving parameters here */
237 double ratio = double(config.sides[0].mfb) /
238 double(config.sides[0].lpb);
239 config.sides[0].lpb += doing.iteration;
240 config.sides[0].mfb = ratio*config.sides[0].lpb;
241 ratio = double(config.sides[1].mfb) /
242 double(config.sides[1].lpb);
243 config.sides[1].lpb += doing.iteration;
244 config.sides[1].mfb = ratio*config.sides[1].lpb;
245 verbose_output_print(0, 1,
246 "# NOTE: current values of lpb/mfb: %d,%d %d,%d\n#\n",
247 config.sides[0].lpb,
248 config.sides[0].mfb,
249 config.sides[1].lpb,
250 config.sides[1].mfb);
251 }
252
253 return config;
254 }/*}}}*/
255
256 void
declare_usage(cxx_param_list & pl)257 siever_config_pool::declare_usage(cxx_param_list & pl)/*{{{*/
258 {
259 param_list_decl_usage(pl, "hint-table", "filename with per-special q sieving data");
260 if (dlp_descent)
261 param_list_decl_usage(pl, "descent-hint-table", "Alias to hint-table");
262 }
263 /*}}}*/
264
siever_config_pool(cxx_param_list & pl)265 siever_config_pool::siever_config_pool(cxx_param_list & pl)/*{{{*/
266 {
267 default_config_ptr = NULL;
268 if (siever_config::parse_default(base, pl))
269 default_config_ptr = &base;
270
271 /* support both, since we've got to realize it's not that much
272 * attached to sieving. */
273 const char * filename = param_list_lookup_string(pl, "hint-table");
274 if (dlp_descent) {
275 const char * filename2 = param_list_lookup_string(pl, "descent-hint-table");
276 if (!filename) {
277 filename = filename2;
278 }
279 if (!filename) {
280 if (!default_config_ptr) {
281 fprintf(stderr,
282 "Error: no default config set, and no hint table either\n");
283 exit(EXIT_FAILURE);
284 }
285 return;
286 }
287 }
288 if (!filename) return;
289
290 char line[1024];
291 FILE * f;
292 f = fopen(filename, "r");
293 if (f == NULL) {
294 fprintf(stderr, "%s: %s\n", filename, strerror(errno));
295 /* There's no point in proceeding, since it would really change
296 * the behaviour of the program to do so */
297 exit(1);
298 }
299 for(;;) {
300 char * x = fgets(line, sizeof(line), f);
301 double t;
302 unsigned long z;
303 /* Tolerate comments and blank lines */
304 if (x == NULL) break;
305 for( ; *x && isspace(*x) ; x++) ;
306 if (*x == '#') continue;
307 if (!*x) continue;
308
309 descent_hint h;
310 siever_config & sc(h);
311
312 /* start with the global defaults */
313 sc = base;
314
315 int side;
316 unsigned int bitsize;
317
318 z = strtoul(x, &x, 10);
319 ASSERT_ALWAYS(z > 0);
320 bitsize = z;
321 switch(*x++) {
322 case '@' :
323 side = strtoul(x, &x, 0);
324 ASSERT_ALWAYS(side < 2);
325 break;
326 default:
327 fprintf(stderr, "%s: parse error at %s\n", filename, line);
328 exit(1);
329 }
330 for( ; *x && isspace(*x) ; x++) ;
331 t = strtod(x, &x); ASSERT_ALWAYS(t >= 0);
332 h.expected_time = t;
333 for( ; *x && isspace(*x) ; x++) ;
334 t = strtod(x, &x); ASSERT_ALWAYS(t >= 0);
335 h.expected_success = t;
336 for( ; *x && isspace(*x) ; x++) ;
337 char letter MAYBE_UNUSED = *x;
338 for( ; *x && !isdigit(*x) ; x++) ;
339 z = strtoul(x, &x, 10); ASSERT_ALWAYS(z > 0);
340 if (letter == 'I') {
341 sc.logA = 2*z-1;
342 } else if (letter == 'A') {
343 sc.logA = z;
344 } else {
345 fprintf(stderr, "%s: parse error (want I= or A=) at %s\n", filename, line);
346 exit(EXIT_FAILURE);
347 }
348
349 for(int s = 0 ; s < 2 ; s++) {
350 for( ; *x && isspace(*x) ; x++) ;
351 z = strtoul(x, &x, 10); ASSERT_ALWAYS(z > 0);
352 sc.sides[s].lim = z;
353 for( ; *x && !isdigit(*x) ; x++) ;
354 z = strtoul(x, &x, 10); ASSERT_ALWAYS(z > 0);
355 sc.sides[s].lpb = z;
356 /* recognize this as a double. If it's < 10, we'll consider
357 * this means lambda */
358 {
359 for( ; *x && !isdigit(*x) ; x++) ;
360 double t = strtod(x, &x); ASSERT_ALWAYS(t > 0);
361 if (t < 10) {
362 sc.sides[s].lambda = t;
363 sc.sides[s].mfb = t * sc.sides[s].lpb;
364 /* Then no "lambda" is allowed */
365 continue;
366 } else {
367 sc.sides[s].mfb = t;
368 }
369 }
370 if (*x == ',') {
371 for( ; *x && !isdigit(*x) ; x++) ;
372 t = strtod(x, &x); ASSERT_ALWAYS(t > 0);
373 sc.sides[s].lambda = t;
374 } else {
375 /* this means "automatic" */
376 sc.sides[s].lambda = 0;
377 }
378 }
379 for( ; *x ; x++) ASSERT_ALWAYS(isspace(*x));
380
381 key_type K(side, bitsize);
382
383 if (hints.find(K) != hints.end()) {
384 fprintf(stderr, "Error: two hints found for %d@%d\n",
385 bitsize, side);
386 exit(EXIT_FAILURE);
387 }
388
389 hints[K] = h;
390 }
391 fclose(f);
392 }/*}}}*/
393
394