Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /*******************************************************************/
16 : /* */
17 : /* RAY CLASS FIELDS */
18 : /* */
19 : /*******************************************************************/
20 : #include "pari.h"
21 : #include "paripriv.h"
22 :
23 : #define DEBUGLEVEL DEBUGLEVEL_bnr
24 :
25 : static GEN
26 1513716 : bnr_get_El(GEN bnr) { return gel(bnr,3); }
27 : static GEN
28 2002457 : bnr_get_U(GEN bnr) { return gel(bnr,4); }
29 : static GEN
30 25396 : bnr_get_Ui(GEN bnr) { return gmael(bnr,4,3); }
31 :
32 : /* faster than Buchray */
33 : GEN
34 35 : bnfnarrow(GEN bnf)
35 : {
36 : GEN nf, cyc, gen, Cyc, Gen, A, GD, v, w, H, invpi, L, R, u, U0, Uoo, archp, sarch;
37 : long r1, j, l, t, RU;
38 : pari_sp av;
39 :
40 35 : bnf = checkbnf(bnf);
41 35 : nf = bnf_get_nf(bnf);
42 35 : r1 = nf_get_r1(nf); if (!r1) return gcopy( bnf_get_clgp(bnf) );
43 :
44 : /* simplified version of nfsign_units; r1 > 0 so bnf.tu = -1 */
45 35 : av = avma; archp = identity_perm(r1);
46 35 : A = bnf_get_logfu(bnf); RU = lg(A)+1;
47 35 : invpi = invr( mppi(nf_get_prec(nf)) );
48 35 : v = cgetg(RU,t_MAT); gel(v, 1) = const_vecsmall(r1, 1); /* nfsign(-1) */
49 98 : for (j=2; j<RU; j++) gel(v,j) = nfsign_from_logarch(gel(A,j-1), invpi, archp);
50 : /* up to here */
51 :
52 35 : v = Flm_image(v, 2); t = lg(v)-1;
53 35 : if (t == r1) { set_avma(av); return gcopy( bnf_get_clgp(bnf) ); }
54 :
55 28 : v = Flm_suppl(v,2); /* v = (sgn(U)|H) in GL_r1(F_2) */
56 28 : H = zm_to_ZM( vecslice(v, t+1, r1) ); /* supplement H of sgn(U) */
57 28 : w = rowslice(Flm_inv(v,2), t+1, r1); /* H*w*z = proj of z on H // sgn(U) */
58 :
59 28 : sarch = nfarchstar(nf, NULL, archp);
60 28 : cyc = bnf_get_cyc(bnf);
61 28 : gen = bnf_get_gen(bnf); l = lg(gen);
62 28 : L = cgetg(l,t_MAT); GD = gmael(bnf,9,3);
63 63 : for (j=1; j<l; j++)
64 : {
65 35 : GEN z = nfsign_from_logarch(gel(GD,j), invpi, archp);
66 35 : gel(L,j) = zc_to_ZC( Flm_Flc_mul(w, z, 2) );
67 : }
68 : /* [cyc, 0; L, 2] = relation matrix for Cl_f */
69 28 : R = shallowconcat(
70 : vconcat(diagonal_shallow(cyc), L),
71 : vconcat(zeromat(l-1, r1-t), scalarmat_shallow(gen_2,r1-t)));
72 28 : Cyc = ZM_snf_group(R, NULL, &u);
73 28 : U0 = rowslice(u, 1, l-1);
74 28 : Uoo = ZM_mul(H, rowslice(u, l, nbrows(u)));
75 28 : l = lg(Cyc); Gen = cgetg(l,t_VEC);
76 91 : for (j = 1; j < l; j++)
77 : {
78 63 : GEN g = gel(U0,j), s = gel(Uoo,j);
79 63 : g = (lg(g) == 1)? gen_1: Q_primpart( idealfactorback(nf,gen,g,0) );
80 63 : if (!ZV_equal0(s))
81 : {
82 28 : GEN a = set_sign_mod_divisor(nf, ZV_to_Flv(s,2), gen_1, sarch);
83 28 : g = is_pm1(g)? a: idealmul(nf, a, g);
84 : }
85 63 : gel(Gen,j) = g;
86 : }
87 28 : return gc_GEN(av, mkvec3(shifti(bnf_get_no(bnf),r1-t), Cyc, Gen));
88 : }
89 :
90 : /********************************************************************/
91 : /** **/
92 : /** REDUCTION MOD IDELE **/
93 : /** **/
94 : /********************************************************************/
95 :
96 : static GEN
97 25515 : compute_fact(GEN nf, GEN U, GEN gen)
98 : {
99 25515 : long i, j, l = lg(U), h = lgcols(U); /* l > 1 */
100 25515 : GEN basecl = cgetg(l,t_VEC), G;
101 :
102 25515 : G = mkvec2(NULL, trivial_fact());
103 54908 : for (j = 1; j < l; j++)
104 : {
105 29393 : GEN z = NULL;
106 98154 : for (i = 1; i < h; i++)
107 : {
108 68761 : GEN g, e = gcoeff(U,i,j); if (!signe(e)) continue;
109 :
110 31801 : g = gel(gen,i);
111 31801 : if (typ(g) != t_MAT)
112 : {
113 20853 : if (z)
114 2247 : gel(z,2) = famat_mulpow_shallow(gel(z,2), g, e);
115 : else
116 18606 : z = mkvec2(NULL, to_famat_shallow(g, e));
117 20853 : continue;
118 : }
119 10948 : gel(G,1) = g;
120 10948 : g = idealpowred(nf,G,e);
121 10948 : z = z? idealmulred(nf,z,g): g;
122 : }
123 29393 : gel(z,2) = famat_reduce(gel(z,2));
124 29393 : gel(basecl,j) = z;
125 : }
126 25515 : return basecl;
127 : }
128 :
129 : static int
130 15449 : too_big(GEN nf, GEN bet)
131 : {
132 15449 : GEN x = nfnorm(nf,bet);
133 15449 : switch (typ(x))
134 : {
135 8974 : case t_INT: return abscmpii(x, gen_1);
136 6475 : case t_FRAC: return abscmpii(gel(x,1), gel(x,2));
137 : }
138 0 : pari_err_BUG("wrong type in too_big");
139 : return 0; /* LCOV_EXCL_LINE */
140 : }
141 :
142 : /* true nf; GTM 193: Algo 4.3.4. Reduce x mod divisor */
143 : static GEN
144 15001 : idealmoddivisor_aux(GEN nf, GEN x, GEN f, GEN sarch)
145 : {
146 15001 : pari_sp av = avma;
147 : GEN a, A;
148 :
149 15001 : if ( is_pm1(gcoeff(f,1,1)) ) /* f = 1 */
150 : {
151 448 : A = idealred(nf, mkvec2(x, gen_1));
152 448 : A = nfinv(nf, gel(A,2));
153 : }
154 : else
155 : {/* given coprime integral ideals x and f (f HNF), compute "small"
156 : * G in x, such that G = 1 mod (f). GTM 193: Algo 4.3.3 */
157 14553 : GEN G = idealaddtoone_raw(nf, x, f);
158 14553 : GEN D = idealaddtoone_i(nf, idealdiv(nf,G,x), f);
159 14553 : A = nfdiv(nf,D,G);
160 : }
161 15001 : if (too_big(nf,A) > 0) return gc_const(av, x);
162 13482 : a = set_sign_mod_divisor(nf, NULL, A, sarch);
163 13482 : if (a != A && too_big(nf,A) > 0) return gc_const(av, x);
164 13482 : return idealmul(nf, a, x);
165 : }
166 :
167 : GEN
168 4214 : idealmoddivisor(GEN bnr, GEN x)
169 : {
170 4214 : GEN nf = bnr_get_nf(bnr), bid = bnr_get_bid(bnr);
171 4214 : return idealmoddivisor_aux(nf, x, bid_get_ideal(bid), bid_get_sarch(bid));
172 : }
173 :
174 : /* v_pr(L0 * cx) */
175 : static long
176 17983 : fast_val(GEN L0, GEN cx, GEN pr)
177 : {
178 17983 : pari_sp av = avma;
179 17983 : long v = typ(L0) == t_INT? 0: ZC_nfval(L0,pr);
180 17983 : if (cx)
181 : {
182 9436 : long w = Q_pval(cx, pr_get_p(pr));
183 9436 : if (w) v += w * pr_get_e(pr);
184 : }
185 17983 : return gc_long(av,v);
186 : }
187 :
188 : /* x coprime to fZ, return y = x mod fZ, y integral */
189 : static GEN
190 4368 : make_integral_Z(GEN x, GEN fZ)
191 : {
192 4368 : GEN d, y = Q_remove_denom(x, &d);
193 4368 : if (d) y = FpC_Fp_mul(y, Fp_inv(d, fZ), fZ);
194 4368 : return y;
195 : }
196 :
197 : /* p pi^(-1) mod f */
198 : static GEN
199 9863 : get_pinvpi(GEN nf, GEN fZ, GEN p, GEN pi, GEN *v)
200 : {
201 9863 : if (!*v) {
202 4368 : GEN invpi = nfinv(nf, pi);
203 4368 : *v = make_integral_Z(RgC_Rg_mul(invpi, p), mulii(p, fZ));
204 : }
205 9863 : return *v;
206 : }
207 : /* uniformizer pi for pr, coprime to F/p */
208 : static GEN
209 10003 : get_pi(GEN F, GEN pr, GEN *v)
210 : {
211 10003 : if (!*v) *v = pr_uniformizer(pr, F);
212 10003 : return *v;
213 : }
214 :
215 : /* true nf */
216 : static GEN
217 32277 : bnr_grp(GEN nf, GEN U, GEN gen, GEN cyc, GEN bid)
218 : {
219 32277 : GEN h = ZV_prod(cyc);
220 : GEN f, fZ, basecl, fa, pr, t, EX, sarch, F, P, vecpi, vecpinvpi;
221 : long i,j,l,lp;
222 :
223 32277 : if (lg(U) == 1) return mkvec3(h, cyc, cgetg(1, t_VEC));
224 25515 : basecl = compute_fact(nf, U, gen); /* generators in factored form */
225 25515 : EX = gel(bid_get_cyc(bid),1); /* exponent of (O/f)^* */
226 25515 : f = bid_get_ideal(bid); fZ = gcoeff(f,1,1);
227 25515 : fa = bid_get_fact(bid);
228 25515 : sarch = bid_get_sarch(bid);
229 25515 : P = gel(fa,1); F = prV_lcm_capZ(P);
230 :
231 25515 : lp = lg(P);
232 25515 : vecpinvpi = cgetg(lp, t_VEC);
233 25515 : vecpi = cgetg(lp, t_VEC);
234 63707 : for (i=1; i<lp; i++)
235 : {
236 38192 : pr = gel(P,i);
237 38192 : gel(vecpi,i) = NULL; /* to be computed if needed */
238 38192 : gel(vecpinvpi,i) = NULL; /* to be computed if needed */
239 : }
240 :
241 25515 : l = lg(basecl);
242 54908 : for (i=1; i<l; i++)
243 : {
244 : GEN p, pi, pinvpi, dmulI, mulI, G, I, A, e, L, newL;
245 : long la, v, k;
246 : pari_sp av;
247 : /* G = [I, A=famat(L,e)] is a generator, I integral */
248 29393 : G = gel(basecl,i);
249 29393 : I = gel(G,1);
250 29393 : A = gel(G,2); L = gel(A,1); e = gel(A,2);
251 : /* if no reduction took place in compute_fact, everybody is still coprime
252 : * to f + no denominators */
253 29393 : if (!I) { gel(basecl,i) = famat_to_nf_moddivisor(nf, L, e, bid); continue; }
254 10787 : if (lg(A) == 1) { gel(basecl,i) = I; continue; }
255 :
256 : /* compute mulI so that mulI * I coprime to f
257 : * FIXME: use idealcoprime ??? (Less efficient. Fix idealcoprime!) */
258 10787 : dmulI = mulI = NULL;
259 25963 : for (j=1; j<lp; j++)
260 : {
261 15176 : pr = gel(P,j);
262 15176 : v = idealval(nf, I, pr);
263 15176 : if (!v) continue;
264 3801 : p = pr_get_p(pr);
265 3801 : pi = get_pi(F, pr, &gel(vecpi,j));
266 3801 : pinvpi = get_pinvpi(nf, fZ, p, pi, &gel(vecpinvpi,j));
267 3801 : t = nfpow_u(nf, pinvpi, (ulong)v);
268 3801 : mulI = mulI? nfmuli(nf, mulI, t): t;
269 3801 : t = powiu(p, v);
270 3801 : dmulI = dmulI? mulii(dmulI, t): t;
271 : }
272 :
273 : /* make all components of L coprime to f.
274 : * Assuming (L^e * I, f) = 1, then newL^e * mulI = L^e */
275 10787 : la = lg(e); newL = cgetg(la, t_VEC);
276 21707 : for (k=1; k<la; k++)
277 : {
278 10920 : GEN cx, LL = nf_to_scalar_or_basis(nf, gel(L,k));
279 10920 : GEN L0 = Q_primitive_part(LL, &cx); /* LL = L0*cx (faster nfval) */
280 28903 : for (j=1; j<lp; j++)
281 : {
282 17983 : pr = gel(P,j);
283 17983 : v = fast_val(L0,cx, pr); /* = val_pr(LL) */
284 17983 : if (!v) continue;
285 6202 : p = pr_get_p(pr);
286 6202 : pi = get_pi(F, pr, &gel(vecpi,j));
287 6202 : if (v > 0)
288 : {
289 6062 : pinvpi = get_pinvpi(nf, fZ, p, pi, &gel(vecpinvpi,j));
290 6062 : t = nfpow_u(nf,pinvpi, (ulong)v);
291 6062 : LL = nfmul(nf, LL, t);
292 6062 : LL = gdiv(LL, powiu(p, v));
293 : }
294 : else
295 : {
296 140 : t = nfpow_u(nf,pi,(ulong)(-v));
297 140 : LL = nfmul(nf, LL, t);
298 : }
299 : }
300 10920 : LL = make_integral(nf,LL,f,P);
301 10920 : gel(newL,k) = typ(LL) == t_INT? LL: FpC_red(LL, fZ);
302 : }
303 :
304 10787 : av = avma;
305 : /* G in nf, = L^e mod f */
306 10787 : G = famat_to_nf_modideal_coprime(nf, newL, e, f, EX);
307 10787 : if (mulI)
308 : {
309 3787 : G = nfmuli(nf, G, mulI);
310 3787 : G = typ(G) == t_COL? ZC_hnfrem(G, ZM_Z_mul(f, dmulI))
311 3787 : : modii(G, mulii(fZ,dmulI));
312 3787 : G = RgC_Rg_div(G, dmulI);
313 : }
314 10787 : G = set_sign_mod_divisor(nf,A,G,sarch);
315 10787 : I = idealmul(nf,I,G);
316 : /* more or less useless, but cheap at this point */
317 10787 : I = idealmoddivisor_aux(nf,I,f,sarch);
318 10787 : gel(basecl,i) = gc_GEN(av, I);
319 : }
320 25515 : return mkvec3(h, cyc, basecl);
321 : }
322 :
323 : /********************************************************************/
324 : /** **/
325 : /** INIT RAY CLASS GROUP **/
326 : /** **/
327 : /********************************************************************/
328 : GEN
329 276095 : bnr_subgroup_check(GEN bnr, GEN H, GEN *pdeg)
330 : {
331 276095 : GEN no = bnr_get_no(bnr);
332 276095 : if (H && isintzero(H)) H = NULL;
333 276093 : if (H)
334 : {
335 117824 : GEN h, cyc = bnr_get_cyc(bnr);
336 117824 : switch(typ(H))
337 : {
338 2548 : case t_INT:
339 2548 : H = scalarmat_shallow(H, lg(cyc)-1);
340 : /* fall through */
341 38710 : case t_MAT:
342 38710 : RgM_check_ZM(H, "bnr_subgroup_check");
343 38710 : H = ZM_hnfmodid(H, cyc);
344 38710 : break;
345 79114 : case t_VEC:
346 79114 : if (char_check(cyc, H)) { H = charker(cyc, H); break; }
347 0 : default: pari_err_TYPE("bnr_subgroup_check", H);
348 : }
349 117824 : h = ZM_det_triangular(H);
350 117823 : if (equalii(h, no)) H = NULL; else no = h;
351 : }
352 276092 : if (pdeg) *pdeg = no;
353 276092 : return H;
354 : }
355 :
356 : /* c a rational content (NULL or t_INT or t_FRAC), return u*c as a ZM/d */
357 : static GEN
358 239868 : ZM_content_mul(GEN u, GEN c, GEN *pd)
359 : {
360 239868 : *pd = gen_1;
361 239868 : if (c)
362 : {
363 164745 : if (typ(c) == t_FRAC) { *pd = gel(c,2); c = gel(c,1); }
364 164745 : if (!is_pm1(c)) u = ZM_Z_mul(u, c);
365 : }
366 239865 : return u;
367 : }
368 :
369 : /* bnr natural generators: bnf gens made coprime to modulus + bid gens.
370 : * Beware: if bnr includes MOD, we may have #El < #bnf.ge*/
371 : static GEN
372 48090 : get_Gen(GEN bnf, GEN bid, GEN El)
373 : {
374 48090 : GEN nf = bnf_get_nf(bnf), gen = bnf_get_gen(bnf), Gen;
375 48090 : long i, l = lg(El);
376 48090 : if (lg(gen) > l) gen = vec_shorten(gen, l-1);
377 48090 : Gen = shallowconcat(gen, bid_get_gen(bid));
378 67060 : for (i = 1; i < l; i++)
379 : {
380 18970 : GEN e = gel(El,i);
381 18970 : if (!isint1(e)) gel(Gen,i) = idealmul(nf, gel(El,i), gel(Gen,i));
382 : }
383 48090 : return Gen;
384 : }
385 :
386 : static GEN
387 269811 : Buchraymod_i(GEN bnf, GEN module, long flag, GEN MOD)
388 : {
389 : GEN nf, cyc0, cyc, gen, Cyc, clg, h, logU, U, Ui, vu;
390 : GEN bid, cycbid, H, El;
391 : long RU, Ri, j, ngen;
392 269811 : const long add_gen = flag & nf_GEN;
393 269811 : const long do_init = flag & nf_INIT;
394 :
395 269811 : if (MOD && typ(MOD) != t_INT)
396 0 : pari_err_TYPE("bnrinit [incorrect cycmod]", MOD);
397 269811 : bnf = checkbnf(bnf);
398 269803 : nf = bnf_get_nf(bnf);
399 269801 : RU = lg(nf_get_roots(nf))-1; /* #K.futu */
400 269801 : El = NULL; /* gcc -Wall */
401 269801 : cyc = cyc0 = bnf_get_cyc(bnf);
402 269799 : gen = bnf_get_gen(bnf); ngen = lg(cyc)-1;
403 :
404 269793 : bid = checkbid_i(module);
405 269796 : if (!bid) bid = Idealstarmod(nf,module,nf_GEN|nf_INIT, MOD);
406 269820 : cycbid = bid_get_cyc(bid);
407 269816 : if (MOD) cyc = ZV_snfclean(ZV_snf_gcd(cyc, MOD));
408 269810 : Ri = lg(cycbid)-1;
409 269810 : if (Ri || add_gen || do_init)
410 : {
411 269811 : GEN fx = bid_get_fact(bid);
412 269808 : long n = Ri? ngen: lg(cyc)-1;
413 269808 : El = cgetg(n+1, t_VEC);
414 306917 : for (j = 1; j <= n; j++)
415 : {
416 37107 : GEN c = idealcoprimefact(nf, gel(gen,j), fx);
417 37106 : gel(El,j) = nf_to_scalar_or_basis(nf,c);
418 : }
419 : }
420 269809 : if (!Ri)
421 : {
422 29939 : GEN no, Gen = add_gen? get_Gen(bnf, bid, El): NULL;
423 29939 : if (MOD) { ngen = lg(cyc)-1; no = ZV_prod(cyc); } else no = bnf_get_no(bnf);
424 29939 : clg = add_gen? mkvec3(no, cyc, Gen): mkvec2(no, cyc);
425 29939 : if (!do_init) return clg;
426 29939 : U = matid(ngen);
427 29939 : U = mkvec3(U, cgetg(1,t_MAT), U);
428 29939 : vu = mkvec3(cgetg(1,t_MAT), matid(RU), gen_1);
429 29939 : return mkvecn(6, bnf, bid, El, U, clg, vu);
430 : }
431 :
432 239870 : logU = ideallog_units0(bnf, bid, MOD);
433 239845 : if (do_init)
434 : { /* (log(Units)|D) * u = (0 | H) */
435 239845 : GEN c1,c2, u,u1,u2, Hi, D = shallowconcat(logU, diagonal_shallow(cycbid));
436 239873 : H = ZM_hnfall_i(D, &u, 1);
437 239859 : u1 = matslice(u, 1,RU, 1,RU);
438 239875 : u2 = matslice(u, 1,RU, RU+1,lg(u)-1);
439 : /* log(Units) (u1|u2) = (0|H) (mod D), H HNF */
440 :
441 239875 : u1 = ZM_lll(u1, 0.99, LLL_INPLACE);
442 239865 : Hi = Q_primitive_part(RgM_inv_upper(H), &c1);
443 239851 : u2 = ZM_mul(ZM_reducemodmatrix(u2,u1), Hi);
444 239861 : u2 = Q_primitive_part(u2, &c2);
445 239861 : u2 = ZM_content_mul(u2, mul_content(c1,c2), &c2);
446 239865 : vu = mkvec3(u2,u1,c2); /* u2/c2 = H^(-1) (mod Im u1) */
447 : }
448 : else
449 : {
450 0 : H = ZM_hnfmodid(logU, cycbid);
451 0 : vu = NULL; /* -Wall */
452 : }
453 239868 : if (!ngen)
454 214052 : h = H;
455 : else
456 : {
457 25816 : GEN L = cgetg(ngen+1, t_MAT), cycgen = bnf_build_cycgen(bnf);
458 52416 : for (j=1; j<=ngen; j++)
459 : {
460 26600 : GEN c = gel(cycgen,j), e = gel(El,j);
461 26600 : if (!equali1(e)) c = famat_mulpow_shallow(c, e, gel(cyc0,j));
462 26600 : gel(L,j) = ideallogmod(nf, c, bid, MOD); /* = log(Gen[j]^cyc[j]) */
463 : }
464 : /* [cyc0, 0; -L, H] = relation matrix for generators Gen of Cl_f */
465 25816 : h = shallowconcat(vconcat(diagonal_shallow(cyc0), ZM_neg(L)),
466 : vconcat(zeromat(ngen, Ri), H));
467 25816 : h = MOD? ZM_hnfmodid(h, MOD): ZM_hnf(h);
468 : }
469 239868 : Cyc = ZM_snf_group(h, &U, &Ui);
470 : /* Gen = clg.gen*U, clg.gen = Gen*Ui */
471 32277 : clg = add_gen? bnr_grp(nf, Ui, get_Gen(bnf, bid, El), Cyc, bid)
472 239875 : : mkvec2(ZV_prod(Cyc), Cyc);
473 239878 : if (!do_init) return clg;
474 239878 : U = mkvec3(vecslice(U, 1,ngen), vecslice(U,ngen+1,lg(U)-1), Ui);
475 239877 : return mkvecn(6, bnf, bid, El, U, clg, vu);
476 : }
477 : GEN
478 41391 : Buchray(GEN bnf, GEN f, long flag)
479 41391 : { return Buchraymod(bnf, f, flag, NULL); }
480 : GEN
481 253805 : Buchraymod(GEN bnf, GEN f, long flag, GEN MOD)
482 : {
483 253805 : pari_sp av = avma;
484 253805 : return gc_GEN(av, Buchraymod_i(bnf, f, flag, MOD));
485 : }
486 : GEN
487 211380 : bnrinitmod(GEN bnf, GEN f, long flag, GEN MOD)
488 : {
489 211380 : switch(flag)
490 : {
491 211334 : case 0: flag = nf_INIT; break;
492 49 : case 1: flag = nf_INIT | nf_GEN; break;
493 0 : default: pari_err_FLAG("bnrinit");
494 : }
495 211383 : return Buchraymod(bnf, f, flag, MOD);
496 : }
497 : GEN
498 0 : bnrinit0(GEN bnf, GEN ideal, long flag)
499 0 : { return bnrinitmod(bnf, ideal, flag, NULL); }
500 :
501 : GEN
502 112 : bnrclassno(GEN bnf,GEN ideal)
503 : {
504 : GEN h, D, bid, cycbid;
505 112 : pari_sp av = avma;
506 :
507 112 : bnf = checkbnf(bnf);
508 112 : h = bnf_get_no(bnf);
509 112 : bid = checkbid_i(ideal);
510 112 : if (!bid) bid = Idealstar(bnf_get_nf(bnf), ideal, nf_INIT);
511 105 : cycbid = bid_get_cyc(bid);
512 105 : if (lg(cycbid) == 1) { set_avma(av); return icopy(h); }
513 84 : D = ideallog_units(bnf, bid); /* (Z_K/f)^* / units ~ Z^n / D */
514 84 : D = ZM_hnfmodid(D,cycbid);
515 84 : return gc_INT(av, mulii(h, ZM_det_triangular(D)));
516 : }
517 : GEN
518 105 : bnrclassno0(GEN A, GEN B, GEN C)
519 : {
520 105 : pari_sp av = avma;
521 105 : GEN h, H = NULL;
522 : /* adapted from ABC_to_bnr, avoid costly bnrinit if possible */
523 105 : if (typ(A) == t_VEC)
524 105 : switch(lg(A))
525 : {
526 14 : case 7: /* bnr */
527 14 : checkbnr(A); H = B;
528 14 : break;
529 91 : case 11: /* bnf */
530 91 : if (!B) pari_err_TYPE("bnrclassno [bnf+missing conductor]",A);
531 91 : if (!C) return bnrclassno(A, B);
532 7 : A = Buchray(A, B, nf_INIT); H = C;
533 7 : break;
534 0 : default: checkbnf(A);/*error*/
535 : }
536 0 : else checkbnf(A);/*error*/
537 :
538 21 : H = bnr_subgroup_check(A, H, &h);
539 21 : if (!H) { set_avma(av); return icopy(h); }
540 14 : return gc_INT(av, h);
541 : }
542 :
543 : /* ZMV_ZCV_mul for two matrices U = [Ux,Uy], it may have more components
544 : * (ignored) and vectors x,y */
545 : static GEN
546 1390930 : ZM2_ZC2_mul(GEN U, GEN x, GEN y)
547 : {
548 1390930 : GEN Ux = gel(U,1), Uy = gel(U,2);
549 1390930 : if (lg(Ux) == 1) return ZM_ZC_mul(Uy,y);
550 163174 : if (lg(Uy) == 1) return ZM_ZC_mul(Ux,x);
551 163174 : return ZC_add(ZM_ZC_mul(Ux,x), ZM_ZC_mul(Uy,y));
552 : }
553 :
554 : GEN
555 1508899 : bnrisprincipalmod(GEN bnr, GEN x, GEN MOD, long flag)
556 : {
557 1508899 : pari_sp av = avma;
558 : GEN E, G, clgp, bnf, nf, bid, ex, cycray, alpha, El;
559 : int trivialbid;
560 :
561 1508899 : checkbnr(bnr);
562 1508900 : El = bnr_get_El(bnr);
563 1508900 : cycray = bnr_get_cyc(bnr);
564 1508900 : if (MOD && flag) pari_err_FLAG("bnrisprincipalmod [MOD!=NULL and flag!=0]");
565 1508900 : if (lg(cycray) == 1 && !(flag & nf_GEN)) return cgetg(1,t_COL);
566 1508739 : if (MOD) cycray = ZV_snf_gcd(cycray, MOD);
567 :
568 1508739 : bnf = bnr_get_bnf(bnr); nf = bnf_get_nf(bnf);
569 1508739 : bid = bnr_get_bid(bnr);
570 1508739 : trivialbid = lg(bid_get_cyc(bid)) == 1;
571 1508739 : if (trivialbid)
572 : {
573 117809 : ex = isprincipal(bnf, x);
574 117809 : setlg(ex, lg(cycray)); /* can happen with MOD */
575 : }
576 : else
577 : {
578 1390930 : GEN v = bnfisprincipal0(bnf, x, nf_FORCE|nf_GENMAT);
579 1390930 : GEN e = gel(v,1), b = gel(v,2);
580 1390930 : long i, j = lg(e);
581 1558832 : for (i = 1; i < j; i++) /* modify b as if bnf.gen were El*bnf.gen */
582 167902 : if (typ(gel(El,i)) != t_INT && signe(gel(e,i))) /* <==> != 1 */
583 31308 : b = famat_mulpow_shallow(b, gel(El,i), negi(gel(e,i)));
584 1390930 : if (!MOD && !(flag & nf_GEN)) MOD = gel(cycray,1);
585 1390930 : ex = ZM2_ZC2_mul(bnr_get_U(bnr), e, ideallogmod(nf, b, bid, MOD));
586 : }
587 1508738 : ex = ZV_ZV_mod(ex, cycray);
588 1508735 : if (!(flag & (nf_GEN|nf_GENMAT))) return gc_upto(av, ex);
589 :
590 : /* compute generator */
591 7049 : E = ZC_neg(ex);
592 7049 : clgp = bnr_get_clgp(bnr);
593 7049 : if (lg(clgp) == 4)
594 21 : G = abgrp_get_gen(clgp);
595 : else
596 : {
597 7028 : G = get_Gen(bnf, bid, El);
598 7028 : E = ZM_ZC_mul(bnr_get_Ui(bnr), E);
599 : }
600 7049 : alpha = isprincipalfact(bnf, x, G, E, nf_GENMAT|nf_GEN_IF_PRINCIPAL|nf_FORCE);
601 7049 : if (alpha == gen_0) pari_err_BUG("isprincipalray");
602 7049 : if (!trivialbid)
603 : {
604 7049 : GEN v = gel(bnr,6), u2 = gel(v,1), u1 = gel(v,2), du2 = gel(v,3);
605 7049 : GEN y = ZM_ZC_mul(u2, ideallog(nf, alpha, bid));
606 7049 : if (!is_pm1(du2)) y = ZC_Z_divexact(y,du2);
607 7049 : y = ZC_reducemodmatrix(y, u1);
608 7049 : if (!ZV_equal0(y))
609 : {
610 4998 : GEN U = shallowcopy(bnf_build_units(bnf));
611 4998 : settyp(U, t_COL);
612 4998 : alpha = famat_div_shallow(alpha, mkmat2(U,y));
613 : }
614 : }
615 7049 : alpha = famat_reduce(alpha);
616 7049 : if (!(flag & nf_GENMAT)) alpha = nffactorback(nf, alpha, NULL);
617 7049 : return gc_GEN(av, mkvec2(ex,alpha));
618 : }
619 :
620 : GEN
621 415000 : bnrisprincipal(GEN bnr, GEN x, long flag)
622 415000 : { return bnrisprincipalmod(bnr, x, NULL, flag); }
623 :
624 : GEN
625 407916 : isprincipalray(GEN bnr, GEN x) { return bnrisprincipal(bnr,x,0); }
626 : GEN
627 0 : isprincipalraygen(GEN bnr, GEN x) { return bnrisprincipal(bnr,x,nf_GEN); }
628 :
629 : /* N! / N^N * (4/pi)^r2 * sqrt(|D|) */
630 : GEN
631 0 : minkowski_bound(GEN D, long N, long r2, long prec)
632 : {
633 0 : pari_sp av = avma;
634 0 : GEN c = divri(mpfactr(N,prec), powuu(N,N));
635 0 : if (r2) c = mulrr(c, powru(divur(4,mppi(prec)), r2));
636 0 : c = mulrr(c, gsqrt(absi_shallow(D),prec));
637 0 : return gc_leaf(av, c);
638 : }
639 :
640 : /* N = [K:Q] > 1, D = disc(K) */
641 : static GEN
642 63 : zimmertbound(GEN D, long N, long R2)
643 : {
644 63 : pari_sp av = avma;
645 : GEN w;
646 :
647 63 : if (N > 20) w = minkowski_bound(D, N, R2, DEFAULTPREC);
648 : else
649 : {
650 63 : const double c[19][11] = {
651 : {/*2*/ 0.6931, 0.45158},
652 : {/*3*/ 1.71733859, 1.37420604},
653 : {/*4*/ 2.91799837, 2.50091538, 2.11943331},
654 : {/*5*/ 4.22701425, 3.75471588, 3.31196660},
655 : {/*6*/ 5.61209925, 5.09730381, 4.60693851, 4.14303665},
656 : {/*7*/ 7.05406203, 6.50550021, 5.97735406, 5.47145968},
657 : {/*8*/ 8.54052636, 7.96438858, 7.40555445, 6.86558259, 6.34608077},
658 : {/*9*/ 10.0630022, 9.46382812, 8.87952524, 8.31139202, 7.76081149},
659 : {/*10*/11.6153797, 10.9966020, 10.3907654, 9.79895170, 9.22232770, 8.66213267},
660 : {/*11*/13.1930961, 12.5573772, 11.9330458, 11.3210061, 10.7222412, 10.1378082},
661 : {/*12*/14.7926394, 14.1420915, 13.5016616, 12.8721114, 12.2542699, 11.6490374,
662 : 11.0573775},
663 : {/*13*/16.4112395, 15.7475710, 15.0929680, 14.4480777, 13.8136054, 13.1903162,
664 : 12.5790381},
665 : {/*14*/18.0466672, 17.3712806, 16.7040780, 16.0456127, 15.3964878, 14.7573587,
666 : 14.1289364, 13.5119848},
667 : {/*15*/19.6970961, 19.0111606, 18.3326615, 17.6620757, 16.9999233, 16.3467686,
668 : 15.7032228, 15.0699480},
669 : {/*16*/21.3610081, 20.6655103, 19.9768082, 19.2953176, 18.6214885, 17.9558093,
670 : 17.2988108, 16.6510652, 16.0131906},
671 :
672 : {/*17*/23.0371259, 22.3329066, 21.6349299, 20.9435607, 20.2591899, 19.5822454,
673 : 18.9131878, 18.2525157, 17.6007672},
674 :
675 : {/*18*/24.7243611, 24.0121449, 23.3056902, 22.6053167, 21.9113705, 21.2242247,
676 : 20.5442836, 19.8719830, 19.2077941, 18.5522234},
677 :
678 : {/*19*/26.4217792, 25.7021950, 24.9879497, 24.2793271, 23.5766321, 22.8801952,
679 : 22.1903709, 21.5075437, 20.8321263, 20.1645647},
680 : {/*20*/28.1285704, 27.4021674, 26.6807314, 25.9645140, 25.2537867, 24.5488420,
681 : 23.8499943, 23.1575823, 22.4719720, 21.7935548, 21.1227537}
682 : };
683 63 : w = mulrr(dbltor(exp(-c[N-2][R2])), gsqrt(absi_shallow(D),DEFAULTPREC));
684 : }
685 63 : return gc_INT(av, ceil_safe(w));
686 : }
687 :
688 : /* return \gamma_n^n if known, an upper bound otherwise */
689 : GEN
690 63 : Hermite_bound(long n, long prec)
691 : {
692 : GEN h,h1;
693 : pari_sp av;
694 :
695 63 : switch(n)
696 : {
697 35 : case 1: return gen_1;
698 14 : case 2: retmkfrac(utoipos(4), utoipos(3));
699 7 : case 3: return gen_2;
700 7 : case 4: return utoipos(4);
701 0 : case 5: return utoipos(8);
702 0 : case 6: retmkfrac(utoipos(64), utoipos(3));
703 0 : case 7: return utoipos(64);
704 0 : case 8: return utoipos(256);
705 0 : case 24: return int2n(48);
706 : }
707 0 : av = avma;
708 0 : h = powru(divur(2,mppi(prec)), n);
709 0 : h1 = sqrr(ggamma(uutoQ(n+4,2),prec));
710 0 : return gc_leaf(av, mulrr(h,h1));
711 : }
712 :
713 : /* 1 if L (= nf != Q) primitive for sure, 0 if MAYBE imprimitive (may have a
714 : * subfield K) */
715 : static long
716 35 : isprimitive(GEN nf)
717 : {
718 35 : long p, i, l, ep, N = nf_get_degree(nf);
719 : GEN D, fa;
720 :
721 35 : p = ucoeff(factoru(N), 1,1); /* smallest prime | N */
722 35 : if (p == N) return 1; /* prime degree */
723 :
724 : /* N = [L:Q] = product of primes >= p, same is true for [L:K]
725 : * d_L = t d_K^[L:K] --> check that some q^p divides d_L */
726 0 : D = nf_get_disc(nf);
727 0 : fa = gel(absZ_factor_limit(D,0),2); /* list of v_q(d_L). Don't check large primes */
728 0 : if (mod2(D)) i = 1;
729 : else
730 : { /* q = 2 */
731 0 : ep = itos(gel(fa,1));
732 0 : if ((ep>>1) >= p) return 0; /* 2 | d_K ==> 4 | d_K */
733 0 : i = 2;
734 : }
735 0 : l = lg(fa);
736 0 : for ( ; i < l; i++)
737 : {
738 0 : ep = itos(gel(fa,i));
739 0 : if (ep >= p) return 0;
740 : }
741 0 : return 1;
742 : }
743 :
744 : static GEN
745 0 : dft_bound(void)
746 : {
747 0 : if (DEBUGLEVEL>1) err_printf("Default bound for regulator: 0.2\n");
748 0 : return dbltor(0.2);
749 : }
750 :
751 : static GEN
752 35 : regulatorbound(GEN bnf)
753 : {
754 : long N, R1, R2, R;
755 : GEN nf, dK, p1, c1;
756 :
757 35 : nf = bnf_get_nf(bnf); N = nf_get_degree(nf);
758 35 : if (!isprimitive(nf)) return dft_bound();
759 :
760 35 : dK = absi_shallow(nf_get_disc(nf));
761 35 : nf_get_sign(nf, &R1, &R2); R = R1+R2-1;
762 35 : c1 = (!R2 && N<12)? int2n(N & (~1UL)): powuu(N,N);
763 35 : if (cmpii(dK,c1) <= 0) return dft_bound();
764 :
765 35 : p1 = sqrr(glog(gdiv(dK,c1),DEFAULTPREC));
766 35 : p1 = divru(gmul2n(powru(divru(mulru(p1,3),N*(N*N-1)-6*R2),R),R2), N);
767 35 : p1 = sqrtr(gdiv(p1, Hermite_bound(R, DEFAULTPREC)));
768 35 : if (DEBUGLEVEL>1) err_printf("Mahler bound for regulator: %Ps\n",p1);
769 35 : return gmax_shallow(p1, dbltor(0.2));
770 : }
771 :
772 : static int
773 70553 : is_unit(GEN M, long r1, GEN x)
774 : {
775 70553 : pari_sp av = avma;
776 70553 : GEN Nx = ground( embed_norm(RgM_zc_mul(M,x), r1) );
777 70553 : return gc_bool(av, is_pm1(Nx));
778 : }
779 :
780 : /* True nf. FIXME: should use smallvectors */
781 : static double
782 42 : minimforunits(GEN nf, long BORNE, ulong w)
783 : {
784 42 : const long prec = MEDDEFAULTPREC;
785 42 : long n, r1, i, j, k, *x, cnt = 0;
786 42 : pari_sp av = avma;
787 : GEN r, M;
788 : double p, norme, normin;
789 : double **q,*v,*y,*z;
790 42 : double eps=0.000001, BOUND = BORNE * 1.00001;
791 :
792 42 : if (DEBUGLEVEL>=2)
793 : {
794 0 : err_printf("Searching minimum of T2-form on units:\n");
795 0 : if (DEBUGLEVEL>2) err_printf(" BOUND = %ld\n",BORNE);
796 : }
797 42 : n = nf_get_degree(nf); r1 = nf_get_r1(nf);
798 42 : minim_alloc(n+1, &q, &x, &y, &z, &v);
799 42 : M = gprec_w(nf_get_M(nf), prec);
800 42 : r = gaussred_from_QR(nf_get_G(nf), prec);
801 231 : for (j=1; j<=n; j++)
802 : {
803 189 : v[j] = gtodouble(gcoeff(r,j,j));
804 651 : for (i=1; i<j; i++) q[i][j] = gtodouble(gcoeff(r,i,j));
805 : }
806 42 : normin = (double)BORNE*(1-eps);
807 42 : k=n; y[n]=z[n]=0;
808 42 : x[n] = (long)(sqrt(BOUND/v[n]));
809 :
810 70553 : for(;;x[1]--)
811 : {
812 : do
813 : {
814 71901 : if (k>1)
815 : {
816 1348 : long l = k-1;
817 1348 : z[l] = 0;
818 5033 : for (j=k; j<=n; j++) z[l] += q[l][j]*x[j];
819 1348 : p = (double)x[k] + z[k];
820 1348 : y[l] = y[k] + p*p*v[k];
821 1348 : x[l] = (long)floor(sqrt((BOUND-y[l])/v[l])-z[l]);
822 1348 : k = l;
823 : }
824 : for(;;)
825 : {
826 73102 : p = (double)x[k] + z[k];
827 73102 : if (y[k] + p*p*v[k] <= BOUND) break;
828 1201 : k++; x[k]--;
829 : }
830 : }
831 71901 : while (k>1);
832 70595 : if (!x[1] && y[1]<=eps) break;
833 :
834 70567 : if (DEBUGLEVEL>8) err_printf(".");
835 70567 : if (++cnt == 5000) return -1.; /* too expensive */
836 :
837 70553 : p = (double)x[1] + z[1]; norme = y[1] + p*p*v[1];
838 70553 : if (is_unit(M, r1, x) && norme < normin)
839 : {
840 : /* exclude roots of unity */
841 56 : if (norme < 2*n)
842 : {
843 42 : GEN t = nfpow_u(nf, zc_to_ZC(x), w);
844 42 : if (typ(t) != t_COL || ZV_isscalar(t)) continue;
845 : }
846 21 : normin = norme*(1-eps);
847 21 : if (DEBUGLEVEL>=2) err_printf("*");
848 : }
849 : }
850 28 : if (DEBUGLEVEL>=2) err_printf("\n");
851 28 : set_avma(av);
852 28 : return normin;
853 : }
854 :
855 : #undef NBMAX
856 : static int
857 1804 : is_zero(GEN x, long bitprec) { return (gexpo(x) < -bitprec); }
858 :
859 : static int
860 1228 : is_complex(GEN x, long bitprec) { return !is_zero(imag_i(x), bitprec); }
861 :
862 : /* assume M_star t_REAL
863 : * FIXME: what does this do ? To be rewritten */
864 : static GEN
865 28 : compute_M0(GEN M_star,long N)
866 : {
867 : long m1,m2,n1,n2,n3,lr,lr1,lr2,i,j,l,vx,vy,vz,vM;
868 : GEN pol,p1,p2,p3,p4,p5,p6,p7,p8,p9,u,v,w,r,r1,r2,M0,M0_pro,S,P,M;
869 : GEN f1,f2,f3,g1,g2,g3,pg1,pg2,pg3,pf1,pf2,pf3,X,Y,Z;
870 28 : long bitprec = 24;
871 :
872 28 : if (N == 2) return gmul2n(sqrr(gacosh(gmul2n(M_star,-1),0)), -1);
873 21 : vx = fetch_var(); X = pol_x(vx);
874 21 : vy = fetch_var(); Y = pol_x(vy);
875 21 : vz = fetch_var(); Z = pol_x(vz);
876 21 : vM = fetch_var(); M = pol_x(vM);
877 :
878 21 : M0 = NULL; m1 = N/3;
879 56 : for (n1=1; n1<=m1; n1++) /* 1 <= n1 <= n2 <= n3 < N */
880 : {
881 35 : m2 = (N-n1)>>1;
882 112 : for (n2=n1; n2<=m2; n2++)
883 : {
884 77 : pari_sp av = avma; n3=N-n1-n2;
885 77 : if (n1==n2 && n1==n3) /* n1 = n2 = n3 = m1 = N/3 */
886 : {
887 7 : p1 = divru(M_star, m1);
888 7 : p4 = sqrtr_abs( mulrr(addsr(1,p1),subrs(p1,3)) );
889 7 : p5 = subrs(p1,1);
890 7 : u = gen_1;
891 7 : v = gmul2n(addrr(p5,p4),-1);
892 7 : w = gmul2n(subrr(p5,p4),-1);
893 7 : M0_pro=gmul2n(mulur(m1,addrr(sqrr(logr_abs(v)),sqrr(logr_abs(w)))), -2);
894 7 : if (DEBUGLEVEL>2)
895 0 : err_printf("[ %ld, %ld, %ld ]: %.28Pg\n",n1,n2,n3,M0_pro);
896 7 : if (!M0 || gcmp(M0_pro,M0) < 0) M0 = M0_pro;
897 : }
898 70 : else if (n1==n2 || n2==n3)
899 42 : { /* n3 > N/3 >= n1 */
900 42 : long k = N - 2*n2;
901 42 : p2 = deg1pol_shallow(stoi(-n2), M_star, vx); /* M* - n2 X */
902 42 : p3 = gmul(powuu(k,k),
903 : gpowgs(gsubgs(RgX_Rg_mul(p2, M_star), k*k), n2));
904 42 : pol = gsub(p3, RgX_mul(monomial(powuu(n2,n2), n2, vx),
905 : gpowgs(p2, N-n2)));
906 42 : r = roots(pol, DEFAULTPREC); lr = lg(r);
907 378 : for (i=1; i<lr; i++)
908 : {
909 : GEN n2S;
910 336 : S = real_i(gel(r,i));
911 336 : if (is_complex(gel(r,i), bitprec) || signe(S) <= 0) continue;
912 :
913 182 : n2S = mulur(n2,S);
914 182 : p4 = subrr(M_star, n2S);
915 182 : P = divrr(mulrr(n2S,p4), subrs(mulrr(M_star,p4),k*k));
916 182 : p5 = subrr(sqrr(S), gmul2n(P,2));
917 182 : if (gsigne(p5) < 0) continue;
918 :
919 140 : p6 = sqrtr(p5);
920 140 : v = gmul2n(subrr(S,p6),-1);
921 140 : if (gsigne(v) <= 0) continue;
922 :
923 126 : u = gmul2n(addrr(S,p6),-1);
924 126 : w = gpow(P, sstoQ(-n2,k), 0);
925 126 : p6 = mulur(n2, addrr(sqrr(logr_abs(u)), sqrr(logr_abs(v))));
926 126 : M0_pro = gmul2n(addrr(p6, mulur(k, sqrr(logr_abs(w)))),-2);
927 126 : if (DEBUGLEVEL>2)
928 0 : err_printf("[ %ld, %ld, %ld ]: %.28Pg\n",n1,n2,n3,M0_pro);
929 126 : if (!M0 || gcmp(M0_pro,M0) < 0) M0 = M0_pro;
930 : }
931 : }
932 : else
933 : {
934 28 : f1 = gsub(gadd(gmulsg(n1,X),gadd(gmulsg(n2,Y),gmulsg(n3,Z))), M);
935 28 : f2 = gmulsg(n1,gmul(Y,Z));
936 28 : f2 = gadd(f2,gmulsg(n2,gmul(X,Z)));
937 28 : f2 = gadd(f2,gmulsg(n3,gmul(X,Y)));
938 28 : f2 = gsub(f2,gmul(M,gmul(X,gmul(Y,Z))));
939 28 : f3 = gsub(gmul(gpowgs(X,n1),gmul(gpowgs(Y,n2),gpowgs(Z,n3))), gen_1);
940 : /* f1 = n1 X + n2 Y + n3 Z - M */
941 : /* f2 = n1 YZ + n2 XZ + n3 XY */
942 : /* f3 = X^n1 Y^n2 Z^n3 - 1*/
943 28 : g1=resultant(f1,f2); g1=primpart(g1);
944 28 : g2=resultant(f1,f3); g2=primpart(g2);
945 28 : g3=resultant(g1,g2); g3=primpart(g3);
946 28 : pf1=gsubst(f1,vM,M_star); pg1=gsubst(g1,vM,M_star);
947 28 : pf2=gsubst(f2,vM,M_star); pg2=gsubst(g2,vM,M_star);
948 28 : pf3=gsubst(f3,vM,M_star); pg3=gsubst(g3,vM,M_star);
949 : /* g3 = Res_Y,Z(f1,f2,f3) */
950 28 : r = roots(pg3,DEFAULTPREC); lr = lg(r);
951 476 : for (i=1; i<lr; i++)
952 : {
953 448 : w = real_i(gel(r,i));
954 448 : if (is_complex(gel(r,i), bitprec) || signe(w) <= 0) continue;
955 140 : p1=gsubst(pg1,vz,w);
956 140 : p2=gsubst(pg2,vz,w);
957 140 : p3=gsubst(pf1,vz,w);
958 140 : p4=gsubst(pf2,vz,w);
959 140 : p5=gsubst(pf3,vz,w);
960 140 : r1 = roots(p1, DEFAULTPREC); lr1 = lg(r1);
961 420 : for (j=1; j<lr1; j++)
962 : {
963 280 : v = real_i(gel(r1,j));
964 280 : if (is_complex(gel(r1,j), bitprec) || signe(v) <= 0
965 280 : || !is_zero(gsubst(p2,vy,v), bitprec)) continue;
966 :
967 164 : p7=gsubst(p3,vy,v);
968 164 : p8=gsubst(p4,vy,v);
969 164 : p9=gsubst(p5,vy,v);
970 164 : r2 = roots(p7, DEFAULTPREC); lr2 = lg(r2);
971 328 : for (l=1; l<lr2; l++)
972 : {
973 164 : u = real_i(gel(r2,l));
974 164 : if (is_complex(gel(r2,l), bitprec) || signe(u) <= 0
975 164 : || !is_zero(gsubst(p8,vx,u), bitprec)
976 164 : || !is_zero(gsubst(p9,vx,u), bitprec)) continue;
977 :
978 164 : M0_pro = mulur(n1, sqrr(logr_abs(u)));
979 164 : M0_pro = gadd(M0_pro, mulur(n2, sqrr(logr_abs(v))));
980 164 : M0_pro = gadd(M0_pro, mulur(n3, sqrr(logr_abs(w))));
981 164 : M0_pro = gmul2n(M0_pro,-2);
982 164 : if (DEBUGLEVEL>2)
983 0 : err_printf("[ %ld, %ld, %ld ]: %.28Pg\n",n1,n2,n3,M0_pro);
984 164 : if (!M0 || gcmp(M0_pro,M0) < 0) M0 = M0_pro;
985 : }
986 : }
987 : }
988 : }
989 77 : if (!M0) set_avma(av); else M0 = gc_GEN(av, M0);
990 : }
991 : }
992 105 : for (i=1;i<=4;i++) (void)delete_var();
993 21 : return M0? M0: gen_0;
994 : }
995 :
996 : static GEN
997 63 : lowerboundforregulator(GEN bnf, GEN units)
998 : {
999 63 : long i, N, R2, RU = lg(units)-1;
1000 : GEN nf, M0, M, G, minunit;
1001 : double bound;
1002 :
1003 63 : if (!RU) return gen_1;
1004 63 : nf = bnf_get_nf(bnf);
1005 63 : N = nf_get_degree(nf);
1006 63 : R2 = nf_get_r2(nf);
1007 :
1008 63 : G = nf_get_G(nf);
1009 63 : minunit = gnorml2(RgM_RgC_mul(G, gel(units,1))); /* T2(units[1]) */
1010 112 : for (i=2; i<=RU; i++)
1011 : {
1012 49 : GEN t = gnorml2(RgM_RgC_mul(G, gel(units,i)));
1013 49 : if (gcmp(t,minunit) < 0) minunit = t;
1014 : }
1015 63 : if (gexpo(minunit) > 30) return NULL;
1016 :
1017 42 : bound = minimforunits(nf, itos(gceil(minunit)), bnf_get_tuN(bnf));
1018 42 : if (bound < 0) return NULL;
1019 28 : if (DEBUGLEVEL>1) err_printf("M* = %Ps\n", dbltor(bound));
1020 28 : M0 = compute_M0(dbltor(bound), N);
1021 28 : if (DEBUGLEVEL>1) err_printf("M0 = %.28Pg\n",M0);
1022 28 : M = gmul2n(divru(gdiv(powrs(M0,RU),Hermite_bound(RU, DEFAULTPREC)),N),R2);
1023 28 : if (cmprr(M, dbltor(0.04)) < 0) return NULL;
1024 28 : M = sqrtr(M);
1025 28 : if (DEBUGLEVEL>1)
1026 0 : err_printf("(lower bound for regulator) M = %.28Pg\n",M);
1027 28 : return M;
1028 : }
1029 :
1030 : /* upper bound for the index of bnf.fu in the full unit group */
1031 : static GEN
1032 63 : bound_unit_index(GEN bnf, GEN units)
1033 : {
1034 63 : pari_sp av = avma;
1035 63 : GEN x = lowerboundforregulator(bnf, units);
1036 63 : if (!x) { set_avma(av); x = regulatorbound(bnf); }
1037 63 : return gc_INT(av, ground(gdiv(bnf_get_reg(bnf), x)));
1038 : }
1039 :
1040 : /* Compute a square matrix of rank #beta attached to a family
1041 : * (P_i), 1<=i<=#beta, of primes s.t. N(P_i) = 1 mod p, and
1042 : * (P_i,beta[j]) = 1 for all i,j. nf = true nf */
1043 : static void
1044 1715 : primecertify(GEN nf, GEN beta, ulong p, GEN bad)
1045 : {
1046 1715 : long lb = lg(beta), rmax = lb - 1;
1047 : GEN M, vQ, L;
1048 : ulong q;
1049 : forprime_t T;
1050 :
1051 1715 : if (p == 2)
1052 49 : L = cgetg(1,t_VECSMALL);
1053 : else
1054 1666 : L = mkvecsmall(p);
1055 1715 : (void)u_forprime_arith_init(&T, 1, ULONG_MAX, 1, p);
1056 1715 : M = cgetg(lb,t_MAT); setlg(M,1);
1057 3577 : while ((q = u_forprime_next(&T)))
1058 : {
1059 : GEN qq, gg, og;
1060 : long lQ, i, j;
1061 : ulong g, m;
1062 3577 : if (!umodiu(bad,q)) continue;
1063 :
1064 3283 : qq = utoipos(q);
1065 3283 : vQ = idealprimedec_limit_f(nf,qq,1);
1066 3283 : lQ = lg(vQ); if (lQ == 1) continue;
1067 :
1068 : /* cf rootsof1_Fl */
1069 2149 : g = pgener_Fl_local(q, L);
1070 2149 : m = (q-1) / p;
1071 2149 : gg = utoipos( Fl_powu(g, m, q) ); /* order p in (Z/q)^* */
1072 2149 : og = mkmat2(mkcol(utoi(p)), mkcol(gen_1)); /* order of g */
1073 :
1074 2149 : if (DEBUGLEVEL>3) err_printf(" generator of (Zk/Q)^*: %lu\n", g);
1075 2849 : for (i = 1; i < lQ; i++)
1076 : {
1077 2415 : GEN C = cgetg(lb, t_VECSMALL);
1078 2415 : GEN Q = gel(vQ,i); /* degree 1 */
1079 2415 : GEN modpr = zkmodprinit(nf, Q);
1080 : long r;
1081 :
1082 6923 : for (j = 1; j < lb; j++)
1083 : {
1084 4508 : GEN t = nf_to_Fp_coprime(nf, gel(beta,j), modpr);
1085 4508 : t = utoipos( Fl_powu(t[2], m, q) );
1086 4508 : C[j] = itou( Fp_log(t, gg, og, qq) ) % p;
1087 : }
1088 2415 : r = lg(M);
1089 2415 : gel(M,r) = C; setlg(M, r+1);
1090 2415 : if (Flm_rank(M, p) != r) { setlg(M,r); continue; }
1091 :
1092 2191 : if (DEBUGLEVEL>2)
1093 : {
1094 0 : if (DEBUGLEVEL>3)
1095 : {
1096 0 : err_printf(" prime ideal Q: %Ps\n",Q);
1097 0 : err_printf(" matrix log(b_j mod Q_i): %Ps\n", M);
1098 : }
1099 0 : err_printf(" new rank: %ld\n",r);
1100 : }
1101 2191 : if (r == rmax) return;
1102 : }
1103 : }
1104 0 : pari_err_BUG("primecertify");
1105 : }
1106 :
1107 : struct check_pr {
1108 : long w; /* #mu(K) */
1109 : GEN mu; /* generator of mu(K) */
1110 : GEN fu;
1111 : GEN cyc;
1112 : GEN cycgen;
1113 : GEN bad; /* p | bad <--> p | some element occurring in cycgen */
1114 : };
1115 :
1116 : static void
1117 1715 : check_prime(ulong p, GEN nf, struct check_pr *S)
1118 : {
1119 1715 : pari_sp av = avma;
1120 1715 : long i,b, lc = lg(S->cyc), lf = lg(S->fu);
1121 1715 : GEN beta = cgetg(lf+lc, t_VEC);
1122 :
1123 1715 : if (DEBUGLEVEL>1) err_printf(" *** testing p = %lu\n",p);
1124 1785 : for (b=1; b<lc; b++)
1125 : {
1126 1484 : if (umodiu(gel(S->cyc,b), p)) break; /* p \nmid cyc[b] */
1127 70 : if (b==1 && DEBUGLEVEL>2) err_printf(" p divides h(K)\n");
1128 70 : gel(beta,b) = gel(S->cycgen,b);
1129 : }
1130 1715 : if (S->w % p == 0)
1131 : {
1132 49 : if (DEBUGLEVEL>2) err_printf(" p divides w(K)\n");
1133 49 : gel(beta,b++) = S->mu;
1134 : }
1135 3787 : for (i=1; i<lf; i++) gel(beta,b++) = gel(S->fu,i);
1136 1715 : setlg(beta, b); /* beta = [cycgen[i] if p|cyc[i], tu if p|w, fu] */
1137 1715 : if (DEBUGLEVEL>3) err_printf(" Beta list = %Ps\n",beta);
1138 1715 : primecertify(nf, beta, p, S->bad); set_avma(av);
1139 1715 : }
1140 :
1141 : static void
1142 63 : init_bad(struct check_pr *S, GEN nf, GEN gen)
1143 : {
1144 63 : long i, l = lg(gen);
1145 63 : GEN bad = gen_1;
1146 :
1147 126 : for (i=1; i < l; i++)
1148 63 : bad = lcmii(bad, gcoeff(gel(gen,i),1,1));
1149 126 : for (i = 1; i < l; i++)
1150 : {
1151 63 : GEN c = gel(S->cycgen,i);
1152 : long j;
1153 63 : if (typ(c) == t_MAT)
1154 : {
1155 63 : GEN g = gel(c,1);
1156 420 : for (j = 1; j < lg(g); j++)
1157 : {
1158 357 : GEN h = idealhnf_shallow(nf, gel(g,j));
1159 357 : bad = lcmii(bad, gcoeff(h,1,1));
1160 : }
1161 : }
1162 : }
1163 63 : S->bad = bad;
1164 63 : }
1165 :
1166 : long
1167 63 : bnfcertify0(GEN bnf, long flag)
1168 : {
1169 63 : pari_sp av = avma;
1170 : long N;
1171 : GEN nf, cyc, B, U;
1172 : ulong bound, p;
1173 : struct check_pr S;
1174 : forprime_t T;
1175 :
1176 63 : bnf = checkbnf(bnf);
1177 63 : nf = bnf_get_nf(bnf);
1178 63 : N = nf_get_degree(nf); if (N==1) return 1;
1179 63 : B = zimmertbound(nf_get_disc(nf), N, nf_get_r2(nf));
1180 63 : if (is_bigint(B))
1181 0 : pari_warn(warner,"Zimmert's bound is large (%Ps), certification will take a long time", B);
1182 63 : if (!is_pm1(nf_get_index(nf)))
1183 : {
1184 42 : GEN D = nf_get_diff(nf), L;
1185 42 : if (DEBUGLEVEL>1) err_printf("**** Testing Different = %Ps\n",D);
1186 42 : L = bnfisprincipal0(bnf, D, nf_FORCE);
1187 42 : if (DEBUGLEVEL>1) err_printf(" is %Ps\n", L);
1188 : }
1189 63 : if (DEBUGLEVEL)
1190 : {
1191 0 : err_printf("PHASE 1 [CLASS GROUP]: are all primes good ?\n");
1192 0 : err_printf(" Testing primes <= %Ps\n", B);
1193 : }
1194 63 : bnftestprimes(bnf, B);
1195 63 : if (flag) return 1;
1196 :
1197 63 : U = bnf_build_units(bnf);
1198 63 : cyc = bnf_get_cyc(bnf);
1199 63 : S.w = bnf_get_tuN(bnf);
1200 63 : S.mu = gel(U,1);
1201 63 : S.fu = vecslice(U,2,lg(U)-1);
1202 63 : S.cyc = cyc;
1203 63 : S.cycgen = bnf_build_cycgen(bnf);
1204 63 : init_bad(&S, nf, bnf_get_gen(bnf));
1205 :
1206 63 : B = bound_unit_index(bnf, S.fu);
1207 63 : if (DEBUGLEVEL)
1208 : {
1209 0 : err_printf("PHASE 2 [UNITS/RELATIONS]: are all primes good ?\n");
1210 0 : err_printf(" Testing primes <= %Ps\n", B);
1211 : }
1212 63 : bound = itou_or_0(B);
1213 63 : if (!bound) pari_err_OVERFLOW("bnfcertify [too many primes to check]");
1214 63 : if (u_forprime_init(&T, 2, bound))
1215 1757 : while ( (p = u_forprime_next(&T)) ) check_prime(p, nf, &S);
1216 63 : if (lg(cyc) > 1)
1217 : {
1218 28 : GEN f = Z_factor(cyc_get_expo(cyc)), P = gel(f,1);
1219 : long i;
1220 28 : if (DEBUGLEVEL>1) err_printf(" Primes dividing h(K)\n\n");
1221 35 : for (i = lg(P)-1; i; i--)
1222 : {
1223 28 : p = itou(gel(P,i)); if (p <= bound) break;
1224 7 : check_prime(p, nf, &S);
1225 : }
1226 : }
1227 63 : return gc_long(av,1);
1228 : }
1229 : long
1230 35 : bnfcertify(GEN bnf) { return bnfcertify0(bnf, 0); }
1231 :
1232 : /*******************************************************************/
1233 : /* */
1234 : /* RAY CLASS FIELDS: CONDUCTORS AND DISCRIMINANTS */
1235 : /* */
1236 : /*******************************************************************/
1237 : static GEN
1238 609 : bnrchar_i(GEN bnr, GEN g, GEN v)
1239 : {
1240 609 : long i, h, l = lg(g), t = typ_NULL;
1241 609 : GEN CH, D, U, U2, H, cycD, dv, dchi, cyc = NULL;
1242 :
1243 609 : if (checkbnr_i(bnr)) { t = typ_BNR; cyc = bnr_get_cyc(bnr); }
1244 14 : else if (checkznstar_i(bnr)) { t = typ_BIDZ; cyc = znstar_get_cyc(bnr); }
1245 0 : else if (typ(bnr) == t_VEC && RgV_is_ZV(bnr)) cyc = bnr;
1246 0 : else pari_err_TYPE("bnrchar", bnr);
1247 609 : switch(typ(g))
1248 : {
1249 : GEN G;
1250 28 : case t_VEC:
1251 28 : G = cgetg(l, t_MAT);
1252 28 : if (t == typ_BNR)
1253 : {
1254 49 : for (i = 1; i < l; i++) gel(G,i) = isprincipalray(bnr, gel(g,i));
1255 14 : cyc = bnr_get_cyc(bnr);
1256 : }
1257 : else
1258 35 : for (i = 1; i < l; i++) gel(G,i) = Zideallog(bnr, gel(g,i));
1259 28 : g = G; break;
1260 581 : case t_MAT:
1261 581 : if (RgM_is_ZM(g)) break;
1262 : default:
1263 0 : pari_err_TYPE("bnrchar",g);
1264 : }
1265 609 : H = ZM_hnfall_i(shallowconcat(g,diagonal_shallow(cyc)), v? &U: NULL, 1);
1266 609 : dv = NULL;
1267 609 : if (v)
1268 : {
1269 42 : GEN w = Q_remove_denom(v, &dv);
1270 42 : if (typ(v)!=t_VEC || lg(v)!=l || !RgV_is_ZV(w)) pari_err_TYPE("bnrchar",v);
1271 42 : if (!dv) v = NULL;
1272 : else
1273 : {
1274 42 : U = rowslice(U, 1, l-1);
1275 42 : w = FpV_red(ZV_ZM_mul(w, U), dv);
1276 140 : for (i = 1; i < l; i++)
1277 105 : if (signe(gel(w,i))) pari_err_TYPE("bnrchar [inconsistent values]",v);
1278 35 : v = vecslice(w,l,lg(w)-1);
1279 : }
1280 : }
1281 : /* chi defined on subgroup H, chi(H[i]) = e(v[i] / dv)
1282 : * unless v = NULL: chi|H = 1*/
1283 602 : h = itos( ZM_det_triangular(H) ); /* #(clgp/H) = number of chars */
1284 602 : if (h == 1) /* unique character, H = Id */
1285 : {
1286 14 : if (v)
1287 14 : v = char_denormalize(cyc,dv,v);
1288 : else
1289 0 : v = zerovec(lg(cyc)-1); /* trivial char */
1290 14 : return mkvec(v);
1291 : }
1292 :
1293 : /* chi defined on a subgroup of index h > 1; U H V = D diagonal,
1294 : * Z^#H / (H) = Z^#H / (D) ~ \oplus (Z/diZ) */
1295 588 : D = ZM_snfall_i(H, &U, NULL, 1);
1296 588 : cycD = cyc_normalize(D); gel(cycD,1) = gen_1; /* cycD[i] = d1/di */
1297 588 : dchi = gel(D,1);
1298 588 : U2 = ZM_diag_mul(cycD, U);
1299 588 : if (v)
1300 : {
1301 21 : GEN Ui = ZM_inv(U, NULL);
1302 21 : GEN Z = hnf_solve(H, ZM_mul_diag(Ui, D));
1303 21 : v = ZV_ZM_mul(ZV_ZM_mul(v, Z), U2);
1304 21 : dchi = mulii(dchi, dv);
1305 21 : U2 = ZM_Z_mul(U2, dv);
1306 : }
1307 588 : CH = cyc2elts(D);
1308 2996 : for (i = 1; i <= h; i++)
1309 : {
1310 2408 : GEN c = zv_ZM_mul(gel(CH,i), U2);
1311 2408 : if (v) c = ZC_add(c, v);
1312 2408 : gel(CH,i) = char_denormalize(cyc, dchi, c);
1313 : }
1314 588 : return CH;
1315 : }
1316 : GEN
1317 609 : bnrchar(GEN bnr, GEN g, GEN v)
1318 : {
1319 609 : pari_sp av = avma;
1320 609 : return gc_GEN(av, bnrchar_i(bnr,g,v));
1321 : }
1322 :
1323 : /* Let bnr1, bnr2 be such that mod(bnr2) | mod(bnr1), compute surjective map
1324 : * p: Cl(bnr1) ->> Cl(bnr2).
1325 : * Write (bnr gens) for the concatenation of the bnf [corrected by El] and bid
1326 : * generators; and bnr.gen for the SNF generators. Then
1327 : * bnr.gen = (bnf.gen*bnr.El | bid.gen) bnr.Ui
1328 : * (bnf.gen*bnr.El | bid.gen) = bnr.gen * bnr.U */
1329 : GEN
1330 15274 : bnrsurjection(GEN bnr1, GEN bnr2)
1331 : {
1332 15274 : GEN bnf = bnr_get_bnf(bnr2), nf = bnf_get_nf(bnf);
1333 15274 : GEN M, U = bnr_get_U(bnr2), bid2 = bnr_get_bid(bnr2);
1334 15274 : GEN gen1 = bid_get_gen(bnr_get_bid(bnr1));
1335 15274 : GEN cyc2 = bnr_get_cyc(bnr2), e2 = cyc_get_expo(cyc2);
1336 15274 : long i, l = lg(bnf_get_cyc(bnf)), lb = lg(gen1);
1337 : /* p(bnr1.gen) = p(bnr1 gens) * bnr1.Ui
1338 : * = (bnr2 gens) * P * bnr1.Ui
1339 : * = bnr2.gen * (bnr2.U * P * bnr1.Ui) */
1340 :
1341 : /* p(bid1.gen) on bid2.gen */
1342 15274 : M = cgetg(lb, t_MAT);
1343 111853 : for (i = 1; i < lb; i++) gel(M,i) = ideallogmod(nf, gel(gen1,i), bid2, e2);
1344 : /* [U[1], U[2]] * [Id, 0; N, M] = [U[1] + U[2]*N, U[2]*M] */
1345 15274 : M = ZM_mul(gel(U,2), M);
1346 15274 : if (l > 1)
1347 : { /* non trivial class group */
1348 : /* p(bnf.gen * bnr1.El) in terms of bnf.gen * bnr2.El and bid2.gen */
1349 861 : GEN El2 = bnr_get_El(bnr2), El1 = bnr_get_El(bnr1);
1350 861 : long ngen2 = lg(bid_get_gen(bid2))-1;
1351 861 : if (!ngen2)
1352 602 : M = gel(U,1);
1353 : else
1354 : {
1355 259 : GEN U1 = gel(U,1), U2 = gel(U,2), T = cgetg(l, t_MAT);
1356 : /* T = U1 + U2 log(El2/El1) */
1357 539 : for (i = 1; i < l; i++)
1358 : { /* bnf gen in bnr1 is bnf.gen * El1 = bnf gen in bnr 2 * El1/El2 */
1359 280 : GEN c = gel(U1,i);
1360 280 : if (typ(gel(El1,i)) != t_INT) /* else El1[i] = 1 => El2[i] = 1 */
1361 : {
1362 98 : GEN z = nfdiv(nf,gel(El1,i),gel(El2,i));
1363 98 : c = ZC_add(c, ZM_ZC_mul(U2, ideallogmod(nf, z, bid2, e2)));
1364 : }
1365 280 : gel(T,i) = c;
1366 : }
1367 259 : M = shallowconcat(T, M);
1368 : }
1369 : }
1370 15274 : M = ZM_ZV_mod(ZM_mul(M, bnr_get_Ui(bnr1)), cyc2);
1371 15274 : return mkvec3(M, bnr_get_cyc(bnr1), cyc2);
1372 : }
1373 :
1374 : /* nchi a normalized character, S a surjective map ; return S(nchi)
1375 : * still normalized wrt the original cyclic structure (S[2]) */
1376 : GEN
1377 1449 : abmap_nchar_image(GEN S, GEN nchi)
1378 : {
1379 1449 : GEN U, M = gel(S,1), Mc = diagonal_shallow(gel(S,3));
1380 1449 : long l = lg(M);
1381 :
1382 1449 : (void)ZM_hnfall_i(shallowconcat(M, Mc), &U, 1); /* identity */
1383 1449 : U = matslice(U,1,l-1, l,lg(U)-1);
1384 1449 : return char_simplify(gel(nchi,1), ZV_ZM_mul(gel(nchi,2), U));
1385 : }
1386 : GEN
1387 1232 : abmap_char_image(GEN S, GEN chi)
1388 : {
1389 1232 : GEN nchi = char_normalize(chi, cyc_normalize(gel(S,2)));
1390 1232 : GEN DC = abmap_nchar_image(S, nchi);
1391 1232 : return char_denormalize(gel(S,3), gel(DC,1), gel(DC,2));
1392 : }
1393 :
1394 : GEN
1395 616 : bnrmap(GEN A, GEN B)
1396 : {
1397 616 : pari_sp av = avma;
1398 : GEN KA, KB, M, c, C;
1399 616 : if ((KA = checkbnf_i(A)))
1400 : {
1401 168 : checkbnr(A); checkbnr(B); KB = bnr_get_bnf(B);
1402 168 : if (!gidentical(KA, KB))
1403 0 : pari_err_TYPE("bnrmap [different fields]", mkvec2(KA,KB));
1404 168 : return gc_GEN(av, bnrsurjection(A,B));
1405 : }
1406 448 : if (lg(A) != 4 || typ(A) != t_VEC) pari_err_TYPE("bnrmap [not a map]", A);
1407 441 : M = gel(A,1); c = gel(A,2); C = gel(A,3);
1408 441 : if (typ(M) != t_MAT || !RgM_is_ZM(M) || typ(c) != t_VEC ||
1409 441 : typ(C) != t_VEC || lg(c) != lg(M) || (lg(M) > 1 && lgcols(M) != lg(C)))
1410 0 : pari_err_TYPE("bnrmap [not a map]", A);
1411 441 : switch(typ(B))
1412 : {
1413 7 : case t_INT: /* subgroup */
1414 7 : B = scalarmat_shallow(B, lg(C)-1);
1415 7 : B = ZM_hnfmodid(B, C); break;
1416 392 : case t_MAT: /* subgroup */
1417 392 : if (!RgM_is_ZM(B)) pari_err_TYPE("bnrmap [not a subgroup]", B);
1418 385 : B = ZM_hnfmodid(B, c); B = abmap_subgroup_image(A, B); break;
1419 21 : case t_VEC: /* character */
1420 21 : if (!char_check(c, B))
1421 14 : pari_err_TYPE("bnrmap [not a character mod mA]", B);
1422 7 : B = abmap_char_image(A, B); break;
1423 21 : case t_COL: /* discrete log mod mA */
1424 21 : if (lg(B) != lg(c) || !RgV_is_ZV(B))
1425 14 : pari_err_TYPE("bnrmap [not a discrete log]", B);
1426 7 : B = ZV_ZV_mod(ZM_ZC_mul(M, B), C);
1427 7 : return gc_upto(av, B);
1428 : }
1429 392 : return gc_GEN(av, B);
1430 : }
1431 :
1432 : /* convert A,B,C to [bnr, H] */
1433 : GEN
1434 273 : ABC_to_bnr(GEN A, GEN B, GEN C, GEN *H, int gen)
1435 : {
1436 273 : if (typ(A) == t_VEC)
1437 273 : switch(lg(A))
1438 : {
1439 119 : case 7: /* bnr */
1440 119 : *H = B; return A;
1441 154 : case 11: /* bnf */
1442 154 : if (!B) pari_err_TYPE("ABC_to_bnr [bnf+missing conductor]",A);
1443 154 : *H = C; return Buchray(A,B, gen? nf_INIT | nf_GEN: nf_INIT);
1444 : }
1445 0 : pari_err_TYPE("ABC_to_bnr",A);
1446 : *H = NULL; return NULL; /* LCOV_EXCL_LINE */
1447 : }
1448 :
1449 : /* OBSOLETE */
1450 : GEN
1451 63 : bnrconductor0(GEN A, GEN B, GEN C, long flag)
1452 : {
1453 63 : pari_sp av = avma;
1454 63 : GEN H, bnr = ABC_to_bnr(A,B,C,&H, 0);
1455 63 : return gc_GEN(av, bnrconductor(bnr, H, flag));
1456 : }
1457 :
1458 : long
1459 35 : bnrisconductor0(GEN A,GEN B,GEN C)
1460 : {
1461 35 : GEN H, bnr = ABC_to_bnr(A,B,C,&H, 0);
1462 35 : return bnrisconductor(bnr, H);
1463 : }
1464 :
1465 : static GEN
1466 613089 : ideallog_to_bnr_i(GEN Ubid, GEN cyc, GEN z)
1467 613089 : { return (lg(Ubid)==1)? zerocol(lg(cyc)-1): ZV_ZV_mod(ZM_ZC_mul(Ubid,z), cyc); }
1468 : /* return bnrisprincipal(bnr, (t)), assuming x = ideallog(t); allow a
1469 : * t_MAT for x, understood as a collection of ideallog(t_i) */
1470 : static GEN
1471 596257 : ideallog_to_bnr(GEN bnr, GEN x)
1472 : {
1473 596257 : GEN U = gel(bnr_get_U(bnr), 2); /* bid part */
1474 596256 : GEN cyc = bnr_get_cyc(bnr);
1475 596270 : if (typ(x) == t_COL) return ideallog_to_bnr_i(U, cyc, x);
1476 985339 : pari_APPLY_same(ideallog_to_bnr_i(U, cyc, gel(x,i)));
1477 : }
1478 : static GEN
1479 484276 : bnr_log_gen_pr(GEN bnr, zlog_S *S, long e, long index)
1480 484276 : { return ideallog_to_bnr(bnr, log_gen_pr(S, index, bnr_get_nf(bnr), e)); }
1481 : static GEN
1482 111991 : bnr_log_gen_arch(GEN bnr, zlog_S *S, long index)
1483 111991 : { return ideallog_to_bnr(bnr, log_gen_arch(S, index)); }
1484 :
1485 : /* A \subset H ? Allow H = NULL = trivial subgroup */
1486 : static int
1487 498875 : contains(GEN H, GEN A)
1488 498875 : { return H? (hnf_solve(H, A) != NULL): gequal0(A); }
1489 :
1490 : /* finite part of the conductor of H is S.P^e2*/
1491 : static GEN
1492 62468 : cond0_e(GEN bnr, GEN H, zlog_S *S)
1493 : {
1494 62468 : long j, k, l = lg(S->k), iscond0 = S->no2;
1495 62468 : GEN e = S->k, e2 = cgetg(l, t_COL);
1496 190620 : for (k = 1; k < l; k++)
1497 : {
1498 162819 : for (j = itos(gel(e,k)); j > 0; j--)
1499 : {
1500 147076 : if (!contains(H, bnr_log_gen_pr(bnr, S, j, k))) break;
1501 34664 : iscond0 = 0;
1502 : }
1503 128149 : gel(e2,k) = utoi(j);
1504 : }
1505 62465 : return iscond0? NULL: e2;
1506 : }
1507 : /* infinite part of the conductor of H in archp form */
1508 : static GEN
1509 62468 : condoo_archp(GEN bnr, GEN H, zlog_S *S)
1510 : {
1511 : long j, k, l;
1512 62468 : GEN archp = S->archp, archp2 = cgetg_copy(archp, &l);
1513 106337 : for (k = j = 1; k < l; k++)
1514 : {
1515 43869 : if (!contains(H, bnr_log_gen_arch(bnr, S, k)))
1516 : {
1517 32158 : archp2[j++] = archp[k];
1518 32158 : continue;
1519 : }
1520 : }
1521 62468 : if (j == l) return NULL;
1522 9618 : setlg(archp2, j); return archp2;
1523 : }
1524 :
1525 : /* true bnr, H subgroup */
1526 : GEN
1527 62468 : bnrconductor_factored_arch(GEN bnr, GEN H, GEN *parch)
1528 : {
1529 62468 : GEN nf = bnr_get_nf(bnr), bid = bnr_get_bid(bnr), e;
1530 : zlog_S S;
1531 :
1532 62468 : init_zlog(&S, bid); e = cond0_e(bnr, H, &S); /* in terms of S.P */
1533 62468 : if (parch)
1534 : {
1535 62468 : GEN archp = condoo_archp(bnr, H, &S);
1536 62468 : *parch = archp? indices_to_vec01(archp, nf_get_r1(nf)): NULL;
1537 : }
1538 62468 : return e? famat_remove_trivial(mkmat2(S.P, e)): NULL;
1539 : }
1540 : /* we return [factor(f0),foo] or NULL if cond(H) = cond(bnr).
1541 : * If flag, return f = [f0,foo] else return [f, factor(f0)]. */
1542 : static GEN
1543 6006 : bnrconductor_factored_i(GEN bnr, GEN H, long flag)
1544 : {
1545 6006 : GEN arch, fa, cond = NULL;
1546 :
1547 6006 : checkbnr(bnr); H = bnr_subgroup_check(bnr, H, NULL);
1548 6006 : fa = bnrconductor_factored_arch(bnr, H, &arch);
1549 6006 : if (!arch)
1550 : {
1551 3801 : GEN mod = bnr_get_mod(bnr);
1552 3801 : if (!fa) cond = mod;
1553 3801 : arch = gel(mod,2);
1554 : }
1555 6006 : if (!cond)
1556 : {
1557 2247 : GEN f0 = fa? factorbackprime(bnr_get_nf(bnr), gel(fa,1), gel(fa,2))
1558 2898 : : bid_get_ideal(bnr_get_bid(bnr));
1559 2898 : cond = mkvec2(f0, arch);
1560 : }
1561 6006 : if (flag) return cond;
1562 623 : if (!fa) fa = bid_get_fact(bnr_get_bid(bnr));
1563 623 : return mkvec2(cond, fa);
1564 : }
1565 : GEN
1566 623 : bnrconductor_factored(GEN bnr, GEN H)
1567 623 : { return bnrconductor_factored_i(bnr, H, 0); }
1568 : GEN
1569 5383 : bnrconductor_raw(GEN bnr, GEN H)
1570 5383 : { return bnrconductor_factored_i(bnr, H, 1); }
1571 :
1572 : /* exponent G/H, assuming H is a left divisor of matdiagonal(G.cyc) */
1573 : static GEN
1574 6041 : quotient_expo(GEN H)
1575 : {
1576 6041 : GEN D = ZM_snf(H); /* structure of G/H */
1577 6041 : return lg(D) == 1? gen_1: gel(D,1);
1578 : }
1579 :
1580 : /* H subgroup or NULL (trivial) */
1581 : GEN
1582 56462 : bnrtoprimitive(GEN bnr, GEN H, GEN mod)
1583 : {
1584 : long fl;
1585 56462 : GEN arch, fa = bnrconductor_factored_arch(bnr, H, &arch);
1586 56462 : if (!fa && !arch) return NULL;
1587 15190 : if (!fa) fa = bid_get_fact(bnr_get_bid(bnr));
1588 15190 : if (!arch) arch = gel(bnr_get_mod(bnr), 2);
1589 15190 : fl = lg(bnr_get_clgp(bnr)) == 4? nf_INIT | nf_GEN: nf_INIT;
1590 15190 : return Buchraymod_i(bnr, mkvec2(fa,arch), fl, mod);
1591 : }
1592 : void
1593 41846 : bnr_sanitize(GEN *pbnr)
1594 : {
1595 41846 : switch(nftyp(*pbnr))
1596 : {
1597 133 : case typ_BNF: *pbnr = Buchray(*pbnr, gen_1, nf_INIT);
1598 41832 : case typ_BNR: break;
1599 14 : default: checkbnr(*pbnr); /* error out */
1600 : }
1601 41832 : }
1602 : void
1603 9247 : bnr_subgroup_sanitize(GEN *pbnr, GEN *pH)
1604 : {
1605 9247 : GEN mod, cyc, bnrc, bnr = *pbnr, H = *pH;
1606 :
1607 9247 : bnr_sanitize(&bnr);
1608 9233 : cyc = bnr_get_cyc(bnr);
1609 9233 : if (!H) mod = cyc_get_expo(cyc);
1610 8827 : else switch(typ(H))
1611 : {
1612 2765 : case t_INT:
1613 2765 : mod = H; H = bnr_subgroup_check(bnr, H, NULL);
1614 2765 : break;
1615 7 : case t_VEC:
1616 7 : if (!char_check(cyc, H))
1617 0 : pari_err_TYPE("bnr_subgroup_sanitize [character]", H);
1618 7 : H = charker(cyc, H); /* character -> subgroup */
1619 7 : mod = quotient_expo(H);
1620 7 : break;
1621 6048 : case t_MAT:
1622 6048 : H = hnfmodid(H, cyc); /* make sure H is a left divisor of Mat(cyc) */
1623 6034 : mod = quotient_expo(H);
1624 6034 : break;
1625 7 : default: pari_err_TYPE("bnr_subroup_sanitize [subgroup]", H);
1626 0 : mod = NULL;
1627 : }
1628 9212 : bnrc = bnrtoprimitive(bnr, H, mod);
1629 9212 : if (!bnrc) bnrc = bnr;
1630 5243 : else if (H)
1631 : {
1632 4578 : GEN map = bnrsurjection(bnr, bnrc);
1633 4578 : H = abmap_subgroup_image(map, H);
1634 : }
1635 9212 : if (!H) H = diagonal_shallow(bnr_get_cyc(bnrc));
1636 9212 : *pbnr = bnrc; *pH = H;
1637 9212 : }
1638 : void
1639 1267 : bnr_char_sanitize(GEN *pbnr, GEN *pchi)
1640 : {
1641 1267 : GEN cyc, bnrc, bnr = *pbnr, chi = *pchi;
1642 :
1643 1267 : bnr_sanitize(&bnr); cyc = bnr_get_cyc(bnr);
1644 1267 : if (!char_check(cyc, chi))
1645 0 : pari_err_TYPE("bnr_char_sanitize [character]", chi);
1646 1267 : bnrc = bnrtoprimitive(bnr, charker(cyc,chi), charorder(cyc,chi));
1647 1267 : if (bnrc)
1648 : {
1649 714 : GEN map = bnrsurjection(bnr, bnrc);
1650 714 : *pbnr = bnrc;
1651 714 : *pchi = abmap_char_image(map, chi);
1652 : }
1653 1267 : }
1654 :
1655 : /* assume lg(CHI) > 1 */
1656 : void
1657 350 : bnr_vecchar_sanitize(GEN *pbnr, GEN *pCHI)
1658 : {
1659 350 : GEN D, nchi, map, H, cyc, o, bnrc, bnr = *pbnr, CHI = *pCHI;
1660 350 : long i, l = lg(CHI);
1661 :
1662 350 : bnr_sanitize(&bnr);
1663 350 : cyc = bnr_get_cyc(bnr); o = gen_1;
1664 1519 : for (i = 1; i < l; i++)
1665 : {
1666 1169 : GEN chi = gel(CHI,i);
1667 1169 : if (!char_check(cyc, chi))
1668 0 : pari_err_TYPE("bnr_char_sanitize [character]", chi);
1669 1169 : o = lcmii(o, charorder(cyc, chi));
1670 : }
1671 350 : H = bnr_subgroup_check(bnr, gel(CHI,1), NULL);
1672 350 : bnrc = bnrtoprimitive(bnr, H, o);
1673 350 : map = bnrc? bnrsurjection(bnr, bnrc): NULL;
1674 350 : if (!bnrc) bnrc = bnr;
1675 350 : D = cyc_normalize(bnr_get_cyc(bnrc));
1676 350 : nchi = cgetg(l, t_VEC);
1677 1512 : for (i = 1; i < l; i++)
1678 : {
1679 1169 : GEN chi = gel(CHI,i);
1680 1169 : if (map) chi = abmap_char_image(map, chi);
1681 1169 : if (i > 1 && !bnrisconductor(bnrc, chi))
1682 7 : pari_err_TYPE("bnr_vecchar_sanitize [different conductors]", CHI);
1683 1162 : gel(nchi, i) = char_renormalize(char_normalize(chi, D), o);
1684 : }
1685 343 : *pbnr = bnrc; *pCHI = mkvec2(o, nchi);
1686 343 : }
1687 :
1688 : /* Given a bnr and a subgroup H0 (possibly given as a
1689 : * character chi, in which case H0 = ker chi) of the ray class group,
1690 : * return [conductor(H0), bnr attached to conductor, image of H0].
1691 : * OBSOLETE: the interface is awkward and inefficient; split the 3 parts
1692 : * in caller 1) bnrtoprimitive with proper MOD depending on 3), 2) bnrmap,
1693 : * 3) image of objects. Don't use this wrapper except for lazy GP use. */
1694 : static GEN
1695 532 : bnrconductormod(GEN bnr, GEN H0, GEN MOD)
1696 : {
1697 532 : GEN H = bnr_subgroup_check(bnr, H0, NULL);
1698 532 : GEN bnrc = bnrtoprimitive(bnr, H, MOD);
1699 532 : int ischi = H0 && typ(H0) == t_VEC; /* character or subgroup ? */
1700 532 : if (!bnrc)
1701 : { /* same conductor */
1702 196 : bnrc = bnr;
1703 196 : if (ischi) H = H0;
1704 : }
1705 : else
1706 : {
1707 336 : if (ischi)
1708 : {
1709 0 : GEN map = bnrsurjection(bnr, bnrc);
1710 0 : H = abmap_char_image(map, H0);
1711 : }
1712 336 : else if (H)
1713 : {
1714 182 : GEN map = bnrsurjection(bnr, bnrc);
1715 182 : H = abmap_subgroup_image(map, H);
1716 : }
1717 : }
1718 532 : if (!H) H = diagonal_shallow(bnr_get_cyc(bnrc));
1719 532 : return mkvec3(bnr_get_mod(bnrc), bnrc, H);
1720 : }
1721 : /* OBSOLETE */
1722 : GEN
1723 616 : bnrconductor(GEN bnr, GEN H, long flag)
1724 : {
1725 616 : pari_sp av = avma;
1726 : GEN v;
1727 616 : if (flag == 0) return bnrconductor_raw(bnr, H);
1728 0 : if (flag > 2 || flag < 0) pari_err_FLAG("bnrconductor");
1729 0 : v = bnrconductormod(bnr, H, NULL); bnr = gel(v,1); H = gel(v,2);
1730 0 : if (flag == 1) gel(v,2) = bnr_get_clgp(gel(v,2));
1731 0 : return gc_GEN(av, v);
1732 : }
1733 :
1734 : long
1735 275525 : bnrisconductor(GEN bnr, GEN H0)
1736 : {
1737 275525 : pari_sp av = avma;
1738 : long j, k, l;
1739 : GEN archp, e, H;
1740 : zlog_S S;
1741 :
1742 275525 : checkbnr(bnr);
1743 275525 : init_zlog(&S, bnr_get_bid(bnr));
1744 275528 : if (!S.no2) return 0;
1745 232184 : H = bnr_subgroup_check(bnr, H0, NULL);
1746 :
1747 232184 : archp = S.archp;
1748 232184 : e = S.k; l = lg(e);
1749 370490 : for (k = 1; k < l; k++)
1750 : {
1751 262505 : j = itos(gel(e,k));
1752 262506 : if (contains(H, bnr_log_gen_pr(bnr, &S, j, k))) return gc_long(av,0);
1753 : }
1754 107985 : l = lg(archp);
1755 137251 : for (k = 1; k < l; k++)
1756 45164 : if (contains(H, bnr_log_gen_arch(bnr, &S, k))) return gc_long(av,0);
1757 92087 : return gc_long(av,1);
1758 : }
1759 :
1760 : /* return the norm group corresponding to the relative extension given by
1761 : * polrel over bnr.bnf, assuming it is abelian and the modulus of bnr is a
1762 : * multiple of the conductor */
1763 : static GEN
1764 826 : rnfnormgroup_i(GEN bnr, GEN polrel)
1765 : {
1766 : long i, j, degrel, degnf, k;
1767 : GEN bnf, index, discnf, nf, G, detG, fa, gdegrel;
1768 : GEN fac, col, cnd;
1769 : forprime_t S;
1770 : ulong p;
1771 :
1772 826 : checkbnr(bnr); bnf = bnr_get_bnf(bnr);
1773 826 : nf = bnf_get_nf(bnf);
1774 826 : cnd = gel(bnr_get_mod(bnr), 1);
1775 826 : polrel = RgX_nffix("rnfnormgroup", nf_get_pol(nf),polrel,1);
1776 826 : if (!gequal1(leading_coeff(polrel)))
1777 0 : pari_err_IMPL("rnfnormgroup for nonmonic polynomials");
1778 :
1779 826 : degrel = degpol(polrel);
1780 826 : if (umodiu(bnr_get_no(bnr), degrel)) return NULL;
1781 : /* degrel-th powers are in norm group */
1782 812 : gdegrel = utoipos(degrel);
1783 812 : G = ZV_snf_gcd(bnr_get_cyc(bnr), gdegrel);
1784 812 : detG = ZV_prod(G);
1785 812 : k = abscmpiu(detG,degrel);
1786 812 : if (k < 0) return NULL;
1787 812 : if (!k) return diagonal(G);
1788 :
1789 266 : G = diagonal_shallow(G);
1790 266 : discnf = nf_get_disc(nf);
1791 266 : index = nf_get_index(nf);
1792 266 : degnf = nf_get_degree(nf);
1793 266 : u_forprime_init(&S, 2, ULONG_MAX);
1794 1582 : while ( (p = u_forprime_next(&S)) )
1795 : {
1796 : long oldf, nfa;
1797 : /* If all pr are unramified and have the same residue degree, p =prod pr
1798 : * and including last pr^f or p^f is the same, but the last isprincipal
1799 : * is much easier! oldf is used to track this */
1800 :
1801 1582 : if (!umodiu(index, p)) continue; /* can't be treated efficiently */
1802 :
1803 : /* primes of degree 1 are enough, and simpler */
1804 1582 : fa = idealprimedec_limit_f(nf, utoipos(p), 1);
1805 1582 : nfa = lg(fa)-1;
1806 1582 : if (!nfa) continue;
1807 : /* all primes above p included ? */
1808 1295 : oldf = (nfa == degnf)? -1: 0;
1809 2471 : for (i=1; i<=nfa; i++)
1810 : {
1811 1442 : GEN pr = gel(fa,i), pp, T, polr, modpr;
1812 : long f, nfac;
1813 : /* if pr (probably) ramified, we have to use all (unramified) P | pr */
1814 1953 : if (idealval(nf,cnd,pr)) { oldf = 0; continue; }
1815 1099 : modpr = zk_to_Fq_init(nf, &pr, &T, &pp); /* T = NULL, pp ignored */
1816 1099 : polr = nfX_to_FqX(polrel, nf, modpr); /* in Fp[X] */
1817 1099 : polr = ZX_to_Flx(polr, p);
1818 1099 : if (!Flx_is_squarefree(polr, p)) { oldf = 0; continue; }
1819 :
1820 1043 : fac = gel(Flx_factor(polr, p), 1);
1821 1043 : f = degpol(gel(fac,1));
1822 1043 : if (f == degrel) continue; /* degrel-th powers already included */
1823 588 : nfac = lg(fac)-1;
1824 : /* check decomposition of pr has Galois type */
1825 1638 : for (j=2; j<=nfac; j++)
1826 1316 : if (degpol(gel(fac,j)) != f) return NULL;
1827 581 : if (oldf < 0) oldf = f; else if (oldf != f) oldf = 0;
1828 :
1829 : /* last prime & all pr^f, pr | p, included. Include p^f instead */
1830 581 : if (oldf && i == nfa && degrel == nfa*f && !umodiu(discnf, p))
1831 0 : pr = utoipos(p);
1832 :
1833 : /* pr^f = N P, P | pr, hence is in norm group */
1834 581 : col = bnrisprincipalmod(bnr,pr,gdegrel,0);
1835 581 : if (f > 1) col = ZC_z_mul(col, f);
1836 581 : G = ZM_hnf(shallowconcat(G, col));
1837 581 : detG = ZM_det_triangular(G);
1838 581 : k = abscmpiu(detG,degrel);
1839 581 : if (k < 0) return NULL;
1840 581 : if (!k) { cgiv(detG); return G; }
1841 : }
1842 : }
1843 0 : return NULL;
1844 : }
1845 : GEN
1846 14 : rnfnormgroup(GEN bnr, GEN polrel)
1847 : {
1848 14 : pari_sp av = avma;
1849 14 : GEN G = rnfnormgroup_i(bnr, polrel);
1850 14 : if (!G) retgc_const(av, cgetg(1, t_MAT));
1851 7 : return gc_upto(av, G);
1852 : }
1853 :
1854 : GEN
1855 0 : nf_deg1_prime(GEN nf)
1856 : {
1857 0 : GEN z, T = nf_get_pol(nf), D = nf_get_disc(nf), f = nf_get_index(nf);
1858 0 : long degnf = degpol(T);
1859 : forprime_t S;
1860 : pari_sp av;
1861 : ulong p;
1862 0 : u_forprime_init(&S, degnf, ULONG_MAX);
1863 0 : av = avma;
1864 0 : while ( (p = u_forprime_next(&S)) )
1865 : {
1866 : ulong r;
1867 0 : if (!umodiu(D, p) || !umodiu(f, p)) continue;
1868 0 : r = Flx_oneroot(ZX_to_Flx(T,p), p);
1869 0 : if (r != p)
1870 : {
1871 0 : z = utoi(Fl_neg(r, p));
1872 0 : z = deg1pol_shallow(gen_1, z, varn(T));
1873 0 : return idealprimedec_kummer(nf, z, 1, utoipos(p));
1874 : }
1875 0 : set_avma(av);
1876 : }
1877 0 : return NULL;
1878 : }
1879 :
1880 : /* Given bnf and T defining an abelian relative extension, compute the
1881 : * corresponding conductor and congruence subgroup. Return
1882 : * [cond,bnr(cond),H] where cond=[ideal,arch] is the conductor. */
1883 : GEN
1884 826 : rnfconductor0(GEN bnf, GEN T, long flag)
1885 : {
1886 826 : pari_sp av = avma;
1887 : GEN P, E, D, nf, module, bnr, H, lim, Tr, MOD;
1888 : long i, l, degT;
1889 :
1890 826 : if (flag < 0 || flag > 2) pari_err_FLAG("rnfconductor");
1891 826 : bnf = checkbnf(bnf); nf = bnf_get_nf(bnf);
1892 812 : Tr = rnfdisc_get_T(nf, T, &lim);
1893 812 : T = nfX_to_monic(nf, Tr, NULL); degT = degpol(T);
1894 812 : if (!lim)
1895 791 : D = rnfdisc_factored(nf, T, NULL);
1896 : else
1897 : {
1898 21 : D = nfX_disc(nf, Q_primpart(Tr));
1899 21 : if (gequal0(D))
1900 0 : pari_err_DOMAIN("rnfconductor","issquarefree(pol)","=",gen_0, Tr);
1901 21 : D = idealfactor_partial(nf, D, lim);
1902 : }
1903 812 : P = gel(D,1); l = lg(P);
1904 812 : E = gel(D,2);
1905 1806 : for (i = 1; i < l; i++) /* cheaply update tame primes */
1906 : { /* v_pr(f) = 1 + \sum_{0 < i < l} g_i/g_0
1907 : <= 1 + max_{i>0} g_i/(g_i-1) \sum_{0 < i < l} g_i -1
1908 : <= 1 + (p/(p-1)) * v_P(e(L/K, pr)), P | pr | p */
1909 994 : GEN pr = gel(P,i), p = pr_get_p(pr), e = gen_1;
1910 994 : ulong q, e0 = itou(gel(E,i));
1911 994 : if (e0 > 1 && cmpiu(p, degT) <= 0)
1912 : {
1913 455 : long v, pp = itou(p);
1914 455 : if ((v = u_lvalrem(degT, pp, &q)))
1915 : { /* e = e_tame * e_wild, e_wild | p^v */
1916 357 : ulong t = ugcdiu(subiu(pr_norm(pr),1), q); /* e_tame | t */
1917 : /* upper bound for 1 + p/(p-1) * v * e(L/Q,p) */
1918 357 : e0 = minuu(e0, 1 + (pp * v * pr_get_e(pr) * upowuu(pp,v) * t) / (pp-1));
1919 357 : e = utoipos(e0);
1920 : }
1921 : }
1922 994 : gel(E,i) = e;
1923 : }
1924 812 : module = mkvec2(D, identity_perm(nf_get_r1(nf)));
1925 812 : MOD = flag? utoipos(degpol(T)): NULL;
1926 812 : bnr = Buchraymod_i(bnf, module, nf_INIT, MOD);
1927 812 : H = rnfnormgroup_i(bnr,T); if (!H) return gc_const(av,gen_0);
1928 1330 : return gc_GEN(av, flag == 2? bnrconductor_factored(bnr, H)
1929 532 : : bnrconductormod(bnr, H, MOD));
1930 : }
1931 : GEN
1932 35 : rnfconductor(GEN bnf, GEN T) { return rnfconductor0(bnf, T, 0); }
1933 :
1934 : static GEN
1935 1568 : prV_norms(GEN x)
1936 2877 : { pari_APPLY_same( pr_norm(gel(x,i))); }
1937 :
1938 : /* Given a number field bnf=bnr[1], a ray class group structure bnr, and a
1939 : * subgroup H (HNF form) of the ray class group, compute [n, r1, dk]
1940 : * attached to H. If flag & rnf_COND, abort (return NULL) if module is not the
1941 : * conductor. If flag & rnf_REL, return relative data, else absolute */
1942 : static GEN
1943 1645 : bnrdisc_i(GEN bnr, GEN H, long flag)
1944 : {
1945 1645 : const long flcond = flag & rnf_COND;
1946 : GEN nf, clhray, E, ED, dk;
1947 : long k, d, l, n, r1;
1948 : zlog_S S;
1949 :
1950 1645 : checkbnr(bnr);
1951 1645 : init_zlog(&S, bnr_get_bid(bnr));
1952 1645 : nf = bnr_get_nf(bnr);
1953 1645 : H = bnr_subgroup_check(bnr, H, &clhray);
1954 1645 : d = itos(clhray);
1955 1645 : if (!H) H = diagonal_shallow(bnr_get_cyc(bnr));
1956 1645 : E = S.k; ED = cgetg_copy(E, &l);
1957 2996 : for (k = 1; k < l; k++)
1958 : {
1959 1365 : long j, e = itos(gel(E,k)), eD = e*d;
1960 1365 : GEN H2 = H;
1961 1582 : for (j = e; j > 0; j--)
1962 : {
1963 1456 : GEN z = bnr_log_gen_pr(bnr, &S, j, k);
1964 : long d2;
1965 1456 : H2 = ZM_hnf(shallowconcat(H2, z));
1966 1456 : d2 = itos( ZM_det_triangular(H2) );
1967 1456 : if (flcond && j==e && d2 == d) return NULL;
1968 1442 : if (d2 == 1) { eD -= j; break; }
1969 217 : eD -= d2;
1970 : }
1971 1351 : gel(ED,k) = utoi(eD); /* v_{P[k]}(relative discriminant) */
1972 : }
1973 1631 : l = lg(S.archp); r1 = nf_get_r1(nf);
1974 1932 : for (k = 1; k < l; k++)
1975 : {
1976 329 : if (!contains(H, bnr_log_gen_arch(bnr, &S, k))) { r1--; continue; }
1977 98 : if (flcond) return NULL;
1978 : }
1979 : /* d = relative degree
1980 : * r1 = number of unramified real places;
1981 : * [P,ED] = factorization of relative discriminant */
1982 1603 : if (flag & rnf_REL)
1983 : {
1984 35 : n = d;
1985 35 : dk = factorbackprime(nf, S.P, ED);
1986 : }
1987 : else
1988 : {
1989 1568 : n = d * nf_get_degree(nf);
1990 1568 : r1= d * r1;
1991 1568 : dk = factorback2(prV_norms(S.P), ED);
1992 1568 : if (((n-r1)&3) == 2) dk = negi(dk); /* (2r2) mod 4 = 2: r2(relext) is odd */
1993 1568 : dk = mulii(dk, powiu(absi_shallow(nf_get_disc(nf)), d));
1994 : }
1995 1603 : return mkvec3(utoipos(n), utoi(r1), dk);
1996 : }
1997 : GEN
1998 1645 : bnrdisc(GEN bnr, GEN H, long flag)
1999 : {
2000 1645 : pari_sp av = avma;
2001 1645 : GEN D = bnrdisc_i(bnr, H, flag);
2002 1645 : return D? gc_GEN(av, D): gc_const(av, gen_0);
2003 : }
2004 : GEN
2005 175 : bnrdisc0(GEN A, GEN B, GEN C, long flag)
2006 : {
2007 175 : GEN H, bnr = ABC_to_bnr(A,B,C,&H, 0);
2008 175 : return bnrdisc(bnr,H,flag);
2009 : }
2010 :
2011 : /* Given a number field bnf=bnr[1], a ray class group structure bnr and a
2012 : * vector chi representing a character on the generators bnr[2][3], compute
2013 : * the conductor of chi. */
2014 : GEN
2015 7 : bnrconductorofchar(GEN bnr, GEN chi)
2016 : {
2017 7 : pari_sp av = avma;
2018 7 : return gc_GEN(av, bnrconductor_raw(bnr, chi));
2019 : }
2020 :
2021 : /* \sum U[i]*y[i], U[i],y[i] ZM, we allow lg(y) > lg(U). */
2022 : static GEN
2023 910 : ZMV_mul(GEN U, GEN y)
2024 : {
2025 910 : long i, l = lg(U);
2026 910 : GEN z = NULL;
2027 910 : if (l == 1) return cgetg(1,t_MAT);
2028 2324 : for (i = 1; i < l; i++)
2029 : {
2030 1442 : GEN u = ZM_mul(gel(U,i), gel(y,i));
2031 1442 : z = z? ZM_add(z, u): u;
2032 : }
2033 882 : return z;
2034 : }
2035 :
2036 : /* t = [bid,U], h = #Cl(K) */
2037 : static GEN
2038 910 : get_classno(GEN t, GEN h)
2039 : {
2040 910 : GEN bid = gel(t,1), m = gel(t,2), cyc = bid_get_cyc(bid), U = bid_get_U(bid);
2041 910 : return mulii(h, ZM_det_triangular(ZM_hnfmodid(ZMV_mul(U,m), cyc)));
2042 : }
2043 :
2044 : static void
2045 28 : chk_listBU(GEN L, const char *s) {
2046 28 : if (typ(L) != t_VEC) pari_err_TYPE(s,L);
2047 28 : if (lg(L) > 1) {
2048 28 : GEN z = gel(L,1);
2049 28 : if (typ(z) != t_VEC) pari_err_TYPE(s,z);
2050 28 : if (lg(z) == 1) return;
2051 28 : z = gel(z,1); /* [bid,U] */
2052 28 : if (typ(z) != t_VEC || lg(z) != 3) pari_err_TYPE(s,z);
2053 28 : checkbid(gel(z,1));
2054 : }
2055 : }
2056 :
2057 : /* Given lists of [bid, unit ideallogs], return lists of ray class numbers */
2058 : GEN
2059 7 : bnrclassnolist(GEN bnf,GEN L)
2060 : {
2061 7 : pari_sp av = avma;
2062 7 : long i, l = lg(L);
2063 : GEN V, h;
2064 :
2065 7 : chk_listBU(L, "bnrclassnolist");
2066 7 : if (l == 1) return cgetg(1, t_VEC);
2067 7 : bnf = checkbnf(bnf);
2068 7 : h = bnf_get_no(bnf);
2069 7 : V = cgetg(l,t_VEC);
2070 392 : for (i = 1; i < l; i++)
2071 : {
2072 385 : GEN v, z = gel(L,i);
2073 385 : long j, lz = lg(z);
2074 385 : gel(V,i) = v = cgetg(lz,t_VEC);
2075 826 : for (j=1; j<lz; j++) gel(v,j) = get_classno(gel(z,j), h);
2076 : }
2077 7 : return gc_GEN(av, V);
2078 : }
2079 :
2080 : static GEN
2081 1484 : Lbnrclassno(GEN L, GEN fac)
2082 : {
2083 1484 : long i, l = lg(L);
2084 2184 : for (i=1; i<l; i++)
2085 2184 : if (gequal(gmael(L,i,1),fac)) return gmael(L,i,2);
2086 0 : pari_err_BUG("Lbnrclassno");
2087 : return NULL; /* LCOV_EXCL_LINE */
2088 : }
2089 :
2090 : static GEN
2091 406 : factordivexact(GEN fa1,GEN fa2)
2092 : {
2093 : long i, j, k, c, l;
2094 : GEN P, E, P1, E1, P2, E2, p1;
2095 :
2096 406 : P1 = gel(fa1,1); E1 = gel(fa1,2); l = lg(P1);
2097 406 : P2 = gel(fa2,1); E2 = gel(fa2,2);
2098 406 : P = cgetg(l,t_COL);
2099 406 : E = cgetg(l,t_COL);
2100 903 : for (c = i = 1; i < l; i++)
2101 : {
2102 497 : j = RgV_isin(P2,gel(P1,i));
2103 497 : if (!j) { gel(P,c) = gel(P1,i); gel(E,c) = gel(E1,i); c++; }
2104 : else
2105 : {
2106 497 : p1 = subii(gel(E1,i), gel(E2,j)); k = signe(p1);
2107 497 : if (k < 0) pari_err_BUG("factordivexact [not exact]");
2108 497 : if (k > 0) { gel(P,c) = gel(P1,i); gel(E,c) = p1; c++; }
2109 : }
2110 : }
2111 406 : setlg(P, c);
2112 406 : setlg(E, c); return mkmat2(P, E);
2113 : }
2114 : /* remove index k */
2115 : static GEN
2116 1169 : factorsplice(GEN fa, long k)
2117 : {
2118 1169 : GEN p = gel(fa,1), e = gel(fa,2), P, E;
2119 1169 : long i, l = lg(p) - 1;
2120 1169 : P = cgetg(l, typ(p));
2121 1169 : E = cgetg(l, typ(e));
2122 1344 : for (i=1; i<k; i++) { P[i] = p[i]; E[i] = e[i]; }
2123 1169 : p++; e++;
2124 1694 : for ( ; i<l; i++) { P[i] = p[i]; E[i] = e[i]; }
2125 1169 : return mkvec2(P,E);
2126 : }
2127 : static GEN
2128 812 : factorpow(GEN fa, long n)
2129 : {
2130 812 : if (!n) return trivial_fact();
2131 812 : return mkmat2(gel(fa,1), gmulsg(n, gel(fa,2)));
2132 : }
2133 : static GEN
2134 1043 : factormul(GEN fa1,GEN fa2)
2135 : {
2136 1043 : GEN p, pnew, e, enew, v, P, y = famat_mul_shallow(fa1,fa2);
2137 : long i, c, lx;
2138 :
2139 1043 : p = gel(y,1); v = indexsort(p); lx = lg(p);
2140 1043 : e = gel(y,2);
2141 1043 : pnew = vecpermute(p, v);
2142 1043 : enew = vecpermute(e, v);
2143 1043 : P = gen_0; c = 0;
2144 2933 : for (i=1; i<lx; i++)
2145 : {
2146 1890 : if (gequal(gel(pnew,i),P))
2147 49 : gel(e,c) = addii(gel(e,c),gel(enew,i));
2148 : else
2149 : {
2150 1841 : c++; P = gel(pnew,i);
2151 1841 : gel(p,c) = P;
2152 1841 : gel(e,c) = gel(enew,i);
2153 : }
2154 : }
2155 1043 : setlg(p, c+1);
2156 1043 : setlg(e, c+1); return y;
2157 : }
2158 :
2159 : static long
2160 168 : get_nz(GEN bnf, GEN ideal, GEN arch, long clhray)
2161 : {
2162 : GEN arch2, mod;
2163 168 : long nz = 0, l = lg(arch), k, clhss;
2164 168 : if (typ(arch) == t_VECSMALL)
2165 14 : arch2 = indices_to_vec01(arch,nf_get_r1(bnf_get_nf(bnf)));
2166 : else
2167 154 : arch2 = leafcopy(arch);
2168 168 : mod = mkvec2(ideal, arch2);
2169 448 : for (k = 1; k < l; k++)
2170 : { /* FIXME: this is wasteful. Use the same algorithm as bnrconductor */
2171 301 : if (signe(gel(arch2,k)))
2172 : {
2173 28 : gel(arch2,k) = gen_0; clhss = itos(bnrclassno(bnf,mod));
2174 28 : gel(arch2,k) = gen_1;
2175 28 : if (clhss == clhray) return -1;
2176 : }
2177 273 : else nz++;
2178 : }
2179 147 : return nz;
2180 : }
2181 :
2182 : static GEN
2183 427 : get_NR1D(long Nf, long clhray, long degk, long nz, GEN fadkabs, GEN idealrel)
2184 : {
2185 : long n, R1;
2186 : GEN dlk;
2187 427 : if (nz < 0) return mkvec3(gen_0,gen_0,gen_0); /*EMPTY*/
2188 406 : n = clhray * degk;
2189 406 : R1 = clhray * nz;
2190 406 : dlk = factordivexact(factorpow(Z_factor(utoipos(Nf)),clhray), idealrel);
2191 : /* r2 odd, set dlk = -dlk */
2192 406 : if (((n-R1)&3)==2) dlk = factormul(to_famat_shallow(gen_m1,gen_1), dlk);
2193 406 : return mkvec3(utoipos(n),
2194 : stoi(R1),
2195 : factormul(dlk,factorpow(fadkabs,clhray)));
2196 : }
2197 :
2198 : /* t = [bid,U], h = #Cl(K) */
2199 : static GEN
2200 469 : get_discdata(GEN t, GEN h)
2201 : {
2202 469 : GEN bid = gel(t,1), fa = bid_get_fact(bid);
2203 469 : GEN P = gel(fa,1), E = vec_to_vecsmall(gel(fa,2));
2204 469 : return mkvec3(mkvec2(P, E), (GEN)itou(get_classno(t, h)), bid_get_mod(bid));
2205 : }
2206 : typedef struct _disc_data {
2207 : long degk;
2208 : GEN bnf, fadk, idealrelinit, V;
2209 : } disc_data;
2210 :
2211 : static GEN
2212 469 : get_discray(disc_data *D, GEN V, GEN z, long N)
2213 : {
2214 469 : GEN idealrel = D->idealrelinit;
2215 469 : GEN mod = gel(z,3), Fa = gel(z,1);
2216 469 : GEN P = gel(Fa,1), E = gel(Fa,2);
2217 469 : long k, nz, clhray = z[2], lP = lg(P);
2218 700 : for (k=1; k<lP; k++)
2219 : {
2220 546 : GEN pr = gel(P,k), p = pr_get_p(pr);
2221 546 : long e, ep = E[k], f = pr_get_f(pr);
2222 546 : long S = 0, norm = N, Npr = upowuu(p[2],f), clhss;
2223 798 : for (e=1; e<=ep; e++)
2224 : {
2225 : GEN fad;
2226 574 : if (e < ep) { E[k] = ep-e; fad = Fa; }
2227 462 : else fad = factorsplice(Fa, k);
2228 574 : norm /= Npr;
2229 574 : clhss = (long)Lbnrclassno(gel(V,norm), fad);
2230 574 : if (e==1 && clhss==clhray) { E[k] = ep; return cgetg(1, t_VEC); }
2231 259 : if (clhss == 1) { S += ep-e+1; break; }
2232 252 : S += clhss;
2233 : }
2234 231 : E[k] = ep;
2235 231 : idealrel = factormul(idealrel, to_famat_shallow(p, utoi(f * S)));
2236 : }
2237 154 : nz = get_nz(D->bnf, gel(mod,1), gel(mod,2), clhray);
2238 154 : return get_NR1D(N, clhray, D->degk, nz, D->fadk, idealrel);
2239 : }
2240 :
2241 : /* Given a list of bids and attached unit log matrices, return the
2242 : * list of discrayabs. Only keep moduli which are conductors. */
2243 : GEN
2244 21 : discrayabslist(GEN bnf, GEN L)
2245 : {
2246 21 : pari_sp av = avma;
2247 21 : long i, l = lg(L);
2248 : GEN nf, V, D, h;
2249 : disc_data ID;
2250 :
2251 21 : chk_listBU(L, "discrayabslist");
2252 21 : if (l == 1) return cgetg(1, t_VEC);
2253 21 : ID.bnf = bnf = checkbnf(bnf);
2254 21 : nf = bnf_get_nf(bnf);
2255 21 : h = bnf_get_no(bnf);
2256 21 : ID.degk = nf_get_degree(nf);
2257 21 : ID.fadk = absZ_factor(nf_get_disc(nf));
2258 21 : ID.idealrelinit = trivial_fact();
2259 21 : V = cgetg(l, t_VEC);
2260 21 : D = cgetg(l, t_VEC);
2261 448 : for (i = 1; i < l; i++)
2262 : {
2263 427 : GEN z = gel(L,i), v, d;
2264 427 : long j, lz = lg(z);
2265 427 : gel(V,i) = v = cgetg(lz,t_VEC);
2266 427 : gel(D,i) = d = cgetg(lz,t_VEC);
2267 896 : for (j=1; j<lz; j++) {
2268 469 : gel(d,j) = get_discdata(gel(z,j), h);
2269 469 : gel(v,j) = get_discray(&ID, D, gel(d,j), i);
2270 : }
2271 : }
2272 21 : return gc_GEN(av, V);
2273 : }
2274 :
2275 : /* a zsimp is [fa, cyc, v]
2276 : * fa: vecsmall factorisation,
2277 : * cyc: ZV (concatenation of (Z_K/pr^k)^* SNFs), the generators
2278 : * are positive at all real places [defined implicitly by weak approximation]
2279 : * v: ZC (log of units on (Z_K/pr^k)^* components) */
2280 : static GEN
2281 28 : zsimp(void)
2282 : {
2283 28 : GEN empty = cgetg(1, t_VECSMALL);
2284 28 : return mkvec3(mkvec2(empty,empty), cgetg(1,t_VEC), cgetg(1,t_MAT));
2285 : }
2286 :
2287 : /* fa a vecsmall factorization, append p^e */
2288 : static GEN
2289 175 : fasmall_append(GEN fa, long p, long e)
2290 : {
2291 175 : GEN P = gel(fa,1), E = gel(fa,2);
2292 175 : retmkvec2(vecsmall_append(P,p), vecsmall_append(E,e));
2293 : }
2294 :
2295 : /* sprk = sprkinit(pr,k), b zsimp with modulus coprime to pr */
2296 : static GEN
2297 518 : zsimpjoin(GEN b, GEN sprk, GEN U_pr, long prcode, long e)
2298 : {
2299 518 : GEN fa, cyc = sprk_get_cyc(sprk);
2300 518 : if (lg(gel(b,2)) == 1) /* trivial group */
2301 343 : fa = mkvec2(mkvecsmall(prcode),mkvecsmall(e));
2302 : else
2303 : {
2304 175 : fa = fasmall_append(gel(b,1), prcode, e);
2305 175 : cyc = shallowconcat(gel(b,2), cyc); /* no SNF ! */
2306 175 : U_pr = vconcat(gel(b,3),U_pr);
2307 : }
2308 518 : return mkvec3(fa, cyc, U_pr);
2309 : }
2310 : /* B a zsimp, sgnU = [cyc[f_oo], sgn_{f_oo}(units)] */
2311 : static GEN
2312 28 : bnrclassno_1(GEN B, ulong h, GEN sgnU)
2313 : {
2314 28 : long lx = lg(B), j;
2315 28 : GEN L = cgetg(lx,t_VEC);
2316 56 : for (j=1; j<lx; j++)
2317 : {
2318 28 : pari_sp av = avma;
2319 28 : GEN b = gel(B,j), cyc = gel(b,2), qm = gel(b,3);
2320 : ulong z;
2321 28 : cyc = shallowconcat(cyc, gel(sgnU,1));
2322 28 : qm = vconcat(qm, gel(sgnU,2));
2323 28 : z = itou( mului(h, ZM_det_triangular(ZM_hnfmodid(qm, cyc))) );
2324 28 : set_avma(av);
2325 28 : gel(L,j) = mkvec2(gel(b,1), mkvecsmall(z));
2326 : }
2327 28 : return L;
2328 : }
2329 :
2330 : static void
2331 1344 : vecselect_p(GEN A, GEN B, GEN p, long init, long lB)
2332 : {
2333 1344 : long i; setlg(B, lB);
2334 2688 : for (i=init; i<lB; i++) B[i] = A[p[i]];
2335 1344 : }
2336 : /* B := p . A = row selection according to permutation p. Treat only lower
2337 : * right corner init x init */
2338 : static void
2339 1022 : rowselect_p(GEN A, GEN B, GEN p, long init)
2340 : {
2341 1022 : long i, lB = lg(A), lp = lg(p);
2342 2436 : for (i=1; i<init; i++) setlg(B[i],lp);
2343 2366 : for ( ; i<lB; i++) vecselect_p(gel(A,i),gel(B,i),p,init,lp);
2344 1022 : }
2345 : static ulong
2346 1022 : hdet(ulong h, GEN m)
2347 : {
2348 1022 : pari_sp av = avma;
2349 1022 : GEN z = mului(h, ZM_det_triangular(ZM_hnf(m)));
2350 1022 : return gc_ulong(av, itou(z));
2351 : }
2352 : static GEN
2353 1106 : bnrclassno_all(GEN B, ulong h, GEN sgnU)
2354 : {
2355 : long lx, k, kk, j, r1, jj, nba, nbarch;
2356 : GEN _2, L, m, H, mm, rowsel;
2357 :
2358 1106 : if (typ(sgnU) == t_VEC) return bnrclassno_1(B,h,sgnU);
2359 1078 : lx = lg(B); if (lx == 1) return B;
2360 :
2361 371 : r1 = nbrows(sgnU); _2 = const_vec(r1, gen_2);
2362 371 : L = cgetg(lx,t_VEC); nbarch = 1L<<r1;
2363 889 : for (j=1; j<lx; j++)
2364 : {
2365 518 : pari_sp av = avma;
2366 518 : GEN b = gel(B,j), cyc = gel(b,2), qm = gel(b,3);
2367 518 : long nc = lg(cyc)-1;
2368 : /* [ qm cyc 0 ]
2369 : * [ sgnU 0 2 ] */
2370 518 : m = ZM_hnfmodid(vconcat(qm, sgnU), shallowconcat(cyc,_2));
2371 518 : mm = RgM_shallowcopy(m);
2372 518 : rowsel = cgetg(nc+r1+1,t_VECSMALL);
2373 518 : H = cgetg(nbarch+1,t_VECSMALL);
2374 1540 : for (k = 0; k < nbarch; k++)
2375 : {
2376 1022 : nba = nc+1;
2377 2366 : for (kk=k,jj=1; jj<=r1; jj++,kk>>=1)
2378 1344 : if (kk&1) rowsel[nba++] = nc + jj;
2379 1022 : setlg(rowsel, nba);
2380 1022 : rowselect_p(m, mm, rowsel, nc+1);
2381 1022 : H[k+1] = hdet(h, mm);
2382 : }
2383 518 : H = gc_leaf(av, H);
2384 518 : gel(L,j) = mkvec2(gel(b,1), H);
2385 : }
2386 371 : return L;
2387 : }
2388 :
2389 : static int
2390 21 : is_module(GEN v)
2391 : {
2392 21 : if (lg(v) != 3 || (typ(v) != t_MAT && typ(v) != t_VEC)) return 0;
2393 21 : return typ(gel(v,1)) == t_VECSMALL && typ(gel(v,2)) == t_VECSMALL;
2394 : }
2395 : GEN
2396 21 : decodemodule(GEN nf, GEN fa)
2397 : {
2398 : long n, nn, k;
2399 21 : pari_sp av = avma;
2400 : GEN G, E, id, pr;
2401 :
2402 21 : nf = checknf(nf);
2403 21 : if (!is_module(fa)) pari_err_TYPE("decodemodule [not a factorization]", fa);
2404 21 : n = nf_get_degree(nf); nn = n*n; id = NULL;
2405 21 : G = gel(fa,1);
2406 21 : E = gel(fa,2);
2407 35 : for (k=1; k<lg(G); k++)
2408 : {
2409 14 : long code = G[k], p = code / nn, j = (code%n)+1;
2410 14 : GEN P = idealprimedec(nf, utoipos(p)), e = stoi(E[k]);
2411 14 : if (lg(P) <= j) pari_err_BUG("decodemodule [incorrect hash code]");
2412 14 : pr = gel(P,j);
2413 14 : id = id? idealmulpowprime(nf,id, pr,e)
2414 14 : : idealpow(nf, pr,e);
2415 : }
2416 21 : if (!id) { set_avma(av); return matid(n); }
2417 14 : return gc_upto(av,id);
2418 : }
2419 :
2420 : /* List of ray class fields. Do all from scratch, bound < 2^30. No subgroups.
2421 : *
2422 : * Output: a vector V, V[k] contains the ideals of norm k. Given such an ideal
2423 : * m, the component is as follows:
2424 : *
2425 : * + if arch = NULL, run through all possible archimedean parts; archs are
2426 : * ordered using inverse lexicographic order, [0,..,0], [1,0,..,0], [0,1,..,0],
2427 : * Component is [m,V] where V is a vector with 2^r1 entries, giving for each
2428 : * arch the triple [N,R1,D], with N, R1, D as in discrayabs; D is in factored
2429 : * form.
2430 : *
2431 : * + otherwise [m,N,R1,D] */
2432 : GEN
2433 28 : discrayabslistarch(GEN bnf, GEN arch, ulong bound)
2434 : {
2435 28 : int allarch = (arch==NULL), flbou = 0;
2436 : long degk, j, k, l, nba, nbarch, r1, c, sqbou;
2437 28 : pari_sp av0 = avma, av, av1;
2438 : GEN nf, p, Z, fa, Disc, U, sgnU, EMPTY, empty, archp;
2439 : GEN res, Ray, discall, idealrel, idealrelinit, fadkabs, BOUND;
2440 : ulong i, h;
2441 : forprime_t S;
2442 :
2443 28 : if (bound == 0)
2444 0 : pari_err_DOMAIN("discrayabslistarch","bound","==",gen_0,utoi(bound));
2445 28 : res = discall = NULL; /* -Wall */
2446 :
2447 28 : bnf = checkbnf(bnf);
2448 28 : nf = bnf_get_nf(bnf);
2449 28 : r1 = nf_get_r1(nf);
2450 28 : degk = nf_get_degree(nf);
2451 28 : fadkabs = absZ_factor(nf_get_disc(nf));
2452 28 : h = itou(bnf_get_no(bnf));
2453 :
2454 28 : if (allarch)
2455 : {
2456 21 : if (r1>15) pari_err_IMPL("r1>15 in discrayabslistarch");
2457 21 : arch = const_vec(r1, gen_1);
2458 : }
2459 7 : else if (lg(arch)-1 != r1)
2460 0 : pari_err_TYPE("Idealstar [incorrect archimedean component]",arch);
2461 28 : U = log_prk_units_init(bnf);
2462 28 : archp = vec01_to_indices(arch);
2463 28 : nba = lg(archp)-1;
2464 28 : sgnU = zm_to_ZM( nfsign_units(bnf, archp, 1) );
2465 28 : if (!allarch) sgnU = mkvec2(const_vec(nba,gen_2), sgnU);
2466 :
2467 28 : empty = cgetg(1,t_VEC);
2468 : /* what follows was rewritten from Ideallist */
2469 28 : BOUND = utoipos(bound);
2470 28 : p = cgetipos(3);
2471 28 : u_forprime_init(&S, 2, bound);
2472 28 : av = avma;
2473 28 : sqbou = (long)sqrt((double)bound) + 1;
2474 28 : Z = const_vec(bound, empty);
2475 28 : gel(Z,1) = mkvec(zsimp());
2476 28 : if (DEBUGLEVEL>1) err_printf("Starting zidealstarunits computations\n");
2477 : /* The goal is to compute Ray (lists of bnrclassno). Z contains "zsimps",
2478 : * simplified bid, from which bnrclassno is easy to compute.
2479 : * Once p > sqbou, delete Z[i] for i > sqbou and compute directly Ray */
2480 28 : Ray = Z;
2481 294 : while ((p[2] = u_forprime_next(&S)))
2482 : {
2483 266 : if (!flbou && p[2] > sqbou)
2484 : {
2485 21 : flbou = 1;
2486 21 : if (DEBUGLEVEL>1) err_printf("\nStarting bnrclassno computations\n");
2487 21 : Z = gc_GEN(av,Z);
2488 21 : Ray = cgetg(bound+1, t_VEC);
2489 889 : for (i=1; i<=bound; i++) gel(Ray,i) = bnrclassno_all(gel(Z,i),h,sgnU);
2490 21 : Z = vecslice(Z, 1, sqbou);
2491 : }
2492 266 : fa = idealprimedec_limit_norm(nf,p,BOUND);
2493 504 : for (j=1; j<lg(fa); j++)
2494 : {
2495 238 : GEN pr = gel(fa,j);
2496 238 : long prcode, f = pr_get_f(pr);
2497 238 : ulong q, Q = upowuu(p[2], f);
2498 :
2499 : /* p, f-1, j-1 as a single integer in "base degk" (f,j <= degk)*/
2500 238 : prcode = (p[2]*degk + f-1)*degk + j-1;
2501 238 : q = Q;
2502 : /* FIXME: if Q = 2, should start at l = 2 */
2503 238 : for (l = 1;; l++) /* Q <= bound */
2504 105 : {
2505 : ulong iQ;
2506 343 : GEN sprk = log_prk_init(nf, pr, l, NULL);
2507 343 : GEN U_pr = log_prk_units(nf, U, sprk);
2508 1582 : for (iQ = Q, i = 1; iQ <= bound; iQ += Q, i++)
2509 : {
2510 1239 : GEN pz, p2, p1 = gel(Z,i);
2511 1239 : long lz = lg(p1);
2512 1239 : if (lz == 1) continue;
2513 :
2514 595 : p2 = cgetg(lz,t_VEC); c = 0;
2515 1113 : for (k=1; k<lz; k++)
2516 : {
2517 658 : GEN z = gel(p1,k), v = gmael(z,1,1); /* primes in zsimp's fact. */
2518 658 : long lv = lg(v);
2519 : /* If z has a power of pr in its modulus, skip it */
2520 658 : if (i != 1 && lv > 1 && v[lv-1] == prcode) break;
2521 518 : gel(p2,++c) = zsimpjoin(z,sprk,U_pr,prcode,l);
2522 : }
2523 595 : setlg(p2, c+1);
2524 595 : pz = gel(Ray,iQ);
2525 595 : if (flbou) p2 = bnrclassno_all(p2,h,sgnU);
2526 595 : if (lg(pz) > 1) p2 = shallowconcat(pz,p2);
2527 595 : gel(Ray,iQ) = p2;
2528 : }
2529 343 : Q = itou_or_0( muluu(Q, q) );
2530 343 : if (!Q || Q > bound) break;
2531 : }
2532 : }
2533 266 : if (gc_needed(av,1))
2534 : {
2535 0 : if(DEBUGMEM>1) pari_warn(warnmem,"[1]: discrayabslistarch");
2536 0 : (void)gc_all(av, flbou? 2: 1, &Z, &Ray);
2537 : }
2538 : }
2539 28 : if (!flbou) /* occurs iff bound = 1,2,4 */
2540 : {
2541 7 : if (DEBUGLEVEL>1) err_printf("\nStarting bnrclassno computations\n");
2542 7 : Ray = cgetg(bound+1, t_VEC);
2543 35 : for (i=1; i<=bound; i++) gel(Ray,i) = bnrclassno_all(gel(Z,i),h,sgnU);
2544 : }
2545 28 : Ray = gc_GEN(av, Ray);
2546 :
2547 28 : if (DEBUGLEVEL>1) err_printf("Starting discrayabs computations\n");
2548 28 : if (allarch) nbarch = 1L<<r1;
2549 : else
2550 : {
2551 7 : nbarch = 1;
2552 7 : discall = cgetg(2,t_VEC);
2553 : }
2554 28 : EMPTY = mkvec3(gen_0,gen_0,gen_0);
2555 28 : idealrelinit = trivial_fact();
2556 28 : av1 = avma;
2557 28 : Disc = const_vec(bound, empty);
2558 924 : for (i=1; i<=bound; i++)
2559 : {
2560 896 : GEN sousdisc, sous = gel(Ray,i);
2561 896 : long ls = lg(sous);
2562 896 : gel(Disc,i) = sousdisc = cgetg(ls,t_VEC);
2563 1442 : for (j=1; j<ls; j++)
2564 : {
2565 546 : GEN b = gel(sous,j), clhrayall = gel(b,2), Fa = gel(b,1);
2566 546 : GEN P = gel(Fa,1), E = gel(Fa,2);
2567 546 : long lP = lg(P), karch;
2568 :
2569 546 : if (allarch) discall = cgetg(nbarch+1,t_VEC);
2570 1596 : for (karch=0; karch<nbarch; karch++)
2571 : {
2572 1050 : long nz, clhray = clhrayall[karch+1];
2573 1050 : if (allarch)
2574 : {
2575 : long ka, k2;
2576 1022 : nba = 0;
2577 2366 : for (ka=karch,k=1; k<=r1; k++,ka>>=1)
2578 1344 : if (ka & 1) nba++;
2579 1918 : for (k2=1,k=1; k<=r1; k++,k2<<=1)
2580 1190 : if (karch&k2 && clhrayall[karch-k2+1] == clhray)
2581 294 : { res = EMPTY; goto STORE; }
2582 : }
2583 756 : idealrel = idealrelinit;
2584 1078 : for (k=1; k<lP; k++) /* cf get_discray */
2585 : {
2586 805 : long e, ep = E[k], pf = P[k] / degk, f = (pf%degk) + 1, S = 0;
2587 805 : ulong normi = i, Npr;
2588 805 : p = utoipos(pf / degk);
2589 805 : Npr = upowuu(p[2],f);
2590 1204 : for (e=1; e<=ep; e++)
2591 : {
2592 : long clhss;
2593 : GEN fad;
2594 910 : if (e < ep) { E[k] = ep-e; fad = Fa; }
2595 707 : else fad = factorsplice(Fa, k);
2596 910 : normi /= Npr;
2597 910 : clhss = Lbnrclassno(gel(Ray,normi),fad)[karch+1];
2598 910 : if (e==1 && clhss==clhray) { E[k] = ep; res = EMPTY; goto STORE; }
2599 427 : if (clhss == 1) { S += ep-e+1; break; }
2600 399 : S += clhss;
2601 : }
2602 322 : E[k] = ep;
2603 322 : idealrel = factormul(idealrel, to_famat_shallow(p, utoi(f * S)));
2604 : }
2605 273 : if (!allarch && nba)
2606 14 : nz = get_nz(bnf, decodemodule(nf,Fa), arch, clhray);
2607 : else
2608 259 : nz = r1 - nba;
2609 273 : res = get_NR1D(i, clhray, degk, nz, fadkabs, idealrel);
2610 1050 : STORE: gel(discall,karch+1) = res;
2611 : }
2612 518 : res = allarch? mkvec2(Fa, discall)
2613 546 : : mkvec4(Fa, gel(res,1), gel(res,2), gel(res,3));
2614 546 : gel(sousdisc,j) = res;
2615 546 : if (gc_needed(av1,1))
2616 : {
2617 : long jj;
2618 0 : if(DEBUGMEM>1) pari_warn(warnmem,"[2]: discrayabslistarch");
2619 0 : for (jj=j+1; jj<ls; jj++) gel(sousdisc,jj) = gen_0; /* dummy */
2620 0 : Disc = gc_GEN(av1, Disc);
2621 0 : sousdisc = gel(Disc,i);
2622 : }
2623 : }
2624 : }
2625 28 : return gc_GEN(av0, Disc);
2626 : }
2627 :
2628 : int
2629 95214 : subgroup_conductor_ok(GEN H, GEN L)
2630 : { /* test conductor */
2631 95214 : long i, l = lg(L);
2632 222388 : for (i = 1; i < l; i++)
2633 180681 : if ( hnf_solve(H, gel(L,i)) ) return 0;
2634 41707 : return 1;
2635 : }
2636 : static GEN
2637 183014 : conductor_elts(GEN bnr)
2638 : {
2639 : long le, la, i, k;
2640 : GEN e, L;
2641 : zlog_S S;
2642 :
2643 183014 : if (!bnrisconductor(bnr, NULL)) return NULL;
2644 56317 : init_zlog(&S, bnr_get_bid(bnr));
2645 56326 : e = S.k; le = lg(e); la = lg(S.archp);
2646 56326 : L = cgetg(le + la - 1, t_VEC);
2647 56327 : i = 1;
2648 129584 : for (k = 1; k < le; k++)
2649 73258 : gel(L,i++) = bnr_log_gen_pr(bnr, &S, itos(gel(e,k)), k);
2650 78955 : for (k = 1; k < la; k++)
2651 22631 : gel(L,i++) = bnr_log_gen_arch(bnr, &S, k);
2652 56324 : return L;
2653 : }
2654 :
2655 : /* Let C a congruence group in bnr, compute its subgroups whose index is
2656 : * described by bound (see subgrouplist) as subgroups of Clk(bnr).
2657 : * Restrict to subgroups having the same conductor as bnr */
2658 : GEN
2659 455 : subgrouplist_cond_sub(GEN bnr, GEN C, GEN bound)
2660 : {
2661 455 : pari_sp av = avma;
2662 : long l, i, j;
2663 455 : GEN D, Mr, U, T, subgrp, L, cyc = bnr_get_cyc(bnr);
2664 :
2665 455 : L = conductor_elts(bnr); if (!L) return cgetg(1,t_VEC);
2666 455 : Mr = diagonal_shallow(cyc);
2667 455 : D = ZM_snfall_i(hnf_solve(C, Mr), &U, NULL, 1);
2668 455 : T = ZM_mul(C, RgM_inv(U));
2669 455 : subgrp = subgrouplist(D, bound);
2670 455 : l = lg(subgrp);
2671 966 : for (i = j = 1; i < l; i++)
2672 : {
2673 511 : GEN H = ZM_hnfmodid(ZM_mul(T, gel(subgrp,i)), cyc);
2674 511 : if (subgroup_conductor_ok(H, L)) gel(subgrp, j++) = H;
2675 : }
2676 455 : setlg(subgrp, j);
2677 455 : return gc_GEN(av, subgrp);
2678 : }
2679 :
2680 : static GEN
2681 182559 : subgroupcond(GEN bnr, GEN indexbound)
2682 : {
2683 182559 : pari_sp av = avma;
2684 182559 : GEN L = conductor_elts(bnr);
2685 :
2686 182551 : if (!L) return cgetg(1, t_VEC);
2687 55863 : L = subgroupcondlist(bnr_get_cyc(bnr), indexbound, L);
2688 55874 : if (indexbound && typ(indexbound) != t_VEC)
2689 : { /* sort by increasing index if not single value */
2690 14 : long i, l = lg(L);
2691 14 : GEN D = cgetg(l,t_VEC);
2692 245 : for (i=1; i<l; i++) gel(D,i) = ZM_det_triangular(gel(L,i));
2693 14 : L = vecreverse( vecpermute(L, indexsort(D)) );
2694 : }
2695 55874 : return gc_GEN(av, L);
2696 : }
2697 :
2698 : GEN
2699 184778 : subgrouplist0(GEN cyc, GEN indexbound, long all)
2700 : {
2701 184778 : if (!all && checkbnr_i(cyc)) return subgroupcond(cyc,indexbound);
2702 2212 : if (typ(cyc) != t_VEC || !RgV_is_ZV(cyc)) cyc = member_cyc(cyc);
2703 2205 : return subgrouplist(cyc,indexbound);
2704 : }
2705 :
2706 : GEN
2707 49 : bnrdisclist0(GEN bnf, GEN L, GEN arch)
2708 : {
2709 49 : if (typ(L)!=t_INT) return discrayabslist(bnf,L);
2710 28 : return discrayabslistarch(bnf,arch,itos(L));
2711 : }
2712 :
2713 : /****************************************************************************/
2714 : /* Galois action on a BNR */
2715 : /****************************************************************************/
2716 : GEN
2717 3094 : bnrautmatrix(GEN bnr, GEN aut)
2718 : {
2719 3094 : pari_sp av = avma;
2720 3094 : GEN bnf = bnr_get_bnf(bnr), nf = bnf_get_nf(bnf), bid = bnr_get_bid(bnr);
2721 3094 : GEN M, Gen = get_Gen(bnf, bid, bnr_get_El(bnr)), cyc = bnr_get_cyc(bnr);
2722 3094 : long i, l = lg(Gen);
2723 :
2724 3094 : M = cgetg(l, t_MAT); aut = nfgaloismatrix(nf, aut);
2725 : /* Gen = clg.gen*U, clg.gen = Gen*Ui */
2726 11676 : for (i = 1; i < l; i++)
2727 8582 : gel(M,i) = isprincipalray(bnr, nfgaloismatrixapply(nf, aut, gel(Gen,i)));
2728 3094 : M = ZM_mul(M, bnr_get_Ui(bnr));
2729 3094 : return gc_GEN(av, ZM_ZV_mod(M, cyc));
2730 : }
2731 :
2732 : GEN
2733 231 : bnrgaloismatrix(GEN bnr, GEN aut)
2734 : {
2735 231 : checkbnr(bnr);
2736 231 : switch (typ(aut))
2737 : {
2738 0 : case t_POL:
2739 : case t_COL:
2740 0 : return bnrautmatrix(bnr, aut);
2741 231 : case t_VEC:
2742 : {
2743 231 : pari_sp av = avma;
2744 231 : long i, l = lg(aut);
2745 : GEN v;
2746 231 : if (l == 9)
2747 : {
2748 7 : GEN g = gal_get_gen(aut);
2749 7 : if (typ(g) == t_VEC) { aut = galoispermtopol(aut, g); l = lg(aut); }
2750 : }
2751 231 : v = cgetg(l, t_VEC);
2752 693 : for(i = 1; i < l; i++) gel(v,i) = bnrautmatrix(bnr, gel(aut,i));
2753 231 : return gc_upto(av, v);
2754 : }
2755 0 : default:
2756 0 : pari_err_TYPE("bnrgaloismatrix", aut);
2757 : return NULL; /*LCOV_EXCL_LINE*/
2758 : }
2759 : }
2760 :
2761 : GEN
2762 3577 : bnrgaloisapply(GEN bnr, GEN mat, GEN x)
2763 : {
2764 3577 : pari_sp av=avma;
2765 : GEN cyc;
2766 3577 : checkbnr(bnr);
2767 3577 : cyc = bnr_get_cyc(bnr);
2768 3577 : if (typ(mat)!=t_MAT || !RgM_is_ZM(mat))
2769 0 : pari_err_TYPE("bnrgaloisapply",mat);
2770 3577 : if (typ(x)!=t_MAT || !RgM_is_ZM(x))
2771 0 : pari_err_TYPE("bnrgaloisapply",x);
2772 3577 : return gc_upto(av, ZM_hnfmodid(ZM_mul(mat, x), cyc));
2773 : }
2774 :
2775 : static GEN
2776 448 : check_bnrgal(GEN bnr, GEN M)
2777 : {
2778 448 : checkbnr(bnr);
2779 448 : if (typ(M)==t_MAT)
2780 0 : return mkvec(M);
2781 448 : else if (typ(M)==t_VEC && lg(M)==9 && typ(gal_get_gen(M))==t_VEC)
2782 : {
2783 224 : pari_sp av = avma;
2784 224 : GEN V = galoispermtopol(M, gal_get_gen(M));
2785 224 : return gc_upto(av, bnrgaloismatrix(bnr, V));
2786 : }
2787 224 : else if (!is_vec_t(typ(M)))
2788 0 : pari_err_TYPE("bnrisgalois",M);
2789 224 : return M;
2790 : }
2791 :
2792 : long
2793 448 : bnrisgalois(GEN bnr, GEN M, GEN H)
2794 : {
2795 448 : pari_sp av = avma;
2796 : long i, l;
2797 448 : if (typ(H)!=t_MAT || !RgM_is_ZM(H))
2798 0 : pari_err_TYPE("bnrisgalois",H);
2799 448 : M = check_bnrgal(bnr, M); l = lg(M);
2800 616 : for (i=1; i<l; i++)
2801 : {
2802 560 : long res = ZM_equal(bnrgaloisapply(bnr,gel(M,i), H), H);
2803 560 : if (!res) return gc_long(av,0);
2804 : }
2805 56 : return gc_long(av,1);
2806 : }
2807 :
2808 : static GEN
2809 14 : bnrlcmcond(GEN bnr1, GEN bnr2)
2810 : {
2811 14 : GEN I1 = bnr_get_bid(bnr1), f1 = bid_get_fact(I1), a1 = bid_get_arch(I1);
2812 14 : GEN I2 = bnr_get_bid(bnr2), f2 = bid_get_fact(I2), a2 = bid_get_arch(I2);
2813 : GEN f, a;
2814 : long i, l;
2815 14 : if (!gidentical(bnr_get_nf(bnr1), bnr_get_nf(bnr2)))
2816 0 : pari_err_TYPE("bnrcompositum[different fields]", mkvec2(bnr1,bnr2));
2817 14 : f = merge_factor(f1, f2, (void*)&cmp_prime_ideal, &cmp_nodata);
2818 14 : a = cgetg_copy(a1, &l);
2819 28 : for (i = 1; i < l; i++)
2820 14 : gel(a,i) = (signe(gel(a1,i)) || signe(gel(a2,i)))? gen_1: gen_0;
2821 14 : return mkvec2(f, a);
2822 : }
2823 : /* H subgroup (of bnr.clgp) in HNF; lift to BNR */
2824 : static GEN
2825 28 : bnrliftsubgroup(GEN BNR, GEN bnr, GEN H)
2826 : {
2827 28 : GEN E = gel(bnrsurjection(BNR, bnr), 1), K = kerint(shallowconcat(E, H));
2828 28 : return ZM_hnfmodid(rowslice(K, 1, lg(E)-1), bnr_get_cyc(BNR));
2829 : }
2830 : GEN
2831 14 : bnrcompositum(GEN fH1, GEN fH2)
2832 : {
2833 14 : pari_sp av = avma;
2834 : GEN bnr1, bnr2, bnr, H1, H2, H, n1, n2, d;
2835 14 : if (typ(fH1) != t_VEC || lg(fH2) != 3) pari_err_TYPE("bnrcompositum", fH1);
2836 14 : if (typ(fH2) != t_VEC || lg(fH2) != 3) pari_err_TYPE("bnrcompositum", fH2);
2837 14 : bnr1 = gel(fH1,1); if (!checkbnr_i(bnr1)) pari_err_TYPE("bnrcompositum",bnr1);
2838 14 : bnr2 = gel(fH2,1); if (!checkbnr_i(bnr2)) pari_err_TYPE("bnrcompositum",bnr2);
2839 14 : H1 = bnr_subgroup_check(bnr1, gel(fH1,2), &n1);
2840 14 : if (!H1) H1 = diagonal_shallow(bnr_get_cyc(bnr1));
2841 14 : H2 = bnr_subgroup_check(bnr2, gel(fH2,2), &n2);
2842 14 : if (!H2) H2 = diagonal_shallow(bnr_get_cyc(bnr2));
2843 14 : bnr = bnrinitmod(bnr_get_bnf(bnr1), bnrlcmcond(bnr1, bnr2), 0, lcmii(n1,n2));
2844 14 : d = bnr_get_cyc(bnr);
2845 14 : H1 = bnrliftsubgroup(bnr, bnr1, H1);
2846 14 : H2 = bnrliftsubgroup(bnr, bnr2, H2);
2847 14 : H = ZM_hnfintersectmod(H1, H2, gel(d,1));
2848 14 : return gc_GEN(av, mkvec2(bnr, ZM_hnfmodid(H, d)));
2849 : }
|
