1 #include "cado.h" // IWYU pragma: keep
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <string.h>
5 #include <math.h>
6 #include <float.h> // DBL_MAX
7 #include <gmp.h>                // for mpz_t, mpz_clear, mpz_init, mpz_submu...
8 #include <stdint.h>             // for uint32_t
9 #include "cado_poly.h"          // for cado_poly_init, cado_poly_read, cado_...
10 #include "macros.h"             // for ASSERT_ALWAYS
11 #include "mpz_poly.h"           // for mpz_poly_s, mpz_poly, mpz_poly_discri...
12 #include "auxiliary.h" /* for common routines with polyselect.c */
13 #include "gmp_aux.h"    // ulong_isprime
14 #include "mod_ul.h"     // modulusul_t
15 #include "modul_poly.h"
16 #include "size_optimization.h"
17 #include "timing.h"             // for seconds
18 
19 /* for the rotation, we try (j*x+k) for |k| <= 2^MAX_k */
20 int MAX_k = 16;
21 
22 // #define KK -16346
23 // #define PP 7
24 
25 /* if gpn and pn are not coprime, let g = gcd(gpn, pn):
26    we must have g*g'*v = r + j*g*p' with g' = gpn/g and p'=pn/g:
27    (a) if g does not divide r, there is no solution
28    (b) if g divides r, then we must have:
29    g'*v = r' + j*p' with r'=r/g, let v0 = r'/g' mod p',
30    then v0, v0+p', ..., v0+(g-1)*p' are the g solutions. */
31 static void
special_update(double * A,long K0,long K1,const residueul_t gpn,modulusul_t Pn,residueul_t r,double alpha_p)32 special_update (double *A, long K0, long K1, const residueul_t gpn,
33 		modulusul_t Pn, residueul_t r, double alpha_p)
34 {
35   modulusul_t Pp;
36   residueul_t rp, gp;
37   int ret;
38   long k0, k;
39   modintul_t gg;
40   unsigned long g, rr, pp, pn;
41 
42   pn = modul_getmod_ul (Pn);
43   modul_gcd (gg, gpn, Pn);
44   g = gg[0];
45   rr = modul_get_ul (r, Pn);
46   if (rr % g == 0)
47     {
48       pp = pn / g;
49       modul_initmod_ul (Pp, pp);
50       modul_init (gp, Pp);
51       modul_init (rp, Pp);
52       modul_set_ul (gp, modul_get_ul (gpn, Pn) / g, Pp);
53       modul_set_ul (rp, rr / g, Pp);
54       ret = modul_inv (gp, gp, Pp);
55       ASSERT_ALWAYS (ret != 0);
56       modul_mul (rp, rp, gp, Pp);
57       k0 = modul_get_ul (rp, Pp);
58       for (k = k0; k <= K1; k += pp)
59 	{
60 	  A[k - K0] -= alpha_p;
61 #if (defined(KK) && defined(PP))
62 	  if (k == KK && (pn % PP) == 0)
63 	    fprintf (stderr, "subtract %f to AA[%d] for pp=%lu, now %f\n",
64 		     alpha_p, KK, pp, A[k - K0]);
65 #endif
66 	}
67       for (k = k0 - pp; k >= K0; k -= pp)
68 	{
69 	  A[k - K0] -= alpha_p;
70 #if (defined(KK) && defined(PP))
71 	  if (k == KK && (pn % PP) == 0)
72 	    fprintf (stderr, "subtract %f to AA[%d] for pp=%lu, now %f\n",
73 		     alpha_p, KK, pp, A[k - K0]);
74 #endif
75 	}
76       modul_clear (gp, Pp);
77       modul_clear (rp, Pp);
78       modul_clearmod (Pp);
79     }
80 }
81 
82 /* f(x) of degree d is the polynomial f0(x) + j*x*g(x)
83    where g(x) = b*x-m.
84    A is a table of K1-K0+1 entries, where A[k-K0] corresponds to the alpha
85    contribution for f(x) + k*g(x)
86    p is the current prime, we consider pn = p^n here */
87 void
update_table(mpz_t * f,int d,mpz_t m,mpz_t b,double * A,long K0,long K1,unsigned long pn,double alpha_p)88 update_table (mpz_t *f, int d, mpz_t m, mpz_t b, double *A, long K0, long K1,
89 	      unsigned long pn, double alpha_p)
90 {
91   long k, k0;
92   modul_poly_t fpn; /* f mod p^n */
93   modulusul_t Pn;
94   residueul_t l, gpn, bpn, r, v;
95   int ret;
96 
97   modul_initmod_ul (Pn, pn);
98   modul_init (gpn, Pn);
99   modul_init (bpn, Pn);
100   modul_init (r, Pn);
101   modul_init (l, Pn);
102   modul_init (v, Pn);
103 
104   modul_poly_init (fpn, d);
105 
106   /* first reduce f(x) and g(x) mod p^n */
107   mpz_poly F;
108   F->coeff = f;
109   F->deg = d;
110   modul_poly_set_mod_raw (fpn, F, Pn);
111 
112   modul_set_ul (gpn, mpz_fdiv_ui (m, pn), Pn);
113   /* invariant: gpn = -g(l) */
114   modul_set_ul (bpn, mpz_fdiv_ui (b, pn), Pn);
115   modul_set_ul (l, 0, Pn);
116 
117   for (;;)
118     {
119       modul_poly_eval (r, fpn, l, Pn);
120       /* invariant: gpn = -g(l) */
121       /* we want r = v*gpn, i.e., v = r/gpn; r is never 0 otherwise f(x) would
122 	 be divisible by p^n, but gpn = g(l) can be 0 */
123       ret = modul_intcmp_ul (gpn, 0);
124       if (ret != 0)
125 	{
126 	  ret = modul_inv (v, gpn, Pn); /* FIXME: use batch inversion */
127 	  if (ret == 0) /* gpn and pn are not coprime */
128 	    special_update (A, K0, K1, gpn, Pn, r, alpha_p);
129 	  else
130 	    {
131 	      modul_mul (v, v, r, Pn);
132 
133 	      k0 = modul_get_ul (v, Pn);
134 	      for (k = k0; k <= K1; k += pn)
135 		{
136 		  A[k - K0] -= alpha_p;
137 #if (defined(KK) && defined(PP))
138 		  if (k == KK && (pn % PP) == 0)
139 		    fprintf (stderr, "subtract %f to AA[%d] for l=%lu, pn=%lu, now %f\n",
140 			     alpha_p, KK, modul_get_ul (l, Pn), pn, A[k - K0]);
141 #endif
142 		}
143 	      for (k = k0 - (long) pn; k >= K0; k -= (long) pn)
144 		{
145 		  A[k - K0] -= alpha_p;
146 #if (defined(KK) && defined(PP))
147 		  if (k == KK && (pn % PP) == 0)
148 		    fprintf (stderr, "subtract %f to AA[%d] for l=%lu, pn=%lu, now %f\n",
149 			     alpha_p, KK, modul_get_ul (l, Pn), pn, A[k - K0]);
150 #endif
151 		}
152 	    }
153 	}
154 
155       /* since g(x) = b*x-m, -g(l) = m-b*l */
156       modul_sub (gpn, gpn, bpn, Pn);
157       modul_add_ul (l, l, 1, Pn);
158       if (modul_intcmp_ul (l, 0) == 0)
159 	break;
160     }
161 
162   modul_clear (gpn, Pn);
163   modul_clear (bpn, Pn);
164   modul_clear (r, Pn);
165   modul_clear (l, Pn);
166   modul_clear (v, Pn);
167   modul_clearmod (Pn);
168   modul_poly_clear (fpn);
169 }
170 
171 /* D <- discriminant (f+k*g), which has degree d */
172 static void
discriminant_k(mpz_t * D,mpz_poly_ptr f,mpz_t m,mpz_t b)173 discriminant_k (mpz_t *D, mpz_poly_ptr f, mpz_t m, mpz_t b)
174 {
175   int i, j, k;
176   uint32_t **M, pivot;
177   int d = f->deg;
178 
179   ASSERT_ALWAYS(d <= 9);
180 
181   /* we first put in D[i] the value of disc(f + i*g) for 0 <= i <= d,
182      thus if disc(f + k*g) = a[d]*k^d + ... + a[0], then
183           D[0] = a[0]
184           D[1] = a[0] + a[1] + ... + a[d]
185           ...
186           D[d] = a[0] + a[1]*d + ... + a[d]*d^d */
187 
188   mpz_poly_discriminant (D[0], f);
189   for (i = 1; i <= d; i++)
190     {
191       /* add b*x - m */
192       mpz_add (f->coeff[1], f->coeff[1], b);
193       mpz_sub (f->coeff[0], f->coeff[0], m);
194       mpz_poly_discriminant (D[i], f);
195     }
196 
197   /* initialize matrix coefficients */
198   M = (uint32_t**) malloc ((d + 1) * sizeof(uint32_t*));
199   for (i = 0; i <= d; i++)
200     {
201       M[i] = (uint32_t*) malloc ((d + 1) * sizeof(uint32_t));
202       M[i][0] = 1;
203       for (j = 1; j <= d; j++)
204         M[i][j] = i * M[i][j-1];
205     }
206 
207   for (j = 0; j < d; j++)
208     {
209       /* invariant: D[i] = M[i][0] * a[0] + ... + M[i][d] * a[d]
210          with M[i][k] = 0 for k < j and k < i */
211       for (i = j + 1; i <= d; i++)
212         {
213           /* eliminate M[i][j] */
214           pivot = M[i][j] / M[j][j];
215           mpz_submul_ui (D[i], D[j], pivot);
216           for (k = j; k <= d; k++)
217             M[i][k] -= pivot * M[j][k];
218         }
219     }
220 
221   /* now we have an upper triangular matrix */
222   for (j = d; j > 0; j--)
223     {
224       for (k = j + 1; k <= d; k++)
225         mpz_submul_ui (D[j], D[k], M[j][k]);
226       ASSERT_ALWAYS(mpz_divisible_ui_p (D[j], M[j][j]));
227       mpz_divexact_ui (D[j], D[j], M[j][j]);
228     }
229 
230   /* restore the original f[] */
231   mpz_submul_ui (f->coeff[1], b, d);
232   mpz_addmul_ui (f->coeff[0], m, d);
233 
234   for (i = 0; i <= d; i++)
235     free (M[i]);
236   free (M);
237 }
238 
239 unsigned long
rotate_area(long K0,long K1,long J0,long J1)240 rotate_area (long K0, long K1, long J0, long J1)
241 {
242   return (unsigned long) (J1 - J0 + 1) * (unsigned long) (K1 - K0 + 1);
243 }
244 
245 /* Use Emmanuel Thome's idea: assuming alpha(f + k*g) admits a Gaussian
246    distribution with mean 'm' and standard deviation 's', then the probability
247    that one polynomial has a value of alpha >= A is
248    1/2*(1 - erf((A-m)/s/sqrt(2))), thus the probability that K polynomials
249    have all their alpha values >= A is [1/2*(1 - erf((A-m)/s/sqrt(2)))]^K.
250    For 'm' and 's' fixed, and a given K, we define the value of A for which
251    this probability is 1/2 to be the expected value of A for K polynomials.
252 
253    We assume lognorm(f + k*g) + E(alpha(f + k*g)) is first decreasing, then
254    increasing, then the optimal K corresponds to the minimum of that function.
255 */
256 void
rotate_bounds(mpz_poly_ptr f,mpz_t b,mpz_t m,long * K0,long * K1,long * J0,long * J1,int verbose)257 rotate_bounds (mpz_poly_ptr f, mpz_t b, mpz_t m, long *K0, long *K1,
258                long *J0, long *J1, int verbose)
259 {
260   /* exp_alpha[i] corresponds to K=2^i polynomials for m=0, s=1
261      f := x -> 1/2*(1 - erf(x/sqrt(2)));
262      seq(fsolve(f(x)^(2^k) = 1/2, x=-log[2](1.0+k)), k=0..30);
263    */
264   double exp_alpha[] = {0.0 /* K=1 */, -0.5449521356 /* 2 */,
265                         -0.9981488825 /* 4 */, -1.385198061 /* 8 */,
266                         -1.723526050 /* 16 */, -2.025111894 /* 32 */,
267                         -2.298313131 /* 64 */, -2.549067054 /* 128 */,
268                         -2.781676726 /* 256 */, -2.999326227 /* 512 */,
269                         -3.204420841 /* 1024 */, -3.398814100 /* 2048 */,
270                         -3.583961388 /* 4096 */, -3.761025425 /* 8192 */,
271                         -3.930949902 /* 16384 */, -4.094511733 /* 32768 */,
272                         -4.252358774 /* 65536 */, -4.405037486 /* 131072 */,
273                         -4.553013560 /* 262144 */, -4.696687518 /* 524288 */,
274                         -4.836406692 /* 2^20 */, -4.972474538 /* 2^21 */,
275                         -5.105157963 /* 2^22 */, -5.234693169 /* 2^23 */,
276                         -5.361290351 /* 2^24 */, -5.485137511 /* 2^25 */,
277                         -5.606403590 /* 2^26 */, -5.725241052 /* 2^27 */,
278                         -5.841788041 /* 2^28 */, -5.956170181 /* 2^29 */,
279                         DBL_MAX};
280   int i;
281   long k0 = 0, j0 = 0;
282   double lognorm, alpha, E0, E, best_E;
283   double skewness = L2_skewness (f, SKEWNESS_DEFAULT_PREC);
284   long jmax = (long) ((double) (1L << MAX_k) / skewness);
285   unsigned long max_area = 1UL << MAX_k;
286 
287 #define MARGIN 0.12 /* we allow a small error margin in the expected lognorm
288                        + alpha values, to get a larger search range */
289 
290   E0 = L2_lognorm (f, skewness);
291 
292   *K0 = -2;
293   *J0 = -16;
294   *J1 = 16;
295 
296   /* look for negative k: -2, -4, -8, ... */
297   best_E = E0;
298   for (i = 1; rotate_area (*K0, -*K0, *J0, *J1) < max_area; i++, *K0 *= 2)
299     {
300       k0 = rotate_aux (f->coeff, b, m, k0, *K0, 0);
301       lognorm = L2_skew_lognorm (f, SKEWNESS_DEFAULT_PREC);
302       alpha = exp_alpha[i];
303       E = lognorm + alpha;
304       if (E < best_E + MARGIN)
305         {
306           if (E < best_E)
307             best_E = E;
308         }
309       else
310         break;
311     }
312   /* go back to k=0 */
313   k0 = rotate_aux (f->coeff, b, m, k0, 0, 0);
314   *K1 = -*K0;
315 
316   /* now try negative j: -1, -3, -7, ... */
317   for (i++; exp_alpha[i] != DBL_MAX && rotate_area (*K0, *K1, *J0, -*J0) < max_area; i++, *J0 = 2 * *J0 - 1)
318     {
319       j0 = rotate_aux (f->coeff, b, m, j0, *J0, 1);
320       lognorm = L2_skew_lognorm (f, SKEWNESS_DEFAULT_PREC);
321       alpha = exp_alpha[i];
322       E = lognorm + alpha;
323       if (E < best_E + MARGIN)
324         {
325           if (E < best_E)
326             best_E = E;
327         }
328       else
329         break;
330       if (1 - 2 * *J0 > jmax)
331         break;
332     }
333   *J1 = -*J0;
334 
335   if (verbose > 0)
336     fprintf (stderr, "# Rotate bounds: %ld <= j <= %ld, %ld <= k <= %ld\n",
337              *J0, *J1, *K0, *K1);
338 
339   /* rotate back to j=k=0 */
340   rotate_aux (f->coeff, b, m, k0, 0, 0);
341   rotate_aux (f->coeff, b, m, j0, 0, 1);
342 }
343 
344 /* Return the smallest value of lognorm + alpha(f + (j*x+k)*(b*x-m)) for
345    j and k small enough such that the norm does not increase too much, and
346    modify f[] accordingly.
347    The parameter "multi" means that several polynomials are wanted. If
348    multi=0 or 1, then only 1 polynomial is returned (classical behavior).
349    Otherwise, multi polynomials are stored in jmin and kmin (that
350    must be initialized arrays with at least multi elements). This option
351    might be useful for Coppersmith variant (a.k.a. MNFS).
352    In the multi case, the smallest of the returned values of lognorm + alpha
353    is returned (and f[] accordingly).
354    Warning: the caller is responsible to update the skewness if needed.
355 
356    FIXME: it would be better to pass the polynomial g to rotate, instead of
357    b and m.
358    */
359 double
rotate(mpz_poly_ptr f,unsigned long alim,mpz_t m,mpz_t b,long * jmin,long * kmin,int multi,int verbose)360 rotate (mpz_poly_ptr f, unsigned long alim, mpz_t m, mpz_t b,
361         long *jmin, long *kmin, int multi, int verbose)
362 {
363   mpz_t *D;
364   long K0, K1, J0, J1, k0, k, i, j, j0, bestk;
365   double *A, alpha, lognorm, best_alpha = DBL_MAX, best_lognorm = DBL_MAX;
366   double corr = 0.0;
367   double alpha0;
368   unsigned long p;
369   double *best_E = NULL; /* set to NULL to avoid warning... */
370   double time_alpha = 0.0, time_norm = 0.0;
371   const int d = f->deg;
372   unsigned long pp;
373   double one_over_pm1, logp, average_alpha = 0.0;
374 
375   /* The code with multi polynomials is here, but not covered at all.
376    * Remove this assert in order to test, and fix the caller which
377    * definitely wants to stick to single-polynomial stuff.
378    */
379   ASSERT_ALWAYS(multi == 0 || multi == 1);
380 
381   /* allocate best_E, to store the best (lognorm+alpha) in multi mode */
382   if (multi > 1)
383     {
384       best_E = (double *) malloc (multi * sizeof (double));
385       for (i = 0; i < multi; ++i)
386         best_E[i] = DBL_MAX;
387     }
388 
389   /* allocate D(k) = disc(f + (j*x+k)*g, x) */
390   D = (mpz_t*) malloc((d+1) * sizeof(mpz_t));
391   for(int i = 0 ; i <= d ; i++)
392       mpz_init(D[i]);
393 
394 
395   /* compute range for k */
396   rotate_bounds (f, b, m, &K0, &K1, &J0, &J1, verbose);
397   ASSERT_ALWAYS(K0 <= 0 && 0 <= K1);
398 
399   /* allocate sieving zone for computing alpha */
400   A = (double*) malloc ((K1 + 1 - K0) * sizeof (double));
401   j0 = k0 = 0; /* the current coefficients f[] correspond to f+(j*x+k)*g */
402 
403   *jmin = *kmin = 0;
404 
405   alpha0 = get_alpha (f, alim); /* value of alpha without rotation */
406 
407   ASSERT_ALWAYS(J0 < 0 && 0 < J1);
408   for (j = 0;;)
409     {
410       /* we consider j=0, 1, ..., J1, then J0, J0+1, ..., -1 */
411       j0 = rotate_aux (f->coeff, b, m, j0, j, 1);
412       /* go back to k=0 for the discriminant */
413       k0 = rotate_aux (f->coeff, b, m, k0, 0, 0);
414       /* D(k) = disc(f + (j*x+k)*g, x) (j is now fixed) */
415       discriminant_k (D, f, m, b);
416 
417       for (k = K0; k <= K1; k++)
418   A[k - K0] = 0.0; /* A[k - K0] will store the value alpha(f + k*g) */
419 
420   for (p = 2; p <= alim; p += 1 + (p & 1))
421     if (ulong_isprime (p))
422       {
423         int i;
424         /* We skip primes which divide all coefficients of f, since then
425            f mod p is zero. This can only happen when p divides N, which is
426            silly in GNFS, but maybe the user is stupid. */
427         for (i = d; i >= 0 && mpz_divisible_ui_p (f->coeff[i], p); i--);
428         if (i < 0)
429           continue;
430 
431   if (k0 != 0)
432     k0 = rotate_aux (f->coeff, b, m, k0, 0, 0);
433 
434   time_alpha -= seconds ();
435 
436   one_over_pm1 = 1.0 / (double) (p - 1);
437   logp = log ((double) p);
438   for (pp = p; pp <= alim; pp *= p)
439     {
440       /* Murphy (page 48) defines cont_p(F) = q_p*p/(p^2-1)
441          = q_p*p/(p+1)*(1/p+1/p^2+...)
442          The contribution for p^k is thus q_p*p/(p+1)/p^k. */
443       alpha = logp / (double) (p + 1) * (double) p / (double) pp;
444       /* the following is the average contribution for a prime not
445          dividing the discriminant, cf alpha.sage, function alpha_p.
446          We take it into account only for p, not for p^2, p^3, ... */
447       if (p == pp)
448         average_alpha += logp * one_over_pm1;
449       /* we do not consider the projective roots here, since the
450          corresponding correction will be considered separately for each
451          row below */
452       /* + alpha_p_projective (f, d, (D->data)[0], p); */
453       update_table (f->coeff, d, m, b, A, K0, K1, pp, alpha);
454     }
455 
456   time_alpha += seconds ();
457       } /* end of loop on primes p */
458 
459   /* determine the best alpha in each row */
460   bestk = K0;
461   for (k = K0 + 1; k <= K1; k++)
462     if (A[k - K0] < A[bestk - K0])
463       bestk = k;
464 
465   /* Correction to apply to the current row (takes into account the projective
466      roots). FIXME: we are lazy here, we should only consider the contribution
467      from the projective roots. */
468   k0 = rotate_aux (f->coeff, b, m, k0, bestk, 0);
469   corr = get_alpha (f, alim) - A[bestk - K0];
470 
471   if (verbose > 1)
472     fprintf (stderr, "# best alpha for j=%ld: k=%ld with %f\n",
473              j, bestk, A[bestk - K0] + corr);
474 
475   /* now finds the best lognorm+alpha */
476   time_norm -= seconds ();
477   for (k = K0; k <= K1; k++)
478     {
479       alpha = A[k - K0] + corr;
480       if (alpha < best_alpha + 2.0)
481         {
482           /* check lognorm + alpha < best_lognorm + best_alpha */
483 
484           /* translate from k0 to k */
485           k0 = rotate_aux (f->coeff, b, m, k0, k, 0);
486           lognorm = L2_skew_lognorm (f, SKEWNESS_DEFAULT_PREC);
487           if (multi <= 1) {
488               if (lognorm + alpha < best_lognorm + best_alpha) {
489                   best_lognorm = lognorm;
490                   best_alpha = alpha;
491                   *kmin = k;
492                   *jmin = j;
493               }
494           } else { /* multi mode */
495               /* Rem: best_lognorm + best_alpha is the worse of the
496                  preselected */
497               double newE = lognorm + alpha;
498               if (newE < best_E[multi-1]) {
499                   int ii;
500                   /* find position; assume list of preselected is sorted */
501                   for (ii = 0; ii < multi; ++ii) {
502                       if (best_E[ii] > newE)
503                           break;
504                   }
505                   ASSERT_ALWAYS(ii < multi);
506                   /* insert */
507                   for (i = multi - 1; i > ii; --i) {
508                       kmin[i] = kmin[i-1];
509                       jmin[i] = jmin[i-1];
510                       best_E[i] = best_E[i-1];
511                   }
512                   kmin[ii] = k;
513                   jmin[ii] = j;
514                   best_E[ii] = newE;
515               }
516           }
517         }
518     }
519   time_norm += seconds ();
520 
521   j++;
522   if (j > J1)
523     j = J0;
524   else if (j == 0)
525     break;
526 
527     } /* end of loop on j */
528 
529   /* we now have f + (j0*x+k0)*(bx-m) and we want f + (jmin*x+kmin)*(bx-m),
530      thus we have to add ((jmin-j0)*x+(kmin-k0)*(bx-m) */
531   /* if you are in multi-mode, we use the best polynomial */
532   rotate_aux (f->coeff, b, m, k0, *kmin, 0);
533   rotate_aux (f->coeff, b, m, j0, *jmin, 1);
534 
535   if ((verbose > 0) && (multi <= 1))
536     {
537       fprintf (stderr, "# Rotate by ");
538       if (*jmin != 0)
539         {
540           if (*jmin == -1)
541             fprintf (stderr, "-");
542           else if (*jmin != 1)
543             fprintf (stderr, "%ld*", *jmin);
544           fprintf (stderr, "x");
545           if (*kmin >= 0)
546             fprintf (stderr, "+");
547         }
548       fprintf (stderr, "%ld: alpha improved from %1.2f to %1.2f (alpha %1.2fs, norm %1.2fs)\n", *kmin, alpha0, best_alpha, time_alpha, time_norm);
549     }
550 
551   if (verbose && (multi > 1)) {
552       fprintf(stderr, "Found the following polynomials  (j, k, E):\n");
553       for (i = 0; i < multi; ++i) {
554           fprintf(stderr, "  %ld\t%ld\t%1.2f\n", jmin[i], kmin[i], best_E[i]);
555       }
556   }
557 
558   free (A);
559 
560   for(int i = 0 ; i <= d ; i++)
561       mpz_clear(D[i]);
562   free(D);
563 
564   {
565       double ret_val = best_lognorm + best_alpha;
566       if (multi>1) {
567           ret_val = best_E[0]; /* we return the smallest */
568           free(best_E);
569       }
570       return ret_val;
571   }
572 }
573 
574 
575 
576 static void
usage_and_die(char * argv0)577 usage_and_die (char *argv0)
578 {
579   fprintf (stderr, "usage: %s [-v] poly kmax\n", argv0);
580   fprintf (stderr, "  apply rotation f += (j*x+k)*g to poly.\n");
581   fprintf (stderr, "  poly: filename of polynomial\n");
582   fprintf (stderr, "  j,k : integers\n");
583   exit (1);
584 }
585 
586 int
main(int argc,char ** argv)587 main (int argc, char **argv)
588 {
589     cado_poly poly;
590     long kmax, jmin, kmin;
591     unsigned long alim = 2000;
592     mpz_t b, m;
593     int argc0 = argc, verbose = 0;
594     char **argv0 = argv;
595 
596     while (argc >= 2 && strcmp (argv[1], "-v") == 0)
597       {
598         argv ++;
599         argc --;
600         verbose ++;
601       }
602 
603     mpz_init(b);
604     mpz_init(m);
605     if (argc != 3)
606         usage_and_die (argv0[0]);
607     cado_poly_init (poly);
608     if (!cado_poly_read(poly, argv[1])) {
609         fprintf(stderr, "Problem when reading file %s\n", argv[1]);
610         usage_and_die (argv0[0]);
611     }
612     kmax = strtol(argv[2], NULL, 10);
613     MAX_k = kmax;
614 
615     poly->skew = L2_skewness (poly->pols[ALG_SIDE], SKEWNESS_DEFAULT_PREC);
616 
617     printf ("Initial polynomial:\n");
618     if (verbose)
619       print_cadopoly_extra (stdout, poly, argc0, argv0, 0);
620     else
621       printf ("skewness=%1.2f, alpha=%1.2f\n", poly->skew,
622               get_alpha (poly->pols[ALG_SIDE], get_alpha_bound ()));
623     size_optimization (poly->pols[ALG_SIDE], poly->pols[RAT_SIDE], poly->pols[ALG_SIDE], poly->pols[RAT_SIDE],
624                        SOPT_DEFAULT_EFFORT, verbose - 1);
625     poly->skew = L2_skewness (poly->pols[ALG_SIDE], SKEWNESS_DEFAULT_PREC);
626 
627     printf ("After norm optimization:\n");
628     if (verbose)
629       print_cadopoly_extra (stdout, poly, argc0, argv0, 0);
630     else
631       printf ("skewness=%1.2f, alpha=%1.2f\n",
632               poly->skew, get_alpha (poly->pols[ALG_SIDE], get_alpha_bound ()));
633 
634     mpz_set (b, poly->pols[RAT_SIDE]->coeff[1]);
635     mpz_neg (m, poly->pols[RAT_SIDE]->coeff[0]);
636     rotate (poly->pols[ALG_SIDE], alim, m, b, &jmin, &kmin, 0, verbose - 1);
637     mpz_set (poly->pols[RAT_SIDE]->coeff[1], b);
638     mpz_neg (poly->pols[RAT_SIDE]->coeff[0], m);
639     /* optimize again, but only translation */
640     mpz_poly_fprintf (stdout, poly->pols[RAT_SIDE]);
641     sopt_local_descent (poly->pols[ALG_SIDE], poly->pols[RAT_SIDE], poly->pols[ALG_SIDE], poly->pols[RAT_SIDE], 1, -1,
642                                           SOPT_DEFAULT_MAX_STEPS, verbose - 1);
643     poly->skew = L2_skewness (poly->pols[ALG_SIDE], SKEWNESS_DEFAULT_PREC);
644     mpz_clear(b);
645     mpz_clear(m);
646 
647     print_cadopoly_extra (stdout, poly, argc0, argv0, 0);
648     return 0;
649 }
650