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