1 /* $Id: buch3.c 7622 2006-01-23 18:46:57Z kb $
2
3 Copyright (C) 2000 The PARI group.
4
5 This file is part of the PARI/GP package.
6
7 PARI/GP is free software; you can redistribute it and/or modify it under the
8 terms of the GNU General Public License as published by the Free Software
9 Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10 ANY WARRANTY WHATSOEVER.
11
12 Check the License for details. You should have received a copy of it, along
13 with the package; see the file 'COPYING'. If not, write to the Free Software
14 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15
16 /*******************************************************************/
17 /* */
18 /* RAY CLASS FIELDS */
19 /* */
20 /*******************************************************************/
21 #include "pari.h"
22 #include "paripriv.h"
23
24 /* Faster than Buchray (because it can use zsignunits: easier zarchstar) */
25 GEN
buchnarrow(GEN bnf)26 buchnarrow(GEN bnf)
27 {
28 GEN nf, cyc, gen, GD, v, invpi, logs, p1, p2, R, basecl, met, u1, archp, clgp;
29 long r1, i, j, ngen, t, lo, c;
30 pari_sp av = avma;
31
32 bnf = checkbnf(bnf);
33 nf = checknf(bnf); r1 = nf_get_r1(nf);
34 clgp = gmael(bnf,8,1);
35 if (!r1) return gcopy(clgp);
36
37 cyc = gel(clgp,2);
38 gen = gel(clgp,3);
39 v = FpM_image(zsignunits(bnf, NULL, 1), gen_2);
40 t = lg(v)-1;
41 if (t == r1) { avma = av; return gcopy(clgp); }
42
43 ngen = lg(gen)-1;
44 p1 = cgetg(ngen+r1-t + 1,t_COL);
45 for (i=1; i<=ngen; i++) p1[i] = gen[i];
46 gen = p1;
47 v = archstar_full_rk(NULL, gmael(nf,5,1), ZM_to_Flm(v, 2), gen + (ngen - t));
48 v = rowslice(v, t+1, r1);
49
50 logs = cgetg(ngen+1,t_MAT);
51 GD = gmael(bnf,9,3); invpi = ginv( mppi(DEFAULTPREC) );
52 archp = perm_identity(r1);
53 for (j=1; j<=ngen; j++)
54 {
55 GEN z = zsign_from_logarch(gel(GD,j), invpi, archp);
56 gel(logs,j) = F2V_red_ip( gmul(v, z) );
57 }
58 /* [ cyc 0 ]
59 * [ logs 2 ] = relation matrix for Cl_f */
60 R = shallowconcat(
61 vconcat(diagonal_i(cyc), logs),
62 vconcat(zeromat(ngen, r1-t), gscalmat(gen_2,r1-t))
63 );
64
65 met = smithrel(R,NULL,&u1);
66 lo = lg(R); c = lg(met);
67 if (DEBUGLEVEL>3) msgtimer("smith/class group");
68
69 basecl = cgetg(c,t_VEC);
70 for (j=1; j<c; j++)
71 {
72 p1 = gcoeff(u1,1,j);
73 p2 = idealpow(nf,gel(gen,1),p1);
74 if (signe(p1) < 0) p2 = Q_primpart(p2);
75 for (i=2; i<lo; i++)
76 {
77 p1 = gcoeff(u1,i,j);
78 if (signe(p1))
79 {
80 p2 = idealmul(nf,p2, idealpow(nf,gel(gen,i),p1));
81 p2 = Q_primpart(p2);
82 }
83 }
84 gel(basecl,j) = p2;
85 }
86 return gerepilecopy(av, mkvec3(shifti(gel(clgp,1), r1-t), met,basecl));
87 }
88
89 /********************************************************************/
90 /** **/
91 /** REDUCTION MOD IDELE **/
92 /** **/
93 /********************************************************************/
94
95 static GEN
compute_fact(GEN nf,GEN u1,GEN gen)96 compute_fact(GEN nf, GEN u1, GEN gen)
97 {
98 GEN G, basecl;
99 long prec,i,j, l = lg(u1), h = lg(u1[1]); /* l > 1 */
100
101 basecl = cgetg(l,t_VEC);
102 prec = nfgetprec(nf);
103 G = cgetg(3,t_VEC);
104 gel(G,2) = cgetg(1,t_MAT);
105
106 for (j=1; j<l; j++)
107 {
108 GEN g,e, z = NULL;
109 for (i=1; i<h; i++)
110 {
111 e = gcoeff(u1,i,j); if (!signe(e)) continue;
112
113 g = gel(gen,i);
114 if (typ(g) != t_MAT)
115 {
116 if (z)
117 gel(z,2) = arch_mul(gel(z,2), to_famat_all(g, e));
118 else
119 z = mkvec2(NULL, to_famat_all(g, e));
120 continue;
121 }
122
123 gel(G,1) = g;
124 g = idealpowred(nf,G,e,prec);
125 z = z? idealmulred(nf,z,g,prec): g;
126 }
127 gel(z,2) = famat_reduce(gel(z,2));
128 gel(basecl,j) = z;
129 }
130 return basecl;
131 }
132
133 /* given two coprime integral ideals x and f (f HNF), compute "small"
134 * non-zero a in x, such that a = 1 mod (f). GTM 193: Algo 4.3.3 */
135 static GEN
redideal(GEN nf,GEN x,GEN f)136 redideal(GEN nf,GEN x,GEN f)
137 {
138 if (gcmp1(gcoeff(f,1,1))) return idealred_elt(nf, x); /* f = 1 */
139 return idealaddtoone_i(nf,x,f); /* a = b mod (x f), != 0 since 1 mod f */
140 }
141
142 static int
too_big(GEN nf,GEN bet)143 too_big(GEN nf, GEN bet)
144 {
145 GEN x = gnorm(coltoalg(nf,bet));
146 switch (typ(x))
147 {
148 case t_INT: return absi_cmp(x, gen_1);
149 case t_FRAC: return absi_cmp(gel(x,1), gel(x,2));
150 }
151 pari_err(bugparier, "wrong type in too_big");
152 return 0; /* not reached */
153 }
154
155 /* GTM 193: Algo 4.3.4. Reduce x mod idele */
156 static GEN
_idealmodidele(GEN nf,GEN x,GEN idele,GEN sarch)157 _idealmodidele(GEN nf, GEN x, GEN idele, GEN sarch)
158 {
159 pari_sp av = avma;
160 GEN a,A,D,G, f = gel(idele,1);
161
162 G = redideal(nf, x, f);
163 D = redideal(nf, idealdiv(nf,G,x), f);
164 A = element_div(nf,D,G);
165 if (too_big(nf,A) > 0) { avma = av; return x; }
166 a = set_sign_mod_idele(nf, NULL, A, idele, sarch);
167 if (a != A && too_big(nf,A) > 0) { avma = av; return x; }
168 return idealmul(nf, a, x);
169 }
170
171 GEN
idealmodidele(GEN bnr,GEN x)172 idealmodidele(GEN bnr, GEN x)
173 {
174 GEN bid = gel(bnr,2), fa2 = gel(bid,4);
175 GEN idele = gel(bid,1);
176 GEN sarch = (GEN)fa2[lg(fa2)-1];
177 return _idealmodidele(checknf(bnr), x, idele, sarch);
178 }
179
180 /* v_pr(L0 * cx). tau = pr[5] or (more efficient) mult. table for pr[5] */
181 static long
fast_val(GEN nf,GEN L0,GEN cx,GEN pr,GEN tau)182 fast_val(GEN nf,GEN L0,GEN cx,GEN pr,GEN tau)
183 {
184 pari_sp av = avma;
185 GEN p = gel(pr,1);
186 long v = int_elt_val(nf,L0,p,tau,NULL);
187 if (cx)
188 {
189 long w = ggval(cx, p);
190 if (w) v += w * itos(gel(pr,3));
191 }
192 avma = av; return v;
193 }
194
195 /* x coprime to fZ, return y = x mod fZ, y integral */
196 static GEN
make_integral_Z(GEN x,GEN fZ)197 make_integral_Z(GEN x, GEN fZ)
198 {
199 GEN d, y = Q_remove_denom(x, &d);
200 if (d) y = FpC_Fp_mul(y, Fp_inv(d, fZ), fZ);
201 return y;
202 }
203
204 /* p pi^(-1) mod f */
205 static GEN
get_pinvpi(GEN nf,GEN fZ,GEN p,GEN pi,GEN * v)206 get_pinvpi(GEN nf, GEN fZ, GEN p, GEN pi, GEN *v)
207 {
208 if (!*v) {
209 GEN invpi = element_inv(nf, pi);
210 *v = make_integral_Z(gmul(p, invpi), mulii(p, fZ));
211 }
212 return *v;
213 }
214 /* p pi^(-1) mod f */
215 static GEN
get_pi(GEN F,GEN pr,GEN * v)216 get_pi(GEN F, GEN pr, GEN *v)
217 {
218 if (!*v) *v = unif_mod_fZ(pr, F);
219 return *v;
220 }
221
222 static GEN
compute_raygen(GEN nf,GEN u1,GEN gen,GEN bid)223 compute_raygen(GEN nf, GEN u1, GEN gen, GEN bid)
224 {
225 GEN f, fZ, basecl, module, fa, fa2, pr, t, EX, sarch, cyc, F;
226 GEN *listpr, *vecpi, *vecpinvpi, *vectau;
227 long i,j,l,lp;
228
229 if (lg(u1) == 1) return cgetg(1, t_VEC);
230
231 /* basecl = generators in factored form */
232 basecl = compute_fact(nf,u1,gen);
233
234 module = gel(bid,1);
235 cyc = gmael(bid,2,2); EX = gel(cyc,1); /* exponent of (O/f)^* */
236 f = gel(module,1); fZ = gcoeff(f,1,1);
237 fa = gel(bid,3);
238 fa2 = gel(bid,4); sarch = (GEN)fa2[lg(fa2)-1];
239 listpr = (GEN*)fa[1]; F = init_unif_mod_fZ((GEN)listpr);
240
241 lp = lg(listpr);
242 vecpinvpi = (GEN*)cgetg(lp, t_VEC);
243 vecpi = (GEN*)cgetg(lp, t_VEC);
244 vectau = (GEN*)cgetg(lp, t_VEC);
245 for (i=1; i<lp; i++)
246 {
247 pr = listpr[i];
248 vecpi[i] = NULL; /* to be computed if needed */
249 vecpinvpi[i] = NULL; /* to be computed if needed */
250 vectau[i] = eltmul_get_table(nf, gel(pr,5));
251 }
252
253 l = lg(basecl);
254 for (i=1; i<l; i++)
255 {
256 GEN p, pi, pinvpi, dmulI, mulI, G, I, A, e, L, newL;
257 long la, v, k;
258 pari_sp av;
259 /* G = [I, A=famat(L,e)] is a generator, I integral */
260 G = gel(basecl,i);
261 I = gel(G,1);
262 A = gel(G,2);
263 L = gel(A,1);
264 e = gel(A,2);
265 /* if no reduction took place in compute_fact, everybody is still coprime
266 * to f + no denominators */
267 if (!I)
268 {
269 gel(basecl,i) = famat_to_nf_modidele(nf, L, e, bid);
270 continue;
271 }
272 if (lg(A) == 1)
273 {
274 gel(basecl,i) = I;
275 continue;
276 }
277
278 /* compute mulI so that mulI * I coprime to f
279 * FIXME: use idealcoprime ??? (Less efficient. Fix idealcoprime!) */
280 dmulI = mulI = NULL;
281 for (j=1; j<lp; j++)
282 {
283 pr = listpr[j];
284 v = idealval(nf, I, pr);
285 if (!v) continue;
286 p = gel(pr,1);
287 pi = get_pi(F, pr, &vecpi[j]);
288 pinvpi = get_pinvpi(nf, fZ, p, pi, &vecpinvpi[j]);
289 t = element_pow(nf, pinvpi, stoi(v));
290 mulI = mulI? element_mul(nf, mulI, t): t;
291 t = powiu(gel(pr,1), v);
292 dmulI = dmulI? mulii(dmulI, t): t;
293 }
294
295 /* make all components of L coprime to f.
296 * Assuming (L^e * I, f) = 1, then newL^e * mulI = L^e */
297 la = lg(e); newL = cgetg(la, t_VEC);
298 for (k=1; k<la; k++)
299 {
300 GEN L0, cx, LL = algtobasis_i(nf, gel(L,k));
301 L0 = Q_primitive_part(LL, &cx); /* LL = L0*cx (faster element_val) */
302 for (j=1; j<lp; j++)
303 {
304 pr = listpr[j];
305 v = fast_val(nf, L0,cx, pr,vectau[j]); /* = val_pr(LL) */
306 if (!v) continue;
307 p = gel(pr,1);
308 pi = get_pi(F, pr, &vecpi[j]);
309 if (v > 0)
310 {
311 pinvpi = get_pinvpi(nf, fZ, p, pi, &vecpinvpi[j]);
312 t = element_pow(nf,pinvpi,stoi(v));
313 LL = element_mul(nf, LL, t);
314 LL = gdiv(LL, powiu(p, v));
315 }
316 else
317 {
318 t = element_pow(nf,pi,stoi(-v));
319 LL = element_mul(nf, LL, t);
320 }
321 }
322 gel(newL,k) = FpC_red(make_integral(nf,LL,f,listpr), fZ);
323 }
324
325 av = avma;
326 /* G in nf, = L^e mod f */
327 G = famat_to_nf_modideal_coprime(nf, newL, e, f, EX);
328 if (mulI)
329 {
330 G = element_muli(nf, G, mulI);
331 G = colreducemodHNF(G, gmul(f, dmulI), NULL);
332 }
333 G = set_sign_mod_idele(nf,A,G,module,sarch);
334 I = idealmul(nf,I,G);
335 if (dmulI) I = gdivexact(I, dmulI);
336 /* more or less useless, but cheap at this point */
337 I = _idealmodidele(nf,I,module,sarch);
338 gel(basecl,i) = gerepilecopy(av, I);
339 }
340 return basecl;
341 }
342
343 /********************************************************************/
344 /** **/
345 /** INIT RAY CLASS GROUP **/
346 /** **/
347 /********************************************************************/
348 static GEN
get_dataunit(GEN bnf,GEN bid)349 get_dataunit(GEN bnf, GEN bid)
350 {
351 GEN D, cyc = gmael(bid,2,2), U = init_units(bnf), nf = gel(bnf,7);
352 long i, l;
353 zlog_S S; init_zlog_bid(&S, bid);
354 D = zsignunits(bnf, S.archp, 1); l = lg(D);
355 for (i = 1; i < l; i++)
356 gel(D,i) = vecmodii(gmul(S.U, zlog(nf, gel(U,i),gel(D,i), &S)), cyc);
357 return shallowconcat(D, diagonal_i(cyc));
358 }
359
360 static GEN
Buchray(GEN bnf,GEN module,long flag)361 Buchray(GEN bnf,GEN module,long flag)
362 {
363 GEN nf, cyc, gen, Gen, u, clg, logs, p1, h, met, u1, u2, U, cycgen;
364 GEN bigres, bid, cycbid, genbid, x, y, funits, H, El;
365 long RU, Ri, j, ngen, lh;
366 const long add_gen = flag & nf_GEN;
367 const long do_init = flag & nf_INIT;
368 pari_sp av = avma;
369
370 bnf = checkbnf(bnf); nf = checknf(bnf);
371 funits = check_units(bnf, "Buchray"); RU = lg(funits);
372 El = Gen = NULL; /* gcc -Wall */
373 bigres = gel(bnf,8);
374 cyc = gmael(bigres,1,2);
375 gen = gmael(bigres,1,3); ngen = lg(cyc)-1;
376
377 bid = Idealstar(nf,module,1);
378 cycbid = gmael(bid,2,2);
379 genbid = gmael(bid,2,3);
380 Ri = lg(cycbid)-1; lh = ngen+Ri;
381
382 x = gmael(bid,1,1);
383 if (Ri || add_gen || do_init)
384 {
385 GEN fx = gel(bid,3);
386 El = cgetg(ngen+1,t_VEC);
387 for (j=1; j<=ngen; j++)
388 {
389 p1 = idealcoprime_fact(nf, gel(gen,j), fx);
390 if (RgV_isscalar(p1)) p1 = gel(p1,1);
391 gel(El,j) = p1;
392 }
393 }
394 if (add_gen)
395 {
396 Gen = cgetg(lh+1,t_VEC);
397 for (j=1; j<=ngen; j++) gel(Gen,j) = idealmul(nf,gel(El,j),gel(gen,j));
398 for ( ; j<=lh; j++) Gen[j] = genbid[j - ngen];
399 }
400 if (!Ri)
401 {
402 clg = cgetg(add_gen? 4: 3,t_VEC);
403 if (add_gen) gel(clg,3) = Gen;
404 gel(clg,1) = gmael(bigres,1,1);
405 gel(clg,2) = cyc;
406 if (!do_init) return gerepilecopy(av,clg);
407 y = cgetg(7,t_VEC);
408 gel(y,1) = bnf;
409 gel(y,2) = bid;
410 gel(y,3) = El;
411 gel(y,4) = matid(ngen);
412 gel(y,5) = clg;
413 gel(y,6) = mkvec2(cgetg(1,t_MAT), matid(RU));
414 return gerepilecopy(av,y);
415 }
416
417 cycgen = check_and_build_cycgen(bnf);
418 /* (log(Units)|D) * u = (0 | H) */
419 H = hnfall_i( get_dataunit(bnf, bid), do_init? &u: NULL, 1);
420 logs = cgetg(ngen+1, t_MAT);
421 /* FIXME: cycgen[j] is not necessarily coprime to bid, but it is made coprime
422 * in famat_zlog using canonical uniformizers [from bid data]: no need to
423 * correct it here. The same ones will be used in bnrisprincipal. Hence
424 * modification by El is useless. */
425 for (j=1; j<=ngen; j++)
426 {
427 p1 = gel(cycgen,j);
428 if (typ(El[j]) != t_INT) /* <==> != 1 */
429 {
430 GEN F = to_famat_all(gel(El,j), gel(cyc,j));
431 p1 = arch_mul(F, p1);
432 }
433 gel(logs,j) = zideallog(nf, p1, bid); /* = log(Gen[j]) */
434 }
435 /* [ cyc 0 ]
436 * [-logs H ] = relation matrix for Cl_f */
437 h = shallowconcat(
438 vconcat(diagonal_i(cyc), gneg_i(logs)),
439 vconcat(zeromat(ngen, Ri), H)
440 );
441 met = smithrel(hnf(h), &U, add_gen? &u1: NULL);
442 clg = cgetg(add_gen? 4: 3, t_VEC);
443 gel(clg,1) = detcyc(met, &j);
444 gel(clg,2) = met;
445 if (add_gen) gel(clg,3) = compute_raygen(nf,u1,Gen,bid);
446 if (!do_init) return gerepilecopy(av, clg);
447
448 u2 = cgetg(Ri+1,t_MAT);
449 u1 = cgetg(RU+1,t_MAT);
450 for (j=1; j<=RU; j++) { u1[j]=u[j]; setlg(u[j],RU+1); }
451 u += RU;
452 for (j=1; j<=Ri; j++) { u2[j]=u[j]; setlg(u[j],RU+1); }
453
454 /* log(Units) U2 = H (mod D)
455 * log(Units) U1 = 0 (mod D) */
456 u1 = lllint_ip(u1,100);
457 u2 = gmul(reducemodinvertible(u2,u1), ginv(H));
458 y = cgetg(7,t_VEC);
459 gel(y,1) = bnf;
460 gel(y,2) = bid;
461 gel(y,3) = El;
462 gel(y,4) = U;
463 gel(y,5) = clg;
464 gel(y,6) = mkvec2(u2,u1);
465 return gerepilecopy(av,y);
466 }
467
468 GEN
buchrayinitgen(GEN bnf,GEN ideal)469 buchrayinitgen(GEN bnf, GEN ideal)
470 { return Buchray(bnf,ideal, nf_INIT | nf_GEN); }
471 GEN
buchrayinit(GEN bnf,GEN ideal)472 buchrayinit(GEN bnf, GEN ideal)
473 { return Buchray(bnf,ideal, nf_INIT); }
474 GEN
buchray(GEN bnf,GEN ideal)475 buchray(GEN bnf, GEN ideal)
476 { return Buchray(bnf,ideal, nf_GEN); }
477
478 GEN
bnrclass0(GEN bnf,GEN ideal,long flag)479 bnrclass0(GEN bnf, GEN ideal, long flag)
480 {
481 switch(flag)
482 {
483 case 0: flag = nf_GEN; break;
484 case 1: flag = nf_INIT; break;
485 case 2: flag = nf_INIT | nf_GEN; break;
486 default: pari_err(flagerr,"bnrclass");
487 }
488 return Buchray(bnf,ideal,flag);
489 }
490
491 GEN
bnrinit0(GEN bnf,GEN ideal,long flag)492 bnrinit0(GEN bnf, GEN ideal, long flag)
493 {
494 switch(flag)
495 {
496 case 0: flag = nf_INIT; break;
497 case 1: flag = nf_INIT | nf_GEN; break;
498 default: pari_err(flagerr,"bnrinit");
499 }
500 return Buchray(bnf,ideal,flag);
501 }
502
503 GEN
bnrclassno(GEN bnf,GEN ideal)504 bnrclassno(GEN bnf,GEN ideal)
505 {
506 GEN nf, h, D, bigres, bid, cycbid;
507 pari_sp av = avma;
508
509 bnf = checkbnf(bnf); nf = gel(bnf,7);
510 bigres = gel(bnf,8); h = gmael(bigres,1,1); /* class number */
511 bid = Idealstar(nf,ideal,0);
512 cycbid = gmael(bid,2,2);
513 if (lg(cycbid) == 1) { avma = av; return icopy(h); }
514 D = get_dataunit(bnf, bid); /* (Z_K/f)^* / units ~ Z^n / D */
515 return gerepileuptoint(av, mulii(h, dethnf_i(hnf(D))));
516 }
517
518 GEN
quick_isprincipalgen(GEN bnf,GEN x)519 quick_isprincipalgen(GEN bnf, GEN x)
520 { /* x \prod g[i]^(-ep[i]) = factorisation of principal ideal */
521 GEN idep, gen = gmael3(bnf,8,1,3), ep = isprincipal(bnf,x);
522 idep = isprincipalfact(bnf, gen, gneg(ep), x, nf_GENMAT|nf_FORCE);
523 return mkvec2(ep, gel(idep,2));
524 }
525
526 GEN
bnrisprincipal(GEN bnr,GEN x,long flag)527 bnrisprincipal(GEN bnr, GEN x, long flag)
528 {
529 long i, j, c;
530 pari_sp av = avma;
531 GEN bnf, nf, bid, U, El, ep, p1, beta, idep, ex, rayclass, divray, genray;
532 GEN alpha;
533
534 checkbnr(bnr); rayclass = gel(bnr,5);
535 divray = gel(rayclass,2); c = lg(divray);
536 ex = cgetg(c,t_COL);
537 if (c == 1 && !(flag & nf_GEN)) return ex;
538
539 bnf = gel(bnr,1); nf = gel(bnf,7);
540 bid = gel(bnr,2);
541 El = gel(bnr,3);
542 U = gel(bnr,4);
543
544 if (typ(x) == t_VEC && lg(x) == 3)
545 { idep = gel(x,2); x = gel(x,1); } /* precomputed */
546 else
547 idep = quick_isprincipalgen(bnf, x);
548 ep = gel(idep,1);
549 beta= gel(idep,2);
550 j = lg(ep);
551 for (i=1; i<j; i++) /* modify beta as if gen -> El.gen (coprime to bid) */
552 if (typ(El[i]) != t_INT && signe(ep[i])) /* <==> != 1 */
553 beta = arch_mul(to_famat_all(gel(El,i), negi(gel(ep,i))), beta);
554 p1 = gmul(U, shallowconcat(ep, zideallog(nf,beta,bid)));
555 ex = vecmodii(p1, divray);
556 if (!(flag & nf_GEN)) return gerepileupto(av, ex);
557
558 /* compute generator */
559 if (lg(rayclass)<=3)
560 pari_err(talker,"please apply bnrinit(,,1) and not bnrinit(,,0)");
561
562 genray = gel(rayclass,3);
563 p1 = isprincipalfact(bnf, genray, gneg(ex), x, nf_GENMAT | nf_FORCE);
564 if (!gcmp0(gel(p1,1))) pari_err(bugparier,"isprincipalray");
565 p1 = gel(p1,2); alpha = factorbackelt(p1, nf, NULL);
566 if (lg(bid[5]) > 1 && lg(gmael(bid,5,1)) > 1)
567 {
568 GEN u = gel(bnr,6), y = gmul(gel(u,1), zideallog(nf, p1, bid));
569 y = reducemodinvertible(y, gel(u,2));
570 alpha = element_div(nf, alpha, factorbackelt(init_units(bnf), y, nf));
571 }
572 return gerepilecopy(av, mkvec2(ex,alpha));
573 }
574
575 GEN
isprincipalray(GEN bnr,GEN x)576 isprincipalray(GEN bnr, GEN x)
577 {
578 return bnrisprincipal(bnr,x,0);
579 }
580
581 GEN
isprincipalraygen(GEN bnr,GEN x)582 isprincipalraygen(GEN bnr, GEN x)
583 {
584 return bnrisprincipal(bnr,x,nf_GEN);
585 }
586
587 /* N! / N^N * (4/pi)^r2 * sqrt(|D|) */
588 GEN
minkowski_bound(GEN D,long N,long r2,long prec)589 minkowski_bound(GEN D, long N, long r2, long prec)
590 {
591 pari_sp av = avma;
592 GEN p1;
593 p1 = gdiv(mpfactr(N,prec), powuu(N,N));
594 p1 = gmul(p1, gpowgs(gdivsg(4,mppi(prec)), r2));
595 p1 = gmul(p1, gsqrt(absi(D),prec));
596 return gerepileupto(av, p1);
597 }
598
599 /* DK = |dK| */
600 static long
zimmertbound(long N,long R2,GEN DK)601 zimmertbound(long N,long R2,GEN DK)
602 {
603 pari_sp av = avma;
604 GEN w;
605 long n;
606
607 if (N < 2) return 1;
608 if (N < 21)
609 {
610 static double c[19][11] = {
611 {/*2*/ 0.6931, 0.45158},
612 {/*3*/ 1.71733859, 1.37420604},
613 {/*4*/ 2.91799837, 2.50091538, 2.11943331},
614 {/*5*/ 4.22701425, 3.75471588, 3.31196660},
615 {/*6*/ 5.61209925, 5.09730381, 4.60693851, 4.14303665},
616 {/*7*/ 7.05406203, 6.50550021, 5.97735406, 5.47145968},
617 {/*8*/ 8.54052636, 7.96438858, 7.40555445, 6.86558259, 6.34608077},
618 {/*9*/ 10.0630022, 9.46382812, 8.87952524, 8.31139202, 7.76081149},
619 {/*10*/11.6153797, 10.9966020, 10.3907654, 9.79895170, 9.22232770, 8.66213267},
620 {/*11*/13.1930961, 12.5573772, 11.9330458, 11.3210061, 10.7222412, 10.1378082},
621 {/*12*/14.7926394, 14.1420915, 13.5016616, 12.8721114, 12.2542699, 11.6490374,
622 11.0573775},
623 {/*13*/16.4112395, 15.7475710, 15.0929680, 14.4480777, 13.8136054, 13.1903162,
624 12.5790381},
625 {/*14*/18.0466672, 17.3712806, 16.7040780, 16.0456127, 15.3964878, 14.7573587,
626 14.1289364, 13.5119848},
627 {/*15*/19.6970961, 19.0111606, 18.3326615, 17.6620757, 16.9999233, 16.3467686,
628 15.7032228, 15.0699480},
629 {/*16*/21.3610081, 20.6655103, 19.9768082, 19.2953176, 18.6214885, 17.9558093,
630 17.2988108, 16.6510652, 16.0131906},
631
632 {/*17*/23.0371259, 22.3329066, 21.6349299, 20.9435607, 20.2591899, 19.5822454,
633 18.9131878, 18.2525157, 17.6007672},
634
635 {/*18*/24.7243611, 24.0121449, 23.3056902, 22.6053167, 21.9113705, 21.2242247,
636 20.5442836, 19.8719830, 19.2077941, 18.5522234},
637
638 {/*19*/26.4217792, 25.7021950, 24.9879497, 24.2793271, 23.5766321, 22.8801952,
639 22.1903709, 21.5075437, 20.8321263, 20.1645647},
640 {/*20*/28.1285704, 27.4021674, 26.6807314, 25.9645140, 25.2537867, 24.5488420,
641 23.8499943, 23.1575823, 22.4719720, 21.7935548, 21.1227537}
642 };
643 w = gmul(dbltor(exp(-c[N-2][R2])), gsqrt(DK,DEFAULTPREC));
644 }
645 else
646 {
647 w = minkowski_bound(DK, N, R2, DEFAULTPREC);
648 }
649 n = itos_or_0( gceil(w) );
650 if (!n) pari_err(talker,"Minkowski bound is too large");
651 if (n > 500000)
652 pari_warn(warner,"large Minkowski bound: certification will be VERY long");
653 avma = av; return n;
654 }
655
656 /* return \gamma_n^n if known, an upper bound otherwise */
657 static GEN
hermiteconstant(long n)658 hermiteconstant(long n)
659 {
660 GEN h,h1;
661 pari_sp av;
662
663 switch(n)
664 {
665 case 1: return gen_1;
666 case 2: return mkfrac(utoipos(4), utoipos(3));
667 case 3: return gen_2;
668 case 4: return utoipos(4);
669 case 5: return utoipos(8);
670 case 6: return mkfrac(utoipos(64), utoipos(3));
671 case 7: return utoipos(64);
672 case 8: return utoipos(256);
673 }
674 av = avma;
675 h = gpowgs(divsr(2,mppi(DEFAULTPREC)), n);
676 h1 = gsqr(ggamma(gdivgs(utoipos(n+4),2),DEFAULTPREC));
677 return gerepileupto(av, gmul(h,h1));
678 }
679
680 /* 1 if L (= nf != Q) primitive for sure, 0 if MAYBE imprimitive (may have a
681 * subfield K) */
682 static long
isprimitive(GEN nf)683 isprimitive(GEN nf)
684 {
685 long p, i, l, ep, N = degpol(nf[1]);
686 GEN d,fa;
687
688 fa = (GEN)factor(utoipos(N))[1]; /* primes | N */
689 p = itos(gel(fa,1)); if (p == N) return 1; /* prime degree */
690
691 /* N = [L:Q] = product of primes >= p, same is true for [L:K]
692 * d_L = t d_K^[L:K] --> check that some q^p divides d_L */
693 d = absi(gel(nf,3));
694 fa = (GEN)auxdecomp(d,0)[2]; /* list of v_q(d_L). Don't check large primes */
695 if (mod2(d)) i = 1;
696 else
697 { /* q = 2 */
698 ep = itos(gel(fa,1));
699 if ((ep>>1) >= p) return 0; /* 2 | d_K ==> 4 | d_K */
700 i = 2;
701 }
702 l = lg(fa);
703 for ( ; i < l; i++)
704 {
705 ep = itos(gel(fa,i));
706 if (ep >= p) return 0;
707 }
708 return 1;
709 }
710
711 static GEN
dft_bound()712 dft_bound()
713 {
714 if (DEBUGLEVEL>1) fprintferr("Default bound for regulator: 0.2\n");
715 return dbltor(0.2);
716 }
717
718 static GEN
regulatorbound(GEN bnf)719 regulatorbound(GEN bnf)
720 {
721 long N, R1, R2, R;
722 GEN nf, dK, p1, c1;
723
724 nf = gel(bnf,7); N = degpol(nf[1]);
725 if (!isprimitive(nf)) return dft_bound();
726
727 dK = absi(gel(nf,3));
728 nf_get_sign(nf, &R1, &R2); R = R1+R2-1;
729 c1 = (!R2 && N<12)? int2n(N & (~1UL)): powuu(N,N);
730 if (cmpii(dK,c1) <= 0) return dft_bound();
731
732 p1 = gsqr(glog(gdiv(dK,c1),DEFAULTPREC));
733 p1 = divrs(gmul2n(gpowgs(divrs(mulrs(p1,3),N*(N*N-1)-6*R2),R),R2), N);
734 p1 = sqrtr(gdiv(p1, hermiteconstant(R)));
735 if (DEBUGLEVEL>1) fprintferr("Mahler bound for regulator: %Z\n",p1);
736 return gmax(p1, dbltor(0.2));
737 }
738
739 /* x given by its embeddings */
740 GEN
norm_by_embed(long r1,GEN x)741 norm_by_embed(long r1, GEN x)
742 {
743 long i, ru = lg(x)-1;
744 GEN p = gel(x,ru);
745 if (r1 == ru)
746 {
747 for (i=ru-1; i>0; i--) p = gmul(p, gel(x,i));
748 return p;
749 }
750 p = gnorm(p);
751 for (i=ru-1; i>r1; i--) p = gmul(p, gnorm(gel(x,i)));
752 for ( ; i>0 ; i--) p = gmul(p, gel(x,i));
753 return p;
754 }
755
756 static int
is_unit(GEN M,long r1,GEN x)757 is_unit(GEN M, long r1, GEN x)
758 {
759 pari_sp av = avma;
760 GEN Nx = ground( norm_by_embed(r1, RgM_zc_mul(M,x)) );
761 int ok = is_pm1(Nx);
762 avma = av; return ok;
763 }
764
765 /* FIXME: should use smallvectors */
766 static GEN
minimforunits(GEN nf,long BORNE,GEN w)767 minimforunits(GEN nf, long BORNE, GEN w)
768 {
769 const long prec = MEDDEFAULTPREC;
770 long n, i, j, k, s, *x, r1, cnt = 0;
771 pari_sp av = avma;
772 GEN u,r,a,M;
773 double p, norme, normin, normax;
774 double **q,*v,*y,*z;
775 double eps=0.000001, BOUND = BORNE * 1.00001;
776
777 if (DEBUGLEVEL>=2)
778 {
779 fprintferr("Searching minimum of T2-form on units:\n");
780 if (DEBUGLEVEL>2) fprintferr(" BOUND = %ld\n",BORNE);
781 flusherr();
782 }
783 r1 = nf_get_r1(nf); n = degpol(nf[1]);
784 minim_alloc(n+1, &q, &x, &y, &z, &v);
785 M = gprec_w(gmael(nf,5,1), prec);
786 a = gmul(gmael(nf,5,2), real_1(prec));
787 r = sqred1_from_QR(a, prec);
788 for (j=1; j<=n; j++)
789 {
790 v[j] = rtodbl(gcoeff(r,j,j));
791 for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j));
792 }
793 normax = 0.; normin = (double)BOUND;
794 s=0; k=n; y[n]=z[n]=0;
795 x[n] = (long)(sqrt(BOUND/v[n]));
796
797 for(;;)
798 {
799 do
800 {
801 if (k>1)
802 {
803 long l = k-1;
804 z[l] = 0;
805 for (j=k; j<=n; j++) z[l] = z[l]+q[l][j]*x[j];
806 p = (double)x[k] + z[k];
807 y[l] = y[k]+p*p*v[k];
808 x[l] = (long)floor(sqrt((BOUND-y[l])/v[l])-z[l]);
809 k = l;
810 }
811 for(;;)
812 {
813 p = (double)x[k] + z[k];
814 if (y[k] + p*p*v[k] <= BOUND) break;
815 k++; x[k]--;
816 }
817 }
818 while (k>1);
819 if (!x[1] && y[1]<=eps) break;
820
821 if (DEBUGLEVEL>8){ fprintferr("."); flusherr(); }
822 if (++cnt == 5000) return NULL; /* too expensive */
823
824 p = (double)x[1] + z[1]; norme = y[1] + p*p*v[1] + eps;
825 if (norme > normax) normax = norme;
826 if (is_unit(M,r1, x)
827 && (norme > 2*n /* exclude roots of unity */
828 || !RgV_isscalar(element_pow(nf, zc_to_ZC(x), w))))
829 {
830 if (norme < normin) normin = norme;
831 if (DEBUGLEVEL>=2) { fprintferr("*"); flusherr(); }
832 }
833 x[k]--;
834 }
835 if (DEBUGLEVEL>=2){ fprintferr("\n"); flusherr(); }
836 avma = av; u = cgetg(4,t_VEC);
837 gel(u,1) = stoi(s<<1);
838 gel(u,2) = dbltor(normax);
839 gel(u,3) = dbltor(normin);
840 return u;
841 }
842
843 #undef NBMAX
844 static int
is_zero(GEN x,long bitprec)845 is_zero(GEN x, long bitprec) { return (gexpo(x) < -bitprec); }
846
847 static int
is_complex(GEN x,long bitprec)848 is_complex(GEN x, long bitprec) { return !is_zero(imag_i(x), bitprec); }
849
850 /* assume M_star t_REAL
851 * FIXME: what does this do ? To be rewritten */
852 static GEN
compute_M0(GEN M_star,long N)853 compute_M0(GEN M_star,long N)
854 {
855 long m1,m2,n1,n2,n3,lr,lr1,lr2,i,j,l,vx,vy,vz,vM;
856 GEN pol,p1,p2,p3,p4,p5,p6,p7,p8,p9,u,v,w,r,r1,r2,M0,M0_pro,S,P,M;
857 GEN f1,f2,f3,g1,g2,g3,pg1,pg2,pg3,pf1,pf2,pf3,X,Y,Z;
858 long bitprec = 24;
859
860 if (N == 2) return gmul2n(gsqr(gach(gmul2n(M_star,-1),0)), -1);
861 vM = fetch_var(); M = pol_x[vM];
862 vz = fetch_var(); Z = pol_x[vz];
863 vy = fetch_var(); Y = pol_x[vy];
864 vx = fetch_var(); X = pol_x[vx];
865
866 M0 = NULL; m1 = N/3;
867 for (n1=1; n1<=m1; n1++)
868 {
869 m2 = (N-n1)>>1;
870 for (n2=n1; n2<=m2; n2++)
871 {
872 pari_sp av = avma; n3=N-n1-n2;
873 if (n1==n2 && n1==n3) /* n1 = n2 = n3 = m1 = N/3 */
874 {
875 p1 = divrs(M_star, m1);
876 p4 = sqrtr_abs( mulrr(addsr(1,p1),subrs(p1,3)) );
877 p5 = subrs(p1,1);
878 u = gen_1;
879 v = gmul2n(addrr(p5,p4),-1);
880 w = gmul2n(subrr(p5,p4),-1);
881 M0_pro=gmul2n(mulsr(m1,addrr(gsqr(logr_abs(v)),gsqr(logr_abs(w)))), -2);
882 if (DEBUGLEVEL>2)
883 {
884 fprintferr("[ %ld, %ld, %ld ]: %Z\n",n1,n2,n3,gprec_w(M0_pro,3));
885 flusherr();
886 }
887 if (!M0 || gcmp(M0_pro,M0) < 0) M0 = M0_pro;
888 }
889 else if (n1==n2 || n2==n3)
890 { /* n3 > N/3 >= n1 */
891 long k = N - 2*n2;
892 p2 = gsub(M_star, gmulgs(X,n2));
893 p3 = gmul(powuu(k,k),
894 gpowgs(gsubgs(gmul(M_star,p2),k*k),n2));
895 pol = gsub(p3, gmul(gmul(powuu(n2,n2),gpowgs(X,n2)),
896 gpowgs(p2, N-n2)));
897 r = roots(pol, DEFAULTPREC); lr = lg(r);
898 for (i=1; i<lr; i++)
899 {
900 S = real_i(gel(r,i));
901 if (is_complex(gel(r,i), bitprec) || signe(S) <= 0) continue;
902
903 p4 = subrr(M_star, mulsr(n2,S));
904 P = divrr(mulrr(mulsr(n2,S),p4), subrs(mulrr(M_star,p4),k*k));
905 p5 = subrr(gsqr(S), gmul2n(P,2));
906 if (gsigne(p5) < 0) continue;
907
908 p6 = sqrtr(p5);
909 v = gmul2n(subrr(S,p6),-1);
910 if (gsigne(v) <= 0) continue;
911
912 u = gmul2n(addrr(S,p6),-1);
913 w = gpow(P, gdivgs(utoineg(n2),k), 0);
914 p6 = mulsr(n2, addrr(gsqr(logr_abs(u)), gsqr(logr_abs(v))));
915 M0_pro = gmul2n(addrr(p6, mulsr(k, gsqr(logr_abs(w)))),-2);
916 if (DEBUGLEVEL>2)
917 {
918 fprintferr("[ %ld, %ld, %ld ]: %Z\n",n1,n2,n3,gprec_w(M0_pro,3));
919 flusherr();
920 }
921 if (!M0 || gcmp(M0_pro,M0) < 0) M0 = M0_pro;
922 }
923 }
924 else
925 {
926 f1 = gsub(gadd(gmulsg(n1,X),gadd(gmulsg(n2,Y),gmulsg(n3,Z))), M);
927 f2 = gmulsg(n1,gmul(Y,Z));
928 f2 = gadd(f2,gmulsg(n2,gmul(X,Z)));
929 f2 = gadd(f2,gmulsg(n3,gmul(X,Y)));
930 f2 = gsub(f2,gmul(M,gmul(X,gmul(Y,Z))));
931 f3 = gsub(gmul(gpowgs(X,n1),gmul(gpowgs(Y,n2),gpowgs(Z,n3))), gen_1);
932 /* f1 = n1 X + n2 Y + n3 Z - M */
933 /* f2 = n1 YZ + n2 XZ + n3 XY */
934 /* f3 = X^n1 Y^n2 Z^n3 - 1*/
935 g1=subres(f1,f2); g1=gdiv(g1,content(g1));
936 g2=subres(f1,f3); g2=gdiv(g2,content(g2));
937 g3=subres(g1,g2); g3=gdiv(g3,content(g3));
938 pf1=gsubst(f1,vM,M_star); pg1=gsubst(g1,vM,M_star);
939 pf2=gsubst(f2,vM,M_star); pg2=gsubst(g2,vM,M_star);
940 pf3=gsubst(f3,vM,M_star); pg3=gsubst(g3,vM,M_star);
941 /* g3 = Res_Y,Z(f1,f2,f3) */
942 r = roots(pg3,DEFAULTPREC); lr = lg(r);
943 for (i=1; i<lr; i++)
944 {
945 w = real_i(gel(r,i));
946 if (is_complex(gel(r,i), bitprec) || signe(w) <= 0) continue;
947 p1=gsubst(pg1,vz,w);
948 p2=gsubst(pg2,vz,w);
949 p3=gsubst(pf1,vz,w);
950 p4=gsubst(pf2,vz,w);
951 p5=gsubst(pf3,vz,w);
952 r1 = roots(p1, DEFAULTPREC); lr1 = lg(r1);
953 for (j=1; j<lr1; j++)
954 {
955 v = real_i(gel(r1,j));
956 if (is_complex(gel(r1,j), bitprec) || signe(v) <= 0
957 || !is_zero(gsubst(p2,vy,v), bitprec)) continue;
958
959 p7=gsubst(p3,vy,v);
960 p8=gsubst(p4,vy,v);
961 p9=gsubst(p5,vy,v);
962 r2 = roots(p7, DEFAULTPREC); lr2 = lg(r2);
963 for (l=1; l<lr2; l++)
964 {
965 u = real_i(gel(r2,l));
966 if (is_complex(gel(r2,l), bitprec) || signe(u) <= 0
967 || !is_zero(gsubst(p8,vx,u), bitprec)
968 || !is_zero(gsubst(p9,vx,u), bitprec)) continue;
969
970 M0_pro = mulsr(n1, gsqr(logr_abs(u)));
971 M0_pro = gadd(M0_pro, mulsr(n2, gsqr(logr_abs(v))));
972 M0_pro = gadd(M0_pro, mulsr(n3, gsqr(logr_abs(w))));
973 M0_pro = gmul2n(M0_pro,-2);
974 if (DEBUGLEVEL>2)
975 {
976 fprintferr("[ %ld, %ld, %ld ]: %Z\n",n1,n2,n3,gprec_w(M0_pro,3));
977 flusherr();
978 }
979 if (!M0 || gcmp(M0_pro,M0) < 0) M0 = M0_pro;
980 }
981 }
982 }
983 }
984 if (!M0) avma = av; else M0 = gerepilecopy(av, M0);
985 }
986 }
987 for (i=1;i<=4;i++) (void)delete_var();
988 return M0? M0: gen_0;
989 }
990
991 static GEN
lowerboundforregulator_i(GEN bnf)992 lowerboundforregulator_i(GEN bnf)
993 {
994 long N,R1,R2,RU,i;
995 GEN nf,M0,M,G,bound,minunit,newminunit;
996 GEN vecminim,p1,pol,y;
997 GEN units = check_units(bnf,"bnfcertify");
998
999 nf = gel(bnf,7); N = degpol(nf[1]);
1000 nf_get_sign(nf, &R1, &R2); RU = R1+R2-1;
1001 if (!RU) return gen_1;
1002
1003 G = gmael(nf,5,2);
1004 units = algtobasis(bnf,units);
1005 minunit = gnorml2(gmul(G, gel(units,1))); /* T2(units[1]) */
1006 for (i=2; i<=RU; i++)
1007 {
1008 newminunit = gnorml2(gmul(G, gel(units,i)));
1009 if (gcmp(newminunit,minunit) < 0) minunit = newminunit;
1010 }
1011 if (gexpo(minunit) > 30) return NULL;
1012
1013 vecminim = minimforunits(nf, itos(gceil(minunit)), gmael3(bnf,8,4,1));
1014 if (!vecminim) return NULL;
1015 bound = gel(vecminim,3);
1016 if (DEBUGLEVEL>1)
1017 {
1018 fprintferr("M* = %Z\n", bound);
1019 if (DEBUGLEVEL>2)
1020 {
1021 pol = gaddgs(gsub(monomial(gen_1,N,0),monomial(bound,1,0)),N-1);
1022 p1 = roots(pol,DEFAULTPREC);
1023 y= real_i((GEN)p1[ N&1? 3: 2]);
1024 M0 = gmul2n(gmulsg(N*(N-1),gsqr(glog(y,DEFAULTPREC))),-2);
1025 fprintferr("pol = %Z\n",pol);
1026 fprintferr("old method: y = %Z, M0 = %Z\n",y,gprec_w(M0,3));
1027 }
1028 }
1029 M0 = compute_M0(bound, N);
1030 if (DEBUGLEVEL>1) { fprintferr("M0 = %Z\n",gprec_w(M0,3)); flusherr(); }
1031 M = gmul2n(gdivgs(gdiv(gpowgs(M0,RU),hermiteconstant(RU)),N),R2);
1032 if (gcmp(M, dbltor(0.04)) < 0) return NULL;
1033 M = gsqrt(M,DEFAULTPREC);
1034 if (DEBUGLEVEL>1)
1035 fprintferr("(lower bound for regulator) M = %Z\n",gprec_w(M,3));
1036 return M;
1037 }
1038
1039 static GEN
lowerboundforregulator(GEN bnf)1040 lowerboundforregulator(GEN bnf)
1041 {
1042 pari_sp av = avma;
1043 GEN x = lowerboundforregulator_i(bnf);
1044 if (!x) { avma = av; x = regulatorbound(bnf); }
1045 return x;
1046 }
1047
1048 /* Compute a square matrix of rank length(beta) associated to a family
1049 * (P_i), 1<=i<=length(beta), of primes s.t. N(P_i) = 1 mod p, and
1050 * (P_i,beta[j]) = 1 for all i,j */
1051 static void
primecertify(GEN bnf,GEN beta,ulong p,GEN bad)1052 primecertify(GEN bnf, GEN beta, ulong p, GEN bad)
1053 {
1054 long i, j, nbcol, lb, nbqq, ra;
1055 GEN nf,mat,gq,LQ,newcol,g,ord,modpr;
1056 ulong q;
1057
1058 ord = NULL; /* gcc -Wall */
1059 nbcol = 0; nf = gel(bnf,7);
1060 lb = lg(beta)-1; mat = cgetg(1,t_MAT); q = 1UL;
1061 for(;;)
1062 {
1063 q += 2*p;
1064 if (!umodiu(bad,q) || !uisprime(q)) continue;
1065
1066 gq = utoipos(q);
1067 LQ = primedec(bnf,gq); nbqq = lg(LQ)-1;
1068 g = NULL;
1069 for (i=1; i<=nbqq; i++)
1070 {
1071 GEN mat1, Q = gel(LQ,i); if (!gcmp1(gel(Q,4))) break;
1072 /* Q has degree 1 */
1073 if (!g)
1074 {
1075 ord = Z_factor( utoipos(q-1) );
1076 g = gener_Fp_local(gq, gel(ord,1)); /* primitive root */
1077 }
1078 modpr = zkmodprinit(nf, Q);
1079 newcol = cgetg(lb+1,t_COL);
1080 for (j=1; j<=lb; j++)
1081 {
1082 GEN t = to_Fp_simple(nf, gel(beta,j), modpr);
1083 gel(newcol,j) = Fp_PHlog(t,g,gq,ord);
1084 }
1085 if (DEBUGLEVEL>3)
1086 {
1087 if (i==1) fprintferr(" generator of (Zk/Q)^*: %Z\n", g);
1088 fprintferr(" prime ideal Q: %Z\n",Q);
1089 fprintferr(" column #%ld of the matrix log(b_j/Q): %Z\n",
1090 nbcol, newcol);
1091 }
1092 mat1 = shallowconcat(mat,newcol); ra = rank(mat1);
1093 if (ra==nbcol) continue;
1094
1095 if (DEBUGLEVEL>2) fprintferr(" new rank: %ld\n",ra);
1096 if (++nbcol == lb) return;
1097 mat = mat1;
1098 }
1099 }
1100 }
1101
1102 static void
check_prime(ulong p,GEN bnf,GEN cyc,GEN cycgen,GEN fu,GEN mu,GEN bad)1103 check_prime(ulong p, GEN bnf, GEN cyc, GEN cycgen, GEN fu, GEN mu, GEN bad)
1104 {
1105 pari_sp av = avma;
1106 long i,b, lc = lg(cyc), w = itos(gel(mu,1)), lf = lg(fu);
1107 GEN beta = cgetg(lf+lc, t_VEC);
1108
1109 if (DEBUGLEVEL>1) fprintferr(" *** testing p = %lu\n",p);
1110 for (b=1; b<lc; b++)
1111 {
1112 if (umodiu(gel(cyc,b), p)) break; /* p \nmid cyc[b] */
1113 if (b==1 && DEBUGLEVEL>2) fprintferr(" p divides h(K)\n");
1114 beta[b] = cycgen[b];
1115 }
1116 if (w % p == 0)
1117 {
1118 if (DEBUGLEVEL>2) fprintferr(" p divides w(K)\n");
1119 beta[b++] = mu[2];
1120 }
1121 for (i=1; i<lf; i++) beta[b++] = fu[i];
1122 setlg(beta, b); /* beta = [cycgen[i] if p|cyc[i], tu if p|w, fu] */
1123 if (DEBUGLEVEL>3) {fprintferr(" Beta list = %Z\n",beta); flusherr();}
1124 primecertify(bnf,beta,p,bad); avma = av;
1125 }
1126
1127 long
certifybuchall(GEN bnf)1128 certifybuchall(GEN bnf)
1129 {
1130 pari_sp av = avma;
1131 long nbgen, i, N, R1, R2;
1132 GEN bad, nf, reg, zu, funits, gen, cycgen, cyc;
1133 byteptr delta = diffptr;
1134 ulong bound, p;
1135
1136 bnf = checkbnf(bnf); nf = gel(bnf,7);
1137 N=degpol(nf[1]); if (N==1) return 1;
1138 nf_get_sign(nf, &R1, &R2);
1139 funits = check_units(bnf,"bnfcertify");
1140 testprimes(bnf, zimmertbound(N,R2,absi(gel(nf,3))));
1141 reg = gmael(bnf,8,2);
1142 cyc = gmael3(bnf,8,1,2); nbgen = lg(cyc)-1;
1143 gen = gmael3(bnf,8,1,3); zu = gmael(bnf,8,4);
1144 bound = itou_or_0( ground(gdiv(reg, lowerboundforregulator(bnf))) );
1145 if (!bound) pari_err(talker,"sorry, too many primes to check");
1146 maxprime_check(bound);
1147 if (DEBUGLEVEL>1)
1148 {
1149 fprintferr("\nPHASE 2: are all primes good ?\n\n");
1150 fprintferr(" Testing primes <= B (= %lu)\n\n",bound); flusherr();
1151 }
1152 cycgen = check_and_build_cycgen(bnf);
1153 for (bad=gen_1,i=1; i<=nbgen; i++)
1154 bad = lcmii(bad, gcoeff(gen[i],1,1));
1155 for (i=1; i<=nbgen; i++)
1156 {
1157 GEN p1 = gel(cycgen,i);
1158 long j;
1159 if (typ(p1) == t_MAT)
1160 {
1161 GEN h, g = gel(p1,1);
1162 for (j = 1; j < lg(g); j++)
1163 {
1164 h = idealhermite(nf, gel(g,j));
1165 bad = lcmii(bad, gcoeff(h,1,1));
1166 }
1167 }
1168 }
1169 /* p | bad <--> p | some element occurring in cycgen[i] */
1170
1171 funits = algtobasis(nf, funits);
1172 zu = mkvec2(gel(zu,1), algtobasis(nf, gel(zu,2)));
1173
1174 for (p = *delta++; p <= bound; ) {
1175 check_prime(p,bnf,cyc,cycgen,funits,zu,bad);
1176 NEXT_PRIME_VIADIFF(p, delta);
1177 }
1178
1179 if (nbgen)
1180 {
1181 GEN f = factor(gel(cyc,1)), f1 = gel(f,1);
1182 long nbf1 = lg(f1);
1183 if (DEBUGLEVEL>1) { fprintferr(" Testing primes | h(K)\n\n"); flusherr(); }
1184 for (i=1; i<nbf1; i++)
1185 {
1186 p = itou(gel(f1,i));
1187 if (p > bound) check_prime(p,bnf,cyc,cycgen,funits,zu,bad);
1188 }
1189 }
1190 avma = av; return 1;
1191 }
1192
1193 /*******************************************************************/
1194 /* */
1195 /* RAY CLASS FIELDS: CONDUCTORS AND DISCRIMINANTS */
1196 /* */
1197 /*******************************************************************/
1198 /* Let bnr1, bnr2 be such that mod(bnr2) | mod(bnr1), compute the
1199 matrix of the surjective map Cl(bnr1) ->> Cl(bnr2) */
1200 GEN
bnrGetSurj(GEN bnr1,GEN bnr2)1201 bnrGetSurj(GEN bnr1, GEN bnr2)
1202 {
1203 long l, i;
1204 GEN M, gen = gmael(bnr1, 5, 3);
1205
1206 l = lg(gen); M = cgetg(l, t_MAT);
1207 for (i = 1; i < l; i++)
1208 gel(M,i) = isprincipalray(bnr2, gel(gen,i));
1209 return M;
1210 }
1211
1212 /* s: <gen> = Cl_f --> Cl_f2 --> 0, H subgroup of Cl_f (generators given as
1213 * HNF on [gen]). Return subgroup s(H) in Cl_f2 */
1214 static GEN
imageofgroup(GEN bnr,GEN bnr2,GEN H)1215 imageofgroup(GEN bnr, GEN bnr2, GEN H)
1216 {
1217 GEN H2, Delta = diagonal_i(gmael(bnr2,5,2)); /* SNF structure of Cl_n */
1218
1219 if (!H) return Delta;
1220 H2 = gmul(bnrGetSurj(bnr, bnr2), H);
1221 return hnf( shallowconcat(H2, Delta) ); /* s(H) in Cl_n */
1222 }
1223
1224 static GEN
args_to_bnr(GEN arg0,GEN arg1,GEN arg2,GEN * subgroup,int gen)1225 args_to_bnr(GEN arg0, GEN arg1, GEN arg2, GEN *subgroup, int gen)
1226 {
1227 GEN bnr,bnf;
1228
1229 if (typ(arg0)!=t_VEC)
1230 pari_err(talker,"neither bnf nor bnr in conductor or discray");
1231 if (!arg1) arg1 = gen_0;
1232 if (!arg2) arg2 = gen_0;
1233
1234 switch(lg(arg0))
1235 {
1236 case 7: /* bnr */
1237 bnr = arg0; (void)checkbnf(gel(bnr,1));
1238 *subgroup = arg1; break;
1239
1240 case 11: /* bnf */
1241 bnf = checkbnf(arg0);
1242 bnr = Buchray(bnf,arg1, gen? nf_INIT | nf_GEN: nf_INIT);
1243 *subgroup = arg2; break;
1244
1245 default: pari_err(talker,"neither bnf nor bnr in conductor or discray");
1246 return NULL; /* not reached */
1247 }
1248 if (!gcmp0(*subgroup))
1249 {
1250 long tx = typ(*subgroup);
1251 if (!is_matvec_t(tx))
1252 pari_err(talker,"bad subgroup in conductor or discray");
1253 }
1254 return bnr;
1255 }
1256
1257 GEN
bnrconductor(GEN arg0,GEN arg1,GEN arg2,GEN all)1258 bnrconductor(GEN arg0,GEN arg1,GEN arg2,GEN all)
1259 {
1260 long flag = all? itos(all): 0;
1261 GEN sub = arg1, bnr = args_to_bnr(arg0,arg1,arg2,&sub, flag > 0);
1262 return conductor(bnr,sub, flag);
1263 }
1264
1265 long
bnrisconductor(GEN arg0,GEN arg1,GEN arg2)1266 bnrisconductor(GEN arg0,GEN arg1,GEN arg2)
1267 {
1268 GEN sub = arg1, bnr = args_to_bnr(arg0,arg1,arg2,&sub, 0);
1269 return itos(conductor(bnr,sub,-1));
1270 }
1271
1272 static GEN
check_subgroup(GEN bnr,GEN H,GEN * clhray,int triv_is_NULL,char * s)1273 check_subgroup(GEN bnr, GEN H, GEN *clhray, int triv_is_NULL, char *s)
1274 {
1275 GEN h, D = NULL;
1276 if (H && gcmp0(H)) H = NULL;
1277 if (H)
1278 {
1279 D = diagonal_i(gmael(bnr,5,2));
1280 H = hnf(H);
1281 if (!hnfdivide(H, D)) pari_err(talker,"incorrect subgroup in %s", s);
1282 h = dethnf_i(H);
1283 if (equalii(h, *clhray)) H = NULL; else *clhray = h;
1284 }
1285 if (!H && !triv_is_NULL) H = D? D: diagonal_i(gmael(bnr,5,2));
1286 return H;
1287 }
1288
1289 /* return bnrisprincipal(bnr, (x)), assuming z = ideallog(x) */
1290 static GEN
ideallog_to_bnr(GEN bnr,GEN z)1291 ideallog_to_bnr(GEN bnr, GEN z)
1292 {
1293 GEN rayclass = gel(bnr,5), U = gel(bnr,4), divray = gel(rayclass,2);
1294 long j, l, lU, lz;
1295 int col;
1296
1297 if (lg(z) == 1) return z;
1298 col = (typ(z) == t_COL); /* else t_MAT */
1299 lz = col? lg(z): lg(z[1]);
1300 lU = lg(U);
1301 if (lz != lU)
1302 {
1303 if (lz == 1) return zerocol(lg(U[1]) - 1); /* lU != 1 */
1304 U = vecslice(U, lU-lz+1, lU-1); /* remove Cl(K) part */
1305 }
1306 z = gmul(U, z);
1307 if (col)
1308 z = vecmodii(z, divray);
1309 else
1310 {
1311 l = lg(z);
1312 for (j = 1; j < l; j++) gel(z,j) = vecmodii(gel(z,j), divray);
1313 }
1314 return z;
1315 }
1316 static GEN
bnr_log_gen_pr(GEN bnr,zlog_S * S,GEN nf,long e,long index)1317 bnr_log_gen_pr(GEN bnr, zlog_S *S, GEN nf, long e, long index)
1318 { return ideallog_to_bnr(bnr, log_gen_pr(S, index, nf, e)); }
1319 static GEN
bnr_log_gen_arch(GEN bnr,zlog_S * S,long index)1320 bnr_log_gen_arch(GEN bnr, zlog_S *S, long index)
1321 { return ideallog_to_bnr(bnr, log_gen_arch(S, index)); }
1322
1323 /* A \subset H ? Allow H = NULL = trivial subgroup */
1324 static int
contains(GEN H,GEN A)1325 contains(GEN H, GEN A)
1326 { return H? (hnf_gauss(H, A) != NULL): gcmp0(A); }
1327
1328 /* (see also Discrayrel). Given a number field bnf=bnr[1], a ray class
1329 * group structure bnr (with generators if all > 0), and a subgroup H of the
1330 * ray class group, compute the conductor of H if all=0. If all > 0, compute
1331 * furthermore the corresponding H' and output
1332 * if all = 1: [[ideal,arch],[hm,cyc,gen],H']
1333 * if all = 2: [[ideal,arch],newbnr,H']
1334 * if all < 0, answer only 1 is module is the conductor, 0 otherwise. */
1335 GEN
conductor(GEN bnr,GEN H0,long all)1336 conductor(GEN bnr, GEN H0, long all)
1337 {
1338 pari_sp av = avma;
1339 long j, k, l;
1340 GEN bnf, nf, bid, ideal, archp, clhray, bnr2, e2, e, mod, H;
1341 int iscond = 1;
1342 zlog_S S;
1343
1344 if (all > 0) checkbnrgen(bnr); else checkbnr(bnr);
1345 bnf = gel(bnr,1);
1346 bid = gel(bnr,2); init_zlog_bid(&S, bid);
1347 clhray = gmael(bnr,5,1);
1348 nf = gel(bnf,7);
1349 H = check_subgroup(bnr, H0, &clhray, 1, "conductor");
1350
1351 archp = S.archp;
1352 e = S.e; l = lg(e);
1353 e2 = cgetg(l, t_COL);
1354 for (k = 1; k < l; k++)
1355 {
1356 for (j = itos(gel(e,k)); j > 0; j--)
1357 {
1358 if (!contains(H, bnr_log_gen_pr(bnr, &S, nf, j, k))) break;
1359 if (all < 0) { avma = av; return gen_0; }
1360 iscond = 0;
1361 }
1362 gel(e2,k) = stoi(j);
1363 }
1364 l = lg(archp);
1365 for (k = 1; k < l; k++)
1366 {
1367 if (!contains(H, bnr_log_gen_arch(bnr, &S, k))) continue;
1368 if (all < 0) { avma = av; return gen_0; }
1369 archp[k] = 0;
1370 iscond = 0;
1371 }
1372 if (all < 0) { avma = av; return gen_1; }
1373 for (j = k = 1; k < l; k++)
1374 if (archp[k]) archp[j++] = archp[k];
1375 setlg(archp, j);
1376 ideal = gequal(e2, e)? gmael(bid,1,1): factorbackprime(nf, S.P, e2);
1377 mod = mkvec2(ideal, perm_to_arch(nf, archp));
1378 if (!all) return gerepilecopy(av, mod);
1379
1380 if (iscond)
1381 {
1382 bnr2 = bnr;
1383 if (!H) H = diagonal_i(gmael(bnr,5,2));
1384 }
1385 else
1386 {
1387 bnr2 = Buchray(bnf, mod, nf_INIT | nf_GEN);
1388 H = imageofgroup(bnr, bnr2, H);
1389 }
1390 return gerepilecopy(av, mkvec3(mod, (all == 1)? gel(bnr2,5): bnr2, H));
1391 }
1392
1393 /* return the norm group corresponding to the relative extension given by
1394 * polrel over bnr.bnf, assuming it is abelian and the modulus of bnr is a
1395 * multiple of the conductor */
1396 GEN
rnfnormgroup(GEN bnr,GEN polrel)1397 rnfnormgroup(GEN bnr, GEN polrel)
1398 {
1399 long i, j, reldeg, nfac, k;
1400 pari_sp av = avma;
1401 GEN bnf,index,discnf,nf,raycl,group,detgroup,fa,greldeg;
1402 GEN famo, fac, col;
1403 byteptr d = diffptr;
1404 ulong p;
1405
1406 checkbnr(bnr); bnf=gel(bnr,1); raycl=gel(bnr,5);
1407 nf=gel(bnf,7);
1408 polrel = fix_relative_pol(nf,polrel,1);
1409 if (typ(polrel)!=t_POL) pari_err(typeer,"rnfnormgroup");
1410 reldeg = degpol(polrel);
1411 /* reldeg-th powers are in norm group */
1412 greldeg = utoipos(reldeg);
1413 group = diagonal_i(FpC_red(gel(raycl,2), greldeg));
1414 for (i=1; i<lg(group); i++)
1415 if (!signe(gcoeff(group,i,i))) gcoeff(group,i,i) = greldeg;
1416 detgroup = dethnf_i(group);
1417 k = cmpiu(detgroup,reldeg);
1418 if (k < 0)
1419 pari_err(talker,"not an Abelian extension in rnfnormgroup?");
1420 if (!k) return gerepilecopy(av, group);
1421
1422 discnf = gel(nf,3);
1423 index = gel(nf,4);
1424 for (p=0 ;;)
1425 {
1426 long oldf = -1, lfa;
1427 /* If all pr are unramified and have the same residue degree, p =prod pr
1428 * and including last pr^f or p^f is the same, but the last isprincipal
1429 * is much easier! oldf is used to track this */
1430
1431 NEXT_PRIME_VIADIFF_CHECK(p,d);
1432 if (!umodiu(index, p)) continue; /* can't be treated efficiently */
1433
1434 fa = primedec(nf, utoipos(p)); lfa = lg(fa)-1;
1435 for (i=1; i<=lfa; i++)
1436 {
1437 GEN pr = gel(fa,i), pp, T, polr, modpr;
1438 long f;
1439
1440 /* primes of degree 1 are enough, and simpler */
1441 if (itos(gel(pr,4)) > 1) break;
1442
1443 modpr = nf_to_ff_init(nf, &pr, &T, &pp);
1444 polr = modprX(polrel, nf, modpr);
1445 /* if pr (probably) ramified, we have to use all (non-ram) P | pr */
1446 if (!FqX_is_squarefree(polr, T,pp)) { oldf = 0; continue; }
1447
1448 famo = FqX_factor(polr, T, pp);
1449 fac = gel(famo,1); f = degpol(gel(fac,1));
1450 nfac = lg(fac)-1;
1451 /* check decomposition of pr has Galois type */
1452 for (j=2; j<=nfac; j++)
1453 if (degpol(fac[j]) != f)
1454 pari_err(talker,"non Galois extension in rnfnormgroup");
1455 if (oldf < 0) oldf = f; else if (oldf != f) oldf = 0;
1456 if (f == reldeg) continue; /* reldeg-th powers already included */
1457
1458 if (oldf && i == lfa && !umodiu(discnf, p)) pr = utoipos(p);
1459
1460 /* pr^f = N P, P | pr, hence is in norm group */
1461 col = gmulsg(f, bnrisprincipal(bnr,pr,0));
1462 group = hnf(shallowconcat(group, col));
1463 detgroup = dethnf_i(group);
1464 k = cmpiu(detgroup,reldeg);
1465 if (k < 0) pari_err(talker,"not an Abelian extension in rnfnormgroup");
1466 if (!k) { cgiv(detgroup); return gerepileupto(av,group); }
1467 }
1468 }
1469 }
1470
1471 static GEN
liftpol(GEN pol,GEN q)1472 liftpol(GEN pol, GEN q)
1473 {
1474 long i, l = lg(pol);
1475 GEN y = cgetg(l, t_POL); y[1] = pol[1];
1476 for (i = 2; i < l; i++)
1477 gel(y,i) = lift_intern(poleval(lift_intern(gel(pol,i)), q));
1478 return y;
1479 }
1480
1481 static int
rnf_is_abelian(GEN nf,GEN pol)1482 rnf_is_abelian(GEN nf, GEN pol)
1483 {
1484 GEN modpr, pr, T, pp, ro, nfL, eq, C, z, a, sig;
1485 long i, j, l, v = varn(nf[1]);
1486 ulong p, k, ka;
1487
1488 eq = rnfequation2(nf,pol);
1489 C = shallowcopy(gel(eq,1)); setvarn(C, v);
1490 a = lift_intern(gel(eq,2)); setvarn(a, v); /* root of nf[1] */
1491 nfL = initalg_i(C, nf_PARTIALFACT, DEFAULTPREC);
1492 z = nfrootsall_and_pr(nfL, liftpol(pol, a));
1493 if (!z) return 0;
1494 ro = gel(z,1); l = lg(ro)-1;
1495 /* small groups are abelian, as are groups of prime order */
1496 if (l < 6 || uisprime(l)) return 1;
1497
1498 pr = gel(z,2);
1499 modpr = nf_to_ff_init(nfL, &pr, &T, &pp);
1500 p = itou(pp);
1501 k = umodiu(gel(eq,3), p);
1502 ka = (k * itou(nf_to_ff(nfL, a, modpr))) % p;
1503 sig= cgetg(l+1, t_VECSMALL);
1504 /* image of c = ro[1] + k a [distinguished root of C] by the l automorphisms
1505 * sig[i]: ro[1] -> ro[i] */
1506 ro = lift_intern(ro);
1507 for (i = 1; i <= l; i++)
1508 sig[i] = Fl_add(ka, itou(nf_to_ff(nfL, gel(ro,i), modpr)), p);
1509 ro = Q_primpart(ro);
1510 for (i=2; i<=l; i++) { /* start at 2, since sig[1] = identity */
1511 gel(ro,i) = ZX_to_Flx(gel(ro,i), p);
1512 for (j=2; j<i; j++)
1513 if (Flx_eval(gel(ro,j), sig[i], p)
1514 != Flx_eval(gel(ro,i), sig[j], p)) return 0;
1515 }
1516 return 1;
1517 }
1518
1519 /* Given bnf and polrel defining an abelian relative extension, compute the
1520 * corresponding conductor and congruence subgroup. Return
1521 * [[ideal,arch],[hm,cyc,gen],group] where [ideal,arch] is the conductor, and
1522 * [hm,cyc,gen] is the corresponding ray class group.
1523 * If flag != 0, check that the extension is abelian */
1524 GEN
rnfconductor(GEN bnf,GEN polrel,long flag)1525 rnfconductor(GEN bnf, GEN polrel, long flag)
1526 {
1527 pari_sp av = avma;
1528 GEN nf, module, bnr, group, p1, pol2;
1529
1530 bnf = checkbnf(bnf); nf = gel(bnf,7);
1531 if (typ(polrel) != t_POL) pari_err(typeer,"rnfconductor");
1532 p1 = unifpol(nf, polrel, t_COL);
1533 pol2 = RgX_rescale(polrel, Q_denom(p1));
1534 if (flag && !rnf_is_abelian(nf, pol2)) { avma = av; return gen_0; }
1535
1536 pol2 = fix_relative_pol(nf, pol2, 1);
1537 module = mkvec2((GEN)rnfdiscf(nf,pol2)[1],
1538 const_vec(nf_get_r1(nf), gen_1));
1539 bnr = Buchray(bnf,module,nf_INIT | nf_GEN);
1540 group = rnfnormgroup(bnr,pol2);
1541 if (!group) { avma = av; return gen_0; }
1542 return gerepileupto(av, conductor(bnr,group,1));
1543 }
1544
1545 /* Given a number field bnf=bnr[1], a ray class group structure bnr, and a
1546 * subgroup H (HNF form) of the ray class group, compute [n, r1, dk]
1547 * associated to H (cf. discrayall). If flcond = 1, abort (return gen_0) if
1548 * module is not the conductor If flrel = 0, compute only N(dk) instead of
1549 * the ideal dk proper */
1550 static GEN
Discrayrel(GEN bnr,GEN H0,long flag)1551 Discrayrel(GEN bnr, GEN H0, long flag)
1552 {
1553 pari_sp av = avma;
1554 long j, k, l, nz, flrel = flag & nf_REL, flcond = flag & nf_COND;
1555 GEN bnf, nf, bid, ideal, archp, clhray, clhss, P, e, dlk, H;
1556 zlog_S S;
1557
1558 checkbnr(bnr);
1559 bnf = gel(bnr,1);
1560 bid = gel(bnr,2); init_zlog_bid(&S, bid);
1561 clhray = gmael(bnr,5,1);
1562 nf = gel(bnf,7);
1563 ideal= gmael(bid,1,1);
1564 H0 = H = check_subgroup(bnr, H0, &clhray, 0, "bnrdiscray");
1565 archp = S.archp;
1566 P = S.P;
1567 e = S.e; l = lg(e);
1568 dlk = flrel? idealpow(nf,ideal,clhray)
1569 : powgi(dethnf_i(ideal),clhray);
1570 for (k = 1; k < l; k++)
1571 {
1572 GEN pr = gel(P,k), sum = gen_0;
1573 long ep = itos(gel(e,k));
1574 for (j = ep; j > 0; j--)
1575 {
1576 GEN z = bnr_log_gen_pr(bnr, &S, nf, j, k);
1577 H = hnf(shallowconcat(H, z));
1578 clhss = dethnf_i(H);
1579 if (flcond && j==ep && equalii(clhss,clhray)) { avma = av; return gen_0; }
1580 if (is_pm1(clhss)) { sum = addis(sum, j); break; }
1581 sum = addii(sum, clhss);
1582 }
1583 dlk = flrel? idealdivpowprime(nf, dlk, pr, sum)
1584 : diviiexact(dlk, powgi(pr_norm(pr),sum));
1585 }
1586 l = lg(archp); nz = nf_get_r1(nf) - (l-1);
1587 for (k = 1; k < l; k++)
1588 {
1589 if (!contains(H0, bnr_log_gen_arch(bnr, &S, k))) continue;
1590 if (flcond) { avma = av; return gen_0; }
1591 nz++;
1592 }
1593 return gerepilecopy(av, mkvec3(clhray, stoi(nz), dlk));
1594 }
1595
1596 static GEN
Discrayabs(GEN bnr,GEN subgroup,long flag)1597 Discrayabs(GEN bnr, GEN subgroup, long flag)
1598 {
1599 pari_sp av = avma;
1600 long clhray, n, R1;
1601 GEN z, p1, D, dk, nf, dkabs;
1602
1603 D = Discrayrel(bnr, subgroup, flag);
1604 if ((flag & nf_REL) || D == gen_0) return D;
1605
1606 nf = checknf(bnr);
1607 dkabs = absi(gel(nf,3));
1608 clhray = itos(gel(D,1)); p1 = powiu(dkabs, clhray);
1609 n = clhray * degpol(nf[1]);
1610 R1= clhray * itos(gel(D,2));
1611 dk = gel(D,3);
1612 if (((n-R1)&3) == 2) dk = negi(dk); /* (2r2) mod 4 = 2 : r2(relext) is odd */
1613 z = cgetg(4,t_VEC);
1614 gel(z,1) = utoipos(n);
1615 gel(z,2) = stoi(R1);
1616 gel(z,3) = mulii(dk,p1); return gerepileupto(av, z);
1617 }
1618
1619 GEN
bnrdisc0(GEN arg0,GEN arg1,GEN arg2,long flag)1620 bnrdisc0(GEN arg0, GEN arg1, GEN arg2, long flag)
1621 {
1622 GEN H, bnr = args_to_bnr(arg0,arg1,arg2,&H, 0);
1623 return Discrayabs(bnr,H,flag);
1624 }
1625 GEN
discrayrel(GEN bnr,GEN H)1626 discrayrel(GEN bnr, GEN H)
1627 { return Discrayrel(bnr,H,nf_REL); }
1628 GEN
discrayrelcond(GEN bnr,GEN H)1629 discrayrelcond(GEN bnr, GEN H)
1630 { return Discrayrel(bnr,H,nf_REL | nf_COND); }
1631 GEN
discrayabs(GEN bnr,GEN H)1632 discrayabs(GEN bnr, GEN H)
1633 { return Discrayabs(bnr,H,0); }
1634 GEN
discrayabscond(GEN bnr,GEN H)1635 discrayabscond(GEN bnr, GEN H)
1636 { return Discrayabs(bnr,H,nf_COND); }
1637
1638 /* chi character of abelian G: chi[i] = chi(z_i), where G = \oplus Z/cyc[i] z_i.
1639 * Return Ker chi [ NULL = trivial subgroup of G ] */
1640 static GEN
KerChar(GEN chi,GEN cyc)1641 KerChar(GEN chi, GEN cyc)
1642 {
1643 long i, l = lg(cyc);
1644 GEN m, U, d1;
1645
1646 if (lg(chi) != l) pari_err(talker,"incorrect character length in KerChar");
1647 if (l == 1) return NULL; /* trivial subgroup */
1648 d1 = gel(cyc,1); m = cgetg(l+1,t_MAT);
1649 for (i=1; i<l; i++)
1650 {
1651 if (typ(chi[i]) != t_INT) pari_err(typeer,"conductorofchar");
1652 gel(m,i) = mkcol(mulii(gel(chi,i), diviiexact(d1, gel(cyc,i))));
1653 }
1654 gel(m,i) = mkcol(d1);
1655 (void)hnfall_i(m, &U, 1);
1656 for (i = 1; i < l; i++) setlg(U[i], l);
1657 setlg(U,l); return U;
1658 }
1659
1660 /* Given a number field bnf=bnr[1], a ray class group structure bnr and a
1661 * vector chi representing a character on the generators bnr[2][3], compute
1662 * the conductor of chi. */
1663 GEN
bnrconductorofchar(GEN bnr,GEN chi)1664 bnrconductorofchar(GEN bnr, GEN chi)
1665 {
1666 pari_sp av = avma; checkbnr(bnr);
1667 return gerepileupto(av, conductor(bnr, KerChar(chi, gmael(bnr,5,2)), 0));
1668 }
1669
1670 /* t = [bid,U], h = #Cl(K) */
1671 static GEN
get_classno(GEN t,GEN h)1672 get_classno(GEN t, GEN h)
1673 {
1674 GEN bid = gel(t,1), cyc = gmael(bid,2,2);
1675 GEN m = shallowconcat(gel(t,2), diagonal_i(cyc));
1676 return mulii(h, dethnf_i(hnf(m)));
1677 }
1678
1679 static void
chk_listBU(GEN L,char * s)1680 chk_listBU(GEN L, char *s) {
1681 if (typ(L) != t_VEC) pari_err(typeer,s);
1682 if (lg(L) > 1) {
1683 GEN z = gel(L,1);
1684 if (typ(z) != t_VEC) pari_err(typeer, s);
1685 if (lg(z) == 1) return;
1686 z = gel(z,1); /* [bid,U] */
1687 if (typ(z) != t_VEC || lg(z) != 3) pari_err(typeer, s);
1688 checkbid(gel(z,1));
1689 }
1690 }
1691
1692 /* Given lists of [zidealstarinit, unit ideallogs], return lists of ray class
1693 * numbers */
1694 GEN
bnrclassnolist(GEN bnf,GEN L)1695 bnrclassnolist(GEN bnf,GEN L)
1696 {
1697 pari_sp av = avma;
1698 long i, j, lz, l = lg(L);
1699 GEN v, z, V, h;
1700
1701 chk_listBU(L, "bnrclassnolist");
1702 if (l == 1) return cgetg(1, t_VEC);
1703 bnf = checkbnf(bnf); h = gmael3(bnf,8,1,1);
1704 V = cgetg(l,t_VEC);
1705 for (i = 1; i < l; i++)
1706 {
1707 z = gel(L,i); lz = lg(z);
1708 gel(V,i) = v = cgetg(lz,t_VEC);
1709 for (j=1; j<lz; j++) gel(v,j) = get_classno(gel(z,j), h);
1710 }
1711 return gerepilecopy(av, V);
1712 }
1713
1714 static GEN
Lbnrclassno(GEN L,GEN fac)1715 Lbnrclassno(GEN L, GEN fac)
1716 {
1717 long i, l = lg(L);
1718 for (i=1; i<l; i++)
1719 if (gequal(gmael(L,i,1),fac)) return gmael(L,i,2);
1720 pari_err(bugparier,"Lbnrclassno");
1721 return NULL; /* not reached */
1722 }
1723
1724 /* returns the first index i<=n such that x=v[i] if it exits, 0 otherwise */
1725 long
isinvector(GEN v,GEN x)1726 isinvector(GEN v, GEN x)
1727 {
1728 long i, l = lg(v);
1729 for (i = 1; i < l; i++)
1730 if (gequal(gel(v,i), x)) return i;
1731 return 0;
1732 }
1733
1734 static GEN
factordivexact(GEN fa1,GEN fa2)1735 factordivexact(GEN fa1,GEN fa2)
1736 {
1737 long i, j, k, c, l;
1738 GEN P, E, P1, E1, P2, E2, p1;
1739
1740 P1 = gel(fa1,1); E1 = gel(fa1,2); l = lg(P1);
1741 P2 = gel(fa2,1); E2 = gel(fa2,2);
1742 P = cgetg(l,t_COL);
1743 E = cgetg(l,t_COL);
1744 for (c = i = 1; i < l; i++)
1745 {
1746 j = isinvector(P2,gel(P1,i));
1747 if (!j) { P[c] = P1[i]; E[c] = E1[i]; c++; }
1748 else
1749 {
1750 p1 = subii(gel(E1,i), gel(E2,j)); k = signe(p1);
1751 if (k < 0) pari_err(talker,"factordivexact is not exact!");
1752 if (k > 0) { P[c] = P1[i]; gel(E,c) = p1; c++; }
1753 }
1754 }
1755 setlg(P, c);
1756 setlg(E, c); return mkmat2(P, E);
1757 }
1758 /* remove index k */
1759 static GEN
factorsplice(GEN fa,long k)1760 factorsplice(GEN fa, long k)
1761 {
1762 GEN p = gel(fa,1), e = gel(fa,2), P, E;
1763 long i, l = lg(p) - 1;
1764 P = cgetg(l, typ(p));
1765 E = cgetg(l, typ(e));
1766 for (i=1; i<k; i++) { P[i] = p[i]; E[i] = e[i]; }
1767 p++; e++;
1768 for ( ; i<l; i++) { P[i] = p[i]; E[i] = e[i]; }
1769 return mkmat2(P,E);
1770 }
1771 static GEN
factorpow(GEN fa,long n)1772 factorpow(GEN fa, long n)
1773 {
1774 if (!n) return trivfact();
1775 return mkmat2(gel(fa,1), gmulsg(n, gel(fa,2)));
1776 }
1777 static GEN
factormul(GEN fa1,GEN fa2)1778 factormul(GEN fa1,GEN fa2)
1779 {
1780 GEN p, pnew, e, enew, v, P, y = concat_factor(fa1,fa2);
1781 long i, c, lx;
1782
1783 p = gel(y,1); v = sindexsort(p); lx = lg(p);
1784 e = gel(y,2);
1785 pnew = vecpermute(p, v);
1786 enew = vecpermute(e, v);
1787 P = gen_0; c = 0;
1788 for (i=1; i<lx; i++)
1789 {
1790 if (gequal(gel(pnew,i),P))
1791 gel(e,c) = addii(gel(e,c),gel(enew,i));
1792 else
1793 {
1794 c++; P = gel(pnew,i);
1795 gel(p,c) = P;
1796 e[c] = enew[i];
1797 }
1798 }
1799 setlg(p, c+1);
1800 setlg(e, c+1); return y;
1801 }
1802
1803
1804 static long
get_nz(GEN bnf,GEN ideal,GEN arch,long clhray)1805 get_nz(GEN bnf, GEN ideal, GEN arch, long clhray)
1806 {
1807 GEN arch2 = shallowcopy(arch), mod = mkvec2(ideal, arch2);
1808 long nz = 0, l = lg(arch), k, clhss;
1809 for (k = 1; k < l; k++)
1810 { /* FIXME: this is wasteful. Use the same algorithm as conductor */
1811 if (signe(arch[k]))
1812 {
1813 gel(arch2,k) = gen_0; clhss = itos(bnrclassno(bnf,mod));
1814 gel(arch2,k) = gen_1;
1815 if (clhss == clhray) return -1;
1816 }
1817 else nz++;
1818 }
1819 return nz;
1820 }
1821
1822 static GEN
get_NR1D(long Nf,long clhray,long degk,long nz,GEN fadkabs,GEN idealrel)1823 get_NR1D(long Nf, long clhray, long degk, long nz, GEN fadkabs, GEN idealrel)
1824 {
1825 long n, R1;
1826 GEN dlk;
1827 if (nz < 0) return NULL;
1828 n = clhray * degk;
1829 R1 = clhray * nz;
1830 dlk = factordivexact(factorpow(factor(utoipos(Nf)),clhray), idealrel);
1831 /* r2 odd, set dlk = -dlk */
1832 if (((n-R1)&3)==2) dlk = factormul(to_famat_all(gen_m1,gen_1), dlk);
1833 return mkvec3(utoipos(n),
1834 stoi(R1),
1835 factormul(dlk,factorpow(fadkabs,clhray)));
1836 }
1837
1838 /* t = [bid,U], h = #Cl(K) */
1839 static GEN
get_discdata(GEN t,GEN h)1840 get_discdata(GEN t, GEN h)
1841 {
1842 GEN bid = gel(t,1), fa = gel(bid,3);
1843 return mkvec3(mkmat2(gel(fa,1), vec_to_vecsmall(gel(fa,2))),
1844 (GEN)itou(get_classno(t, h)),
1845 gel(bid,1));
1846 }
1847 typedef struct _disc_data {
1848 long degk;
1849 GEN bnf, fadk, idealrelinit, V;
1850 } disc_data;
1851
1852 static GEN
get_discray(disc_data * D,GEN V,GEN x,GEN z,long N)1853 get_discray(disc_data *D, GEN V, GEN x, GEN z, long N)
1854 {
1855 GEN idealrel = D->idealrelinit;
1856 GEN mod = gel(z,3), Fa = gel(z,1);
1857 GEN P = gel(Fa,1), E = gel(Fa,2);
1858 long k, nz, clhray = z[2], lP = lg(P);
1859 for (k=1; k<lP; k++)
1860 {
1861 GEN pr = gel(P,k), p = gel(pr,1);
1862 long e, ep = E[k], f = itos(gel(pr,4));
1863 long S = 0, norm = N, Npr, clhss;
1864 Npr = itos(powiu(p,f));
1865 for (e=1; e<=ep; e++)
1866 {
1867 GEN fad;
1868 if (e < ep) { E[k] = ep-e; fad = Fa; }
1869 else fad = factorsplice(Fa, k);
1870 norm /= Npr;
1871 clhss = (long)Lbnrclassno(gel(V,norm), fad);
1872 if (e==1 && clhss==clhray) { E[k] = ep; return cgetg(1, t_VEC); }
1873 if (clhss == 1) { S += ep-e+1; break; }
1874 S += clhss;
1875 }
1876 E[k] = ep;
1877 idealrel = factormul(idealrel, to_famat_all(p, utoi(f * S)));
1878 }
1879 nz = get_nz(D->bnf, gel(mod,1), gel(mod,2), clhray);
1880 return get_NR1D(N, clhray, D->degk, nz, D->fadk, idealrel);
1881 }
1882
1883 /* Given a list of bids and associated unit log matrices, return the
1884 * list of discrayabs. Only keep moduli which are conductors. */
1885 GEN
discrayabslist(GEN bnf,GEN L)1886 discrayabslist(GEN bnf, GEN L)
1887 {
1888 pari_sp av = avma;
1889 long i, j, lz, l = lg(L);
1890 GEN nf, v, z, V, D, d, h;
1891 disc_data ID;
1892
1893 chk_listBU(L, "discrayabslist");
1894 if (l == 1) return cgetg(1, t_VEC);
1895 ID.bnf = bnf = checkbnf(bnf);
1896 nf = gel(bnf,7);
1897 h = gmael3(bnf,8,1,1);
1898 ID.degk = degpol(nf[1]);
1899 ID.fadk = factor(absi(gel(nf,3)));
1900 ID.idealrelinit = trivfact();
1901 V = cgetg(l, t_VEC);
1902 D = cgetg(l, t_VEC);
1903 for (i = 1; i < l; i++)
1904 {
1905 z = gel(L,i); lz = lg(z);
1906 gel(V,i) = v = cgetg(lz,t_VEC);
1907 gel(D,i) = d = cgetg(lz,t_VEC);
1908 for (j=1; j<lz; j++) {
1909 gel(d,j) = get_discdata(gel(z,j), h);
1910 gel(v,j) = get_discray(&ID, D, gel(z,j), gel(d,j), i);
1911 }
1912 }
1913 return gerepilecopy(av, V);
1914 }
1915
1916 /* BIG VECTOR:
1917 * Interface: a container v whose length is arbitrary (< 2^30), bigel(v,i)
1918 * refers to the i-th component. It is an lvalue.
1919 *
1920 * Implementation: a vector v whose components have exactly 2^LGVINT entries
1921 * but for the last one which is allowed to be shorter. v[i][j]
1922 * (where j<=2^LGVINT) is understood as component number I = (i-1)*2^LGVINT+j
1923 * in a unique huge vector V. */
1924
1925 #define SHLGVINT 15
1926 #define LGVINT (1L << SHLGVINT)
1927 #define vext0(i) ((((i)-1)>>SHLGVINT)+1)
1928 #define vext1(i) ((i)&(LGVINT-1))
1929 #define bigel(v,i) gmael((v), vext0(i), vext1(i))
1930
1931 /* allocate an extended vector (t_VEC of t_VEC) for N _true_ components */
1932 static GEN
bigcgetvec(long N)1933 bigcgetvec(long N)
1934 {
1935 long i, nv = vext0(N);
1936 GEN v = cgetg(nv+1,t_VEC);
1937 for (i=1; i<nv; i++) gel(v,i) = cgetg(LGVINT+1,t_VEC);
1938 gel(v,nv) = cgetg(vext1(N)+1,t_VEC); return v;
1939 }
1940
1941 static GEN
zsimp(GEN bid,GEN embunit)1942 zsimp(GEN bid, GEN embunit)
1943 {
1944 GEN empty = cgetg(1, t_VECSMALL);
1945 return mkvec4(mkmat2(empty,empty), gmael(bid,2,2),
1946 gel(bid,5), embunit);
1947 }
1948
1949 static GEN
zsimpjoin(GEN b,GEN bid,GEN embunit,long prcode,long e)1950 zsimpjoin(GEN b, GEN bid, GEN embunit, long prcode, long e)
1951 {
1952 long i, l1, l2, nbgen, c;
1953 pari_sp av = avma;
1954 GEN fa, U, U1, U2, cyc1, cyc2, u1u2, D;
1955
1956 fa = gel(b,1);
1957 U1 = gel(b,3); cyc1 = gel(b,2); l1 = lg(cyc1);
1958 U2 = gel(bid,5); cyc2 = gmael(bid,2,2); l2 = lg(cyc2);
1959 nbgen = l1+l2-2;
1960 if (nbgen)
1961 {
1962 u1u2 = matsnf0(diagonal_i(shallowconcat(cyc1,cyc2)), 1 | 4); /* all && clean */
1963 U = gel(u1u2,1);
1964 D = gel(u1u2,3);
1965 U = shallowconcat(
1966 l1==1 ? zeromat(nbgen, lg(U1)-1): gmul(vecslice(U, 1, l1-1), U1),
1967 l1>nbgen? zeromat(nbgen, lg(U2)-1): gmul(vecslice(U, l1, nbgen), U2)
1968 );
1969 }
1970 else
1971 {
1972 c = lg(U1)+lg(U2)-1; U = cgetg(c,t_MAT);
1973 for (i=1; i<c; i++) gel(U,i) = cgetg(1,t_COL);
1974 D = cgetg(1,t_MAT);
1975 }
1976 return gerepilecopy(av, mkvec4(
1977 mkmat2(vecsmall_append(gel(fa,1), prcode),
1978 vecsmall_append(gel(fa,2), e)),
1979 mattodiagonal_i(D),
1980 U,
1981 vconcat(gel(b,4),embunit)
1982 ));
1983 }
1984
1985 static GEN
bnrclassnointern(GEN B,GEN h)1986 bnrclassnointern(GEN B, GEN h)
1987 {
1988 long lx = lg(B), j;
1989 GEN b, m, qm, L = cgetg(lx,t_VEC);
1990 for (j=1; j<lx; j++)
1991 {
1992 b = gel(B,j); qm = gmul(gel(b,3),gel(b,4));
1993 m = shallowconcat(qm, diagonal_i(gel(b,2)));
1994 gel(L,j) = mkvec2(gel(b,1),
1995 mkvecsmall( itou( mulii(h, dethnf_i(hnf(m))) ) ));
1996 }
1997 return L;
1998 }
1999
2000 static GEN
bnrclassnointernarch(GEN B,GEN h,GEN matU)2001 bnrclassnointernarch(GEN B, GEN h, GEN matU)
2002 {
2003 long lx, nc, k, kk, j, r1, jj, nba, nbarch;
2004 GEN _2, b, qm, L, cyc, m, H, mm, rowsel;
2005
2006 if (!matU) return bnrclassnointern(B,h);
2007 lx = lg(B); if (lx == 1) return B;
2008
2009 r1 = lg(matU[1])-1; _2 = const_vec(r1, gen_2);
2010 L = cgetg(lx,t_VEC); nbarch = 1<<r1;
2011 for (j=1; j<lx; j++)
2012 {
2013 b = gel(B,j); qm = gmul(gel(b,3),gel(b,4));
2014 cyc = gel(b,2); nc = lg(cyc)-1;
2015 /* [ qm cyc 0 ]
2016 * [ matU 0 2 ] */
2017 m = shallowconcat(vconcat(qm, matU),
2018 diagonal_i(shallowconcat(cyc, _2)));
2019 m = hnf(m); mm = shallowcopy(m);
2020 H = cgetg(nbarch+1,t_VECSMALL);
2021 rowsel = cgetg(nc+r1+1,t_VECSMALL);
2022 for (k = 0; k < nbarch; k++)
2023 {
2024 nba = nc+1;
2025 for (kk=k,jj=1; jj<=r1; jj++,kk>>=1)
2026 if (kk&1) rowsel[nba++] = nc + jj;
2027 setlg(rowsel, nba);
2028 rowselect_p(m, mm, rowsel, nc+1);
2029 H[k+1] = itou( mulii(h, dethnf_i(hnf(mm))) );
2030 }
2031 gel(L,j) = mkvec2(gel(b,1), H);
2032 }
2033 return L;
2034 }
2035
2036 GEN
decodemodule(GEN nf,GEN fa)2037 decodemodule(GEN nf, GEN fa)
2038 {
2039 long n, nn, k;
2040 pari_sp av = avma;
2041 GEN G, E, id, pr;
2042
2043 nf = checknf(nf);
2044 if (typ(fa)!=t_MAT || lg(fa)!=3)
2045 pari_err(talker,"not a factorisation in decodemodule");
2046 n = degpol(nf[1]); nn = n*n; id = NULL;
2047 G = gel(fa,1);
2048 E = gel(fa,2);
2049 for (k=1; k<lg(G); k++)
2050 {
2051 long code = itos(gel(G,k)), p = code / nn, j = (code%n)+1;
2052 GEN P = primedec(nf, utoipos(p)), e = gel(E,k);
2053 if (lg(P) <= j) pari_err(talker, "incorrect hash code in decodemodule");
2054 pr = gel(P,j);
2055 id = id? idealmulpowprime(nf,id, pr,e)
2056 : idealpow(nf, pr,e);
2057 }
2058 if (!id) { avma = av; return matid(n); }
2059 return gerepileupto(av,id);
2060 }
2061
2062 /* List of ray class fields. Do all from scratch, bound < 2^30. No subgroups.
2063 *
2064 * Output: a "big vector" V (cf bigcgetvec). V[k] is a vector indexed by
2065 * the ideals of norm k. Given such an ideal m, the component is as follows:
2066 *
2067 * + if arch = NULL, run through all possible archimedean parts; archs are
2068 * ordered using inverse lexicographic order, [0,..,0], [1,0,..,0], [0,1,..,0],
2069 * Component is [m,V] where V is a vector with 2^r1 entries, giving for each
2070 * arch the triple [N,R1,D], with N, R1, D as in discrayabs; D is in factored
2071 * form.
2072 *
2073 * + otherwise [m,N,R1,D] */
2074 GEN
discrayabslistarch(GEN bnf,GEN arch,long bound)2075 discrayabslistarch(GEN bnf, GEN arch, long bound)
2076 {
2077 byteptr dif = diffptr + 1;
2078 int allarch = (arch==NULL), flbou = 0;
2079 long degk, i, j, k, sqbou, l, nba, nbarch, ii, r1, c;
2080 pari_sp av0 = avma, av, av1, lim;
2081 GEN nf, p, Z, fa, ideal, bidp, matarchunit, Disc, U, sgnU, EMPTY;
2082 GEN res, embunit, h, Ray, discall, idealrel, idealrelinit, fadkabs;
2083
2084 if (bound <= 0) pari_err(talker,"non-positive bound in Discrayabslist");
2085 res = discall = NULL; /* -Wall */
2086
2087 bnf = checkbnf(bnf);
2088 nf = gel(bnf,7); r1 = nf_get_r1(nf);
2089 degk = degpol(nf[1]);
2090 fadkabs = factor(absi(gel(nf,3)));
2091 h = gmael3(bnf,8,1,1);
2092 U = init_units(bnf);
2093 sgnU = zsignunits(bnf, NULL, 1);
2094
2095 if (allarch) arch = const_vec(r1, gen_1);
2096 bidp = Idealstar(nf, mkvec2(gen_1, arch), 0);
2097 if (allarch) {
2098 matarchunit = zlog_units(nf, U, sgnU, bidp);
2099 bidp = Idealstar(nf,matid(degk),0);
2100 if (r1>15) pari_err(talker,"r1>15 in discrayabslistarch");
2101 nba = r1;
2102 } else {
2103 matarchunit = (GEN)NULL;
2104 for (nba=0,i=1; i<=r1; i++) if (signe(arch[i])) nba++;
2105 }
2106
2107 /* what follows was rewritten from Ideallist */
2108 p = utoipos(2);
2109 av = avma; lim = stack_lim(av,1);
2110 sqbou = (long)sqrt((double)bound) + 1;
2111 Z = bigcgetvec(bound);
2112 for (i=2; i<=bound; i++) bigel(Z,i) = cgetg(1,t_VEC);
2113 embunit = zlog_units(nf, U, sgnU, bidp);
2114 bigel(Z,1) = mkvec(zsimp(bidp,embunit));
2115 if (DEBUGLEVEL>1) fprintferr("Starting zidealstarunits computations\n");
2116 maxprime_check((ulong)bound);
2117 /* The goal is to compute Ray (lists of bnrclassno). Z contains "zsimps",
2118 * simplified zidealstarinit, from which bnrclassno is easy to compute.
2119 * Once p > sqbou, delete Z[i] for i > sqbou and compute directly Ray */
2120 Ray = Z;
2121 while (p[2] <= bound)
2122 {
2123 if (!flbou && p[2] > sqbou)
2124 {
2125 GEN z;
2126 flbou = 1;
2127 if (DEBUGLEVEL>1) fprintferr("\nStarting bnrclassno computations\n");
2128 Z = gerepilecopy(av,Z); av1 = avma;
2129 Ray = bigcgetvec(bound);
2130 for (i=1; i<=bound; i++)
2131 bigel(Ray,i) = bnrclassnointernarch(bigel(Z,i),h,matarchunit);
2132 Ray = gerepilecopy(av1,Ray);
2133 z = bigcgetvec(sqbou);
2134 for (i=1; i<=sqbou; i++) bigel(z,i) = bigel(Z,i);
2135 Z = z;
2136 }
2137 fa = primedec(nf,p);
2138 for (j=1; j<lg(fa); j++)
2139 {
2140 GEN pr = gel(fa,j);
2141 long prcode, q, f = itos(gel(pr,4)), Q = itos_or_0(powiu(p,f));
2142 if (!Q || Q > bound) continue;
2143
2144 /* p, f-1, j-1 as a single integer in "base degk" (f,j <= degk)*/
2145 prcode = (p[2]*degk + f-1)*degk + j-1;
2146 q = Q; ideal = pr;
2147 for (l=1;; l++) /* Q <= bound */
2148 {
2149 bidp = Idealstar(nf,ideal,0);
2150 embunit = zlog_units_noarch(nf, U, bidp);
2151 for (i=Q; i<=bound; i+=Q)
2152 {
2153 GEN pz, p2, p1 = bigel(Z,i/Q);
2154 long lz = lg(p1);
2155 if (lz == 1) continue;
2156
2157 p2 = cgetg(lz,t_VEC); c = 0;
2158 for (k=1; k<lz; k++)
2159 {
2160 GEN z = gel(p1,k), v = gmael(z,1,1); /* primes in zsimp's fact. */
2161 long lv = lg(v);
2162 /* If z has a power of pr in its modulus, skip it */
2163 if (Q != i && lv > 1 && v[lv-1] == prcode) break;
2164 gel(p2,++c) = zsimpjoin(z,bidp,embunit,prcode,l);
2165 }
2166
2167 setlg(p2, c+1);
2168 pz = bigel(Ray,i);
2169 if (flbou) p2 = bnrclassnointernarch(p2,h,matarchunit);
2170 if (lg(pz) > 1) p2 = shallowconcat(pz,p2);
2171 bigel(Ray,i) = p2;
2172 }
2173 Q = itos_or_0( mulss(Q, q) );
2174 if (!Q || Q > bound) break;
2175
2176 ideal = idealmul(nf,ideal,pr);
2177 }
2178 }
2179 if (low_stack(lim, stack_lim(av,1)))
2180 {
2181 if(DEBUGMEM>1) pari_warn(warnmem,"[1]: discrayabslistarch");
2182 gerepileall(av, flbou? 2: 1, &Z, &Ray);
2183 }
2184 NEXT_PRIME_VIADIFF(p[2], dif);
2185 }
2186 if (!flbou) /* occurs iff bound = 1,2,4 */
2187 {
2188 if (DEBUGLEVEL>1) fprintferr("\nStarting bnrclassno computations\n");
2189 Ray = bigcgetvec(bound);
2190 for (i=1; i<=bound; i++)
2191 bigel(Ray,i) = bnrclassnointernarch(bigel(Z,i),h,matarchunit);
2192 }
2193 Ray = gerepilecopy(av, Ray);
2194
2195 if (DEBUGLEVEL>1) fprintferr("Starting discrayabs computations\n");
2196 if (allarch) nbarch = 1<<r1;
2197 else
2198 {
2199 nbarch = 1;
2200 discall = cgetg(2,t_VEC);
2201 }
2202 EMPTY = mkvec3(gen_0,gen_0,gen_0);
2203 idealrelinit = trivfact();
2204 av1 = avma; lim = stack_lim(av1,1);
2205 Disc = bigcgetvec(bound);
2206 for (i=1; i<=bound; i++) bigel(Disc,i) = cgetg(1,t_VEC);
2207 for (ii=1; ii<=bound; ii++)
2208 {
2209 GEN sous, sousdisc;
2210 long ls;
2211 i = ii;
2212 sous = bigel(Ray,i);
2213 ls = lg(sous); bigel(Disc,ii) = sousdisc = cgetg(ls,t_VEC);
2214 for (j=1; j<ls; j++)
2215 {
2216 GEN b = gel(sous,j), clhrayall = gel(b,2), Fa = gel(b,1);
2217 GEN P = gel(Fa,1), E = gel(Fa,2);
2218 long lP = lg(P), karch;
2219
2220 if (allarch) discall = cgetg(nbarch+1,t_VEC);
2221 for (karch=0; karch<nbarch; karch++)
2222 {
2223 long nz, clhray = clhrayall[karch+1];
2224 if (allarch)
2225 {
2226 long ka, k2;
2227 nba = 0;
2228 for (ka=karch,k=1; k<=r1; k++,ka>>=1)
2229 if (ka & 1) nba++;
2230 for (k2=1,k=1; k<=r1; k++,k2<<=1)
2231 if (karch&k2 && clhrayall[karch-k2+1] == clhray)
2232 { res = EMPTY; goto STORE; }
2233 }
2234 idealrel = idealrelinit;
2235 for (k=1; k<lP; k++) /* cf get_discray */
2236 {
2237 long e, ep = E[k], pf = P[k] / degk, f = (pf%degk) + 1;
2238 long S = 0, normi = i, Npr, clhss;
2239 p = utoipos(pf / degk);
2240 Npr = itos(powiu(p,f));
2241 for (e=1; e<=ep; e++)
2242 {
2243 GEN fad;
2244 if (e < ep) { E[k] = ep-e; fad = Fa; }
2245 else fad = factorsplice(Fa, k);
2246 normi /= Npr;
2247 clhss = Lbnrclassno(bigel(Ray,normi),fad)[karch+1];
2248 if (e==1 && clhss==clhray) { E[k] = ep; res = EMPTY; goto STORE; }
2249 if (clhss == 1) { S += ep-e+1; break; }
2250 S += clhss;
2251 }
2252 E[k] = ep;
2253 idealrel = factormul(idealrel, to_famat_all(p, utoi(f * S)));
2254 }
2255 if (!allarch && nba)
2256 nz = get_nz(bnf, decodemodule(nf,Fa), arch, clhray);
2257 else
2258 nz = r1 - nba;
2259 res = get_NR1D(i, clhray, degk, nz, fadkabs, idealrel);
2260 STORE: gel(discall,karch+1) = res;
2261 }
2262 res = allarch? mkvec2(Fa, discall)
2263 : mkvec4(Fa, gel(res,1), gel(res,2), gel(res,3));
2264 gel(sousdisc,j) = res;
2265 if (low_stack(lim, stack_lim(av1,1)))
2266 {
2267 if(DEBUGMEM>1) pari_warn(warnmem,"[2]: discrayabslistarch");
2268 Disc = gerepilecopy(av1, Disc);
2269 }
2270 }
2271 }
2272 return gerepilecopy(av0, Disc);
2273 }
2274 GEN
discrayabslistlong(GEN bnf,long bound)2275 discrayabslistlong(GEN bnf, long bound) {
2276 GEN nf = checknf(bnf);
2277 long r1 = nf_get_r1(nf);
2278 return discrayabslistarch(bnf,zerovec(r1),bound);
2279 }
2280
2281 static GEN
subgroupcond(GEN bnr,GEN indexbound)2282 subgroupcond(GEN bnr, GEN indexbound)
2283 {
2284 pari_sp av = avma;
2285 long i, k, l, le, la;
2286 GEN e, li, p1, lidet, perm, L, nf = checknf(bnr);
2287 zlog_S S;
2288
2289 checkbnr(bnr); init_zlog_bid(&S, gel(bnr,2));
2290 e = S.e; le = lg(e); la = lg(S.archp);
2291 L = cgetg(le + la - 1, t_VEC);
2292 i = 1;
2293 for (k = 1; k < le; k++)
2294 gel(L,i++) = bnr_log_gen_pr(bnr, &S, nf, itos(gel(e,k)), k);
2295 for (k = 1; k < la; k++)
2296 gel(L,i++) = bnr_log_gen_arch(bnr, &S, k);
2297 li = subgroupcondlist(gmael(bnr,5,2), indexbound, L);
2298 l = lg(li);
2299 /* sort by increasing index */
2300 lidet = cgetg(l,t_VEC);
2301 for (i=1; i<l; i++) gel(lidet,i) = dethnf_i(gel(li,i));
2302 perm = sindexsort(lidet); p1 = li; li = cgetg(l,t_VEC);
2303 for (i=1; i<l; i++) li[i] = p1[perm[l-i]];
2304 return gerepilecopy(av,li);
2305 }
2306
2307 GEN
subgrouplist0(GEN bnr,GEN indexbound,long all)2308 subgrouplist0(GEN bnr, GEN indexbound, long all)
2309 {
2310 if (typ(bnr)!=t_VEC) pari_err(typeer,"subgrouplist");
2311 if (lg(bnr)!=1 && typ(bnr[1])!=t_INT)
2312 {
2313 if (!all) return subgroupcond(bnr,indexbound);
2314 checkbnr(bnr); bnr = gmael(bnr,5,2);
2315 }
2316 return subgrouplist(bnr,indexbound);
2317 }
2318
2319 GEN
bnrdisclist0(GEN bnf,GEN L,GEN arch)2320 bnrdisclist0(GEN bnf, GEN L, GEN arch)
2321 {
2322 if (typ(L)!=t_INT) return discrayabslist(bnf,L);
2323 return discrayabslistarch(bnf,arch,itos(L));
2324 }
2325