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, °F, &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, °F, &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, °F, &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 : }
|
