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