OpenXM_contrib/pari-2.2/src/basemath/buch3.c - diff

Return to buch3.c CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath

Diff for /OpenXM_contrib/pari-2.2/src/basemath/Attic/buch3.c between version 1.1 and 1.2

version 1.1, 2001/10/02 11:17:03

version 1.2, 2002/09/11 07:26:50

Line 21 Foundation, Inc., 59 Temple Place - Suite 330, Boston,

#include "pari.h"

#include "parinf.h"

extern GEN concatsp3(GEN x, GEN y, GEN z);

extern void testprimes(GEN bnf, long bound);

extern GEN Fp_PHlog(GEN a, GEN g, GEN p, GEN ord);

extern GEN FqX_factor(GEN x, GEN T, GEN p);

extern GEN anti_unif_mod_f(GEN nf, GEN pr, GEN sqf);

extern GEN arch_mul(GEN x, GEN y);

extern GEN check_and_build_cycgen(GEN bnf);

extern GEN colreducemodHNF(GEN x, GEN y, GEN *Q);

extern GEN detcyc(GEN cyc);

extern GEN famat_reduce(GEN fa);

extern GEN famat_to_nf_modideal_coprime(GEN nf, GEN g, GEN e, GEN id);

extern GEN famat_to_nf_modidele(GEN nf, GEN g, GEN e, GEN bid);

extern GEN gmul_mat_smallvec(GEN x, GEN y);

extern GEN idealaddtoone_i(GEN nf, GEN x, GEN y);

extern GEN idealprodprime(GEN nf, GEN L);

extern GEN ideleaddone_aux(GEN nf,GEN x,GEN ideal);

extern GEN isprincipalfact(GEN bnf,GEN P, GEN e, GEN C, long flag);

extern GEN logunitmatrix(GEN nf,GEN funits,GEN racunit,GEN bid);

extern GEN vconcat(GEN Q1, GEN Q2);

extern GEN make_integral(GEN nf, GEN L0, GEN f, GEN *listpr, GEN *ptd1);

extern void minim_alloc(long n, double ***q, GEN *x, double **y, double **z, double **v);

extern GEN sqred1_from_QR(GEN x, long prec);

extern GEN subgroupcondlist(GEN cyc, GEN bound, GEN listKer);

extern GEN to_Fp_simple(GEN nf, GEN x, GEN ffproj);

extern GEN to_famat_all(GEN x, GEN y);

extern GEN trivfact(void);

extern GEN arch_mul(GEN x, GEN y);

extern GEN unif_mod_f(GEN nf, GEN pr, GEN sqf);

extern GEN isprincipalfact(GEN bnf,GEN P, GEN e, GEN C, long flag);

extern GEN vconcat(GEN Q1, GEN Q2);

extern GEN idealaddtoone_i(GEN nf, GEN x, GEN y);

extern long FqX_is_squarefree(GEN P, GEN T, GEN p);

extern GEN ideallllred_elt(GEN nf, GEN I);

extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *newx);

extern void minim_alloc(long n, double ***q, GEN *x, double **y, double **z, double **v);

extern void rowselect_p(GEN A, GEN B, GEN p, long init);

/* U W V = D, Ui = U^(-1) */

GEN

compute_class_number(GEN W, GEN *D,GEN *Ui,GEN *V)

{

GEN S = smith2(W);

*Ui= ginv((GEN)S[1]);

if (V) *V = (GEN)S[2];

*D = (GEN)S[3];

if (DEBUGLEVEL>=4) msgtimer("smith/class group");

return dethnf_i(*D);

}

/* FIXME: obsolete, see zarchstar (which is much slower unfortunately). */

static GEN

get_full_rank(GEN nf, GEN v, GEN _0, GEN _1, GEN gen, long ngen, long rankmax)

Line 63 get_full_rank(GEN nf, GEN v, GEN _0, GEN _1, GEN gen,

Line 66 get_full_rank(GEN nf, GEN v, GEN _0, GEN _1, GEN gen,

limr = (limr-1)>>1;

for (k=rr; k<=limr; k++)

{

long av1=avma;

gpmem_t av1=avma;

alpha = gzero;

for (kk=k,i=1; i<=N; i++,kk/=rr)

{

Line 75 get_full_rank(GEN nf, GEN v, GEN _0, GEN _1, GEN gen,

Line 78 get_full_rank(GEN nf, GEN v, GEN _0, GEN _1, GEN gen,

vecsign[i] = (gsigne(gsubst(alpha,va,(GEN)rac[i])) > 0)? (long)_0

: (long)_1;

v1 = concatsp(v, vecsign);

if (rank(v1) == rankinit) avma=av1;

if (rank(v1) == rankinit) { avma = av1; continue; }

else

{

v=v1; rankinit++; ngen++;

gen[ngen] = (long) alpha;

if (rankinit == rankmax) return ginv(v); /* full rank */

vecsign = cgetg(rankmax+1,t_COL);

}

Line 91 GEN

Line 93 GEN

buchnarrow(GEN bnf)

{

GEN nf,_0,_1,cyc,gen,v,matsign,arch,cycgen,logs;

GEN dataunit,p1,p2,h,clh,basecl,met,u1;

GEN dataunit,p1,p2,h,basecl,met,u1;

long r1,R,i,j,ngen,sizeh,t,lo,c;

ulong av = avma;

gpmem_t av = avma;

bnf = checkbnf(bnf);

nf = checknf(bnf); r1 = nf_get_r1(nf);

Line 134 buchnarrow(GEN bnf)

Line 136 buchnarrow(GEN bnf)

vconcat(diagonal(cyc), logs),

vconcat(zeromat(ngen, r1-t), gscalmat(gdeux,r1-t))

);

clh = compute_class_number(h,&met,&u1,NULL);

lo = lg(met)-1;

lo = lg(h)-1;

for (c=1; c<=lo; c++)

met = smithrel(h,NULL,&u1);

if (gcmp1(gcoeff(met,c,c))) break;

c = lg(met);

if (DEBUGLEVEL>3) msgtimer("smith/class group");

u1 = reducemodmatrix(u1,h);

basecl = cgetg(c,t_VEC);

for (j=1; j<c; j++)

{

Line 152 buchnarrow(GEN bnf)

Line 154 buchnarrow(GEN bnf)

if (signe(p1))

{

p2 = idealmul(nf,p2, idealpow(nf,(GEN)gen[i],p1));

p1 = content(p2); if (!gcmp1(p1)) p2 = gdiv(p2,p1);

p2 = primpart(p2);

}

basecl[j] = (long)p2;

}

v = cgetg(4,t_VEC);

v[1] = lcopy(clh); setlg(met,c);

v[1] = (long)dethnf_i(h);

v[2] = (long)mattodiagonal(met);

v[2] = (long)met;

v[3] = lcopy(basecl); return gerepileupto(av, v);

v[3] = (long)basecl; return gerepilecopy(av, v);

}

/* given two coprime ideals x (integral) and id, compute alpha in x,

/********************************************************************/

* alpha = 1 mod (id), with x/alpha nearly reduced.

/** **/

/** REDUCTION MOD IDELE **/

/** **/

/********************************************************************/

static GEN

findalpha(GEN nf,GEN x,GEN id)

compute_fact(GEN nf, GEN u1, GEN gen)

{

GEN p1, y;

GEN G, basecl;

GEN alp = idealaddtoone_i(nf,x,id);

long prec,i,j, l = lg(u1), h = lg(u1[1]); /* l > 1 */

y = ideallllred_elt(nf, idealmullll(nf,x,id));

basecl = cgetg(l,t_VEC);

p1 = ground(element_div(nf,alp,y));

prec = nfgetprec(nf);

alp = gsub(alp, element_mul(nf,p1,y));

G = cgetg(3,t_VEC);

return gcmp0(alp)? y: alp;

G[2] = lgetg(1,t_MAT);

for (j=1; j<l; j++)

{

GEN g,e, z = NULL;

for (i=1; i<h; i++)

{

e = gcoeff(u1,i,j); if (!signe(e)) continue;

g = (GEN)gen[i];

if (typ(g) != t_MAT)

{

if (z)

z[2] = (long)arch_mul((GEN)z[2], to_famat_all(g, e));

else

{

z = cgetg(3,t_VEC);

z[1] = 0;

z[2] = (long)to_famat_all(g, e);

}

continue;

}

G[1] = (long)g;

g = idealpowred(nf,G,e,prec);

z = z? idealmulred(nf,z,g,prec): g;

}

z[2] = (long)famat_reduce((GEN)z[2]);

basecl[j] = (long)z;

}

return basecl;

}

/* given two coprime integral ideals x and f (f HNF), compute "small"

* non-zero a in x, such that a = 1 mod (f). GTM 193: Algo 4.3.3 */

static GEN

redideal(GEN nf,GEN x,GEN f)

{

GEN q, y, b;

if (gcmp1(gcoeff(f,1,1))) return idealred_elt(nf, x); /* f = 1 */

b = idealaddtoone_i(nf,x,f); /* a = b mod (x f) */

y = idealred_elt(nf, idealmullll(nf,x,f));

q = ground(element_div(nf,b,y));

return gsub(b, element_mul(nf,q,y)); /* != 0 since = 1 mod f */

}

static int

too_big(GEN nf, GEN bet)

{

Line 191 too_big(GEN nf, GEN bet)

Line 241 too_big(GEN nf, GEN bet)

return 0; /* not reached */

}

/* GTM 193: Algo 4.3.4. Reduce x mod idele */

static GEN

idealmodidele(GEN nf, GEN x, GEN ideal, GEN sarch, GEN arch)

_idealmodidele(GEN nf, GEN x, GEN idele, GEN sarch)

{

long av = avma,i,l;

gpmem_t av = avma;

GEN p1,p2,alp,bet,b;

GEN a,A,D,G, f = (GEN)idele[1];

nf=checknf(nf); alp=findalpha(nf,x,ideal);

G = redideal(nf, x, f);

p1=idealdiv(nf,alp,x);

D = redideal(nf, idealdiv(nf,G,x), f);

bet = element_div(nf,findalpha(nf,p1,ideal),alp);

A = element_div(nf,D,G);

if (too_big(nf,bet) > 0) { avma=av; return x; }

if (too_big(nf,A) > 0) { avma = av; return x; }

p1=(GEN)sarch[2]; l=lg(p1);

a = set_sign_mod_idele(nf, NULL, A, idele, sarch);

if (l > 1)

if (a != A && too_big(nf,A) > 0) { avma = av; return x; }

{

return idealmul(nf, a, x);

b=bet; p2=lift_intern(gmul((GEN)sarch[3],zsigne(nf,bet,arch)));

for (i=1; i<l; i++)

if (signe(p2[i])) bet = element_mul(nf,bet,(GEN)p1[i]);

if (b != bet && too_big(nf,bet) > 0) { avma=av; return x; }

}

return idealmul(nf,bet,x);

}

static GEN

GEN

idealmulmodidele(GEN nf,GEN x,GEN y, GEN ideal,GEN sarch,GEN arch)

idealmodidele(GEN bnr, GEN x)

{

return idealmodidele(nf,idealmul(nf,x,y),ideal,sarch,arch);

GEN bid = (GEN)bnr[2], fa2 = (GEN)bid[4];

GEN idele = (GEN)bid[1];

GEN sarch = (GEN)fa2[lg(fa2)-1];

return _idealmodidele(checknf(bnr), x, idele, sarch);

}

/* assume n > 0 */

/* v_pr(L0 * cx). tau = pr[5] or (more efficient) mult. table for pr[5] */

/* FIXME: should compute x^n = a I using idealred, then reduce a mod idele */

static long

fast_val(GEN nf,GEN L0,GEN cx,GEN pr,GEN tau)

{

GEN p = (GEN)pr[1];

long v = int_elt_val(nf,L0,p,tau,NULL);

if (cx)

{

long w = ggval(cx, p);

if (w) v += w * itos((GEN)pr[3]);

}

return v;

}

static GEN

idealpowmodidele(GEN nf,GEN x,GEN n, GEN ideal,GEN sarch,GEN arch)

compute_raygen(GEN nf, GEN u1, GEN gen, GEN bid)

{

long i,m,av=avma;

GEN f, fZ, basecl, module, fa, fa2, *listpr, *listep, *vecinvpi, *vectau;

GEN y;

GEN pr, t, sqf, EX, sarch, cyc;

ulong j;

long i,j,l,lp;

if (cmpis(n, 16) < 0)

/* basecl = generators in factored form */

basecl = compute_fact(nf,u1,gen);

module = (GEN)bid[1];

cyc = gmael(bid,2,2); EX = (GEN)cyc[1]; /* exponent of (O/f)^* */

f = (GEN)module[1]; fZ = gcoeff(f,1,1);

fa = (GEN)bid[3];

fa2 = (GEN)bid[4]; sarch = (GEN)fa2[lg(fa2)-1];

listpr = (GEN*)fa[1];

listep = (GEN*)fa[2];

lp = lg(listpr);

/* sqf = squarefree kernel of f */

sqf = lp <= 2? NULL: idealprodprime(nf, (GEN)listpr);

vecinvpi = (GEN*)cgetg(lp, t_VEC);

vectau = (GEN*)cgetg(lp, t_VEC);

for (i=1; i<lp; i++)

{

if (gcmp1(n)) return x;

pr = listpr[i];

x = idealpow(nf,x,n);

vecinvpi[i] = NULL; /* to be computed if needed */

x = idealmodidele(nf,x,ideal,sarch,arch);

vectau[i] = eltmul_get_table(nf, (GEN)pr[5]);

return gerepileupto(av,x);

}

i = lgefint(n)-1; m=n[i]; j=HIGHBIT;

l = lg(basecl);

while ((m&j)==0) j>>=1;

for (i=1; i<l; i++)

y = x;

for (j>>=1; j; j>>=1)

{

y = idealmul(nf,y,y);

GEN d, invpi, mulI, G, I, A, e, L, newL;

if (m&j) y = idealmul(nf,y,x);

long la, v, k;

y = idealmodidele(nf,y,ideal,sarch,arch);

/* G = [I, A=famat(L,e)] is a generator, I integral */

}

G = (GEN)basecl[i];

for (i--; i>=2; i--)

I = (GEN)G[1];

for (m=n[i],j=HIGHBIT; j; j>>=1)

A = (GEN)G[2];

L = (GEN)A[1];

e = (GEN)A[2];

/* if no reduction took place in compute_fact, everybody is still coprime

* to f + no denominators */

if (!I)

{

y = idealmul(nf,y,y);

basecl[i] = (long)famat_to_nf_modidele(nf, L, e, bid);

if (m&j) y = idealmul(nf,y,x);

continue;

y = idealmodidele(nf,y,ideal,sarch,arch);

}

return gerepileupto(av,y);

if (lg(A) == 1)

{

basecl[i] = (long)I;

continue;

}

/* compute mulI so that mulI * I coprime to f

* FIXME: use idealcoprime ??? (Less efficient. Fix idealcoprime!) */

mulI = NULL;

for (j=1; j<lp; j++)

{

invpi = vecinvpi[j];

pr = listpr[j];

v = idealval(nf, I, pr);

if (v) {

if (!invpi) invpi = vecinvpi[j] = anti_unif_mod_f(nf, pr, sqf);

t = element_pow(nf,invpi,stoi(v));

mulI = mulI? element_mul(nf, mulI, t): t;

}

/* make all components of L coprime to f.

* Assuming (L^e * I, f) = 1, then newL^e * mulI = L^e */

la = lg(e); newL = cgetg(la, t_VEC);

for (k=1; k<la; k++)

{

GEN L0, cx, LL = (GEN)L[k];

if (typ(LL) != t_COL) LL = algtobasis(nf, LL);

L0 = Q_primitive_part(LL, &cx); /* LL = L0*cx (faster element_val) */

for (j=1; j<lp; j++)

{

invpi = vecinvpi[j];

pr = listpr[j];

v = fast_val(nf, L0,cx, pr,vectau[j]); /* = val_pr(LL) */

if (v) {

if (!invpi) invpi = vecinvpi[j] = anti_unif_mod_f(nf, pr, sqf);

t = element_pow(nf,invpi,stoi(v));

LL = element_mul(nf, LL, t);

}

LL = make_integral(nf, LL, f, listpr, NULL);

newL[k] = (long)FpV_red(LL, fZ);

}

/* G in nf, = L^e mod f */

G = famat_to_nf_modideal_coprime(nf, newL, gmod(e,EX), f);

d = NULL;

if (mulI)

{

GEN mod;

/* mulI integral, d | mulI * I, sqf(dO_K) | f */

mulI = make_integral(nf,mulI,f,listpr, &d);

/* reduce mulI mod (f * (mulI,d)) */

if (!d) mod = f;

else

{

mod = hnfmodid(eltmul_get_table(nf,mulI), d); /* = (mulI,d) */

mod = idealmul(nf,f, mod);

}

mulI = colreducemodHNF(mulI, mod, NULL);

G = element_muli(nf, G, mulI);

/* L^e * I = (G/d * I) mod f */

G = colreducemodHNF(G, mod, NULL);

}

G = set_sign_mod_idele(nf,A,G,module,sarch);

I = idealmul(nf,I,G);

if (d) I = Q_div_to_int(I,d);

/* more or less useless, but cheap at this point */

I = _idealmodidele(nf,I,module,sarch);

basecl[i] = (long)I;

}

return basecl;

}

/********************************************************************/

/** **/

/** INIT RAY CLASS GROUP **/

/** **/

/********************************************************************/

static GEN

buchrayall(GEN bnf,GEN module,long flag)

{

GEN nf,cyc,gen,genplus,fa2,sarch,hmatu,u,clg,logs;

GEN nf,cyc,gen,genplus,hmatu,u,clg,logs;

GEN dataunit,p1,p2,h,clh,genray,met,u1,u2,u1old,cycgen;

GEN dataunit,p1,h,met,u1,u2,U,cycgen;

GEN racunit,bigres,bid,cycbid,genbid,x,y,funits,hmat,vecel;

long RU,Ri,i,j,ngen,lh,lo,c,av=avma;

long RU, Ri, j, ngen, lh;

const int add_gen = flag & nf_GEN;

const int do_init = flag & nf_INIT;

gpmem_t av=avma;

bnf = checkbnf(bnf); nf = checknf(bnf);

funits = check_units(bnf, "buchrayall"); RU = lg(funits);

Line 275 buchrayall(GEN bnf,GEN module,long flag)

Line 438 buchrayall(GEN bnf,GEN module,long flag)

Ri = lg(cycbid)-1; lh = ngen+Ri;

x = idealhermite(nf,module);

if (Ri || flag & (nf_INIT|nf_GEN))

if (Ri || add_gen || do_init)

{

vecel = cgetg(ngen+1,t_VEC);

for (j=1; j<=ngen; j++)

Line 285 buchrayall(GEN bnf,GEN module,long flag)

Line 448 buchrayall(GEN bnf,GEN module,long flag)

vecel[j]=(long)p1;

}

if (flag & nf_GEN)

if (add_gen)

{

genplus = cgetg(lh+1,t_VEC);

for (j=1; j<=ngen; j++)

genplus[j] = (long) idealmul(nf,(GEN)vecel[j],(GEN)gen[j]);

genplus[j] = (long)idealmul(nf,(GEN)vecel[j],(GEN)gen[j]);

for ( ; j<=lh; j++)

genplus[j] = genbid[j-ngen];

}

if (!Ri)

{

if (!(flag & nf_GEN)) clg = cgetg(3,t_VEC);

clg = cgetg(add_gen? 4: 3,t_VEC);

else

if (add_gen) clg[3] = (long)genplus;

{ clg = cgetg(4,t_VEC); clg[3] = (long)genplus; }

clg[1] = mael(bigres,1,1);

clg[2] = (long)cyc;

if (!(flag & nf_INIT)) return gerepilecopy(av,clg);

if (!do_init) return gerepilecopy(av,clg);

y = cgetg(7,t_VEC);

y[1] = lcopy(bnf);

y[1] = (long)bnf;

y[2] = lcopy(bid);

y[2] = (long)bid;

y[3] = lcopy(vecel);

y[3] = (long)vecel;

y[4] = (long)idmat(ngen);

y[5] = lcopy(clg); u = cgetg(3,t_VEC);

y[5] = (long)clg; u = cgetg(3,t_VEC);

y[6] = (long)u;

u[1] = lgetg(1,t_MAT);

u[2] = (long)idmat(RU);

return gerepileupto(av,y);

return gerepilecopy(av,y);

}

fa2 = (GEN)bid[4]; sarch = (GEN)fa2[lg(fa2)-1];

cycgen = check_and_build_cycgen(bnf);

dataunit = cgetg(RU+1,t_MAT); racunit = gmael(bigres,4,2);

Line 340 buchrayall(GEN bnf,GEN module,long flag)

Line 501 buchrayall(GEN bnf,GEN module,long flag)

vconcat(diagonal(cyc), gneg_i(logs)),

vconcat(zeromat(ngen, Ri), hmat)

);

clh = compute_class_number(h,&met,&u1,NULL);

met = smithrel(hnf(h), &U, add_gen? &u1: NULL);

u1old = u1; lo = lg(met)-1;

clg = cgetg(add_gen? 4: 3, t_VEC);

for (c=1; c<=lo; c++)

clg[1] = (long)detcyc(met);

if (gcmp1(gcoeff(met,c,c))) break;

clg[2] = (long)met;

if (add_gen) clg[3] = (long)compute_raygen(nf,u1,genplus,bid);

if (!do_init) return gerepilecopy(av, clg);

if (flag & nf_GEN)

{

GEN Id = idmat(degpol(nf[1])), arch = (GEN)module[2];

u1 = reducemodmatrix(u1,h);

genray = cgetg(c,t_VEC);

for (j=1; j<c; j++)

{

GEN *op, minus = Id, plus = Id;

long av1 = avma, s;

for (i=1; i<=lo; i++)

{

p1 = gcoeff(u1,i,j);

if (!(s = signe(p1))) continue;

if (s > 0) op = + else { op = − p1 = negi(p1); }

p1 = idealpowmodidele(nf,(GEN)genplus[i],p1,x,sarch,arch);

*op = *op==Id? p1

: idealmulmodidele(nf,*op,p1,x,sarch,arch);

}

if (minus == Id) p1 = plus;

else

{

p2 = ideleaddone_aux(nf,minus,module);

p1 = idealdivexact(nf,idealmul(nf,p2,plus),minus);

p1 = idealmodidele(nf,p1,x,sarch,arch);

}

genray[j]=lpileupto(av1,p1);

}

clg = cgetg(4,t_VEC); clg[3] = lcopy(genray);

} else clg = cgetg(3,t_VEC);

clg[1] = licopy(clh); setlg(met,c);

clg[2] = (long)mattodiagonal(met);

if (!(flag & nf_INIT)) return gerepileupto(av,clg);

u2 = cgetg(Ri+1,t_MAT);

u1 = cgetg(RU+1,t_MAT); u = (GEN)hmatu[2];

for (j=1; j<=RU; j++) { u1[j]=u[j]; setlg(u[j],RU+1); }

u += RU;

for (j=1; j<=Ri; j++) { u2[j]=u[j]; setlg(u[j],RU+1); }

p1 = lllint(u1); p2 = ginv(hmat);

y = cgetg(7,t_VEC);

y[1] = lcopy(bnf);

y[1] = (long)bnf;

y[2] = lcopy(bid);

y[2] = (long)bid;

y[3] = lcopy(vecel);

y[3] = (long)vecel;

y[4] = linv(u1old);

y[4] = (long)U;

y[5] = lcopy(clg); u = cgetg(3,t_VEC);

y[5] = (long)clg; u = cgetg(3,t_VEC);

y[6] = (long)u;

u[1] = lmul(u2,p2);

u[1] = lmul(u2, ginv(hmat));

u[2] = lmul(u1,p1);

u[2] = (long)lllint_ip(u1,100);

return gerepileupto(av,y);

return gerepilecopy(av,y);

}

GEN

Line 445 rayclassno(GEN bnf,GEN ideal)

Line 574 rayclassno(GEN bnf,GEN ideal)

{

GEN nf,h,dataunit,racunit,bigres,bid,cycbid,funits,H;

long RU,i;

ulong av = avma;

gpmem_t av = avma;

bnf = checkbnf(bnf); nf = (GEN)bnf[7];

bigres = (GEN)bnf[8]; h = gmael(bigres,1,1); /* class number */

Line 470 quick_isprincipalgen(GEN bnf, GEN x)

Line 599 quick_isprincipalgen(GEN bnf, GEN x)

GEN z = cgetg(3, t_VEC), gen = gmael3(bnf,8,1,3);

GEN idep, ep = isprincipal(bnf,x);

/* x \prod g[i]^(-ep[i]) = factorisation of principal ideal */

idep = isprincipalfact(bnf, gen, gneg(ep), x, nf_GENMAT | nf_GEN);

idep = isprincipalfact(bnf, gen, gneg(ep), x, nf_GENMAT|nf_FORCE);

z[1] = (long)ep;

z[2] = idep[2]; return z;

}

Line 478 quick_isprincipalgen(GEN bnf, GEN x)

Line 607 quick_isprincipalgen(GEN bnf, GEN x)

GEN

isprincipalrayall(GEN bnr, GEN x, long flag)

{

long av=avma,i,j,c;

long i, j, c;

GEN bnf,nf,bid,matu,vecel,ep,p1,beta,idep,y,rayclass;

gpmem_t av=avma;

GEN bnf,nf,bid,matu,vecel,ep,p1,beta,idep,ex,rayclass;

GEN divray,genray,alpha,alphaall,racunit,res,funit;

checkbnr(bnr);

Line 488 isprincipalrayall(GEN bnr, GEN x, long flag)

Line 618 isprincipalrayall(GEN bnr, GEN x, long flag)

vecel=(GEN)bnr[3];

matu =(GEN)bnr[4];

rayclass=(GEN)bnr[5];

divray = (GEN)rayclass[2]; c = lg(divray)-1;

ex = cgetg(c+1,t_COL);

if (c == 0 && !(flag & nf_GEN)) return ex;

if (typ(x) == t_VEC && lg(x) == 3)

{ idep = (GEN)x[2]; x = (GEN)x[1]; } /* precomputed */

Line 495 isprincipalrayall(GEN bnr, GEN x, long flag)

Line 628 isprincipalrayall(GEN bnr, GEN x, long flag)

idep = quick_isprincipalgen(bnf, x);

ep = (GEN)idep[1];

beta= (GEN)idep[2];

c = lg(ep);

j = lg(ep);

for (i=1; i<c; i++) /* modify beta as if gen -> vecel.gen (coprime to bid) */

for (i=1; i<j; i++) /* modify beta as if gen -> vecel.gen (coprime to bid) */

if (typ(vecel[i]) != t_INT && signe(ep[i])) /* <==> != 1 */

beta = arch_mul(to_famat_all((GEN)vecel[i], negi((GEN)ep[i])), beta);

p1 = gmul(matu, concatsp(ep, zideallog(nf,beta,bid)));

divray = (GEN)rayclass[2]; c = lg(divray);

for (i=1; i<=c; i++)

y = cgetg(c,t_COL);

ex[i] = lmodii((GEN)p1[i],(GEN)divray[i]);

for (i=1; i<c; i++)

if (!(flag & nf_GEN)) return gerepileupto(av, ex);

y[i] = lmodii((GEN)p1[i],(GEN)divray[i]);

if (!(flag & nf_GEN)) return gerepileupto(av, y);

/* compute generator */

if (lg(rayclass)<=3)

Line 512 isprincipalrayall(GEN bnr, GEN x, long flag)

Line 643 isprincipalrayall(GEN bnr, GEN x, long flag)

genray = (GEN)rayclass[3];

/* TODO: should be using nf_GENMAT and function should return a famat */

alphaall = isprincipalfact(bnf, genray, gneg(y), x, nf_GEN | nf_FORCE);

alphaall = isprincipalfact(bnf, genray, gneg(ex), x, nf_GEN | nf_FORCE);

if (!gcmp0((GEN)alphaall[1])) err(bugparier,"isprincipalray (bug1)");

res = (GEN)bnf[8];

Line 533 isprincipalrayall(GEN bnr, GEN x, long flag)

Line 664 isprincipalrayall(GEN bnr, GEN x, long flag)

alpha = gdiv(alpha,p1);

}

p1 = cgetg(4,t_VEC);

p1[1] = lcopy(y);

p1[1] = lcopy(ex);

p1[2] = (long)algtobasis(nf,alpha);

p1[3] = lcopy((GEN)alphaall[3]);

return gerepileupto(av, p1);

Line 542 isprincipalrayall(GEN bnr, GEN x, long flag)

Line 673 isprincipalrayall(GEN bnr, GEN x, long flag)

GEN

isprincipalray(GEN bnr, GEN x)

{

return isprincipalrayall(bnr,x,nf_REGULAR);

return isprincipalrayall(bnr,x,0);

}

GEN

Line 551 isprincipalraygen(GEN bnr, GEN x)

Line 682 isprincipalraygen(GEN bnr, GEN x)

return isprincipalrayall(bnr,x,nf_GEN);

}

/* N! / N^N * (4/pi)^r2 * sqrt(|D|) */

GEN

minkowski_bound(GEN D, long N, long r2, long prec)

{

long av = avma;

gpmem_t av = avma;

GEN p1;

p1 = gdiv(mpfactr(N,prec), gpowgs(stoi(N),N));

p1 = gmul(p1, gpowgs(gdivsg(4,mppi(prec)), r2));

Line 566 minkowski_bound(GEN D, long N, long r2, long prec)

Line 698 minkowski_bound(GEN D, long N, long r2, long prec)

static long

zimmertbound(long N,long R2,GEN DK)

{

long av = avma;

gpmem_t av = avma;

GEN w;

if (N < 2) return 1;

if (N < 21)

{

static double c[21][11] = {

static double c[19][11] = {

{/*2*/ 0.6931, 0.45158},

{/*3*/ 1.71733859, 1.37420604},

{/*4*/ 2.91799837, 2.50091538, 2.11943331},

Line 619 zimmertbound(long N,long R2,GEN DK)

Line 751 zimmertbound(long N,long R2,GEN DK)

avma = av; return itos(w);

}

/* all primes up to Minkowski bound factor on factorbase ? */

static void

testprime(GEN bnf, long minkowski)

{

ulong av = avma;

long pp,i,nbideal,k,pmax;

GEN f,p1,vectpp,fb,dK, nf=checknf(bnf);

byteptr delta = diffptr;

if (DEBUGLEVEL>1)

fprintferr("PHASE 1: check primes to Zimmert bound = %ld\n\n",minkowski);

f=(GEN)nf[4]; dK=(GEN)nf[3];

if (!gcmp1(f))

{

GEN different = gmael(nf,5,5);

if (DEBUGLEVEL>1)

fprintferr("**** Testing Different = %Z\n",different);

p1 = isprincipalall(bnf,different,nf_FORCE);

if (DEBUGLEVEL>1) fprintferr(" is %Z\n",p1);

}

fb=(GEN)bnf[5];

p1 = gmael(fb, lg(fb)-1, 1); /* largest p in factorbase */

pp = 0; pmax = is_bigint(p1)? VERYBIGINT: itos(p1);

if ((ulong)minkowski > maxprime()) err(primer1);

while (pp < minkowski)

{

pp += *delta++;

if (DEBUGLEVEL>1) fprintferr("*** p = %ld\n",pp);

vectpp=primedec(bnf,stoi(pp)); nbideal=lg(vectpp)-1;

/* loop through all P | p if ramified, all but one otherwise */

if (!smodis(dK,pp)) nbideal++;

for (i=1; i<nbideal; i++)

{

GEN P = (GEN)vectpp[i];

if (DEBUGLEVEL>1)

fprintferr(" Testing P = %Z\n",P);

if (cmpis(idealnorm(bnf,P), minkowski) < 1)

{

if (pp <= pmax && (k = tablesearch(fb, P, cmp_prime_ideal)))

{

if (DEBUGLEVEL>1) fprintferr(" #%ld in factor base\n",k);

}

else

{

p1 = isprincipal(bnf,P);

if (DEBUGLEVEL>1) fprintferr(" is %Z\n",p1);

}

else if (DEBUGLEVEL>1)

fprintferr(" Norm(P) > Zimmert bound\n");

}

avma = av;

}

if (DEBUGLEVEL>1) { fprintferr("End of PHASE 1.\n\n"); flusherr(); }

}

/* return \gamma_n^n if known, an upper bound otherwise */

static GEN

hermiteconstant(long n)

{

GEN h,h1;

long av;

gpmem_t av;

switch(n)

{

Line 694 hermiteconstant(long n)

Line 770 hermiteconstant(long n)

case 8: return stoi(256);

}

av = avma;

h = gpuigs(divsr(2,mppi(DEFAULTPREC)), n);

h = gpowgs(divsr(2,mppi(DEFAULTPREC)), n);

h1 = gsqr(ggamma(gdivgs(stoi(n+4),2),DEFAULTPREC));

return gerepileupto(av, gmul(h,h1));

}

Line 731 isprimitive(GEN nf)

Line 807 isprimitive(GEN nf)

}

static GEN

dft_bound()

{

if (DEBUGLEVEL>1) fprintferr("Default bound for regulator: 0.2\n");

return dbltor(0.2);

}

static GEN

regulatorbound(GEN bnf)

{

long N,R1,R2,R;

long N, R1, R2, R;

GEN nf,dKa,bound,p1,c1;

GEN nf, dKa, p1, c1;

nf = (GEN)bnf[7]; N = degpol(nf[1]);

bound = dbltor(0.2);

if (!isprimitive(nf)) return dft_bound();

if (!isprimitive(nf))

{

if (DEBUGLEVEL>1) fprintferr("Default bound for regulator: 0.2\n");

return bound;

}

dKa = absi((GEN)nf[3]);

R1 = itos(gmael(nf,2,1));

nf_get_sign(nf, &R1, &R2); R = R1+R2-1;

R2 = itos(gmael(nf,2,2)); R = R1+R2-1;

if (!R2 && N<12) c1 = gpowgs(stoi(4),N>>1); else c1 = gpowgs(stoi(N),N);

if (!R2 && N<12) c1 = gpuigs(stoi(4),N>>1); else c1 = gpuigs(stoi(N),N);

if (cmpii(dKa,c1) <= 0) return dft_bound();

if (cmpii(dKa,c1) <= 0)

{

if (DEBUGLEVEL>1) fprintferr("Default bound for regulator: 0.2\n");

return bound;

}

p1 = gsqr(glog(gdiv(dKa,c1),DEFAULTPREC));

p1 = gdivgs(gmul2n(gpuigs(gdivgs(gmulgs(p1,3),N*(N*N-1)-6*R2),R),R2),N);

p1 = divrs(gmul2n(gpowgs(divrs(mulrs(p1,3),N*(N*N-1)-6*R2),R),R2), N);

p1 = gsqrt(gdiv(p1, hermiteconstant(R)), DEFAULTPREC);

p1 = mpsqrt(gdiv(p1, hermiteconstant(R)));

if (gcmp(p1,bound) > 0) bound = p1;

if (DEBUGLEVEL>1) fprintferr("Mahler bound for regulator: %Z\n",p1);

return bound;

return gmax(p1, dbltor(0.2));

}

/* x given by its embeddings */

Line 780 norm_by_embed(long r1, GEN x)

Line 854 norm_by_embed(long r1, GEN x)

static int

is_unit(GEN M, long r1, GEN x)

{

long av = avma;

gpmem_t av = avma;

GEN Nx = ground( norm_by_embed(r1, gmul_mat_smallvec(M,x)) );

int ok = is_pm1(Nx);

avma = av; return ok;

}

#define NBMAX 5000

/* FIXME: should use smallvectors */

static GEN

minimforunits(GEN nf, long BORNE, long stockmax)

minimforunits(GEN nf, long BORNE, GEN w)

{

const long prec = MEDDEFAULTPREC;

long av = avma,n,i,j,k,s,norme,normax,*x,cmpt,r1;

long n, i, j, k, s, *x, r1, cnt = 0;

GEN u,r,S,a,M,p1;

gpmem_t av = avma;

double p;

GEN u,r,a,M;

double p, norme, normin, normax;

double **q,*v,*y,*z;

double eps=0.000001, BOUND = BORNE + eps;

double eps=0.000001, BOUND = BORNE * 1.00001;

if (DEBUGLEVEL>=2)

{

Line 804 minimforunits(GEN nf, long BORNE, long stockmax)

Line 878 minimforunits(GEN nf, long BORNE, long stockmax)

if (DEBUGLEVEL>2) fprintferr(" BOUND = %ld\n",BOUND);

flusherr();

}

r1 = nf_get_r1(nf);

r1 = nf_get_r1(nf); n = degpol(nf[1]);

a = gmael(nf,5,3); n = lg(a);

minim_alloc(n+1, &q, &x, &y, &z, &v);

minim_alloc(n, &q, &x, &y, &z, &v);

M = gprec_w(gmael(nf,5,1), prec);

n--;

a = gmul(gmael(nf,5,2), realun(prec));

u = lllgram(a,prec);

r = sqred1_from_QR(a, prec);

M = gmul(gmael(nf,5,1), u); /* embeddings of T2-reduced basis */

M = gprec_w(M, prec);

a = gmul(qf_base_change(a,u,1), realun(prec));

r = sqred1(a);

for (j=1; j<=n; j++)

{

v[j] = rtodbl(gcoeff(r,j,j));

for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j));

}

normax=0;

normax = 0.; normin = (double)BOUND;

S = cgetg(stockmax+1,t_MAT);

s=0; k=n; y[n]=z[n]=0;

s=0; k=n; cmpt=0; y[n]=z[n]=0;

x[n]=(long)(sqrt(BOUND/v[n]));

for(;;)

Line 846 minimforunits(GEN nf, long BORNE, long stockmax)

Line 915 minimforunits(GEN nf, long BORNE, long stockmax)

}

while (k>1);

if (!x[1] && y[1]<=eps) break;

if (++cmpt > NBMAX) { av=avma; return NULL; }

if (DEBUGLEVEL>8){ fprintferr("."); flusherr(); }

p = x[1]+z[1]; norme = (long)(y[1] + p*p*v[1] + eps);

if (++cnt == 5000) return NULL; /* too expensive */

p = x[1]+z[1]; norme = y[1] + p*p*v[1] + eps;

if (norme > normax) normax = norme;

if (is_unit(M,r1, x))

if (is_unit(M,r1, x)

&& (norme > 2*n /* exclude roots of unity */

|| !isnfscalar(element_pow(nf, small_to_col(x), w))))

{

if (norme < normin) normin = norme;

if (DEBUGLEVEL>=2) { fprintferr("*"); flusherr(); }

cmpt = 0;

if (++s <= stockmax)

{

p1 = cgetg(n+1,t_COL);

for (i=1; i<=n; i++) p1[i]=lstoi(x[i]);

S[s] = lmul(u,p1);

}

x[k]--;

}

if (DEBUGLEVEL>=2){ fprintferr("\n"); flusherr(); }

k = (s<stockmax)? s:stockmax; setlg(S,k+1);

avma = av; u = cgetg(4,t_VEC);

S = gerepilecopy(av, S);

u = cgetg(4,t_VEC);

u[1] = lstoi(s<<1);

u[2] = lstoi(normax);

u[2] = (long)dbltor(normax);

u[3] = (long)S; return u;

u[3] = (long)dbltor(normin);

return u;

}

#undef NBMAX

Line 901 compute_M0(GEN M_star,long N)

Line 966 compute_M0(GEN M_star,long N)

m2 = (N-n1)>>1;

for (n2=n1; n2<=m2; n2++)

{

long av = avma; n3=N-n1-n2; prec=PREC;

gpmem_t av = avma; n3=N-n1-n2; prec=PREC;

if (n1==n2 && n1==n3) /* n1 = n2 = n3 = m1 = N/3 */

{

p1=gdivgs(M_star,m1);

Line 924 compute_M0(GEN M_star,long N)

Line 989 compute_M0(GEN M_star,long N)

{ /* n3 > N/3 >= n1 */

k = n2; kk = N-2*k;

p2=gsub(M_star,gmulgs(X,k));

p3=gmul(gpuigs(stoi(kk),kk),gpuigs(gsubgs(gmul(M_star,p2),kk*kk),k));

p3=gmul(gpowgs(stoi(kk),kk),gpowgs(gsubgs(gmul(M_star,p2),kk*kk),k));

pol=gsub(p3,gmul(gmul(gpuigs(stoi(k),k),gpuigs(X,k)),gpuigs(p2,N-k)));

pol=gsub(p3,gmul(gmul(gpowgs(stoi(k),k),gpowgs(X,k)),gpowgs(p2,N-k)));

prec=gprecision(pol); if (!prec) prec = MEDDEFAULTPREC;

r=roots(pol,prec); lr = lg(r);

for (i=1; i<lr; i++)

Line 943 compute_M0(GEN M_star,long N)

Line 1008 compute_M0(GEN M_star,long N)

if (gsigne(v) <= 0) continue;

u=gmul2n(gadd(S,p6),-1);

w=gpui(P,gdivgs(stoi(-k),kk),prec);

w=gpow(P,gdivgs(stoi(-k),kk),prec);

p6=gmulsg(k,gadd(gsqr(glog(u,prec)),gsqr(glog(v,prec))));

M0_pro=gmul2n(gadd(p6,gmulsg(kk,gsqr(glog(w,prec)))),-2);

if (DEBUGLEVEL>2)

Line 961 compute_M0(GEN M_star,long N)

Line 1026 compute_M0(GEN M_star,long N)

f2 = gadd(f2,gmulsg(n2,gmul(X,Z)));

f2 = gadd(f2,gmulsg(n3,gmul(X,Y)));

f2 = gsub(f2,gmul(M,gmul(X,gmul(Y,Z))));

f3 = gsub(gmul(gpuigs(X,n1),gmul(gpuigs(Y,n2),gpuigs(Z,n3))), gun);

f3 = gsub(gmul(gpowgs(X,n1),gmul(gpowgs(Y,n2),gpowgs(Z,n3))), gun);

/* f1 = n1 X + n2 Y + n3 Z - M */

/* f2 = n1 YZ + n2 XZ + n3 XY */

/* f3 = X^n1 Y^n2 Z^n3 - 1*/

Line 1019 compute_M0(GEN M_star,long N)

Line 1084 compute_M0(GEN M_star,long N)

if (!M0) avma = av; else M0 = gerepilecopy(av, M0);

}

for (i=1;i<=4;i++) delete_var();

for (i=1;i<=4;i++) (void)delete_var();

return M0? M0: gzero;

}

static GEN

lowerboundforregulator_i(GEN bnf)

{

long N,R1,R2,RU,i,nbrootsofone,nbmin;

long N,R1,R2,RU,i;

GEN rootsofone,nf,M0,M,m,col,T2,bound,minunit,newminunit;

GEN nf,M0,M,G,bound,minunit,newminunit;

GEN vecminim,colalg,p1,pol,y;

GEN vecminim,p1,pol,y;

GEN units = check_units(bnf,"bnfcertify");

rootsofone=gmael(bnf,8,4); nbrootsofone=itos((GEN)rootsofone[1]);

nf = (GEN)bnf[7]; N = degpol(nf[1]);

nf=(GEN)bnf[7]; T2=gmael(nf,5,3); N=degpol(nf[1]);

nf_get_sign(nf, &R1, &R2); RU = R1+R2-1;

R1=itos(gmael(nf,2,1)); R2=itos(gmael(nf,2,2)); RU=R1+R2-1;

if (!RU) return gun;

if (RU==0) return gun;

units=algtobasis(bnf,units); minunit=qfeval(T2,(GEN)units[1]);

G = gmael(nf,5,2);

units = algtobasis(bnf,units);

minunit = gnorml2(gmul(G, (GEN)units[1])); /* T2(units[1]) */

for (i=2; i<=RU; i++)

{

newminunit=qfeval(T2,(GEN)units[i]);

newminunit = gnorml2(gmul(G, (GEN)units[i]));

if (gcmp(newminunit,minunit)<0) minunit=newminunit;

if (gcmp(newminunit,minunit) < 0) minunit = newminunit;

}

if (gcmpgs(minunit,1000000000)>0) return NULL;

if (gexpo(minunit) > 30) return NULL;

vecminim = minimforunits(nf,itos(gceil(minunit)),10000);

vecminim = minimforunits(nf, itos(gceil(minunit)), gmael3(bnf,8,4,1));

if (!vecminim) return NULL;

m=(GEN)vecminim[3]; nbmin=lg(m)-1;

bound = (GEN)vecminim[3];

if (nbmin==10000) return NULL;

i = gexpo(gsub(bound,minunit));

bound=gaddgs(minunit,1);

if (i > -5) err(bugparier,"lowerboundforregulator");

for (i=1; i<=nbmin; i++)

{

col=(GEN)m[i]; colalg=basistoalg(nf,col);

if (!gcmp1(lift_intern(gpuigs(colalg,nbrootsofone))))

{

newminunit=qfeval(T2,col);

if (gcmp(newminunit,bound)<0) bound=newminunit;

}

if (gcmp(bound,minunit)>0) err(talker,"bug in lowerboundforregulator");

if (DEBUGLEVEL>1)

{

fprintferr("M* = %Z\n",gprec_w(bound,3));

fprintferr("M* = %Z\n", bound);

if (DEBUGLEVEL>2)

{

p1=polx[0]; pol=gaddgs(gsub(gpuigs(p1,N),gmul(bound,p1)),N-1);

p1=polx[0]; pol=gaddgs(gsub(gpowgs(p1,N),gmul(bound,p1)),N-1);

p1 = roots(pol,DEFAULTPREC);

if (N&1) y=greal((GEN)p1[3]); else y=greal((GEN)p1[2]);

M0 = gmul2n(gmulsg(N*(N-1),gsqr(glog(y,DEFAULTPREC))),-2);

fprintferr("pol = %Z\n",pol);

fprintferr("old method: y = %Z, M0 = %Z\n",y,gprec_w(M0,3));

}

flusherr();

}

M0 = compute_M0(bound,N);

if (DEBUGLEVEL>1) { fprintferr("M0 = %Z\n",gprec_w(M0,3)); flusherr(); }

M = gmul2n(gdivgs(gdiv(gpuigs(M0,RU),hermiteconstant(RU)),N),R2);

M = gmul2n(gdivgs(gdiv(gpowgs(M0,RU),hermiteconstant(RU)),N),R2);

if (gcmp(M,dbltor(0.04))<0) return NULL;

if (gcmp(M, dbltor(0.04)) < 0) return NULL;

M = gsqrt(M,DEFAULTPREC);

if (DEBUGLEVEL>1)

{

fprintferr("(lower bound for regulator) M = %Z\n",gprec_w(M,3));

flusherr();

}

return M;

}

static GEN

lowerboundforregulator(GEN bnf)

{

long av = avma;

gpmem_t av = avma;

GEN x = lowerboundforregulator_i(bnf);

if (!x) { avma = av; x = regulatorbound(bnf); }

return x;

}

extern GEN to_Fp_simple(GEN x, GEN prh);

extern GEN Fp_PHlog(GEN a, GEN g, GEN p, GEN ord);

/* Compute a square matrix of rank length(beta) associated to a family

* (P_i), 1<=i<=length(beta), of primes s.t. N(P_i) = 1 mod pp, and

* (P_i,beta[j]) = 1 for all i,j */

Line 1105 static void

Line 1155 static void

primecertify(GEN bnf,GEN beta,long pp,GEN big)

{

long i,j,qq,nbcol,lb,nbqq,ra,N;

GEN nf,mat,mat1,qgen,decqq,newcol,Qh,Q,g,ord;

GEN nf,mat,mat1,qgen,decqq,newcol,Q,g,ord,modpr;

ord = NULL; /* gcc -Wall */

nbcol = 0; nf = (GEN)bnf[7]; N = degpol(nf[1]);

Line 1126 primecertify(GEN bnf,GEN beta,long pp,GEN big)

Line 1176 primecertify(GEN bnf,GEN beta,long pp,GEN big)

g = lift_intern(gener(qgen)); /* primitive root */

ord = decomp(stoi(qq-1));

}

Qh = prime_to_ideal(nf,Q);

modpr = zkmodprinit(nf, Q);

newcol = cgetg(lb+1,t_COL);

for (j=1; j<=lb; j++)

{

GEN t = to_Fp_simple((GEN)beta[j], Qh);

GEN t = to_Fp_simple(nf, (GEN)beta[j], modpr);

newcol[j] = (long)Fp_PHlog(t,g,qgen,ord);

}

if (DEBUGLEVEL>3)

Line 1143 primecertify(GEN bnf,GEN beta,long pp,GEN big)

Line 1193 primecertify(GEN bnf,GEN beta,long pp,GEN big)

mat1 = concatsp(mat,newcol); ra = rank(mat1);

if (ra==nbcol) continue;

if (DEBUGLEVEL>2) fprintferr(" new rank: %ld\n\n",ra);

if (DEBUGLEVEL>2) fprintferr(" new rank: %ld\n",ra);

if (++nbcol == lb) return;

mat = mat1;

}

Line 1153 primecertify(GEN bnf,GEN beta,long pp,GEN big)

Line 1203 primecertify(GEN bnf,GEN beta,long pp,GEN big)

static void

check_prime(long p, GEN bnf, GEN cyc, GEN cycgen, GEN fu, GEN mu, GEN big)

{

ulong av = avma;

gpmem_t av = avma;

long i,b, lc = lg(cyc), w = itos((GEN)mu[1]), lf = lg(fu);

GEN beta = cgetg(lf+lc, t_VEC);

Line 1178 check_prime(long p, GEN bnf, GEN cyc, GEN cycgen, GEN

Line 1228 check_prime(long p, GEN bnf, GEN cyc, GEN cycgen, GEN

long

certifybuchall(GEN bnf)

{

ulong av = avma;

gpmem_t av = avma;

long nbgen,i,j,p,N,R1,R2,R,bound;

GEN big,nf,reg,rootsofone,funits,gen,p1,gbound,cycgen,cyc;

byteptr delta = diffptr;

bnf = checkbnf(bnf); nf = (GEN)bnf[7];

N=degpol(nf[1]); if (N==1) return 1;

R1=itos(gmael(nf,2,1)); R2=itos(gmael(nf,2,2)); R=R1+R2-1;

nf_get_sign(nf, &R1, &R2); R = R1+R2-1;

funits = check_units(bnf,"bnfcertify");

testprime(bnf, zimmertbound(N,R2,absi((GEN)nf[3])));

testprimes(bnf, zimmertbound(N,R2,absi((GEN)nf[3])));

reg = gmael(bnf,8,2);

cyc = gmael3(bnf,8,1,2); nbgen = lg(cyc)-1;

gen = gmael3(bnf,8,1,3); rootsofone = gmael(bnf,8,4);

Line 1224 certifybuchall(GEN bnf)

Line 1274 certifybuchall(GEN bnf)

rootsofone = dummycopy(rootsofone);

rootsofone[2] = (long)algtobasis(nf, (GEN)rootsofone[2]);

for (p = *delta++; p <= bound; p += *delta++)

for (p = *delta++; p <= bound; ) {

check_prime(p,bnf,cyc,cycgen,funits,rootsofone,big);

NEXT_PRIME_VIADIFF(p, delta);

}

if (nbgen)

{

Line 1266 imageofgroup0(GEN gen,GEN bnr,GEN H)

Line 1318 imageofgroup0(GEN gen,GEN bnr,GEN H)

static GEN

imageofgroup(GEN gen,GEN bnr,GEN H)

{

ulong av = avma;

gpmem_t av = avma;

return gerepileupto(av,imageofgroup0(gen,bnr,H));

}

Line 1327 bnrisconductor(GEN arg0,GEN arg1,GEN arg2)

Line 1379 bnrisconductor(GEN arg0,GEN arg1,GEN arg2)

return itos(conductor(bnr,sub,-1));

}

GEN

isprincipalray_init(GEN bnf, GEN x)

{

GEN z = cgetg(3,t_VEC);

z[2] = (long)quick_isprincipalgen(bnf, x);

z[1] = (long)x; return z;

}

/* special for isprincipalrayall */

static GEN

getgen(GEN bnf, GEN gen)

{

long i,l = lg(gen);

GEN p1, g = cgetg(l, t_VEC);

GEN g = cgetg(l, t_VEC);

for (i=1; i<l; i++)

{

g[i] = (long)isprincipalray_init(bnf, (GEN)gen[i]);

p1 = cgetg(3,t_VEC); g[i] = (long)p1;

p1[1] = (long)gen[i];

p1[2] = (long)quick_isprincipalgen(bnf, (GEN)gen[i]);

}

return g;

}

Line 1352 getgen(GEN bnf, GEN gen)

Line 1408 getgen(GEN bnf, GEN gen)

GEN

conductor(GEN bnr, GEN H, long all)

{

ulong av = avma;

gpmem_t av = avma;

long r1,j,k,ep;

GEN bnf,nf,gen,bid,ideal,arch,p1,clhray,clhss,fa,arch2,bnr2,P,ex,mod;

Line 1363 conductor(GEN bnr, GEN H, long all)

Line 1419 conductor(GEN bnr, GEN H, long all)

nf = (GEN)bnf[7]; r1 = nf_get_r1(nf);

ideal= gmael(bid,1,1);

arch = gmael(bid,1,2);

if (gcmp0(H)) H = NULL;

if (H && gcmp0(H)) H = NULL;

else

if (H)

{

p1 = gauss(H, diagonal(gmael(bnr,5,2)));

if (!gcmp1(denom(p1))) err(talker,"incorrect subgroup in conductor");

Line 1384 conductor(GEN bnr, GEN H, long all)

Line 1440 conductor(GEN bnr, GEN H, long all)

ep = (all>=0)? itos((GEN)ex[k]): 1;

for (j=1; j<=ep; j++)

{

mod[1] = (long)idealdivexact(nf,ideal,pr);

mod[1] = (long)idealdivpowprime(nf,ideal,pr,gun);

clhss = orderofquotient(bnf,mod,H,gen);

if (!egalii(clhss,clhray)) break;

if (all < 0) { avma = av; return gzero; }

Line 1414 conductor(GEN bnr, GEN H, long all)

Line 1470 conductor(GEN bnr, GEN H, long all)

/* etant donne un bnr et un polynome relatif, trouve le groupe des normes

correspondant a l'extension relative en supposant qu'elle est abelienne

et que le module donnant bnr est multiple du conducteur. Verifie que

et que le module donnant bnr est multiple du conducteur. */

l'extension est bien abelienne (sous GRH) si rnf != NULL, dans ce cas

GEN

rnf est le rnf de l'extension relative. */

rnfnormgroup(GEN bnr, GEN polrel)

static GEN

rnfnormgroup0(GEN bnr, GEN polrel, GEN rnf)

{

long av=avma,i,j,reldeg,sizemat,p,pmax,nfac,k;

long i, j, reldeg, sizemat, p, nfac, k;

GEN bnf,polreldisc,discnf,nf,raycl,group,detgroup,fa,greldeg;

gpmem_t av = avma;

GEN contreld,primreld,reldisc,famo,ep,fac,col,p1,bd,upnf;

GEN bnf,index,discnf,nf,raycl,group,detgroup,fa,greldeg;

byteptr d = diffptr + 1; /* start at p = 2 */

GEN famo,ep,fac,col;

byteptr d = diffptr;

checkbnr(bnr); bnf=(GEN)bnr[1]; raycl=(GEN)bnr[5];

nf=(GEN)bnf[7];

polrel = fix_relative_pol(nf,polrel,1);

if (typ(polrel)!=t_POL) err(typeer,"rnfnormgroup");

reldeg=degpol(polrel);

reldeg = degpol(polrel);

/* reldeg-th powers are in norm group */

greldeg = stoi(reldeg);

group = diagonal(gmod((GEN)raycl[2], greldeg));

Line 1437 rnfnormgroup0(GEN bnr, GEN polrel, GEN rnf)

Line 1492 rnfnormgroup0(GEN bnr, GEN polrel, GEN rnf)

if (!signe(gcoeff(group,i,i))) coeff(group,i,i) = (long)greldeg;

detgroup = dethnf_i(group);

k = cmpis(detgroup,reldeg);

if (k<0)

if (k < 0)

{

if (rnf) return NULL;

err(talker,"not an Abelian extension in rnfnormgroup?");

}

if (!k) return gerepilecopy(av, group);

if (!rnf && !k) return gerepilecopy(av, group);

polreldisc=discsr(polrel);

if (rnf)

{

reldisc=gmael(rnf,3,1);

upnf=nfinit0(gmael(rnf,11,1),1,DEFAULTPREC);

}

else

{

reldisc = idealhermite(nf,polreldisc);

upnf = NULL;

}

reldisc = idealmul(nf, reldisc, gmael3(bnr,2,1,1));

contreld= content(reldisc);

primreld= gcmp1(contreld)? reldisc: gdiv(reldisc, contreld);

discnf = (GEN)nf[3];

k = degpol(nf[1]);

index = (GEN)nf[4];

bd = gmulsg(k, glog(absi(discnf), DEFAULTPREC));

bd = gadd(bd,glog(mpabs(det(reldisc)),DEFAULTPREC));

p1 = dbltor(reldeg * k * 2.5 + 5);

bd = gfloor(gsqr(gadd(gmulsg(4,bd),p1)));

pmax = is_bigint(bd)? 0: itos(bd);

if (rnf)

{

if (DEBUGLEVEL) fprintferr("rnfnormgroup: bound for primes = %Z\n", bd);

if (!pmax) err(warner,"rnfnormgroup: prime bound too large, can't certify");

}

sizemat=lg(group)-1;

for (p=2; !pmax || p < pmax; p += *d++)

for (p=0 ;;)

{

long oldf = -1, lfa;

/* If all pr are unramified and have the same residue degree, p =prod pr

* and including last pr^f or p^f is the same, but the last isprincipal

* is much easier! oldf is used to track this */

if (!*d) err(primer1);

NEXT_PRIME_VIADIFF_CHECK(p,d);

if (!smodis(contreld,p)) continue; /* all pr|p ramified */

if (!smodis(index, p)) continue; /* can't be treated efficiently */

fa = primedec(nf,stoi(p)); lfa = lg(fa)-1;

fa = primedec(nf, stoi(p)); lfa = lg(fa)-1;

for (i=1; i<=lfa; i++)

{

GEN pr = (GEN)fa[i];

GEN pr = (GEN)fa[i], pp, T, polr;

GEN modpr = nf_to_ff_init(nf, &pr, &T, &pp);

long f;

polr = modprX(polrel, nf, modpr);

/* if pr (probably) ramified, we have to use all (non-ram) P | pr */

if (!FqX_is_squarefree(polr, T,pp)) { oldf = 0; continue; }

famo = FqX_factor(polr, T, pp);

fac = (GEN)famo[1]; f = degpol((GEN)fac[1]);

ep = (GEN)famo[2]; nfac = lg(ep)-1;

/* check decomposition of pr has Galois type */

if (element_val(nf,polreldisc,pr) != 0)

for (j=1; j<=nfac; j++)

{

/* if pr ramified, we will have to use all (non-ram) P | pr */

if (!gcmp1((GEN)ep[j])) err(bugparier,"rnfnormgroup");

if (idealval(nf,primreld,pr)!=0) { oldf = 0; continue; }

if (degpol(fac[j]) != f)

err(talker,"non Galois extension in rnfnormgroup");

famo=idealfactor(upnf,rnfidealup(rnf,pr));

ep=(GEN)famo[2];

fac=(GEN)famo[1];

nfac=lg(ep)-1;

f = itos(gmael(fac,1,4));

for (j=1; j<=nfac; j++)

{

if (!gcmp1((GEN)ep[j])) err(bugparier,"rnfnormgroup");

if (itos(gmael(fac,j,4)) != f)

{

if (rnf) return NULL;

err(talker,"non Galois extension in rnfnormgroup");

}

else

{

famo=nffactormod(nf,polrel,pr);

ep=(GEN)famo[2];

fac=(GEN)famo[1];

nfac=lg(ep)-1;

f = degpol((GEN)fac[1]);

for (j=1; j<=nfac; j++)

{

if (!gcmp1((GEN)ep[j])) err(bugparier,"rnfnormgroup");

if (degpol(fac[j]) != f)

{

if (rnf) return NULL;

err(talker,"non Galois extension in rnfnormgroup");

}

if (oldf < 0) oldf = f; else if (oldf != f) oldf = 0;

if (f == reldeg) continue; /* reldeg-th powers already included */

if (oldf && i == lfa && !smodis(discnf, p)) pr = stoi(p);

/* pr^f = N P, P | pr, hence is in norm group */

col = gmulsg(f, isprincipalrayall(bnr,pr,nf_REGULAR));

col = gmulsg(f, isprincipalrayall(bnr,pr,0));

group = hnf(concatsp(group, col));

detgroup = dethnf_i(group);

k = cmpis(detgroup,reldeg);

if (k < 0)

{

if (rnf) return NULL;

err(talker,"not an Abelian extension in rnfnormgroup");

}

if (!k) { cgiv(detgroup); return gerepileupto(av,group); }

if (!rnf && !k) { cgiv(detgroup); return gerepileupto(av,group); }

}

if (k>0) err(bugparier,"rnfnormgroup");

cgiv(detgroup); return gerepileupto(av,group);

}

GEN

static int

rnfnormgroup(GEN bnr, GEN polrel)

rnf_is_abelian(GEN nf, GEN pol)

{

return rnfnormgroup0(bnr,polrel,NULL);

GEN ro, rores, nfL, L, mod, d;

long i,j, l;

nf = checknf(nf);

L = rnfequation(nf,pol);

mod = dummycopy(L); setvarn(mod, varn(nf[1]));

nfL = _initalg(mod, nf_PARTIALFACT, DEFAULTPREC);

ro = nfroots(nfL, L);

l = lg(ro)-1;

if (l != degpol(pol)) return 0;

if (isprime(stoi(l))) return 1;

ro = Q_remove_denom(ro, &d);

if (!d) rores = ro;

else

{

rores = cgetg(l+1, t_VEC);

for (i=1; i<=l; i++)

rores[i] = (long)rescale_pol((GEN)ro[i], d);

}

/* assume roots are sorted by increasing degree */

for (i=1; i<l; i++)

for (j=1; j<i; j++)

{

GEN a = RX_RXQ_compo((GEN)rores[j], (GEN)ro[i], mod);

GEN b = RX_RXQ_compo((GEN)rores[i], (GEN)ro[j], mod);

if (d) a = gmul(a, gpowgs(d, degpol(ro[i]) - degpol(ro[j])));

if (!gegal(a, b)) return 0;

}

return 1;

}

/* Etant donne un bnf et un polynome relatif polrel definissant une extension

Line 1562 rnfnormgroup(GEN bnr, GEN polrel)

Line 1588 rnfnormgroup(GEN bnr, GEN polrel)

a l'extension definie par polrel sous la forme

[[ideal,arch],[hm,cyc,gen],group] ou [ideal,arch] est le conducteur, et

[hm,cyc,gen] est le groupe de classes de rayon correspondant. Verifie

(sous GRH) que polrel definit bien une extension abelienne si flag != 0 */

que polrel definit bien une extension abelienne si flag != 0 */

GEN

rnfconductor(GEN bnf, GEN polrel, long flag)

{

long av=avma,tetpil,R1,i,v;

long R1, i;

GEN nf,module,rnf,arch,bnr,group,p1,pol2;

gpmem_t av = avma;

GEN nf,module,arch,bnr,group,p1,pol2;

bnf = checkbnf(bnf); nf=(GEN)bnf[7];

bnf = checkbnf(bnf); nf = (GEN)bnf[7];

if (typ(polrel)!=t_POL) err(typeer,"rnfconductor");

module=cgetg(3,t_VEC); R1=itos(gmael(nf,2,1));

p1=unifpol(nf, polrel, 0);

v=varn(polrel);

p1=unifpol((GEN)bnf[7],polrel,0);

p1=denom(gtovec(p1));

pol2=gsubst(polrel,v,gdiv(polx[v],p1));

pol2=rescale_pol(polrel, p1);

pol2=gmul(pol2,gpuigs(p1,degpol(pol2)));

if (flag && !rnf_is_abelian(bnf, pol2)) { avma = av; return gzero; }

if (flag)

{

module = cgetg(3,t_VEC);

rnf=rnfinitalg(bnf,pol2,DEFAULTPREC);

module[1] = rnfdiscf(nf,pol2)[1];

module[1]=mael(rnf,3,1);

R1 = nf_get_r1(nf); arch = cgetg(R1+1,t_VEC);

}

module[2] = (long)arch; for (i=1; i<=R1; i++) arch[i]=un;

else

bnr = buchrayall(bnf,module,nf_INIT | nf_GEN);

{

group = rnfnormgroup(bnr,pol2);

rnf=NULL;

module[1]=rnfdiscf(nf,pol2)[1];

}

arch=cgetg(R1+1,t_VEC);

module[2]=(long)arch; for (i=1; i<=R1; i++) arch[i]=un;

bnr=buchrayall(bnf,module,nf_INIT | nf_GEN);

group=rnfnormgroup0(bnr,pol2,rnf);

if (!group) { avma = av; return gzero; }

tetpil=avma;

return gerepileupto(av, conductor(bnr,group,1));

return gerepile(av,tetpil,conductor(bnr,group,1));

}

/* Given a number field bnf=bnr[1], a ray class group structure

Line 1604 rnfconductor(GEN bnf, GEN polrel, long flag)

Line 1621 rnfconductor(GEN bnf, GEN polrel, long flag)

static GEN

discrayrelall(GEN bnr, GEN H, long flag)

{

ulong av = avma;

gpmem_t av = avma;

long r1,j,k,ep, nz, flrel = flag & nf_REL, flcond = flag & nf_COND;

GEN bnf,nf,gen,bid,ideal,arch,p1,clhray,clhss,fa,arch2,idealrel,P,ex,mod,dlk;

Line 1639 discrayrelall(GEN bnr, GEN H, long flag)

Line 1656 discrayrelall(GEN bnr, GEN H, long flag)

mod[1] = (long)ideal;

for (j=1; j<=ep; j++)

{

mod[1] = (long)idealdivexact(nf,(GEN)mod[1],pr);

mod[1] = (long)idealdivpowprime(nf,(GEN)mod[1],pr,gun);

clhss = orderofquotient(bnf,mod,H,gen);

if (flcond && j==1 && egalii(clhss,clhray)) { avma = av; return gzero; }

if (is_pm1(clhss)) { S = addis(S, ep-j+1); break; }

S = addii(S, clhss);

}

idealrel = flrel? idealmul(nf,idealrel, idealpow(nf,pr, S))

idealrel = flrel? idealmulpowprime(nf,idealrel, pr, S)

: mulii(idealrel, powgi(idealnorm(nf,pr),S));

}

dlk = flrel? idealdivexact(nf,idealpow(nf,ideal,clhray), idealrel)

Line 1673 discrayrelall(GEN bnr, GEN H, long flag)

Line 1690 discrayrelall(GEN bnr, GEN H, long flag)

static GEN

discrayabsall(GEN bnr, GEN subgroup,long flag)

{

ulong av = avma;

gpmem_t av = avma;

long degk,clhray,n,R1;

GEN z,p1,D,dk,nf,dkabs,bnf;

Line 1717 discrayrelcond(GEN bnr, GEN H)

Line 1734 discrayrelcond(GEN bnr, GEN H)

GEN

discrayabs(GEN bnr, GEN H)

{

return discrayabsall(bnr,H,nf_REGULAR);

return discrayabsall(bnr,H,0);

}

GEN

Line 1732 discrayabscond(GEN bnr, GEN H)

Line 1749 discrayabscond(GEN bnr, GEN H)

GEN

bnrconductorofchar(GEN bnr, GEN chi)

{

ulong av = avma;

gpmem_t av = avma;

long nbgen,i;

GEN p1,m,U,d1,cyc;

GEN m,U,d1,cyc;

checkbnrgen(bnr);

cyc = gmael(bnr,5,2); nbgen = lg(cyc)-1;

Line 1746 bnrconductorofchar(GEN bnr, GEN chi)

Line 1763 bnrconductorofchar(GEN bnr, GEN chi)

for (i=1; i<=nbgen; i++)

{

if (typ(chi[i]) != t_INT) err(typeer,"conductorofchar");

p1 = cgetg(2,t_COL); m[i] = (long)p1;

m[i] = (long)_col(mulii((GEN)chi[i], divii(d1, (GEN)cyc[i])));

p1[1] = lmulii((GEN)chi[i], divii(d1, (GEN)cyc[i]));

}

p1 = cgetg(2,t_COL); m[i] = (long)p1;

m[i] = (long)_col(d1);

p1[1] = (long)d1; U = (GEN)hnfall(m)[2];

U = (GEN)hnfall(m)[2];

setlg(U,nbgen+1);

for (i=1; i<=nbgen; i++) setlg(U[i],nbgen+1); /* U = Ker chi */

return gerepileupto(av, conductor(bnr,U,0));

Line 1761 bnrconductorofchar(GEN bnr, GEN chi)

Line 1777 bnrconductorofchar(GEN bnr, GEN chi)

GEN

rayclassnolist(GEN bnf,GEN lists)

{

ulong av = avma;

gpmem_t av = avma;

long i,j,lx,ly;

GEN blist,ulist,Llist,h,b,u,L,m;

Line 1875 factorpow(GEN fa,long n)

Line 1891 factorpow(GEN fa,long n)

GEN

discrayabslist(GEN bnf,GEN lists)

{

ulong av = avma;

gpmem_t av = avma;

long ii,jj,i,j,k,clhss,ep,clhray,lx,ly,r1,degk,nz;

long n,R1,lP;

GEN hlist,blist,dlist,nf,dkabs,b,h,d;

Line 2016 zsimp(GEN bid, GEN matunit)

Line 2032 zsimp(GEN bid, GEN matunit)

static GEN

zsimpjoin(GEN bidsimp, GEN bidp, GEN dummyfa, GEN matunit)

{

long i,l1,l2,nbgen,c, av = avma;

long i, l1, l2, nbgen, c;

gpmem_t av = avma;

GEN U,U1,U2,cyc1,cyc2,cyc,u1u2,met, y = cgetg(5,t_VEC);

y[1] = (long)vconcat((GEN)bidsimp[1],dummyfa);

Line 2068 rayclassnointern(GEN blist, GEN h)

Line 2085 rayclassnointern(GEN blist, GEN h)

return L;

}

void rowselect_p(GEN A, GEN B, GEN p, long init);

static GEN

rayclassnointernarch(GEN blist, GEN h, GEN matU)

{

Line 2113 rayclassnointernarch(GEN blist, GEN h, GEN matU)

Line 2128 rayclassnointernarch(GEN blist, GEN h, GEN matU)

GEN

decodemodule(GEN nf, GEN fa)

{

long n,k,j,fauxpr,av=avma;

long n, k, j, fauxpr;

gpmem_t av=avma;

GEN g,e,id,pr;

nf = checknf(nf);

Line 2127 decodemodule(GEN nf, GEN fa)

Line 2143 decodemodule(GEN nf, GEN fa)

fauxpr = itos((GEN)g[k]);

j = (fauxpr%n)+1; fauxpr /= n*n;

pr = (GEN)primedec(nf,stoi(fauxpr))[j];

id = idealmul(nf,id, idealpow(nf,pr,(GEN)e[k]));

id = idealmulpowprime(nf,id, pr,(GEN)e[k]);

}

return gerepileupto(av,id);

}

Line 2157 discrayabslistarchintern(GEN bnf, GEN arch, long bound

Line 2173 discrayabslistarchintern(GEN bnf, GEN arch, long bound

long degk,i,j,k,p2s,lfa,lp1,sqbou,cex, allarch;

long ffs,fs,resp,flbou,nba, k2,karch,kka,nbarch,jjj,jj,square;

long ii2,ii,ly,clhray,lP,ep,S,clhss,normps,normi,nz,r1,R1,n,c;

ulong q, av0 = avma, av,av1,lim;

ulong q;

gpmem_t av0 = avma, av, av1, lim;

GEN nf,p,z,p1,p2,p3,fa,pr,normp,ideal,bidp,z2,matarchunit;

GEN funits,racunit,embunit,sous,clh,sousray,raylist;

GEN clhrayall,discall,faall,Id,idealrel,idealrelinit;

Line 2170 discrayabslistarchintern(GEN bnf, GEN arch, long bound

Line 2187 discrayabslistarchintern(GEN bnf, GEN arch, long bound

mod = Id = dlk = ideal = clhrayall = discall = faall = NULL;

/* ce qui suit recopie d'assez pres ideallistzstarall */

if (DEBUGLEVEL>2) timer2();

if (DEBUGLEVEL>2) (void)timer2();

bnf = checkbnf(bnf); flbou=0;

nf = (GEN)bnf[7]; r1 = nf_get_r1(nf);

degk = degpol(nf[1]);

Line 2219 discrayabslistarchintern(GEN bnf, GEN arch, long bound

Line 2236 discrayabslistarchintern(GEN bnf, GEN arch, long bound

if (bound > (long)maxprime()) err(primer1);

for (p[2]=0; p[2]<=bound; )

{

p[2] += *ptdif++;

NEXT_PRIME_VIADIFF(p[2], ptdif);

if (!flbou && p[2]>sqbou)

{

if (DEBUGLEVEL>1) fprintferr("\nStarting rayclassno computations\n");

Line 2387 discrayabslistarchintern(GEN bnf, GEN arch, long bound

Line 2404 discrayabslistarchintern(GEN bnf, GEN arch, long bound

if (!allarch && nba)

{

p1 = (GEN)primedec(nf,gprime)[ffs%degk+1];

ideal = idealmul(nf,ideal,idealpow(nf,p1,(GEN)ex[k]));

ideal = idealmulpowprime(nf,ideal,p1,(GEN)ex[k]);

}

S=0; clhss=0;

normi = ii; normps= itos(gpuigs(gprime,resp));

normi = ii; normps= itos(gpowgs(gprime,resp));

for (j=1; j<=ep; j++)

{

GEN fad, fad1, fad2;

Line 2485 GEN

Line 2502 GEN

discrayabslistarchsquare(GEN bnf, GEN arch, long bound)

{ return discrayabslistarchintern(bnf,arch,bound, -1); }

static GEN

/* Let bnr1, bnr2 be such that mod(bnr2) | mod(bnr1), compute the

computehuv(GEN bnr,GEN id, GEN arch)

matrix of the surjective map Cl(bnr1) ->> Cl(bnr2) */

GEN

bnrGetSurj(GEN bnr1, GEN bnr2)

{

long i,nbgen, av = avma;

long nbg, i;

GEN bnf,bnrnew,listgen,P,U,DC;

GEN gen, M;

GEN newmod=cgetg(3,t_VEC);

newmod[1]=(long)id;

newmod[2]=(long)arch;

bnf=(GEN)bnr[1];

gen = gmael(bnr1, 5, 3);

bnrnew=buchrayall(bnf,newmod,nf_INIT);

nbg = lg(gen) - 1;

listgen=gmael(bnr,5,3); nbgen=lg(listgen);

P=cgetg(nbgen,t_MAT);

M = cgetg(nbg + 1, t_MAT);

for (i=1; i<nbgen; i++)

for (i = 1; i <= nbg; i++)

P[i] = (long)isprincipalray(bnrnew,(GEN)listgen[i]);

M[i] = (long)isprincipalray(bnr2, (GEN)gen[i]);

DC=diagonal(gmael(bnrnew,5,2));

U=(GEN)hnfall(concatsp(P,DC))[2];

return M;

setlg(U,nbgen); for (i=1; i<nbgen; i++) setlg(U[i], nbgen);

return gerepileupto(av, hnf(concatsp(U,diagonal(gmael(bnr,5,2)))));

}

/* 0 if hinv*list[i] has a denominator for all i, 1 otherwise */

/* Kernel of the above map */

static int

static GEN

hnflistdivise(GEN list,GEN h)

bnrGetKer(GEN bnr, GEN mod2)

{

long av = avma, i, I = lg(list);

long i, n;

GEN hinv = ginv(h);

gpmem_t av = avma;

GEN P,U, bnr2 = buchrayall(bnr,mod2,nf_INIT);

for (i=1; i<I; i++)

P = bnrGetSurj(bnr, bnr2); n = lg(P);

if (gcmp1(denom(gmul(hinv,(GEN)list[i])))) break;

P = concatsp(P, diagonal(gmael(bnr2,5,2)));

avma = av; return i < I;

U = (GEN)hnfall(P)[2]; setlg(U,n);

for (i=1; i<n; i++) setlg(U[i], n);

U = concatsp(U, diagonal(gmael(bnr, 5,2)));

return gerepileupto(av, hnf(U));

}

static GEN

subgroupcond(GEN bnr, long indexbound)

subgroupcond(GEN bnr, GEN indexbound)

{

ulong av = avma;

gpmem_t av = avma;

long i,j,lgi,lp;

GEN li,p1,lidet,perm,nf,bid,ideal,arch,primelist,listkernels;

GEN li,p1,lidet,perm,nf,bid,ideal,arch,arch2,primelist,listKer;

GEN mod = cgetg(3, t_VEC);

checkbnrgen(bnr); bid=(GEN)bnr[2];

ideal=gmael(bid,1,1);

arch =gmael(bid,1,2); nf=gmael(bnr,1,7);

primelist=gmael(bid,3,1); lp=lg(primelist)-1;

listkernels=cgetg(lp+lg(arch),t_VEC);

mod[2] = (long)arch;

listKer=cgetg(lp+lg(arch),t_VEC);

for (i=1; i<=lp; )

{

p1=idealdiv(nf,ideal,(GEN)primelist[i]);

mod[1] = (long)idealdivpowprime(nf,ideal,(GEN)primelist[i],gun);

listkernels[i++]=(long)computehuv(bnr,p1,arch);

listKer[i++] = (long)bnrGetKer(bnr,mod);

}

mod[1] = (long)ideal; arch2 = dummycopy(arch);

mod[2] = (long)arch2;

for (j=1; j<lg(arch); j++)

if (signe((GEN)arch[j]))

{

p1=dummycopy(arch); p1[j]=zero;

arch2[j] = zero;

listkernels[i++]=(long)computehuv(bnr,ideal,p1);

listKer[i++] = (long)bnrGetKer(bnr,mod);

arch2[j] = un;

}

setlg(listkernels,i);

setlg(listKer,i);

li=subgrouplist(gmael(bnr,5,2),indexbound);

li = subgroupcondlist(gmael(bnr,5,2),indexbound,listKer);

lgi=lg(li);

lgi = lg(li);

for (i=1,j=1; j<lgi; j++)

if (!hnflistdivise(listkernels, (GEN)li[j])) li[i++] = li[j];

/* sort by increasing index */

lgi = i; setlg(li,i); lidet=cgetg(lgi,t_VEC);

lidet = cgetg(lgi,t_VEC);

for (i=1; i<lgi; i++) lidet[i]=(long)dethnf_i((GEN)li[i]);

for (i=1; i<lgi; i++) lidet[i] = (long)dethnf_i((GEN)li[i]);

perm = sindexsort(lidet); p1=li; li=cgetg(lgi,t_VEC);

perm = sindexsort(lidet); p1 = li; li = cgetg(lgi,t_VEC);

for (i=1; i<lgi; i++) li[i] = p1[perm[lgi-i]];

return gerepilecopy(av,li);

}

GEN

subgrouplist0(GEN bnr, long indexbound, long all)

subgrouplist0(GEN bnr, GEN indexbound, long all)

{

if (typ(bnr)!=t_VEC) err(typeer,"subgrouplist");

if (lg(bnr)!=1 && typ(bnr[1])!=t_INT)

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>