1 /* dlpolyselect - discrete logarithm polynomial selection via Joux-Lercier's
2 algorithm
3
4 Example with Apache 512-bit key (Table 1 of
5 https://weakdh.org/imperfect-forward-secrecy-ccs15.pdf):
6
7 # we first search a degree (3,2) pair
8 $ ./dlpolyselect -df 3 -dg 2 -N 8372421755538377327377912526045445423027732035562313241965800453667849685158691589507936013805295187219621475007123900107532269487803598942841993804845107 -bound 8
9 ...
10 c3: 4
11 c2: 3
12 c1: -6
13 c0: -8
14 Y2: 430506467754036192750431535685941771886700493958222
15 Y1: 687078639517265884307374391507473735402330150842399
16 Y0: -533201855944366475689054404823391134825953676216635
17 skew: 1.23
18 # f lognorm 1.28, alpha 0.62, score 1.90
19 # g lognorm 116.59, alpha -2.18, score 114.41
20 # f+g score 116.31
21
22 # then we compare with a degree (4,3) pair
23 $ ./dlpolyselect -df 4 -dg 3 -N 8372421755538377327377912526045445423027732035562313241965800453667849685158691589507936013805295187219621475007123900107532269487803598942841993804845107 -bound 8
24 ...
25 c4: 6
26 c3: 1
27 c2: -8
28 c1: -3
29 c0: 6
30 Y3: -64071264306884991611859009153886700616
31 Y2: -115884379190374676852348454783130883186
32 Y1: 382823080720299801084267253734861739469
33 Y0: -14529492984288436819691699253895818591
34 skew: 1.05
35 # f lognorm 1.15, alpha -0.04, score 1.12
36 # g lognorm 87.27, alpha -4.07, score 83.20
37 # f+g score 84.32
38
39 To recover the polynomial pair used for the DLP768 record
40 (https://listserv.nodak.edu/cgi-bin/wa.exe?A2=NMBRTHRY;a0c66b63.1606 and
41 http://eprint.iacr.org/2017/067.pdf):
42
43 # -modr r -modm m enables one to only consider polynomials f
44 # with index r mod m in Stage 1 (default is r=0 and m=1).
45 # This is useful to distribute the work among several machines.
46 $ ./dlpolyselect -N 1219344858334286932696341909195796109526657386154251328029273656175766870980306505584577389125860826715201547225794072935883258868036433287217994721542199148182841505800433148410869683590659346847659519108393837414567892730579162319 -df 4 -dg 3 -bound 140 -modr 66234274 -modm 100000000
47 ...
48 c4: 55
49 c3: 5
50 c2: -86
51 c1: 34
52 c0: -140
53 Y3: 277260730400349522890422618473498148528706115003337935150
54 Y2: 217583293626947899787577441128333027617541095004734736415
55 Y1: -1937981312833038778565617469829395544065255938015920309679
56 Y0: -370863403886416141150505523919527677231932618184100095924
57 skew: 1.37
58 # f lognorm 3.93, alpha -1.96, score 1.97
59 # g lognorm 130.19, alpha -5.47, score 124.73
60 # f+g score 126.69
61
62 */
63
64 #include "cado.h" // IWYU pragma: keep
65 #include <float.h> // for DBL_MAX
66 #include <stdlib.h>
67 #include <math.h> // pow
68 #include <time.h>
69 #include <limits.h> // for ULONG_MAX
70 #include <stdio.h> // for fprintf, printf, stderr, fflush, stdout
71 #include <string.h> // for strcmp
72 #include <gmp.h> // for mpz_t, mpz_clear, mpz_init, gmp_printf, mpz_...
73 #include "cado_poly.h" // for cado_poly_fprintf_MurphyE, cado_poly
74 #include "macros.h" // for ASSERT_ALWAYS, ASSERT
75 #include "omp_proxy.h" // IWYU pragma: keep
76 #include "auxiliary.h"
77 #include "gcd.h" // gcd_uint64
78 #include "lll.h" // mat_Z, LLL
79 #include "mpz_poly.h"
80 #include "murphyE.h"
81 #include "rootfinder.h"
82 #include "timing.h" // for seconds
83 #include "usp.h" // numberOfRealRoots
84
85 /* We assume a difference <= ALPHA_BOUND_GUARD between alpha computed
86 with ALPHA_BOUND_SMALL and ALPHA_BOUND. In practice the largest value
87 observed is 0.79. */
88 #define ALPHA_BOUND_GUARD 1.0
89
90 /* global variables */
91 double best_score_f = DBL_MAX, worst_score_f = DBL_MIN, sum_score_f = 0;
92 unsigned long f_candidate = 0; /* number of irreducibility tests */
93 unsigned long f_irreducible = 0; /* number of irreducible polynomials */
94 double max_guard = DBL_MIN;
95 int skewed = 0; /* boolean (use skewed polynomial) */
96 unsigned long *count = NULL; /* use only with -skewed */
97 int easySM = 0; /* see the -easySM option */
98 int rrf = -1; /* see the -rrf option */
99 int rrg = -1; /* see the -rrg option */
100
101 /* MPZ_POLY_TIMINGS is defined (or maybe not) in utils/mpz_poly.h
102 */
103 #ifndef MPZ_POLY_H_
104 #error "please include mpz_poly.h first"
105 #endif
106 #ifdef MPZ_POLY_TIMINGS
107 static double timer[4] = {0.0, };
108 #endif
109
110 double Bf = 0.0, Bg = 0.0, Area = 0.0;
111 double bestE = 0.0; /* best Murphy-E so far */
112 int opt_flag = 0; /* 0: optimize "simple" E
113 1: optimize Muprhy E */
114
115 #define TIMER_ROOTS 0
116 #define TIMER_IRRED 1
117 #define TIMER_LLL 2
118 #define TIMER_MURPHYE 3
119
120 #ifdef MPZ_POLY_TIMINGS
121 #define START_TIMER double t = seconds_thread ()
122 #define END_TIMER(x) add_timer (x, seconds_thread () - t)
123
124 /* flag=0: computing roots of f mod p
125 1: checking irreducibility of f
126 2: LLL reduction
127 3: computing MurphyE */
128 static void
add_timer(int flag,double t)129 add_timer (int flag, double t)
130 {
131 #ifdef HAVE_OPENMP
132 #pragma omp critical
133 #endif
134 timer[flag] += t;
135 }
136 #else
137 #define START_TIMER
138 #define END_TIMER(x)
139 #endif
140
141
142 static int
check_SM(mpz_poly ff,mpz_t ell)143 check_SM (mpz_poly ff, mpz_t ell)
144 {
145 if (ff->deg <= 2) {
146 return 1;
147 }
148 if (ff->deg > 4) {
149 fprintf(stderr, "Not implemented\n");
150 ASSERT_ALWAYS(0);
151 }
152 // in degree 3 and 4, the minimum number of SMs is 1. We
153 // check that we have at least one root mod ell.
154 gmp_randstate_t rstate;
155 gmp_randinit_default(rstate);
156 int nr = mpz_poly_roots_mpz(NULL, ff, ell, rstate);
157 gmp_randclear(rstate);
158 return (nr >= 1);
159 }
160
161
162 /*
163 Print two nonlinear poly info. Return non-zero for record polynomials.
164 */
165 static int
print_nonlinear_poly_info(mpz_poly ff,double alpha_f,mpz_poly gg,int format,mpz_t n,mpz_t ell)166 print_nonlinear_poly_info (mpz_poly ff, double alpha_f, mpz_poly gg,
167 int format, mpz_t n, mpz_t ell)
168 {
169 unsigned int i;
170 double skew, logmu[2], alpha_g_approx, alpha_g, score, score_approx;
171 int df = ff->deg;
172 mpz_t *f = ff->coeff;
173 int dg = gg->deg;
174 mpz_t *g = gg->coeff;
175 static double best_score = DBL_MAX;
176 /* the coefficients of g are O(n^(1/df)) */
177
178 /* we use the skewness minimizing the sum of lognorms */
179 skew = L2_combined_skewness2 (ff, gg, SKEWNESS_DEFAULT_PREC);
180 logmu[1] = L2_lognorm (gg, skew);
181 logmu[0] = L2_lognorm (ff, skew);
182 /* first estimate alpha with a small bound */
183 alpha_g_approx = get_alpha (gg, ALPHA_BOUND_SMALL);
184
185 score_approx = logmu[1] + alpha_g_approx + logmu[0] + alpha_f;
186
187 if (score_approx >= best_score + ALPHA_BOUND_GUARD)
188 return 0;
189
190 /* now get a more precise alpha value */
191 alpha_g = get_alpha (gg, get_alpha_bound ());
192
193 score = logmu[1] + alpha_g + logmu[0] + alpha_f;
194 if (score_approx - score > max_guard)
195 #ifdef HAVE_OPENMP
196 #pragma omp critical
197 #endif
198 max_guard = score_approx - score;
199
200 double E = 0.0;
201 if (opt_flag == 0)
202 {
203 if (score >= best_score)
204 return 0; /* only print record scores */
205 }
206 else /* optimize Murphy-E */
207 {
208 if (score >= best_score + 1.0) /* the guard 1.0 seems good in
209 practice */
210 return 0;
211
212 /* compute Murphy-E */
213 cado_poly p;
214 START_TIMER;
215 p->pols[ALG_SIDE]->coeff = f;
216 p->pols[ALG_SIDE]->deg = df;
217 p->pols[RAT_SIDE]->coeff = g;
218 p->pols[RAT_SIDE]->deg = dg;
219 p->skew = skew;
220 E = MurphyE (p, Bf, Bg, Area, MURPHY_K, get_alpha_bound ());
221 END_TIMER (TIMER_MURPHYE);
222 if (E <= bestE)
223 return 0;
224 }
225
226 /* Possibly check the number of roots mod ell of f and g, assuming that
227 * they have the minimum number of real roots. */
228 if (easySM) {
229 if (! check_SM(ff, ell))
230 return 0;
231 if (! check_SM(gg, ell))
232 return 0;
233 }
234 bestE = E;
235
236 #ifdef HAVE_OPENMP
237 #pragma omp critical
238 #endif
239 {
240 if (score < best_score)
241 best_score = score;
242
243 if (format == 1)
244 gmp_printf ("n: %Zd\n", n);
245 else
246 gmp_printf ("N %Zd\n", n);
247
248 if (format == 1) {
249 for (i = df + 1; i -- != 0; )
250 gmp_printf ("c%u: %Zd\n", i, f[i]);
251 }
252 else {
253 for (i = df + 1; i -- != 0; )
254 gmp_printf ("X%u %Zd\n", i, f[i]);
255 }
256 if (format == 1) {
257 for (i = dg + 1; i -- != 0; )
258 gmp_printf ("Y%u: %Zd\n", i, g[i]);
259 }
260 else {
261 for (i = dg + 1; i -- != 0; )
262 gmp_printf ("Y%u %Zd\n", i, g[i]);
263 }
264 printf ("skew: %1.2f\n", skew);
265 int nr = numberOfRealRoots ((const mpz_t *) f, df, 0, 0, NULL);
266 printf ("# f lognorm %1.2f, alpha %1.2f, score %1.2f, %d rroot(s)\n",
267 logmu[0], alpha_f, logmu[0] + alpha_f, nr);
268 nr = numberOfRealRoots ((const mpz_t *) g, dg, 0, 0, NULL);
269 printf ("# g lognorm %1.2f, alpha %1.2f, score %1.2f, %d rroot(s)\n",
270 logmu[1], alpha_g, logmu[1] + alpha_g, nr);
271 printf ("# f+g score %1.2f\n", score);
272 if (opt_flag)
273 cado_poly_fprintf_MurphyE (stdout, E, Bf, Bg, Area, "");
274
275 printf ("\n");
276 fflush (stdout);
277 }
278 return 1;
279 }
280
281 /* return the number of polynomials we look for with -skewed:
282 * we compute the maximal skewness s = B^(-2/d)
283 * the coefficients of degree >= d/2 are bounded by B
284 * the coefficients of degree i < d/2 are bounded by B*s^(d/2-i) */
285 static unsigned long
get_maxtries(unsigned int B,unsigned int d)286 get_maxtries (unsigned int B, unsigned int d)
287 {
288 unsigned long maxtries = 1, c;
289 unsigned int i;
290 double skew;
291
292 ASSERT_ALWAYS(d >= 2);
293
294 skew = pow ((double) B, 2.0 / (double) d);
295 count = malloc ((d + 1) * sizeof (unsigned long));
296 for (i = 0; i <= d; i++)
297 {
298 if (i == d)
299 c = B; /* coefficient in 1..B */
300 else if (i == d - 1)
301 c = B + 1; /* coefficient in 0..B */
302 else if (2 * i >= d)
303 c = 2 * B + 1; /* coefficient in -B..B */
304 else /* i < d/2 */
305 {
306 c = B * (unsigned long) pow (skew, (double) (d - 2 * i) / 2.0);
307 if (i > 0)
308 c = 2 * c + 1;
309 else
310 c = 2 * c;
311 }
312 count[i] = c;
313 if (maxtries >= ULONG_MAX / count[i])
314 {
315 fprintf (stderr, "Error, too large -bound option\n");
316 exit (1);
317 }
318 maxtries *= count[i];
319 }
320 return maxtries;
321 }
322
323 static int
generate_f(mpz_t * f,unsigned int d,unsigned long idx,unsigned int bound)324 generate_f (mpz_t *f, unsigned int d, unsigned long idx, unsigned int bound)
325 {
326 unsigned int i, j;
327 int ok = 1;
328 int *a;
329
330 a = malloc ((d + 1) * sizeof (int));
331
332 /* the coefficient of degree j can take count[j] different values */
333 for (j = 0; j <= d; j++)
334 {
335 long k;
336 a[j] = idx % count[j];
337 idx = idx / count[j];
338 if (j == d)
339 a[j] ++; /* coefficient in [1..B] */
340 else if (j == 0) /* coefficient in [-k..k] except 0 */
341 {
342 k = count[j] / 2; /* count[j]=2*k */
343 /* 0..k-1 -> -k..-1
344 k..2k-1 -> 1..k */
345 a[j] = (a[j] < k) ? a[j] - k : a[j] - (k - 1);
346 }
347 else if (j < d - 1)
348 {
349 k = count[j] / 2; /* count[j]=2*k+1 */
350 a[j] -= k;
351 }
352 }
353 ASSERT_ALWAYS(idx == 0);
354
355 /* Check if polynomial agrees with maximal skewness: for i > d/2 and j < d/2,
356 the line going through |a[i]| and |a[j]| should not exceed bound at d/2 */
357 for (i = d / 2 + 1; i <= d && ok; i++)
358 {
359 unsigned ai = abs (a[i]);
360 if (ai == 0)
361 continue;
362 for (j = 0; j < (d + 1) / 2 && ok; j++)
363 {
364 unsigned aj = abs (a[j]);
365 double s = pow ((double) aj / (double) ai, 1.0 / (double) (i - j));
366 double mid = (double) ai * pow (s, (double) (2 * i - d) / 2.0);
367 ok = mid <= (double) bound;
368 }
369 }
370
371 /* Check if the reverse polynomial has smaller rank.
372 This test discards about 7.7% of the polynomials for d=4 and bound=6. */
373 if (ok)
374 for (j = 0; 2 * j < d; j++)
375 if (abs (a[d-j]) != abs(a[j]))
376 {
377 ok = abs (a[d-j]) < abs(a[j]);
378 break;
379 }
380
381 /* Since f(x) is equivalent to f(-x), if a[d-1]=0, then the largest a[d-3],
382 a[d-5], ... that is non zero should be positive. Discards 13.5% of the
383 remaining polynomials. */
384 if (ok)
385 for (j = d + 1; j >= 2;)
386 {
387 j -= 2;
388 if (a[j] != 0)
389 {
390 ok = a[j] > 0;
391 break;
392 }
393 }
394
395 /* Check if +-1 is a root of f. Discards 4.3% of the remaining polynomials. */
396 if (ok)
397 {
398 int value_one = 0, value_minus_one = 0;
399 for (j = 0; j <= d; j++)
400 {
401 value_one += a[j];
402 value_minus_one += (j & 1) ? -a[j] : a[j];
403 }
404 ok = value_one != 0 && value_minus_one != 0;
405 }
406
407 /* Content test. Discards 2.3% of the remaining polynomials. */
408 if (ok)
409 {
410 unsigned long g = a[d];
411 for (j = 0; j < d; j++)
412 g = gcd_int64 (g, a[j]);
413 ok = g == 1;
414 }
415
416 if (ok)
417 for (j = 0; j <= d; j++)
418 mpz_set_si (f[j], a[j]);
419
420 free (a);
421
422 return ok;
423 }
424
425 /*
426 Generate polynomial f(x) of degree d with rank 'idx',
427 with coefficients in [-bound, bound].
428 The coefficient of degree d can be taken in [1, bound],
429 due to the symmetry f(x) -> -f(x).
430 The coefficient of degree d-1 can be taken in [0, bound],
431 due to the symmetry f(x) -> f(-x).
432 The coefficient of degree 0 should not be 0.
433 Thus there are bound*(bound+1)*(2*bound+1)^(d-2)*(2*bound) possible values.
434 Return 0 if the corresponding polynomial is not irreducible.
435 Return !=0 if the poly is valid.
436 */
437 static int
polygen_JL_f(int d,unsigned int bound,mpz_t * f,unsigned long idx)438 polygen_JL_f (int d, unsigned int bound, mpz_t *f, unsigned long idx)
439 {
440 int ok = 1;
441 /* compute polynomial of index idx and check it is irreducible */
442 mpz_poly ff;
443 ff->deg = d;
444 ff->coeff = f;
445 //mpz_poly_init(ff, d);
446 if (!skewed) {
447 int max_abs_coeffs;
448 unsigned long next_counter;
449 ok = mpz_poly_setcoeffs_counter(ff, &max_abs_coeffs, &next_counter,
450 d, idx, bound);
451 } else {
452 ok = generate_f (f, d, idx, bound);
453 }
454
455 // to be compatible with the previous version, the count on the number of valid polys was before
456 // the content test.
457 // Now the content test is included in mpz_poly_setcoeffs_counter()
458 // so the number of candidates will be lower than before.
459 // to get the previous value, use
460 // #define POLY_CONTENT -6
461 // if ((ok == 1) || ((ok == 0) && (max_abs_coeffs == POLY_CONTENT))){
462 // #undef POLY_CONTENT
463 if (ok)
464 #ifdef HAVE_OPENMP
465 #pragma omp critical
466 #endif
467 f_candidate ++;
468
469 /* irreducibility test */
470 if (ok)
471 {
472 START_TIMER;
473 ok = mpz_poly_squarefree_p (ff);
474 /* check number of real roots */
475 if (ok && (easySM || rrf != -1))
476 {
477 int nr = numberOfRealRoots ((const mpz_t *) ff->coeff, ff->deg, 0.0, 0, NULL);
478 if (easySM)
479 /* check that the number of real roots is minimal */
480 ok = nr == (ff->deg & 1);
481 else
482 ok = nr == rrf;
483 }
484 if (ok)
485 ok = mpz_poly_is_irreducible_z (ff);
486 END_TIMER (TIMER_IRRED);
487 }
488
489 return ok;
490 }
491
492 /* Generate polynomial g(x) of degree dg, given root 'root' of f mod N.
493 It might be better to take into account the skewness of f in the LLL
494 lattice, but experimentally this does not give better results (probably
495 because LLL is not very sensible to a small change of the skewness).
496 kN is the product k*N, where k is the multiplier. */
497 static void
polygen_JL_g(mpz_t kN,int dg,mat_Z g,mpz_t root,double skew_f)498 polygen_JL_g (mpz_t kN, int dg, mat_Z g, mpz_t root, double skew_f)
499 {
500 int i, j;
501 mpz_t a, b, det, r;
502 unsigned long skew, skew_powi;
503
504 skew = skew_f < 0.5 ? 1 : round (skew_f);
505
506 mpz_init (det);
507 mpz_init_set_ui (a, 1);
508 mpz_init_set_ui (b, 1);
509 mpz_init_set (r, root);
510 for (i = 0; i <= dg + 1; i ++) {
511 for (j = 0; j <= dg + 1; j ++) {
512 mpz_set_ui (g.coeff[i][j], 0);
513 }
514 }
515
516 for (i = skew_powi = 1; i <= dg + 1; i++) {
517 for (j = 1; j <= dg + 1; j++) {
518 if (i == 1)
519 {
520 if (j == 1)
521 mpz_set (g.coeff[j][i], kN);
522 else
523 {
524 mpz_neg (g.coeff[j][i], r);
525 mpz_mul (r, r, root);
526 }
527 }
528 else if (i == j)
529 {
530 ASSERT_ALWAYS((double) skew_powi * (double) skew < (double) ULONG_MAX);
531 skew_powi *= skew;
532 mpz_set_ui (g.coeff[j][i], skew_powi);
533 }
534 }
535 }
536
537 START_TIMER;
538 LLL (det, g, NULL, a, b);
539 END_TIMER (TIMER_LLL);
540
541 /* divide row i back by skew^i */
542 skew_powi = 1;
543 for (i = 2; i <= dg + 1; i++)
544 {
545 skew_powi *= skew;
546 for (j = 1; j <= dg + 1; j++)
547 {
548 ASSERT_ALWAYS (mpz_divisible_ui_p (g.coeff[j][i], skew_powi));
549 mpz_divexact_ui (g.coeff[j][i], g.coeff[j][i], skew_powi);
550 }
551 }
552
553 mpz_clear (det);
554 mpz_clear (a);
555 mpz_clear (b);
556 mpz_clear (r);
557 }
558
559 /* JL method to generate degree d and d-1 polynomials.
560 Given irreducible polynomial f of degree df, find roots of f mod n,
561 and for each root, use Joux-Lercier method to find good polynomials g. */
562 static void
polygen_JL2(mpz_t n,unsigned int df,unsigned int dg,unsigned long nb_comb,mpz_poly f,long bound2,mpz_t ell)563 polygen_JL2 (mpz_t n,
564 unsigned int df, unsigned int dg,
565 unsigned long nb_comb, mpz_poly f, long bound2, mpz_t ell)
566 {
567 unsigned int i, j, nr, format = 1;
568 mpz_t *rf, c;
569 mat_Z g;
570 mpz_poly *v, u;
571 long *a;
572 double alpha_f;
573 gmp_randstate_t rstate;
574 gmp_randinit_default(rstate);
575
576 ASSERT_ALWAYS (df >= 3);
577 mpz_init (c);
578 rf = (mpz_t *) malloc (df * sizeof(mpz_t));
579 for (i = 0; i < df; i ++)
580 mpz_init (rf[i]);
581 g.coeff = (mpz_t **) malloc ((dg + 2) * sizeof(mpz_t*));
582 g.NumRows = g.NumCols = dg + 1;
583 for (i = 0; i <= dg + 1; i ++) {
584 g.coeff[i] = (mpz_t *) malloc ((dg + 2) * sizeof(mpz_t));
585 for (j = 0; j <= dg + 1; j ++) {
586 mpz_init (g.coeff[i][j]);
587 }
588 }
589 v = malloc ((dg + 1) * sizeof (mpz_poly));
590 for (j = 0; j <= dg; j++)
591 {
592 v[j]->deg = dg;
593 v[j]->coeff = g.coeff[j + 1] + 1;
594 }
595 mpz_poly_init (u, dg);
596 u->deg = dg;
597 a = malloc ((dg + 1) * sizeof (long));
598
599 /* compute roots of the polynomial f */
600 START_TIMER;
601 nr = mpz_poly_roots_mpz (rf, f, n, rstate);
602 END_TIMER (TIMER_ROOTS);
603 ASSERT(nr <= df);
604
605 /* update the best and worst score for f (FIXME: even if f has no roots?) */
606 double skew_f, lognorm_f, score_f;
607 alpha_f = get_alpha (f, get_alpha_bound ());
608 skew_f = L2_skewness (f, SKEWNESS_DEFAULT_PREC);
609 lognorm_f = L2_lognorm (f, skew_f);
610 score_f = lognorm_f + alpha_f;
611
612 if (score_f < best_score_f)
613 #ifdef HAVE_OPENMP
614 #pragma omp critical
615 #endif
616 best_score_f = score_f;
617
618 if (score_f > worst_score_f)
619 #ifdef HAVE_OPENMP
620 #pragma omp critical
621 #endif
622 worst_score_f = score_f;
623
624 #ifdef HAVE_OPENMP
625 #pragma omp critical
626 #endif
627 {
628 f_irreducible ++;
629 sum_score_f += score_f;
630 }
631
632 /* for each root of f mod n, generate the corresponding g */
633 for (i = 0; i < nr; i ++) {
634 /* generate g of degree dg */
635 polygen_JL_g (n, dg, g, rf[i], skew_f);
636
637 /* we skip idx = 0 which should correspond to c[0] = ... = c[dg] = 0 */
638 for (unsigned long idx = 1; idx < nb_comb; idx ++)
639 {
640 unsigned long k = idx;
641
642 /* compute first index */
643 a[0] = k % (bound2 + 1);
644 k = k / (bound2 + 1);
645 for (j = 1; j <= dg; j++)
646 {
647 a[j] = k % (2 * bound2 + 1);
648 k = k / (2 * bound2 + 1);
649 a[j] = (a[j] <= bound2) ? a[j] : a[j] - (2 * bound2 + 1);
650 }
651 ASSERT_ALWAYS(k == 0);
652 mpz_poly_mul_si (u, v[0], a[0]);
653 for (j = 1; j <= dg; j++)
654 mpz_poly_addmul_si (u, v[j], a[j]);
655
656 /* adjust degree of u */
657 for (u->deg = dg; u->deg >= 0 && mpz_cmp_ui (u->coeff[u->deg], 0)
658 == 0; u->deg--);
659
660 #if 0
661 /* If u is not square-free or irreducible, skip it. However, this
662 test is very expensive, and non-irreducible polynomials should
663 not happen in practice for large input N, thus we disable. */
664 if (mpz_cmp_ui (u->coeff[0], 0) == 0 || !mpz_poly_squarefree_p (u)
665 || !mpz_poly_is_irreducible_z (u))
666 continue;
667 #endif
668 /* check the number of real roots of g */
669 if (easySM || rrg != -1)
670 {
671 int nr = numberOfRealRoots ((const mpz_t *) u->coeff, u->deg, 0.0, 0, NULL);
672 int ok;
673 if (easySM)
674 ok = nr == (u->deg & 1);
675 else
676 ok = nr == rrg;
677 if (! ok)
678 continue; // skip this g
679 }
680
681 if (print_nonlinear_poly_info (f, alpha_f, u, format, n, ell))
682 {
683 #if 0 /* print coefficients of record combination */
684 for (j = 0; j <= dg; j++)
685 printf ("%ld ", a[j]);
686 printf ("\n");
687 #endif
688 }
689 }
690 }
691
692 /* clear */
693 gmp_randclear(rstate);
694 free (a);
695 mpz_poly_clear (u);
696 free (v);
697 for (i = 0; i < df; i ++)
698 mpz_clear (rf[i]);
699 for (i = 0; i <= dg + 1; i ++) {
700 for (j = 0; j <= dg + 1; j ++)
701 mpz_clear (g.coeff[i][j]);
702 free(g.coeff[i]);
703 }
704 free (g.coeff);
705 free (rf);
706 mpz_clear (c);
707 }
708
709 /* JL method to generate d and d-1 polynomial.
710 Generate polynomial f of degree df with |f[i]| <= bound and index 'idx'. */
711 static void
polygen_JL1(mpz_t n,unsigned int df,unsigned int dg,unsigned int bound,unsigned long idx,unsigned long nb_comb,unsigned int bound2,mpz_t ell)712 polygen_JL1 (mpz_t n,
713 unsigned int df, unsigned int dg, unsigned int bound,
714 unsigned long idx, unsigned long nb_comb, unsigned int bound2,
715 mpz_t ell)
716 {
717 unsigned int i;
718 mpz_t *f;
719 mpz_poly ff;
720 int irred;
721
722 ASSERT_ALWAYS (df >= 3);
723 f = (mpz_t *) malloc ((df + 1)*sizeof(mpz_t));
724 for (i = 0; i <= df; i ++)
725 mpz_init (f[i]);
726
727 ff->deg = df;
728 ff->coeff = f;
729
730 /* generate f of degree d with small coefficients */
731 irred = polygen_JL_f (df, bound, f, idx);
732 if (irred)
733 polygen_JL2 (n, df, dg, nb_comb, ff, bound2, ell);
734 /* clear */
735 for (i = 0; i <= df; i ++)
736 mpz_clear (f[i]);
737 free (f);
738 }
739
740 static void
usage()741 usage ()
742 {
743 fprintf (stderr, "./dlpolyselect -N xxx -df xxx -dg xxx -bound xxx [-modr xxx] [-modm xxx] [-t xxx] [-easySM <ell>] [-skewed] [-rrf nnn]\n");
744 fprintf (stderr, "Mandatory parameters:\n");
745 fprintf (stderr, " -N xxx input number\n");
746 fprintf (stderr, " -df xxx degree of polynomial f\n");
747 fprintf (stderr, " -dg xxx degree of polynomial g\n");
748 fprintf (stderr, " -bound xxx bound for absolute value of coefficients of f\n");
749 fprintf (stderr, "Optional parameters:\n");
750 fprintf (stderr, " -modr r -modm m processes only polynomials of index r mod m\n");
751 fprintf (stderr, " -t nnn uses n threads\n");
752 fprintf (stderr, " -easySM ell generates polynomials with minimal number of SMs mod ell\n");
753 fprintf (stderr, " -skewed search for skewed polynomials\n");
754 fprintf (stderr, " -rrf nnn f should have nnn real roots\n");
755 fprintf (stderr, " -rrg nnn g should have nnn real roots\n");
756 fprintf (stderr, " -Bf nnn sieving bound for f\n");
757 fprintf (stderr, " -Bg nnn sieving bound for g\n");
758 fprintf (stderr, " -area nnn sieving area\n");
759 exit (1);
760 }
761
762 int
main(int argc,char * argv[])763 main (int argc, char *argv[])
764 {
765 int i;
766 mpz_t N;
767 mpz_t ell;
768 unsigned int df = 0, dg = 0;
769 int nthreads = 1;
770 unsigned int bound = 4; /* bound on the coefficients of f */
771 unsigned long maxtries;
772 double t;
773 unsigned long modr = 0, modm = 1;
774
775 t = seconds ();
776 mpz_init (N);
777 mpz_init (ell);
778
779 /* printf command-line */
780 printf ("#");
781 for (i = 0; i < argc; i++)
782 printf (" %s", argv[i]);
783 printf ("\n");
784 fflush (stdout);
785
786 /* parsing */
787 while (argc >= 2 && argv[1][0] == '-')
788 {
789 if (argc >= 3 && strcmp (argv[1], "-N") == 0) {
790 mpz_set_str (N, argv[2], 10);
791 argv += 2;
792 argc -= 2;
793 }
794 else if (argc >= 3 && strcmp (argv[1], "-easySM") == 0) {
795 mpz_set_str (ell, argv[2], 10);
796 easySM = 1;
797 argv += 2;
798 argc -= 2;
799 }
800 else if (argc >= 3 && strcmp (argv[1], "-df") == 0) {
801 df = atoi (argv[2]);
802 argv += 2;
803 argc -= 2;
804 }
805 else if (argc >= 3 && strcmp (argv[1], "-dg") == 0) {
806 dg = atoi (argv[2]);
807 argv += 2;
808 argc -= 2;
809 }
810 else if (argc >= 3 && strcmp (argv[1], "-bound") == 0) {
811 bound = atoi (argv[2]);
812 argv += 2;
813 argc -= 2;
814 }
815 else if (argc >= 3 && strcmp (argv[1], "-modm") == 0) {
816 modm = strtoul (argv[2], NULL, 10);
817 argv += 2;
818 argc -= 2;
819 }
820 else if (argc >= 3 && strcmp (argv[1], "-modr") == 0) {
821 modr = strtoul (argv[2], NULL, 10);
822 argv += 2;
823 argc -= 2;
824 }
825 else if (argc >= 3 && strcmp (argv[1], "-t") == 0) {
826 nthreads = atoi (argv[2]);
827 argv += 2;
828 argc -= 2;
829 }
830 /* if rrf = -1 (default), f might have any number of real roots,
831 otherwise it should have exactly 'rrf' real roots */
832 else if (argc >= 3 && strcmp (argv[1], "-rrf") == 0) {
833 rrf = atoi (argv[2]);
834 argv += 2;
835 argc -= 2;
836 }
837 /* if rrg = -1 (default), g might have any number of real roots,
838 otherwise it should have exactly 'rrg' real roots */
839 else if (argc >= 3 && strcmp (argv[1], "-rrg") == 0) {
840 rrg = atoi (argv[2]);
841 argv += 2;
842 argc -= 2;
843 }
844 else if (argc >= 2 && strcmp (argv[1], "-skewed") == 0) {
845 skewed = 1;
846 argv += 1;
847 argc -= 1;
848 }
849 else if (argc >= 3 && strcmp (argv[1], "-Bf") == 0) {
850 Bf = atof (argv[2]);
851 argv += 2;
852 argc -= 2;
853 }
854 else if (argc >= 3 && strcmp (argv[1], "-Bg") == 0) {
855 Bg = atof (argv[2]);
856 argv += 2;
857 argc -= 2;
858 }
859 else if (argc >= 3 && strcmp (argv[1], "-area") == 0) {
860 Area = atof (argv[2]);
861 argv += 2;
862 argc -= 2;
863 }
864 else {
865 fprintf (stderr, "Invalid option: %s\n", argv[1]);
866 usage();
867 exit (1);
868 }
869 }
870
871 if (mpz_cmp_ui (N, 0) <= 0) {
872 fprintf (stderr, "Error, missing input number (-N option)\n");
873 usage ();
874 }
875
876 if (df == 0) {
877 fprintf (stderr, "Error, missing degree (-df option)\n");
878 usage ();
879 }
880
881 if (dg == 0 || dg >= df) {
882 fprintf (stderr, "Error, missing or erroneous degree (-dg option): ");
883 fprintf (stderr, "one should have dg < df.\n");
884 usage ();
885 }
886
887 opt_flag = Bf != 0 && Bg != 0 && Area != 0;
888
889 srand (time (NULL));
890
891 ASSERT_ALWAYS (bound >= 1);
892
893 if (skewed)
894 maxtries = get_maxtries (bound, df);
895 else {
896 /* check modm has no common factor with B, B+1, 2B+1 and 2B to avoid
897 a bias between classes mod 'modm' */
898 if (gcd_uint64 (modm, 2 * bound) != 1 ||
899 gcd_uint64 (modm, bound + 1) != 1 ||
900 gcd_uint64 (modm, 2 * bound + 1) != 1)
901 {
902 fprintf (stderr, "Error, modm should be coprime to "
903 "2*bound*(bound+1)*(2*bound+1)\n");
904 exit (1);
905 }
906 double maxtries_double = (double) bound;
907 maxtries_double *= (double) (bound + 1);
908 maxtries_double *= pow ((double) (2 * bound + 1), (double) (df - 2));
909 maxtries_double *= (double) (2 * bound);
910 if (maxtries_double >= (double) ULONG_MAX)
911 maxtries = ULONG_MAX;
912 else
913 maxtries = (unsigned long) maxtries_double;
914 }
915
916 unsigned int bound2 = 1; /* bound on the coefficients of linear
917 combinations from the LLL short vectors */
918 double nb_comb_f;
919 unsigned long nb_comb;
920 /* we try all combinations u = c[0]*v[0] + ... + c[dg]*v[dg] with
921 -bound2 <= c[j] <= bound2, except 0 <= c[0] <= bound2 since u and -u
922 are equivalent, and except c[0] = ... = c[dg] = 0, this gives a total
923 of (bound2+1)*(2*bound2+1)^dg-1 */
924 nb_comb_f = (double) (bound2 + 1) * pow ((double) (2 * bound2 + 1),
925 (double) dg);
926 nb_comb = (nb_comb_f > (double) ULONG_MAX) ? ULONG_MAX
927 : (unsigned long) nb_comb_f;
928
929 printf ("# will generate about %lu polynomials\n", maxtries / modm);
930
931 #ifdef HAVE_OPENMP
932 omp_set_num_threads (nthreads);
933 #pragma omp parallel for schedule(dynamic)
934 #else
935 if (nthreads > 1) {
936 fprintf(stderr, "Warning: openmp unavailable, -t ignored\n");
937 }
938 #endif
939 for (unsigned long c = modr; c < maxtries; c += modm)
940 polygen_JL1 (N, df, dg, bound, c, nb_comb, bound2, ell);
941
942 t = seconds () - t;
943
944 printf ("# found %lu irreducible f out of %lu candidates out of %lu\n",
945 f_irreducible, f_candidate, maxtries / modm);
946 printf ("# best f-score %1.2f, av. %1.2f, worst %1.2f, max alpha-guard %1.2f\n",
947 best_score_f, sum_score_f / f_irreducible, worst_score_f,
948 max_guard);
949 if (max_guard > ALPHA_BOUND_GUARD)
950 printf ("# Warning: max_guard > ALPHA_BOUND_GUARD, might "
951 "have missed some polynomials\n");
952 printf ("# Time %.2fs", t);
953 #ifdef MPZ_POLY_TIMINGS
954 printf (" (roots %.2fs, irred %.2fs, lll %.2fs, MurphyE %.2fs)",
955 timer[TIMER_ROOTS], timer[TIMER_IRRED], timer[TIMER_LLL],
956 timer[TIMER_MURPHYE]);
957 printf ("\n#");
958 print_timings_pow_mod_f_mod_p();
959 #endif
960 printf ("\n");
961
962 mpz_clear (N);
963 mpz_clear (ell);
964 if (skewed)
965 free (count);
966
967 return 0;
968 }
969