Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - bb_group.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.1 lcov report (development 30702-bddb8d6928) Lines: 554 589 94.1 %
Date: 2026-02-23 02:23:56 Functions: 37 37 100.0 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.16