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