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