Line data Source code
1 : /* Copyright (C) 2013 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_fflog
19 :
20 : /* Let [ be the following order on Fp: 0 [ p-1 [ 1 [ p-2 [ 2 .. [ p\2
21 : and [[ the lexicographic extension of [ to Fp[T]. Compute the
22 : isomorphism (Fp[X], [[) -> (N,<) on P */
23 :
24 : static long
25 391126 : Flx_cindex(GEN P, ulong p)
26 : {
27 391126 : long d = degpol(P), i;
28 391126 : ulong s = 0, p2 = (p-1)>>1;
29 1776942 : for (i = 0; i <= d; ++i)
30 : {
31 1385816 : ulong x = P[d-i+2];
32 1385816 : if (x<=p2) x = 2*x; else x = 1+2*(p-1-x);
33 1385816 : s = p*s+x;
34 : }
35 391126 : return s;
36 : }
37 :
38 : /* Compute the polynomial immediately after t for the [[ order */
39 :
40 : static void
41 423626 : Flx_cnext(GEN t, ulong p)
42 : {
43 : long i;
44 423626 : long p2 = p>>1;
45 540130 : for(i=2;;i++)
46 540130 : if (t[i]==p2)
47 116504 : t[i]=0;
48 : else
49 : {
50 423626 : t[i] = t[i]<p2 ? p-1-t[i]: p-t[i];
51 423626 : break;
52 : }
53 423626 : }
54 :
55 : static int
56 28 : has_deg1_auto(GEN T, ulong p, ulong pi)
57 : {
58 28 : long i, n = degpol(T);
59 28 : GEN a = polx_Flx(get_Flx_var(T));
60 672 : for (i=1; i<n; i++)
61 : {
62 644 : a = Flxq_powu_pre(a, p, T, p, pi);
63 644 : if (degpol(a)==1) return 1;
64 : }
65 28 : return 0;
66 : }
67 :
68 : static void
69 1057 : smallirred_Flx_next(GEN a, ulong p, ulong pi)
70 : {
71 : do
72 : {
73 : long i;
74 1449 : for(i=2;;i++)
75 1449 : if (++uel(a,i)==p) uel(a,i)=0;
76 1057 : else break;
77 1057 : } while (!Flx_is_irred(a, p) || has_deg1_auto(a,p,pi) );
78 28 : }
79 :
80 :
81 : /* For p = 1 [3], we do not always get full rank using the default model
82 : * this is a table of working models to use
83 : * this save time */
84 : static GEN
85 28 : smallirred_table(ulong p, ulong n)
86 : {
87 28 : switch(p)
88 : {
89 0 : case 7:
90 0 : if (n==21) return mkvecsmall3(3,5,1);
91 0 : if (n==31) return mkvecsmall3(5,1,1);
92 0 : break;
93 0 : case 13:
94 0 : if (n==19) return mkvecsmall3(8,0,1);
95 0 : if (n==23) return mkvecsmall3(2,1,1);
96 0 : if (n==37) return mkvecsmall3(3,1,1);
97 0 : break;
98 0 : case 19:
99 0 : if (n==17) return mkvecsmall2(12,2);
100 0 : if (n==29) return mkvecsmall2(7,2);
101 0 : if (n==32) return mkvecsmall3(15,1,1);
102 0 : if (n==38) return mkvecsmall3(6,1,1);
103 0 : break;
104 0 : case 31:
105 0 : if (n==21) return mkvecsmall2(6,4);
106 0 : if (n==23) return mkvecsmall3(4,1,1);
107 0 : if (n==24) return mkvecsmall4(9,2,1,1);
108 0 : if (n==26) return mkvecsmall4(18,0,0,1);
109 0 : if (n==27) return mkvecsmall3(13,3,1);
110 0 : break;
111 0 : case 37:
112 0 : if (n==22) return mkvecsmall2(15,13);
113 0 : if (n==24) return mkvecsmall2(11,4);
114 0 : if (n==29) return mkvecsmall3(20,5,1);
115 0 : break;
116 0 : case 43:
117 0 : if (n==23) return mkvecsmall3(12,1,1);
118 0 : if (n==26) return mkvecsmall4(15,0,0,1);
119 0 : if (n==27) return mkvecsmall3(1,0,25);
120 0 : break;
121 0 : case 61:
122 0 : if (n==26) return mkvecsmall3(37,1,1);
123 0 : if (n==29) return mkvecsmall3(8,2,1);
124 0 : break;
125 0 : case 67:
126 0 : if (n==26) return mkvecsmall4(60,0,0,1);
127 0 : if (n==27) return mkvecsmall2(7,7);
128 0 : break;
129 0 : case 73:
130 0 : if (n==26) return mkvecsmall3(50,1,1);
131 0 : if (n==27) return mkvecsmall3(14,0,1);
132 0 : break;
133 0 : case 79:
134 0 : if (n==23) return mkvecsmall3(36,1,1);
135 0 : if (n==26) return mkvecsmall4(22,0,0,1);
136 0 : if (n==27) return mkvecsmall3(66,0,2);
137 0 : if (n==28) return mkvecsmall3(13,1,1);
138 0 : if (n==29) return mkvecsmall4(1,0,0,1);
139 0 : break;
140 0 : case 103:
141 0 : if (n==23) return mkvecsmall2(10,2);
142 0 : break;
143 0 : case 153:
144 0 : if (n==23) return mkvecsmall2(33,1);
145 0 : break;
146 : }
147 28 : return NULL;
148 : }
149 :
150 : /* Avoid automorphisms of degree 1 */
151 : static GEN
152 28 : smallirred_Flx(ulong p, ulong n, long sv, ulong pi)
153 : {
154 28 : GEN a = zero_zv(n+2);
155 28 : GEN b = smallirred_table(p, n);
156 28 : a[1] = sv; a[3] = 1; a[n+2] = 1;
157 28 : if (b)
158 : {
159 0 : long i, lb = lg(b);
160 0 : for (i = 1 ; i < lb; i++) a[i+1] = b[i];
161 : }
162 : else
163 28 : smallirred_Flx_next(a, p, pi);
164 28 : return a;
165 : }
166 :
167 : struct Flxq_log_rel
168 : {
169 : long nbrel;
170 : GEN rel;
171 : long nb;
172 : long r, off, nbmax, nbexp;
173 : ulong nbtest;
174 : };
175 :
176 : static GEN
177 4647 : cindex_Flx(long c, long d, ulong p, long v)
178 : {
179 4647 : GEN P = cgetg(d+3, t_VECSMALL);
180 : long i;
181 4647 : P[1] = v;
182 31758 : for (i = 0; i <= d; ++i)
183 : {
184 27111 : ulong x = c%p;
185 27111 : P[i+2] = (x&1) ? p-1-(x>>1) : x>>1;
186 27111 : c/=p;
187 : }
188 4647 : return Flx_renormalize(P, d+3);
189 : }
190 :
191 : static GEN
192 10933 : factorel(GEN h, ulong p)
193 : {
194 10933 : GEN F = Flx_factor(h, p);
195 10934 : GEN F1 = gel(F, 1), F2 = gel(F, 2);
196 10934 : long i, l1 = lg(F1)-1;
197 10934 : GEN p2 = cgetg(l1+1, t_VECSMALL);
198 10934 : GEN e2 = cgetg(l1+1, t_VECSMALL);
199 51716 : for (i = 1; i <= l1; ++i)
200 : {
201 40782 : p2[i] = Flx_cindex(gel(F1, i), p);
202 40782 : e2[i] = F2[i];
203 : }
204 10934 : return mkmat2(p2, e2);
205 : }
206 :
207 : static long
208 74256 : Flx_addifsmooth3(pari_sp *av, struct Flxq_log_rel *r, GEN h, long u, long v, long w, ulong p)
209 : {
210 74256 : long off = r->off;
211 74256 : r->nbtest++;
212 74256 : if (Flx_is_smooth(h, r->r, p))
213 : {
214 5670 : GEN z = factorel(h, p);
215 5670 : if (v<0)
216 1225 : z = mkmat2(vecsmall_append(gel(z,1),off+u),vecsmall_append(gel(z,2),-1));
217 : else
218 13335 : z = famatsmall_reduce(mkmat2(
219 4445 : vecsmall_concat(gel(z,1),mkvecsmall3(off+u,off+v,off+w)),
220 4445 : vecsmall_concat(gel(z,2),mkvecsmall3(-1,-1,-1))));
221 5670 : gel(r->rel,++r->nbrel) = gc_GEN(*av,z);
222 5670 : if (DEBUGLEVEL && (r->nbrel&511UL)==0)
223 0 : err_printf("%ld%% ",r->nbrel*100/r->nbexp);
224 5670 : *av = avma;
225 68586 : } else set_avma(*av);
226 74256 : return r->nbrel==r->nb || r->nbrel==r->nbmax;
227 : }
228 :
229 : static void
230 423661 : Flx_renormalize_inplace(GEN x, long lx)
231 : {
232 : long i;
233 4695440 : for (i = lx-1; i>1; i--)
234 4692785 : if (x[i]) break;
235 423661 : setlg(x, i+1);
236 423814 : }
237 :
238 : /* Let T*X^e=C^3-R
239 : * a+b+c = 0
240 : * (C+a)*(C+b)*(C+c) = C^3+ (a*b+a*c+b*c)*C+a*b*c
241 : * = R + (a*b+a*c+b*c)*C+a*b*c
242 : * = R + (a*b-c^2)*C+a*b*c */
243 : static void
244 14 : Flxq_log_cubic(struct Flxq_log_rel *r, GEN C, GEN R, ulong p, ulong pi)
245 : {
246 14 : long l = lg(C);
247 14 : GEN a = zero_zv(l); /*We allocate one extra word to catch overflow*/
248 14 : GEN b = zero_zv(l);
249 14 : pari_sp av = avma;
250 : long i,j,k;
251 2800 : for(i=0; ; i++, Flx_cnext(a, p))
252 : {
253 2800 : Flx_renormalize_inplace(a, l+1);
254 2800 : r->nb++;
255 2800 : if (Flx_addifsmooth3(&av, r, Flx_add(a, C, p), i, -1, -1, p)) return;
256 29400 : for(j=2; j<=l; j++) b[j] = 0;
257 353129 : for(j=0; j<=i; j++, Flx_cnext(b, p))
258 : {
259 : GEN h,c;
260 : GEN pab,pabc,pabc2;
261 350343 : Flx_renormalize_inplace(b, l+1);
262 350343 : c = Flx_neg(Flx_add(a,b,p),p);
263 350343 : k = Flx_cindex(c, p);
264 350343 : if (k > j) continue;
265 71456 : pab = Flx_mul_pre(a, b, p, pi);
266 71456 : pabc = Flx_mul_pre(pab,c,p,pi);
267 71456 : pabc2= Flx_sub(pab,Flx_sqr_pre(c,p,pi),p);
268 71456 : h = Flx_add(R,Flx_add(Flx_mul_pre(C,pabc2,p,pi),pabc,p), p);
269 71456 : h = Flx_normalize(h, p);
270 71456 : if (Flx_addifsmooth3(&av, r, h, i, j, k, p)) return;
271 : }
272 : }
273 : }
274 :
275 : static GEN
276 58 : Flxq_log_find_rel(GEN b, long r, GEN T, ulong p, ulong pi, GEN *g, long *e)
277 : {
278 58 : pari_sp av = avma;
279 : while (1)
280 758 : {
281 : GEN M, z;
282 816 : *g = Flxq_mul_pre(*g, b, T, p, pi); (*e)++;
283 816 : M = Flx_halfgcd_all_pre(*g,T,p,pi,&z,NULL);
284 816 : if (Flx_is_smooth_pre(gcoeff(M,1,1), r, p, pi)
285 216 : && Flx_is_smooth_pre(z, r, p, pi))
286 : {
287 58 : GEN F = factorel(z, p);
288 58 : GEN G = factorel(gcoeff(M,1,1), p);
289 58 : GEN rel = mkmat2(vecsmall_concat(gel(F, 1),gel(G, 1)),
290 58 : vecsmall_concat(gel(F, 2),zv_neg(gel(G, 2))));
291 58 : return gc_all(av,2,&rel,g);
292 : }
293 758 : if (gc_needed(av,2))
294 : {
295 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flxq_log_find_rel");
296 0 : *g = gc_GEN(av, *g);
297 : }
298 : }
299 : }
300 :
301 : /* Generalised Odlyzko formulae ( EUROCRYPT '84, LNCS 209, pp. 224-314, 1985. ) */
302 : /* Return the number of monic, k smooth, degree n polynomials for k=1..r */
303 : static GEN
304 1036 : smoothness_vec(ulong p, long r, long n)
305 : {
306 : long i,j,k;
307 1036 : GEN R = cgetg(r+1, t_VEC), pp = utoipos(p);
308 1036 : GEN V = cgetg(n+1, t_VEC);
309 12768 : for (j = 1; j <= n; ++j)
310 11732 : gel(V, j) = binomialuu(p+j-1,j);
311 1036 : gel(R, 1) = gel(V, n);
312 3304 : for (k = 2; k <= r; ++k)
313 : {
314 2268 : GEN W = cgetg(n+1, t_VEC);
315 2268 : GEN Ik = ffnbirred(pp, k);
316 25872 : for (j = 1; j <= n; ++j)
317 : {
318 23604 : long l = j/k;
319 23604 : GEN s = gen_0;
320 23604 : pari_sp av2 = avma;
321 23604 : if (l*k == j)
322 : {
323 7812 : s = binomial(addiu(Ik,l-1), l);
324 7812 : l--;
325 : }
326 85372 : for (i = 0; i <= l; ++i)
327 61768 : s = addii(s, mulii(gel(V, j-k*i), binomial(addis(Ik,i-1), i)));
328 23604 : gel(W, j) = gc_INT(av2, s);
329 : }
330 2268 : V = W;
331 2268 : gel(R, k) = gel(V, n);
332 : }
333 1036 : return R;
334 : }
335 :
336 : /* Solve N^2*pr/6 + N*prC = N+fb
337 : N^2*pr/6 + N*(prC-1) -fb = 0
338 : */
339 :
340 : static GEN
341 868 : smooth_cost(GEN fb, GEN pr, GEN prC)
342 : {
343 868 : GEN a = gdivgu(pr,6);
344 868 : GEN b = gsubgs(prC,1);
345 868 : GEN c = gneg(fb);
346 868 : GEN vD = gsqrt(gsub(gsqr(b),gmul2n(gmul(a,c),2)),BIGDEFAULTPREC);
347 868 : return ceil_safe(gdiv(gsub(vD,b),gmul2n(a,1)));
348 : }
349 :
350 : /* Return best choice of r.
351 : We loop over d until there is sufficiently many triples (a,b,c) (a+b+c=0)
352 : of degree <=d with respect to the probability of smoothness of (a*b-c^2)*C
353 : */
354 :
355 : static GEN
356 28 : smooth_best(long p, long n, long *pt_r, long *pt_nb)
357 : {
358 28 : pari_sp av = avma, av2;
359 28 : GEN bestc = NULL, pp = utoipos(p);
360 28 : long bestr = 0, bestFB = 0;
361 28 : long r,d, dC = (n+2)/3;
362 196 : for (r = 1; r < dC; ++r)
363 : {
364 168 : GEN fb = ffsumnbirred(pp, r);
365 168 : GEN smoothC = smoothness_vec(p,r,dC);
366 168 : GEN prC = gdiv(gel(smoothC,r), powuu(p,dC));
367 168 : ulong rels = 0;
368 168 : av2 = avma;
369 938 : for(d=0; d<dC && rels < ULONG_MAX; d++)
370 : {
371 : GEN c;
372 868 : long dt = dC+2*d;
373 868 : GEN smooth = smoothness_vec(p,r,dt);
374 868 : GEN pr = gdiv(gel(smooth,r), powuu(p,dt));
375 868 : GEN FB = addii(fb,powuu(p,d));
376 868 : GEN N = smooth_cost(subiu(FB,rels),pr,prC);
377 868 : GEN Nmax = powuu(p,d+1);
378 868 : if (gcmp(N,Nmax) >= 0)
379 : {
380 770 : rels = itou_or_0(addui(rels, gceil(gmul(gdivgu(sqri(Nmax),6),pr))));
381 770 : if (!rels) rels = ULONG_MAX;
382 770 : set_avma(av2);
383 770 : continue;
384 : }
385 98 : c = gdivgu(addii(powuu(p,2*d),sqri(N)),6);
386 98 : FB = addii(FB,N);
387 98 : if ((!bestc || gcmp(gmul2n(c,r), gmul2n(bestc,bestr)) < 0))
388 : {
389 42 : if (DEBUGLEVEL)
390 0 : err_printf("r=%ld d=%ld fb=%Ps early rels=%lu P=%.5Pe -> C=%.5Pe \n",
391 : r, dt, FB, rels, pr, c);
392 42 : bestc = c;
393 42 : bestr = r;
394 42 : bestFB = itos_or_0(FB);
395 : }
396 98 : break;
397 : }
398 : }
399 28 : *pt_r=bestr;
400 28 : *pt_nb=bestFB;
401 28 : return bestc ? gc_upto(av, gceil(bestc)): NULL;
402 : }
403 :
404 : static GEN
405 28 : check_kernel(long r, GEN M, long nbi, long nbrow, GEN T, ulong p, ulong pi, GEN m)
406 : {
407 28 : pari_sp av = avma;
408 28 : long N = 3*upowuu(p, r);
409 28 : GEN K = FpMs_leftkernel_elt(M, nbrow, m);
410 28 : long i, f=0, tbs;
411 28 : long lm = lgefint(m), u=1;
412 : GEN tab, g;
413 28 : GEN q = powuu(p,degpol(T));
414 28 : GEN idx = diviiexact(subiu(q,1),m);
415 : pari_timer ti;
416 28 : if (DEBUGLEVEL) timer_start(&ti);
417 224 : while (signe(gel(K,u))==0)
418 196 : u++;
419 28 : K = FpC_Fp_mul(K, Fp_inv(gel(K, u), m), m);
420 28 : g = Flxq_pow_pre(cindex_Flx(u, r, p, T[1]), idx, T, p, pi);
421 28 : tbs = maxss(1, expu(nbi/expi(m)));
422 28 : tab = Flxq_pow_init_pre(g, q, tbs, T, p, pi);
423 28 : setlg(K, N);
424 46662 : for (i=1; i<N; i++)
425 : {
426 46634 : GEN k = gel(K,i);
427 46634 : pari_sp av = avma;
428 51027 : long t = signe(k) && Flx_equal(Flxq_pow_table_pre(tab, k, T, p, pi),
429 4393 : Flxq_pow_pre(cindex_Flx(i,r,p,T[1]), idx, T, p, pi));
430 46634 : set_avma(av);
431 46634 : if (!t)
432 42241 : gel(K,i) = cgetineg(lm);
433 : else
434 4393 : f++;
435 : }
436 28 : if (DEBUGLEVEL) timer_printf(&ti,"found %ld/%ld logs", f, nbi);
437 28 : if (f < maxss(3,maxss(p/2,nbi/p))) return NULL; /* Not enough logs found */
438 28 : return gc_GEN(av, K);
439 : }
440 :
441 : static GEN
442 28 : Flxq_log_rec(GEN W, GEN a, long r, GEN T, ulong p, ulong pi, GEN m)
443 : {
444 28 : long AV = 0, u = 1;
445 28 : GEN g = a, b;
446 : pari_timer ti;
447 280 : while (!equali1(gel(W,u)))
448 252 : u++;
449 28 : b = cindex_Flx(u, r, p, T[1]);
450 : while(1)
451 1 : {
452 : long i, l;
453 : GEN V, F, E, Ao;
454 29 : timer_start(&ti);
455 29 : V = Flxq_log_find_rel(b, r, T, p, pi, &g, &AV);
456 29 : if (DEBUGLEVEL>1) timer_printf(&ti,"%ld-smooth element",r);
457 29 : F = gel(V,1); E = gel(V,2);
458 29 : l = lg(F);
459 29 : Ao = gen_0;
460 256 : for(i=1; i<l; i++)
461 : {
462 228 : GEN R = gel(W,F[i]);
463 228 : if (signe(R)<=0)
464 1 : break;
465 227 : Ao = Fp_add(Ao, mulis(R, E[i]), m);
466 : }
467 29 : if (i==l) return subis(Ao,AV);
468 : }
469 : }
470 :
471 : static int
472 308 : Flxq_log_use_index_cubic(GEN m, GEN T0, ulong p)
473 : {
474 308 : pari_sp av = avma;
475 308 : long n = get_Flx_degree(T0), r, nb, e = expi(m);
476 308 : if (e >= 27 && 10*(n*expu(p)+6)<=e*e)
477 : {
478 14 : GEN cost_rho = sqrti(shifti(m,2));
479 14 : GEN cost = smooth_best(p, n, &r, &nb);
480 14 : int use = (cost && gcmp(cost,cost_rho)<0);
481 14 : set_avma(av);
482 14 : return use;
483 294 : } else return 0;
484 : }
485 :
486 : static GEN
487 14 : Flxq_log_index_cubic(GEN a0, GEN b0, GEN m, GEN T0, ulong p)
488 : {
489 14 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
490 14 : long n = get_Flx_degree(T0), r, nb;
491 14 : pari_sp av = avma;
492 : struct Flxq_log_rel rel;
493 : long nbi;
494 : GEN W, M, S, T, a, b, Ao, Bo, e, C, R;
495 : pari_timer ti;
496 14 : GEN cost = smooth_best(p, n, &r, &nb);
497 14 : GEN cost_rho = sqrti(shifti(m,2));
498 14 : if (!cost || gcmp(cost,cost_rho)>=0) return gc_NULL(av);
499 14 : nbi = itos(ffsumnbirred(stoi(p), r));
500 14 : if (DEBUGLEVEL)
501 : {
502 0 : err_printf("Size FB=%ld, looking for %ld relations, %Ps tests needed\n", nbi, nb,cost);
503 0 : timer_start(&ti);
504 : }
505 14 : T = smallirred_Flx(p,n,get_Flx_var(T0), pi);
506 : for(;;)
507 : {
508 14 : S = Flx_ffisom(T0,T,p);
509 14 : a = Flx_Flxq_eval_pre(a0, S, T, p, pi);
510 14 : b = Flx_Flxq_eval_pre(b0, S, T, p, pi);
511 14 : C = Flx_shift(pol1_Flx(get_Flx_var(T)), (n+2)/3);
512 14 : R = Flxq_powu_pre(C,3,T,p,pi);
513 14 : if (DEBUGLEVEL)
514 0 : timer_printf(&ti," model change: %Ps",Flx_to_ZX(T));
515 14 : rel.nbmax=2*nb;
516 14 : M = cgetg(rel.nbmax+1, t_VEC);
517 14 : rel.rel = M;
518 14 : rel.nbrel = 0; rel.r = r; rel.off = 3*upowuu(p,r);
519 14 : rel.nb = nbi; rel.nbexp = nb; rel.nbtest=0;
520 14 : Flxq_log_cubic(&rel, C, R, p, pi);
521 14 : setlg(M,1+rel.nbrel);
522 14 : if (DEBUGLEVEL)
523 : {
524 0 : err_printf("\n");
525 0 : timer_printf(&ti," %ld relations, %ld generators (%ld tests)",rel.nbrel,rel.nb,rel.nbtest);
526 : }
527 14 : W = check_kernel(r, M, nbi, rel.off + rel.nb - nbi, T, p, pi, m);
528 14 : if (W) break;
529 0 : if (DEBUGLEVEL) timer_start(&ti);
530 0 : smallirred_Flx_next(T,p, pi);
531 : }
532 14 : if (DEBUGLEVEL) timer_start(&ti);
533 14 : Ao = Flxq_log_rec(W, a, r, T, p, pi, m);
534 14 : if (DEBUGLEVEL) timer_printf(&ti,"smooth element");
535 14 : Bo = Flxq_log_rec(W, b, r, T, p, pi, m);
536 14 : if (DEBUGLEVEL) timer_printf(&ti,"smooth generator");
537 14 : e = Fp_div(Ao, Bo, m);
538 14 : if (!Flx_equal(Flxq_pow_pre(b0, e, T0, p, pi), a0)) pari_err_BUG("Flxq_log");
539 14 : return gc_upto(av, e);
540 : }
541 :
542 21876 : INLINE GEN Flx_frob(GEN u, ulong p) { return Flx_inflate(u, p); }
543 :
544 : static GEN
545 32291 : rel_Coppersmith(long r, GEN u, GEN v, long h, GEN R, long d, ulong p, ulong pi)
546 : {
547 : GEN a, b, F, G, M;
548 32291 : if (degpol(Flx_gcd_pre(u,v,p,pi))) return NULL;
549 32267 : a = Flx_add(Flx_shift(u, h), v, p);
550 32377 : if (lgpol(a)==0 || !Flx_is_smooth_pre(a, r, p, pi)) return NULL;
551 8360 : b = Flx_add(Flx_mul_pre(R, Flx_frob(u, p), p, pi),
552 : Flx_shift(Flx_frob(v, p),d), p);
553 8415 : if (!Flx_is_smooth_pre(b, r, p, pi)) return NULL;
554 2523 : F = factorel(a, p); G = factorel(b, p);
555 5048 : M = mkmat2(vecsmall_concat(gel(F, 1), vecsmall_append(gel(G, 1), 2*p)),
556 5048 : vecsmall_concat(zv_z_mul(gel(F, 2),p), vecsmall_append(zv_neg(gel(G, 2)),d)));
557 2524 : return famatsmall_reduce(M);
558 : }
559 :
560 : GEN
561 1355 : Flxq_log_Coppersmith_worker(GEN u, long i, GEN V, GEN R)
562 : {
563 1355 : long r = V[1], h = V[2], d = V[3], p = V[4], pi = V[5], dT = V[6];
564 1355 : pari_sp ltop = avma;
565 1355 : GEN v = zero_zv(dT+2);
566 1355 : GEN L = cgetg(2*i+1, t_VEC);
567 1354 : pari_sp av = avma;
568 : long j;
569 1354 : long nbtest=0, rel = 1;
570 1354 : ulong lu = Flx_lead(u), lv;
571 70482 : for (j=1; j<=i; j++)
572 : {
573 : GEN z;
574 69127 : Flx_cnext(v, p);
575 69136 : Flx_renormalize_inplace(v, dT+2);
576 69278 : lv = Flx_lead(v);
577 69263 : set_avma(av);
578 69256 : if (lu != 1 && lv != 1) continue;
579 39994 : if (degpol(Flx_gcd_pre(u, v, p, pi))!=0) continue;
580 26915 : if (lu==1)
581 : {
582 14936 : z = rel_Coppersmith(r, u, v, h, R, d, p, pi);
583 14956 : nbtest++;
584 14956 : if (z) { gel(L, rel++) = z; av = avma; }
585 : }
586 26935 : if (i==j) continue;
587 26850 : if (lv==1)
588 : {
589 17359 : z = rel_Coppersmith(r, v, u, h, R, d, p, pi);
590 17406 : nbtest++;
591 17406 : if (z) { gel(L, rel++) = z; av = avma; }
592 : }
593 : }
594 1355 : setlg(L,rel);
595 1355 : return gc_GEN(ltop, mkvec2(stoi(nbtest), L));
596 : }
597 :
598 : static GEN
599 14 : Flxq_log_Coppersmith(long nbrel, long r, GEN T, ulong p, ulong pi)
600 : {
601 : pari_sp av;
602 14 : long dT = degpol(T);
603 14 : long h = dT/p, d = dT-(h*p);
604 14 : GEN R = Flx_sub(Flx_shift(pol1_Flx(T[1]), dT), T, p);
605 14 : GEN u = zero_zv(dT+2);
606 : GEN done;
607 14 : long nbtest = 0, rel = 0;
608 14 : GEN M = cgetg(nbrel+1, t_VEC);
609 14 : long i = 1;
610 14 : GEN worker = snm_closure(is_entry("_Flxq_log_Coppersmith_worker"),
611 : mkvec2(mkvecsmalln(6, r,h,d,p,pi,dT), R));
612 : struct pari_mt pt;
613 14 : long running, pending = 0, stop=0;
614 14 : if (DEBUGLEVEL) err_printf("Coppersmith (R = %ld): ",degpol(R));
615 14 : mt_queue_start(&pt, worker);
616 14 : av = avma;
617 1399 : while ((running = !stop) || pending)
618 : {
619 : GEN L;
620 : long l, j;
621 1385 : Flx_cnext(u, p);
622 1385 : Flx_renormalize_inplace(u, dT+2);
623 1385 : mt_queue_submit(&pt, 0, running ? mkvec2(u, stoi(i)): NULL);
624 1385 : done = mt_queue_get(&pt, NULL, &pending);
625 1385 : if (!done) continue;
626 1355 : L = gel(done, 2); nbtest += itos(gel(done,1));
627 1355 : l = lg(L);
628 1355 : if (l > 1)
629 : {
630 3392 : for (j=1; j<l; j++)
631 : {
632 2466 : if (rel>nbrel) break;
633 2429 : gel(M,++rel) = gel(L,j);
634 2429 : if (DEBUGLEVEL && (rel&511UL)==0)
635 0 : err_printf("%ld%%[%ld] ",rel*100/nbrel,i);
636 : }
637 963 : av = avma;
638 : }
639 392 : else set_avma(av);
640 1355 : if (rel>nbrel) stop = 1;
641 1355 : i++;
642 : }
643 14 : mt_queue_end(&pt);
644 14 : if (DEBUGLEVEL) err_printf(": %ld tests\n", nbtest);
645 14 : return M;
646 : }
647 :
648 : static GEN Flxq_log_Coppersmith_d(GEN W, GEN g, long r, GEN T, ulong p, ulong pi, GEN mo);
649 :
650 : static GEN
651 50 : Flxq_log_from_rel(GEN W, GEN rel, long r, GEN T, ulong p, ulong pi, GEN m)
652 : {
653 50 : pari_sp av = avma;
654 50 : GEN F = gel(rel,1), E = gel(rel,2), o = gen_0;
655 50 : long i, l = lg(F);
656 356 : for(i=1; i<l; i++)
657 : {
658 306 : GEN R = gel(W, F[i]);
659 306 : if (signe(R)==0) /* Already failed */
660 0 : return NULL;
661 306 : else if (signe(R)<0) /* Not yet tested */
662 : {
663 1 : setsigne(gel(W,F[i]),0);
664 1 : R = Flxq_log_Coppersmith_d(W, cindex_Flx(F[i],r,p,T[1]), r, T, p, pi, m);
665 1 : if (!R) return NULL;
666 : }
667 306 : o = Fp_add(o, mulis(R, E[i]), m);
668 : }
669 50 : return gc_INT(av, o);
670 : }
671 :
672 : static GEN
673 51 : Flxq_log_Coppersmith_d(GEN W, GEN g, long r, GEN T, ulong p, ulong pi, GEN mo)
674 : {
675 51 : pari_sp av = avma, av2;
676 51 : long dg = degpol(g), k = r-1, m = maxss((dg-k)/2,0);
677 51 : long i, j, l = dg-m, N;
678 51 : GEN v = cgetg(k+m+1,t_MAT);
679 51 : long dT = degpol(T);
680 51 : long h = dT/p, d = dT-h*p;
681 51 : GEN R = Flx_rem_pre(Flx_shift(pol1_Flx(T[1]), dT), T, p, pi);
682 51 : GEN z = Flx_rem_pre(Flx_shift(pol1_Flx(T[1]), h), g, p, pi);
683 350 : for(i=1; i<=k+m; i++)
684 : {
685 299 : gel(v,i) = Flx_to_Flv(Flx_shift(z,-l),m);
686 299 : z = Flx_rem_pre(Flx_shift(z,1),g,p,pi);
687 : }
688 51 : v = Flm_ker(v,p);
689 286 : for(i=1; i<=k; i++)
690 235 : gel(v,i) = Flv_to_Flx(gel(v,i),T[1]);
691 51 : N = upowuu(p,k);
692 51 : av2 = avma;
693 2538 : for (i=1; i<N; i++)
694 : {
695 : GEN p0,q,qh,a,b;
696 2537 : ulong el = i;
697 2537 : set_avma(av2);
698 2537 : q = pol0_Flx(T[1]);
699 14713 : for (j=1; j<=k; j++)
700 : {
701 12176 : ulong r = el % p;
702 12176 : el /= p;
703 12176 : if (r) q = Flx_add(q, Flx_Fl_mul(gel(v,j), r, p), p);
704 : }
705 2537 : qh = Flx_shift(q, h);
706 2537 : p0 = Flx_rem_pre(qh, g, p, pi);
707 2537 : b = Flx_sub(Flx_mul_pre(R, Flx_frob(q, p), p, pi),
708 : Flx_shift(Flx_frob(p0, p), d), p);
709 2537 : if (lgpol(b)==0 || !Flx_is_smooth_pre(b, r, p, pi)) continue;
710 59 : a = Flx_div_pre(Flx_sub(qh, p0, p), g, p, pi);
711 59 : if (degpol(Flx_gcd_pre(a, q, p, pi)) && degpol(Flx_gcd_pre(a, p0, p, pi)))
712 0 : continue;
713 59 : if (!(lgpol(a)==0 || !Flx_is_smooth_pre(a, r, p, pi)))
714 : {
715 50 : GEN F = factorel(b, p);
716 50 : GEN G = factorel(a, p);
717 50 : GEN FG = vecsmall_concat(vecsmall_append(gel(F, 1), 2*p), gel(G, 1));
718 50 : GEN E = vecsmall_concat(vecsmall_append(gel(F, 2), -d),
719 50 : zv_z_mul(gel(G, 2),-p));
720 50 : GEN R = famatsmall_reduce(mkmat2(FG, E));
721 50 : GEN l = Flxq_log_from_rel(W, R, r, T, p, pi, mo);
722 50 : if (!l) continue;
723 50 : l = Fp_divu(l,p,mo);
724 50 : if (dg <= r)
725 : {
726 1 : long idx = Flx_cindex(g, p);
727 1 : affii(l, gel(W, idx));
728 1 : if (DEBUGLEVEL>1) err_printf("Found %lu\n", idx);
729 : }
730 50 : return gc_INT(av, l);
731 : }
732 : }
733 1 : set_avma(av);
734 1 : return NULL;
735 : }
736 :
737 : static GEN
738 28 : Flxq_log_Coppersmith_rec(GEN W, long r2, GEN a, long r, GEN T, ulong p, ulong pi, GEN m)
739 : {
740 28 : GEN b = polx_Flx(T[1]);
741 28 : long AV = 0;
742 28 : GEN g = a, bad = pol0_Flx(T[1]);
743 : pari_timer ti;
744 : while(1)
745 1 : {
746 : long i, l;
747 : GEN V, F, E, Ao;
748 29 : timer_start(&ti);
749 29 : V = Flxq_log_find_rel(b, r2, T, p, pi, &g, &AV);
750 29 : if (DEBUGLEVEL>1) timer_printf(&ti,"%ld-smooth element",r2);
751 29 : F = gel(V,1); E = gel(V,2);
752 29 : l = lg(F);
753 29 : Ao = gen_0;
754 225 : for(i=1; i<l; i++)
755 : {
756 197 : GEN Fi = cindex_Flx(F[i], r2, p, T[1]);
757 : GEN R;
758 197 : if (degpol(Fi) <= r)
759 : {
760 147 : if (signe(gel(W,F[i]))==0)
761 0 : break;
762 147 : else if (signe(gel(W,F[i]))<0)
763 : {
764 0 : setsigne(gel(W,F[i]),0);
765 0 : R = Flxq_log_Coppersmith_d(W,Fi,r,T,p,pi,m);
766 : } else
767 147 : R = gel(W,F[i]);
768 : }
769 : else
770 : {
771 50 : if (Flx_equal(Fi,bad)) break;
772 50 : R = Flxq_log_Coppersmith_d(W,Fi,r,T,p,pi,m);
773 50 : if (!R) bad = Fi;
774 : }
775 197 : if (!R) break;
776 196 : Ao = Fp_add(Ao, mulis(R, E[i]), m);
777 : }
778 29 : if (i==l) return subis(Ao,AV);
779 : }
780 : }
781 :
782 : static GEN
783 14 : Flxq_log_index_Coppersmith(GEN a0, GEN b0, GEN m, GEN T0, ulong p)
784 : {
785 14 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
786 14 : pari_sp av = avma;
787 14 : GEN M, S, a, b, Ao=NULL, Bo=NULL, W, e;
788 : pari_timer ti;
789 14 : double rf = p ==3 ? 1.2 : .9;
790 14 : long n = degpol(T0), r = (long) sqrt(n*rf);
791 : GEN T;
792 14 : long r2 = 3*r/2;
793 14 : long nbi = itos(ffsumnbirred(utoipos(p), r)), nbrel=nbi*5/4;
794 14 : if (DEBUGLEVEL)
795 : {
796 0 : err_printf("Coppersmith: Parameters r=%ld r2=%ld\n", r,r2);
797 0 : err_printf("Coppersmith: Size FB=%ld rel. needed=%ld\n", nbi, nbrel);
798 0 : timer_start(&ti);
799 : }
800 14 : T = smallirred_Flx(p,n,get_Flx_var(T0), pi);
801 14 : S = Flx_ffisom(T0,T,p);
802 14 : a = Flx_Flxq_eval_pre(a0, S, T, p, pi);
803 14 : b = Flx_Flxq_eval_pre(b0, S, T, p, pi);
804 14 : if (DEBUGLEVEL) timer_printf(&ti,"model change");
805 14 : M = Flxq_log_Coppersmith(nbrel, r, T, p, pi);
806 14 : if (DEBUGLEVEL) timer_printf(&ti,"relations");
807 14 : W = check_kernel(r, M, nbi, 3*upowuu(p,r), T, p, pi, m);
808 14 : timer_start(&ti);
809 14 : Ao = Flxq_log_Coppersmith_rec(W, r2, a, r, T, p, pi, m);
810 14 : if (DEBUGLEVEL) timer_printf(&ti,"smooth element");
811 14 : Bo = Flxq_log_Coppersmith_rec(W, r2, b, r, T, p, pi, m);
812 14 : if (DEBUGLEVEL) timer_printf(&ti,"smooth generator");
813 14 : e = Fp_div(Ao, Bo, m);
814 14 : if (!Flx_equal(Flxq_pow_pre(b0,e,T0,p,pi), a0)) pari_err_BUG("Flxq_log");
815 14 : return gc_upto(av, e);
816 : }
817 :
818 : GEN
819 28 : Flxq_log_index(GEN a, GEN b, GEN m, GEN T, ulong p)
820 : {
821 28 : long d = get_Flx_degree(T);
822 28 : if (p==3 || (p==5 && d>41))
823 14 : return Flxq_log_index_Coppersmith(a, b, m, T, p);
824 14 : else return Flxq_log_index_cubic(a, b, m, T, p);
825 : }
826 :
827 : int
828 164210 : Flxq_log_use_index(GEN m, GEN T, ulong p)
829 : {
830 164210 : long d = get_Flx_degree(T);
831 164210 : if (p==3 || (p==5 && d>41))
832 24269 : return 1;
833 139941 : else if (d<=4 || d==6)
834 139633 : return 0;
835 : else
836 308 : return Flxq_log_use_index_cubic(m, T, p);
837 : }
|
