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 - pclgp.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27775-aca467eab2) Lines: 1832 2430 75.4 %
Date: 2022-07-03 07:33:15 Functions: 155 187 82.9 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000, 2012  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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : #define DEBUGLEVEL DEBUGLEVEL_subcyclo
      18             : 
      19             : /* written by Takashi Fukuda */
      20             : 
      21             : #define onevec(x) const_vec(x,gen_1)
      22             : #define nullvec() cgetg(1, t_VEC)
      23             : #define order_f_x(f, x) Fl_order(x%f, eulerphiu(f), f)
      24             : 
      25             : #define USE_MLL         (1L<<0)
      26             : #define NO_PLUS_PART    (1L<<1)
      27             : #define NO_MINUS_PART   (1L<<2)
      28             : #define SKIP_PROPER     (1L<<3)
      29             : #define SAVE_MEMORY     (1L<<4)
      30             : #define USE_FULL_EL     (1L<<5)
      31             : #define USE_BASIS       (1L<<6)
      32             : #define USE_FACTOR      (1L<<7)
      33             : #define USE_GALOIS_POL  (1L<<8)
      34             : #define USE_F           (1L<<9)
      35             : 
      36             : static ulong
      37      178002 : _get_d(GEN H) { return umael(H, 2, 1);}
      38             : static ulong
      39      119820 : _get_f(GEN H) { return umael(H, 2, 2);}
      40             : static ulong
      41       30843 : _get_h(GEN H) { return umael(H, 2, 3);}
      42             : static long
      43       28392 : _get_s(GEN H) { return umael(H, 2, 4);}
      44             : static long
      45       53957 : _get_g(GEN H) { return umael(H, 2, 5);}
      46             : static GEN
      47       30787 : _get_H(GEN H) { return gel(H, 3);}
      48             : static ulong
      49      112755 : K_get_d(GEN K) { return _get_d(gel(K,1)); }
      50             : static ulong
      51       82965 : K_get_f(GEN K) { return _get_f(gel(K,1)); }
      52             : static ulong
      53       17424 : K_get_h(GEN K) { return _get_h(gel(K,1)); }
      54             : static long
      55           0 : K_get_s(GEN K) { return _get_s(gel(K,1)); }
      56             : static ulong
      57       17102 : K_get_g(GEN K) { return _get_g(gel(K,1)); }
      58             : static GEN
      59       17368 : K_get_H(GEN K) { return _get_H(gel(K,1)); }
      60             : static ulong
      61       59314 : K_get_dchi(GEN K) { return gel(K,6)[1]; }
      62             : static ulong
      63       23660 : K_get_nconj(GEN K) { return gel(K,6)[2]; }
      64             : 
      65             : /* temporary dummy implementation of factcylo */
      66             : static GEN
      67         308 : Flx_factcyclo(ulong n, ulong p, ulong m, ulong flag)
      68             : {
      69             :   (void) m; (void) flag;
      70         308 :   return gel(Flx_factor(ZX_to_Flx(polcyclo(n,0), p), p), 1);
      71             : }
      72             : 
      73             : static GEN
      74          42 : FpX_factcyclo(ulong n, GEN p, ulong m, ulong flag)
      75             : {
      76             :   (void) m; (void) flag;
      77          42 :   return gel(FpX_factor(polcyclo(n,0), p), 1);
      78             : }
      79             : 
      80             : /* G=<s> is a cyclic group of order n and t=s^(-1).
      81             :  *  convert sum_i a_i*s^i to sum_i b_i*t^i */
      82             : static GEN
      83          14 : Flx_recip1_inplace(GEN x, long pn)
      84             : {
      85          14 :   long i, lx = lg(x);
      86          14 :   if(lx-2 != pn) /* This case scarcely occurs */
      87             :   {
      88           0 :     long ly = pn+2;
      89           0 :     GEN y = const_vecsmall(ly, 0);
      90           0 :     y[1] = x[1];y[2] = x[2];
      91           0 :     for(i=3;i<lx;i++) y[ly+2-i] = x[i];
      92           0 :     return Flx_renormalize(y, ly);
      93             :   }
      94             :   else /* almost all cases */
      95             :   {
      96          14 :     long t, mid = (lx+1)>>1;
      97        7168 :     for(i=3;i<=mid;i++)
      98             :     {
      99        7154 :       t = x[i];x[i] = x[lx+2-i];x[lx+2-i] = t;
     100             :     }
     101          14 :     return Flx_renormalize(x, lx);
     102             :   }
     103             : }
     104             : 
     105             : /* Return h^degpol(P) P(x / h) */
     106             : static GEN
     107          14 : Flx_rescale_inplace(GEN P, ulong h, ulong p)
     108             : {
     109          14 :   long i, l = lg(P);
     110          14 :   ulong hi = h;
     111       14322 :   for (i=l-2; i>=2; i--)
     112             :   {
     113       14322 :     P[i] = Fl_mul(P[i], hi, p);
     114       14322 :     if (i == 2) break;
     115       14308 :     hi = Fl_mul(hi,h, p);
     116             :   }
     117          14 :   return P;
     118             : }
     119             : 
     120             : static GEN
     121          14 : zx_to_Flx_inplace(GEN x, ulong p)
     122             : {
     123          14 :   long i, lx = lg(x);
     124       14350 :   for (i=2; i<lx; i++) uel(x,i) = umodsu(x[i], p);
     125          14 :   return Flx_renormalize(x, lx);
     126             : }
     127             : 
     128             : /* zero pol of n components (i.e. deg=n-1). need to pass to ZX_renormalize */
     129             : INLINE GEN
     130       53263 : pol_zero(long n)
     131             : {
     132             :   long i;
     133       53263 :   GEN p = cgetg(n+2, t_POL);
     134       53263 :   p[1] = evalsigne(1) | evalvarn(0);
     135     1799812 :   for (i = 2; i < n+2; i++) gel(p, i) = gen_0;
     136       53263 :   return p;
     137             : }
     138             : 
     139             : /* e[i+1] = L*i + K for i >= n; determine K,L and reduce n if possible */
     140             : static GEN
     141          21 : vecsmall2vec2(GEN e, long n)
     142             : {
     143          21 :   long L = e[n+1] - e[n], K = e[n+1] - L*n;
     144          42 :   n--; while (n >= 0 && e[n+1] - L*n == K) n--;
     145          21 :   if (n < 0) e = nullvec(); else { setlg(e, n+2); e = zv_to_ZV(e); }
     146          21 :   return mkvec3(utoi(L), stoi(K), e); /* L >= 0 */
     147             : }
     148             : 
     149             : /* z=zeta_{p^n}; return k s.t. (z-1)^k || f(z) assuming deg(f)<phi(p^n) */
     150             : static long
     151          42 : zx_p_val(GEN f, ulong p, ulong n)
     152             : {
     153          42 :   pari_sp av = avma;
     154          42 :   ulong x = zx_lval(f, p);
     155          42 :   if (x) { f = zx_z_divexact(f, upowuu(p, x)); x *= (p-1)*upowuu(p, n-1); }
     156          42 :   x += Flx_val(Flx_translate1(zx_to_Flx(f, p), p));
     157          42 :   return gc_long(av, x);
     158             : }
     159             : 
     160             : static long
     161         315 : ZX_p_val(GEN f, ulong p, ulong n)
     162             : {
     163         315 :   pari_sp av = avma;
     164         315 :   ulong x = ZX_lval(f, p);
     165         315 :   if (x) { f = ZX_Z_divexact(f, powuu(p, x)); x *= (p-1)*upowuu(p, n-1); }
     166         315 :   x += Flx_val(Flx_translate1(ZX_to_Flx(f, p), p));
     167         315 :   return gc_long(av, x);
     168             : }
     169             : 
     170             : static GEN
     171          35 : set_A(GEN B, int *chi)
     172             : {
     173          35 :   long a, i, j, B1 = B[1], l = lg(B);
     174          35 :   GEN A = cgetg(l, t_VECSMALL);
     175     1687014 :   for (a = 0, j = 1; j < B1; j++) a += chi[j];
     176          35 :   A[1] = a;
     177         714 :   for (i = 2; i < l; i++)
     178             :   {
     179         679 :     long Bi = B[i];
     180    28159243 :     for (a = A[i-1], j = B[i-1]; j < Bi; j++) a += chi[j];
     181         679 :     A[i] = a;
     182             :   }
     183          35 :   return A;
     184             : }
     185             : 
     186             : /* g_n(a)=g_n(b) <==> a^2=b^2 mod 2^(n+2) <==> a=b,-b mod 2^(n+2)
     187             :  * g_n(a)=g_n(1+q0)^k <==> a=x(1+q0)^k x=1,-1
     188             :  * gam[1+a]=k, k<0  ==> g_n(a)=0
     189             :  *             k>=0 ==> g_n(a)^(-1)=gamma^k, gamma=g_n(1+q0) */
     190             : static GEN
     191          14 : set_gam2(long q01, long n)
     192             : {
     193             :   long i, x, x1, pn, pn2;
     194             :   GEN gam;
     195          14 :   pn = (1L<<n);
     196          14 :   pn2 = (pn<<2);
     197          14 :   gam = const_vecsmall(pn2, -1);
     198          14 :   x=Fl_inv(q01, pn2); x1=1;
     199       14350 :   for (i=0; i<pn; i++)
     200             :   {
     201       14336 :     gam[1+x1] = gam[1+Fl_neg(x1, pn2)] = i;
     202       14336 :     x1 = Fl_mul(x1, x, pn2);
     203             :   }
     204          14 :   return gam;
     205             : }
     206             : 
     207             : /* g_n(a)=g_n(b) <==> a^(p-1)=b^(p-1) mod p^(n+1) <==> a=xb x=<g^(p^n)>
     208             :  * g_n(a)=g_n(1+q0)^k <==> a=x(1+q0)^k x=<g^(p^n)>
     209             :  * gam[1+a]=k, k<0  ==> g_n(a)=0
     210             :  *             k>=0 ==> g_n(a)^(-1)=gamma^k, gamma=g_n(1+q0) */
     211             : static GEN
     212         476 : set_gam(long q01, long p, long n)
     213             : {
     214             :   long i, j, g, g1, x, x1, p1, pn, pn1;
     215             :   GEN A, gam;
     216         476 :   p1 = p-1; pn = upowuu(p, n); pn1 = p*pn;
     217         476 :   gam = const_vecsmall(pn1, -1);
     218         476 :   g = pgener_Zl(p); g1 = Fl_powu(g, pn, pn1);
     219         476 :   A = Fl_powers(g1, p1-1, pn1);  /* A[1+i]=g^(i*p^n) mod p^(n+1), 0<=i<=p-2 */
     220         476 :   x = Fl_inv(q01, pn1); x1 = 1;
     221      568694 :   for (i=0; i<pn; i++)
     222             :   {
     223     2086126 :     for (j=1; j<=p1; j++) gam[1+Fl_mul(x1, A[j], pn1)] = i;
     224      568218 :     x1 = Fl_mul(x1, x, pn1);
     225             :   }
     226         476 :   return gam;
     227             : }
     228             : 
     229             : /* k=Q(sqrt(m)), A_n=p-class gr. of k_n, |A_n|=p^(e_n)
     230             :  * return e_n-e_(n-1)
     231             :  * essential assumption : m is not divisible by p
     232             :  * Gold, Acta Arith. XXVI (1974), p.22 formula (3) */
     233             : static long
     234          35 : ediff(ulong p, long m, ulong n, int *chi)
     235             : {
     236          35 :   pari_sp av = avma;
     237             :   long j, lx, *px;
     238             :   ulong i, d, s, y, g, p1, pn, pn1, pn_1, phipn, phipn1;
     239             :   GEN A, B, x, gs, cs;
     240             : 
     241          35 :   d=((m-1)%4==0)?labs(m):4*labs(m);
     242          35 :   p1=p-1; pn_1=upowuu(p, n-1); pn=p*pn_1; pn1=p*pn; phipn=p1*pn_1; phipn1=p1*pn;
     243          35 :   lx=2*p1*phipn;
     244          35 :   y=Fl_inv(pn1%d, d); g=pgener_Zl(p);  /* pn1 may > d */
     245          35 :   cs = cgetg(2+phipn, t_VECSMALL); cs[1] = evalvarn(0);
     246          35 :   x = cgetg(1+lx, t_VECSMALL);
     247          35 :   gs = Fl_powers(g, phipn1-1, pn1); /* gs[1+i]=g^i(mod p^(n+1)), 0<=i<p^(n+1) */
     248             : 
     249         105 :   for (px=x,i=0; i<p1; i++)
     250             :   {
     251          70 :     long ipn=i*pn+1,ipnpn=ipn+phipn;
     252         546 :     for (s=0; s<phipn; s++)
     253             :     {
     254         476 :       *++px = (y*gs[s+ipn])%d;  /* gs[s+ipn] may > d */
     255         476 :       *++px = (y*gs[(s%pn_1)+ipnpn])%d;
     256             :     }
     257             :   }
     258          35 :   B = vecsmall_uniq(x);
     259          35 :   A = set_A(B, chi);
     260         273 :   for (s=0; s<phipn; s++)
     261             :   {
     262         238 :     long a=0, ipn=1, spn1=s%pn_1;
     263         714 :     for (i=0; i<p1; i++)
     264             :     {
     265         476 :       if ((j=zv_search(B, (y*gs[s+ipn])%d))<=0)
     266           0 :         pari_err_BUG("zv_search failed\n");
     267         476 :       a+=A[j];
     268         476 :       if ((j=zv_search(B, (y*gs[spn1+ipn+phipn])%d))<=0)
     269           0 :         pari_err_BUG("zv_search failed\n");
     270         476 :       a-=A[j];
     271         476 :       ipn+=pn;
     272             :     }
     273         238 :     cs[2+s] = a;
     274             :   }
     275          35 :   cs = zx_renormalize(cs, lg(cs));
     276          35 :   y = (lg(cs)==3) ? phipn*z_lval(cs[2], p) : zx_p_val(cs, p, n);
     277          35 :   return gc_long(av, y);
     278             : }
     279             : 
     280             : static GEN
     281           0 : quadteichstk(GEN Chi, int *chi, GEN Gam, long p, long m, long n)
     282             : {
     283           0 :   GEN Gam1 = Gam+1, xi;
     284             :   long i, j, j0, d, f0, pn, pn1, deg, pn1d;
     285             : 
     286           0 :   d = ((m&3)==1)?m:m<<2;
     287           0 :   f0 = ulcm(p, d)/p;
     288           0 :   pn = upowuu(p, n); pn1 = p*pn; pn1d = pn1%d;
     289           0 :   xi = cgetg(pn+2, t_POL); xi[1] = evalsigne(1) | evalvarn(0);
     290           0 :   for (i=0; i<pn; i++) gel(xi, 2+i) = const_vecsmall(p, 0);
     291           0 :   for (j=1; j<pn1; j++)
     292             :   {
     293             :     long jp, ipn1d, *xij0;
     294           0 :     if ((j0 = Gam1[j])<0) continue;
     295           0 :     jp = j%p; ipn1d = j%d; xij0 = gel(xi, 2+j0)+2;
     296           0 :     for (i=1; i<f0; i++)
     297             :     {
     298             :       int sgn;
     299           0 :       if ((ipn1d += pn1d) >= d) ipn1d -= d;
     300           0 :       if ((sgn = chi[ipn1d])==0) continue;
     301           0 :       deg = Chi[jp];  /* jp!=0 because j0>=0 */
     302           0 :       if (sgn>0) xij0[deg] += i;
     303           0 :       else xij0[deg] -= i;
     304             :     }
     305             :   }
     306           0 :   for (i=0; i<pn; i++) gel(xi, 2+i) = zx_renormalize(gel(xi, 2+i), p+1);
     307           0 :   return FlxX_renormalize(xi, pn+2);  /* zxX_renormalize does not exist */
     308             : }
     309             : 
     310             : #ifdef DEBUG_QUADSTK
     311             : /* return f0*xi_n */
     312             : static GEN
     313             : quadstkp_by_def(int *chi, GEN gam, long n, long p, long f, long f0)
     314             : {
     315             :   long i, a, a1, pn, pn1, qn;
     316             :   GEN x, x2, gam1 = gam+1;
     317             :   pn = upowuu(p, n); pn1 = p*pn; qn = f0*pn1;
     318             :   x = const_vecsmall(pn+1, 0); x2 = x+2;
     319             :   for (a=1; a<qn; a++)
     320             :   {
     321             :     int sgn;
     322             :     if ((a1=gam1[a%pn1])<0 || (sgn=chi[a%f])==0) continue;
     323             :     if (sgn>0) x2[a1]+=a;
     324             :     else x2[a1]-=a;
     325             :   }
     326             :   for (i=0; i<pn; i++)
     327             :   {
     328             :     if (x2[i]%pn1) pari_err_BUG("stickel. ele. is not integral.\n");
     329             :     else x2[i]/=pn1;
     330             :   }
     331             :   return zx_renormalize(x, pn+2);
     332             : }
     333             : #endif
     334             : 
     335             : /* f!=p
     336             :  * xi_n = f0^(-1)*
     337             :  *   sum_{0<=j<pn1,(j,p)=1}(Q_n/Q,j)^(-1)*(sum_{0<=i<f0}i*chi^(-1)(pn1*i+j)) */
     338             : static GEN
     339          14 : quadstkp1(int *chi, GEN gam, long n, long p, long f, long f0)
     340             : {
     341             :   long i, j, j0, pn, pn1, pn1f, den;
     342             :   GEN x, x2;
     343          14 :   pn = upowuu(p, n); pn1 = p*pn; pn1f = pn1%f;
     344          14 :   x = const_vecsmall(pn+1, 0); x2 = x+2;
     345          14 :   if (f==3) den = (chi[p%f]>0)?f0<<1:2;
     346          14 :   else if (f==4) den = (chi[p%f]>0)?f0<<1:f0;
     347          14 :   else den = f0<<1;
     348     1653372 :   for (j=1; j<pn1; j++)
     349             :   {
     350             :     long ipn1;
     351     1653358 :     if (j%p==0) continue;
     352     1102248 :     j0 = gam[1+j]; ipn1 = j%f;
     353   263437272 :     for (i=1; i<f0; i++)
     354             :     {
     355             :       int sgn;
     356   262335024 :       if ((ipn1+=pn1f)>=f) ipn1-=f;
     357   262335024 :       if ((sgn = chi[ipn1])>0) x2[j0]+=i;
     358   131716319 :       else if (sgn<0) x2[j0]-=i;
     359             :     }
     360             :   }
     361      551138 :   for (i=0; i<pn; i++)
     362             :   {
     363      551124 :     if (x2[i]%den) pari_err_BUG("stickel. ele. is not integral.\n");
     364      551124 :     else x2[i]/=den;
     365             :   }
     366          14 :   return zx_renormalize(x, pn+2);
     367             : }
     368             : 
     369             : /* f==p */
     370             : static GEN
     371           0 : quadstkp2(int *chi, GEN gam, long n, long p)
     372             : {
     373             :   long a, a1, i, pn, pn1, amodp;
     374           0 :   GEN x, x2, gam1 = gam+1;
     375           0 :   pn = upowuu(p, n); pn1 = p*pn;
     376           0 :   x = const_vecsmall(pn+1, 0); x2 = x+2;
     377           0 :   for (a=1,amodp=0; a<pn1; a++)
     378             :   {
     379             :     int sgn;
     380           0 :     if (++amodp==p) {amodp = 0; continue; }
     381           0 :     if ((sgn = chi[amodp])==0) continue;
     382           0 :     a1=gam1[a];
     383           0 :     if (sgn>0) x2[a1]+=a;
     384           0 :     else x2[a1]-=a;
     385             :   }
     386           0 :   for (i=0; i<pn; i++)
     387             :   {
     388           0 :     if (x2[i]%pn1) pari_err_BUG("stickel. ele. is not integral.\n");
     389           0 :     else x2[i]/=pn1;
     390             :   }
     391           0 :   return zx_renormalize(x, pn+2);
     392             : }
     393             : 
     394             : /*  p>=3
     395             :  *  f = conductor of Q(sqrt(m))
     396             :  *  q0 = lcm(f,p) = f0*p
     397             :  *  qn = q0*p^n = f0*p^(n+1)
     398             :  *  xi_n = qn^(-1)*sum_{1<=a<=qn,(a,qn)=1} a*chi(a)^(-1)*(Q_n/Q,a)^(-1) */
     399             : static GEN
     400          14 : quadstkp(long p, long m, long n, int *chi)
     401             : {
     402             :   long f, f0, pn, pn1, q0;
     403             :   GEN gam;
     404          14 :   f = ((m-1)%4==0)?labs(m):4*labs(m);
     405          14 :   pn = upowuu(p, n); pn1 = p*pn;
     406          14 :   if (f % p) { q0 = f * p; f0 = f; } else { q0 = f; f0 = f / p; }
     407          14 :   gam = set_gam((1+q0)%pn1, p, n);
     408             : #ifdef DEBUG_QUADSTK
     409             :   return quadstkp_by_def(chi, gam, n, p, f, f0);
     410             : #else
     411          14 :   return (f0!=1)?quadstkp1(chi, gam, n, p, f, f0):quadstkp2(chi, gam, n, p);
     412             : #endif
     413             : }
     414             : 
     415             : /* p=2 */
     416             : static GEN
     417          14 : quadstk2(long m, long n, int *chi)
     418             : {
     419             :   long i, j, j0, f, f0, pn, pn1, pn2, pn2f, q0;
     420             :   GEN x, x2, gam;
     421          14 :   f = ((m-1)%4==0)?labs(m):4*labs(m);
     422          14 :   pn = 1L<<n; pn1 = pn<<1; pn2 = pn1<<1; pn2f = pn2%f;
     423          14 :   q0 = (f&1)?f*4:f;
     424          14 :   f0 = (f&1)?f:f/4;
     425          14 :   x = const_vecsmall(pn+1, 0); x2 = x+2;
     426          14 :   gam = set_gam2((1+q0)%pn2, n);
     427       57344 :   for (j=1; j<pn2; j++)
     428             :   {
     429             :     long ipn2;
     430       57330 :     if (!(j&1)) continue;
     431       28672 :     j0 = gam[1+j];
     432       28672 :     ipn2 = j%f;
     433             :     /* for (i=1; i<f0; i++) x2[j0]+=i*chi[(i*pn2+j)%f]; */
     434     1691648 :     for (i=1; i<f0; i++)
     435             :     {
     436             :       int sgn;
     437     1662976 :       if ((ipn2+=pn2f)>=f) ipn2-=f;
     438     1662976 :       if ((sgn=chi[ipn2])>0) x2[j0]+=i;
     439      845390 :       else if (sgn<0) x2[j0]-=i;
     440             :     }
     441             :   }
     442       14350 :   for (f0<<=1, i=0; i<pn; i++)
     443             :   {
     444       14336 :     if (x2[i]%f0) pari_err_BUG("stickel. ele. is not integral.\n");
     445       14336 :     else x2[i]/=f0;
     446             :   }
     447          14 :   return zx_renormalize(x, pn+2);
     448             : }
     449             : 
     450             : /* Chin is a generator of the group of the characters of G(Q_n/Q).
     451             :  * chin[1+a]=k, k<0  ==> Chin(a)=0
     452             :  *              k>=0 ==> Chin(a)=zeta_{p^n}^k */
     453             : static GEN
     454          21 : set_chin(long p, long n)
     455             : {
     456          21 :   long i, j, x = 1, g, gpn, pn, pn1;
     457             :   GEN chin, chin1;
     458          21 :   pn = upowuu(p, n); pn1 = p*pn;
     459          21 :   chin = const_vecsmall(pn1, -1); chin1 = chin+1;
     460          21 :   g = pgener_Zl(p); gpn = Fl_powu(g, pn, pn1);
     461         294 :   for (i=0; i<pn; i++)
     462             :   {
     463         273 :     long y = x;
     464         819 :     for (j=1; j<p; j++)
     465             :     {
     466         546 :       chin1[y] = i;
     467         546 :       y = Fl_mul(y, gpn, pn1);
     468             :     }
     469         273 :     x = Fl_mul(x, g, pn1);
     470             :   }
     471          21 :   return chin;
     472             : }
     473             : 
     474             : /* k=Q(sqrt(m)), A_n=p-class gr. of k_n, |A_n|=p^(e_n), p|m
     475             :  * return e_n-e_(n-1).
     476             :  * There is an another method using the Stickelberger element based on
     477             :  * Coates-Lichtenbaum, Ann. Math. vol.98 No.3 (1973), 498-550, Lemma 2.15.
     478             :  * If kro(m,p)!=1, then orders of two groups coincide.
     479             :  * ediff_ber is faster than the Stickelberger element. */
     480             : static long
     481          21 : ediff_ber(ulong p, long m, ulong n, int *chi)
     482             : {
     483          21 :   pari_sp av = avma;
     484             :   long a, d, e, x, y, pn, pn1, qn1;
     485          21 :   GEN B, B2, chin = set_chin(p, n)+1;
     486             : 
     487          21 :   d = ((m-1)%4==0)?labs(m):4*labs(m);
     488          21 :   pn = upowuu(p, n); pn1 = p*pn; qn1 = (d*pn)>>1;
     489          21 :   B = const_vecsmall(pn+1, 0); B2 = B+2;
     490   522105969 :   for (a=x=y=1; a <= qn1; a++) /* x=a%d, y=a%pn1 */
     491             :   {
     492   522105948 :     int sgn = chi[x];
     493   522105948 :     if (sgn)
     494             :     {
     495   172937856 :       long k = chin[y];
     496   172937856 :       if (k >= 0) { if (sgn > 0) B2[k]++; else B2[k]--; }
     497             :     }
     498   522105948 :     if (++x == d) x = 0;
     499   522105948 :     if (++y == pn1) y = 0;
     500             :   }
     501          21 :   B = zx_renormalize(B, pn+2);
     502           7 :   e = (n==1)? zx_p_val(B, p, n)
     503          21 :             : ZX_p_val(ZX_rem(zx_to_ZX(B), polcyclo(pn, 0)), p, n);
     504          21 :   if (p==3 && chi[2] < 0) e--;  /* 2 is a primitive root of 3^n (n>=1) */
     505          21 :   return gc_long(av, e);
     506             : }
     507             : 
     508             : #ifdef DEBUG
     509             : /* slow */
     510             : static int*
     511             : set_quad_chi_1(long m)
     512             : {
     513             :   long a, d, f;
     514             :   int *chi;
     515             :   d=((m-1)%4==0)?m:4*m; f=labs(d);
     516             :   chi= (int*)stack_calloc(sizeof(int)*f);
     517             :   for (a=1; a<f; a++) chi[a]=kross(d, a);
     518             :   return chi;
     519             : }
     520             : #endif
     521             : 
     522             : /* chi[a]=kross(d, a)   0<=a<=f-1
     523             :  * d=discriminant of Q(sqrt(m)), f=abs(d)
     524             :  *
     525             :  * Algorithm: m=-p1*p2*...*pr ==> kross(d,gi)=-1 (1<=i<=r), gi=proot(pi)
     526             :  * set_quad_chi_1(m)=set_quad_chi_2(m) for all square-free m s.t. |m|<10^5. */
     527             : static int*
     528          49 : set_quad_chi_2(long m)
     529             : {
     530          49 :   long d = (m-1) % 4? 4*m: m, f = labs(d);
     531          49 :   GEN fa = factoru(f), P = gel(fa, 1), E = gel(fa,2), u, v;
     532          49 :   long i, j, np, nm, l = lg(P);
     533          49 :   int *chi = (int*)stack_calloc(sizeof(int)*f);
     534          49 :   pari_sp av = avma;
     535          49 :   int *plus = (int*)stack_calloc(sizeof(int)*f), *p0 = plus;
     536          49 :   int *minus = (int*)stack_calloc(sizeof(int)*f), *p1 = minus;
     537             : 
     538          49 :   u = cgetg(32, t_VECSMALL);
     539          49 :   v = cgetg(32, t_VECSMALL);
     540         147 :   for (i = 1; i < l; i++)
     541             :   {
     542          98 :     ulong p = upowuu(P[i], E[i]);
     543          98 :     u[i] = p * Fl_inv(p, f / p);
     544          98 :     v[i] = Fl_sub(1, u[i], f);
     545             :   }
     546          49 :   if (E[1]==2)       /* f=4*(-m) */
     547             :   {
     548          14 :     *p0++ = Fl_add(v[1], u[1], f);
     549          14 :     *p1++ = Fl_add(Fl_mul(3, v[1], f), u[1], f);
     550          14 :     i = 2;
     551             :   }
     552          35 :   else if (E[1]==3)  /* f=8*(-m) */
     553             :   {
     554             :     ulong a;
     555           7 :     *p0++ = Fl_add(v[1], u[1], f);
     556           7 :     a = Fl_add(Fl_mul(3, v[1], f), u[1], f);
     557           7 :     if (kross(d, a) > 0) *p0++ = a; else *p1++ = a;
     558           7 :     a = Fl_add(Fl_mul(5, v[1], f), u[1], f);
     559           7 :     if (kross(d, a) > 0) *p0++ = a; else *p1++ = a;
     560           7 :     a = Fl_add(Fl_mul(7, v[1], f), u[1], f);
     561           7 :     if (kross(d, a) > 0) *p0++ = a; else *p1++ = a;
     562           7 :     i = 2;
     563             :   }
     564             :   else              /* f=-m */
     565          28 :   {*p0++ = 1; i = 1; }
     566         126 :   for (; i < l; i++)
     567             :   {
     568          77 :     ulong gn, g = pgener_Fl(P[i]);
     569          77 :     gn = g = Fl_add(Fl_mul(g, v[i], f), u[i], f);
     570          77 :     np = p0-plus;
     571          77 :     nm = p1-minus;
     572             :     for (;;)
     573             :     {
     574     4802406 :       for (j = 0; j < np; j++) *p1++ = Fl_mul(plus[j], gn, f);
     575     4799655 :       for (j = 0; j < nm; j++) *p0++ = Fl_mul(minus[j], gn, f);
     576        7679 :       gn = Fl_mul(gn, g, f); if (gn == 1) break;
     577     4763745 :       for (j= 0; j < np; j++) *p0++ = Fl_mul(plus[j], gn, f);
     578     4761022 :       for (j = 0; j < nm; j++) *p1++ = Fl_mul(minus[j], gn, f);
     579        7602 :       gn = Fl_mul(gn, g, f); if (gn == 1) break;
     580             :     }
     581             :   }
     582          49 :   np = p0-plus;
     583          49 :   nm = p1-minus;
     584     9548224 :   for (i = 0; i < np; i++) chi[plus[i]] = 1;
     585     9548224 :   for (i = 0; i < nm; i++) chi[minus[i]] = -1;
     586          49 :   set_avma(av); return chi;
     587             : }
     588             : 
     589             : static long
     590        8995 : srh_x(GEN T, long n, long x)
     591             : {
     592       30086 :   for (; x<n; x++) if (!T[x]) return x;
     593         623 :   return -1;
     594             : }
     595             : 
     596             : /* G is a cyclic group of order d. hat(G)=<chi>.
     597             :  * chi, chi^p, ... , chi^(p^(d_chi-1)) are conjugate.
     598             :  * {chi^j | j in C} are repre. of Q_p-congacy classes of inj. chars.
     599             :  *
     600             :  * C is a set of representatives of H/<p>, where H=(Z/dZ)^* */
     601             : static GEN
     602        1134 : set_C(long p, long d, long d_chi, long n_conj)
     603             : {
     604        1134 :   long i, j, x, y, pmodd = p%d;
     605        1134 :   GEN T = const_vecsmall(d, 0)+1;
     606        1134 :   GEN C = cgetg(1+n_conj, t_VECSMALL);
     607        1134 :   if (n_conj==1) { C[1] = 1; return C; }
     608        9618 :   for (i=0, x=1; x >= 0; x = srh_x(T, d, x))
     609             :   {
     610        8995 :     if (cgcd(x, d)==1) C[++i] = x;
     611       40929 :     for (j=0, y=x; j<d_chi; j++) T[y = Fl_mul(y, pmodd, d)] = 1;
     612             :   }
     613         623 :   return C;
     614             : }
     615             : 
     616             : static GEN
     617         343 : FpX_one_cyclo(long n, GEN p)
     618             : {
     619         343 :   if (lgefint(p)==3)
     620         301 :     return Flx_to_ZX(gel(Flx_factcyclo(n, p[2], 1, 0), 1));
     621             :   else
     622          42 :     return gel(FpX_factcyclo(n, p, 1, 0), 1);
     623             : }
     624             : 
     625             : static void
     626       17094 : Flx_red_inplace(GEN x, ulong p)
     627             : {
     628       17094 :   long i, l = lg(x);
     629      274540 :   for (i=2; i<l; i++) x[i] = uel(x, i)%p;
     630       17094 :   Flx_renormalize(x, l);
     631       17094 : }
     632             : 
     633             : /* x[i], T[i] < pn */
     634             : static GEN
     635       39046 : Flxq_xi_conj(GEN x, GEN T, long j, long d, long pn)
     636             : {
     637       39046 :   long i, deg = degpol(x);
     638       39046 :   GEN z = const_vecsmall(d+1, 0);
     639     1032304 :   for (i=0; i<=deg; i++) z[2+Fl_mul(i, j, d)] = x[2+i];
     640       39046 :   return Flx_rem(Flx_renormalize(z, d+2), T, pn);
     641             : }
     642             : 
     643             : static GEN
     644         966 : FlxqX_xi_conj(GEN x, GEN T, long j, long d, long pn)
     645             : {
     646         966 :   long i, l = lg(x);
     647             :   GEN z;
     648         966 :   z = cgetg(l, t_POL); z[1] = evalsigne(1) | evalvarn(0);
     649       40012 :   for (i=2; i<l; i++) gel(z, i) = Flxq_xi_conj(gel(x, i), T, j, d, pn);
     650         966 :   return z;
     651             : }
     652             : 
     653             : static GEN
     654           0 : FlxqX_xi_norm(GEN x, GEN T, long p, long d, long pn)
     655             : {
     656           0 :   long i, d_chi = degpol(T);
     657           0 :   GEN z = x, z1 = x;
     658           0 :   for (i=1; i<d_chi; i++)
     659             :   {
     660           0 :     z1 = FlxqX_xi_conj(z1, T, p, d, pn);
     661           0 :     z = FlxqX_mul(z, z1, T, pn);
     662             :   }
     663           0 :   return z;
     664             : }
     665             : 
     666             : /* assume 0 <= x[i], y[j] <= m-1 */
     667             : static GEN
     668          15 : FpV_shift_add(GEN x, GEN y, GEN m, long start, long end)
     669             : {
     670             :   long i, j;
     671      222300 :   for (i=start, j=1; i<=end; i++, j++)
     672             :   {
     673      222285 :     GEN z = addii(gel(x, i), gel(y, j));
     674      222285 :     gel(x, i) = (cmpii(z, m)>=0)?subii(z, m):z;
     675             :   }
     676          15 :   return x;
     677             : }
     678             : 
     679             : /* assume 0 <= x[i], y[j] <= m-1 */
     680             : static GEN
     681          10 : FpV_shift_sub(GEN x, GEN y, GEN m, long start, long end)
     682             : {
     683             :   long i, j;
     684      112430 :   for (i=start, j=1; i<=end; i++, j++)
     685             :   {
     686      112420 :     GEN z = subii(gel(x, i), gel(y, j));
     687      112420 :     gel(x, i) = (signe(z)<0)?addii(z, m):z;
     688             :   }
     689          10 :   return x;
     690             : }
     691             : 
     692             : /* assume 0 <= x[i], y[j] <= m-1 */
     693             : static GEN
     694         173 : Flv_shift_add(GEN x, GEN y, ulong m, long start, long end)
     695             : {
     696             :   long i, j;
     697     2320113 :   for (i=start, j=1; i<=end; i++, j++)
     698             :   {
     699     2319940 :     ulong xi = x[i], yj = y[j];
     700     2319940 :     x[i] = Fl_add(xi, yj, m);
     701             :   }
     702         173 :   return x;
     703             : }
     704             : 
     705             : /* assume 0 <= x[i], y[j] <= m-1 */
     706             : static GEN
     707         102 : Flv_shift_sub(GEN x, GEN y, ulong m, long start, long end)
     708             : {
     709             :   long i, j;
     710     1165182 :   for (i=start, j=1; i<=end; i++, j++)
     711             :   {
     712     1165080 :     ulong xi = x[i], yj = y[j];
     713     1165080 :     x[i] = Fl_sub(xi, yj, m);
     714             :   }
     715         102 :   return x;
     716             : }
     717             : 
     718             : /* return 0 if p|x. else return 1 */
     719             : INLINE long
     720         896 : Flx_divcheck(GEN x, ulong p)
     721             : {
     722         896 :   long i, l = lg(x);
     723         910 :   for (i=2; i<l; i++) if (uel(x, i)%p) return 1;
     724         448 :   return 0;
     725             : }
     726             : 
     727             : static long
     728         448 : FlxX_weier_deg(GEN pol, long p)
     729             : {
     730         448 :   long i, l = lg(pol);
     731         896 :   for (i=2; i<l && Flx_divcheck(gel(pol, i), p)==0; i++);
     732         448 :   return (i<l)?i-2:-1;
     733             : }
     734             : 
     735             : static long
     736        1582 : Flx_weier_deg(GEN pol, long p)
     737             : {
     738        1582 :   long i, l = lg(pol);
     739        3997 :   for (i=2; i<l && pol[i]%p==0; i++);
     740        1582 :   return (i<l)?i-2:-1;
     741             : }
     742             : 
     743             : static GEN
     744         308 : Flxn_shift_mul(GEN g, long n, GEN p, long d, long m)
     745             : {
     746         308 :   return Flx_shift(Flxn_mul(g, p, d, m), n);
     747             : }
     748             : 
     749             : INLINE long
     750        1057 : deg_trunc(long lam, long p, long n, long pn)
     751             : {
     752             :   long r, x, d;
     753        1260 :   for (r=1,x=p; x<lam; r++) x *= p;  /* r is min int s.t. lam<=p^r */
     754        1057 :   if ((d = (n-r+2)*lam+1)>=pn) d = pn;
     755        1057 :   return d;
     756             : }
     757             : 
     758             : /*  Flx_translate1_basecase(g, pn) becomes slow when degpol(g)>1000.
     759             :  *  So I wrote Flxn_translate1().
     760             :  *  I need lambda to truncate pol.
     761             :  *  But I need to translate T --> 1+T to know lambda.
     762             :  *  Though the code has a little overhead, it is still fast. */
     763             : static GEN
     764         756 : Flxn_translate1(GEN g, long p, long n)
     765             : {
     766             :   long i, j, d, lam, pn, start;
     767         756 :   if (n==1) start = 3;
     768          70 :   else if (n==2) start = 9;
     769          70 :   else start = 10;
     770         756 :   pn = upowuu(p, n);
     771         756 :   for (lam=start; lam; lam<<=1)  /* least upper bound is 3 */
     772             :   {
     773             :     GEN z;
     774         756 :     d = deg_trunc(lam, p, n, pn);
     775         756 :     z = const_vecsmall(d+1, 0);  /* z[2],...,z[d+1] <--> a_0,...,a_{d-1} */
     776       44954 :     for (i=degpol(g); i>=0; i--)
     777             :     {
     778     1683486 :       for (j=d+1; j>2; j--) z[j] = Fl_add(z[j], z[j-1], pn);  /* z = z*(1+T) */
     779       44198 :       z[2] = Fl_add(z[2], g[2+i], pn);
     780             :     }
     781         756 :     z = Flx_renormalize(z, d+2);
     782         756 :     if (Flx_weier_deg(z, p) <= lam) return z;
     783             :   }
     784             :   return NULL; /*LCOV_EXCL_LINE*/
     785             : }
     786             : 
     787             : static GEN
     788         224 : FlxXn_translate1(GEN g, long p, long n)
     789             : {
     790             :   long i, j, d, lam, pn, start;
     791             :   GEN z;
     792         224 :   if (n==1) start = 3;
     793           0 :   else if (n==2) start = 9;
     794           0 :   else start = 10;
     795         224 :   pn = upowuu(p, n);
     796         224 :   for (lam=start; lam; lam<<=1)  /* least upper bound is 3 */
     797             :   {
     798         224 :     d = deg_trunc(lam, p, n, pn);
     799         224 :     z = const_vec(d+1, pol0_Flx(0));  /* z[2],...,z[d+1] <--> a_0,...,a_{d-1} */
     800         224 :     settyp(z, t_POL); z[1] = evalsigne(1) | evalvarn(0);
     801        9408 :     for (i=degpol(g); i>=0; i--)
     802             :     {
     803       64288 :       for (j=d+1; j>2; j--) gel(z, j) = Flx_add(gel(z, j), gel(z, j-1), pn);
     804        9184 :       gel(z, 2) = Flx_add(gel(z, 2), gel(g, 2+i), pn);
     805             :     }
     806         224 :     z = FlxX_renormalize(z, d+2);
     807         224 :     if (FlxX_weier_deg(z, p) <= lam) return z;
     808             :   }
     809             :   return NULL; /*LCOV_EXCL_LINE*/
     810             : }
     811             : 
     812             : /* lam < 0 => error (lambda can't be determined)
     813             :  * lam = 0 => return 1
     814             :  * lam > 0 => return dist. poly. of degree lam. */
     815             : static GEN
     816          84 : Flxn_Weierstrass_prep(GEN g, long p, long n, long d_chi)
     817             : {
     818          84 :   long i, r0, d, dg = degpol(g), lam, pn, t;
     819             :   ulong lam0;
     820             :   GEN U, UINV, P, PU, g0, g1, gp, gU;
     821          84 :   if ((lam = Flx_weier_deg(g, p))==0) return(pol1_Flx(0));
     822          77 :   else if (lam<0)
     823           0 :     pari_err(e_MISC,"Flxn_Weierstrass_prep: precision too low. Increase n!");
     824          77 :   lam0 = lam/d_chi;
     825          77 :   pn = upowuu(p, n);
     826          77 :   d = deg_trunc(lam, p, n, pn);
     827          77 :   if (d>dg) d = dg;
     828          77 :   if (d<=lam) d=1+lam;
     829         140 :   for (r0=1; upowuu(p, r0)<lam0; r0++);
     830          77 :   g = Flxn_red(g, d);
     831          77 :   t = Fl_inv(g[2+lam], pn);
     832          77 :   g = Flx_Fl_mul(g, t, pn);  /* normalized so as g[2+lam]=1 */
     833          77 :   U = Flx_shift(g, -lam);
     834          77 :   UINV = Flxn_inv(U, d, pn);
     835          77 :   P = zx_z_divexact(Flxn_red(g, lam), p);  /* assume g[i] <= LONG_MAX */
     836          77 :   PU = Flxn_mul(P, UINV, d, pn);
     837          77 :   gU = Flxn_mul(g, UINV, d, pn);
     838          77 :   g0 = pol1_Flx(0);
     839          77 :   g1 = pol1_Flx(0);
     840         385 :   for (i=1; i<n; i++)
     841             :   {
     842         308 :     g1 = Flxn_shift_mul(g1, -lam, PU, d, pn);
     843         308 :     gp = Flx_Fl_mul(g1, upowuu(p, i), pn);
     844         308 :     g0 = (i&1)?Flx_sub(g0, gp, pn):Flx_add(g0, gp, pn);
     845             :   }
     846          77 :   g0 = Flxn_mul(g0, gU, lam+1, pn);
     847          77 :   g0 = Flx_red(g0, upowuu(p, (p==2)?n-r0:n+1-r0));
     848          77 :   return g0;
     849             : }
     850             : 
     851             : /* xi_n and Iwasawa pol. for Q(sqrt(m)) and p
     852             :  *
     853             :  * (flag&1)!=0 ==> output xi_n
     854             :  * (flag&2)!=0 ==> output power series
     855             :  * (flag&4)!=0 ==> output Iwasawa polynomial */
     856             : static GEN
     857          14 : imagquadstkpol(long p, long m, long n)
     858             : {
     859          14 :   long pn = upowuu(p, n);
     860             :   GEN pol, stk, stk2;
     861             :   int *chi;
     862          14 :   if (p==2 && (m==-1 || m==-2 || m==-3 || m==-6)) return nullvec();
     863          14 :   if (p==3 && m==-3) return nullvec();
     864          14 :   if (p==2 && m%2==0) m /= 2;
     865          14 :   chi = set_quad_chi_2(m);
     866          14 :   stk = (p==2)? quadstk2(m, n, chi): quadstkp(p, m, n, chi);
     867          14 :   stk2 = zx_to_Flx(stk, pn);
     868          14 :   pol = Flxn_Weierstrass_prep(zlx_translate1(stk2, p, n), p, n, 1);
     869          14 :   return degpol(pol)? mkvec(Flx_to_ZX(pol)): nullvec();
     870             : }
     871             : 
     872             : /* a mod p == g^i mod p ==> omega(a)=zeta_(p-1)^(-i)
     873             :  *  Chi[g^i mod p]=i (0 <= i <= p-2) */
     874             : static GEN
     875           0 : get_teich(long p, long g)
     876             : {
     877           0 :   long i, gi = 1, p1 = p-1;
     878           0 :   GEN Chi = cgetg(p, t_VECSMALL);
     879           0 :   for (i=0; i<p1; i++) { Chi[gi] = i; gi = Fl_mul(gi, g, p); }
     880           0 :   return Chi;
     881             : }
     882             : 
     883             : /* Ichimura-Sumida criterion for Greenberg conjecture for real quadratic field.
     884             :  * chi: character of Q(sqrt(m)), omega: Teichmuller character mod p or 4.
     885             :  * Get Stickelberger element from chi^* = omega*chi^(-1) and convert it to
     886             :  * power series by the correspondence (Q_n/Q,1+q0)^(-1) <-> (1+T)(1+q0)^(-1) */
     887             : static GEN
     888          14 : realquadstkpol(long p, long m, long n)
     889             : {
     890             :   int *chi;
     891          14 :   long pnm1 = upowuu(p, n-1),pn = p*pnm1, pn1 = p*pn, d, q0;
     892             :   GEN stk, ser, pol;
     893          14 :   if (m==1) pari_err_DOMAIN("quadstkpol", "m", "=", gen_1, gen_1);
     894          14 :   if (p==2 && (m&1)==0) m>>=1;
     895          14 :   d = ((m&3)==1)?m:m<<2;
     896          14 :   q0 = ulcm((p==2)?4:p, d);
     897          14 :   if (p==2)
     898             :   {
     899          14 :     chi = set_quad_chi_2(-m);
     900          14 :     stk = quadstk2(-m, n, chi);
     901          14 :     stk = zx_to_Flx_inplace(stk, pn);
     902             :   }
     903           0 :   else if (p==3 && m%3==0 && kross(-m/3,3)==1)
     904           0 :   {
     905           0 :     long m3 = m/3;
     906           0 :     chi = set_quad_chi_2(-m3);
     907           0 :     stk = quadstkp(3, -m3, n, chi);
     908           0 :     stk = zx_to_Flx_inplace(stk, pn);
     909             :   }
     910             :   else
     911             :   {
     912           0 :     long g = pgener_Zl(p);
     913           0 :     long x = Fl_powu(Fl_inv(g, p), pnm1, pn);
     914           0 :     GEN Chi = get_teich(p, g);
     915           0 :     GEN Gam = set_gam((1+q0)%pn1, p, n);
     916           0 :     chi = set_quad_chi_2(m);
     917           0 :     stk = quadteichstk(Chi, chi, Gam, p, m, n);  /* exact */
     918           0 :     stk = zxX_to_FlxX(stk, pn);  /* approx. */
     919           0 :     stk = FlxY_evalx(stk, x, pn);
     920             :   }
     921          14 :   stk = Flx_rescale_inplace(Flx_recip1_inplace(stk, pn), (1+q0)%pn, pn);
     922          14 :   ser = Flxn_translate1(stk, p, n);
     923          14 :   pol = Flxn_Weierstrass_prep(ser, p, n, 1);
     924          14 :   return degpol(pol)? mkvec(Flx_to_ZX(pol)): nullvec();
     925             : }
     926             : 
     927             : /* m > 0 square-free. lambda_2(Q(sqrt(-m)))
     928             :  * Kida, Tohoku Math. J. vol.31 (1979), 91-96, Theorem 1. */
     929             : static GEN
     930           0 : quadlambda2(long m)
     931             : {
     932             :   long i, l, L;
     933             :   GEN P;
     934           0 :   if ((m&1)==0) m >>= 1;  /* lam_2(Q(sqrt(-m)))=lam_2(Q(sqrt(-2*m))) */
     935           0 :   if (m <= 3) return mkvecs(0);
     936           0 :   P = gel(factoru(m), 1); l = lg(P);
     937           0 :   for (L = -1,i = 1; i < l; i++) L += 1L << (-3 + vals(P[i]-1) + vals(P[i]+1));
     938           0 :   return mkvecs(L);
     939             : }
     940             : 
     941             : /* Iwasawa lambda invariant of Q(sqrt(m)) (m<0) for p
     942             :  * |A_n|=p^(e[n])
     943             :  * kross(m,p)!=1 : e[n]-e[n-1]<eulerphi(p^n)  ==> lambda=e[n]-e[n-1]
     944             :  * kross(m,p)==1 : e[n]-e[n-1]<=eulerphi(p^n) ==> lambda=e[n]-e[n-1]
     945             :  * Gold, Acta Arith. XXVI (1974), p.25, Cor. 3
     946             :  * Gold, Acta Arith. XXVI (1975), p.237, Cor. */
     947             : static GEN
     948          21 : quadlambda(long p, long m)
     949             : {
     950             :   long flag, n, phipn;
     951          21 :   GEN e = cgetg(31, t_VECSMALL);
     952             :   int *chi;
     953          21 :   if (m>0) pari_err_IMPL("plus part of lambda invariant in quadlambda()");
     954          21 :   if (p==2) return quadlambda2(-m);
     955          21 :   if (p==3 && m==-3) return mkvec3(gen_0, gen_0, nullvec());
     956          21 :   flag = kross(m, p);
     957          21 :   e[1] = Z_lval(quadclassno(quaddisc(stoi(m))), p);
     958          21 :   if (flag!=1 && e[1]==0) return mkvec3(gen_0, gen_0, nullvec());
     959          21 :   chi = set_quad_chi_2(m);
     960          21 :   phipn = p-1;  /* phipn=phi(p^n) */
     961          56 :   for (n=1; n; n++, phipn *= p)
     962             :   {
     963          56 :     long L = flag? ediff(p, m, n, chi): ediff_ber(p, m, n, chi);
     964          56 :     e[n+1] = e[n] + L;
     965          56 :     if ((flag!=1 && (L < phipn))|| (flag==1 && (L <= phipn))) break;
     966             :   }
     967          21 :   return vecsmall2vec2(e, n);
     968             : }
     969             : 
     970             : /* factor n-th cyclotomic polynomial mod p^r and return a minimal
     971             :  *  polynomial of zeta_n over Q_p.
     972             :  *  phi(n)=deg*n_conj, n_conj == 1 <=> polcyclo(n) is irred mod p. */
     973             : static GEN
     974         945 : set_minpol(ulong n, GEN p, ulong r, long n_conj)
     975             : {
     976             :   GEN z, v, pol, pr;
     977             :   pari_timer ti;
     978         945 :   if (umodiu(p, n)==1) /* zeta_n in Z_p, faster than polcyclo() */
     979             :   {
     980         420 :     GEN prm1 = powiu(p, r-1), pr = mulii(prm1, p); /* pr=p^r */
     981         420 :     GEN prn = diviuexact(subii(pr, prm1), n);      /* prn=phi(p^r)/n */
     982         420 :     z = Fp_pow(pgener_Fp(p), prn, pr);
     983         420 :     return deg1pol_shallow(gen_1, Fp_neg(z, pr), 0);
     984             :   }
     985         525 :   pr = powiu(p, r);
     986         525 :   pol = polcyclo(n, 0);
     987         525 :   if (n_conj==1) return FpX_red(pol, pr);
     988         343 :   if (DEBUGLEVEL>3) timer_start(&ti);
     989         343 :   z = FpX_one_cyclo(n, p);
     990         343 :   if (DEBUGLEVEL>3) timer_printf(&ti, "FpX_one_cyclo:n=%ld  ", n);
     991         343 :   v = ZpX_liftfact(pol, mkvec2(z, FpX_div(pol, z, p)), pr, p, r);
     992         343 :   return gel(v, 1);
     993             : }
     994             : 
     995             : static GEN
     996          91 : set_minpol_teich(ulong g_K, GEN p, ulong r)
     997             : {
     998          91 :   GEN prm1 = powiu(p, r-1), pr = mulii(prm1, p), z;
     999          91 :   z = Fp_pow(Fp_inv(utoi(g_K), p), prm1, pr);
    1000          91 :   return deg1pol_shallow(gen_1, Fp_neg(z, pr), 0);
    1001             : }
    1002             : 
    1003             : static long
    1004       18963 : srh_1(GEN H)
    1005             : {
    1006       18963 :   GEN bits = gel(H, 3);
    1007       18963 :   ulong f = bits[1];
    1008       18963 :   return F2v_coeff(bits, f-1);
    1009             : }
    1010             : 
    1011             : /* (1/f)sum_{1<=a<=f}a*chi^{-1}(a) = -(1/(2-chi(a)))sum_{1<=a<=f/2} chi^{-1}(a)
    1012             :  *  does not overflow */
    1013             : static GEN
    1014       13146 : zx_ber_num(GEN Chi, long f, long d)
    1015             : {
    1016       13146 :   long i, f2 = f>>1;
    1017       13146 :   GEN x = const_vecsmall(d+1, 0), x2 = x+2;
    1018    51965081 :   for (i = 1; i <= f2; i++)
    1019    51951935 :     if (Chi[i] >= 0) x2[Chi[i]] ++;
    1020       13146 :   return zx_renormalize(x, d+2);
    1021             : }
    1022             : 
    1023             : /* x a zx
    1024             :  * zx_ber_num is O(f). ZX[FpX,Flx]_ber_conj is O(d). Sometimes d<<f. */
    1025             : static GEN
    1026       26257 : ZX_ber_conj(GEN x, long j, long d)
    1027             : {
    1028       26257 :   long i, deg = degpol(x);
    1029       26257 :   GEN y = pol_zero(d), x2 = x+2, y2 = y+2;
    1030      818202 :   for (i=0; i<=deg; i++) gel(y2, Fl_mul(i, j, d)) = stoi(x2[i]);
    1031       26257 :   return ZX_renormalize(y, d+2);
    1032             : }
    1033             : 
    1034             : /* x a zx */
    1035             : static GEN
    1036         252 : FpX_ber_conj(GEN x, long j, long d, GEN p)
    1037             : {
    1038         252 :   long i, deg = degpol(x);
    1039         252 :   GEN y = pol_zero(d), x2 = x+2, y2 = y+2;
    1040         756 :   for (i=0; i<=deg; i++) gel(y2, Fl_mul(i, j, d)) = modsi(x2[i], p);
    1041         252 :   return FpX_renormalize(y, d+2);
    1042             : }
    1043             : 
    1044             : /* x a zx */
    1045             : static GEN
    1046       21756 : Flx_ber_conj(GEN x, long j, long d, ulong p)
    1047             : {
    1048       21756 :   long i, deg = degpol(x);
    1049       21756 :   GEN y = const_vecsmall(d+1, 0), x2 = x+2, y2 = y+2;
    1050     1076565 :   for (i=0; i<=deg; i++) y2[Fl_mul(i, j, d)] = umodsu(x2[i], p);
    1051       21756 :   return Flx_renormalize(y, d+2);
    1052             : }
    1053             : 
    1054             : static GEN
    1055       26257 : ZX_ber_den(GEN Chi, long j, long d)
    1056             : {
    1057       26257 :   GEN x = pol_zero(d), x2 = x+2;
    1058       26257 :   if (Chi[2]>=0) gel(x2, Fl_neg(Fl_mul(Chi[2], j, d), d)) = gen_1;
    1059       26257 :   gel(x2, 0) = subiu(gel(x2, 0), 2);
    1060       26257 :   return ZX_renormalize(x, d+2);
    1061             : }
    1062             : 
    1063             : static GEN
    1064       14490 : Flx_ber_den(GEN Chi, long j, long d, ulong p)
    1065             : {
    1066       14490 :   GEN x = const_vecsmall(d+1, 0), x2 = x+2;
    1067       14490 :   if (Chi[2]>=0) x2[Fl_neg(Fl_mul(Chi[2], j, d), d)] = 1;
    1068       14490 :   x2[0] = Fl_sub(x2[0], 2, p);
    1069       14490 :   return Flx_renormalize(x, d+2);
    1070             : }
    1071             : 
    1072             : /* x is ZX of deg <= d-1 */
    1073             : static GEN
    1074         196 : ber_conj(GEN x, long k, long d)
    1075             : {
    1076         196 :   long i, deg = degpol(x);
    1077         196 :   GEN z = pol_zero(d);
    1078         196 :   if (k==1)
    1079           0 :     for (i=0; i<=deg; i++) gel(z, 2+i) = gel(x, 2+i);
    1080             :   else
    1081      122206 :     for (i=0; i<=deg; i++) gel(z, 2+Fl_mul(i, k, d)) = gel(x, 2+i);
    1082         196 :   return ZX_renormalize(z, d+2);
    1083             : }
    1084             : 
    1085             : /* The computation is fast when p^n and el=1+k*f*p^n are less than 2^64
    1086             :  *  for m <= n <= M
    1087             :  *  We believe M>=3 is enough when f%p=0 and M>=2 is enough for other case
    1088             :  *  because we expect that p^2 does not divide |A_{K,psi}| for a large p.
    1089             :  *  FIXME: M should be set according to p and f. */
    1090             : static void
    1091         196 : set_p_f(GEN pp, ulong f, long *pm, long *pM)
    1092             : {
    1093         196 :   ulong p = itou_or_0(pp);
    1094         196 :   if (!p || p >= 2000000) { *pm=2; *pM = dvdui(f, pp)? 3: 2; }
    1095         189 :   else if (p == 3)      { *pm=5; *pM=20; }
    1096         119 :   else if (p == 5)      { *pm=5; *pM=13; }
    1097          56 :   else if (p == 7)      { *pm=5; *pM=11; }
    1098          42 :   else if (p == 11)     { *pm=5; *pM=9; }
    1099          35 :   else if (p == 13)     { *pm=5; *pM=8; }
    1100          28 :   else if (p < 400)     { *pm=5; *pM=7; }
    1101           0 :   else if (p < 5000)    { *pm=3; *pM=5; }
    1102           0 :   else if (p < 50000)   { *pm=2; *pM=4; }
    1103           0 :   else                  { *pm=2; *pM=3; }
    1104         196 : }
    1105             : 
    1106             : static GEN
    1107       18795 : subgp2ary(GEN H, long n)
    1108             : {
    1109       18795 :   GEN v = gel(H, 3), w = cgetg(n+1, t_VECSMALL);
    1110       18795 :   long i, j, f = v[1];
    1111   399982464 :   for (i = 1, j = 0; i <= f; i++)
    1112   399963669 :     if (F2v_coeff(v,i)) w[++j] = i;
    1113       18795 :   return w;
    1114             : }
    1115             : 
    1116             : static GEN
    1117       19124 : Flv_FlvV_factorback(GEN g, GEN x, ulong q)
    1118       90216 : { pari_APPLY_ulong(Flv_factorback(g, gel(x,i), q)) }
    1119             : 
    1120             : /* lift chi character on G/H to character on G */
    1121             : static GEN
    1122       18795 : zncharlift(GEN chi, GEN ncycGH, GEN U, GEN cycG)
    1123             : {
    1124       18795 :   GEN nchi = char_normalize(chi, ncycGH);
    1125       18795 :   GEN c = ZV_ZM_mul(gel(nchi, 2), U), d = gel(nchi, 1);
    1126       18795 :   return char_denormalize(cycG, d, c);
    1127             : }
    1128             : 
    1129             : /* 0 <= c[i] < d, i=1..r; (c[1],...,c[r], d) = 1; find e[i] such that
    1130             :  * sum e[i]*c[i] = 1 mod d */
    1131             : static GEN
    1132       18795 : Flv_extgcd(GEN c, ulong d)
    1133             : {
    1134       18795 :   long i, j, u, f, l = lg(c);
    1135       18795 :   GEN e = zero_zv(l-1);
    1136       18795 :   if (l == 1) return e;
    1137       46053 :   for (f = d, i = 1; f != 1 && i < l; i++)
    1138             :   {
    1139       27258 :     f = cbezout(f, itou(gel(c,i)), &u, &e[i]);
    1140       27258 :     if (!e[i]) continue;
    1141       25004 :     e[i] = umodsu(e[i], d);
    1142       25004 :     u = umodsu(u, d);
    1143       32998 :     if (u != 1) for (j = 1; j < i; j++) e[j] = Fl_mul(e[j], u, d);
    1144             :   }
    1145       18795 :   return e;
    1146             : }
    1147             : 
    1148             : /* f!=p; return exact xi. */
    1149             : static GEN
    1150         462 : get_xi_1(GEN Chi, GEN Gam, long p, long f, long n, long d, ulong pm)
    1151             : {
    1152         462 :   GEN Gam1 = Gam+1, xi;
    1153             :   long i, j, j0, f0, pn, pn1, deg, pn1f;
    1154             : 
    1155         462 :   f0 = (f%p)?f:f/p;
    1156         462 :   pn = upowuu(p, n); pn1 = p*pn; pn1f = pn1%f;
    1157         462 :   xi = cgetg(pn+2, t_POL); xi[1] = evalsigne(1) | evalvarn(0);
    1158       17556 :   for (i=0; i<pn; i++) gel(xi, 2+i) = const_vecsmall(d+1, 0);
    1159      432754 :   for (j=1; j<pn1; j++)
    1160             :   {
    1161             :     long ipn1,*xij0;
    1162      432292 :     if ((j0 = Gam1[j])<0) continue;
    1163      415660 :     ipn1 = j%f; xij0 = gel(xi, 2+j0)+2;
    1164  4027401588 :     for (i=1; i<f0; i++)
    1165             :     {
    1166  4026985928 :       if ((ipn1 += pn1f) >= f) ipn1 -= f;
    1167  4026985928 :       if (ipn1==0 || (deg = Chi[ipn1])<0) continue;
    1168  2203029612 :       xij0[deg] += i;
    1169             :     }
    1170             :   }
    1171       17556 :   for (i=0; i<pn; i++) Flx_red_inplace(gel(xi, 2+i), pm);
    1172         462 :   return FlxX_renormalize(xi, pn+2);
    1173             : }
    1174             : 
    1175             : /* f=p; return p^(n+1)*xi mod pm. */
    1176             : static GEN
    1177           0 : get_xi_2(GEN Chi, GEN Gam, long p, long f, long n, long d, ulong pm)
    1178             : {
    1179             :   long a, amodf, i, j0, pn, pn1, deg;
    1180           0 :   GEN Gam1 = Gam+1, xi;
    1181             : 
    1182           0 :   pn = upowuu(p, n); pn1 = p*pn;
    1183           0 :   xi = cgetg(pn+2, t_POL); xi[1] = evalsigne(1) | evalvarn(0);
    1184           0 :   for (i=0; i<pn; i++) gel(xi, 2+i) = const_vecsmall(d+1, 0);
    1185           0 :   for (a=1,amodf=0; a<pn1; a++)  /* xi is exact */
    1186             :   {
    1187           0 :     if (++amodf==f) amodf = 0;
    1188           0 :     if ((j0=Gam1[a])<0 || amodf==0 || (deg=Chi[amodf])<0) continue;
    1189           0 :     mael(xi, 2+j0, 2+deg) += a;
    1190             :   }
    1191           0 :   for (i=0; i<pn; i++) Flx_red_inplace(gel(xi, 2+i), pm);
    1192           0 :   return FlxX_renormalize(xi, pn+2);
    1193             : }
    1194             : 
    1195             : static GEN
    1196          56 : pol_chi_xi(GEN K, long p, long j, long n)
    1197             : {
    1198          56 :   pari_sp av = avma;
    1199          56 :   GEN MinPol2 = gel(K, 7), xi = gel(K, 8);
    1200          56 :   long d = K_get_d(K), f = K_get_f(K), d_chi = K_get_dchi(K);
    1201          56 :   long wd, minpolpow = (f==p)?2*n+1:n, pm = upowuu(p, minpolpow);
    1202             :   GEN ser, pol, xi_conj;
    1203             :   pari_timer ti;
    1204             : 
    1205             :   /* xi is FlxX mod p^m, MinPol2 is Flx mod p^m, xi_conj is FlxqX. */
    1206          56 :   xi_conj = FlxqX_xi_conj(xi, MinPol2, j, d, pm);
    1207          56 :   if (d_chi==1)  /* d_chi==1 if f==p */
    1208             :   {
    1209          56 :     xi_conj = FlxX_to_Flx(xi_conj);
    1210          56 :     if (f==p) xi_conj = zx_z_divexact(xi_conj, upowuu(p, n+1));
    1211             :   }
    1212             :   /* Now xi_conj is mod p^n */
    1213          56 :   if (DEBUGLEVEL>1) timer_start(&ti);
    1214          56 :   ser = (d_chi==1) ? Flxn_translate1(xi_conj, p, n)
    1215          56 :     : FlxXn_translate1(xi_conj, p, n);
    1216          56 :   if (DEBUGLEVEL>1) timer_printf(&ti, "Flx%sn_translate1",(d_chi==1)?"":"X");
    1217          56 :   wd = (d_chi==1)?Flx_weier_deg(ser, p):FlxX_weier_deg(ser, p);
    1218          56 :   if (wd<0) pari_err(e_MISC,"pol_chi_xi: precision too low. Increase n!\n");
    1219          56 :   else if (wd==0) return pol_1(0);
    1220             :   /* wd>0, convert to dist. poly. */
    1221          56 :   if (d_chi>1)  /* f!=p. minpolpow==n */
    1222             :   {
    1223           0 :     ser = FlxqX_xi_norm(ser, MinPol2, p, d, upowuu(p, n));
    1224           0 :     ser = FlxX_to_Flx(ser);
    1225             :   }
    1226          56 :   pol = Flx_to_ZX(Flxn_Weierstrass_prep(ser, p, n, d_chi));
    1227          56 :   setvarn(pol, fetch_user_var("T"));
    1228             : #ifdef DEBUG
    1229             :   if (wd>0 && d_chi>1)
    1230             :     err_printf("(wd,d_chi,p,f,d,j,H)=(%ld,%ld,%ld,%ld,%ld,%ld,%Ps)\n",
    1231             :         wd,d_chi,p,f,d,j,gmael3(K, 1, 1, 1));
    1232             : #endif
    1233          56 :   return gerepilecopy(av, pol);
    1234             : }
    1235             : 
    1236             : /* return 0 if lam_psi (psi=chi^j) is determined to be zero.
    1237             :  * else return -1.
    1238             :  * If psi(p)!=1, then N_{Q(zeta_d)/Q}(1-psi(p))!=0 (mod p) */
    1239             : static long
    1240       14504 : lam_chi_ber(GEN K, long p, long j)
    1241             : {
    1242       14504 :   pari_sp av = avma;
    1243       14504 :   GEN B1, B2, Chi = gel(K, 2), MinPol2 = gel(K, 7), B_num = gel(K, 8);
    1244       14504 :   long x, p2 = p*p, d = K_get_d(K), f = K_get_f(K);
    1245             : 
    1246       14504 :   if (f == d+1 && p == f && j == 1) return 0;  /* Teichmuller */
    1247             : 
    1248       14490 :   B1 = Flx_rem(Flx_ber_conj(B_num, j, d, p2), MinPol2, p2);
    1249       14490 :   B2 = Flx_rem(Flx_ber_den(Chi, j, d, p2), MinPol2, p2);
    1250       14490 :   if (degpol(B1)<0 || degpol(B2)<0)
    1251          49 :     return gc_long(av, -1); /* 0 mod p^2 */
    1252       14441 :   x = zx_lval(B1, p) - zx_lval(B2, p);
    1253       14441 :   if (x<0) pari_err_BUG("subcycloiwasawa [Bernoulli number]");
    1254       14441 :   return gc_long(av, x==0 ? 0: -1);
    1255             : }
    1256             : 
    1257             : static long
    1258         910 : lam_chi_xi(GEN K, long p, long j, long n)
    1259             : {
    1260         910 :   pari_sp av = avma;
    1261         910 :   GEN xi_conj, z, MinPol2 = gel(K, 7), xi = gel(K, 8);
    1262         910 :   long d = K_get_d(K), f = K_get_f(K), d_chi = K_get_dchi(K);
    1263         910 :   long wd, minpolpow = (f==p)?n+2:1, pm = upowuu(p, minpolpow);
    1264             : 
    1265             :   /* xi is FlxX mod p^m, MinPol2 is Flx mod p^m, xi_conj is FlxqX. */
    1266         910 :   xi_conj = FlxqX_xi_conj(xi, MinPol2, j, d, pm);
    1267         910 :   if (d_chi==1)  /* d_chi==1 if f==p */
    1268             :   {
    1269         686 :     xi_conj = FlxX_to_Flx(xi_conj);
    1270         686 :     if (f==p) xi_conj = zx_z_divexact(xi_conj, upowuu(p, n+1));
    1271             :   }
    1272             :   /* Now xi_conj is mod p^n */
    1273         686 :   z = (d_chi==1) ? Flxn_translate1(xi_conj, p, n)
    1274         910 :     : FlxXn_translate1(xi_conj, p, n);
    1275         910 :   wd = (d_chi==1)?Flx_weier_deg(z, p):FlxX_weier_deg(z, p);
    1276             : #ifdef DEBUG
    1277             :   if (wd>0 && d_chi>1)
    1278             :     err_printf("(wd,d_chi,p,f,d,j,H)=(%ld,%ld,%ld,%ld,%ld,%ld,%Ps)\n",
    1279             :         wd,d_chi,p,f,d,j,gmael3(K, 1, 1, 1));
    1280             : #endif
    1281         910 :   return gc_long(av, wd<0 ? -1 : wd*d_chi);
    1282             : }
    1283             : 
    1284             : /* K = [H1, Chi, Minpol, C, [d_chi, n_conj]] */
    1285             : static GEN
    1286          56 : imag_cyc_pol(GEN K, long p, long n)
    1287             : {
    1288          56 :   pari_sp av = avma;
    1289          56 :   GEN Chi = gel(K, 2), MinPol = gel(K, 3), C = gel(K, 4), MinPol2;
    1290          56 :   long d_K = K_get_d(K), f = K_get_f(K), n_conj = K_get_nconj(K);
    1291          56 :   long i, q0, pn1, pM, pmodf = p%f, n_done = 0;
    1292          56 :   GEN z = nullvec(), Gam, xi, Lam, K2;
    1293             : 
    1294          56 :   Lam = const_vecsmall(n_conj, -1);
    1295          56 :   if (pmodf==0 || Chi[pmodf]) /* mark trivial chi-part using Bernoulli number */
    1296             :   {
    1297          14 :     MinPol2 = ZX_to_Flx(MinPol, p*p);  /* p^2 for B_{1,chi} */
    1298          14 :     K2 = shallowconcat(K, mkvec2(MinPol2, zx_ber_num(Chi, f, d_K)));
    1299          42 :     for (i=1; i<=n_conj; i++)
    1300          28 :       if ((Lam[i] = lam_chi_ber(K2, p, C[i])) == 0) n_done++;
    1301          14 :     if (n_conj==n_done) return gerepilecopy(av, z); /* all chi-parts trivial */
    1302             :   }
    1303          49 :   q0 = (f%p)? f*p: f;
    1304          49 :   pn1 = upowuu(p, n+1);
    1305          49 :   Gam = set_gam((1+q0)%pn1, p, n);
    1306          49 :   pM = upowuu(p, (f==p)? 2*n+1: n);
    1307          49 :   MinPol2 = ZX_to_Flx(MinPol, pM);
    1308           0 :   xi = (f==p)? get_xi_2(Chi, Gam, p, f, n, d_K, pM)
    1309          49 :              : get_xi_1(Chi, Gam, p, f, n, d_K, pM);
    1310          49 :   K2 = shallowconcat(K, mkvec2(MinPol2, xi));
    1311         105 :   for (i=1; i<=n_conj; i++)
    1312             :   {
    1313             :     GEN z1;
    1314          56 :     if (Lam[i]>=0) continue;
    1315          56 :     z1 = pol_chi_xi(K2, p, C[i], n);
    1316          56 :     if (degpol(z1)) z = vec_append(z, z1);  /* degpol(z1) may be zero */
    1317             :   }
    1318          49 :   return gerepilecopy(av, z);
    1319             : }
    1320             : 
    1321             : /* K is an imaginary cyclic extension of degree d contained in Q(zeta_f)
    1322             :  * H is the subgr of G=(Z/fZ)^* corresponding to K
    1323             :  * h=|H|, d*h=phi(f)
    1324             :  * G/H=<g> i.e. g^d \in H
    1325             :  * d_chi=[Q_p(zeta_d):Q_p], i.e. p^d_chi=1 (mod d)
    1326             :  * An inj. char. of G(K/Q) is automatically imaginary.
    1327             :  *
    1328             :  * G(K/Q)=G/H=<g>, chi:G(K/Q) -> overline{Q_p} s.t. chi(g)=zeta_d^(-1)
    1329             :  * Chi[a]=k, k<0  => chi(a)=0
    1330             :  *           k>=0 => chi(a)=zeta_d^(-k)
    1331             :  * psi=chi^j, j in C : repre. of inj. odd char.
    1332             :  * psi(p)==1 <=> chi(p)^j==0 <=> j*Chi[p]=0 (mod d) <=> Chi[p]==0
    1333             :  *
    1334             :  * K = [H1, Chi, Minpol, C, [d_chi, n_conj]] */
    1335             : static long
    1336        3262 : imag_cyc_lam(GEN K, long p)
    1337             : {
    1338        3262 :   pari_sp av = avma;
    1339        3262 :   GEN Chi = gel(K, 2), MinPol = gel(K, 3), C = gel(K, 4), MinPol2;
    1340        3262 :   long d_K = K_get_d(K), f = K_get_f(K), n_conj = K_get_nconj(K);
    1341        3262 :   long i, q0, n, pmodf = p%f, n_done = 0;
    1342             :   ulong pn1, pM;
    1343        3262 :   GEN p0 = utoi(p), Gam, Lam, xi, K2;
    1344             : 
    1345        3262 :   q0 = (f%p)? f*p: f;
    1346        3262 :   Lam = const_vecsmall(n_conj, -1);
    1347        3262 :   if (pmodf==0 || Chi[pmodf])  /* 1st trial is Bernoulli number */
    1348             :   {
    1349        3052 :     MinPol2 = ZX_to_Flx(MinPol, p*p);  /* p^2 for B_{1,chi} */
    1350        3052 :     K2 = shallowconcat(K, mkvec2(MinPol2, zx_ber_num(Chi, f, d_K)));
    1351       17528 :     for (i=1; i<=n_conj; i++)
    1352       14476 :       if ((Lam[i] = lam_chi_ber(K2, p, C[i])) == 0) n_done++;
    1353        3052 :     if (n_conj==n_done) return gc_long(av, 0);  /* all chi-parts trivial */
    1354             :   }
    1355         413 :   pM = pn1 = p;
    1356         413 :   for (n=1; n>=0; n++)  /* 2nd trial is Stickelberger element */
    1357             :   {
    1358         413 :     pn1 *= p; /* p^(n+1) */
    1359         413 :     if (f == p)
    1360             :     { /* do not use set_minpol: it returns a new pol for each call */
    1361           0 :       GEN fac, cofac, v, pol = polcyclo(d_K, 0);
    1362           0 :       pM = pn1 * p; /* p^(n+2) */
    1363           0 :       fac = FpX_red(MinPol, p0); cofac = FpX_div(pol, fac, p0);
    1364           0 :       v = ZpX_liftfact(pol, mkvec2(fac, cofac), utoipos(pM), p0, n+2);
    1365           0 :       MinPol2 = gel(v, 1);
    1366             :     }
    1367         413 :     Gam = set_gam((1+q0)%pn1, p, n);
    1368         413 :     MinPol2 = ZX_to_Flx(MinPol, pM);
    1369           0 :     xi = (f==p)? get_xi_2(Chi, Gam, p, f, n, d_K, pM)
    1370         413 :                : get_xi_1(Chi, Gam, p, f, n, d_K, pM);
    1371         413 :     K2 = shallowconcat(K, mkvec2(MinPol2, xi));
    1372        2205 :     for (i=1; i<=n_conj; i++)
    1373        1792 :       if (Lam[i]<0 && (Lam[i] = lam_chi_xi(K2, p, C[i], n)) >= 0) n_done++;
    1374         413 :     if (n_conj==n_done) break;
    1375             :   }
    1376         413 :   return gc_long(av, zv_sum(Lam));
    1377             : }
    1378             : static GEN
    1379         329 : GHinit(long f, GEN HH, GEN *pcycGH)
    1380             : {
    1381         329 :   GEN G = znstar0(utoipos(f), 1);
    1382         329 :   GEN U, Ui, cycG, cycGH, ncycGH, gG, gGH, vChar, vH1, P, gH = gel(HH, 1);
    1383         329 :   long i, expG, n_f, lgH = lg(gH); /* gens. of H */
    1384         329 :   P = cgetg(lgH, t_MAT);
    1385         805 :   for (i = 1; i < lgH; i++) gel(P,i) = Zideallog(G, utoi(gH[i]));
    1386             : 
    1387             :   /* group structure of G/H */
    1388         329 :   cycG = znstar_get_cyc(G);
    1389         329 :   expG = itou(gel(cycG, 1));
    1390             :   /* gG generators of G, gGH generators of G/H: gGH = g.Ui, g = gGH.U */
    1391         329 :   cycGH = ZM_snf_group(hnfmodid(P, cycG), &U, &Ui);
    1392         329 :   ncycGH = cyc_normalize(cycGH);
    1393         329 :   gG = ZV_to_Flv(znstar_get_gen(G), f); /* gens. of G */
    1394             :   /* generators of G/H */
    1395         329 :   gGH = Flv_FlvV_factorback(gG, ZM_to_Flm(Ui, expG), f);
    1396         329 :   vChar = chargalois(cycGH, NULL); n_f = lg(vChar)-2;
    1397         329 :   vH1 = cgetg(n_f+1, t_VEC);
    1398       19124 :   for (i = 1; i <= n_f; i++)
    1399             :   { /* skip trivial character */
    1400       18795 :     GEN chi = gel(vChar,i+1), nchi = char_normalize(chi, ncycGH);
    1401       18795 :     GEN chiG, E, H1, C = gel(nchi, 2);
    1402       18795 :     long e, he, gen, d = itou(gel(nchi, 1));
    1403             :     /* chi(prod g[i]^e[i]) = e(sum e[i]*C[i] / d), chi has order d = #(G/H1)*/
    1404       18795 :     E = Flv_extgcd(C, d); /* \sum C[i]*E[i] = 1 in Z/dZ */
    1405             : 
    1406       18795 :     chiG = zncharlift(chi, ncycGH, U, cycG);
    1407       18795 :     H1 =  charker(cycG, chiG); /* H1 < G with G/H1 cyclic */
    1408       18795 :     e = itou( zncharconductor(G, chiG) ); /* cond H1 = cond chi */
    1409       18795 :     H1 = Flv_FlvV_factorback(zv_to_Flv(gG, e), ZM_to_Flm(H1, expG), e);
    1410       18795 :     gen = Flv_factorback(zv_to_Flv(gGH, e), E, e);
    1411       18795 :     H1 = znstar_generate(e, H1); /* G/H1 = <gen>, chi(gen) = e(1/d) */
    1412       18795 :     he = eulerphiu(e) / d;
    1413             :     /* G/H1 = <gen> cyclic of index d, e = cond(H1) */
    1414       18795 :     gel(vH1,i) = mkvec3(H1, mkvecsmall5(d,e,he,srh_1(H1), gen),
    1415             :                         subgp2ary(H1, he));
    1416             :   }
    1417         329 :   if (pcycGH) *pcycGH = cycGH;
    1418         329 :   return vH1;
    1419             : }
    1420             : 
    1421             : /* aH=g^iH ==> chi(a)=zeta_n^(-i); Chi[g^iH]=i; Chi[0] is never accessed */
    1422             : static GEN
    1423       13419 : get_chi(GEN H1)
    1424             : {
    1425       13419 :   GEN H = _get_H(H1);
    1426       13419 :   long i, j, gi, d = _get_d(H1), f = _get_f(H1), h = _get_h(H1), g = _get_g(H1);
    1427       13419 :   GEN Chi = const_vecsmall(f-1, -1);
    1428             : 
    1429     5584159 :   for (j=1; j<=h; j++) Chi[H[j]] = 0; /* i = 0 */
    1430      396592 :   for (i = 1, gi = g; i < d; i++)
    1431             :   {
    1432    40492081 :     for (j=1; j<=h; j++) Chi[Fl_mul(gi, H[j], f)] = i;
    1433      383173 :     gi = Fl_mul(gi, g, f);
    1434             :   }
    1435       13419 :   return Chi;
    1436             : }
    1437             : 
    1438             : static void
    1439          14 : errpdiv(const char *f, GEN p, long d)
    1440             : {
    1441          14 :   pari_err_DOMAIN(f, "p", "divides",
    1442          14 :                   strtoGENstr(stack_sprintf("[F:Q] = %ld", d)), p);
    1443           0 : }
    1444             : /* p odd doesn't divide degF; return lambda invariant if n==0 and
    1445             :  * iwasawa polynomials if n>=1 */
    1446             : static GEN
    1447          49 : abeliwasawa(long p, long f, GEN HH, long degF, long n)
    1448             : {
    1449          49 :   long lam = 0, i, n_f;
    1450          49 :   GEN vH1, vData, z = nullvec(), p0 = utoi(p) ;
    1451             : 
    1452          49 :   vH1 = GHinit(f, HH, NULL); n_f = lg(vH1)-1;
    1453          49 :   vData = const_vec(degF, NULL);
    1454        6608 :   for (i=1; i<=n_f; i++) /* prescan. set Teichmuller */
    1455             :   {
    1456        6573 :     GEN H1 = gel(vH1,i);
    1457        6573 :     long d_K = _get_d(H1), f_K = _get_f(H1), g_K = _get_g(H1);
    1458             : 
    1459        6573 :     if (f_K == d_K+1 && p == f_K)  /* found K=Q(zeta_p) */
    1460             :     {
    1461          14 :       long d_chi = 1, n_conj = eulerphiu(d_K);
    1462          14 :       GEN C = set_C(p, d_K, d_chi, n_conj);
    1463          14 :       long minpow = n? 2*n+1: 2;
    1464          14 :       GEN MinPol = set_minpol_teich(g_K, p0, minpow);
    1465          14 :       gel(vData, d_K) = mkvec4(MinPol, C, NULL, mkvecsmall2(d_chi, n_conj));
    1466          14 :       break;
    1467             :     }
    1468             :   }
    1469             : 
    1470        6664 :   for (i=1; i<=n_f; i++)
    1471             :   {
    1472        6615 :     GEN H1 = gel(vH1,i), z1, Chi, K;
    1473        6615 :     long d_K = _get_d(H1), s = _get_s(H1);
    1474             : 
    1475        6615 :     if (s) continue;  /* F is real */
    1476             : #ifdef DEBUG
    1477             :     err_printf("  handling %s cyclic subfield K, deg(K)=%ld, cond(K)=%ld\n",
    1478             :         s? "a real": "an imaginary", d_K, _get_f(H1));
    1479             : #endif
    1480        3318 :     if (!gel(vData, d_K))
    1481             :     {
    1482         126 :       long d_chi = order_f_x(d_K, p), n_conj = eulerphiu(d_K)/d_chi;
    1483         126 :       GEN C = set_C(p, d_K, d_chi, n_conj);
    1484         126 :       long minpow = n? n+1: 2;
    1485         126 :       GEN MinPol = set_minpol(d_K, p0, minpow, n_conj);
    1486         126 :       gel(vData, d_K) = mkvec4(MinPol, C, NULL, mkvecsmall2(d_chi, n_conj));
    1487             :     }
    1488        3318 :     Chi = get_chi(H1);
    1489        3318 :     K = shallowconcat(mkvec2(H1, Chi), gel(vData, d_K));
    1490        3318 :     if (n==0) lam += imag_cyc_lam(K, p);
    1491          56 :     else if (lg(z1 = imag_cyc_pol(K, p, n)) > 1) z = shallowconcat(z, z1);
    1492             :   }
    1493          49 :   return n? z: mkvecs(lam);
    1494             : }
    1495             : 
    1496             : static GEN
    1497          77 : ary2mat(GEN x, long n)
    1498             : {
    1499             :   long i, j;
    1500          77 :   GEN z = cgetg(n+1,t_MAT);
    1501         182 :   for (i=1; i<=n; i++)
    1502             :   {
    1503         105 :     gel(z,i) = cgetg(n+1,t_COL);
    1504         280 :     for (j=1; j<=n; j++) gmael(z, i, j) = utoi(x[(i-1)*n+j-1]);
    1505             :   }
    1506          77 :   return z;
    1507             : }
    1508             : 
    1509             : static long
    1510           0 : is_cyclic(GEN x)
    1511             : {
    1512           0 :   GEN y = gel(x, 2);
    1513           0 :   long i, l = lg(y), n = 0;
    1514           0 :   for (i = 1; i < l; i++) if (signe(gel(y,i))) n++;
    1515           0 :   return n <= 1;
    1516             : }
    1517             : 
    1518             : static GEN
    1519          77 : make_p_part(GEN y, ulong p, long d_pow)
    1520             : {
    1521          77 :   long i, l = lg(y);
    1522          77 :   GEN z = cgetg(l, t_VECSMALL);
    1523         182 :   for (i = 1; i < l; i++) z[i] = signe(gel(y,i))? Z_lval(gel(y,i), p): d_pow;
    1524          77 :   return z;
    1525             : }
    1526             : 
    1527             : static GEN
    1528          77 : structure_MLL(GEN y, long d_pow)
    1529             : {
    1530          77 :   long y0, i, l = lg(y);
    1531          77 :   GEN x = gen_0, E = cgetg(l, t_VEC);
    1532         182 :   for (i = 1; i < l; i++)
    1533             :   {
    1534         105 :     if ((y0 = d_pow-y[i]) < 0) y0 = 0;
    1535         105 :     x = addiu(x, y0);
    1536         105 :     gel(E, l-i) = utoi(y0);
    1537             :   }
    1538          77 :   return mkvec2(x, E);
    1539             : }
    1540             : 
    1541             : static long
    1542          14 : find_del_el(GEN *oldgr, GEN newgr, long n, long n_el, long d_chi)
    1543             : {
    1544          14 :   if (n_el==1) return 1;
    1545          14 :   if (equalis(gmael(newgr, 2, 1), n_el)) return n;
    1546          14 :   if (cmpii(gel(*oldgr, 1), gel(newgr, 1)) >= 0) return n;
    1547          14 :   if (n > 1 && is_cyclic(newgr)) { *oldgr = newgr; return n-1; }
    1548          14 :   if (n == n_el) return n;
    1549          14 :   if (cmpis(gel(newgr, 1), n*d_chi) < 0) return n;
    1550          14 :   return 0;
    1551             : }
    1552             : 
    1553             : static GEN
    1554           7 : subgr2vecsmall(GEN H, long h, long f)
    1555             : {
    1556             :   long i;
    1557           7 :   GEN z = const_vecsmall(f-1, 0); /* f=lg(z) */
    1558        2023 :   for (i=1; i<=h; i++) z[H[i]] = 1; /* H[i]!=0 */
    1559           7 :   return z;
    1560             : }
    1561             : 
    1562             : /* K is the subfield of Q(zeta_f) with degree d corresponding to the subgroup
    1563             :  * H in (Z/fZ)^*; for a divisor e of f, zeta_e \in K <=> H \subset He. */
    1564             : static long
    1565         119 : root_of_1(long f, GEN H)
    1566             : {
    1567         119 :   GEN g = gel(H, 1); /* generators */
    1568         119 :   long e = f, i, l = lg(g);
    1569         119 :   for (i = 1; i < l; i++)
    1570             :   {
    1571          98 :     e = cgcd(e, g[i] - 1);
    1572          98 :     if (e <= 2) return 2;
    1573             :   }
    1574          21 :   return odd(e)? (e<<1): e;
    1575             : }
    1576             : 
    1577             : static long
    1578         259 : find_ele(GEN H)
    1579             : {
    1580         259 :   long i, f=lg(H);
    1581      369852 :   for (i=1; i<f; i++) if (H[i]) return i;
    1582           7 :   return 0;
    1583             : }
    1584             : 
    1585             : static void
    1586         252 : delete_ele(GEN H, long j, long el)
    1587             : {
    1588         252 :   long f = lg(H), x = 1;
    1589        2016 :   do H[Fl_mul(j,x,f)] = 0;
    1590        2016 :   while ((x=Fl_mul(x,el,f))!=1);
    1591         252 : }
    1592             : 
    1593             : static GEN
    1594           7 : get_coset(GEN H, long h, long f, long el)
    1595             : {
    1596           7 :   long i, j, k = h/order_f_x(f, el);
    1597           7 :   GEN H2, coset = const_vecsmall(k, 0);
    1598           7 :   H2 = subgr2vecsmall(H, h, f);
    1599         259 :   for (i=0; (j=find_ele(H2))>0; i++)
    1600             :   {
    1601         252 :     coset[1+i] = j;
    1602         252 :     delete_ele(H2, j, el);
    1603             :   }
    1604           7 :   if (i != k) pari_err_BUG("failed to find coset\n");
    1605           7 :   return coset;
    1606             : }
    1607             : 
    1608             : static long
    1609        3024 : srh_pol(GEN xpows, GEN vn, GEN pols, long el, long k, long f)
    1610             : {
    1611        3024 :   pari_sp av = avma;
    1612        3024 :   long i, j, l = lg(pols), d = degpol(gel(pols, 1));
    1613        3024 :   GEN pol = gel(pols, 1);
    1614             : 
    1615      654696 :   for (i=1; i<l; i++)
    1616             :   {
    1617             :     GEN x, y, z;
    1618      654696 :     if (vn[i]==0) continue;
    1619      331170 :     y = gel(pols, vn[i]);
    1620      331170 :     z = pol0_Flx(0);
    1621     3311700 :     for (j=0; j<=d; j++)
    1622     2980530 :       z = Flx_add(z, Flx_Fl_mul(gel(xpows, 1+Fl_mul(j, k, f)), y[2+j], el), el);
    1623      331170 :     x = Flx_rem(z, pol, el);
    1624      331170 :     if (lg(x)==2)
    1625        3024 :     {vn[i] = 0; return gc_long(av, i); }  /* pols[i] is min pol of zeta_f^k */
    1626             :   }
    1627           0 :   pari_err_BUG("subcyclopclgp [minimal polinomial]");
    1628           0 :   return 0;  /* to suppress warning */
    1629             : }
    1630             : 
    1631             : /* e_chi[i mod dK] = chi(i*j), i = 0..dK-1; beware: e_chi is translated ! */
    1632             : static GEN
    1633       35857 : get_e_chi(GEN K, ulong j, ulong d, ulong *pdK)
    1634             : {
    1635       35857 :   ulong i, dK = K_get_d(K);
    1636       35857 :   GEN TR = gel(K,4) + 2, e_chi = cgetg(dK+1, t_VECSMALL) + 1;
    1637       35857 :   if (j == 1)
    1638      288267 :     for (i = 0; i < dK; i++) e_chi[i] = umodiu(gel(TR, i), d);
    1639             :   else
    1640      981980 :     for (i = 0; i < dK; i++) e_chi[i] = umodiu(gel(TR, Fl_mul(i, j, dK)), d);
    1641       35857 :   *pdK = dK; return e_chi;
    1642             : }
    1643             : static GEN
    1644        5145 : get_e_chi_ll(GEN K, ulong j, GEN d)
    1645             : {
    1646        5145 :   ulong i, dK = umael3(K, 1, 2, 1);
    1647        5145 :   GEN TR = gel(K,4) + 2, e_chi = cgetg(dK+1, t_VEC) + 1;
    1648      246785 :   for (i = 0; i < dK; i++) gel(e_chi,i) = modii(gel(TR, Fl_mul(i, j, dK)), d);
    1649        5145 :   return e_chi;
    1650             : }
    1651             : 
    1652             : /* el=1 (mod f) */
    1653             : static long
    1654           0 : chk_el_real_f(GEN K, ulong p, ulong d_pow, ulong el, ulong j0)
    1655             : {
    1656           0 :   pari_sp av = avma;
    1657           0 :   GEN H = K_get_H(K);
    1658           0 :   ulong d_K, f = K_get_f(K), h = K_get_h(K), g_K = K_get_g(K);
    1659           0 :   ulong  i, j, gi, d = upowuu(p, d_pow), dp = d*p;
    1660           0 :   ulong g_el, z_f, flag = 0, el1f = (el-1)/f, el1dp = (el-1)/dp;
    1661           0 :   GEN e_chi = get_e_chi(K, j0, dp, &d_K);
    1662           0 :   GEN vz_f, xi_el = cgetg(d_K+1, t_VECSMALL)+1;
    1663             : 
    1664           0 :   g_el = pgener_Fl(el);
    1665           0 :   z_f = Fl_powu(g_el, el1f, el);
    1666           0 :   vz_f = Fl_powers(z_f, f-1, el)+1;
    1667             : 
    1668           0 :   for (gi = 1, i = 0; i < d_K; i++)
    1669             :   {
    1670           0 :     ulong x = 1;
    1671           0 :     for (j = 1; j <= h; j++)
    1672             :     {
    1673           0 :       ulong y = Fl_mul(H[j], gi, f);
    1674           0 :       x = Fl_mul(x, vz_f[y]-1, el); /* vz_f[y] = z_f^y  */
    1675             :     }
    1676           0 :     gi = Fl_mul(gi, g_K, f);
    1677           0 :     xi_el[i] = x;  /* xi_el[i] = xi^{g_K^i} mod el */
    1678             :   }
    1679           0 :   for (i=0; i<d_K; i++)
    1680             :   {
    1681           0 :     ulong x = 1;
    1682           0 :     for (j=0; j<d_K; j++)
    1683           0 :       x = Fl_mul(x, Fl_powu(xi_el[j], e_chi[(i+j)%d_K], el), el);
    1684           0 :     if ((x = Fl_powu(x, el1dp, el))!=1) flag = 1;
    1685           0 :     if (Fl_powu(x, p, el)!=1) return gc_long(av,0);
    1686             :   }
    1687           0 :   return gc_long(av, flag?1:0);
    1688             : }
    1689             : 
    1690             : /* For a cyclic field K contained in Q(zeta_f),
    1691             :  * computes minimal polynomial T of theta=Tr_{Q(zeta_f)/K}(zeta_f) over Q
    1692             :  * and conjugates of theta */
    1693             : static GEN
    1694          14 : minpol_theta(GEN K)
    1695             : {
    1696          14 :   GEN HH = gmael3(K,1,1,1);
    1697          14 :   return galoissubcyclo(utoi(K_get_f(K)), HH, 0, 0);
    1698             : }
    1699             : 
    1700             : /*  xi[1+i] = i-th conj of xi = Tr_{Q(zeta_f)/K}(1-zeta_f).
    1701             :  * |1-(cos(x)+i*sin(x))|^2 = 2(1-cos(x)) */
    1702             : static GEN
    1703           0 : xi_approx(GEN K, long prec)
    1704             : {
    1705           0 :   pari_sp av = avma;
    1706           0 :   GEN H = K_get_H(K);
    1707           0 :   long d_K = K_get_d(K), f = K_get_f(K), h = K_get_h(K), g_K = K_get_g(K);
    1708           0 :   GEN xi = cgetg(d_K+1, t_COL), vz_f = grootsof1(f, prec);
    1709           0 :   long i, j, g = 1, h2 = h>>1;
    1710           0 :   for (i=1; i<=d_K; i++)
    1711             :   {
    1712           0 :     GEN y = real_1(prec);
    1713           0 :     for (j=1; j<=h2; j++)
    1714             :     {
    1715           0 :       GEN z = gmael(vz_f, 1+Fl_mul(H[j], g, f), 1);
    1716           0 :       y = mulrr(y, shiftr(subsr(1, z), 1));
    1717             :     }
    1718           0 :     gel(xi, i) = y;
    1719           0 :     g = Fl_mul(g, g_K, f);
    1720             :   }
    1721           0 :   return gerepilecopy(av, xi);
    1722             : }
    1723             : 
    1724             : static GEN
    1725          47 : theta_xi_el(GEN K, ulong el)
    1726             : {
    1727          47 :   pari_sp av = avma;
    1728          47 :   GEN H = K_get_H(K);
    1729          47 :   ulong d_K = K_get_d(K), f = K_get_f(K), h = K_get_h(K), g_K = K_get_g(K);
    1730          47 :   GEN theta = cgetg(d_K+1, t_VECSMALL), xi = cgetg(d_K+1, t_VECSMALL), vz_f;
    1731          47 :   ulong i, j, g = 1, x, y, g_el, z_f;
    1732             : 
    1733          47 :   g_el = pgener_Fl(el);
    1734          47 :   z_f = Fl_powu(g_el, (el-1)/f, el);
    1735          47 :   vz_f = Fl_powers(z_f, f-1, el);
    1736        1109 :   for (i=1; i<=d_K; i++)
    1737             :   {
    1738        1062 :     x = 0; y = 1;
    1739       28254 :     for (j=1; j<=h; j++)
    1740             :     {
    1741       27192 :       ulong z = vz_f[1+Fl_mul(H[j], g, f)];
    1742       27192 :       x = Fl_add(x, z, el);
    1743       27192 :       y = Fl_mul(y, z-1, el);
    1744             :     }
    1745        1062 :     theta[i] = x;
    1746        1062 :     xi[i] = y;
    1747        1062 :     g = Fl_mul(g, g_K, f);
    1748             :   }
    1749          47 :   return gerepilecopy(av, mkvec2(theta, xi));
    1750             : }
    1751             : 
    1752             : static GEN
    1753          47 : make_Xi(GEN xi, long d)
    1754             : {
    1755             :   long i, j;
    1756          47 :   GEN Xi = cgetg(d+1, t_MAT);
    1757        1109 :   for (j=0; j<d; j++)
    1758             :   {
    1759        1062 :     GEN x = cgetg(d+1, t_VECSMALL);
    1760        1062 :     gel(Xi, 1+j) = x;
    1761       27714 :     for (i=0; i<d; i++) x[1+i] = xi[1+(i+j)%d];
    1762             :   }
    1763          47 :   return Xi;
    1764             : }
    1765             : 
    1766             : static GEN
    1767          47 : make_Theta(GEN theta, ulong d, ulong el)
    1768             : {
    1769             :   ulong i;
    1770          47 :   GEN Theta = cgetg(d+1, t_MAT);
    1771        1109 :   for (i=1; i<=d; i++) gel(Theta, i) = Fl_powers(theta[i], d-1, el);
    1772          47 :   return Flm_inv(Theta, el);
    1773             : }
    1774             : 
    1775             : static GEN
    1776          47 : Xi_el(GEN K, GEN tInvA, ulong el)
    1777             : {
    1778          47 :   pari_sp av = avma;
    1779          47 :   ulong d_K = K_get_d(K);
    1780          47 :   GEN tx = theta_xi_el(K, el), Theta, Xi, X;
    1781             : 
    1782          47 :   if ((Theta = make_Theta(gel(tx, 1), d_K, el))==NULL) return NULL;
    1783          47 :   Xi = make_Xi(gel(tx, 2), d_K);
    1784          47 :   X = Flm_mul(Flm_mul(Xi, Theta, el), ZM_to_Flm(tInvA, el), el);
    1785          47 :   return gerepilecopy(av, X);
    1786             : }
    1787             : 
    1788             : static GEN
    1789           0 : pol_xi_el(GEN K, ulong el)
    1790             : {
    1791           0 :   pari_sp av = avma;
    1792           0 :   ulong d_K = K_get_d(K), f = K_get_f(K), h = K_get_h(K), g_K = K_get_g(K);
    1793           0 :   GEN H = K_get_H(K), xi = cgetg(d_K+1, t_VECSMALL), vz_f;
    1794           0 :   ulong i, j, g = 1, y, g_el, z_f;
    1795             : 
    1796           0 :   g_el = pgener_Fl(el);
    1797           0 :   z_f = Fl_powu(g_el, (el-1)/f, el);
    1798           0 :   vz_f = Fl_powers(z_f, f-1, el);
    1799           0 :   for (i=1; i<=d_K; i++)
    1800             :   {
    1801           0 :     y = 1;
    1802           0 :     for (j=1; j<=h; j++)
    1803             :     {
    1804           0 :       ulong z = vz_f[1+Fl_mul(H[j], g, f)];
    1805           0 :       y = Fl_mul(y, z-1, el);
    1806             :     }
    1807           0 :     xi[i] = y;
    1808           0 :     g = Fl_mul(g, g_K, f);
    1809             :   }
    1810           0 :   return gerepilecopy(av, Flv_roots_to_pol(xi, el, 0));
    1811             : }
    1812             : 
    1813             : /* theta[1+i] = i-th conj of theta; xi[1+i] = i-th conj of xi. */
    1814             : static GEN
    1815          14 : theta_xi_approx(GEN K, long prec)
    1816             : {
    1817          14 :   pari_sp av = avma;
    1818          14 :   GEN H = K_get_H(K);
    1819          14 :   long d_K = K_get_d(K), f = K_get_f(K), h = K_get_h(K), g_K = K_get_g(K);
    1820          14 :   GEN theta = cgetg(d_K+1, t_VEC), xi = cgetg(d_K+1, t_VEC);
    1821          14 :   GEN vz_f = grootsof1(f, prec);
    1822          14 :   long i, j, g = 1, h2 = h>>1;
    1823             : 
    1824         238 :   for (i=1; i<=d_K; i++)
    1825             :   {
    1826         224 :     GEN x = real_0(prec), y = real_1(prec);
    1827        5068 :     for (j=1; j<=h2; j++)
    1828             :     {
    1829        4844 :       GEN z = gmael(vz_f, 1+Fl_mul(H[j], g, f), 1);
    1830        4844 :       x = addrr(x, z);
    1831        4844 :       y = mulrr(y, shiftr(subsr(1, z), 1));
    1832             :     }
    1833         224 :     gel(theta, i) = shiftr(x, 1);
    1834         224 :     gel(xi, i) = y;
    1835         224 :     g = Fl_mul(g, g_K, f);
    1836             :   }
    1837          14 :   return gerepilecopy(av, mkvec2(theta, xi));
    1838             : }
    1839             : 
    1840             : static GEN
    1841          14 : bound_coeff_xi(GEN K, GEN tInvA)
    1842             : {
    1843          14 :   pari_sp av = avma;
    1844          14 :   long d_K = K_get_d(K), prec = MEDDEFAULTPREC, i;
    1845          14 :   GEN tInvV, R = cgetg(d_K+1, t_MAT), theta_xi = theta_xi_approx(K, prec);
    1846          14 :   GEN theta = gel(theta_xi, 1), xi = gel(theta_xi, 2), x1, x2, bound;
    1847             : 
    1848         238 :   for (i=1; i<=d_K; i++)
    1849             :   {
    1850         224 :     GEN z = gpowers(gel(theta, i), d_K-1);
    1851         224 :     settyp(z, t_COL);
    1852         224 :     gel(R, i) = z;
    1853             :   }
    1854          14 :   tInvV = RgM_mul(RgM_inv(R), tInvA);
    1855          14 :   x1 = gsupnorm(tInvV, prec); x2 = gsupnorm(xi, prec);
    1856          14 :   bound = mulrs(mulrr(x1, x2), 3*d_K);
    1857          14 :   return gerepilecopy(av, bound);
    1858             : }
    1859             : 
    1860             : static GEN
    1861          14 : get_Xi(GEN K, GEN tInvA)
    1862             : {
    1863          14 :   pari_sp av = avma;
    1864             :   GEN M0, XI, EL, Xi;
    1865          14 :   ulong f = K_get_f(K), el, e, n, i;
    1866             :   forprime_t T0;
    1867             : 
    1868          14 :   M0 = bound_coeff_xi(K, tInvA);
    1869          14 :   e = expo(M0)+1; n = e/(BITS_IN_LONG-1); n++;
    1870          14 :   EL = cgetg(1+n, t_VECSMALL);
    1871          14 :   XI = cgetg(1+n, t_VEC);
    1872          14 :   u_forprime_arith_init(&T0, LONG_MAX, ULONG_MAX, 1, f);
    1873             : 
    1874          61 :   for (i=1; i<=n; i++)
    1875             :   {
    1876          47 :     el = u_forprime_next(&T0);
    1877          47 :     while ((Xi=Xi_el(K, tInvA, el))==NULL) el = u_forprime_next(&T0);
    1878          47 :     gel(XI, i) = Xi;
    1879          47 :     EL[i] = el;
    1880             :   }
    1881          14 :   return gerepileupto(av, nmV_chinese_center(XI, EL, NULL));
    1882             : }
    1883             : 
    1884             : /* K is a cyclic field of conductor f with degree d=d_K
    1885             :  * xi = Norm_{Q(zeta_f)/K}(1-zeta_f)
    1886             :  * 1: T, min poly of a=Tr_{Q(zeta_f)/K}(zeta_f) over Q
    1887             :  * 2: B, power basis of K with respect to a
    1888             :  * 3: A, rational matrix s.t. t(v_1,...v_d) = A*t(1,a,...,a^{d-1})
    1889             :  * 4: Xi, integer matrix s.t. t(xi^(1),...,xi^(d)) = Xi*t(v_1,...,v_d) */
    1890             : static GEN
    1891          14 : xi_data_basis(GEN K)
    1892             : {
    1893          14 :   pari_sp av = avma;
    1894          14 :   GEN T = minpol_theta(K);
    1895             :   GEN InvA, A, M, Xi, A_den;
    1896          14 :   GEN D, B = nfbasis(T, &D);
    1897             :   pari_timer ti;
    1898          14 :   if (DEBUGLEVEL>1) timer_start(&ti);
    1899          14 :   A = RgXV_to_RgM(B, lg(B)-1);
    1900          14 :   M = gmael(A, 1, 1);
    1901          14 :   if (!equali1(M)) A = RgM_Rg_div(A, M);
    1902          14 :   InvA = QM_inv(A);
    1903          14 :   A = Q_remove_denom(A, &A_den);
    1904          14 :   if (A_den==NULL) A_den = gen_1;
    1905          14 :   Xi = get_Xi(K, shallowtrans(InvA));
    1906          14 :   if (DEBUGLEVEL>1) timer_printf(&ti, "xi_data_basis");
    1907          14 :   return gerepilecopy(av, mkvec5(T, B, shallowtrans(A), Xi, A_den));
    1908             : }
    1909             : 
    1910             : /* When factorization of polcyclo mod el is difficult, one can try to
    1911             :  * check the condition of el using an integral basis of K.
    1912             :  * This is useful when d_K is small. */
    1913             : static long
    1914          14 : chk_el_real_basis(GEN K, long p, long d_pow, long el, long j0)
    1915             : {
    1916          14 :   pari_sp av = avma;
    1917          14 :   GEN xi = gel(K, 7), T = gel(xi, 1), A = gel(xi, 3), Xi = gel(xi, 4);
    1918          14 :   GEN A_den = gel(xi, 5);
    1919          14 :   ulong i, j, x, found = 0;
    1920             :   GEN v_el, xi_el;
    1921             :   GEN e_chi, xi_e_chi;
    1922             :   ulong d_K, d, dp, el1dp;
    1923             : 
    1924          14 :   if (dvdiu(A_den, el)) return 0;
    1925             : 
    1926          14 :   d = upowuu(p, d_pow); dp = d*p; el1dp = (el-1)/dp;
    1927          14 :   e_chi = get_e_chi(K, j0, dp, &d_K);
    1928          14 :   xi_e_chi = cgetg(d_K+1, t_VECSMALL)+1;
    1929             : 
    1930          14 :   if (DEBUGLEVEL>1) err_printf("chk_el_real_basis: d_K=%ld el=%ld\n",d_K,el);
    1931          14 :   A = ZM_to_Flm(A, el);
    1932          14 :   A = Flm_Fl_mul(A, Fl_inv(umodiu(A_den, el), el), el);
    1933          14 :   x = Flx_oneroot_split(ZX_to_Flx(T, el), el);
    1934          14 :   v_el = Flm_Flc_mul(A, Fl_powers(x, d_K-1, el), el);
    1935          14 :   xi_el = Flm_Flc_mul(ZM_to_Flm(Xi, el), v_el, el);
    1936          14 :   if (DEBUGLEVEL>2) err_printf("el=%ld xi_el=%Ps\n", el, xi_el);
    1937         238 :   for (i=0; i<d_K; i++)
    1938             :   {
    1939         224 :     long z = 1;
    1940        5208 :     for (j=0; j<d_K; j++)
    1941        4984 :       z = Fl_mul(z, Fl_powu(xi_el[1+j], e_chi[(i+j)%d_K], el), el);
    1942         224 :     xi_e_chi[i] = z;
    1943             :   }
    1944          14 :   if (DEBUGLEVEL>2) err_printf("xi_e_chi=%Ps\n", xi_e_chi-1);
    1945         238 :   for (i=0; i<d_K; i++)
    1946             :   {
    1947         224 :     long x = Fl_powu(xi_e_chi[i], el1dp, el);
    1948         224 :     if (x!=1) found = 1;
    1949         224 :     if (Fl_powu(x, p, el)!=1) return gc_long(av, 0);
    1950             :   }
    1951          14 :   return gc_long(av, found);
    1952             : }
    1953             : 
    1954             : static GEN
    1955           0 : bound_pol_xi(GEN K)
    1956             : {
    1957           0 :   pari_sp av = avma;
    1958           0 :   GEN xi = xi_approx(K, MEDDEFAULTPREC);
    1959           0 :   GEN M = real_1(MEDDEFAULTPREC), one = rtor(dbltor(1.0001), MEDDEFAULTPREC);
    1960           0 :   long i, n = lg(xi);
    1961             : 
    1962           0 :   for (i=1; i<n; i++) M = mulrr(M, addrr(one, gel(xi, i)));
    1963           0 :   M = mulrs(M, 3);
    1964           0 :   return gerepilecopy(av, M);
    1965             : }
    1966             : 
    1967             : static GEN
    1968           0 : minpol_xi(GEN K)
    1969             : {
    1970           0 :   pari_sp av = avma;
    1971             :   GEN M0, POL, EL;
    1972           0 :   ulong f = K_get_f(K), el, e, n, i;
    1973             :   forprime_t T0;
    1974             : 
    1975           0 :   M0 = bound_pol_xi(K);
    1976           0 :   e = expo(M0)+1; n = e/(BITS_IN_LONG-1); n++;
    1977           0 :   EL = cgetg(1+n, t_VECSMALL);
    1978           0 :   POL = cgetg(1+n, t_VEC);
    1979           0 :   u_forprime_arith_init(&T0, LONG_MAX, ULONG_MAX, 1, f);
    1980           0 :   for (i=1; i<=n; i++)
    1981             :   {
    1982           0 :     el = u_forprime_next(&T0);
    1983           0 :     gel(POL, i) = pol_xi_el(K, el);
    1984           0 :     EL[i] = el;
    1985             :   }
    1986           0 :   return gerepileupto(av, nxV_chinese_center(POL, EL, NULL));
    1987             : }
    1988             : 
    1989             : static long
    1990           0 : find_conj_el(GEN K, GEN pol, GEN Den)
    1991             : {
    1992           0 :   pari_sp av = avma;
    1993           0 :   GEN H = K_get_H(K);
    1994           0 :   ulong d_K = K_get_d(K), f = K_get_f(K), h = K_get_h(K), g = K_get_g(K);
    1995           0 :   ulong j, k, el, g_el, z_f, xi = 1, xi_g = 1;
    1996           0 :   GEN T = NULL, vz_f;
    1997             : 
    1998           0 :   for (el=f+1; el; el+=f)
    1999           0 :     if (uisprime(el) && dvdiu(Den, el)==0)
    2000             :     {
    2001           0 :       T = ZX_to_Flx(pol, el);
    2002           0 :       T = Flx_Fl_mul(T, Fl_inv(umodiu(Den, el), el), el);
    2003           0 :       break;
    2004             :     }
    2005           0 :   g_el = pgener_Fl(el);
    2006           0 :   z_f = Fl_powu(g_el, (el-1)/f, el);
    2007           0 :   vz_f = Fl_powers(z_f, f-1, el);
    2008           0 :   for (j=1; j<=h; j++)
    2009           0 :     xi = Fl_mul(xi, vz_f[1+H[j]]-1, el);
    2010           0 :   for (j=1; j<=h; j++)
    2011           0 :     xi_g = Fl_mul(xi_g, vz_f[1+Fl_mul(H[j], g, f)]-1, el);
    2012           0 :   for (k=1; k<=d_K; k++)
    2013             :   {
    2014           0 :     xi = Flx_eval(T, xi, el);
    2015           0 :     if (xi == xi_g) break;
    2016             :   }
    2017           0 :   if (xi != xi_g) pari_err_BUG("find_conj_el");
    2018           0 :   return gc_long(av, k);
    2019             : }
    2020             : 
    2021             : /* G = H_1*H_2*...*H_m is cyclic of order n, H_i=<perms[i]>
    2022             :  * G is not necessarily a direct product.
    2023             :  * If p^e || n, then p^e || |H_i| for some i.
    2024             :  * return a generator of G. */
    2025             : static GEN
    2026           0 : find_gen(GEN perms, long n)
    2027             : {
    2028           0 :   GEN fa = factoru(n), P = gel(fa, 1), E = gel(fa, 2);
    2029           0 :   long i, j, l = lg(perms), r = lg(P);
    2030           0 :   GEN gen = cgetg(r, t_VEC), orders = cgetg(l, t_VECSMALL), perm;
    2031             : 
    2032           0 :   for (i=1; i<l; i++) orders[i] = perm_orderu(gel(perms, i));
    2033           0 :   for (i=1; i<r; i++)
    2034             :   {
    2035           0 :     long pe = upowuu(P[i], E[i]);
    2036           0 :     for (j=1; j<l; j++) if (orders[j]%pe==0) break;
    2037           0 :     gel(gen, i) = perm_powu(gel(perms, j), orders[j]/pe);
    2038             :   }
    2039           0 :   perm = gel(gen, 1);
    2040           0 :   for (i=2; i<l; i++) perm = perm_mul(perm, gel(gen, i));
    2041           0 :   return perm;
    2042             : }
    2043             : 
    2044             : /* R is the roots of T. R[1+i] = R[1]^(g_K^i), 0 <= i <= d_K-1
    2045             :  * 1: min poly T of xi over Q
    2046             :  * 2: F(x)\in Q[x] s.t. xi^(g_K)=F(xi) */
    2047             : static GEN
    2048           0 : xi_data_galois(GEN K)
    2049             : {
    2050           0 :   pari_sp av = avma;
    2051             :   GEN T, G, perms, perm, pol, pol2, Den;
    2052           0 :   ulong k, d_K = K_get_d(K);
    2053             :   pari_timer ti;
    2054             : 
    2055           0 :   if (DEBUGLEVEL>1) timer_start(&ti);
    2056           0 :   T = minpol_xi(K);
    2057           0 :   if (DEBUGLEVEL>1) timer_printf(&ti, "minpol_xi");
    2058           0 :   G = galoisinit(T, NULL);
    2059           0 :   if (DEBUGLEVEL>1) timer_printf(&ti, "galoisinit");
    2060           0 :   perms = gal_get_gen(G);
    2061           0 :   perm = (lg(perms)==2)?gel(perms, 1):find_gen(perms, d_K);
    2062           0 :   if (DEBUGLEVEL>1) timer_start(&ti);
    2063           0 :   pol = galoispermtopol(G, perm);
    2064           0 :   if (DEBUGLEVEL>1) timer_printf(&ti, "galoispermtopol");
    2065           0 :   pol = Q_remove_denom(pol, &Den);
    2066           0 :   if (Den==NULL) Den = gen_1;
    2067           0 :   k = find_conj_el(K, pol, Den);
    2068           0 :   if (DEBUGLEVEL>1) timer_printf(&ti,"find_conj");
    2069           0 :   pol2 = galoispermtopol(G, perm_powu(perm, k));
    2070           0 :   pol2 = Q_remove_denom(pol2, &Den);
    2071           0 :   if (Den==NULL) Den = gen_1;
    2072           0 :   return gerepilecopy(av, mkvec3(T, pol2, Den));
    2073             : }
    2074             : 
    2075             : /* If g(X)\in Q[X] s.t. g(xi)=xi^{g_K} was found,
    2076             :  * then we fix an integer x_0 s.t. xi=x_0 (mod el) and construct x_i
    2077             :  * s.t. xi^{g_K^i}=x_i (mod el) using g(X). */
    2078             : static long
    2079           0 : chk_el_real_galois(GEN K, long p, long d_pow, long el, long j0)
    2080             : {
    2081           0 :   pari_sp av = avma;
    2082           0 :   GEN xi = gel(K, 7), T = gel(xi, 1), F = gel(xi, 2), Den = gel(xi, 3);
    2083             :   GEN Fel, xi_el, xi_e_chi, e_chi;
    2084           0 :   ulong i, j, x, found = 0;
    2085             :   ulong d_K, d, dp, el1dp;
    2086             : 
    2087           0 :   if (dvdiu(Den, el)) return 0;
    2088             : 
    2089           0 :   d = upowuu(p, d_pow); dp = d*p; el1dp = (el-1)/dp;
    2090           0 :   e_chi = get_e_chi(K, j0, dp, &d_K);
    2091           0 :   xi_el = cgetg(d_K+1, t_VECSMALL)+1;
    2092           0 :   xi_e_chi = cgetg(d_K+1, t_VECSMALL)+1;
    2093             : 
    2094           0 :   if (DEBUGLEVEL>1) err_printf("chk_el_real_galois: d_K=%ld el=%ld\n",d_K,el);
    2095           0 :   Fel = ZX_to_Flx(F, el);
    2096           0 :   Fel = Flx_Fl_mul(Fel, Fl_inv(umodiu(Den, el), el), el);
    2097           0 :   x = Flx_oneroot_split(ZX_to_Flx(T, el), el);
    2098           0 :   for (i=0; i<d_K; i++) { xi_el[i] = x; x = Flx_eval(Fel, x, el); }
    2099           0 :   if (DEBUGLEVEL>2) err_printf("el=%ld xi_el=%Ps\n", el, xi_el-1);
    2100           0 :   for (i=0; i<d_K; i++)
    2101             :   {
    2102           0 :     long z = 1;
    2103           0 :     for (j=0; j<d_K; j++)
    2104           0 :       z = Fl_mul(z, Fl_powu(xi_el[j], e_chi[(i+j)%d_K], el), el);
    2105           0 :     xi_e_chi[i] = z;
    2106             :   }
    2107           0 :   if (DEBUGLEVEL>2) err_printf("xi_e_chi=%Ps\n", xi_e_chi-1);
    2108           0 :   for (i=0; i<d_K; i++)
    2109             :   {
    2110           0 :     long x = Fl_powu(xi_e_chi[i], el1dp, el);
    2111           0 :     if (x!=1) found = 1;
    2112           0 :     if (Fl_powu(x, p, el)!=1) return gc_long(av, 0);
    2113             :   }
    2114           0 :   return gc_long(av, found);
    2115             : }
    2116             : 
    2117             : /* checks the condition of el using the irreducible polynomial G_K(X) of zeta_f
    2118             :  * over K. G_K(X) mod el is enough for our purpose and it is obtained by
    2119             :  * factoring polcyclo(f) mod el */
    2120             : static long
    2121           7 : chk_el_real_factor(GEN K, long p, long d_pow, long el, long j0)
    2122             : {
    2123           7 :   pari_sp av = avma;
    2124           7 :   GEN H = K_get_H(K);
    2125           7 :   ulong f = K_get_f(K), h = K_get_h(K), g_K = K_get_g(K);
    2126           7 :   ulong i, j, k, d_K, d = upowuu(p, d_pow), dp = d*p, found = 0;
    2127             :   GEN pols, coset, vn_g, polnum, xpows, G_K;
    2128           7 :   ulong el1dp = (el-1)/dp, n_coset, n_g, gi;
    2129           7 :   GEN e_chi = get_e_chi(K, j0, dp, &d_K);
    2130             :   pari_timer ti;
    2131             : 
    2132           7 :   if (DEBUGLEVEL>1) err_printf("chk_el_real_factor: f=%ld el=%ld\n",f,el);
    2133           7 :   coset = get_coset(H, h, f, el);
    2134           7 :   if (DEBUGLEVEL>1)
    2135             :   {
    2136           0 :     timer_start(&ti);
    2137           0 :     err_printf("factoring polyclo(d) (mod %ld)\n",f, el);
    2138             :   }
    2139           7 :   pols = Flx_factcyclo(f, el, 0, 0);
    2140           7 :   if (DEBUGLEVEL>1) timer_printf(&ti,"Flx_factcyclo(%lu,%lu)",f,el);
    2141           7 :   n_coset = lg(coset)-1;
    2142           7 :   n_g = lg(pols)-1;
    2143           7 :   vn_g = identity_perm(n_g);
    2144             : 
    2145           7 :   polnum = const_vec(d_K, NULL);
    2146          91 :   for (i=1; i<=d_K; i++) gel(polnum, i) = const_vecsmall(n_coset, 0);
    2147           7 :   xpows = Flxq_powers(polx_Flx(0), f-1, gel(pols, 1), el);
    2148          91 :   for (gi=1,i=1; i<=d_K; i++)
    2149             :   {
    2150        3108 :     for (j=1; j<=n_coset; j++)
    2151             :     {
    2152        3024 :       long x, conj = Fl_mul(gi, coset[j], f);
    2153        3024 :       x = srh_pol(xpows, vn_g, pols, el, conj, f);
    2154        3024 :       gel(polnum, i)[j] = x;
    2155             :     }
    2156          84 :     gi = Fl_mul(gi, g_K, f);
    2157             :   }
    2158           7 :   G_K = const_vec(d_K, NULL);
    2159          91 :   for (i=1; i<=d_K; i++)
    2160             :   {
    2161          84 :     GEN z = pol1_Flx(0);
    2162        3108 :     for (j=1; j<=n_coset; j++) z = Flx_mul(z, gel(pols, gel(polnum,i)[j]), el);
    2163          84 :     gel(G_K, i) = z;
    2164             :   }
    2165           7 :   if (DEBUGLEVEL>2) err_printf("G_K(x)=%Ps\n",Flx_to_ZX(gel(G_K, 1)));
    2166          91 :   for (k=0; k<d_K; k++)
    2167             :   {
    2168          84 :     long x = 1;
    2169        1092 :     for (i = 0; i < d_K; i++)
    2170             :     {
    2171             :       long x0, t;
    2172        1008 :       x0 = Flx_eval(gel(G_K, 1+i), 1, el);
    2173        1008 :       t = Fl_powu(x0, e_chi[(i+k)%d_K], el);
    2174        1008 :       x = Fl_mul(x, t, el);
    2175             :     }
    2176          84 :     x = Fl_powu(x, el1dp, el);
    2177          84 :     if (x!=1) found = 1;
    2178          84 :     if (Fl_powu(x, p, el)!=1) return gc_long(av, 0);
    2179             :   }
    2180           7 :   return gc_long(av, found);
    2181             : }
    2182             : 
    2183             : static long
    2184          21 : chk_el_real_chi(GEN K, ulong p, ulong d_pow, ulong el, ulong j0, long flag)
    2185             : {
    2186          21 :   ulong f = K_get_f(K);
    2187             : 
    2188          21 :   if (el%f == 1)
    2189           0 :     return chk_el_real_f(K, p, d_pow, el, j0); /* must be faster */
    2190          21 :   if (flag&USE_BASIS)
    2191          14 :     return chk_el_real_basis(K, p, d_pow, el, j0);
    2192           7 :   if (flag&USE_GALOIS_POL)
    2193           0 :     return chk_el_real_galois(K, p, d_pow, el, j0);
    2194           7 :   return chk_el_real_factor(K, p, d_pow, el, j0);
    2195             : }
    2196             : 
    2197             : static long
    2198         616 : chk_ell_real(GEN K, long d2, GEN ell, long j0)
    2199             : {
    2200         616 :   pari_sp av = avma;
    2201         616 :   GEN H = K_get_H(K);
    2202         616 :   ulong f = K_get_f(K), h = K_get_h(K), g_K = K_get_g(K);
    2203             :   ulong d_K, i, j, gi;
    2204         616 :   GEN e_chi = get_e_chi(K, j0, d2, &d_K);
    2205         616 :   GEN g_ell, z_f, vz_f, xi_el = cgetg(d_K+1, t_VEC)+1;
    2206         616 :   GEN ell_1 = subiu(ell,1), ell1d2 = diviuexact(ell_1, d2);
    2207             : 
    2208         616 :   g_ell = pgener_Fp(ell);
    2209         616 :   z_f = Fp_pow(g_ell, diviuexact(subiu(ell, 1), f), ell);
    2210         616 :   vz_f = Fp_powers(z_f, f-1, ell)+1;
    2211       12964 :   for (gi=1, i=0; i<d_K; i++)
    2212             :   {
    2213       12348 :     GEN x = gen_1;
    2214     1074948 :     for (j = 1; j <= h; j++)
    2215             :     {
    2216     1062600 :       ulong y = Fl_mul(H[j], gi, f);
    2217     1062600 :       x = Fp_mul(x, subiu(gel(vz_f, y), 1), ell); /* vz_f[y] = z_f^y  */
    2218             :     }
    2219       12348 :     gi = Fl_mul(gi, g_K, f);
    2220       12348 :     gel(xi_el, i) = x;  /* xi_el[i]=xi^{g_K^i} */
    2221             :   }
    2222        1421 :   for (i=0; i<d_K; i++)
    2223             :   {
    2224        1386 :     GEN x = gen_1;
    2225       31976 :     for (j=0; j<d_K; j++)
    2226       30590 :       x = Fp_mul(x, Fp_powu(gel(xi_el, j), e_chi[(i+j)%d_K], ell), ell);
    2227        1386 :     x = Fp_pow(x, ell1d2, ell);
    2228        1386 :     if (!equali1(x)) return gc_long(av, 0);
    2229             :   }
    2230          35 :   return gc_long(av, 1);
    2231             : }
    2232             : 
    2233             : static GEN
    2234          21 : next_el_real(GEN K, long p, long d_pow, GEN elg, long j0, long flag)
    2235             : {
    2236          21 :   GEN Chi = gel(K, 2);
    2237          21 :   ulong f = K_get_f(K), h = K_get_h(K), d = upowuu(p, d_pow), d2 = d*d;
    2238          21 :   ulong D = (flag & USE_F)? d2*f: d2<<1, el = elg[1] + D;
    2239             : 
    2240             :   /* O(el*h) -> O(el*log(el)) by FFT */
    2241          21 :   if (1000 < h && el < h) { el = (h/D)*D+1; if (el < h) el += D; }
    2242          21 :   if (flag&USE_F)  /* el=1 (mod f) */
    2243             :   {
    2244           0 :     for (;; el += D)
    2245           0 :       if (uisprime(el) && chk_el_real_f(K, p, d_pow, el, j0)) break;
    2246             :   }
    2247             :   else
    2248             :   {
    2249         413 :     for (;; el += D)
    2250         455 :       if (Chi[el%f]==0 && uisprime(el) &&
    2251          42 :           chk_el_real_chi(K, p, d_pow, el, j0, flag)) break;
    2252             :   }
    2253          21 :   return mkvecsmall2(el, pgener_Fl(el));
    2254             : }
    2255             : 
    2256             : static GEN
    2257          35 : next_ell_real(GEN K, GEN ellg, long d2, GEN df0l, long j0)
    2258             : {
    2259          35 :   GEN ell = gel(ellg, 1);
    2260        7602 :   for (ell = addii(ell, df0l);; ell = addii(ell, df0l))
    2261        7602 :     if (BPSW_psp(ell) && chk_ell_real(K, d2, ell, j0))
    2262          35 :       return mkvec2(ell, pgener_Fp(ell));
    2263             : }
    2264             : 
    2265             : /* #velg >= n */
    2266             : static long
    2267           0 : delete_el(GEN velg, long n)
    2268             : {
    2269             :   long i, l;
    2270           0 :   if (DEBUGLEVEL>1) err_printf("deleting %ld ...\n", gmael(velg, n, 1));
    2271           0 :   for (l = lg(velg)-1; l >= 1; l--) if (gel(velg, l)) break;
    2272           0 :   for (i = n; i < l; i++) gel(velg, i) = gel(velg, i+1);
    2273           0 :   return l;
    2274             : }
    2275             : 
    2276             : /* velg has n components */
    2277             : static GEN
    2278          21 : set_ell_real(GEN K, GEN velg, long n, long d_chi, long d2, long f0, long j0)
    2279             : {
    2280          21 :   long i, n_ell = n*d_chi;
    2281          21 :   GEN z = cgetg(n_ell + 1, t_VEC);
    2282          21 :   GEN df0l = muluu(d2, f0), ellg = mkvec2(gen_1, gen_1);
    2283          42 :   for (i=1; i<=n; i++) df0l = muliu(df0l, gel(velg, i)[1]);
    2284          56 :   for (i=1; i<=n_ell; i++) ellg = gel(z, i)= next_ell_real(K, ellg, d2, df0l, j0);
    2285          21 :   return z;
    2286             : }
    2287             : 
    2288             : static GEN
    2289         182 : G_K_vell(GEN K, GEN vellg, ulong gk)
    2290             : {
    2291         182 :   pari_sp av = avma;
    2292         182 :   GEN H = K_get_H(K);
    2293         182 :   ulong f = K_get_f(K), h = K_get_h(K);
    2294         182 :   GEN z_f, vz_f, A, P, M, z =  cgetg(h+1, t_VEC);
    2295         182 :   ulong i, lv = lg(vellg);
    2296             : 
    2297         182 :   A=cgetg(lv, t_VEC);
    2298         182 :   P=cgetg(lv, t_VEC);
    2299         728 :   for (i=1; i<lv; i++)
    2300             :   {
    2301         546 :     GEN ell = gmael(vellg, i, 1), g_ell = gmael(vellg, i, 2);
    2302         546 :     gel(A, i) = Fp_pow(g_ell, diviuexact(subiu(ell, 1), f), ell);
    2303         546 :     gel(P, i) = ell;
    2304             :   }
    2305         182 :   z_f = ZV_chinese(A, P, &M);
    2306         182 :   vz_f = Fp_powers(z_f, f-1, M)+1;
    2307        3822 :   for (i=1; i<=h; i++) gel(z, i) = gel(vz_f, Fl_mul(H[i], gk, f));
    2308         182 :   return gerepilecopy(av, FpV_roots_to_pol(z, M, 0));
    2309             : }
    2310             : 
    2311             : /* f=cond(K), M=product of ell in vell, G(K/Q)=<g_K>
    2312             :  * G_K[1+i]=minimal polynomial of zeta_f^{g_k^i} over K mod M, 0 <= i < d_K */
    2313             : static GEN
    2314           7 : make_G_K(GEN K, GEN vellg)
    2315             : {
    2316           7 :   ulong d_K = K_get_d(K), f = K_get_f(K), g_K = K_get_g(K);
    2317           7 :   GEN G_K = cgetg(d_K+1, t_VEC);
    2318           7 :   ulong i, g = 1;
    2319             : 
    2320         189 :   for (i=0; i<d_K; i++)
    2321             :   {
    2322         182 :     gel(G_K, 1+i) = G_K_vell(K, vellg, g);
    2323         182 :     g = Fl_mul(g, g_K, f);
    2324             :   }
    2325           7 :   return G_K;
    2326             : }
    2327             : 
    2328             : static GEN
    2329          12 : G_K_p(GEN K, GEN ellg, ulong gk)
    2330             : {
    2331          12 :   pari_sp av = avma;
    2332          12 :   ulong i, f = K_get_f(K), h = K_get_h(K);
    2333          12 :   GEN ell = gel(ellg, 1), g_ell = gel(ellg, 2);
    2334          12 :   GEN H = K_get_H(K), z_f, vz_f, z = cgetg(h+1, t_VEC);
    2335             : 
    2336          12 :   z_f = Fp_pow(g_ell, diviuexact(subiu(ell, 1), f), ell);
    2337          12 :   vz_f = Fp_powers(z_f, f-1, ell)+1;
    2338        3468 :   for (i=1; i<=h; i++) gel(z, i) = gel(vz_f, Fl_mul(H[i], gk, f));
    2339          12 :   return gerepilecopy(av, FpV_roots_to_pol(z, ell, 0));
    2340             : }
    2341             : 
    2342             : static GEN
    2343         114 : G_K_l(GEN K, GEN ellg, ulong gk)
    2344             : {
    2345         114 :   pari_sp av = avma;
    2346         114 :   ulong ell = itou(gel(ellg, 1)), g_ell = itou(gel(ellg, 2));
    2347         114 :   ulong f = K_get_f(K), h = K_get_h(K), i, z_f;
    2348         114 :   GEN H = K_get_H(K), vz_f, z = cgetg(h+1, t_VEC);
    2349             : 
    2350         114 :   z_f = Fl_powu(g_ell, (ell-1) / f, ell);
    2351         114 :   vz_f = Fl_powers(z_f, f-1, ell)+1;
    2352       26898 :   for (i=1; i<=h; i++) z[i] = vz_f[Fl_mul(H[i], gk, f)];
    2353         114 :   return gerepilecopy(av, Flv_roots_to_pol(z, ell, 0));
    2354             : }
    2355             : 
    2356             : static GEN
    2357           6 : vz_vell(long d, GEN vellg, GEN *pM)
    2358             : {
    2359           6 :   long i, l = lg(vellg);
    2360           6 :   GEN A = cgetg(l, t_VEC), P = cgetg(l, t_VEC), z;
    2361             : 
    2362          18 :   for (i = 1; i < l; i++)
    2363             :   {
    2364          12 :     GEN ell = gmael(vellg, i, 1), g_ell = gmael(vellg, i, 2);
    2365          12 :     gel(A, i) = Fp_pow(g_ell, diviuexact(subiu(ell, 1), d), ell);
    2366          12 :     gel(P, i) = ell;
    2367             :   }
    2368           6 :   z = ZV_chinese(A, P, pM); return Fp_powers(z, d-1, *pM);
    2369             : }
    2370             : 
    2371             : static GEN
    2372           0 : D_xi_el_vell_FFT(GEN K, GEN elg, GEN vellg, ulong d, ulong j0, GEN vG_K)
    2373             : {
    2374           0 :   pari_sp av = avma;
    2375           0 :   ulong d_K, h = K_get_h(K), d_chi = K_get_dchi(K);
    2376           0 :   ulong el = elg[1], g_el = elg[2], el_1 = el-1;
    2377             :   ulong i, j, i2, k, dwel;
    2378           0 :   GEN u = cgetg(el+2, t_POL) , v = cgetg(h+3, t_POL);
    2379           0 :   GEN w = cgetg(el+1, t_VEC), ww;
    2380           0 :   GEN M, vz_el, G_K, z = const_vec(d_chi, gen_1);
    2381           0 :   GEN e_chi = get_e_chi(K, j0, d, &d_K);
    2382             : 
    2383           0 :   vz_el = vz_vell(el, vellg, &M);
    2384           0 :   u[1] = evalsigne(1) | evalvarn(0);
    2385           0 :   v[1] = evalsigne(1) | evalvarn(0);
    2386             : 
    2387           0 :   for (i=i2=0; i<el; i++)
    2388             :   {
    2389           0 :     ulong j2 = i2?el-i2:i2; /* i2=(i*i)%el */
    2390           0 :     gel(u, 2+i) = gel(vz_el, 1+j2);
    2391           0 :     if ((i2+=i+i+1)>=el) i2%=el;
    2392             :   }
    2393           0 :   for (k=0; k<d_K; k++)
    2394             :   {
    2395           0 :     pari_sp av = avma;
    2396             :     pari_timer ti;
    2397             :     long gd, gi;
    2398           0 :     GEN x1 = gen_1;
    2399           0 :     G_K = gel(vG_K, 1+k);
    2400           0 :     for (i=i2=0; i<=h; i++)
    2401             :     {
    2402           0 :       gel(v, 2+i) = Fp_mul(gel(G_K, 2+i), gel(vz_el, 1+i2), M);
    2403           0 :       if ((i2+=i+i+1)>=el) i2%=el;
    2404             :     }
    2405           0 :     if (DEBUGLEVEL>2) timer_start(&ti);
    2406           0 :     ww = ZX_mul(u, v);
    2407           0 :     if (DEBUGLEVEL>2)
    2408           0 :       timer_printf(&ti, "ZX_mul:%ld/%ld h*el=%ld*%ld", k, d_K, h, el);
    2409           0 :     dwel = degpol(ww)-el;
    2410           0 :     for (i=0; i<=dwel; i++) gel(w, 1+i) = addii(gel(ww, 2+i), gel(ww, 2+i+el));
    2411           0 :     for (; i<el; i++) gel(w, 1+i) = gel(ww, 2+i);
    2412           0 :     for (i=i2=1; i<el; i++)  /* w[i]=G_K(z_el^(2*i)) */
    2413             :     {
    2414           0 :       gel(w, i) = Fp_mul(gel(w, 1+i), gel(vz_el, 1+i2), M);
    2415           0 :       if ((i2+=i+i+1)>=el) i2%=el;
    2416             :     }
    2417           0 :     gd = Fl_powu(g_el, d, el);  /* a bit faster */
    2418           0 :     gi = g_el;
    2419           0 :     for (i=1; i<d; i++)
    2420             :     {
    2421           0 :       GEN xi = gen_1;
    2422           0 :       long gdi = gi;
    2423           0 :       for (j=0; i+j<el_1; j+=d)
    2424             :       {
    2425           0 :         xi = Fp_mul(xi, gel(w, (gdi+gdi)%el), M);
    2426           0 :         gdi = Fl_mul(gdi, gd, el);
    2427             :       }
    2428           0 :       gi = Fl_mul(gi, g_el, el);
    2429           0 :       xi = Fp_powu(xi, i, M);
    2430           0 :       x1 = Fp_mul(x1, xi, M);
    2431             :     }
    2432           0 :     for (i=1; i<=d_chi; i++)
    2433             :     {
    2434           0 :       GEN x2 = Fp_powu(x1, e_chi[(k+i-1)%d_K], M);
    2435           0 :       gel(z, i) = Fp_mul(gel(z, i), x2, M);
    2436             :     }
    2437           0 :     z = gerepilecopy(av, z);
    2438             :   }
    2439           0 :   return gerepilecopy(av, z);
    2440             : }
    2441             : 
    2442             : static GEN
    2443           0 : D_xi_el_vell(GEN K, GEN elg, GEN vellg, ulong d, ulong j0)
    2444             : {
    2445           0 :   pari_sp av = avma;
    2446           0 :   GEN H = K_get_H(K);
    2447           0 :   ulong f = K_get_f(K), h = K_get_h(K), g_K = K_get_g(K);
    2448             :   GEN z_f, z_el, vz_f, vz_el;
    2449           0 :   ulong el = elg[1], g_el = elg[2], el_1 = el-1;
    2450           0 :   ulong i, j, k, d_K, lv = lg(vellg), d_chi = K_get_dchi(K);
    2451           0 :   GEN A, B, P, M, z = const_vec(d_chi, gen_1);
    2452           0 :   GEN e_chi = get_e_chi(K, j0, d, &d_K);
    2453             : 
    2454           0 :   A=cgetg(lv, t_VEC);
    2455           0 :   B=cgetg(lv, t_VEC);
    2456           0 :   P=cgetg(lv, t_VEC);
    2457           0 :   for (i = 1; i < lv; i++)
    2458             :   {
    2459           0 :     GEN ell = gmael(vellg, i, 1), g_ell = gmael(vellg, i, 2);
    2460           0 :     GEN ell_1 = subiu(ell, 1);
    2461           0 :     gel(A, i) = Fp_pow(g_ell, diviuexact(ell_1, f), ell);
    2462           0 :     gel(B, i) = Fp_pow(g_ell, diviuexact(ell_1, el), ell);
    2463           0 :     gel(P, i) = ell;
    2464             :   }
    2465           0 :   z_f = ZV_chinese(A, P, &M);
    2466           0 :   z_el = ZV_chinese(B, P, NULL);
    2467           0 :   vz_f = Fp_powers(z_f, f-1, M);
    2468           0 :   vz_el = Fp_powers(z_el, el-1, M);
    2469           0 :   for (k = 0; k < d_K; k++)
    2470             :   {
    2471           0 :     pari_sp av = avma;
    2472           0 :     GEN x0 = gen_1;
    2473           0 :     long gk = Fl_powu(g_K, k, f);
    2474           0 :     for (i=1; i<el_1; i++)
    2475             :     {
    2476           0 :       long gi = Fl_powu(g_el, i, el);
    2477           0 :       GEN x1 = gen_1;
    2478           0 :       GEN x2 = gel(vz_el, 1+gi);
    2479           0 :       for (j=1; j<=h; j++)
    2480             :       {
    2481           0 :         long y = Fl_mul(H[j], gk, f);
    2482           0 :         x1 = Fp_mul(x1, Fp_sub(x2, gel(vz_f, 1+y), M), M);
    2483             :       }
    2484           0 :       x1 = Fp_powu(x1, i%d, M);
    2485           0 :       x0 = Fp_mul(x0, x1, M);
    2486             :     }
    2487           0 :     for (i=1; i<=d_chi; i++)
    2488             :     {
    2489           0 :       GEN x2 = Fp_powu(x0, e_chi[(k+i-1)%d_K], M);
    2490           0 :       gel(z, i) = Fp_mul(gel(z, i), x2, M);
    2491             :     }
    2492           0 :     z = gerepilecopy(av, z);
    2493             :   }
    2494           0 :   return gerepilecopy(av, z);
    2495             : }
    2496             : 
    2497             : static GEN
    2498          34 : D_xi_el_Flx_mul(GEN K, GEN elg, GEN ellg, GEN vG_K, ulong d, ulong j0)
    2499             : {
    2500          34 :   pari_sp av = avma;
    2501          34 :   ulong d_K, f = K_get_f(K), h = K_get_h(K), g_K = K_get_g(K);
    2502          34 :   ulong el = elg[1], g_el = elg[2], el_1 = el-1, d_chi = K_get_dchi(K);
    2503          34 :   ulong ell = itou(gel(ellg, 1)), g_ell = itou(gel(ellg, 2)), z_el;
    2504          34 :   GEN u = cgetg(el+2, t_VECSMALL), v = cgetg(h+3, t_VECSMALL);
    2505          34 :   GEN w = cgetg(el+1, t_VECSMALL), ww;
    2506          34 :   GEN vz_el, G_K, z = const_vecsmall(d_chi, 1);
    2507          34 :   GEN e_chi = get_e_chi(K, j0, d, &d_K);
    2508             :   ulong i, j, i2, k, dwel;
    2509             : 
    2510          34 :   u[1] = evalvarn(0);
    2511          34 :   v[1] = evalvarn(0);
    2512          34 :   z_el = Fl_powu(g_ell, (ell - 1) / el, ell);
    2513          34 :   vz_el = Fl_powers(z_el, el_1, ell)+1;
    2514             : 
    2515      700376 :   for (i=i2=0; i<el; i++)
    2516             :   {
    2517      700342 :     ulong j2 = i2?el-i2:i2;
    2518      700342 :     u[2+i] = vz_el[j2];
    2519      700342 :     if ((i2+=i+i+1)>=el) i2%=el;  /* i2=(i*i)%el */
    2520             :   }
    2521         694 :   for (k=0; k<d_K; k++)
    2522             :   {
    2523         660 :     pari_sp av = avma;
    2524             :     pari_timer ti;
    2525         660 :     ulong gk = Fl_powu(g_K, k, f);
    2526         660 :     long gd, gi, x1 = 1;
    2527         660 :     if (DEBUGLEVEL>2) timer_start(&ti);
    2528         660 :     G_K = (vG_K==NULL)?G_K_l(K, ellg, gk):ZX_to_Flx(gel(vG_K, 1+k), ell);
    2529         660 :     if (DEBUGLEVEL>2) timer_printf(&ti, "G_K_l");
    2530       39024 :     for (i=i2=0; i<=h; i++)
    2531             :     {
    2532       38364 :       v[2+i] = Fl_mul(G_K[2+i], vz_el[i2], ell);
    2533       38364 :       if ((i2+=i+i+1)>=el) i2%=el;  /* i2=(i*i)%el */
    2534             :     }
    2535         660 :     if (DEBUGLEVEL>2) timer_start(&ti);
    2536         660 :     ww = Flx_mul(u, v, ell);
    2537         660 :     if (DEBUGLEVEL>2)
    2538           0 :       timer_printf(&ti, "Flx_mul:%ld/%ld h*el=%ld*%ld", k, d_K, h, el);
    2539         660 :     dwel=degpol(ww)-el; /* dwel=h-1 */
    2540       38364 :     for (i=0; i<=dwel; i++) w[1+i] = Fl_add(ww[2+i], ww[2+i+el], ell);
    2541     8388480 :     for (; i<el; i++) w[1+i] = ww[2+i];
    2542     8425524 :     for (i=i2=1; i<el; i++)  /* w[i]=G_K(z_el^(2*i)) */
    2543             :     {
    2544     8424864 :       w[i] = Fl_mul(w[1+i], vz_el[i2], ell);
    2545     8424864 :       if ((i2+=i+i+1)>=el) i2%=el;  /* i2=(i*i)%el */
    2546             :     }
    2547         660 :     gd = Fl_powu(g_el, d, el);  /* a bit faster */
    2548         660 :     gi = g_el;
    2549        4596 :     for (i=1; i<d; i++)
    2550             :     {
    2551        3936 :       long xi = 1, gdi = gi;
    2552     8163696 :       for (j=0; i+j<el_1; j+=d)
    2553             :       {
    2554     8159760 :         xi = Fl_mul(xi, w[(gdi+gdi)%el], ell);
    2555     8159760 :         gdi = Fl_mul(gdi, gd, el);
    2556             :       }
    2557        3936 :       gi = Fl_mul(gi, g_el, el);
    2558        3936 :       xi = Fl_powu(xi, i, ell);
    2559        3936 :       x1 = Fl_mul(x1, xi, ell);
    2560             :     }
    2561        2412 :     for (i=1; i<=d_chi; i++)
    2562             :     {
    2563        1752 :       long x2 = Fl_powu(x1, e_chi[(k+i-1)%d_K], ell);
    2564        1752 :       z[i] = Fl_mul(z[i], x2, ell);
    2565             :     }
    2566         660 :     set_avma(av);
    2567             :   }
    2568          34 :   return gerepilecopy(av, Flv_to_ZV(z));
    2569             : }
    2570             : 
    2571             : static GEN
    2572          35 : D_xi_el_ZX_mul(GEN K, GEN elg, GEN ellg, GEN vG_K, ulong d, ulong j0)
    2573             : {
    2574          35 :   pari_sp av = avma;
    2575          35 :   GEN ell = gel(ellg,1), g_ell, u, v, w, ww, z_el, vz_el, G_K, z, e_chi;
    2576             :   ulong d_K, f, h, g_K, el, g_el, el_1, d_chi, i, j, i2, k, dwel;
    2577             : 
    2578          35 :   if (lgefint(ell) == 3) return D_xi_el_Flx_mul(K, elg, ellg, vG_K, d, j0);
    2579           1 :   f = K_get_f(K); h = K_get_h(K); g_K = K_get_g(K);
    2580           1 :   el = elg[1]; g_el = elg[2]; el_1 = el-1; d_chi = K_get_dchi(K);
    2581           1 :   g_ell = gel(ellg, 2);
    2582           1 :   z = const_vec(d_chi, gen_1);
    2583           1 :   e_chi = get_e_chi(K, j0, d, &d_K);
    2584             : 
    2585           1 :   u = cgetg(el+2,t_POL); u[1] = evalsigne(1) | evalvarn(0);
    2586           1 :   v = cgetg(h+3, t_POL); v[1] = evalsigne(1) | evalvarn(0);
    2587           1 :   w = cgetg(el+1, t_VEC);
    2588           1 :   z_el = Fp_pow(g_ell, diviuexact(subiu(ell, 1), el), ell);
    2589           1 :   vz_el = Fp_powers(z_el, el_1, ell)+1;
    2590             : 
    2591      114998 :   for (i=i2=0; i<el; i++)
    2592             :   {
    2593      114997 :     ulong j2 = i2?el-i2:i2; /* i2=(i*i)%el */
    2594      114997 :     gel(u, 2+i) = gel(vz_el, j2);
    2595      114997 :     if ((i2+=i+i+1)>=el) i2%=el;
    2596             :   }
    2597          13 :   for (k=0; k<d_K; k++)
    2598             :   {
    2599          12 :     pari_sp av = avma;
    2600             :     pari_timer ti;
    2601          12 :     long gd, gi, gk = Fl_powu(g_K, k, f);
    2602          12 :     GEN x1 = gen_1;
    2603          12 :     if (DEBUGLEVEL>2) timer_start(&ti);
    2604          12 :     G_K = (vG_K==NULL) ? G_K_p(K, ellg, gk):RgX_to_FpX(gel(vG_K, 1+k), ell);
    2605          12 :     if (DEBUGLEVEL>2) timer_printf(&ti, "G_K_p");
    2606        3480 :     for (i=i2=0; i<=h; i++)
    2607             :     {
    2608        3468 :       gel(v, 2+i) = Fp_mul(gel(G_K, 2+i), gel(vz_el, i2), ell);
    2609        3468 :       if ((i2+=i+i+1)>=el) i2%=el;
    2610             :     }
    2611          12 :     if (DEBUGLEVEL>2) timer_start(&ti);
    2612          12 :     ww = ZX_mul(u, v);
    2613          12 :     if (DEBUGLEVEL>2)
    2614           0 :       timer_printf(&ti, "ZX_mul:%ld/%ld h*el=%ld*%ld", k, d_K, h, el);
    2615          12 :     dwel = degpol(ww)-el;
    2616        3468 :     for (i=0; i<=dwel; i++) gel(w, 1+i) = addii(gel(ww, 2+i), gel(ww, 2+i+el));
    2617     1376520 :     for (; i<el; i++) gel(w, 1+i) = gel(ww, 2+i);
    2618     1379964 :     for (i=i2=1; i<el; i++)  /* w[i]=G_K(z_el^(2*i)) */
    2619             :     {
    2620     1379952 :       gel(w, i) = Fp_mul(gel(w, 1+i), gel(vz_el, i2), ell);
    2621     1379952 :       if ((i2+=i+i+1)>=el) i2%=el;
    2622             :     }
    2623          12 :     gd = Fl_powu(g_el, d, el);  /* a bit faster */
    2624          12 :     gi = g_el;
    2625         444 :     for (i=1; i<d; i++)
    2626             :     {
    2627         432 :       GEN xi = gen_1;
    2628         432 :       long gdi = gi;
    2629     1343088 :       for (j=0; i+j<el_1; j+=d)
    2630             :       {
    2631     1342656 :         xi = Fp_mul(xi, gel(w, (gdi+gdi)%el), ell);
    2632     1342656 :         gdi = Fl_mul(gdi, gd, el);
    2633             :       }
    2634         432 :       gi = Fl_mul(gi, g_el, el);
    2635         432 :       xi = Fp_powu(xi, i, ell);
    2636         432 :       x1 = Fp_mul(x1, xi, ell);
    2637             :     }
    2638          24 :     for (i=1; i<=d_chi; i++)
    2639             :     {
    2640          12 :       GEN x2 = Fp_powu(x1, e_chi[(k+i-1)%d_K], ell);
    2641          12 :       gel(z, i) = Fp_mul(gel(z, i), x2, ell);
    2642             :     }
    2643          12 :     z = gerepilecopy(av, z);
    2644             :   }
    2645           1 :   return gerepilecopy(av, z);
    2646             : }
    2647             : 
    2648             : static GEN
    2649           0 : D_xi_el_ss(GEN K, GEN elg, GEN ellg, ulong d, ulong j0)
    2650             : {
    2651           0 :   pari_sp av = avma;
    2652           0 :   GEN H = K_get_H(K);
    2653           0 :   ulong d_K, f = K_get_f(K), h = K_get_h(K), g_K = K_get_g(K);
    2654           0 :   ulong el = elg[1], g_el = elg[2], el_1 = el-1;
    2655           0 :   ulong ell = itou(gel(ellg, 1)), g_ell = itou(gel(ellg, 2));
    2656           0 :   ulong i, j, k, gk, z_f, z_el, d_chi = K_get_dchi(K);
    2657           0 :   GEN vz_f, vz_el, z = const_vecsmall(d_chi, 1);
    2658           0 :   GEN e_chi = get_e_chi(K, j0, d, &d_K);
    2659             : 
    2660           0 :   z_f = Fl_powu(g_ell, (ell - 1) / f, ell);
    2661           0 :   z_el = Fl_powu(g_ell, (ell - 1) / el, ell);
    2662           0 :   vz_f = Fl_powers(z_f, f-1, ell)+1;
    2663           0 :   vz_el = Fl_powers(z_el, el-1, ell)+1;
    2664           0 :   gk = 1; /* g_K^k */
    2665           0 :   for (k = 0; k < d_K; k++)
    2666             :   {
    2667           0 :     ulong x0 = 1, gi = g_el; /* g_el^i */
    2668           0 :     for (i = 1; i < el_1; i++)
    2669             :     {
    2670           0 :       ulong x1 = 1, x2 = vz_el[gi];
    2671           0 :       for (j=1; j<=h; j++)
    2672             :       {
    2673           0 :         ulong y = Fl_mul(H[j], gk, f);
    2674           0 :         x1 = Fl_mul(x1, Fl_sub(x2, vz_f[y], ell), ell);
    2675             :       }
    2676           0 :       x1 = Fl_powu(x1, i%d, ell);
    2677           0 :       x0 = Fl_mul(x0, x1, ell);
    2678           0 :       gi = Fl_mul(gi, g_el, el);
    2679             :     }
    2680           0 :     for (i = 1; i <= d_chi; i++)
    2681             :     {
    2682           0 :       ulong x2 = Fl_powu(x0, e_chi[(k+i-1)%d_K], ell);
    2683           0 :       z[i] = Fl_mul(z[i], x2, ell);
    2684             :     }
    2685           0 :     gk = Fl_mul(gk, g_K, f);
    2686             :   }
    2687           0 :   return gerepileupto(av, Flv_to_ZV(z));
    2688             : }
    2689             : 
    2690             : static GEN
    2691           0 : D_xi_el_sl(GEN K, GEN elg, GEN ellg, ulong d, ulong j0)
    2692             : {
    2693           0 :   pari_sp av = avma;
    2694           0 :   GEN ell = gel(ellg, 1), H;
    2695             :   GEN g_ell, ell_1, z_f, z_el, vz_f, vz_el, z, e_chi;
    2696             :   ulong d_K, f, h, g_K, el, g_el, el_1, d_chi, i, j, k, gk;
    2697             : 
    2698           0 :   if (lgefint(ell) == 3) return D_xi_el_ss(K, elg, ellg, d, j0);
    2699           0 :   H = K_get_H(K);
    2700           0 :   f = K_get_f(K); h = K_get_h(K); g_K = K_get_g(K);
    2701           0 :   el = elg[1]; g_el = elg[2]; el_1 = el-1; d_chi = K_get_dchi(K);
    2702           0 :   g_ell = gel(ellg, 2); ell_1 = subiu(ell, 1);
    2703           0 :   z = const_vec(d_chi, gen_1);
    2704           0 :   e_chi = get_e_chi(K, j0, d, &d_K);
    2705             : 
    2706           0 :   z_f = Fp_pow(g_ell, diviuexact(ell_1, f), ell);
    2707           0 :   z_el = Fp_pow(g_ell, diviuexact(ell_1, el), ell);
    2708           0 :   vz_f = Fp_powers(z_f, f-1, ell) + 1;
    2709           0 :   vz_el = Fp_powers(z_el, el-1, ell) + 1;
    2710           0 :   gk = 1; /* g_K^k */
    2711           0 :   for (k = 0; k < d_K; k++)
    2712             :   {
    2713           0 :     pari_sp av2 = avma;
    2714           0 :     GEN x0 = gen_1;
    2715           0 :     ulong gi = g_el; /* g_el^i */
    2716           0 :     for (i = 1; i < el_1; i++)
    2717             :     {
    2718           0 :       pari_sp av3 = avma;
    2719           0 :       GEN x1 = gen_1, x2 = gel(vz_el, gi);
    2720           0 :       for (j = 1; j <= h; j++)
    2721             :       {
    2722           0 :         ulong y = Fl_neg(Fl_mul(H[j], gk, f), f);
    2723           0 :         x1 = Fp_mul(x1, Fp_sub(x2, gel(vz_f, y), ell), ell);
    2724             :       }
    2725           0 :       x1 = Fp_powu(x1, i%d, ell);
    2726           0 :       x0 = gerepileuptoint(av3, Fp_mul(x0, x1, ell));
    2727           0 :       gi = Fl_mul(gi, g_el, el);
    2728             :     }
    2729           0 :     for (i=1; i<=d_chi; i++)
    2730             :     {
    2731           0 :       GEN x2 = Fp_powu(x0, e_chi[(k+i-1)%d_K], ell);
    2732           0 :       gel(z, i) = Fp_mul(gel(z, i), x2, ell);
    2733             :     }
    2734           0 :     if (k == d_K-1) break;
    2735           0 :     z = gerepilecopy(av2, z);
    2736           0 :     gk = Fl_mul(gk, g_K, f);
    2737             :   }
    2738           0 :   return gerepilecopy(av, z);
    2739             : }
    2740             : 
    2741             : static long
    2742         175 : get_y(GEN z, GEN ellg, long d)
    2743             : {
    2744         175 :   GEN ell = gel(ellg, 1), g_ell = gel(ellg, 2);
    2745         175 :   GEN elld = diviuexact(subiu(ell, 1), d);
    2746         175 :   GEN g_elld = Fp_pow(g_ell, elld, ell);
    2747         175 :   GEN x = Fp_pow(modii(z, ell), elld, ell);
    2748             :   long k;
    2749      211778 :   for (k=0; k<d; k++)
    2750             :   {
    2751      211778 :     if (equali1(x)) break;
    2752      211603 :     x = Fp_mul(x, g_elld, ell);
    2753             :   }
    2754         175 :   if (k==0) k=d;
    2755         161 :   else if (d<=k) pari_err_BUG("subcyclopclgp [MLL]");
    2756         175 :   return k;
    2757             : }
    2758             : 
    2759             : static void
    2760           0 : real_MLLn(long *y, GEN K, ulong p, ulong d_pow, ulong n,
    2761             :     GEN velg, GEN vellg, GEN vG_K, ulong j0)
    2762             : {
    2763           0 :   pari_sp av = avma;
    2764           0 :   ulong i, j, k, d = upowuu(p, d_pow), h = gmael(K, 1, 2)[3];
    2765           0 :   ulong row = lg(vellg)-1;
    2766           0 :   for (i=1; i<=n; i++)
    2767             :   {
    2768           0 :     GEN elg = gel(velg, i), z;
    2769           0 :     ulong el = elg[1], nz;
    2770             :     pari_timer ti;
    2771           0 :     if (DEBUGLEVEL>1) timer_start(&ti);
    2772           0 :     z = (h<el) ? D_xi_el_vell_FFT(K, elg, vellg, d, j0, vG_K)
    2773           0 :                : D_xi_el_vell(K, elg, vellg, d, j0);
    2774           0 :     if (DEBUGLEVEL>1) timer_printf(&ti, "subcyclopclgp:[D_xi_el]");
    2775           0 :     if (DEBUGLEVEL>2) err_printf("z=%Ps\n", z);
    2776           0 :     nz = lg(z)-1;
    2777           0 :     for (k = 1; k <= nz; k++)
    2778           0 :       for (j=1; j<=row; j++)
    2779           0 :         y[(j-1)*row+(i-1)*nz+k-1] = get_y(gel(z, k), gel(vellg, j), d);
    2780           0 :     set_avma(av);
    2781             :   }
    2782           0 : }
    2783             : 
    2784             : static void
    2785          14 : real_MLL1(long *y, GEN K, ulong p, ulong d_pow, GEN velg, GEN vellg, ulong j0)
    2786             : {
    2787          14 :   ulong h = gmael(K, 1, 2)[3], d = upowuu(p, d_pow);
    2788          14 :   GEN elg = gel(velg, 1), ellg = gel(vellg, 1), z;
    2789          14 :   ulong el = elg[1];
    2790             :   pari_timer ti;
    2791             : 
    2792          14 :   if (DEBUGLEVEL>2) timer_start(&ti);
    2793          14 :   z = h < el? D_xi_el_ZX_mul(K, elg, ellg, NULL, d, j0)
    2794          14 :             : D_xi_el_sl(K, elg, ellg, d, j0);
    2795          14 :   if (DEBUGLEVEL>2) timer_printf(&ti, "subcyclopclgp:[D_xi_el]");
    2796          14 :   if (DEBUGLEVEL>2) err_printf("z=%Ps\n", z);
    2797          14 :   y[0] = get_y(gel(z, 1), ellg, d);
    2798          14 : }
    2799             : 
    2800             : static void
    2801           7 : real_MLL(long *y, GEN K, long p, long d_pow, long n,
    2802             :     GEN velg, GEN vellg, GEN vG_K, long j0)
    2803             : {
    2804           7 :   long i, j, row = lg(vellg)-1;
    2805           7 :   ulong k, h = gmael(K, 1, 2)[3], d = upowuu(p, d_pow);
    2806          28 :   for (j=1; j<=row; j++)
    2807             :   {
    2808          21 :     GEN ellg = gel(vellg, j);
    2809          42 :     for (i=1; i<=n; i++)
    2810             :     {
    2811          21 :       pari_sp av2 = avma;
    2812          21 :       GEN elg = gel(velg, i), z;
    2813          21 :       ulong el = elg[1], nz;
    2814             :       pari_timer ti;
    2815          21 :       if (DEBUGLEVEL>2) timer_start(&ti);
    2816          21 :       z = h < el? D_xi_el_ZX_mul(K, elg, ellg, vG_K, d, j0)
    2817          21 :                 : D_xi_el_sl(K, elg, ellg, d, j0);
    2818          21 :       if (DEBUGLEVEL>2) timer_printf(&ti, "subcyclopclgp:[D_xi_el]");
    2819          21 :       if (DEBUGLEVEL>3) err_printf("z=%Ps\n", z);
    2820          21 :       nz = lg(z)-1;
    2821          84 :       for (k = 1; k <= nz; k++)
    2822          63 :         y[(j-1)*row+(i-1)*nz+k-1] = get_y(gel(z, k), ellg, d);
    2823          21 :       set_avma(av2);
    2824             :     }
    2825             :   }
    2826           7 : }
    2827             : 
    2828             : static long
    2829          21 : use_basis(long d_K, long f) { return (d_K<=10 || (d_K<=30 && f<=5000)); }
    2830             : 
    2831             : static long
    2832           7 : use_factor(ulong f)
    2833           7 : { GEN fa = factoru(f), P = gel(fa, 1); return (P[lg(P)-1]<500); }
    2834             : 
    2835             : /* group structure, destroy gr */
    2836             : static GEN
    2837          63 : get_str(GEN gr)
    2838             : {
    2839          63 :   GEN z = gel(gr,2);
    2840          63 :   long i, j, l = lg(z);
    2841         154 :   for (i = j = 1; i < l; i++)
    2842          91 :     if (lgefint(gel(z, i)) > 2) gel(z,j++) = gel(z,i);
    2843          63 :   setlg(z, j); return z;
    2844             : }
    2845             : 
    2846             : static GEN
    2847          21 : cyc_real_MLL(GEN K, ulong p, long d_pow, long j0, long flag)
    2848             : {
    2849          21 :   ulong d_K = K_get_d(K), f = K_get_f(K), d_chi = K_get_dchi(K);
    2850          21 :   ulong n, n0 = 1, f0, n_el = d_pow, d = upowuu(p, d_pow), rank = n_el*d_chi;
    2851          21 :   GEN velg = const_vec(n_el, NULL), vellg = NULL;
    2852          21 :   GEN oldgr = mkvec2(gen_0, NULL), newgr = mkvec2(gen_0, NULL);
    2853          21 :   long *y0 = (long*)stack_calloc(sizeof(long)*rank*rank);
    2854             : 
    2855          21 :   if (DEBUGLEVEL>1)
    2856           0 :     err_printf("cyc_real_MLL:p=%ld d_pow=%ld deg(K)=%ld cond(K)=%ld g_K=%ld\n",
    2857             :         p, d_pow, d_K, f, K_get_g(K));
    2858          21 :   gel(K, 2) = get_chi(gel(K,1));
    2859          21 :   if (f-1 <= (d_K<<1)) flag |= USE_F;
    2860          21 :   else if (use_basis(d_K, f)) flag |= USE_BASIS;
    2861           7 :   else if (use_factor(f)) flag |= USE_FACTOR;
    2862           0 :   else flag |= USE_GALOIS_POL;
    2863          21 :   if (flag&USE_BASIS) K = vec_append(K, xi_data_basis(K));
    2864           7 :   else if (flag&USE_GALOIS_POL) K = vec_append(K, xi_data_galois(K));
    2865          21 :   f0 = f%p?f:f/p;
    2866          21 :   gel(velg, 1) = next_el_real(K, p, d_pow, mkvecsmall2(1, 1), j0, flag);
    2867          21 :   if (flag&USE_FULL_EL)
    2868             :   {
    2869           0 :     for (n=2; n<=n_el; n++)
    2870           0 :       gel(velg, n) = next_el_real(K, p, d_pow, gel(velg, n+1), j0, flag);
    2871           0 :     n0 = n_el;
    2872             :   }
    2873             : 
    2874          21 :   for (n=n0; n<=n_el; n++) /* loop while structure is unknown */
    2875             :   {
    2876          21 :     pari_sp av2 = avma;
    2877             :     long n_ell, m, M;
    2878             :     GEN y;
    2879             :     pari_timer ti;
    2880          21 :     if (DEBUGLEVEL>2) timer_start(&ti);
    2881          21 :     vellg = set_ell_real(K, velg, n, d_chi, d*d, f0, j0);
    2882          21 :     n_ell = lg(vellg) -1; /* equal to n*d_chi */
    2883          21 :     if (DEBUGLEVEL>2) timer_printf(&ti, "set_ell_real");
    2884          21 :     if (DEBUGLEVEL>3) err_printf("vel=%Ps\nvell=%Ps\n", velg, vellg);
    2885          21 :     if (n_ell==1)
    2886          14 :       real_MLL1(y0, K, p, d_pow, velg, vellg, j0);
    2887             :     else
    2888             :     {
    2889             :       GEN vG_K;
    2890           7 :       if (DEBUGLEVEL>2) timer_start(&ti);
    2891           7 :       vG_K = make_G_K(K, vellg);
    2892           7 :       if (DEBUGLEVEL>2) timer_printf(&ti, "make_G_K");
    2893           7 :       if (lgefint(gmael(vellg, n_ell, 1))<=3 || (flag&SAVE_MEMORY))
    2894           7 :         real_MLL(y0, K, p, d_pow, n, velg, vellg, vG_K, j0);
    2895             :       else
    2896           0 :         real_MLLn(y0, K, p, d_pow, n, velg, vellg, vG_K, j0);
    2897             :     }
    2898          21 :     set_avma(av2);
    2899          21 :     y = ary2mat(y0, n_ell);
    2900          21 :     if (DEBUGLEVEL>3) err_printf("y=%Ps\n", y);
    2901          21 :     y = ZM_snf(y);
    2902          21 :     if (DEBUGLEVEL>3) err_printf("y=%Ps\n", y);
    2903          21 :     y = make_p_part(y, p, d_pow);
    2904          21 :     if (DEBUGLEVEL>3) err_printf("y=%Ps\n", y);
    2905          21 :     newgr = structure_MLL(y, d_pow);
    2906          21 :     if (DEBUGLEVEL>3)
    2907           0 :       err_printf("d_pow=%ld d_chi=%ld old=%Ps new=%Ps\n",d_pow,d_chi,oldgr,newgr);
    2908          21 :     if (equalsi(d_pow*d_chi, gel(newgr, 1))) break;
    2909           0 :     if ((m = find_del_el(&oldgr, newgr, n, n_el, d_chi)))
    2910           0 :     { M = m = delete_el(velg, m); n--; }
    2911             :     else
    2912           0 :     { M = n+1; m = n; }
    2913           0 :     gel(velg, M) = next_el_real(K, p, d_pow, gel(velg, m), j0, flag);
    2914             :   }
    2915          21 :   return get_str(newgr);
    2916             : }
    2917             : 
    2918             : static GEN
    2919           0 : cyc_buch(long dK, GEN p, long d_pow)
    2920             : {
    2921           0 :   GEN z = Buchquad(stoi(dK), 0.0, 0.0, 0), cyc = gel(z,2);
    2922           0 :   long i, l = lg(cyc);
    2923           0 :   if (Z_pval(gel(z,1), p) != d_pow) pari_err_BUG("subcyclopclgp [Buchquad]");
    2924           0 :   for (i = 1; i < l; i++)
    2925             :   {
    2926           0 :     long x = Z_pval(gel(cyc, i), p); if (!x) break;
    2927           0 :     gel(cyc, i) = utoipos(x);
    2928             :   }
    2929           0 :   setlg(cyc, i); return cyc;
    2930             : }
    2931             : 
    2932             : static void
    2933           0 : verbose_output(GEN K, GEN p, long pow, long j)
    2934             : {
    2935           0 :   long d = K_get_d(K), f = K_get_f(K), s = K_get_s(K), d_chi = K_get_dchi(K);
    2936           0 :   err_printf("|A_K_psi|=%Ps^%ld, psi=chi^%ld, d_psi=%ld, %s,\n\
    2937             :     [K:Q]=%ld, [f,H]=[%ld, %Ps]\n",
    2938           0 :     p,pow*d_chi,j,d_chi,s?"real":"imaginary",d,f,zv_to_ZV(gmael3(K,1,1,1)));
    2939           0 : }
    2940             : 
    2941             : static int
    2942       35091 : cyc_real_pre(GEN K, GEN xi, ulong p, ulong j, long el)
    2943             : {
    2944       35091 :   pari_sp av = avma;
    2945       35091 :   ulong i, d_K, x = 1;
    2946       35091 :   GEN e_chi = get_e_chi(K, j, p, &d_K);
    2947             : 
    2948       35091 :   xi++;
    2949     1254827 :   for (i = 0; i < d_K; i++) x = Fl_mul(x, Fl_powu(xi[i], e_chi[i], el), el);
    2950       35091 :   return gc_ulong(av, Fl_powu(x, (el-1)/p, el));
    2951             : }
    2952             : 
    2953             : /* return vec[-1,[],0], vec[0,[],0], vec[1,[1],0], vec[2,[1,1],0] etc */
    2954             : static GEN
    2955       26579 : cyc_real_ss(GEN K, GEN xi, ulong p, long j, long pow, long el, ulong pn, long flag)
    2956             : {
    2957       26579 :   ulong d_chi = K_get_dchi(K);
    2958       26579 :   if (cyc_real_pre(K, xi, pn, j, el) == 1) return NULL; /* not determined */
    2959       21273 :   if (--pow==0) return mkvec3(gen_0, nullvec(), gen_0); /* trivial */
    2960         105 :   if (DEBUGLEVEL) verbose_output(K, utoi(p), pow, j);
    2961         105 :   if (flag&USE_MLL)
    2962             :   {
    2963          21 :     pari_sp av = avma;
    2964          21 :     GEN gr = (K_get_d(K) == 2)? cyc_buch(K_get_f(K), utoi(p), pow)
    2965          21 :                                : cyc_real_MLL(K, p, pow, j, flag);
    2966          21 :     return gerepilecopy(av, mkvec3(utoipos(d_chi * pow), gr, gen_0));
    2967             :   }
    2968          84 :   if (pow==1) return mkvec3(utoi(d_chi), onevec(d_chi), gen_0);
    2969          21 :   return mkvec3(utoi(pow*d_chi), nullvec(), gen_0);
    2970             : }
    2971             : 
    2972             : static GEN
    2973        5145 : cyc_real_ll(GEN K, GEN xi, GEN p, long j, long pow, GEN el, GEN pn, long flag)
    2974             : {
    2975        5145 :   pari_sp av = avma;
    2976        5145 :   ulong i, d_K = K_get_d(K), d_chi = K_get_dchi(K);
    2977        5145 :   GEN e_chi = get_e_chi_ll(K, j, pn), x = gen_1;
    2978             : 
    2979        5145 :   xi++;
    2980      246785 :   for (i = 0; i < d_K; i++)
    2981      241640 :     x = Fp_mul(x, Fp_pow(gel(xi, i), gel(e_chi, i), el), el);
    2982        5145 :   x = Fp_pow(x, diviiexact(subiu(el, 1), pn), el); /* x = x^(el-1)/pn mod el */
    2983        5145 :   set_avma(av); if (equali1(x)) return NULL; /* not determined */
    2984        5145 :   if (--pow==0) return mkvec3(gen_0, nullvec(), gen_0); /* trivial */
    2985           0 :   if (DEBUGLEVEL) verbose_output(K, p, pow, j);
    2986           0 :   if (flag&USE_MLL)
    2987           0 :     pari_err_IMPL(stack_sprintf("flag=%ld for large prime", USE_MLL));
    2988           0 :   if (pow==1) return mkvec3(utoi(d_chi), onevec(d_chi), gen_0);
    2989           0 :   return mkvec3(utoi(pow*d_chi), nullvec(), gen_0);
    2990             : }
    2991             : 
    2992             : /* xi[1+i] = xi^(g^i), 0 <= i <= d-1 */
    2993             : static GEN
    2994       13797 : xi_conj_s(GEN K, ulong el)
    2995             : {
    2996       13797 :   pari_sp av = avma;
    2997       13797 :   GEN  H = K_get_H(K);
    2998       13797 :   long d = K_get_d(K), f = K_get_f(K), h = K_get_h(K), g = K_get_g(K);
    2999       13797 :   long i, gi = 1, z = Fl_powu(pgener_Fl(el), (el-1)/f, el);
    3000       13797 :   GEN vz = Fl_powers(z, f, el)+1, xi = cgetg(d+1, t_VECSMALL);
    3001             : 
    3002      372953 :   for (i=1; i<=d; i++)
    3003             :   {
    3004      359156 :     long j, x = 1;
    3005   170777852 :     for (j=1; j<=h; j++)
    3006   170418696 :       x = Fl_mul(x, vz[Fl_mul(H[j], gi, f)]-1, el);
    3007      359156 :     xi[i] = x;
    3008      359156 :     gi = Fl_mul(gi, g, f);
    3009             :   }
    3010       13797 :   return gerepilecopy(av, xi);
    3011             : }
    3012             : 
    3013             : static GEN
    3014        1673 : xi_conj_l(GEN K, GEN el)
    3015             : {
    3016        1673 :   pari_sp av = avma;
    3017        1673 :   GEN  H = K_get_H(K);
    3018        1673 :   long d = K_get_d(K), f = K_get_f(K), h = K_get_h(K), g = K_get_g(K);
    3019        1673 :   long i, gi = 1;
    3020        1673 :   GEN z = Fp_pow(pgener_Fp(el), diviuexact(subiu(el, 1), f), el);
    3021        1673 :   GEN vz = Fp_powers(z, f, el)+1, xi = cgetg(d+1, t_VEC);
    3022             : 
    3023       70462 :   for (i=1; i<=d; i++)
    3024             :   {
    3025             :     long j;
    3026       68789 :     GEN x = gen_1;
    3027     6756043 :     for (j=1; j<=h; j++)
    3028     6687254 :       x = Fp_mul(x, subiu(gel(vz, Fl_mul(H[j], gi, f)), 1), el);
    3029       68789 :     gel(xi, i) = x;
    3030       68789 :     gi = Fl_mul(gi, g, f);
    3031             :   }
    3032        1673 :   return gerepilecopy(av, xi);
    3033             : }
    3034             : 
    3035             : static GEN
    3036       10164 : pclgp_cyc_real(GEN K, GEN p, long max_pow, long flag)
    3037             : {
    3038       10164 :   const long NUM_EL = 20;
    3039       10164 :   GEN C = gel(K, 5);
    3040       10164 :   long f_K = K_get_f(K), n_conj = K_get_nconj(K);
    3041       10164 :   long i, pow, n_el, n_done = 0;
    3042       10164 :   GEN gr = nullvec(), Done = const_vecsmall(n_conj, 0), xi;
    3043       10164 :   long first = 1;
    3044             : 
    3045       10290 :   for (pow=1; pow<=max_pow; pow++)
    3046             :   {
    3047       10290 :     GEN pn = powiu(p, pow), fpn = muliu(pn, f_K), el = addiu(fpn, 1);
    3048      109753 :     for (n_el = 0; n_el < NUM_EL; el = addii(el, fpn))
    3049             :     {
    3050             :       ulong uel;
    3051      109627 :       if (!BPSW_psp(el)) continue;
    3052       15470 :       n_el++; uel = itou_or_0(el);
    3053       15470 :       if (uel)
    3054             :       {
    3055       13797 :         xi = xi_conj_s(K, uel);
    3056       13797 :         if (first && n_conj > 10) /* mark trivial chi-part */
    3057             :         {
    3058        8715 :           for (i = 1; i <= n_conj; i++)
    3059             :           {
    3060        8512 :             if (cyc_real_pre(K, xi, p[2], C[i], uel) == 1) continue;
    3061        8260 :             Done[i] = 1;
    3062        8260 :             if (++n_done == n_conj) return gr;
    3063             :           }
    3064         203 :           first = 0; continue;
    3065             :         }
    3066             :       }
    3067             :       else
    3068        1673 :         xi = xi_conj_l(K, el);
    3069       55356 :       for (i = 1; i <= n_conj; i++)
    3070             :       {
    3071             :         GEN z;
    3072       50253 :         if (Done[i]) continue;
    3073       31724 :         if (uel)
    3074       26579 :           z = cyc_real_ss(K, xi, p[2], C[i], pow, uel, itou(pn), flag);
    3075             :         else
    3076        5145 :           z = cyc_real_ll(K, xi, p, C[i], pow, el, pn, flag);
    3077       31724 :         if (!z) continue;
    3078       26418 :         Done[i] = 1;
    3079       26418 :         if (!isintzero(gel(z, 1))) gr = vec_append(gr, z);
    3080       26418 :         if (++n_done == n_conj) return gr;
    3081             :       }
    3082             :     }
    3083             :   }
    3084           0 :   pari_err_BUG("pclgp_cyc_real: max_pow is not enough");
    3085             :   return NULL; /*LCOV_EXCL_LINE*/
    3086             : }
    3087             : 
    3088             : /* return (el, g_el) */
    3089             : static GEN
    3090          56 : next_el_imag(GEN elg, long f)
    3091             : {
    3092          56 :   long el = elg[1];
    3093          56 :   if (odd(f)) f<<=1;
    3094         140 :   while (!uisprime(el+=f));
    3095          56 :   return mkvecsmall2(el, pgener_Fl(el));
    3096             : }
    3097             : 
    3098             : /* return (ell, g_ell) */
    3099             : static GEN
    3100          70 : next_ell_imag(GEN ellg, GEN df0l)
    3101             : {
    3102          70 :   GEN ell = gel(ellg, 1);
    3103         770 :   while (!BPSW_psp(ell = addii(ell, df0l)));
    3104          70 :   return mkvec2(ell, pgener_Fp(ell));
    3105             : }
    3106             : 
    3107             : static GEN
    3108          56 : set_ell_imag(GEN velg, long n, long d_chi, GEN df0)
    3109             : {
    3110          56 :   long i, n_ell = n*d_chi;
    3111          56 :   GEN z = cgetg(n_ell + 1, t_VEC);
    3112          56 :   GEN df0l = shifti(df0, 1), ellg = mkvec2(gen_1, gen_1);
    3113         126 :   for (i=1; i<=n; i++) df0l = muliu(df0l, gel(velg, i)[1]);
    3114         126 :   for (i=1; i<=n_ell; i++) ellg = gel(z, i)= next_ell_imag(ellg, df0l);
    3115          56 :   return z;
    3116             : }
    3117             : 
    3118             : /* U(X)=u(x)+u(X)*X^f+...+f(X)*X^((m-1)f) or u(x)-u(X)*X^f+...
    3119             :  * U(X)V(X)=u(X)V(X)(1+X^f+...+X^((m-1)f))
    3120             :  *         =w_0+w_1*X+...+w_{f+el-3}*X^(f+el-3)
    3121             :  * w_i (1 <= i <= f+el-2) are needed.
    3122             :  * w_{f+el-2}=0 if el-1 == f.
    3123             :  * W_i = w_i + w_{i+el-1} (1 <= i <= f-1). */
    3124             : static GEN
    3125          85 : gauss_Flx_mul(ulong f, GEN elg, GEN ellg)
    3126             : {
    3127          85 :   pari_sp av = avma;
    3128          85 :   ulong el = elg[1], g_el= elg[2];
    3129          85 :   ulong el_1 = el-1, f2 = f<<1, lv = el_1, lu = f, m = el_1/f;
    3130          85 :   ulong ell = itou(gel(ellg, 1)), g_ell = itou(gel(ellg, 2));
    3131          85 :   ulong z_2f = Fl_powu(g_ell, (ell - 1) / f2, ell);
    3132          85 :   ulong z_el = Fl_powu(g_ell, (ell - 1) / el, ell);
    3133             :   ulong i, i2, gi;
    3134          85 :   GEN W = cgetg(f+1, t_VECSMALL), vz_2f, vz_el;
    3135          85 :   GEN u = cgetg(lu+2, t_VECSMALL), v = cgetg(lv+2, t_VECSMALL), w0;
    3136             : 
    3137          85 :   u[1] = evalsigne(1);
    3138          85 :   v[1] = evalsigne(1);
    3139          85 :   vz_2f = Fl_powers(z_2f, f2-1, ell);
    3140          85 :   vz_el = Fl_powers(z_el, el_1, ell);
    3141      528459 :   for (i=i2=0; i<lu; i++)
    3142             :   {
    3143             :     long j2; /* i2=(i*i)%f2, gi=g_el^i */
    3144      528374 :     j2 = i2?f2-i2:i2;
    3145      528374 :     u[2+i] = vz_2f[1+j2];
    3146      528374 :     if ((i2+=i+i+1)>=f2) i2-=f2; /* same as i2%=f2 */
    3147             :   }
    3148     1600537 :   for (gi=1,i=i2=0; i<lv; i++)
    3149             :   {
    3150     1600452 :     v[2+i] = Fl_mul(vz_2f[1+i2], vz_el[1+gi], ell);
    3151     1600452 :     gi = Fl_mul(gi, g_el, el);
    3152     1600452 :     if ((i2+=i+i+1)>=f2) i2%=f2; /* i2-=f2 does not work */
    3153             :   }
    3154          85 :   w0 = Flx_mul(u, v, ell) + 1;
    3155          85 :   if (m==1)
    3156             :   { /* el_1=f */
    3157           0 :     for (i=1; i<f; i++) W[i] = Fl_add(w0[i], w0[i+lv], ell);
    3158           0 :     W[f] = w0[f];
    3159             :   }
    3160             :   else
    3161             :   {
    3162          85 :     ulong start = 1+f, end = f+el-1;
    3163          85 :     GEN w = cgetg(end+1, t_VECSMALL);
    3164     2128826 :     for (i=1; i<end; i++) w[i] = w0[i];
    3165          85 :     w[end] = 0;
    3166         360 :     for (i=1; i<m; i++, start+=f)
    3167         550 :       w = both_odd(f,i)? Flv_shift_sub(w, w0, ell, start, end)
    3168         275 :                        : Flv_shift_add(w, w0, ell, start, end);
    3169      528459 :     for (i=0; i<f; i++) W[1+i] = Fl_add(w[1+i], w[1+i+lv], ell);
    3170             :   }
    3171      528374 :   for (i=i2=1; i<f; i++)
    3172             :   {
    3173      528289 :     W[i]=Fl_mul(W[1+i], vz_2f[1+i2], ell);
    3174      528289 :     if ((i2+=i+i+1)>=f2) i2%=f2;
    3175             :   }
    3176             :   /* W[r]=tau_{LL}^{sigma_r}, 1<= r <= f-1 */
    3177          85 :   return gerepilecopy(av, Flv_to_ZV(W));
    3178             : }
    3179             : 
    3180             : static GEN
    3181          90 : gauss_ZX_mul(ulong f, GEN elg, GEN ellg)
    3182             : {
    3183          90 :   pari_sp av = avma, av2;
    3184             :   ulong el, g_el, el_1, f2, lv, lu, m, i, i2, gi;
    3185          90 :   GEN  ell = gel(ellg, 1), g_ell, ell_1, z_2f, z_el, W, vz_2f, vz_el, u, v, w0;
    3186             : 
    3187          90 :   if (lgefint(ell) == 3) return gauss_Flx_mul(f, elg, ellg);
    3188           5 :   g_ell = gel(ellg, 2); ell_1 = subiu(ell, 1);
    3189           5 :   el = elg[1]; g_el = elg[2]; el_1 = el-1;
    3190           5 :   f2 = f<<1; lv=el_1; lu=f; m=el_1/f;
    3191           5 :   z_2f = Fp_pow(g_ell, diviuexact(ell_1, f2), ell);
    3192           5 :   z_el = Fp_pow(g_ell, diviuexact(ell_1, el), ell);
    3193           5 :   W = cgetg(f+1, t_VEC);
    3194           5 :   u = cgetg(lu+2, t_POL); u[1] = evalsigne(1) | evalvarn(0);
    3195           5 :   v = cgetg(lv+2, t_POL); v[1] = evalsigne(1) | evalvarn(0);
    3196           5 :   vz_2f = Fp_powers(z_2f, f2-1, ell);
    3197           5 :   vz_el = Fp_powers(z_el, el_1, ell);
    3198       35264 :   for (gi=1,i=i2=0; i<lu; i++)
    3199             :   {
    3200             :     long j2; /* i2=(i*i)%f2, gi=g_el^i */
    3201       35259 :     j2 = i2?f2-i2:i2;
    3202       35259 :     gel(u, 2+i) = gel(vz_2f, 1+j2);
    3203       35259 :     if ((i2+=i+i+1)>=f2) i2-=f2;
    3204             :   }
    3205       82787 :   for (gi=1,i=i2=0; i<lv; i++)
    3206             :   {
    3207       82782 :     gel(v, 2+i) = Fp_mul(gel(vz_2f, 1+i2), gel(vz_el, 1+gi), ell);
    3208       82782 :     gi = Fl_mul(gi, g_el, el);
    3209       82782 :     if ((i2+=i+i+1)>=f2) i2%=f2;
    3210             :   }
    3211           5 :   w0 = FpX_mul(u, v, ell) + 1; av2 = avma;
    3212           5 :   if (m==1)
    3213             :   {
    3214           0 :     for (i=1; i < f; i++) gel(W,i) = Fp_add(gel(w0, i), gel(w0, i+lv), ell);
    3215           0 :     gel(W, f) = gel(w0, f);
    3216             :   }
    3217             :   else
    3218             :   {
    3219           5 :     ulong start = 1+f, end = f+el-1;
    3220           5 :     GEN w = cgetg(end+1, t_VEC);
    3221      118041 :     for (i=1; i<end; i++) gel(w, i) = gel(w0, i);
    3222           5 :     gel(w, end) = gen_0;
    3223          15 :     for (i=1; i<m; i++, start+=f)
    3224             :     {
    3225          13 :       w = both_odd(f,i)? FpV_shift_sub(w, w0, ell, start, end)
    3226          10 :                        : FpV_shift_add(w, w0, ell, start, end);
    3227          10 :       if ((i & 7) == 0) w = gerepilecopy(av2, w);
    3228             :     }
    3229       35264 :     for (i = 1; i <= f; i++) gel(W, i) = addii(gel(w, i), gel(w, i+lv));
    3230             :   }
    3231       35259 :   for (i = i2 = 1; i < f; i++)
    3232             :   {
    3233       35254 :     gel(W, i) = Fp_mul(gel(W, 1+i), gel(vz_2f, 1+i2), ell);
    3234       35254 :     if ((i2+=i+i+1) >= f2) i2 %= f2;
    3235             :   }
    3236           5 :   return gerepilecopy(av, W);  /* W[r]=tau_{LL}^{sigma_r}, 1<= r <= f-1 */
    3237             : }
    3238             : 
    3239             : /* fast but consumes memory */
    3240             : static GEN
    3241           4 : gauss_el_vell(ulong f, GEN elg, GEN vellg, GEN vz_2f)
    3242             : {
    3243           4 :   pari_sp av = avma, av2;
    3244           4 :   ulong el = elg[1], g_el = elg[2], el_1 = el-1;
    3245           4 :   ulong lv=el_1, f2=f<<1, lu=f, m=el_1/f;
    3246           4 :   GEN W = cgetg(f+1, t_VEC), vz_el, u, v, w0, M;
    3247             :   ulong i, i2, gi;
    3248             : 
    3249           4 :   vz_el = vz_vell(el, vellg, &M);
    3250           4 :   u = cgetg(lu+2, t_POL); u[1] = evalsigne(1) | evalvarn(0);
    3251           4 :   v = cgetg(lv+2, t_POL); v[1] = evalsigne(1) | evalvarn(0);
    3252       25554 :   for (i=i2=0; i<lu; i++)
    3253             :   {
    3254             :     long j2; /* i2=(i*i)%f2, gi=g_el^i */
    3255       25550 :     j2 = i2?f2-i2:i2;
    3256       25550 :     gel(u, 2+i) = gel(vz_2f, 1+j2);
    3257       25550 :     if ((i2+=i+i+1)>=f2) i2%=f2;
    3258             :   }
    3259       86874 :   for (gi=1,i=i2=0; i<lv; i++)
    3260             :   {
    3261       86870 :     gel(v, 2+i) = Fp_mul(gel(vz_2f, 1+i2), gel(vz_el, 1+gi), M);
    3262       86870 :     gi = Fl_mul(gi, g_el, el);
    3263       86870 :     if ((i2+=i+i+1)>=f2) i2%=f2;
    3264             :   }
    3265           4 :   w0 = FpX_mul(u, v, M) + 1; av2 = avma;
    3266           4 :   if (m==1)
    3267             :   { /* el_1=f */
    3268           0 :     for (i=1; i < f; i++) gel(W,i) = Fp_add(gel(w0, i), gel(w0, i+lv), M);
    3269           0 :     gel(W, f) = gel(w0, f);
    3270             :   }
    3271             :   else
    3272             :   {
    3273           4 :     ulong start = 1+f, end = f+el-1;
    3274           4 :     GEN w = cgetg(end+1, t_VEC);
    3275      112420 :     for (i=1; i<end; i++) gel(w, i) = gel(w0, i);
    3276           4 :     gel(w, end) = gen_0;
    3277          19 :     for (i=1; i<m; i++, start+=f)
    3278             :     {
    3279          22 :       w = both_odd(f,i)? FpV_shift_sub(w, w0, M, start, end)
    3280          15 :                        : FpV_shift_add(w, w0, M, start, end);
    3281          15 :       if ((i & 7) == 0) w = gerepilecopy(av2, w);
    3282             :     }
    3283       25554 :     for (i = 1; i <= f; i++) gel(W, i) = Fp_add(gel(w, i), gel(w, i+lv), M);
    3284             :   }
    3285       25550 :   for (i = i2 = 1; i < f; i++)
    3286             :   {
    3287       25546 :     gel(W, i) = Fp_mul(gel(W, 1+i), gel(vz_2f, 1+i2), M);
    3288       25546 :     if ((i2+=i+i+1) >= f2) i2 %= f2;
    3289             :   }
    3290           4 :   return gerepilecopy(av, W);  /* W[r]=tau_{LL}^{sigma_r}, 1<= r <= f-1 */
    3291             : }
    3292             : 
    3293             : static GEN
    3294          94 : norm_chi(GEN K, GEN TAU, ulong p, long d_pow, GEN ell, long j0)
    3295             : {
    3296          94 :   pari_sp av = avma;
    3297          94 :   GEN H = K_get_H(K);
    3298          94 :   ulong d_K, f_K = K_get_f(K), h = K_get_h(K), g_K = K_get_g(K);
    3299          94 :   ulong i, j, gi, pd = upowuu(p, d_pow), d_chi = K_get_dchi(K);
    3300          94 :   GEN z = const_vec(d_chi, gen_1);
    3301          94 :   GEN e_chi = get_e_chi(K, j0, pd, &d_K);
    3302             : 
    3303        1420 :   for (gi=1, i=0; i<d_K; i++)
    3304             :   {
    3305        1326 :     GEN y = gen_1;
    3306      230862 :     for (j=1; j<=h; j++)
    3307      229536 :       y = Fp_mul(y, gel(TAU, Fl_mul(gi, H[j], f_K)), ell);
    3308        1326 :     gi = Fl_mul(gi, g_K, f_K);
    3309        2652 :     for (j=1; j<=d_chi; j++)
    3310             :     {
    3311        1326 :       GEN y2 = Fp_powu(y, e_chi[(i+j-1)%d_K], ell);
    3312        1326 :       gel(z, j) = Fp_mul(gel(z, j), y2, ell);
    3313             :     }
    3314             :   }
    3315          94 :   return gerepilecopy(av, z);
    3316             : }
    3317             : 
    3318             : static void
    3319           2 : imag_MLLn(long *y, GEN K, ulong p, long d_pow, long n,
    3320             :     GEN velg, GEN vellg, long j0)
    3321             : {
    3322           2 :   long f = K_get_f(K), d = upowuu(p, d_pow), row = lg(vellg)-1, i, j, k, nz;
    3323           2 :   GEN g, z, M, vz_2f = vz_vell(f << 1, vellg, &M);
    3324           6 :   for (i=1; i<=n; i++)
    3325             :   {
    3326           4 :     pari_sp av = avma;
    3327           4 :     GEN elg = gel(velg, i);
    3328           4 :     if (DEBUGLEVEL>1) err_printf("(f,el-1)=(%ld,%ld*%ld)\n", f,(elg[1]-1)/f,f);
    3329           4 :     g = gauss_el_vell(f, elg, vellg, vz_2f);
    3330           4 :     z = norm_chi(K, g, p, d_pow, M, j0);
    3331           4 :     nz = lg(z)-1;
    3332           8 :     for (k = 1; k <= nz; k++)
    3333          12 :       for (j = 1; j <= row; j++)
    3334           8 :         y[(j-1)*row+(i-1)*nz+k-1] = get_y(gel(z, k), gel(vellg, j), d);
    3335           4 :     set_avma(av);
    3336             :   }
    3337           2 : }
    3338             : 
    3339             : static void
    3340          42 : imag_MLL1(long *y, GEN K, ulong p, long d_pow, GEN velg, GEN vellg, long j0)
    3341             : {
    3342          42 :   long f = K_get_f(K), d = upowuu(p, d_pow);
    3343          42 :   GEN elg = gel(velg, 1), ellg = gel(vellg, 1), ell = gel(ellg, 1), g, z;
    3344             : 
    3345          42 :   if (DEBUGLEVEL>1) err_printf("(f,el-1)=(%ld,%ld*%ld)\n", f, (elg[1]-1)/f, f);
    3346          42 :   g = gauss_ZX_mul(f, elg, ellg);
    3347          42 :   z = norm_chi(K, g, p, d_pow, ell, j0);
    3348          42 :   y[0] = get_y(gel(z, 1), ellg, d);
    3349          42 : }
    3350             : 
    3351             : static void
    3352          12 : imag_MLL(long *y, GEN K, ulong p, long d_pow, long n, GEN velg, GEN vellg,
    3353             :     long j0)
    3354             : {
    3355          12 :   pari_sp av = avma;
    3356          12 :   long i, j, f = K_get_f(K), d = upowuu(p, d_pow), row = lg(vellg)-1;
    3357             : 
    3358          36 :   for (j=1; j<=row; j++)
    3359             :   {
    3360          24 :     GEN ellg = gel(vellg, j), ell = gel(ellg, 1);
    3361          72 :     for (i=1; i<=n; i++)
    3362             :     {
    3363          48 :       GEN elg = gel(velg, i), g, z;
    3364             :       ulong k, nz;
    3365          48 :       if (DEBUGLEVEL>1) err_printf("(f,el-1)=(%ld,%ld*%ld)\n",f,(elg[1]-1)/f,f);
    3366          48 :       g = gauss_ZX_mul(f, elg, ellg);
    3367          48 :       z = norm_chi(K, g, p, d_pow, ell, j0);
    3368          48 :       nz = lg(z)-1;
    3369          96 :       for (k = 1; k <= nz; k++)
    3370          48 :         y[(j-1)*row+(i-1)*nz+k-1] = get_y(gel(z, k), ellg, d);
    3371          48 :       set_avma(av);
    3372             :     }
    3373             :   }
    3374          12 : }
    3375             : 
    3376             : /* return an upper bound >= 0 if one was found, otherwise return -1.
    3377             :  * set chi-part to be (1) if chi is Teichmuller character.
    3378             :  * B_{1,omega^(-1)} is not p-adic integer. */
    3379             : static GEN
    3380          42 : cyc_imag_MLL(GEN K, ulong p, long d_pow, long j, long flag)
    3381             : {
    3382          42 :   long f = K_get_f(K), d_chi = K_get_dchi(K);
    3383          42 :   long n, n0 = 1, n_el = d_pow, d = upowuu(p, d_pow), rank = n_el*d_chi;
    3384          42 :   GEN df0, velg = const_vec(n_el, NULL), vellg = NULL;
    3385          42 :   GEN oldgr = mkvec2(gen_0, NULL), newgr = mkvec2(gen_0, NULL);
    3386          42 :   long *y0 = (long*)stack_calloc(sizeof(long)*rank*rank);
    3387             : 
    3388          42 :   if (DEBUGLEVEL>1)
    3389           0 :     err_printf("cyc_imag_MLL:p=%ld d_pow=%ld deg(K)=%ld cond(K)=%ld avma=%ld\n",
    3390             :         p, d_pow, K_get_d(K), f, avma);
    3391          42 :   df0 = muluu(d, f%p?f:f/p);
    3392          42 :   gel(velg, 1) = next_el_imag(mkvecsmall2(1, 1), f);
    3393          42 :   if (flag&USE_FULL_EL)
    3394             :   {
    3395           0 :     for (n=2; n<=n_el; n++) gel(velg, n) = next_el_imag(gel(velg, n-1), f);
    3396           0 :     n0 = n_el;
    3397             :   }
    3398          56 :   for (n=n0; n<=n_el; n++) /* loop while structure is unknown */
    3399             :   {
    3400          56 :     pari_sp av2 = avma;
    3401             :     pari_timer ti;
    3402             :     long n_ell, m, M;
    3403             :     GEN y;
    3404          56 :     vellg = set_ell_imag(velg, n, d_chi, df0);
    3405          56 :     n_ell = lg(vellg)-1;  /* equal to n*d_chi */
    3406          56 :     if (DEBUGLEVEL>2) err_printf("velg=%Ps\nvellg=%Ps\n", velg, vellg);
    3407          56 :     if (DEBUGLEVEL>2) timer_start(&ti);
    3408          56 :     if (n_ell==1)
    3409          42 :       imag_MLL1(y0, K, p, d_pow, velg, vellg, j);
    3410          14 :     else if (lgefint(gmael(vellg, n, 1))<=3 || (flag&SAVE_MEMORY))
    3411          12 :       imag_MLL(y0, K, p, d_pow, n, velg, vellg, j);
    3412             :     else
    3413           2 :       imag_MLLn(y0, K, p, d_pow, n, velg, vellg, j);
    3414          56 :     set_avma(av2);
    3415          56 :     if (DEBUGLEVEL>2) timer_printf(&ti, "gauss sum");
    3416          56 :     y = ary2mat(y0, n_ell);
    3417          56 :     if (DEBUGLEVEL>3) err_printf("y=%Ps\n", y);
    3418          56 :     y = ZM_snf(y);
    3419          56 :     if (DEBUGLEVEL>3) err_printf("y=%Ps\n", y);
    3420          56 :     y = make_p_part(y, p, d_pow);
    3421          56 :     if (DEBUGLEVEL>3) err_printf("y=%Ps\n", y);
    3422          56 :     newgr = structure_MLL(y, d_pow);
    3423          56 :     if (DEBUGLEVEL>3)
    3424           0 :       err_printf("d_pow=%ld d_chi=%ld old=%Ps new=%Ps\n",d_pow,d_chi,oldgr,newgr);
    3425          56 :     if (equalsi(d_pow*d_chi, gel(newgr, 1))) break;
    3426          14 :     if ((m = find_del_el(&oldgr, newgr, n, n_el, d_chi)))
    3427           0 :     { M = m = delete_el(velg, m); n--; }
    3428             :     else
    3429          14 :     { M = n+1; m = n; }
    3430          14 :     gel(velg, M) = next_el_imag(gel(velg, m), f);
    3431             :   }
    3432          42 :   return get_str(newgr);
    3433             : }
    3434             : 
    3435             : /* When |A_psi|=p^e, A_psi=(p^e1,...,p^er) (psi=chi^j),
    3436             :  *  return vec[e, [e1, ... ,er], 1].
    3437             :  * If gr str is not determined, return vec[e, [], 1].
    3438             :  * If |A_chi|=1, return vec[0, [], 1].
    3439             :  * If |A_chi|=p, return vec[1, [1], 1].
    3440             :  * If e is not determined, return vec[-1, [], 1].
    3441             :  * If psi is Teichmuller, return vec[0, [], 1].
    3442             :  * B_{1,omega^(-1)} is not p-adic integer. */
    3443             : static GEN
    3444       26334 : cyc_imag(GEN K, GEN B, GEN p, long j, GEN powp, long flag)
    3445             : {
    3446       26334 :   pari_sp av = avma;
    3447       26334 :   GEN MinPol = gel(K, 3), Chi = gel(K, 2), B1, B2, gr;
    3448       26334 :   long x, d_K = K_get_d(K), f_K = K_get_f(K), d_chi = K_get_dchi(K);
    3449             : 
    3450       26334 :   if (f_K == d_K+1 && equaliu(p, f_K) && j == 1) /* Teichmuller */
    3451          77 :     return mkvec3(gen_0, nullvec(), gen_1);
    3452       26257 :   B1 = FpX_rem(ZX_ber_conj(B, j, d_K), MinPol, powp);
    3453       26257 :   B2 = FpX_rem(ZX_ber_den(Chi, j, d_K), MinPol, powp);
    3454       26257 :   if (degpol(B1)<0 || degpol(B2)<0)
    3455             :   {
    3456           0 :     set_avma(av);
    3457           0 :     return mkvec3(gen_m1, nullvec(), gen_1); /* B=0(mod p^pow) */
    3458             :   }
    3459       26257 :   x = ZX_pval(B1, p) - ZX_pval(B2, p);
    3460       26257 :   set_avma(av);
    3461       26257 :   if (x<0) pari_err_BUG("subcyclopclgp [Bernoulli number]");
    3462       26257 :   if (DEBUGLEVEL && x) verbose_output(K, p, x, j);
    3463       26257 :   if (x==0) return mkvec3(gen_0, nullvec(), gen_1); /* trivial */
    3464         588 :   if (x==1) return mkvec3(utoi(d_chi), onevec(d_chi), gen_1);
    3465         140 :   if ((flag&USE_MLL)==0) return mkvec3(utoi(x*d_chi), nullvec(), gen_1);
    3466          42 :   gr = d_K == 2? cyc_buch(-f_K, p, x): cyc_imag_MLL(K, itou(p), x, j, flag);
    3467          42 :   return gerepilecopy(av, mkvec3(utoipos(d_chi * x), gr, gen_1));
    3468             : }
    3469             : 
    3470             : /* handle representatives of all injective characters, d_chi=[Q_p(zeta_d):Q_p],
    3471             :  * d=d_K */
    3472             : static GEN
    3473       10080 : pclgp_cyc_imag(GEN K, GEN p, long start_pow, long max_pow, long flag)
    3474             : {
    3475       10080 :   GEN C = gel(K, 5), Chi = gel(K, 2);
    3476       10080 :   long n_conj = K_get_nconj(K), d_K = K_get_d(K), f_K = K_get_f(K);
    3477       10080 :   long i, pow, n_done = 0;
    3478       10080 :   GEN gr = nullvec(), Done = const_vecsmall(n_conj, 0);
    3479       10080 :   GEN B = zx_ber_num(Chi, f_K, d_K), B_num;
    3480             : 
    3481       10080 :   if (lgefint(p)==3 && n_conj>10) /* mark trivial chi-part by pre-calculation */
    3482             :   {
    3483         595 :     ulong up = itou(p);
    3484         595 :     GEN minpol = ZX_to_Flx(gel(K, 3), up);
    3485        7350 :     for (i=1; i<=n_conj; i++)
    3486             :     {
    3487        7168 :       pari_sp av = avma;
    3488             :       long degB;
    3489        7168 :       B_num = Flx_rem(Flx_ber_conj(B, C[i], d_K, up), minpol, up);
    3490        7168 :       degB = degpol(B_num);
    3491        7168 :       set_avma(av);
    3492        7168 :       if (degB<0) continue;
    3493        6937 :       Done[i] = 1;
    3494        6937 :       if (++n_done == n_conj) return gr;
    3495             :     }
    3496             :   }
    3497        9667 :   for (pow = start_pow; pow<=max_pow; pow++)
    3498             :   {
    3499        9667 :     GEN powp = powiu(p, pow);
    3500       27503 :     for (i = 1; i <= n_conj; i++)
    3501             :     {
    3502             :       GEN z;
    3503       27503 :       if (Done[i]) continue;
    3504       26334 :       z = cyc_imag(K, B, p, C[i], powp, flag);
    3505       26334 :       if (equalim1(gel(z, 1))) continue;
    3506       26334 :       Done[i] = 1;
    3507       26334 :       if (!isintzero(gel(z, 1))) gr = vec_append(gr, z);
    3508       26334 :       if (++n_done == n_conj) return gr;
    3509             :     }
    3510             :   }
    3511           0 :   pari_err_BUG("pclgp_cyc_imag: max_pow is not enough");
    3512             :   return NULL; /*LCOV_EXCL_LINE*/
    3513             : }
    3514             : 
    3515             : static GEN
    3516         392 : gather_part(GEN g, long sgn)
    3517             : {
    3518         392 :   long i, j, l = lg(g), ord = 0, flag = 1;
    3519         392 :   GEN z2 = cgetg(l, t_VEC);
    3520        1778 :   for (i = j = 1; i < l; i++)
    3521             :   {
    3522        1386 :     GEN t = gel(g,i);
    3523        1386 :     if (equaliu(gel(t, 3), sgn))
    3524             :     {
    3525         693 :       ord += itou(gel(t, 1));
    3526         693 :       if (lg(gel(t, 2)) == 1) flag = 0;
    3527         693 :       gel(z2, j++) = gel(t, 2);
    3528             :     }
    3529             :   }
    3530         392 :   if (flag==0 || ord==0) z2 = nullvec();
    3531             :   else
    3532             :   {
    3533         126 :     setlg(z2, j); z2 = shallowconcat1(z2);
    3534         126 :     ZV_sort_inplace(z2); vecreverse_inplace(z2);
    3535             :   }
    3536         392 :   return mkvec2(utoi(ord), z2);
    3537             : }
    3538             : 
    3539             : #ifdef DEBUG
    3540             : static void
    3541             : handling(GEN K)
    3542             : {
    3543             :   long d_K = K_get_d(K), f_K = K_get_f(K), s_K = K_get_s(K), g_K = K_get_g(K);
    3544             :   long d_chi = K_get_dchi(K);
    3545             :   err_printf("  handling %s cyclic subfield K,\
    3546             :       deg(K)=%ld, cond(K)=%ld g_K=%ld d_chi=%ld H=%Ps\n",
    3547             :       s_K? "a real": "an imaginary",d_K,f_K,g_K,d_chi,zv_to_ZV(gmael3(K,1,1,1)));
    3548             : }
    3549             : #endif
    3550             : 
    3551             : /* HH a t_VECSMALL listing group generators
    3552             :  * Aoki and Fukuda, LNCS vol.4076 (2006), 56-74. */
    3553             : static GEN
    3554         161 : pclgp(GEN p0, long f, GEN HH, long degF, long flag)
    3555             : {
    3556             :   long start_pow, max_pow, ip, lp, i, n_f;
    3557         161 :   GEN vH1, z, vData, cycGH, vp = typ(p0) == t_INT? mkvec(p0): p0;
    3558             : 
    3559         161 :   vH1 = GHinit(f, HH, &cycGH); n_f = lg(vH1)-1;
    3560             : #ifdef DEBUG
    3561             :   err_printf("F is %s, deg(F)=%ld, ", srh_1(HH)? "real": "imaginary", degF);
    3562             :   err_printf("cond(F)=%ld, G(F/Q)=%Ps\n",f, cycGH);
    3563             :   err_printf("F has %ld cyclic subfield%s except for Q.\n", n_f,n_f>1?"s":"");
    3564             : #endif
    3565             : 
    3566         161 :   lp = lg(vp); z = cgetg(lp, t_MAT);
    3567         357 :   for (ip = 1; ip < lp; ip++)
    3568             :   {
    3569         196 :     pari_sp av = avma;
    3570         196 :     long n_sub=0, n_chi=0;
    3571         196 :     GEN gr=nullvec(), p = gel(vp, ip), zi;
    3572             :     /* find conductor e of cyclic subfield K and set the subgroup HE of (Z/eZ)^*
    3573             :      * corresponding to K */
    3574         196 :     set_p_f(p, f, &start_pow, &max_pow);
    3575         196 :     vData = const_vec(degF, NULL);
    3576             : 
    3577       16982 :     for (i=1; i<=n_f; i++) /* prescan. set Teichmuller */
    3578             :     {
    3579       16863 :       GEN H1 = gel(vH1, i);
    3580       16863 :       long d_K = _get_d(H1), f_K = _get_f(H1), g_K = _get_g(H1);
    3581             : 
    3582       16863 :       if (f_K == d_K+1 && equaliu(p, f_K)) /* found K=Q(zeta_p) */
    3583             :       {
    3584             :         pari_timer ti;
    3585          77 :         GEN pnmax = powiu(p, max_pow), vNewton, C, MinPol;
    3586          77 :         long d_chi = 1, n_conj = eulerphiu(d_K);
    3587          77 :         ulong pmodd = umodiu(p, d_K);
    3588             : 
    3589          77 :         C = set_C(pmodd, d_K, d_chi, n_conj);
    3590          77 :         MinPol = set_minpol_teich(g_K, p, max_pow);
    3591          77 :         if (DEBUGLEVEL>3) timer_start(&ti);
    3592          77 :         vNewton = FpX_Newton(MinPol, d_K+1, pnmax);
    3593          77 :         if (DEBUGLEVEL>3)
    3594           0 :           timer_printf(&ti, "FpX_Newton: teich: %ld %ld", degpol(MinPol), d_K);
    3595          77 :         gel(vData, d_K) = mkvec4(MinPol, vNewton, C,
    3596             :                                  mkvecsmall2(d_chi, n_conj));
    3597          77 :         break;
    3598             :       }
    3599             :     }
    3600             : 
    3601       20440 :     for (i=1; i<=n_f; i++) /* handle all cyclic K */
    3602             :     {
    3603       20244 :       GEN H1 = gel(vH1, i), K, z1, Chi;
    3604       20244 :       long d_K = _get_d(H1), s_K = _get_s(H1);
    3605             :       pari_sp av2;
    3606             : 
    3607       20244 :       if ((flag&SKIP_PROPER) && degF != d_K) continue;
    3608       20244 :       if (!gel(vData, d_K))
    3609             :       {
    3610             :         pari_timer ti;
    3611         819 :         GEN pnmax = powiu(p, max_pow), vNewton, C, MinPol;
    3612         819 :         ulong pmodd = umodiu(p, d_K);
    3613         819 :         long d_chi = order_f_x(d_K, pmodd), n_conj = eulerphiu(d_K)/d_chi;
    3614             : 
    3615         819 :         C = set_C(pmodd, d_K, d_chi, n_conj);
    3616         819 :         MinPol = set_minpol(d_K, p, max_pow, n_conj);
    3617         819 :         if (DEBUGLEVEL>3) timer_start(&ti);
    3618             :         /* vNewton[2+i] = vNewton[2+i+d_K]. We need vNewton[2+i] for
    3619             :          * 0 <= i < d_K. But vNewton[2+d_K-1] may be 0 and will be deleted.
    3620             :          * So we need vNewton[2+d_K] not to delete vNewton[2+d_K-1]. */
    3621         819 :         vNewton = FpX_Newton(MinPol, d_K+1, pnmax);
    3622         819 :         if (DEBUGLEVEL>3)
    3623           0 :           timer_printf(&ti, "FpX_Newton: %ld %ld", degpol(MinPol), d_K);
    3624         819 :         gel(vData, d_K) = mkvec4(MinPol, vNewton, C,
    3625             :                                  mkvecsmall2(d_chi, n_conj));
    3626             :       }
    3627       20244 :       av2 = avma;
    3628       20244 :       Chi = s_K? NULL: get_chi(H1);
    3629       20244 :       K = shallowconcat(mkvec2(H1, Chi), gel(vData, d_K));
    3630             : #ifdef DEBUG
    3631             :       handling(K);
    3632             : #endif
    3633       20244 :       if (s_K && !(flag&NO_PLUS_PART))
    3634       10164 :         z1 = pclgp_cyc_real(K, p, max_pow, flag);
    3635       10080 :       else if (!s_K && !(flag&NO_MINUS_PART))
    3636       10080 :         z1 = pclgp_cyc_imag(K, p, start_pow, max_pow, flag);
    3637           0 :       else { set_avma(av2); continue; }
    3638       20244 :       n_sub++; n_chi += gmael(vData, d_K, 4)[2]; /* += n_conj */
    3639       20244 :       if (lg(z1) == 1) set_avma(av2);
    3640         658 :       else gr = gerepilecopy(av2, shallowconcat(gr, z1));
    3641             :     }
    3642         196 :     zi = mkcol(p);
    3643         196 :     zi = vec_append(zi, (flag&NO_PLUS_PART)?nullvec():gather_part(gr, 0));
    3644         196 :     zi = vec_append(zi, (flag&NO_MINUS_PART)?nullvec():gather_part(gr, 1));
    3645         196 :     zi = shallowconcat(zi, mkcol3(cycGH, utoi(n_sub), utoi(n_chi)));
    3646         196 :     gel(z, ip) = gerepilecopy(av, zi);
    3647             :   }
    3648         161 :   return typ(p0) == t_INT? shallowtrans(gel(z,1)): shallowtrans(z);
    3649             : }
    3650             : 
    3651             : static GEN
    3652         413 : reduce_gcd(GEN x1, GEN x2)
    3653             : {
    3654         413 :   GEN d = gcdii(x1, x2);
    3655         413 :   x1 = diviiexact(x1, d);
    3656         413 :   x2 = diviiexact(x2, d);
    3657         413 :   return mkvec2(x1, x2);
    3658             : }
    3659             : 
    3660             : /* norm of x0 (= pol of zeta_d with deg <= d-1) by g of order n
    3661             :  * x0^{1+g+g^2+...+g^(n-1)} */
    3662             : static GEN
    3663          49 : ber_norm_cyc(GEN x0, long g, long n, long d)
    3664             : {
    3665          49 :   pari_sp av = avma;
    3666          49 :   long i, ei, di, fi = 0, l = ulogint(n, 2);
    3667          49 :   GEN xi = x0;
    3668          49 :   ei = 1L << l; di = n / ei;
    3669         203 :   for (i = 1; i <= l; i++)
    3670             :   {
    3671         154 :     if (odd(di)) fi += ei;
    3672         154 :     ei = 1L << (l-i); di = n / ei;
    3673         154 :     xi = ZX_mod_Xnm1(ZX_mul(xi, ber_conj(xi, Fl_powu(g, ei, d), d)), d);
    3674         154 :     if (odd(di))
    3675          42 :       xi = ZX_mod_Xnm1(ZX_mul(xi, ber_conj(x0, Fl_powu(g, fi, d), d)), d);
    3676             :   }
    3677          49 :   return gerepilecopy(av, xi);
    3678             : }
    3679             : 
    3680             : /* x0 a ZX of deg < d */
    3681             : static GEN
    3682          21 : ber_norm_by_cyc(GEN x0, long d, GEN MinPol)
    3683             : {
    3684          21 :   pari_sp av=avma;
    3685          21 :   GEN x = x0, z = znstar(utoi(d)), cyc = gel(z, 2), gen = gel(z, 3);
    3686          21 :   long i, l = lg(cyc);
    3687             :   pari_timer ti;
    3688             : 
    3689          21 :   if (DEBUGLEVEL>1) timer_start(&ti);
    3690          70 :   for (i = 1; i < l; i++)
    3691          49 :     x = ber_norm_cyc(x, itou(gmael(gen, i, 2)), itou(gel(cyc, i)), d);
    3692          21 :   if (DEBUGLEVEL>1) timer_printf(&ti, "ber_norm_by_cyc [ber_norm_cyc]");
    3693          21 :   x = ZX_rem(x, MinPol);  /* slow */
    3694          21 :   if (DEBUGLEVEL>1) timer_printf(&ti, "ber_norm_by_cyc [ZX_rem]");
    3695          21 :   if (lg(x) != 3) pari_err_BUG("subcyclohminus [norm of Bernoulli number]");
    3696          21 :   return gerepilecopy(av, gel(x, 2));
    3697             : }
    3698             : 
    3699             : /* MinPol = polcyclo(d_K, 0).
    3700             :  * MinPol = fac*cofac (mod p).
    3701             :  * B is zv.
    3702             :  * K : H1, MinPol, [fac, cofac], C, [d_chi, n_conj] */
    3703             : static long
    3704          98 : ber_norm_by_val(GEN K, GEN B, GEN p)
    3705             : {
    3706          98 :   pari_sp av = avma;
    3707          98 :   GEN MinPol = gel(K, 2), C = gel(K, 4);
    3708          98 :   GEN vfac = gel(K, 3), fac = gel(vfac, 1), cofac = gel(vfac, 2);
    3709          98 :   long d_chi = K_get_dchi(K), n_conj = K_get_nconj(K), d_K = K_get_d(K);
    3710          98 :   long i, r, n_done = 0, x = 0, dcofac = degpol(cofac);
    3711             :   GEN pr, Done;
    3712             : 
    3713          98 :   Done = const_vecsmall(n_conj, 0);
    3714          98 :   if (lgefint(p)==3)
    3715             :   { /* mark trivial chi-part by pre-calculation */
    3716          98 :     ulong up = itou(p);
    3717          98 :     GEN facs = ZX_to_Flx(fac, up);
    3718         196 :     for (i = 1; i <= n_conj; i++)
    3719             :     {
    3720          98 :       pari_sp av2 = avma;
    3721          98 :       GEN B_conj = Flx_rem(Flx_ber_conj(B, C[i], d_K, up), facs, up);
    3722          98 :       long degB = degpol(B_conj);
    3723          98 :       set_avma(av2); if (degB < 0) continue;
    3724           0 :       Done[i] = 1; if (++n_done == n_conj) return gc_long(av, x);
    3725             :     }
    3726             :   }
    3727             :   else
    3728             :   {
    3729           0 :     for (i = 1; i <= n_conj; i++)
    3730             :     {
    3731           0 :       pari_sp av2 = avma;
    3732           0 :       GEN B_conj = FpX_rem(FpX_ber_conj(B, C[i], d_K, p), fac, p);
    3733           0 :       long degB = degpol(B_conj);
    3734           0 :       set_avma(av2); if (degB < 0) continue;
    3735           0 :       Done[i] = 1; if (++n_done == n_conj) return gc_long(av, x);
    3736             :     }
    3737             :   }
    3738         252 :   for (pr = p, r = 2; r; r <<= 1)
    3739             :   {
    3740             :     GEN polr;
    3741         252 :     pr = sqri(pr); /* p^r */
    3742         252 :     polr = (dcofac==0)? FpX_red(MinPol, pr)
    3743         252 :                       : gel(ZpX_liftfact(MinPol, vfac, pr, p, r), 1);
    3744         406 :     for (i = 1; i <= n_conj; i++)
    3745             :     {
    3746         252 :       pari_sp av2 = avma;
    3747             :       GEN B_conj;
    3748             :       long degB;
    3749         252 :       if (Done[i]) continue;
    3750         252 :       B_conj = FpX_rem(FpX_ber_conj(B, C[i], d_K, pr), polr, pr);
    3751         252 :       degB = degpol(B_conj);
    3752         252 :       set_avma(av2); if (degB < 0) continue;
    3753          98 :       x += d_chi * ZX_pval(B_conj, p);
    3754          98 :       Done[i] = 1; if (++n_done == n_conj) return gc_long(av, x);
    3755             :     }
    3756             :   }
    3757             :   pari_err_BUG("ber_norm_by_val"); return 0;/*LCOV_EXCL_LINE*/
    3758             : }
    3759             : 
    3760             : /* n > 2, p = odd prime not dividing n, e > 0, pe = p^e; d = n*p^e
    3761             :  * return generators of the subgroup H of (Z/dZ)^* corresponding to
    3762             :  * Q(zeta_{p^e}): H = {1<=a<=d | gcd(a,n)=1, a=1(mod p^e)} */
    3763             : static GEN
    3764           0 : znstar_subgr(ulong n, ulong pe, ulong d)
    3765             : {
    3766           0 :   GEN z = znstar(utoi(n)), g = gel(z, 3), G;
    3767           0 :   long i, l = lg(g);
    3768           0 :   G = cgetg(l, t_VECSMALL);
    3769           0 :   for (i=1; i<l; i++) G[i] = u_chinese_coprime(itou(gmael(g,i,2)), 1, n, pe, d);
    3770           0 :   return mkvec2(gel(z,2), G);
    3771             : }
    3772             : 
    3773             : /* K is a cyclic extension of degree n*p^e (n>=4 is even).
    3774             :  * x a ZX of deg < n*p^e. */
    3775             : static long
    3776           0 : ber_norm_with_val(GEN x, long n, ulong p, ulong e)
    3777             : {
    3778           0 :   pari_sp av = avma;
    3779           0 :   long i, j, r, degx, pe = upowuu(p, e), d = n*pe;
    3780           0 :   GEN z, gr, gen, y = cgetg(pe+2, t_POL), MinPol = polcyclo(n, 0);
    3781           0 :   y[1] = evalsigne(1) | evalvarn(0);
    3782           0 :   z = znstar_subgr(n, pe, d);
    3783           0 :   gr = gel(z, 1); gen = gel(z, 2); r = lg(gr)-1;
    3784           0 :   for (i=1; i<=r; i++)
    3785           0 :     x = ber_norm_cyc(x, itou(gel(gen, i)), itou(gel(gr, i)), d);
    3786           0 :   degx = degpol(x);
    3787           0 :   for (j=0; j<pe; j++)
    3788             :   {
    3789           0 :     GEN t = pol_zero(n), z;
    3790           0 :     long a = j; /* a=i*pe+j */
    3791           0 :     for (i=0; i<n; i++)
    3792             :     {
    3793           0 :       if (a>degx) break;
    3794           0 :       gel(t, 2+a%n) = gel(x, 2+a);
    3795           0 :       a += pe;
    3796             :     }
    3797           0 :     z = ZX_rem(ZX_renormalize(t, 2+n), MinPol);
    3798           0 :     if (degpol(z)<0) gel(y, 2+j) = gen_0;
    3799           0 :     else if (degpol(z)==0) gel(y, 2+j) = gel(z, 2);
    3800           0 :     else pari_err_BUG("ber_norm_subgr");
    3801             :   }
    3802           0 :   y = ZX_renormalize(y, pe+2);
    3803           0 :   if (e>1) y = ZX_rem(y, polcyclo(pe, 0));
    3804           0 :   return gc_long(av,  ZX_p_val(y, p, e));
    3805             : }
    3806             : 
    3807             : /* K is a cyclic extension of degree 2*p^e. x a ZX of deg < 2*p^e. In most
    3808             :  * cases, deg(x)=2*p^e-1. But deg(x) can be any value in [0, 2*p^e-1]. */
    3809             : static long
    3810         301 : ber_norm_with_val2(GEN x, ulong p, ulong e)
    3811             : {
    3812         301 :   pari_sp av = avma;
    3813         301 :   long i, d = degpol(x), pe = upowuu(p, e);
    3814         301 :   GEN y = pol_zero(pe);
    3815         301 :   if (d == 2*pe-1)
    3816             :   {
    3817       38416 :     for (i = 0; i < pe; i++)
    3818       76230 :       gel(y, 2+i) = odd(i)? subii(gel(x, 2+i+pe), gel(x, 2+i))
    3819       38115 :                           : subii(gel(x, 2+i), gel(x, 2+i+pe));
    3820             :   }
    3821             :   else
    3822             :   {
    3823           0 :     for (i = 0; i < pe && i <= d; i++)
    3824           0 :       gel(y, 2+i) = odd(i)? negi(gel(x, 2+i)): gel(x, 2+i);
    3825           0 :     for (i = pe; i <= d; i++)
    3826           0 :       gel(y, 2+i-pe) = odd(i)? subii(gel(y, 2+i-pe), gel(x, 2+i))
    3827           0 :                              : addii(gel(y, 2+i-pe), gel(x, 2+i));
    3828             :   }
    3829         301 :   y = ZX_renormalize(y, 2+pe);
    3830         301 :   if (e > 1) y = ZX_rem(y, polcyclo(pe, 0));
    3831         301 :   return gc_long(av, ZX_p_val(y, p, e));
    3832             : }
    3833             : 
    3834             : /* K : H1, MinPol, [fac, cofac], C, [d_chi, n_conj] */
    3835             : static GEN
    3836         812 : ber_cyc5(GEN K, GEN p)
    3837             : {
    3838         812 :   pari_sp av = avma;
    3839         812 :   GEN MinPol = gel(K, 2), H = K_get_H(K);
    3840         812 :   long d = K_get_d(K), f = K_get_f(K), h = K_get_h(K), g = K_get_g(K);
    3841         812 :   GEN x, x1, x2, y, B = const_vecsmall(d+1, 0);
    3842         812 :   long i, j, gi, e, f2 = f>>1, dMinPol = degpol(MinPol), chi2 = -1, *B2 = B+2;
    3843             : 
    3844             :   /* get_chi inlined here to save memory */
    3845    18111989 :   for (j=1; j<=h; j++) /* i = 0 */
    3846             :   {
    3847    18111177 :     if (H[j] == 2) chi2 = 0;
    3848    18111177 :     if (H[j] <= f2) B2[0]++; /* Chi[H[j]] = 0 */
    3849             :   }
    3850       97314 :   for (i = 1, gi = g; i < d; i++)
    3851             :   {
    3852    93017085 :     for (j=1; j<=h; j++)
    3853             :     {
    3854    92920583 :       long t = Fl_mul(gi, H[j], f); /* Chi[t] = i */
    3855    92920583 :       if (t == 2) chi2 = i;
    3856    92920583 :       if (t <= f2) B2[i]++;
    3857             :     }
    3858       96502 :     gi = Fl_mul(gi, g, f);
    3859             :   }
    3860         812 :   y = zx_to_ZX(zx_renormalize(B, d+2));
    3861             : 
    3862         812 :   if (p)
    3863             :   {
    3864             :     ulong n;
    3865         399 :     e = u_pvalrem(d, p, &n);
    3866         399 :     if (e == 0)
    3867          98 :       x1 = utoi(ber_norm_by_val(K, B, p));
    3868         301 :     else if (n > 2)
    3869           0 :       x1 = utoi(ber_norm_with_val(y, n, itou(p), e));
    3870             :     else
    3871         301 :       x1 = utoi(ber_norm_with_val2(y, itou(p), e));
    3872             :   }
    3873             :   else
    3874             :   {
    3875         413 :     if (dMinPol > 100)
    3876          21 :       x1 = ber_norm_by_cyc(y, d, MinPol);
    3877             :     else
    3878         392 :       x1 = ZX_resultant(MinPol, ZX_rem(y, MinPol));
    3879             :   }
    3880             : 
    3881         812 :   if (chi2 < 0) /* chi2 = Chi[2] */
    3882           0 :     x2 = shifti(gen_1, 2*dMinPol);
    3883         812 :   else if (chi2 == 0)
    3884          21 :     x2 = shifti(gen_1, dMinPol);
    3885             :   else
    3886             :   {
    3887         791 :     long e = d/ugcd(chi2, d);
    3888         791 :     x2 = powiu(polcyclo_eval(e, gen_2), eulerphiu(d)/eulerphiu(e));
    3889         791 :     x2 = shifti(x2, dMinPol);
    3890             :   }
    3891         812 :   if (p) x = stoi(itou(x1)-Z_pval(x2, p)); else x = reduce_gcd(x1, x2);
    3892         812 :   return gerepilecopy(av, x);
    3893             : }
    3894             : 
    3895             : /*  Hirabayashi-Yoshino, Manuscripta Math. vol.60, 423-436 (1988), Theorem 1
    3896             :  *
    3897             :  *  F is a subfield of Q(zeta_f)
    3898             :  *  f=p^a => Q=1
    3899             :  *  If F=Q(zeta_f), Q=1 <=> f=p^a
    3900             :  *  If f=4*p^a, p^a*q^b (p,q are odd primes), Q=2 <=> [Q(zeta_f):F] is odd */
    3901             : static long
    3902          21 : unit_index(ulong d, ulong f)
    3903             : {
    3904             :   ulong r, d_f;
    3905          21 :   GEN fa = factoru(f), P = gel(fa, 1), E = gel(fa, 2); r = lg(P)-1;
    3906          21 :   if (r==1) return 1;  /* f=P^a */
    3907           7 :   d_f = eulerphiu_fact(fa);
    3908           7 :   if (d==d_f) return 2;  /* F=Q(zeta_f) */
    3909           0 :   if (r==2 && ((P[1]==2 && E[1]==2) || P[1]>2)) return odd(d_f/d)?2:1;
    3910           0 :   return 0;
    3911             : }
    3912             : 
    3913             : /* Compute relative class number h of the subfield K of Q(zeta_f)
    3914             :  * corresponding to the subgroup HH of (Z/fZ)^*.
    3915             :  * If p!=NULL, then return valuation(h,p). */
    3916             : static GEN
    3917         119 : rel_class_num(long f, GEN HH, long degF, GEN p)
    3918             : {
    3919             :   long i, n_f, W, Q;
    3920         119 :   GEN vH1, vData, x, z = gen_1, z1 = gen_0, z2 = mkvec2(gen_1, gen_1);
    3921             : 
    3922         119 :   vH1 = GHinit(f, HH, NULL); n_f = lg(vH1)-1;
    3923         119 :   vData = const_vec(degF, NULL);
    3924        1652 :   for (i=1; i<=n_f; i++)
    3925             :   {
    3926        1533 :     GEN H1 = gel(vH1, i), K;
    3927        1533 :     long d_K = _get_d(H1), s = _get_s(H1);
    3928             : 
    3929        1533 :     if (s) continue;  /* F is real */
    3930             : #ifdef DEBUG
    3931             :     err_printf("  handling %s cyclic subfield K, deg(K)=%ld, cond(K)=%ld\n",
    3932             :         s? "a real": "an imaginary", d_K, _get_f(H1));
    3933             : #endif
    3934         812 :     if (!gel(vData, d_K))
    3935             :     {
    3936             :       GEN C, MinPol, fac, cofac;
    3937             :       ulong d_chi, n_conj;
    3938         497 :       MinPol = polcyclo(d_K,0);
    3939         497 :       if (p && umodui(d_K, p))
    3940          98 :       {
    3941          98 :         ulong pmodd = umodiu(p, d_K);
    3942          98 :         GEN MinPol_p = FpX_red(MinPol, p);
    3943          98 :         d_chi = order_f_x(d_K, pmodd);
    3944          98 :         n_conj = eulerphiu(d_K)/d_chi;
    3945          98 :         if (n_conj==1) fac = MinPol_p;  /* polcyclo(d_K) is irred mod p */
    3946           0 :         else fac = FpX_one_cyclo(d_K, p);
    3947          98 :         cofac = FpX_div(MinPol_p, fac, p);
    3948          98 :         C = set_C(pmodd, d_K, d_chi, n_conj);
    3949             :       }
    3950             :       else
    3951             :       {
    3952         399 :         fac = cofac = C = NULL;
    3953         399 :         d_chi = n_conj = 0;
    3954             :       }
    3955         497 :       gel(vData, d_K) = mkvec5(MinPol, mkvec2(fac, cofac), C,
    3956             :                                NULL, mkvecsmall2(d_chi, n_conj));
    3957             :     }
    3958         812 :     K = vec_prepend(gel(vData, d_K), H1);
    3959         812 :     z = ber_cyc5(K, p);
    3960         812 :     if (p) z1 = addii(z1, z);
    3961             :     else
    3962             :     {
    3963         413 :       gel(z2, 1) = mulii(gel(z2, 1), gel(z, 1));
    3964         413 :       gel(z2, 2) = mulii(gel(z2, 2), gel(z, 2));
    3965             :     }
    3966             :   }
    3967         119 :   W = root_of_1(f, HH);
    3968         119 :   if (p) return addiu(z1, z_pval(W, p));
    3969          21 :   Q = unit_index(degF, f);
    3970          21 :   x = dvmdii(muliu(gel(z2,1), 2 * W), gel(z2,2), &z1);
    3971          21 :   if (signe(z1)) pari_err_BUG("subcyclohminus [norm of Bernoulli number]");
    3972          21 :   if (!Q && mpodd(x)) Q = 2; /* FIXME: can this happen ? */
    3973          21 :   if (Q == 1) x = shifti(x, -1);
    3974          21 :   return mkvec2(x, utoi(Q));
    3975             : }
    3976             : 
    3977             : static void
    3978         336 : checkp(const char *fun, long degF, GEN p)
    3979             : {
    3980         336 :   if (!BPSW_psp(p)) pari_err_PRIME(fun, p);
    3981         329 :   if (equaliu(p, 2)) pari_err_DOMAIN(fun,"p","=", gen_2, p);
    3982         315 :   if (degF && dvdsi(degF, p)) errpdiv(fun, p, degF);
    3983         301 : }
    3984             : 
    3985             : /* if flag is set, handle quadratic fields specially (don't set H) */
    3986             : static long
    3987         427 : subcyclo_init(const char *fun, GEN FH, long *pdegF, GEN *pH, long flag)
    3988             : {
    3989         427 :   long f = 0, degF = 2;
    3990         427 :   GEN F = NULL, H = NULL;
    3991         427 :   if (typ(FH) == t_POL)
    3992             :   {
    3993          70 :     degF = degpol(FH);
    3994          70 :     if (degF < 1 || !RgX_is_ZX(FH)) pari_err_TYPE(fun, FH);
    3995          70 :     if (flag && degF == 2)
    3996             :     {
    3997          49 :       F = coredisc(ZX_disc(FH));
    3998          49 :       if (is_bigint(F))
    3999           0 :         pari_err_IMPL(stack_sprintf("conductor f > %lu in %s", ULONG_MAX, fun));
    4000          49 :       f = itos(F); if (f == 1) degF = 1;
    4001             :     }
    4002             :     else
    4003             :     {
    4004          21 :       GEN z, bnf = Buchall(pol_x(fetch_var()), 0, DEFAULTPREC);
    4005          21 :       z = rnfconductor(bnf, FH); H = gel(z,3);
    4006          21 :       f = subcyclo_nH(fun, gel(z,2), &H);
    4007          21 :       delete_var();
    4008          21 :       H = znstar_generate(f, H); /* group elements */
    4009             :     }
    4010             :   }
    4011             :   else
    4012             :   {
    4013         357 :     long l = lg(FH), fH;
    4014         357 :     if (typ(FH) == t_INT) F = FH;
    4015         273 :     else if (typ(FH) == t_VEC && (l == 2 || l == 3))
    4016             :     {
    4017         273 :       F = gel(FH, 1);
    4018         273 :       if (l == 3) H = gel(FH, 2);
    4019             :     }
    4020           0 :     else pari_err_TYPE(fun, FH);
    4021         357 :     f = subcyclo_nH(fun, F, &H);
    4022         350 :     H = znstar_generate(f, H); /* group elements */
    4023         350 :     fH = znstar_conductor(H);
    4024         350 :     if (fH == 1) degF = 1;
    4025             :     else
    4026             :     {
    4027         350 :       if (fH != f) { H = znstar_reduce_modulus(H, fH); f = fH; }
    4028         350 :       degF = eulerphiu(f) / zv_prod(gel(H, 2));
    4029             :     }
    4030             :   }
    4031         420 :   *pH = H; *pdegF = degF; return f;
    4032             : }
    4033             : 
    4034             : GEN
    4035         210 : subcyclopclgp(GEN FH, GEN p, long flag)
    4036             : {
    4037         210 :   pari_sp av = avma;
    4038             :   GEN H;
    4039         210 :   long degF, f = subcyclo_init("subcyclopclgp", FH, &degF, &H, 0);
    4040         203 :   if (typ(p) == t_VEC)
    4041             :   {
    4042          28 :     long i, l = lg(p);
    4043          77 :     for (i = 1; i < l; i++) checkp("subcyclopclgp", degF, gel(p, i));
    4044          14 :     if (f == 1) { set_avma(av); return const_vec(l-1, nullvec()); }
    4045             :   }
    4046             :   else
    4047             :   {
    4048         175 :     checkp("subcyclopclgp", degF, p);
    4049         154 :     if (f == 1) { set_avma(av); return nullvec(); }
    4050             :   }
    4051         168 :   if (flag >= USE_BASIS) pari_err_FLAG("subcyclopclgp");
    4052         161 :   return gerepilecopy(av, pclgp(p, f, H, degF, flag));
    4053             : }
    4054             : 
    4055             : static GEN
    4056          98 : subcycloiwasawa_i(GEN FH, GEN P, long n)
    4057             : {
    4058             :   long B, p, f, degF;
    4059             :   GEN H;
    4060          98 :   const char *fun = "subcycloiwasawa";
    4061             : 
    4062          98 :   if (typ(P) != t_INT) pari_err_TYPE(fun, P);
    4063          98 :   if (n < 0) pari_err_DOMAIN(fun, "n", "<", gen_0, stoi(n));
    4064          98 :   B = 1L << (BITS_IN_LONG/4);
    4065          98 :   if (is_bigint(P) || cmpiu(P, B) > 0)
    4066           0 :     pari_err_IMPL(stack_sprintf("prime p > %ld in %s", B, fun));
    4067          98 :   p = itos(P);
    4068          98 :   if (p <= 1 || !uisprime(p)) pari_err_PRIME(fun, P);
    4069          98 :   f = subcyclo_init(fun, FH, &degF, &H, 1);
    4070          98 :   if (degF == 1) return NULL;
    4071          98 :   if (degF == 2)
    4072             :   {
    4073          49 :     long m = ((f & 3) == 0)? f / 4: f;
    4074          49 :     if (H && !srh_1(H)) m = -m;
    4075          49 :     if (!n) return quadlambda(p, m);
    4076          28 :     return m < 0? imagquadstkpol(p, m, n): realquadstkpol(p, m, n);
    4077             :   }
    4078          49 :   if (p == 2) pari_err_DOMAIN(fun, "p", "=", gen_2, gen_2);
    4079          49 :   if (srh_1(H)) return NULL;
    4080          49 :   if (degF % p == 0) errpdiv("abeliwasawa", P, degF);
    4081          49 :   return abeliwasawa(p, f, H, degF, n);
    4082             : }
    4083             : GEN
    4084          98 : subcycloiwasawa(GEN FH, GEN P, long n)
    4085             : {
    4086          98 :   pari_sp av = avma;
    4087          98 :   GEN z = subcycloiwasawa_i(FH, P, n);
    4088          98 :   if (!z) { set_avma(av); return n? nullvec(): mkvec(gen_0); }
    4089          98 :   return gerepilecopy(av, z);
    4090             : }
    4091             : 
    4092             : GEN
    4093         119 : subcyclohminus(GEN FH, GEN P)
    4094             : {
    4095         119 :   const char *fun = "subcyclohminus";
    4096         119 :   pari_sp av = avma;
    4097             :   GEN H;
    4098         119 :   long degF, f = subcyclo_init(fun, FH, &degF, &H, 0);
    4099         119 :   if (P)
    4100             :   {
    4101          98 :     if (typ(P) != t_INT) pari_err_TYPE(fun, P);
    4102          98 :     if (isintzero(P)) P = NULL; else checkp(fun, 0, P);
    4103             :   }
    4104         119 :   if (degF == 1 ||  srh_1(H) == 1) return gen_1;
    4105         119 :   return gerepilecopy(av, rel_class_num(f, H, degF, P));
    4106             : }
    4107             : 
    4108             : /* y in H <=> H[y]==1; return n s.t. x^n in H; if H=0, then return 0 */
    4109             : static long
    4110           0 : order_H_x(GEN H, long x, GEN div_H)
    4111             : {
    4112           0 :   long i, y=x, f=H[1], h = lg(div_H);
    4113           0 :   y = Fl_powu(x, div_H[1], f);
    4114           0 :   if (F2v_coeff(H, y)) return div_H[1];
    4115           0 :   for (i=2; i<h; i++)
    4116             :   {
    4117           0 :     y = Fl_mul(y, Fl_powu(x, div_H[i]-div_H[i-1], f), f);
    4118           0 :     if (F2v_coeff(H, y)) return div_H[i];
    4119             :   }
    4120           0 :   pari_err_BUG("znsubgroupgeneratos [order_H_x]");
    4121             :   return 0;/*LCOV_EXCL_LINE*/
    4122             : }
    4123             : 
    4124             : /* find y=xh (h\in H) s.t. y^o=1 */
    4125             : static long
    4126           0 : adjust_H(GEN H, long x, long o, long flag)
    4127             : {
    4128           0 :   ulong y, h, f=H[1];
    4129           0 :   if (flag==0) return x;
    4130           0 :   y = Fl_inv(Fl_powu(x, o, f), f);
    4131           0 :   for (h=1; h<f; h++)
    4132           0 :     if (F2v_coeff(H, h) && y==Fl_powu(h, o, f)) return Fl_mul(x, h, f);
    4133           0 :   pari_err_BUG("znsubgroupgenerators [adjust_H]");
    4134             :   return 1;/*LCOV_EXCL_LINE*/
    4135             : }
    4136             : 
    4137             : static long
    4138           0 : max_order_ele(GEN H1, GEN H2, GEN div_H, long flag)
    4139             : {
    4140           0 :   long i, n=F2v_hamming(H2), max=0, max_i=0, f=H2[1];
    4141           0 :   for (i=1; i<f; i++)
    4142             :   {
    4143             :     long x;
    4144           0 :     if (F2v_coeff(H2, i)==0) continue;
    4145           0 :     x = order_H_x(H1, i, div_H);
    4146           0 :     if (x==n) return adjust_H(H1, i, n, flag);
    4147           0 :     else if (x>max) { max = x; max_i = i; }
    4148             :   }
    4149           0 :   return adjust_H(H1, max_i, max, flag);
    4150             : }
    4151             : 
    4152             : static void
    4153           0 : enlarge_H(GEN H, long x, GEN div_H)
    4154             : {
    4155           0 :   pari_sp av = avma;
    4156           0 :   long i, j, y=x, o, l=lg(H), f=H[1];
    4157             :   GEN H1;
    4158           0 :   if ((o=order_H_x(H, x, div_H))==1) return;
    4159           0 :   H1 = vecsmall_copy(H);
    4160           0 :   for (i=1; i<f; i++)
    4161             :   {
    4162           0 :     if (F2v_coeff(H, i)==0) continue;
    4163           0 :     for (y=x, j=1; j<o; j++)
    4164             :     {
    4165           0 :       F2v_set(H1, Fl_mul(i, y, f));
    4166           0 :       y = Fl_mul(y, x, f);
    4167             :     }
    4168             :   }
    4169           0 :   for (i=2; i<l; i++) H[i] = H1[i];
    4170           0 :   set_avma(av);
    4171             : }
    4172             : 
    4173             : /* H F2v subgroup of (Z/fZ)^*: 1<= a <= f, a in H <=> F2v_coeff(H, a)=1
    4174             :  * Return generators */
    4175             : GEN
    4176           0 : znsubgroupgenerators(GEN H, long flag)
    4177             : {
    4178             :   long a, f, g;
    4179           0 :   pari_sp av = avma;
    4180           0 :   GEN v, H1, H2, z = const_vecsmall(0,0), div_H;
    4181           0 :   if (typ(H)==t_VEC) H = ZV_to_F2v(H);
    4182           0 :   else if (typ(H)==t_VECSMALL) H = Flv_to_F2v(H);
    4183           0 :   else pari_err_TYPE("znsubgroupgenerators", H);
    4184           0 :   f = H[1]; z = cgetg(1,t_VECSMALL);
    4185           0 :   div_H = divisorsu(F2v_hamming(H));
    4186           0 :   H1 = zero_F2v(f); F2v_set(H1, 1);
    4187           0 :   H2 = H;
    4188           0 :   for (a = 2; a <= f; a++)
    4189           0 :     if (F2v_coeff(H, a) && ugcd(a, f) != 1) F2v_set(H, a);
    4190           0 :   while ((g=max_order_ele(H1, H2, div_H, flag)))
    4191             :   {
    4192           0 :     z = vecsmall_append(z, g);
    4193           0 :     enlarge_H(H1, g, div_H);
    4194           0 :     F2v_negimply_inplace(H2, H1);  /* H2 = H2 - H1 */
    4195             :   }
    4196           0 :   v = cgetg(3, t_VEC);
    4197           0 :   gel(v,1) = utoi(f);
    4198           0 :   gel(v,2) = zv_to_ZV(z); return gerepileupto(av, v);
    4199             : }

Generated by: LCOV version 1.13