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