CoCalc -- bb

Testing latest pari + WASM + node.js... and it works?! Wow.
pari / 2021-07-01 / pari-2.14.0.alpha / src / basemath / bb_group.c
²⁸⁴⁸⁶ views
License: GPL3
ubuntu2004
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
#define DEBUGLEVEL DEBUGLEVEL_bb_group
25

26
/***********************************************************************/
27
/**                                                                   **/
28
/**                    POWERING                                       **/
29
/**                                                                   **/
30
/***********************************************************************/
31

32
/* return (n>>(i+1-l)) & ((1<<l)-1) */
33
static ulong
34
int_block(GEN n, long i, long l)
35
{
36
  long q = divsBIL(i), r = remsBIL(i)+1, lr;
37
  GEN nw = int_W(n, q);
38
  ulong w = (ulong) *nw, w2;
39
  if (r>=l) return (w>>(r-l))&((1UL<<l)-1);
40
  w &= (1UL<<r)-1; lr = l-r;
41
  w2 = (ulong) *int_precW(nw); w2 >>= (BITS_IN_LONG-lr);
42
  return (w<<lr)|w2;
43
}
44

45
/* assume n != 0, t_INT. Compute x^|n| using sliding window powering */
46
static GEN
47
sliding_window_powu(GEN x, ulong n, long e, void *E, GEN (*sqr)(void*,GEN),
48
                                                     GEN (*mul)(void*,GEN,GEN))
49
{
50
  pari_sp av;
51
  long i, l = expu(n), u = (1UL<<(e-1));
52
  long w, v;
53
  GEN tab = cgetg(1+u, t_VEC);
54
  GEN x2 = sqr(E, x), z = NULL, tw;
55
  gel(tab, 1) = x;
56
  for (i=2; i<=u; i++) gel(tab,i) = mul(E, gel(tab,i-1), x2);
57
  av = avma;
58
  while (l>=0)
59
  {
60
    if (e > l+1) e = l+1;
61
    w = (n>>(l+1-e)) & ((1UL<<e)-1); v = vals(w); l-=e;
62
    tw = gel(tab, 1+(w>>(v+1)));
63
    if (z)
64
    {
65
      for (i=1; i<=e-v; i++) z = sqr(E, z);
66
      z = mul(E, z, tw);
67
    } else z = tw;
68
    for (i=1; i<=v; i++) z = sqr(E, z);
69
    while (l>=0)
70
    {
71
      if (gc_needed(av,1))
72
      {
73
        if (DEBUGMEM>1) pari_warn(warnmem,"sliding_window_powu (%ld)", l);
74
        z = gerepilecopy(av, z);
75
      }
76
      if (n&(1UL<<l)) break;
77
      z = sqr(E, z); l--;
78
    }
79
  }
80
  return z;
81
}
82

83
/* assume n != 0, t_INT. Compute x^|n| using sliding window powering */
84
static GEN
85
sliding_window_pow(GEN x, GEN n, long e, void *E, GEN (*sqr)(void*,GEN),
86
                                                  GEN (*mul)(void*,GEN,GEN))
87
{
88
  pari_sp av;
89
  long i, l = expi(n), u = (1UL<<(e-1));
90
  long w, v;
91
  GEN tab = cgetg(1+u, t_VEC);
92
  GEN x2 = sqr(E, x), z = NULL, tw;
93
  gel(tab, 1) = x;
94
  for (i=2; i<=u; i++) gel(tab,i) = mul(E, gel(tab,i-1), x2);
95
  av = avma;
96
  while (l>=0)
97
  {
98
    if (e > l+1) e = l+1;
99
    w = int_block(n,l,e); v = vals(w); l-=e;
100
    tw = gel(tab, 1+(w>>(v+1)));
101
    if (z)
102
    {
103
      for (i=1; i<=e-v; i++) z = sqr(E, z);
104
      z = mul(E, z, tw);
105
    } else z = tw;
106
    for (i=1; i<=v; i++) z = sqr(E, z);
107
    while (l>=0)
108
    {
109
      if (gc_needed(av,1))
110
      {
111
        if (DEBUGMEM>1) pari_warn(warnmem,"sliding_window_pow (%ld)", l);
112
        z = gerepilecopy(av, z);
113
      }
114
      if (int_bit(n,l)) break;
115
      z = sqr(E, z); l--;
116
    }
117
  }
118
  return z;
119
}
120

121
/* assume n != 0, t_INT. Compute x^|n| using leftright binary powering */
122
static GEN
123
leftright_binary_powu(GEN x, ulong n, void *E, GEN (*sqr)(void*,GEN),
124
                                              GEN (*mul)(void*,GEN,GEN))
125
{
126
  pari_sp av = avma;
127
  GEN  y;
128
  int j;
129

130
  if (n == 1) return x;
131
  y = x; j = 1+bfffo(n);
132
  /* normalize, i.e set highest bit to 1 (we know n != 0) */
133
  n<<=j; j = BITS_IN_LONG-j;
134
  /* first bit is now implicit */
135
  for (; j; n<<=1,j--)
136
  {
137
    y = sqr(E,y);
138
    if (n & HIGHBIT) y = mul(E,y,x); /* first bit set: multiply by base */
139
    if (gc_needed(av,1))
140
    {
141
      if (DEBUGMEM>1) pari_warn(warnmem,"leftright_powu (%d)", j);
142
      y = gerepilecopy(av, y);
143
    }
144
  }
145
  return y;
146
}
147

148
GEN
149
gen_powu_i(GEN x, ulong n, void *E, GEN (*sqr)(void*,GEN),
150
                                    GEN (*mul)(void*,GEN,GEN))
151
{
152
  long l;
153
  if (n == 1) return x;
154
  l = expu(n);
155
  if (l<=8)
156
    return leftright_binary_powu(x, n, E, sqr, mul);
157
  else
158
    return sliding_window_powu(x, n, l<=24? 2: 3, E, sqr, mul);
159
}
160

161
GEN
162
gen_powu(GEN x, ulong n, void *E, GEN (*sqr)(void*,GEN),
163
                                  GEN (*mul)(void*,GEN,GEN))
164
{
165
  pari_sp av = avma;
166
  if (n == 1) return gcopy(x);
167
  return gerepilecopy(av, gen_powu_i(x,n,E,sqr,mul));
168
}
169

170
GEN
171
gen_pow_i(GEN x, GEN n, void *E, GEN (*sqr)(void*,GEN),
172
                                 GEN (*mul)(void*,GEN,GEN))
173
{
174
  long l, e;
175
  if (lgefint(n)==3) return gen_powu_i(x, uel(n,2), E, sqr, mul);
176
  l = expi(n);
177
  if      (l<=64)  e = 3;
178
  else if (l<=160) e = 4;
179
  else if (l<=384) e = 5;
180
  else if (l<=896) e = 6;
181
  else             e = 7;
182
  return sliding_window_pow(x, n, e, E, sqr, mul);
183
}
184

185
GEN
186
gen_pow(GEN x, GEN n, void *E, GEN (*sqr)(void*,GEN),
187
                               GEN (*mul)(void*,GEN,GEN))
188
{
189
  pari_sp av = avma;
190
  return gerepilecopy(av, gen_pow_i(x,n,E,sqr,mul));
191
}
192

193
/* assume n > 0. Compute x^n using left-right binary powering */
194
GEN
195
gen_powu_fold_i(GEN x, ulong n, void *E, GEN  (*sqr)(void*,GEN),
196
                                         GEN (*msqr)(void*,GEN))
197
{
198
  pari_sp av = avma;
199
  GEN y;
200
  int j;
201

202
  if (n == 1) return x;
203
  y = x; j = 1+bfffo(n);
204
  /* normalize, i.e set highest bit to 1 (we know n != 0) */
205
  n<<=j; j = BITS_IN_LONG-j;
206
  /* first bit is now implicit */
207
  for (; j; n<<=1,j--)
208
  {
209
    if (n & HIGHBIT) y = msqr(E,y); /* first bit set: multiply by base */
210
    else y = sqr(E,y);
211
    if (gc_needed(av,1))
212
    {
213
      if (DEBUGMEM>1) pari_warn(warnmem,"gen_powu_fold (%d)", j);
214
      y = gerepilecopy(av, y);
215
    }
216
  }
217
  return y;
218
}
219
GEN
220
gen_powu_fold(GEN x, ulong n, void *E, GEN (*sqr)(void*,GEN),
221
                                       GEN (*msqr)(void*,GEN))
222
{
223
  pari_sp av = avma;
224
  if (n == 1) return gcopy(x);
225
  return gerepilecopy(av, gen_powu_fold_i(x,n,E,sqr,msqr));
226
}
227

228
/* assume N != 0, t_INT. Compute x^|N| using left-right binary powering */
229
GEN
230
gen_pow_fold_i(GEN x, GEN N, void *E, GEN (*sqr)(void*,GEN),
231
                                      GEN (*msqr)(void*,GEN))
232
{
233
  long ln = lgefint(N);
234
  if (ln == 3) return gen_powu_fold_i(x, N[2], E, sqr, msqr);
235
  else
236
  {
237
    GEN nd = int_MSW(N), y = x;
238
    ulong n = *nd;
239
    long i;
240
    int j;
241
    pari_sp av = avma;
242

243
    if (n == 1)
244
      j = 0;
245
    else
246
    {
247
      j = 1+bfffo(n); /* < BIL */
248
      /* normalize, i.e set highest bit to 1 (we know n != 0) */
249
      n <<= j; j = BITS_IN_LONG - j;
250
    }
251
    /* first bit is now implicit */
252
    for (i=ln-2;;)
253
    {
254
      for (; j; n<<=1,j--)
255
      {
256
        if (n & HIGHBIT) y = msqr(E,y); /* first bit set: multiply by base */
257
        else y = sqr(E,y);
258
        if (gc_needed(av,1))
259
        {
260
          if (DEBUGMEM>1) pari_warn(warnmem,"gen_pow_fold (%d)", j);
261
          y = gerepilecopy(av, y);
262
        }
263
      }
264
      if (--i == 0) return y;
265
      nd = int_precW(nd);
266
      n = *nd; j = BITS_IN_LONG;
267
    }
268
  }
269
}
270
GEN
271
gen_pow_fold(GEN x, GEN n, void *E, GEN (*sqr)(void*,GEN),
272
                                    GEN (*msqr)(void*,GEN))
273
{
274
  pari_sp av = avma;
275
  return gerepilecopy(av, gen_pow_fold_i(x,n,E,sqr,msqr));
276
}
277

278
GEN
279
gen_pow_init(GEN x, GEN n, long k, void *E, GEN (*sqr)(void*,GEN), GEN (*mul)(void*,GEN,GEN))
280
{
281
  long i, j, l = expi(n);
282
  long m = 1UL<<(k-1);
283
  GEN x2 = sqr(E, x), y = gcopy(x);
284
  GEN R = cgetg(m+1, t_VEC);
285
  for(i = 1; i <= m; i++)
286
  {
287
    GEN C = cgetg(l+1, t_VEC);
288
    gel(C,1) = y;
289
    for(j = 2; j <= l; j++)
290
      gel(C,j) = sqr(E, gel(C,j-1));
291
    gel(R,i) = C;
292
    y = mul(E, y, x2);
293
  }
294
  return R;
295
}
296

297
GEN
298
gen_pow_table(GEN R, GEN n, void *E, GEN (*one)(void*), GEN (*mul)(void*,GEN,GEN))
299
{
300
  long e = expu(lg(R)-1) + 1;
301
  long l = expi(n);
302
  long i, w;
303
  GEN z = one(E), tw;
304
  for(i=0; i<=l; )
305
  {
306
    if (int_bit(n, i)==0) { i++; continue; }
307
    if (i+e-1>l) e = l+1-i;
308
    w = int_block(n,i+e-1,e);
309
    tw = gmael(R, 1+(w>>1), i+1);
310
    z = mul(E, z, tw);
311
    i += e;
312
  }
313
  return z;
314
}
315

316
GEN
317
gen_powers(GEN x, long l, int use_sqr, void *E, GEN (*sqr)(void*,GEN),
318
                                      GEN (*mul)(void*,GEN,GEN), GEN (*one)(void*))
319
{
320
  long i;
321
  GEN V = cgetg(l+2,t_VEC);
322
  gel(V,1) = one(E); if (l==0) return V;
323
  gel(V,2) = gcopy(x); if (l==1) return V;
324
  gel(V,3) = sqr(E,x);
325
  if (use_sqr)
326
    for(i = 4; i < l+2; i++)
327
      gel(V,i) = (i&1)? sqr(E,gel(V, (i+1)>>1))
328
                      : mul(E,gel(V, i-1),x);
329
  else
330
    for(i = 4; i < l+2; i++)
331
      gel(V,i) = mul(E,gel(V,i-1),x);
332
  return V;
333
}
334

335
GEN
336
producttree_scheme(long n)
337
{
338
  GEN v, w;
339
  long i, j, k, u, l;
340
  if (n<=2) return mkvecsmall(n);
341
  u = expu(n-1);
342
  v = cgetg(n+1,t_VECSMALL);
343
  w = cgetg(n+1,t_VECSMALL);
344
  v[1] = n; l = 1;
345
  for (i=1; i<=u; i++)
346
  {
347
    for(j=1, k=1; j<=l; j++, k+=2)
348
    {
349
      long vj = v[j], v2 = vj>>1;
350
      w[k]    = vj-v2;
351
      w[k+1]  = v2;
352
    }
353
    swap(v,w); l<<=1;
354
  }
355
  fixlg(v, l+1); set_avma((pari_sp)v); return v;
356
}
357

358
GEN
359
gen_product(GEN x, void *E, GEN (*mul)(void *,GEN,GEN))
360
{
361
  pari_sp av;
362
  long i, k, l = lg(x);
363
  pari_timer ti;
364
  GEN y, v;
365

366
  if (l <= 2) return l == 1? gen_1: gcopy(gel(x,1));
367
  y = cgetg(l, t_VEC); av = avma;
368
  v = producttree_scheme(l-1);
369
  l = lg(v);
370
  if (DEBUGLEVEL>7) timer_start(&ti);
371
  for (k = i = 1; k < l; i += v[k++])
372
    gel(y,k) = v[k]==1? gel(x,i): mul(E, gel(x,i), gel(x,i+1));
373
  while (k > 2)
374
  {
375
    long n = k - 1;
376
    if (DEBUGLEVEL>7) timer_printf(&ti,"gen_product: remaining objects %ld",n);
377
    for (k = i = 1; i < n; i += 2) gel(y,k++) = mul(E, gel(y,i), gel(y,i+1));
378
    if (gc_needed(av,1)) gerepilecoeffs(av, y+1, k-1);
379
  }
380
  return gel(y,1);
381
}
382

383
/***********************************************************************/
384
/**                                                                   **/
385
/**                    DISCRETE LOGARITHM                             **/
386
/**                                                                   **/
387
/***********************************************************************/
388
static GEN
389
iter_rho(GEN x, GEN g, GEN q, GEN A, ulong h, void *E, const struct bb_group *grp)
390
{
391
  GEN a = gel(A,1), b = gel(A,2), c = gel(A,3);
392
  switch((h | grp->hash(a)) % 3UL)
393
  {
394
    case 0: return mkvec3(grp->pow(E,a,gen_2), Fp_mulu(b,2,q), Fp_mulu(c,2,q));
395
    case 1: return mkvec3(grp->mul(E,a,x), addiu(b,1), c);
396
    case 2: return mkvec3(grp->mul(E,a,g), b, addiu(c,1));
397
  }
398
  return NULL; /* LCOV_EXCL_LINE */
399
}
400

401
/*Generic Pollard rho discrete log algorithm*/
402
static GEN
403
gen_Pollard_log(GEN x, GEN g, GEN q, void *E, const struct bb_group *grp)
404
{
405
  pari_sp av=avma;
406
  GEN A, B, l, sqrt4q = sqrti(shifti(q,4));
407
  ulong i, h = 0, imax = itou_or_0(sqrt4q);
408
  if (!imax) imax = ULONG_MAX;
409
  do {
410
 rho_restart:
411
    A = B = mkvec3(x,gen_1,gen_0);
412
    i=0;
413
    do {
414
      if (i>imax)
415
      {
416
        h++;
417
        if (DEBUGLEVEL)
418
          pari_warn(warner,"changing Pollard rho hash seed to %ld",h);
419
        goto rho_restart;
420
      }
421
      A = iter_rho(x, g, q, A, h, E, grp);
422
      B = iter_rho(x, g, q, B, h, E, grp);
423
      B = iter_rho(x, g, q, B, h, E, grp);
424
      if (gc_needed(av,2))
425
      {
426
        if(DEBUGMEM>1) pari_warn(warnmem,"gen_Pollard_log");
427
        gerepileall(av, 2, &A, &B);
428
      }
429
      i++;
430
    } while (!grp->equal(gel(A,1), gel(B,1)));
431
    gel(A,2) = modii(gel(A,2), q);
432
    gel(B,2) = modii(gel(B,2), q);
433
    h++;
434
  } while (equalii(gel(A,2), gel(B,2)));
435
  l = Fp_div(Fp_sub(gel(B,3), gel(A,3),q),Fp_sub(gel(A,2), gel(B,2), q), q);
436
  return gerepileuptoint(av, l);
437
}
438

439
/* compute a hash of g^(i-1), 1<=i<=n. Return [sorted hash, perm, g^-n] */
440
GEN
441
gen_Shanks_init(GEN g, long n, void *E, const struct bb_group *grp)
442
{
443
  GEN p1 = g, G, perm, table = cgetg(n+1,t_VECSMALL);
444
  pari_sp av=avma;
445
  long i;
446
  table[1] = grp->hash(grp->pow(E,g,gen_0));
447
  for (i=2; i<=n; i++)
448
  {
449
    table[i] = grp->hash(p1);
450
    p1 = grp->mul(E,p1,g);
451
    if (gc_needed(av,2))
452
    {
453
      if(DEBUGMEM>1) pari_warn(warnmem,"gen_Shanks_log, baby = %ld", i);
454
      p1 = gerepileupto(av, p1);
455
    }
456
  }
457
  G = gerepileupto(av, grp->pow(E,p1,gen_m1)); /* g^-n */
458
  perm = vecsmall_indexsort(table);
459
  table = vecsmallpermute(table,perm);
460
  return mkvec4(table,perm,g,G);
461
}
462
/* T from gen_Shanks_init(g,n). Return v < n*N such that x = g^v or NULL */
463
GEN
464
gen_Shanks(GEN T, GEN x, ulong N, void *E, const struct bb_group *grp)
465
{
466
  pari_sp av=avma;
467
  GEN table = gel(T,1), perm = gel(T,2), g = gel(T,3), G = gel(T,4);
468
  GEN p1 = x;
469
  long n = lg(table)-1;
470
  ulong k;
471
  for (k=0; k<N; k++)
472
  { /* p1 = x G^k, G = g^-n */
473
    long h = grp->hash(p1), i = zv_search(table, h);
474
    if (i)
475
    {
476
      do i--; while (i && table[i] == h);
477
      for (i++; i <= n && table[i] == h; i++)
478
      {
479
        GEN v = addiu(muluu(n,k), perm[i]-1);
480
        if (grp->equal(grp->pow(E,g,v),x)) return gerepileuptoint(av,v);
481
        if (DEBUGLEVEL)
482
          err_printf("gen_Shanks_log: false positive %lu, %lu\n", k,h);
483
      }
484
    }
485
    p1 = grp->mul(E,p1,G);
486
    if (gc_needed(av,2))
487
    {
488
      if(DEBUGMEM>1) pari_warn(warnmem,"gen_Shanks_log, k = %lu", k);
489
      p1 = gerepileupto(av, p1);
490
    }
491
  }
492
  return NULL;
493
}
494
/* Generic Shanks baby-step/giant-step algorithm. Return log_g(x), ord g = q.
495
 * One-shot: use gen_Shanks_init/log if many logs are desired; early abort
496
 * if log < sqrt(q) */
497
static GEN
498
gen_Shanks_log(GEN x, GEN g, GEN q, void *E, const struct bb_group *grp)
499
{
500
  pari_sp av=avma, av1;
501
  long lbaby, i, k;
502
  GEN p1, table, giant, perm, ginv;
503
  p1 = sqrti(q);
504
  if (abscmpiu(p1,LGBITS) >= 0)
505
    pari_err_OVERFLOW("gen_Shanks_log [order too large]");
506
  lbaby = itos(p1)+1; table = cgetg(lbaby+1,t_VECSMALL);
507
  ginv = grp->pow(E,g,gen_m1);
508
  av1 = avma;
509
  for (p1=x, i=1;;i++)
510
  {
511
    if (grp->equal1(p1)) { set_avma(av); return stoi(i-1); }
512
    table[i] = grp->hash(p1); if (i==lbaby) break;
513
    p1 = grp->mul(E,p1,ginv);
514
    if (gc_needed(av1,2))
515
    {
516
      if(DEBUGMEM>1) pari_warn(warnmem,"gen_Shanks_log, baby = %ld", i);
517
      p1 = gerepileupto(av1, p1);
518
    }
519
  }
520
  p1 = giant = gerepileupto(av1, grp->mul(E,x,grp->pow(E, p1, gen_m1)));
521
  perm = vecsmall_indexsort(table);
522
  table = vecsmallpermute(table,perm);
523
  av1 = avma;
524
  for (k=1; k<= lbaby; k++)
525
  {
526
    long h = grp->hash(p1), i = zv_search(table, h);
527
    if (i)
528
    {
529
      while (table[i] == h && i) i--;
530
      for (i++; i <= lbaby && table[i] == h; i++)
531
      {
532
        GEN v = addiu(mulss(lbaby-1,k),perm[i]-1);
533
        if (grp->equal(grp->pow(E,g,v),x)) return gerepileuptoint(av,v);
534
        if (DEBUGLEVEL)
535
          err_printf("gen_Shanks_log: false positive %ld, %lu\n", k,h);
536
      }
537
    }
538
    p1 = grp->mul(E,p1,giant);
539
    if (gc_needed(av1,2))
540
    {
541
      if(DEBUGMEM>1) pari_warn(warnmem,"gen_Shanks_log, k = %ld", k);
542
      p1 = gerepileupto(av1, p1);
543
    }
544
  }
545
  set_avma(av); return cgetg(1, t_VEC); /* no solution */
546
}
547

548
/*Generic discrete logarithme in a group of prime order p*/
549
GEN
550
gen_plog(GEN x, GEN g, GEN p, void *E, const struct bb_group *grp)
551
{
552
  if (grp->easylog)
553
  {
554
    GEN e = grp->easylog(E, x, g, p);
555
    if (e) return e;
556
  }
557
  if (grp->equal1(x)) return gen_0;
558
  if (grp->equal(x,g)) return gen_1;
559
  if (expi(p)<32) return gen_Shanks_log(x,g,p,E,grp);
560
  return gen_Pollard_log(x, g, p, E, grp);
561
}
562

563
GEN
564
get_arith_ZZM(GEN o)
565
{
566
  if (!o) return NULL;
567
  switch(typ(o))
568
  {
569
    case t_INT:
570
      if (signe(o) > 0) return mkvec2(o, Z_factor(o));
571
      break;
572
    case t_MAT:
573
      if (is_Z_factorpos(o)) return mkvec2(factorback(o), o);
574
      break;
575
    case t_VEC:
576
      if (lg(o) == 3 && signe(gel(o,1)) > 0 && is_Z_factorpos(gel(o,2))) return o;
577
      break;
578
  }
579
  pari_err_TYPE("generic discrete logarithm (order factorization)",o);
580
  return NULL; /* LCOV_EXCL_LINE */
581
}
582
GEN
583
get_arith_Z(GEN o)
584
{
585
  if (!o) return NULL;
586
  switch(typ(o))
587
  {
588
    case t_INT:
589
      if (signe(o) > 0) return o;
590
      break;
591
    case t_MAT:
592
      o = factorback(o);
593
      if (typ(o) == t_INT && signe(o) > 0) return o;
594
      break;
595
    case t_VEC:
596
      if (lg(o) != 3) break;
597
      o = gel(o,1);
598
      if (typ(o) == t_INT && signe(o) > 0) return o;
599
      break;
600
  }
601
  pari_err_TYPE("generic discrete logarithm (order factorization)",o);
602
  return NULL; /* LCOV_EXCL_LINE */
603
}
604

605
/* Generic Pohlig-Hellman discrete logarithm: smallest integer n >= 0 such that
606
 * g^n=a. Assume ord(g) | ord; grp->easylog() is an optional trapdoor
607
 * function that catches easy logarithms */
608
GEN
609
gen_PH_log(GEN a, GEN g, GEN ord, void *E, const struct bb_group *grp)
610
{
611
  pari_sp av = avma;
612
  GEN v, ginv, fa, ex;
613
  long i, j, l;
614

615
  if (grp->equal(g, a)) /* frequent special case */
616
    return grp->equal1(g)? gen_0: gen_1;
617
  if (grp->easylog)
618
  {
619
    GEN e = grp->easylog(E, a, g, ord);
620
    if (e) return e;
621
  }
622
  v = get_arith_ZZM(ord);
623
  ord= gel(v,1);
624
  fa = gel(v,2);
625
  ex = gel(fa,2);
626
  fa = gel(fa,1); l = lg(fa);
627
  ginv = grp->pow(E,g,gen_m1);
628
  v = cgetg(l, t_VEC);
629
  for (i = 1; i < l; i++)
630
  {
631
    GEN q = gel(fa,i), qj, gq, nq, ginv0, a0, t0;
632
    long e = itos(gel(ex,i));
633
    if (DEBUGLEVEL>5)
634
      err_printf("Pohlig-Hellman: DL mod %Ps^%ld\n",q,e);
635
    qj = new_chunk(e+1);
636
    gel(qj,0) = gen_1;
637
    gel(qj,1) = q;
638
    for (j = 2; j <= e; j++) gel(qj,j) = mulii(gel(qj,j-1), q);
639
    t0 = diviiexact(ord, gel(qj,e));
640
    a0 = grp->pow(E, a, t0);
641
    ginv0 = grp->pow(E, ginv, t0); /* ord(ginv0) | q^e */
642
    if (grp->equal1(ginv0)) { gel(v,i) = mkintmod(gen_0, gen_1); continue; }
643
    do gq = grp->pow(E,g, mulii(t0, gel(qj,--e))); while (grp->equal1(gq));
644
    /* ord(gq) = q */
645
    nq = gen_0;
646
    for (j = 0;; j++)
647
    { /* nq = sum_{i<j} b_i q^i */
648
      GEN b = grp->pow(E,a0, gel(qj,e-j));
649
      /* cheap early abort: wrong local order */
650
      if (j == 0 && !grp->equal1(grp->pow(E,b,q))) {
651
        set_avma(av); return cgetg(1, t_VEC);
652
      }
653
      b = gen_plog(b, gq, q, E, grp);
654
      if (typ(b) != t_INT) { set_avma(av); return cgetg(1, t_VEC); }
655
      nq = addii(nq, mulii(b, gel(qj,j)));
656
      if (j == e) break;
657

658
      a0 = grp->mul(E,a0, grp->pow(E,ginv0, b));
659
      ginv0 = grp->pow(E,ginv0, q);
660
    }
661
    gel(v,i) = mkintmod(nq, gel(qj,e+1));
662
  }
663
  return gerepileuptoint(av, lift(chinese1_coprime_Z(v)));
664
}
665

666
/***********************************************************************/
667
/**                                                                   **/
668
/**                    ORDER OF AN ELEMENT                            **/
669
/**                                                                   **/
670
/***********************************************************************/
671
/*Find the exact order of a assuming a^o==1*/
672
GEN
673
gen_order(GEN a, GEN o, void *E, const struct bb_group *grp)
674
{
675
  pari_sp av = avma;
676
  long i, l;
677
  GEN m;
678

679
  m = get_arith_ZZM(o);
680
  if (!m) pari_err_TYPE("gen_order [missing order]",a);
681
  o = gel(m,1);
682
  m = gel(m,2); l = lgcols(m);
683
  for (i = l-1; i; i--)
684
  {
685
    GEN t, y, p = gcoeff(m,i,1);
686
    long j, e = itos(gcoeff(m,i,2));
687
    if (l == 2) {
688
      t = gen_1;
689
      y = a;
690
    } else {
691
      t = diviiexact(o, powiu(p,e));
692
      y = grp->pow(E, a, t);
693
    }
694
    if (grp->equal1(y)) o = t;
695
    else {
696
      for (j = 1; j < e; j++)
697
      {
698
        y = grp->pow(E, y, p);
699
        if (grp->equal1(y)) break;
700
      }
701
      if (j < e) {
702
        if (j > 1) p = powiu(p, j);
703
        o = mulii(t, p);
704
      }
705
    }
706
  }
707
  return gerepilecopy(av, o);
708
}
709

710
/*Find the exact order of a assuming a^o==1, return [order,factor(order)] */
711
GEN
712
gen_factored_order(GEN a, GEN o, void *E, const struct bb_group *grp)
713
{
714
  pari_sp av = avma;
715
  long i, l, ind;
716
  GEN m, F, P;
717

718
  m = get_arith_ZZM(o);
719
  if (!m) pari_err_TYPE("gen_factored_order [missing order]",a);
720
  o = gel(m,1);
721
  m = gel(m,2); l = lgcols(m);
722
  P = cgetg(l, t_COL); ind = 1;
723
  F = cgetg(l, t_COL);
724
  for (i = l-1; i; i--)
725
  {
726
    GEN t, y, p = gcoeff(m,i,1);
727
    long j, e = itos(gcoeff(m,i,2));
728
    if (l == 2) {
729
      t = gen_1;
730
      y = a;
731
    } else {
732
      t = diviiexact(o, powiu(p,e));
733
      y = grp->pow(E, a, t);
734
    }
735
    if (grp->equal1(y)) o = t;
736
    else {
737
      for (j = 1; j < e; j++)
738
      {
739
        y = grp->pow(E, y, p);
740
        if (grp->equal1(y)) break;
741
      }
742
      gel(P,ind) = p;
743
      gel(F,ind) = utoipos(j);
744
      if (j < e) {
745
        if (j > 1) p = powiu(p, j);
746
        o = mulii(t, p);
747
      }
748
      ind++;
749
    }
750
  }
751
  setlg(P, ind); P = vecreverse(P);
752
  setlg(F, ind); F = vecreverse(F);
753
  return gerepilecopy(av, mkvec2(o, mkmat2(P,F)));
754
}
755

756
/* E has order o[1], ..., or o[#o], draw random points until all solutions
757
 * but one are eliminated */
758
GEN
759
gen_select_order(GEN o, void *E, const struct bb_group *grp)
760
{
761
  pari_sp ltop = avma, btop;
762
  GEN lastgood, so, vo;
763
  long lo = lg(o), nbo=lo-1;
764
  if (nbo == 1) return icopy(gel(o,1));
765
  so = ZV_indexsort(o); /* minimize max( o[i+1] - o[i] ) */
766
  vo = zero_zv(lo);
767
  lastgood = gel(o, so[nbo]);
768
  btop = avma;
769
  for(;;)
770
  {
771
    GEN lasto = gen_0;
772
    GEN P = grp->rand(E), t = mkvec(gen_0);
773
    long i;
774
    for (i = 1; i < lo; i++)
775
    {
776
      GEN newo = gel(o, so[i]);
777
      if (vo[i]) continue;
778
      t = grp->mul(E,t, grp->pow(E, P, subii(newo,lasto)));/*P^o[i]*/
779
      lasto = newo;
780
      if (!grp->equal1(t))
781
      {
782
        if (--nbo == 1) { set_avma(ltop); return icopy(lastgood); }
783
        vo[i] = 1;
784
      }
785
      else
786
        lastgood = lasto;
787
    }
788
    set_avma(btop);
789
  }
790
}
791

792
/*******************************************************************/
793
/*                                                                 */
794
/*                          n-th ROOT                              */
795
/*                                                                 */
796
/*******************************************************************/
797
/* Assume l is prime. Return a generator of the l-th Sylow and set *zeta to an element
798
 * of order l.
799
 *
800
 * q = l^e*r, e>=1, (r,l)=1
801
 * UNCLEAN */
802
static GEN
803
gen_lgener(GEN l, long e, GEN r,GEN *zeta, void *E, const struct bb_group *grp)
804
{
805
  const pari_sp av1 = avma;
806
  GEN m, m1;
807
  long i;
808
  for (;; set_avma(av1))
809
  {
810
    m1 = m = grp->pow(E, grp->rand(E), r);
811
    if (grp->equal1(m)) continue;
812
    for (i=1; i<e; i++)
813
    {
814
      m = grp->pow(E,m,l);
815
      if (grp->equal1(m)) break;
816
    }
817
    if (i==e) break;
818
  }
819
  *zeta = m; return m1;
820
}
821

822
/* Let G be a cyclic group of order o>1. Returns a (random) generator */
823

824
GEN
825
gen_gener(GEN o, void *E, const struct bb_group *grp)
826
{
827
  pari_sp ltop = avma, av;
828
  long i, lpr;
829
  GEN F, N, pr, z=NULL;
830
  F = get_arith_ZZM(o);
831
  N = gel(F,1); pr = gel(F,2); lpr = lgcols(pr);
832
  av = avma;
833

834
  for (i = 1; i < lpr; i++)
835
  {
836
    GEN l = gcoeff(pr,i,1);
837
    long e = itos(gcoeff(pr,i,2));
838
    GEN r = diviiexact(N,powis(l,e));
839
    GEN zetan, zl = gen_lgener(l,e,r,&zetan,E,grp);
840
    z = i==1 ? zl: grp->mul(E,z,zl);
841
    if (gc_needed(av,2))
842
    { /* n can have lots of prime factors*/
843
      if(DEBUGMEM>1) pari_warn(warnmem,"gen_gener");
844
      z = gerepileupto(av, z);
845
    }
846
  }
847
  return gerepileupto(ltop, z);
848
}
849

850
/* solve x^l = a , l prime in G of order q.
851
 *
852
 * q =  (l^e)*r, e >= 1, (r,l) = 1
853
 * y is not an l-th power, hence generates the l-Sylow of G
854
 * m = y^(q/l) != 1 */
855
static GEN
856
gen_Shanks_sqrtl(GEN a, GEN l, long e, GEN r, GEN y, GEN m,void *E,
857
                 const struct bb_group *grp)
858
{
859
  pari_sp av = avma;
860
  long k;
861
  GEN p1, u1, u2, v, w, z, dl;
862

863
  (void)bezout(r,l,&u1,&u2);
864
  v = grp->pow(E,a,u2);
865
  w = grp->pow(E,v,l);
866
  w = grp->mul(E,w,grp->pow(E,a,gen_m1));
867
  while (!grp->equal1(w))
868
  {
869
    k = 0;
870
    p1 = w;
871
    do
872
    {
873
      z = p1; p1 = grp->pow(E,p1,l);
874
      k++;
875
    } while(!grp->equal1(p1));
876
    if (k==e) return gc_NULL(av);
877
    dl = gen_plog(z,m,l,E,grp);
878
    if (typ(dl) != t_INT) return gc_NULL(av);
879
    dl = negi(dl);
880
    p1 = grp->pow(E, grp->pow(E,y, dl), powiu(l,e-k-1));
881
    m = grp->pow(E,m,dl);
882
    e = k;
883
    v = grp->mul(E,p1,v);
884
    y = grp->pow(E,p1,l);
885
    w = grp->mul(E,y,w);
886
    if (gc_needed(av,1))
887
    {
888
      if(DEBUGMEM>1) pari_warn(warnmem,"gen_Shanks_sqrtl");
889
      gerepileall(av,4, &y,&v,&w,&m);
890
    }
891
  }
892
  return gerepilecopy(av, v);
893
}
894
/* Return one solution of x^n = a in a cyclic group of order q
895
 *
896
 * 1) If there is no solution, return NULL.
897
 *
898
 * 2) If there is a solution, there are exactly m of them [m = gcd(q-1,n)].
899
 * If zetan!=NULL, *zetan is set to a primitive m-th root of unity so that
900
 * the set of solutions is { x*zetan^k; k=0..m-1 }
901
 */
902
GEN
903
gen_Shanks_sqrtn(GEN a, GEN n, GEN q, GEN *zetan, void *E, const struct bb_group *grp)
904
{
905
  pari_sp ltop = avma;
906
  GEN m, u1, u2, z;
907
  int is_1;
908

909
  if (is_pm1(n))
910
  {
911
    if (zetan) *zetan = grp->pow(E,a,gen_0);
912
    return signe(n) < 0? grp->pow(E,a,gen_m1): gcopy(a);
913
  }
914
  is_1 = grp->equal1(a);
915
  if (is_1 && !zetan) return gcopy(a);
916

917
  m = bezout(n,q,&u1,&u2);
918
  z = grp->pow(E,a,gen_0);
919
  if (!is_pm1(m))
920
  {
921
    GEN F = Z_factor(m);
922
    long i, j, e;
923
    GEN r, zeta, y, l;
924
    pari_sp av1 = avma;
925
    for (i = nbrows(F); i; i--)
926
    {
927
      l = gcoeff(F,i,1);
928
      j = itos(gcoeff(F,i,2));
929
      e = Z_pvalrem(q,l,&r);
930
      y = gen_lgener(l,e,r,&zeta,E,grp);
931
      if (zetan) z = grp->mul(E,z, grp->pow(E,y,powiu(l,e-j)));
932
      if (!is_1) {
933
        do
934
        {
935
          a = gen_Shanks_sqrtl(a,l,e,r,y,zeta,E,grp);
936
          if (!a) return gc_NULL(ltop);
937
        } while (--j);
938
      }
939
      if (gc_needed(ltop,1))
940
      { /* n can have lots of prime factors*/
941
        if(DEBUGMEM>1) pari_warn(warnmem,"gen_Shanks_sqrtn");
942
        gerepileall(av1, zetan? 2: 1, &a, &z);
943
      }
944
    }
945
  }
946
  if (!equalii(m, n))
947
    a = grp->pow(E,a,modii(u1,q));
948
  if (zetan)
949
  {
950
    *zetan = z;
951
    gerepileall(ltop,2,&a,zetan);
952
  }
953
  else /* is_1 is 0: a was modified above -> gerepileupto valid */
954
    a = gerepileupto(ltop, a);
955
  return a;
956
}
957

958
/*******************************************************************/
959
/*                                                                 */
960
/*               structure of groups with pairing                  */
961
/*                                                                 */
962
/*******************************************************************/
963

964
static GEN
965
ellgroup_d2(GEN N, GEN d)
966
{
967
  GEN r = gcdii(N, d);
968
  GEN F1 = gel(Z_factor(r), 1);
969
  long i, j, l1 = lg(F1);
970
  GEN F = cgetg(3, t_MAT);
971
  gel(F,1) = cgetg(l1, t_COL);
972
  gel(F,2) = cgetg(l1, t_COL);
973
  for (i = 1, j = 1; i < l1; ++i)
974
  {
975
    long v = Z_pval(N, gel(F1, i));
976
    if (v<=1) continue;
977
    gcoeff(F, j  , 1) = gel(F1, i);
978
    gcoeff(F, j++, 2) = stoi(v);
979
  }
980
  setlg(F[1],j); setlg(F[2],j);
981
  return j==1 ? NULL : mkvec2(factorback(F), F);
982
}
983

984
GEN
985
gen_ellgroup(GEN N, GEN d, GEN *pt_m, void *E, const struct bb_group *grp,
986
             GEN pairorder(void *E, GEN P, GEN Q, GEN m, GEN F))
987
{
988
  pari_sp av = avma;
989
  GEN N0, N1, F;
990
  if (pt_m) *pt_m = gen_1;
991
  if (is_pm1(N)) return cgetg(1,t_VEC);
992
  F = ellgroup_d2(N, d);
993
  if (!F) {set_avma(av); return mkveccopy(N);}
994
  N0 = gel(F,1); N1 = diviiexact(N, N0);
995
  while(1)
996
  {
997
    pari_sp av2 = avma;
998
    GEN P, Q, d, s, t, m;
999

1000
    P = grp->pow(E,grp->rand(E), N1);
1001
    s = gen_order(P, F, E, grp); if (equalii(s, N0)) {set_avma(av); return mkveccopy(N);}
1002

1003
    Q = grp->pow(E,grp->rand(E), N1);
1004
    t = gen_order(Q, F, E, grp); if (equalii(t, N0)) {set_avma(av); return mkveccopy(N);}
1005

1006
    m = lcmii(s, t);
1007
    d = pairorder(E, P, Q, m, F);
1008
    /* structure is [N/d, d] iff m d == N0. Note that N/d = N1 m */
1009
    if (is_pm1(d) && equalii(m, N0)) {set_avma(av); return mkveccopy(N);}
1010
    if (equalii(mulii(m, d), N0))
1011
    {
1012
      GEN g = mkvec2(mulii(N1,m), d);
1013
      if (pt_m) *pt_m = m;
1014
      gerepileall(av,pt_m?2:1,&g,pt_m);
1015
      return g;
1016
    }
1017
    set_avma(av2);
1018
  }
1019
}
1020

1021
GEN
1022
gen_ellgens(GEN D1, GEN d2, GEN m, void *E, const struct bb_group *grp,
1023
             GEN pairorder(void *E, GEN P, GEN Q, GEN m, GEN F))
1024
{
1025
  pari_sp ltop = avma, av;
1026
  GEN F, d1, dm;
1027
  GEN P, Q, d, s;
1028
  F = get_arith_ZZM(D1);
1029
  d1 = gel(F, 1), dm =  diviiexact(d1,m);
1030
  av = avma;
1031
  do
1032
  {
1033
    set_avma(av);
1034
    P = grp->rand(E);
1035
    s = gen_order(P, F, E, grp);
1036
  } while (!equalii(s, d1));
1037
  av = avma;
1038
  do
1039
  {
1040
    set_avma(av);
1041
    Q = grp->rand(E);
1042
    d = pairorder(E, grp->pow(E, P, dm), grp->pow(E, Q, dm), m, F);
1043
  } while (!equalii(d, d2));
1044
  return gerepilecopy(ltop, mkvec2(P,Q));
1045
}
1046

1047
Product

Resources

Company
