1 /* Copyright (C) 2000-2004 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 /***********************************************************************/
16 /** **/
17 /** GENERIC ALGORITHMS ON BLACKBOX GROUP **/
18 /** **/
19 /***********************************************************************/
20 #include "pari.h"
21 #include "paripriv.h"
22 #undef pow /* AIX: pow(a,b) is a macro, wrongly expanded on grp->pow(a,b,c) */
23
24 /***********************************************************************/
25 /** **/
26 /** POWERING **/
27 /** **/
28 /***********************************************************************/
29
30 /* return (n>>(i+1-l)) & ((1<<l)-1) */
31 static ulong
int_block(GEN n,long i,long l)32 int_block(GEN n, long i, long l)
33 {
34 long q = divsBIL(i), r = remsBIL(i)+1, lr;
35 GEN nw = int_W(n, q);
36 ulong w = (ulong) *nw, w2;
37 if (r>=l) return (w>>(r-l))&((1UL<<l)-1);
38 w &= (1UL<<r)-1; lr = l-r;
39 w2 = (ulong) *int_precW(nw); w2 >>= (BITS_IN_LONG-lr);
40 return (w<<lr)|w2;
41 }
42
43 /* assume n != 0, t_INT. Compute x^|n| using sliding window powering */
44 static GEN
sliding_window_powu(GEN x,ulong n,long e,void * E,GEN (* sqr)(void *,GEN),GEN (* mul)(void *,GEN,GEN))45 sliding_window_powu(GEN x, ulong n, long e, void *E, GEN (*sqr)(void*,GEN),
46 GEN (*mul)(void*,GEN,GEN))
47 {
48 pari_sp av;
49 long i, l = expu(n), u = (1UL<<(e-1));
50 long w, v;
51 GEN tab = cgetg(1+u, t_VEC);
52 GEN x2 = sqr(E, x), z = NULL, tw;
53 gel(tab, 1) = x;
54 for (i=2; i<=u; i++) gel(tab,i) = mul(E, gel(tab,i-1), x2);
55 av = avma;
56 while (l>=0)
57 {
58 if (e > l+1) e = l+1;
59 w = (n>>(l+1-e)) & ((1UL<<e)-1); v = vals(w); l-=e;
60 tw = gel(tab, 1+(w>>(v+1)));
61 if (z)
62 {
63 for (i=1; i<=e-v; i++) z = sqr(E, z);
64 z = mul(E, z, tw);
65 } else z = tw;
66 for (i=1; i<=v; i++) z = sqr(E, z);
67 while (l>=0)
68 {
69 if (gc_needed(av,1))
70 {
71 if (DEBUGMEM>1) pari_warn(warnmem,"sliding_window_powu (%ld)", l);
72 z = gerepilecopy(av, z);
73 }
74 if (n&(1UL<<l)) break;
75 z = sqr(E, z); l--;
76 }
77 }
78 return z;
79 }
80
81 /* assume n != 0, t_INT. Compute x^|n| using sliding window powering */
82 static GEN
sliding_window_pow(GEN x,GEN n,long e,void * E,GEN (* sqr)(void *,GEN),GEN (* mul)(void *,GEN,GEN))83 sliding_window_pow(GEN x, GEN n, long e, void *E, GEN (*sqr)(void*,GEN),
84 GEN (*mul)(void*,GEN,GEN))
85 {
86 pari_sp av;
87 long i, l = expi(n), u = (1UL<<(e-1));
88 long w, v;
89 GEN tab = cgetg(1+u, t_VEC);
90 GEN x2 = sqr(E, x), z = NULL, tw;
91 gel(tab, 1) = x;
92 for (i=2; i<=u; i++) gel(tab,i) = mul(E, gel(tab,i-1), x2);
93 av = avma;
94 while (l>=0)
95 {
96 if (e > l+1) e = l+1;
97 w = int_block(n,l,e); v = vals(w); l-=e;
98 tw = gel(tab, 1+(w>>(v+1)));
99 if (z)
100 {
101 for (i=1; i<=e-v; i++) z = sqr(E, z);
102 z = mul(E, z, tw);
103 } else z = tw;
104 for (i=1; i<=v; i++) z = sqr(E, z);
105 while (l>=0)
106 {
107 if (gc_needed(av,1))
108 {
109 if (DEBUGMEM>1) pari_warn(warnmem,"sliding_window_pow (%ld)", l);
110 z = gerepilecopy(av, z);
111 }
112 if (int_bit(n,l)) break;
113 z = sqr(E, z); l--;
114 }
115 }
116 return z;
117 }
118
119 /* assume n != 0, t_INT. Compute x^|n| using leftright binary powering */
120 static GEN
leftright_binary_powu(GEN x,ulong n,void * E,GEN (* sqr)(void *,GEN),GEN (* mul)(void *,GEN,GEN))121 leftright_binary_powu(GEN x, ulong n, void *E, GEN (*sqr)(void*,GEN),
122 GEN (*mul)(void*,GEN,GEN))
123 {
124 pari_sp av = avma;
125 GEN y;
126 int j;
127
128 if (n == 1) return x;
129 y = x; j = 1+bfffo(n);
130 /* normalize, i.e set highest bit to 1 (we know n != 0) */
131 n<<=j; j = BITS_IN_LONG-j;
132 /* first bit is now implicit */
133 for (; j; n<<=1,j--)
134 {
135 y = sqr(E,y);
136 if (n & HIGHBIT) y = mul(E,y,x); /* first bit set: multiply by base */
137 if (gc_needed(av,1))
138 {
139 if (DEBUGMEM>1) pari_warn(warnmem,"leftright_powu (%d)", j);
140 y = gerepilecopy(av, y);
141 }
142 }
143 return y;
144 }
145
146 GEN
gen_powu_i(GEN x,ulong n,void * E,GEN (* sqr)(void *,GEN),GEN (* mul)(void *,GEN,GEN))147 gen_powu_i(GEN x, ulong n, void *E, GEN (*sqr)(void*,GEN),
148 GEN (*mul)(void*,GEN,GEN))
149 {
150 long l;
151 if (n == 1) return x;
152 l = expu(n);
153 if (l<=8)
154 return leftright_binary_powu(x, n, E, sqr, mul);
155 else
156 return sliding_window_powu(x, n, l<=24? 2: 3, E, sqr, mul);
157 }
158
159 GEN
gen_powu(GEN x,ulong n,void * E,GEN (* sqr)(void *,GEN),GEN (* mul)(void *,GEN,GEN))160 gen_powu(GEN x, ulong n, void *E, GEN (*sqr)(void*,GEN),
161 GEN (*mul)(void*,GEN,GEN))
162 {
163 pari_sp av = avma;
164 if (n == 1) return gcopy(x);
165 return gerepilecopy(av, gen_powu_i(x,n,E,sqr,mul));
166 }
167
168 GEN
gen_pow_i(GEN x,GEN n,void * E,GEN (* sqr)(void *,GEN),GEN (* mul)(void *,GEN,GEN))169 gen_pow_i(GEN x, GEN n, void *E, GEN (*sqr)(void*,GEN),
170 GEN (*mul)(void*,GEN,GEN))
171 {
172 long l, e;
173 if (lgefint(n)==3) return gen_powu_i(x, uel(n,2), E, sqr, mul);
174 l = expi(n);
175 if (l<=64) e = 3;
176 else if (l<=160) e = 4;
177 else if (l<=384) e = 5;
178 else if (l<=896) e = 6;
179 else e = 7;
180 return sliding_window_pow(x, n, e, E, sqr, mul);
181 }
182
183 GEN
gen_pow(GEN x,GEN n,void * E,GEN (* sqr)(void *,GEN),GEN (* mul)(void *,GEN,GEN))184 gen_pow(GEN x, GEN n, void *E, GEN (*sqr)(void*,GEN),
185 GEN (*mul)(void*,GEN,GEN))
186 {
187 pari_sp av = avma;
188 return gerepilecopy(av, gen_pow_i(x,n,E,sqr,mul));
189 }
190
191 /* assume n > 0. Compute x^n using left-right binary powering */
192 GEN
gen_powu_fold_i(GEN x,ulong n,void * E,GEN (* sqr)(void *,GEN),GEN (* msqr)(void *,GEN))193 gen_powu_fold_i(GEN x, ulong n, void *E, GEN (*sqr)(void*,GEN),
194 GEN (*msqr)(void*,GEN))
195 {
196 pari_sp av = avma;
197 GEN y;
198 int j;
199
200 if (n == 1) return x;
201 y = x; j = 1+bfffo(n);
202 /* normalize, i.e set highest bit to 1 (we know n != 0) */
203 n<<=j; j = BITS_IN_LONG-j;
204 /* first bit is now implicit */
205 for (; j; n<<=1,j--)
206 {
207 if (n & HIGHBIT) y = msqr(E,y); /* first bit set: multiply by base */
208 else y = sqr(E,y);
209 if (gc_needed(av,1))
210 {
211 if (DEBUGMEM>1) pari_warn(warnmem,"gen_powu_fold (%d)", j);
212 y = gerepilecopy(av, y);
213 }
214 }
215 return y;
216 }
217 GEN
gen_powu_fold(GEN x,ulong n,void * E,GEN (* sqr)(void *,GEN),GEN (* msqr)(void *,GEN))218 gen_powu_fold(GEN x, ulong n, void *E, GEN (*sqr)(void*,GEN),
219 GEN (*msqr)(void*,GEN))
220 {
221 pari_sp av = avma;
222 if (n == 1) return gcopy(x);
223 return gerepilecopy(av, gen_powu_fold_i(x,n,E,sqr,msqr));
224 }
225
226 /* assume N != 0, t_INT. Compute x^|N| using left-right binary powering */
227 GEN
gen_pow_fold_i(GEN x,GEN N,void * E,GEN (* sqr)(void *,GEN),GEN (* msqr)(void *,GEN))228 gen_pow_fold_i(GEN x, GEN N, void *E, GEN (*sqr)(void*,GEN),
229 GEN (*msqr)(void*,GEN))
230 {
231 long ln = lgefint(N);
232 if (ln == 3) return gen_powu_fold_i(x, N[2], E, sqr, msqr);
233 else
234 {
235 GEN nd = int_MSW(N), y = x;
236 ulong n = *nd;
237 long i;
238 int j;
239 pari_sp av = avma;
240
241 if (n == 1)
242 j = 0;
243 else
244 {
245 j = 1+bfffo(n); /* < BIL */
246 /* normalize, i.e set highest bit to 1 (we know n != 0) */
247 n <<= j; j = BITS_IN_LONG - j;
248 }
249 /* first bit is now implicit */
250 for (i=ln-2;;)
251 {
252 for (; j; n<<=1,j--)
253 {
254 if (n & HIGHBIT) y = msqr(E,y); /* first bit set: multiply by base */
255 else y = sqr(E,y);
256 if (gc_needed(av,1))
257 {
258 if (DEBUGMEM>1) pari_warn(warnmem,"gen_pow_fold (%d)", j);
259 y = gerepilecopy(av, y);
260 }
261 }
262 if (--i == 0) return y;
263 nd = int_precW(nd);
264 n = *nd; j = BITS_IN_LONG;
265 }
266 }
267 }
268 GEN
gen_pow_fold(GEN x,GEN n,void * E,GEN (* sqr)(void *,GEN),GEN (* msqr)(void *,GEN))269 gen_pow_fold(GEN x, GEN n, void *E, GEN (*sqr)(void*,GEN),
270 GEN (*msqr)(void*,GEN))
271 {
272 pari_sp av = avma;
273 return gerepilecopy(av, gen_pow_fold_i(x,n,E,sqr,msqr));
274 }
275
276 GEN
gen_pow_init(GEN x,GEN n,long k,void * E,GEN (* sqr)(void *,GEN),GEN (* mul)(void *,GEN,GEN))277 gen_pow_init(GEN x, GEN n, long k, void *E, GEN (*sqr)(void*,GEN), GEN (*mul)(void*,GEN,GEN))
278 {
279 long i, j, l = expi(n);
280 long m = 1UL<<(k-1);
281 GEN x2 = sqr(E, x), y = gcopy(x);
282 GEN R = cgetg(m+1, t_VEC);
283 for(i = 1; i <= m; i++)
284 {
285 GEN C = cgetg(l+1, t_VEC);
286 gel(C,1) = y;
287 for(j = 2; j <= l; j++)
288 gel(C,j) = sqr(E, gel(C,j-1));
289 gel(R,i) = C;
290 y = mul(E, y, x2);
291 }
292 return R;
293 }
294
295 GEN
gen_pow_table(GEN R,GEN n,void * E,GEN (* one)(void *),GEN (* mul)(void *,GEN,GEN))296 gen_pow_table(GEN R, GEN n, void *E, GEN (*one)(void*), GEN (*mul)(void*,GEN,GEN))
297 {
298 long e = expu(lg(R)-1) + 1;
299 long l = expi(n);
300 long i, w;
301 GEN z = one(E), tw;
302 for(i=0; i<=l; )
303 {
304 if (int_bit(n, i)==0) { i++; continue; }
305 if (i+e-1>l) e = l+1-i;
306 w = int_block(n,i+e-1,e);
307 tw = gmael(R, 1+(w>>1), i+1);
308 z = mul(E, z, tw);
309 i += e;
310 }
311 return z;
312 }
313
314 GEN
gen_powers(GEN x,long l,int use_sqr,void * E,GEN (* sqr)(void *,GEN),GEN (* mul)(void *,GEN,GEN),GEN (* one)(void *))315 gen_powers(GEN x, long l, int use_sqr, void *E, GEN (*sqr)(void*,GEN),
316 GEN (*mul)(void*,GEN,GEN), GEN (*one)(void*))
317 {
318 long i;
319 GEN V = cgetg(l+2,t_VEC);
320 gel(V,1) = one(E); if (l==0) return V;
321 gel(V,2) = gcopy(x); if (l==1) return V;
322 gel(V,3) = sqr(E,x);
323 if (use_sqr)
324 for(i = 4; i < l+2; i++)
325 gel(V,i) = (i&1)? sqr(E,gel(V, (i+1)>>1))
326 : mul(E,gel(V, i-1),x);
327 else
328 for(i = 4; i < l+2; i++)
329 gel(V,i) = mul(E,gel(V,i-1),x);
330 return V;
331 }
332
333 GEN
producttree_scheme(long n)334 producttree_scheme(long n)
335 {
336 GEN v, w;
337 long i, j, k, u, l;
338 if (n<=2) return mkvecsmall(n);
339 u = expu(n-1);
340 v = cgetg(n+1,t_VECSMALL);
341 w = cgetg(n+1,t_VECSMALL);
342 v[1] = n; l = 1;
343 for (i=1; i<=u; i++)
344 {
345 for(j=1, k=1; j<=l; j++, k+=2)
346 {
347 long vj = v[j], v2 = vj>>1;
348 w[k] = vj-v2;
349 w[k+1] = v2;
350 }
351 swap(v,w); l<<=1;
352 }
353 fixlg(v, l+1); set_avma((pari_sp)v); return v;
354 }
355
356 GEN
gen_product(GEN x,void * E,GEN (* mul)(void *,GEN,GEN))357 gen_product(GEN x, void *E, GEN (*mul)(void *,GEN,GEN))
358 {
359 pari_sp av;
360 long i, k, l = lg(x);
361 pari_timer ti;
362 GEN y, v;
363
364 if (l <= 2) return l == 1? gen_1: gcopy(gel(x,1));
365 y = cgetg(l, t_VEC); av = avma;
366 v = producttree_scheme(l-1);
367 l = lg(v);
368 if (DEBUGLEVEL>7) timer_start(&ti);
369 for (k = i = 1; k < l; i += v[k++])
370 gel(y,k) = v[k]==1? gel(x,i): mul(E, gel(x,i), gel(x,i+1));
371 while (k > 2)
372 {
373 long n = k - 1;
374 if (DEBUGLEVEL>7) timer_printf(&ti,"gen_product: remaining objects %ld",n);
375 for (k = i = 1; i < n; i += 2) gel(y,k++) = mul(E, gel(y,i), gel(y,i+1));
376 if (gc_needed(av,1)) gerepilecoeffs(av, y+1, k-1);
377 }
378 return gel(y,1);
379 }
380
381 /***********************************************************************/
382 /** **/
383 /** DISCRETE LOGARITHM **/
384 /** **/
385 /***********************************************************************/
386 static GEN
iter_rho(GEN x,GEN g,GEN q,GEN A,ulong h,void * E,const struct bb_group * grp)387 iter_rho(GEN x, GEN g, GEN q, GEN A, ulong h, void *E, const struct bb_group *grp)
388 {
389 GEN a = gel(A,1), b = gel(A,2), c = gel(A,3);
390 switch((h | grp->hash(a)) % 3UL)
391 {
392 case 0: return mkvec3(grp->pow(E,a,gen_2), Fp_mulu(b,2,q), Fp_mulu(c,2,q));
393 case 1: return mkvec3(grp->mul(E,a,x), addiu(b,1), c);
394 case 2: return mkvec3(grp->mul(E,a,g), b, addiu(c,1));
395 }
396 return NULL; /* LCOV_EXCL_LINE */
397 }
398
399 /*Generic Pollard rho discrete log algorithm*/
400 static GEN
gen_Pollard_log(GEN x,GEN g,GEN q,void * E,const struct bb_group * grp)401 gen_Pollard_log(GEN x, GEN g, GEN q, void *E, const struct bb_group *grp)
402 {
403 pari_sp av=avma;
404 GEN A, B, l, sqrt4q = sqrti(shifti(q,4));
405 ulong i, h = 0, imax = itou_or_0(sqrt4q);
406 if (!imax) imax = ULONG_MAX;
407 do {
408 rho_restart:
409 A = B = mkvec3(x,gen_1,gen_0);
410 i=0;
411 do {
412 if (i>imax)
413 {
414 h++;
415 if (DEBUGLEVEL)
416 pari_warn(warner,"changing Pollard rho hash seed to %ld",h);
417 goto rho_restart;
418 }
419 A = iter_rho(x, g, q, A, h, E, grp);
420 B = iter_rho(x, g, q, B, h, E, grp);
421 B = iter_rho(x, g, q, B, h, E, grp);
422 if (gc_needed(av,2))
423 {
424 if(DEBUGMEM>1) pari_warn(warnmem,"gen_Pollard_log");
425 gerepileall(av, 2, &A, &B);
426 }
427 i++;
428 } while (!grp->equal(gel(A,1), gel(B,1)));
429 gel(A,2) = modii(gel(A,2), q);
430 gel(B,2) = modii(gel(B,2), q);
431 h++;
432 } while (equalii(gel(A,2), gel(B,2)));
433 l = Fp_div(Fp_sub(gel(B,3), gel(A,3),q),Fp_sub(gel(A,2), gel(B,2), q), q);
434 return gerepileuptoint(av, l);
435 }
436
437 /* compute a hash of g^(i-1), 1<=i<=n. Return [sorted hash, perm, g^-n] */
438 GEN
gen_Shanks_init(GEN g,long n,void * E,const struct bb_group * grp)439 gen_Shanks_init(GEN g, long n, void *E, const struct bb_group *grp)
440 {
441 GEN p1 = g, G, perm, table = cgetg(n+1,t_VECSMALL);
442 pari_sp av=avma;
443 long i;
444 table[1] = grp->hash(grp->pow(E,g,gen_0));
445 for (i=2; i<=n; i++)
446 {
447 table[i] = grp->hash(p1);
448 p1 = grp->mul(E,p1,g);
449 if (gc_needed(av,2))
450 {
451 if(DEBUGMEM>1) pari_warn(warnmem,"gen_Shanks_log, baby = %ld", i);
452 p1 = gerepileupto(av, p1);
453 }
454 }
455 G = gerepileupto(av, grp->pow(E,p1,gen_m1)); /* g^-n */
456 perm = vecsmall_indexsort(table);
457 table = vecsmallpermute(table,perm);
458 return mkvec4(table,perm,g,G);
459 }
460 /* T from gen_Shanks_init(g,n). Return v < n*N such that x = g^v or NULL */
461 GEN
gen_Shanks(GEN T,GEN x,ulong N,void * E,const struct bb_group * grp)462 gen_Shanks(GEN T, GEN x, ulong N, void *E, const struct bb_group *grp)
463 {
464 pari_sp av=avma;
465 GEN table = gel(T,1), perm = gel(T,2), g = gel(T,3), G = gel(T,4);
466 GEN p1 = x;
467 long n = lg(table)-1;
468 ulong k;
469 for (k=0; k<N; k++)
470 { /* p1 = x G^k, G = g^-n */
471 long h = grp->hash(p1), i = zv_search(table, h);
472 if (i)
473 {
474 do i--; while (i && table[i] == h);
475 for (i++; i <= n && table[i] == h; i++)
476 {
477 GEN v = addiu(muluu(n,k), perm[i]-1);
478 if (grp->equal(grp->pow(E,g,v),x)) return gerepileuptoint(av,v);
479 if (DEBUGLEVEL)
480 err_printf("gen_Shanks_log: false positive %lu, %lu\n", k,h);
481 }
482 }
483 p1 = grp->mul(E,p1,G);
484 if (gc_needed(av,2))
485 {
486 if(DEBUGMEM>1) pari_warn(warnmem,"gen_Shanks_log, k = %lu", k);
487 p1 = gerepileupto(av, p1);
488 }
489 }
490 return NULL;
491 }
492 /* Generic Shanks baby-step/giant-step algorithm. Return log_g(x), ord g = q.
493 * One-shot: use gen_Shanks_init/log if many logs are desired; early abort
494 * if log < sqrt(q) */
495 static GEN
gen_Shanks_log(GEN x,GEN g,GEN q,void * E,const struct bb_group * grp)496 gen_Shanks_log(GEN x, GEN g, GEN q, void *E, const struct bb_group *grp)
497 {
498 pari_sp av=avma, av1;
499 long lbaby, i, k;
500 GEN p1, table, giant, perm, ginv;
501 p1 = sqrti(q);
502 if (abscmpiu(p1,LGBITS) >= 0)
503 pari_err_OVERFLOW("gen_Shanks_log [order too large]");
504 lbaby = itos(p1)+1; table = cgetg(lbaby+1,t_VECSMALL);
505 ginv = grp->pow(E,g,gen_m1);
506 av1 = avma;
507 for (p1=x, i=1;;i++)
508 {
509 if (grp->equal1(p1)) { set_avma(av); return stoi(i-1); }
510 table[i] = grp->hash(p1); if (i==lbaby) break;
511 p1 = grp->mul(E,p1,ginv);
512 if (gc_needed(av1,2))
513 {
514 if(DEBUGMEM>1) pari_warn(warnmem,"gen_Shanks_log, baby = %ld", i);
515 p1 = gerepileupto(av1, p1);
516 }
517 }
518 p1 = giant = gerepileupto(av1, grp->mul(E,x,grp->pow(E, p1, gen_m1)));
519 perm = vecsmall_indexsort(table);
520 table = vecsmallpermute(table,perm);
521 av1 = avma;
522 for (k=1; k<= lbaby; k++)
523 {
524 long h = grp->hash(p1), i = zv_search(table, h);
525 if (i)
526 {
527 while (table[i] == h && i) i--;
528 for (i++; i <= lbaby && table[i] == h; i++)
529 {
530 GEN v = addiu(mulss(lbaby-1,k),perm[i]-1);
531 if (grp->equal(grp->pow(E,g,v),x)) return gerepileuptoint(av,v);
532 if (DEBUGLEVEL)
533 err_printf("gen_Shanks_log: false positive %ld, %lu\n", k,h);
534 }
535 }
536 p1 = grp->mul(E,p1,giant);
537 if (gc_needed(av1,2))
538 {
539 if(DEBUGMEM>1) pari_warn(warnmem,"gen_Shanks_log, k = %ld", k);
540 p1 = gerepileupto(av1, p1);
541 }
542 }
543 set_avma(av); return cgetg(1, t_VEC); /* no solution */
544 }
545
546 /*Generic discrete logarithme in a group of prime order p*/
547 GEN
gen_plog(GEN x,GEN g,GEN p,void * E,const struct bb_group * grp)548 gen_plog(GEN x, GEN g, GEN p, void *E, const struct bb_group *grp)
549 {
550 if (grp->easylog)
551 {
552 GEN e = grp->easylog(E, x, g, p);
553 if (e) return e;
554 }
555 if (grp->equal1(x)) return gen_0;
556 if (grp->equal(x,g)) return gen_1;
557 if (expi(p)<32) return gen_Shanks_log(x,g,p,E,grp);
558 return gen_Pollard_log(x, g, p, E, grp);
559 }
560
561 GEN
get_arith_ZZM(GEN o)562 get_arith_ZZM(GEN o)
563 {
564 if (!o) return NULL;
565 switch(typ(o))
566 {
567 case t_INT:
568 if (signe(o) > 0) return mkvec2(o, Z_factor(o));
569 break;
570 case t_MAT:
571 if (is_Z_factorpos(o)) return mkvec2(factorback(o), o);
572 break;
573 case t_VEC:
574 if (lg(o) == 3 && signe(gel(o,1)) > 0 && is_Z_factorpos(gel(o,2))) return o;
575 break;
576 }
577 pari_err_TYPE("generic discrete logarithm (order factorization)",o);
578 return NULL; /* LCOV_EXCL_LINE */
579 }
580 GEN
get_arith_Z(GEN o)581 get_arith_Z(GEN o)
582 {
583 if (!o) return NULL;
584 switch(typ(o))
585 {
586 case t_INT:
587 if (signe(o) > 0) return o;
588 break;
589 case t_MAT:
590 o = factorback(o);
591 if (typ(o) == t_INT && signe(o) > 0) return o;
592 break;
593 case t_VEC:
594 if (lg(o) != 3) break;
595 o = gel(o,1);
596 if (typ(o) == t_INT && signe(o) > 0) return o;
597 break;
598 }
599 pari_err_TYPE("generic discrete logarithm (order factorization)",o);
600 return NULL; /* LCOV_EXCL_LINE */
601 }
602
603 /* Generic Pohlig-Hellman discrete logarithm: smallest integer n >= 0 such that
604 * g^n=a. Assume ord(g) | ord; grp->easylog() is an optional trapdoor
605 * function that catches easy logarithms */
606 GEN
gen_PH_log(GEN a,GEN g,GEN ord,void * E,const struct bb_group * grp)607 gen_PH_log(GEN a, GEN g, GEN ord, void *E, const struct bb_group *grp)
608 {
609 pari_sp av = avma;
610 GEN v, ginv, fa, ex;
611 long i, j, l;
612
613 if (grp->equal(g, a)) /* frequent special case */
614 return grp->equal1(g)? gen_0: gen_1;
615 if (grp->easylog)
616 {
617 GEN e = grp->easylog(E, a, g, ord);
618 if (e) return e;
619 }
620 v = get_arith_ZZM(ord);
621 ord= gel(v,1);
622 fa = gel(v,2);
623 ex = gel(fa,2);
624 fa = gel(fa,1); l = lg(fa);
625 ginv = grp->pow(E,g,gen_m1);
626 v = cgetg(l, t_VEC);
627 for (i = 1; i < l; i++)
628 {
629 GEN q = gel(fa,i), qj, gq, nq, ginv0, a0, t0;
630 long e = itos(gel(ex,i));
631 if (DEBUGLEVEL>5)
632 err_printf("Pohlig-Hellman: DL mod %Ps^%ld\n",q,e);
633 qj = new_chunk(e+1);
634 gel(qj,0) = gen_1;
635 gel(qj,1) = q;
636 for (j = 2; j <= e; j++) gel(qj,j) = mulii(gel(qj,j-1), q);
637 t0 = diviiexact(ord, gel(qj,e));
638 a0 = grp->pow(E, a, t0);
639 ginv0 = grp->pow(E, ginv, t0); /* ord(ginv0) | q^e */
640 if (grp->equal1(ginv0)) { gel(v,i) = mkintmod(gen_0, gen_1); continue; }
641 do gq = grp->pow(E,g, mulii(t0, gel(qj,--e))); while (grp->equal1(gq));
642 /* ord(gq) = q */
643 nq = gen_0;
644 for (j = 0;; j++)
645 { /* nq = sum_{i<j} b_i q^i */
646 GEN b = grp->pow(E,a0, gel(qj,e-j));
647 /* cheap early abort: wrong local order */
648 if (j == 0 && !grp->equal1(grp->pow(E,b,q))) {
649 set_avma(av); return cgetg(1, t_VEC);
650 }
651 b = gen_plog(b, gq, q, E, grp);
652 if (typ(b) != t_INT) { set_avma(av); return cgetg(1, t_VEC); }
653 nq = addii(nq, mulii(b, gel(qj,j)));
654 if (j == e) break;
655
656 a0 = grp->mul(E,a0, grp->pow(E,ginv0, b));
657 ginv0 = grp->pow(E,ginv0, q);
658 }
659 gel(v,i) = mkintmod(nq, gel(qj,e+1));
660 }
661 return gerepileuptoint(av, lift(chinese1_coprime_Z(v)));
662 }
663
664 /***********************************************************************/
665 /** **/
666 /** ORDER OF AN ELEMENT **/
667 /** **/
668 /***********************************************************************/
669 /*Find the exact order of a assuming a^o==1*/
670 GEN
gen_order(GEN a,GEN o,void * E,const struct bb_group * grp)671 gen_order(GEN a, GEN o, void *E, const struct bb_group *grp)
672 {
673 pari_sp av = avma;
674 long i, l;
675 GEN m;
676
677 m = get_arith_ZZM(o);
678 if (!m) pari_err_TYPE("gen_order [missing order]",a);
679 o = gel(m,1);
680 m = gel(m,2); l = lgcols(m);
681 for (i = l-1; i; i--)
682 {
683 GEN t, y, p = gcoeff(m,i,1);
684 long j, e = itos(gcoeff(m,i,2));
685 if (l == 2) {
686 t = gen_1;
687 y = a;
688 } else {
689 t = diviiexact(o, powiu(p,e));
690 y = grp->pow(E, a, t);
691 }
692 if (grp->equal1(y)) o = t;
693 else {
694 for (j = 1; j < e; j++)
695 {
696 y = grp->pow(E, y, p);
697 if (grp->equal1(y)) break;
698 }
699 if (j < e) {
700 if (j > 1) p = powiu(p, j);
701 o = mulii(t, p);
702 }
703 }
704 }
705 return gerepilecopy(av, o);
706 }
707
708 /*Find the exact order of a assuming a^o==1, return [order,factor(order)] */
709 GEN
gen_factored_order(GEN a,GEN o,void * E,const struct bb_group * grp)710 gen_factored_order(GEN a, GEN o, void *E, const struct bb_group *grp)
711 {
712 pari_sp av = avma;
713 long i, l, ind;
714 GEN m, F, P;
715
716 m = get_arith_ZZM(o);
717 if (!m) pari_err_TYPE("gen_factored_order [missing order]",a);
718 o = gel(m,1);
719 m = gel(m,2); l = lgcols(m);
720 P = cgetg(l, t_COL); ind = 1;
721 F = cgetg(l, t_COL);
722 for (i = l-1; i; i--)
723 {
724 GEN t, y, p = gcoeff(m,i,1);
725 long j, e = itos(gcoeff(m,i,2));
726 if (l == 2) {
727 t = gen_1;
728 y = a;
729 } else {
730 t = diviiexact(o, powiu(p,e));
731 y = grp->pow(E, a, t);
732 }
733 if (grp->equal1(y)) o = t;
734 else {
735 for (j = 1; j < e; j++)
736 {
737 y = grp->pow(E, y, p);
738 if (grp->equal1(y)) break;
739 }
740 gel(P,ind) = p;
741 gel(F,ind) = utoipos(j);
742 if (j < e) {
743 if (j > 1) p = powiu(p, j);
744 o = mulii(t, p);
745 }
746 ind++;
747 }
748 }
749 setlg(P, ind); P = vecreverse(P);
750 setlg(F, ind); F = vecreverse(F);
751 return gerepilecopy(av, mkvec2(o, mkmat2(P,F)));
752 }
753
754 /* E has order o[1], ..., or o[#o], draw random points until all solutions
755 * but one are eliminated */
756 GEN
gen_select_order(GEN o,void * E,const struct bb_group * grp)757 gen_select_order(GEN o, void *E, const struct bb_group *grp)
758 {
759 pari_sp ltop = avma, btop;
760 GEN lastgood, so, vo;
761 long lo = lg(o), nbo=lo-1;
762 if (nbo == 1) return icopy(gel(o,1));
763 so = ZV_indexsort(o); /* minimize max( o[i+1] - o[i] ) */
764 vo = zero_zv(lo);
765 lastgood = gel(o, so[nbo]);
766 btop = avma;
767 for(;;)
768 {
769 GEN lasto = gen_0;
770 GEN P = grp->rand(E), t = mkvec(gen_0);
771 long i;
772 for (i = 1; i < lo; i++)
773 {
774 GEN newo = gel(o, so[i]);
775 if (vo[i]) continue;
776 t = grp->mul(E,t, grp->pow(E, P, subii(newo,lasto)));/*P^o[i]*/
777 lasto = newo;
778 if (!grp->equal1(t))
779 {
780 if (--nbo == 1) { set_avma(ltop); return icopy(lastgood); }
781 vo[i] = 1;
782 }
783 else
784 lastgood = lasto;
785 }
786 set_avma(btop);
787 }
788 }
789
790 /*******************************************************************/
791 /* */
792 /* n-th ROOT */
793 /* */
794 /*******************************************************************/
795 /* Assume l is prime. Return a generator of the l-th Sylow and set *zeta to an element
796 * of order l.
797 *
798 * q = l^e*r, e>=1, (r,l)=1
799 * UNCLEAN */
800 static GEN
gen_lgener(GEN l,long e,GEN r,GEN * zeta,void * E,const struct bb_group * grp)801 gen_lgener(GEN l, long e, GEN r,GEN *zeta, void *E, const struct bb_group *grp)
802 {
803 const pari_sp av1 = avma;
804 GEN m, m1;
805 long i;
806 for (;; set_avma(av1))
807 {
808 m1 = m = grp->pow(E, grp->rand(E), r);
809 if (grp->equal1(m)) continue;
810 for (i=1; i<e; i++)
811 {
812 m = grp->pow(E,m,l);
813 if (grp->equal1(m)) break;
814 }
815 if (i==e) break;
816 }
817 *zeta = m; return m1;
818 }
819
820 /* Let G be a cyclic group of order o>1. Returns a (random) generator */
821
822 GEN
gen_gener(GEN o,void * E,const struct bb_group * grp)823 gen_gener(GEN o, void *E, const struct bb_group *grp)
824 {
825 pari_sp ltop = avma, av;
826 long i, lpr;
827 GEN F, N, pr, z=NULL;
828 F = get_arith_ZZM(o);
829 N = gel(F,1); pr = gel(F,2); lpr = lgcols(pr);
830 av = avma;
831
832 for (i = 1; i < lpr; i++)
833 {
834 GEN l = gcoeff(pr,i,1);
835 long e = itos(gcoeff(pr,i,2));
836 GEN r = diviiexact(N,powis(l,e));
837 GEN zetan, zl = gen_lgener(l,e,r,&zetan,E,grp);
838 z = i==1 ? zl: grp->mul(E,z,zl);
839 if (gc_needed(av,2))
840 { /* n can have lots of prime factors*/
841 if(DEBUGMEM>1) pari_warn(warnmem,"gen_gener");
842 z = gerepileupto(av, z);
843 }
844 }
845 return gerepileupto(ltop, z);
846 }
847
848 /* solve x^l = a , l prime in G of order q.
849 *
850 * q = (l^e)*r, e >= 1, (r,l) = 1
851 * y is not an l-th power, hence generates the l-Sylow of G
852 * m = y^(q/l) != 1 */
853 static GEN
gen_Shanks_sqrtl(GEN a,GEN l,long e,GEN r,GEN y,GEN m,void * E,const struct bb_group * grp)854 gen_Shanks_sqrtl(GEN a, GEN l, long e, GEN r, GEN y, GEN m,void *E,
855 const struct bb_group *grp)
856 {
857 pari_sp av = avma;
858 long k;
859 GEN p1, u1, u2, v, w, z, dl;
860
861 (void)bezout(r,l,&u1,&u2);
862 v = grp->pow(E,a,u2);
863 w = grp->pow(E,v,l);
864 w = grp->mul(E,w,grp->pow(E,a,gen_m1));
865 while (!grp->equal1(w))
866 {
867 k = 0;
868 p1 = w;
869 do
870 {
871 z = p1; p1 = grp->pow(E,p1,l);
872 k++;
873 } while(!grp->equal1(p1));
874 if (k==e) return gc_NULL(av);
875 dl = gen_plog(z,m,l,E,grp);
876 if (typ(dl) != t_INT) return gc_NULL(av);
877 dl = negi(dl);
878 p1 = grp->pow(E, grp->pow(E,y, dl), powiu(l,e-k-1));
879 m = grp->pow(E,m,dl);
880 e = k;
881 v = grp->mul(E,p1,v);
882 y = grp->pow(E,p1,l);
883 w = grp->mul(E,y,w);
884 if (gc_needed(av,1))
885 {
886 if(DEBUGMEM>1) pari_warn(warnmem,"gen_Shanks_sqrtl");
887 gerepileall(av,4, &y,&v,&w,&m);
888 }
889 }
890 return gerepilecopy(av, v);
891 }
892 /* Return one solution of x^n = a in a cyclic group of order q
893 *
894 * 1) If there is no solution, return NULL.
895 *
896 * 2) If there is a solution, there are exactly m of them [m = gcd(q-1,n)].
897 * If zetan!=NULL, *zetan is set to a primitive m-th root of unity so that
898 * the set of solutions is { x*zetan^k; k=0..m-1 }
899 */
900 GEN
gen_Shanks_sqrtn(GEN a,GEN n,GEN q,GEN * zetan,void * E,const struct bb_group * grp)901 gen_Shanks_sqrtn(GEN a, GEN n, GEN q, GEN *zetan, void *E, const struct bb_group *grp)
902 {
903 pari_sp ltop = avma;
904 GEN m, u1, u2, z;
905 int is_1;
906
907 if (is_pm1(n))
908 {
909 if (zetan) *zetan = grp->pow(E,a,gen_0);
910 return signe(n) < 0? grp->pow(E,a,gen_m1): gcopy(a);
911 }
912 is_1 = grp->equal1(a);
913 if (is_1 && !zetan) return gcopy(a);
914
915 m = bezout(n,q,&u1,&u2);
916 z = grp->pow(E,a,gen_0);
917 if (!is_pm1(m))
918 {
919 GEN F = Z_factor(m);
920 long i, j, e;
921 GEN r, zeta, y, l;
922 pari_sp av1 = avma;
923 for (i = nbrows(F); i; i--)
924 {
925 l = gcoeff(F,i,1);
926 j = itos(gcoeff(F,i,2));
927 e = Z_pvalrem(q,l,&r);
928 y = gen_lgener(l,e,r,&zeta,E,grp);
929 if (zetan) z = grp->mul(E,z, grp->pow(E,y,powiu(l,e-j)));
930 if (!is_1) {
931 do
932 {
933 a = gen_Shanks_sqrtl(a,l,e,r,y,zeta,E,grp);
934 if (!a) return gc_NULL(ltop);
935 } while (--j);
936 }
937 if (gc_needed(ltop,1))
938 { /* n can have lots of prime factors*/
939 if(DEBUGMEM>1) pari_warn(warnmem,"gen_Shanks_sqrtn");
940 gerepileall(av1, zetan? 2: 1, &a, &z);
941 }
942 }
943 }
944 if (!equalii(m, n))
945 a = grp->pow(E,a,modii(u1,q));
946 if (zetan)
947 {
948 *zetan = z;
949 gerepileall(ltop,2,&a,zetan);
950 }
951 else /* is_1 is 0: a was modified above -> gerepileupto valid */
952 a = gerepileupto(ltop, a);
953 return a;
954 }
955
956 /*******************************************************************/
957 /* */
958 /* structure of groups with pairing */
959 /* */
960 /*******************************************************************/
961
962 static GEN
ellgroup_d2(GEN N,GEN d)963 ellgroup_d2(GEN N, GEN d)
964 {
965 GEN r = gcdii(N, d);
966 GEN F1 = gel(Z_factor(r), 1);
967 long i, j, l1 = lg(F1);
968 GEN F = cgetg(3, t_MAT);
969 gel(F,1) = cgetg(l1, t_COL);
970 gel(F,2) = cgetg(l1, t_COL);
971 for (i = 1, j = 1; i < l1; ++i)
972 {
973 long v = Z_pval(N, gel(F1, i));
974 if (v<=1) continue;
975 gcoeff(F, j , 1) = gel(F1, i);
976 gcoeff(F, j++, 2) = stoi(v);
977 }
978 setlg(F[1],j); setlg(F[2],j);
979 return j==1 ? NULL : mkvec2(factorback(F), F);
980 }
981
982 GEN
gen_ellgroup(GEN N,GEN d,GEN * pt_m,void * E,const struct bb_group * grp,GEN pairorder (void * E,GEN P,GEN Q,GEN m,GEN F))983 gen_ellgroup(GEN N, GEN d, GEN *pt_m, void *E, const struct bb_group *grp,
984 GEN pairorder(void *E, GEN P, GEN Q, GEN m, GEN F))
985 {
986 pari_sp av = avma;
987 GEN N0, N1, F;
988 if (pt_m) *pt_m = gen_1;
989 if (is_pm1(N)) return cgetg(1,t_VEC);
990 F = ellgroup_d2(N, d);
991 if (!F) {set_avma(av); return mkveccopy(N);}
992 N0 = gel(F,1); N1 = diviiexact(N, N0);
993 while(1)
994 {
995 pari_sp av2 = avma;
996 GEN P, Q, d, s, t, m;
997
998 P = grp->pow(E,grp->rand(E), N1);
999 s = gen_order(P, F, E, grp); if (equalii(s, N0)) {set_avma(av); return mkveccopy(N);}
1000
1001 Q = grp->pow(E,grp->rand(E), N1);
1002 t = gen_order(Q, F, E, grp); if (equalii(t, N0)) {set_avma(av); return mkveccopy(N);}
1003
1004 m = lcmii(s, t);
1005 d = pairorder(E, P, Q, m, F);
1006 /* structure is [N/d, d] iff m d == N0. Note that N/d = N1 m */
1007 if (is_pm1(d) && equalii(m, N0)) {set_avma(av); return mkveccopy(N);}
1008 if (equalii(mulii(m, d), N0))
1009 {
1010 GEN g = mkvec2(mulii(N1,m), d);
1011 if (pt_m) *pt_m = m;
1012 gerepileall(av,pt_m?2:1,&g,pt_m);
1013 return g;
1014 }
1015 set_avma(av2);
1016 }
1017 }
1018
1019 GEN
gen_ellgens(GEN D1,GEN d2,GEN m,void * E,const struct bb_group * grp,GEN pairorder (void * E,GEN P,GEN Q,GEN m,GEN F))1020 gen_ellgens(GEN D1, GEN d2, GEN m, void *E, const struct bb_group *grp,
1021 GEN pairorder(void *E, GEN P, GEN Q, GEN m, GEN F))
1022 {
1023 pari_sp ltop = avma, av;
1024 GEN F, d1, dm;
1025 GEN P, Q, d, s;
1026 F = get_arith_ZZM(D1);
1027 d1 = gel(F, 1), dm = diviiexact(d1,m);
1028 av = avma;
1029 do
1030 {
1031 set_avma(av);
1032 P = grp->rand(E);
1033 s = gen_order(P, F, E, grp);
1034 } while (!equalii(s, d1));
1035 av = avma;
1036 do
1037 {
1038 set_avma(av);
1039 Q = grp->rand(E);
1040 d = pairorder(E, grp->pow(E, P, dm), grp->pow(E, Q, dm), m, F);
1041 } while (!equalii(d, d2));
1042 return gerepilecopy(ltop, mkvec2(P,Q));
1043 }
1044