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