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