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

Return to buch2.c CVS log

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

Diff for /OpenXM_contrib/pari-2.2/src/basemath/Attic/buch2.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 gscal(GEN x,GEN y);

extern GEN nfbasic_to_nf(nfbasic_t *T, GEN ro, long prec);

extern GEN get_nfindex(GEN bas);

extern GEN sqred1_from_QR(GEN x, long prec);

extern GEN computeGtwist(GEN nf, GEN vdir);

extern GEN famat_mul(GEN f, GEN g);

extern GEN famat_to_nf(GEN nf, GEN f);

extern GEN famat_to_arch(GEN nf, GEN fa, long prec);

extern GEN to_famat_all(GEN x, GEN y);

extern int addcolumntomatrix(GEN V, GEN invp, GEN L);

extern double check_bach(double cbach, double B);

extern GEN gmul_mat_smallvec(GEN x, GEN y);

extern GEN gmul_mati_smallvec(GEN x, GEN y);

extern GEN get_arch(GEN nf,GEN x,long prec);

extern GEN get_arch_real(GEN nf,GEN x,GEN *emb,long prec);

extern GEN get_roots(GEN x,long r1,long ru,long prec);

extern GEN get_roots(GEN x,long r1,long prec);

extern void get_nf_matrices(GEN nf, long small);

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

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

extern GEN init_idele(GEN nf);

extern GEN norm_by_embed(long r1, GEN x);

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

extern GEN idealmulpowprime(GEN nf, GEN x, GEN vp, GEN n);

extern GEN arch_mul(GEN x, GEN y);

extern GEN vecdiv(GEN x, GEN y);

extern void wr_rel(GEN col);

extern GEN vecmul(GEN x, GEN y);

extern void dbg_rel(long s, GEN col);

extern GEN mul_real(GEN x, GEN y);

#define SFB_MAX 2

#define SFB_STEP 2

Line 48 static const int CBUCHG = (1<<RANDOM_BITS) - 1;

Line 53 static const int CBUCHG = (1<<RANDOM_BITS) - 1;

static const int randshift = BITS_IN_RANDOM-1 - RANDOM_BITS;

#undef RANDOM_BITS

static long KC,KCZ,KCZ2,MAXRELSUP;

/* used by factor[elt|gen|gensimple] to return factorizations of smooth elts

static long primfact[500],expoprimfact[500];

* HACK: MAX_FACTOR_LEN never checked, we assume default value is enough

static GEN *idealbase, vectbase, FB, numFB, powsubFB, numideal;

* (since elts have small norm) */

static long primfact[500], exprimfact[500];

/* FB[i] i-th rational prime used in factor base

/* a factor base contains only non-inert primes

* numFB[i] index k such that FB[k]=i (0 if i is not prime)

* KC = # of P in factor base (p <= n, NP <= n2)

* KC2= # of P assumed to generate class group (NP <= n2)

* vectbase vector of all ideals in FB

* KCZ = # of rational primes under ideals counted by KC

* vecbase o subFB part of FB used to build random relations

* KCZ2= same for KC2 */

* powsubFB array lg(subFB) x (CBUCHG+1) all powers up to CBUCHG

typedef struct {

GEN FB; /* FB[i] = i-th rational prime used in factor base */

* idealbase[i] prime ideals above i in FB

GEN LP; /* vector of all prime ideals in FB */

* numideal[i] index k such that idealbase[k] = i.

GEN *LV; /* LV[p] = vector of P|p, NP <= n2

* isclone() is set for LV[p] iff all P|p are in FB

* mat all relations found (as long integers, not reduced)

* LV[i], i not prime or i > n2, is undefined! */

* cglob lg(mat) = total number of relations found so far

GEN iLP; /* iLP[p] = i such that LV[p] = [LP[i],...] */

long KC, KCZ, KCZ2;

* Use only non-inert primes, coprime to discriminant index F:

GEN subFB; /* LP o subFB = part of FB used to build random relations */

* KC = number of prime ideals in factor base (norm < Bach cst)

GEN powsubFB; /* array of [P^0,...,P^CBUCHG], P in LP o subFB */

* KC2= number of prime ideals assumed to generate class group (>= KC)

GEN perm; /* permutation of LP used to represent relations [updated by

hnfspec/hnfadd: dense rows come first] */

* KCZ = number of rational primes under ideals counted by KC

} FB_t;

* KCZ2= same for KC2

/* x a t_VECSMALL */

static long

ccontent(GEN x)

{

long i, l = lg(x), s=labs(x[1]);

long i, l = lg(x), s = labs(x[1]);

for (i=2; i<l && s>1; i++) s = cgcd(x[i],s);

for (i=2; i<l && s!=1; i++) s = cgcd(x[i],s);

return s;

}

static void

desallocate(GEN **M)

desallocate(GEN **M, GEN *first_nz)

{

GEN *m = *M;

long i;

if (m)

{

for (i=lg(m)-1; i; i--) free(m[i]);

free(m); *M = NULL;

free((void*)*M); *M = NULL;

free((void*)*first_nz); *first_nz = NULL;

}

/* Return the list of indexes or the primes chosen for the subFB.

GEN

* Fill vperm (if !=0): primes ideals sorted by increasing norm (except the

cgetalloc(GEN x, size_t l, long t)

* ones in subFB come first [dense rows come first for hnfspec])

* ss = number of rational primes whose divisors are all in FB

static GEN

subFBgen(long N,long m,long minsFB,GEN vperm, long *ptss)

{

long av = avma,i,j, lv=lg(vectbase),s=0,s1=0,n=0,ss=0,z=0;

x = (GEN)gprealloc((void*)x, l * sizeof(long));

GEN y1,y2,subFB,perm,perm1,P,Q;

x[0] = evaltyp(t) | evallg(l); return x;

double prod;

}

(void)new_chunk(lv); /* room for subFB */

static void

y1 = cgetg(lv,t_COL);

reallocate(long max, GEN *M, GEN *first_nz)

y2 = cgetg(lv,t_COL);

{

for (i=1,P=(GEN)vectbase[i];;P=Q)

size_t l = max+1;

{ /* we'll sort ideals by norm (excluded ideals = "zero") */

*M = cgetalloc(*M, l, t_VEC);

long e = itos((GEN)P[3]);

*first_nz = cgetalloc(*first_nz, l, t_VECSMALL);

long ef= e*itos((GEN)P[4]);

}

s1 += ef;

/* don't take P|p if P ramified, or all other Q|p already there */

y2[i] = (long)powgi((GEN)P[1],(GEN)P[4]);

static int

/* take only unramified ideals */

ok_subFB(FB_t *F, long t, GEN D)

if (e>1) { y1[i]=zero; s=0; z++; } else { y1[i]=y2[i]; s += ef; }

{

GEN LP, P = (GEN)F->LP[t];

long p = itos((GEN)P[1]);

LP = F->LV[p];

return smodis(D,p) && (!isclone(LP) || t != F->iLP[p] + lg(LP)-1);

}

i++; Q = (GEN)vectbase[i];

/* set subFB, reset F->powsubFB

if (i == lv || !egalii((GEN)P[1], (GEN)Q[1]))

* Fill F->perm (if != NULL): primes ideals sorted by increasing norm (except

{ /* don't take all P above a given p (delete the last one) */

* the ones in subFB come first [dense rows for hnfspec]) */

if (s == N) { y1[i-1]=zero; z++; }

static int

if (s1== N) ss++;

subFBgen(FB_t *F, GEN nf, double PROD, long minsFB)

if (i == lv) break;

{

s = s1 = 0;

GEN y, perm, yes, no, D = (GEN)nf[3];

}

long i, j, k, iyes, ino, lv = F->KC + 1;

double prod;

const int init = (F->perm == NULL);

gpmem_t av;

if (init)

{

F->LP = cgetg(lv, t_VEC);

F->perm = cgetg(lv, t_VECSMALL);

}

if (z+minsFB >= lv) return NULL;

prod = 1.0;

av = avma;

perm = sindexsort(y1) + z; /* skip "zeroes" (excluded ideals) */

y = cgetg(lv,t_COL); /* Norm P */

for(;;)

for (k=0, i=1; i <= F->KCZ; i++)

{

if (++n > minsFB && (z+n >= lv || prod > m + 0.5)) break;

GEN LP = F->LV[F->FB[i]];

prod *= gtodouble((GEN)y1[perm[n]]);

long l = lg(LP);

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

{

GEN P = (GEN)LP[j];

k++;

y[k] = (long)powgi((GEN)P[1], (GEN)P[4]);

F->LP[k] = (long)P; /* noop if init = 1 */

}

if (prod < m) return NULL;

/* perm sorts LP by increasing norm */

n--;

perm = sindexsort(y);

no = cgetg(lv, t_VECSMALL); ino = 1;

yes = cgetg(lv, t_VECSMALL); iyes = 1;

prod = 1.0;

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

{

long t = perm[i];

if (!ok_subFB(F, t, D)) { no[ino++] = t; continue; }

/* take the first (non excluded) n ideals (wrt norm), put them first, and

yes[iyes++] = t;

* sort the rest by increasing norm */

prod *= (double)itos((GEN)y[t]);

for (j=1; j<=n; j++) y2[perm[j]] = zero;

if (iyes > minsFB && prod > PROD) break;

perm1 = sindexsort(y2); avma = av;

}

if (i == lv) return 0;

subFB = cgetg(n+1, t_VECSMALL);

avma = av; /* HACK: gcopy(yes) still safe */

if (vperm)

if (init)

{

for (j=1; j<=n; j++) vperm[j] = perm[j];

GEN p = F->perm;

for ( ; j<lv; j++) vperm[j] = perm1[j];

for (j=1; j<iyes; j++) p[j] = yes[j];

for (i=1; i<ino; i++, j++) p[j] = no[i];

for ( ; j<lv; j++) p[j] = perm[j];

}

for (j=1; j<=n; j++) subFB[j] = perm[j];

setlg(yes, iyes);

F->subFB = gcopy(yes);

F->powsubFB = NULL; return 1;

}

static int

subFBgen_increase(FB_t *F, GEN nf, long step)

{

GEN yes, D = (GEN)nf[3];

long i, iyes, lv = F->KC + 1, minsFB = lg(F->subFB)-1 + step;

if (DEBUGLEVEL)

yes = cgetg(minsFB+1, t_VECSMALL); iyes = 1;

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

{

if (DEBUGLEVEL>3)

long t = F->perm[i];

{

if (!ok_subFB(F, t, D)) continue;

fprintferr("\n***** IDEALS IN FACTORBASE *****\n\n");

for (i=1; i<=KC; i++) fprintferr("no %ld = %Z\n",i,vectbase[i]);

yes[iyes++] = t;

fprintferr("\n***** IDEALS IN SUB FACTORBASE *****\n\n");

if (iyes > minsFB) break;

outerr(vecextract_p(vectbase,subFB));

fprintferr("\n***** INITIAL PERMUTATION *****\n\n");

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

}

msgtimer("sub factorbase (%ld elements)",n);

}

powsubFB = NULL;

if (i == lv) return 0;

*ptss = ss; return subFB;

F->subFB = yes;

F->powsubFB = NULL; return 1;

}

static GEN

mulred(GEN nf,GEN x, GEN I, long prec)

{

ulong av = avma;

gpmem_t av = avma;

GEN y = cgetg(3,t_VEC);

y[1] = (long)idealmulh(nf,I,(GEN)x[1]);

Line 185 mulred(GEN nf,GEN x, GEN I, long prec)

Line 219 mulred(GEN nf,GEN x, GEN I, long prec)

/* Compute powers of prime ideals (P^0,...,P^a) in subFB (a > 1)

* powsubFB[j][i] contains P_i^j in LLL form + archimedean part in

* MULTIPLICATIVE form (logs are expensive)

* MULTIPLICATIVE form (logs are expensive) */

static void

powsubFBgen(GEN nf,GEN subFB,long a,long prec)

powsubFBgen(FB_t *F, GEN nf, long a, long prec)

{

long i,j, n = lg(subFB), RU = lg(nf[6]);

long i,j, n = lg(F->subFB), RU = lg(nf[6]);

GEN *pow, arch0 = cgetg(RU,t_COL);

for (i=1; i<RU; i++) arch0[i] = un; /* 0 in multiplicative notations */

if (DEBUGLEVEL) fprintferr("Computing powers for sub-factor base:\n");

powsubFB = cgetg(n, t_VEC);

F->powsubFB = cgetg(n, t_VEC);

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

{

GEN vp = (GEN)vectbase[subFB[i]];

GEN vp = (GEN)F->LP[ F->subFB[i] ];

GEN z = cgetg(3,t_VEC); z[1]=vp[1]; z[2]=vp[2];

pow = (GEN*)cgetg(a+1,t_VEC); powsubFB[i] = (long)pow;

pow = (GEN*)cgetg(a+1,t_VEC); F->powsubFB[i] = (long)pow;

pow[1]=cgetg(3,t_VEC);

pow[1][1] = (long)z;

pow[1][2] = (long)arch0;

Line 215 powsubFBgen(GEN nf,GEN subFB,long a,long prec)

Line 248 powsubFBgen(GEN nf,GEN subFB,long a,long prec)

if (DEBUGLEVEL) msgtimer("powsubFBgen");

}

/* Compute FB, numFB, idealbase, vectbase, numideal.

/* Compute FB, LV, iLP + KC*. Reset perm

* n2: bound for norm of tested prime ideals (includes be_honest())

* n : bound for prime ideals used to build relations (Bach cst) ( <= n2 )

* n : bound for p, such that P|p (NP <= n2) used to build relations

* Return prod_{p<=n2} (1-1/p) / prod_{Norm(P)<=n2} (1-1/Norm(P)),

* close to residue of zeta_K at 1 = 2^r1 (2pi)^r2 h R / (w D)

* close to residue of zeta_K at 1 = 2^r1 (2pi)^r2 h R / (w D) */

static GEN

FBgen(GEN nf,long n2,long n)

FBgen(FB_t *F, GEN nf,long n2,long n)

{

byteptr delta=diffptr;

byteptr delta = diffptr;

long KC2,i,j,k,p,lon,ip,nor, N = degpol(nf[1]);

long i, p, ip;

GEN p2,p1,NormP,lfun;

GEN prim, Res;

long prim[] = { evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),0 };

numFB = cgetg(n2+1,t_VECSMALL);

if (maxprime() < n2) err(primer1);

FB = cgetg(n2+1,t_VECSMALL);

F->FB = cgetg(n2+1, t_VECSMALL);

numideal = cgetg(n2+1,t_VECSMALL);

F->iLP = cgetg(n2+1, t_VECSMALL);

idealbase= (GEN*)cgetg(n2+1,t_VEC);

F->LV = (GEN*)new_chunk(n2+1);

lfun=realun(DEFAULTPREC);

Res = realun(DEFAULTPREC);

p=*delta++; i=0; ip=0; KC=0;

prim = icopy(gun);

while (p<=n2)

i = ip = 0;

F->KC = F->KCZ = 0;

for (p = 0;;) /* p <= n2 */

{

long av = avma, av1;

gpmem_t av = avma, av1;

if (DEBUGLEVEL>=2) { fprintferr(" %ld",p); flusherr(); }

long k, l;

prim[2] = p; p1 = primedec(nf,prim); lon=lg(p1);

GEN P, a, b;

av1 = avma;

divrsz(mulsr(p-1,lfun),p,lfun);

NEXT_PRIME_VIADIFF(p, delta);

if (itos(gmael(p1,1,4)) == N) /* p inert */

if (!F->KC && p > n) { F->KCZ = i; F->KC = ip; }

if (p > n2) break;

if (DEBUGLEVEL>1) { fprintferr(" %ld",p); flusherr(); }

prim[2] = p; P = primedec(nf,prim); l = lg(P);

/* a/b := (p-1)/p * prod_{P|p, NP <= n2} NP / (NP-1) */

av1 = avma; a = b = NULL;

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

{

NormP = gpowgs(prim,N);

GEN NormP = powgi(prim, gmael(P,k,4));

if (!is_bigint(NormP) && (nor=NormP[2]) <= n2)

long nor;

divrsz(mulsr(nor,lfun),nor-1, lfun);

if (is_bigint(NormP) || (nor=NormP[2]) > n2) break;

avma = av1;

}

else

{

numideal[p]=ip;

i++; numFB[p]=i; FB[i]=p;

for (k=1; k<lon; k++,ip++)

{

NormP = powgi(prim,gmael(p1,k,4));

if (is_bigint(NormP) || (nor=NormP[2]) > n2) break;

divrsz(mulsr(nor,lfun),nor-1, lfun);

if (a) { a = mului(nor, a); b = mului(nor-1, b); }

}

else { a = utoi(nor / p); b = utoi((nor-1) / (p-1)); }

/* keep all ideals with Norm <= n2 */

avma = av1;

if (k == lon)

setisclone(p1); /* flag it: all prime divisors in FB */

else

{ setlg(p1,k); p1 = gerepilecopy(av,p1); }

idealbase[i] = p1;

}

if (!*delta) err(primer1);

if (a) affrr(divri(mulir(a,Res),b), Res);

p += *delta++;

else affrr(divrs(mulsr(p-1,Res),p), Res);

if (KC == 0 && p>n) { KCZ=i; KC=ip; }

avma = av1;

}

if (l == 2 && itos(gmael(P,1,3)) == 1) continue; /* p inert */

if (!KC) return NULL;

KCZ2=i; KC2=ip; MAXRELSUP = min(50,4*KC);

setlg(FB,KCZ2);

/* keep non-inert ideals with Norm <= n2 */

setlg(numFB,KCZ2);

if (k == l)

setlg(numideal,KCZ2);

setisclone(P); /* flag it: all prime divisors in FB */

setlg(idealbase,KCZ2);

else

vectbase=cgetg(KC+1,t_COL);

{ setlg(P,k); P = gerepilecopy(av,P); }

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

F->FB[++i]= p;

{

F->LV[p] = P;

p1 = idealbase[i]; k=lg(p1);

F->iLP[p] = ip; ip += k-1;

p2 = vectbase + numideal[FB[i]];

for (j=1; j<k; j++) p2[j]=p1[j];

}

if (! F->KC) return NULL;

setlg(F->FB, F->KCZ+1); F->KCZ2 = i;

if (DEBUGLEVEL)

{

if (DEBUGLEVEL>1) fprintferr("\n");

if (DEBUGLEVEL>6)

{

fprintferr("########## FACTORBASE ##########\n\n");

fprintferr("KC2=%ld, KC=%ld, KCZ=%ld, KCZ2=%ld, MAXRELSUP=%ld\n",

fprintferr("KC2=%ld, KC=%ld, KCZ=%ld, KCZ2=%ld\n",

KC2, KC, KCZ, KCZ2, MAXRELSUP);

ip, F->KC, F->KCZ, F->KCZ2);

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

for (i=1; i<=F->KCZ; i++) fprintferr("++ LV[%ld] = %Z",i,F->LV[F->FB[i]]);

fprintferr("++ idealbase[%ld] = %Z",i,idealbase[i]);

}

msgtimer("factor base");

}

return lfun;

F->perm = NULL; return Res;

}

/* can we factor I / m ? (m pseudo minimum, computed in ideallllredpart1) */

/* SMOOTH IDEALS */

static long

static void

factorgen(GEN nf,GEN idealvec,long kcz,long limp)

store(long i, long e)

{

long i,j,n1,ip,v,p,k,lo,ifinal;

primfact[++primfact[0]] = i; /* index */

GEN x,q,r,P,p1,listexpo;

exprimfact[primfact[0]] = e; /* exponent */

GEN I = (GEN)idealvec[1];

}

GEN m = (GEN)idealvec[2];

GEN Nm= absi( subres(gmul((GEN)nf[7],m), (GEN)nf[1]) ); /* |Nm| */

x = divii(Nm, dethnf_i(I)); /* m in I, so NI | Nm */

/* divide out x by all P|p, where x as in can_factor(). k = v_p(Nx) */

if (is_pm1(x)) { primfact[0]=0; return 1; }

static int

listexpo = new_chunk(kcz+1);

divide_p_elt(FB_t *F, long p, long k, GEN nf, GEN m)

for (i=1; ; i++)

{

GEN P, LP = F->LV[p];

long j, v, l = lg(LP), ip = F->iLP[p];

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

{

p=FB[i]; q=dvmdis(x,p,&r);

P = (GEN)LP[j];

for (k=0; !signe(r); k++) { x=q; q=dvmdis(x,p,&r); }

v = int_elt_val(nf, m, (GEN)P[1], (GEN)P[5], NULL); /* v_P(m) */

listexpo[i] = k;

if (!v) continue;

if (cmpis(q,p)<=0) break;

store(ip + j, v); /* v = v_P(m) > 0 */

if (i==kcz) return 0;

k -= v * itos((GEN)P[4]);

if (!k) return 1;

}

if (cmpis(x,limp) > 0) return 0;

return 0;

}

ifinal = i; lo = 0;

static int

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

divide_p_id(FB_t *F, long p, long k, GEN nf, GEN I)

{

GEN P, LP = F->LV[p];

long j, v, l = lg(LP), ip = F->iLP[p];

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

{

k = listexpo[i];

P = (GEN)LP[j];

if (k)

v = idealval(nf,I, P);

{

if (!v) continue;

p = FB[i]; p1 = idealbase[numFB[p]];

store(ip + j, v); /* v = v_P(I) > 0 */

n1 = lg(p1); ip = numideal[p];

k -= v * itos((GEN)P[4]);

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

if (!k) return 1;

{

P = (GEN)p1[j];

v = idealval(nf,I, P) - element_val2(nf,m,Nm, P);

if (v) /* hence < 0 */

{

primfact[++lo]=ip+j; expoprimfact[lo]=v;

k += v * itos((GEN)P[4]);

if (!k) break;

}

if (k) return 0;

}

if (is_pm1(x)) { primfact[0]=lo; return 1; }

return 0;

}

p = itos(x); p1 = idealbase[numFB[p]];

static int

n1 = lg(p1); ip = numideal[p];

divide_p_quo(FB_t *F, long p, long k, GEN nf, GEN I, GEN m)

for (k=1,j=1; j<n1; j++)

{

GEN P, LP = F->LV[p];

long j, v, l = lg(LP), ip = F->iLP[p];

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

{

P = (GEN)p1[j];

P = (GEN)LP[j];

v = idealval(nf,I, P) - element_val2(nf,m,Nm, P);

v = int_elt_val(nf, m, (GEN)P[1], (GEN)P[5], NULL); /* v_P(m) */

if (v)

if (!v) continue;

{

v = idealval(nf,I, P) - v;

primfact[++lo]=ip+j; expoprimfact[lo]=v;

if (!v) continue;

k += v*itos((GEN)P[4]);

store(ip + j, v); /* v = v_P(I / m) < 0 */

if (!k) { primfact[0]=lo; return 1; }

k += v * itos((GEN)P[4]);

}

if (!k) return 1;

}

return 0;

}

/* can we factor x ? Nx = norm(x) */

/* is *N > 0 a smooth rational integer wrt F ? (put the exponents in *ex) */

static long

static int

factorelt(GEN nf,GEN cbase,GEN x,GEN Nx,long kcz,long limp)

smooth_int(FB_t *F, GEN *N, GEN *ex)

{

long i,j,n1,ip,v,p,k,lo,ifinal;

GEN q, r, FB = F->FB;

GEN q,r,P,p1,listexpo;

const long KCZ = F->KCZ;

const long limp = FB[KCZ]; /* last p in FB */

long i, p, k;

if (is_pm1(Nx)) { primfact[0]=0; return 1; }

*ex = new_chunk(KCZ+1);

listexpo = new_chunk(kcz+1);

for (i=1; ; i++)

{

p=FB[i]; q=dvmdis(Nx,p,&r);

p = FB[i]; q = dvmdis(*N,p,&r);

for (k=0; !signe(r); k++) { Nx=q; q=dvmdis(Nx,p,&r); }

for (k=0; !signe(r); k++) { *N = q; q = dvmdis(*N, p, &r); }

listexpo[i] = k;

(*ex)[i] = k;

if (cmpis(q,p)<=0) break;

if (cmpis(q,p) <= 0) break;

if (i==kcz) return 0;

if (i == KCZ) return 0;

}

if (cmpis(Nx,limp) > 0) return 0;

(*ex)[0] = i;

return (cmpis(*N,limp) <= 0);

}

if (cbase) x = gmul(cbase,x);

static int

ifinal=i; lo = 0;

divide_p(FB_t *F, long p, long k, GEN nf, GEN I, GEN m)

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

{

if (!m) return divide_p_id (F,p,k,nf,I);

k = listexpo[i];

if (!I) return divide_p_elt(F,p,k,nf,m);

if (k)

return divide_p_quo(F,p,k,nf,I,m);

{

}

p = FB[i]; p1 = idealbase[numFB[p]];

n1 = lg(p1); ip = numideal[p];

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

{

P = (GEN)p1[j];

v = int_elt_val(nf,x,(GEN)P[1],(GEN)P[5], NULL, k);

if (v)

{

primfact[++lo]=ip+j; expoprimfact[lo]=v;

k -= v * itos((GEN)P[4]);

if (!k) break;

}

if (k) return 0;

}

if (is_pm1(Nx)) { primfact[0]=lo; return 1; }

p = itos(Nx); p1 = idealbase[numFB[p]];

n1 = lg(p1); ip = numideal[p];

/* Let x = m if I == NULL,

for (k=1,j=1; j<n1; j++)

* I if m == NULL,

* I/m otherwise.

* Can we factor x ? N = Norm x > 0 */

static long

can_factor(FB_t *F, GEN nf, GEN I, GEN m, GEN N)

{

GEN ex;

long i;

primfact[0] = 0;

if (is_pm1(N)) return 1;

if (!smooth_int(F, &N, &ex)) return 0;

for (i=1; i<=ex[0]; i++)

if (ex[i] && !divide_p(F, F->FB[i], ex[i], nf, I, m)) return 0;

return is_pm1(N) || divide_p(F, itos(N), 1, nf, I, m);

}

/* can we factor I/m ? [m in I from pseudomin] */

static long

factorgen(FB_t *F, GEN nf, GEN I, GEN m)

{

GEN Nm = absi( subres(gmul((GEN)nf[7],m), (GEN)nf[1]) ); /* |Nm| */

GEN N = diviiexact(Nm, dethnf_i(I)); /* N(m / I) */

return can_factor(F, nf, I, m, N);

}

/* FUNDAMENTAL UNITS */

/* a, m real. Return (Re(x) + a) + I * (Im(x) % m) */

static GEN

addRe_modIm(GEN x, GEN a, GEN m)

{

GEN re, im, z;

if (typ(x) == t_COMPLEX)

{

P = (GEN)p1[j];

re = gadd((GEN)x[1], a);

v = int_elt_val(nf,x,(GEN)P[1],(GEN)P[5], NULL, k);

im = gmod((GEN)x[2], m);

if (v)

if (gcmp0(im)) z = re;

else

{

primfact[++lo]=ip+j; expoprimfact[lo]=v;

z = cgetg(3,t_COMPLEX);

k -= v*itos((GEN)P[4]);

z[1] = (long)re;

if (!k) { primfact[0]=lo; return 1; }

z[2] = (long)im;

}

return 0;

else

z = gadd(x, a);

return z;

}

/* clean archimedean components */

static GEN

cleancol(GEN x,long N,long PRECREG)

cleanarch(GEN x, long N, long prec)

{

long i,j,av,tetpil,tx=typ(x),R1,RU;

long i, R1, RU, tx = typ(x);

GEN s,s2,re,pi4,im,y;

GEN s, y, pi2;

if (tx==t_MAT)

if (tx == t_MAT)

{

y=cgetg(lg(x),tx);

y = cgetg(lg(x), tx);

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

for (i=1; i < lg(x); i++)

y[j]=(long)cleancol((GEN)x[j],N,PRECREG);

y[i] = (long)cleanarch((GEN)x[i], N, prec);

return y;

}

if (!is_vec_t(tx)) err(talker,"not a vector/matrix in cleancol");

if (!is_vec_t(tx)) err(talker,"not a vector/matrix in cleanarch");

av = avma; RU=lg(x)-1; R1 = (RU<<1)-N;

RU = lg(x)-1; R1 = (RU<<1)-N;

re = greal(x); s=(GEN)re[1]; for (i=2; i<=RU; i++) s=gadd(s,(GEN)re[i]);

s = gdivgs(sum(greal(x), 1, RU), -N); /* -log |norm(x)| / N */

im = gimag(x); s = gdivgs(s,-N);

y = cgetg(RU+1,tx);

s2 = (N>R1)? gmul2n(s,1): NULL;

pi2 = Pi2n(1, prec);

pi4=gmul2n(mppi(PRECREG),2);

for (i=1; i<=R1; i++) y[i] = (long)addRe_modIm((GEN)x[i], s, pi2);

tetpil=avma; y=cgetg(RU+1,tx);

if (i <= RU)

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

{

GEN p1=cgetg(3,t_COMPLEX); y[i]=(long)p1;

GEN pi4 = Pi2n(2, prec), s2 = gmul2n(s, 1);

p1[1] = ladd((GEN)re[i], (i<=R1)?s:s2);

for ( ; i<=RU; i++) y[i] = (long)addRe_modIm((GEN)x[i], s2, pi4);

p1[2] = lmod((GEN)im[i], pi4);

}

return gerepile(av,tetpil,y);

return y;

}

#define RELAT 0

enum { RELAT, LARGE, PRECI };

#define LARGE 1

#define PRECI 2

static GEN

not_given(long av, long flun, long reason)

not_given(gpmem_t av, long fl, long reason)

{

if (labs(flun)==2)

if (! (fl & nf_FORCE))

{

char *s;

switch(reason)

{

case RELAT: s = "not enough relations for fundamental units"; break;

case LARGE: s = "fundamental units too large"; break;

case PRECI: s = "insufficient precision for fundamental units"; break;

default: s = "unknown problem with fundamental units";

}

err(warner,"%s, not given",s);

}

avma=av; return cgetg(1,t_MAT);

avma = av; return cgetg(1,t_MAT);

}

/* check whether exp(x) will get too big */

static long

expgexpo(GEN x)

{

long i,j,e, E = -HIGHEXPOBIT;

long i,j,e, E = - (long)HIGHEXPOBIT;

GEN p1;

for (i=1; i<lg(x); i++)

Line 527 gauss_realimag(GEN x, GEN y)

Line 563 gauss_realimag(GEN x, GEN y)

y = split_realimag(y,r1,r2); return gauss(M, y);

}

static GEN

GEN

getfu(GEN nf,GEN *ptxarch,GEN reg,long flun,long *pte,long prec)

getfu(GEN nf,GEN *ptA,GEN reg,long fl,long *pte,long prec)

{

long av = avma,e,i,j,R1,RU,N=degpol(nf[1]);

long e, i, j, R1, RU, N=degpol(nf[1]);

GEN p1,p2,u,y,matep,s,xarch,vec;

gpmem_t av = avma;

GEN *gptr[2];

GEN p1,p2,u,y,matep,s,A,vec;

if (DEBUGLEVEL) fprintferr("\n#### Computing fundamental units\n");

R1 = itos(gmael(nf,2,1)); RU = (N+R1)>>1;

if (RU==1) { *pte=BIGINT; return cgetg(1,t_VEC); }

*pte = 0; xarch = *ptxarch;

*pte = 0; A = *ptA;

if (gexpo(reg) < -8) return not_given(av,flun,RELAT);

matep = cgetg(RU,t_MAT);

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

{

s = gzero; for (i=1; i<=RU; i++) s = gadd(s,greal(gcoeff(xarch,i,j)));

s = gzero; for (i=1; i<=RU; i++) s = gadd(s,greal(gcoeff(A,i,j)));

s = gdivgs(s, -N);

p1=cgetg(RU+1,t_COL); matep[j]=(long)p1;

for (i=1; i<=R1; i++) p1[i] = ladd(s, gcoeff(xarch,i,j));

for (i=1; i<=R1; i++) p1[i] = ladd(s, gcoeff(A,i,j));

for ( ; i<=RU; i++) p1[i] = ladd(s, gmul2n(gcoeff(xarch,i,j),-1));

for ( ; i<=RU; i++) p1[i] = ladd(s, gmul2n(gcoeff(A,i,j),-1));

}

if (prec <= 0) prec = gprecision(xarch);

if (prec <= 0) prec = gprecision(A);

u = lllintern(greal(matep),1,prec);

u = lllintern(greal(matep),100,1,prec);

if (!u) return not_given(av,flun,PRECI);

if (!u) return not_given(av,fl,PRECI);

p1 = gmul(matep,u);

if (expgexpo(p1) > 20) { *pte = BIGINT; return not_given(av,flun,LARGE); }

if (expgexpo(p1) > 20) { *pte = BIGINT; return not_given(av,fl,LARGE); }

matep = gexp(p1,prec);

y = grndtoi(gauss_realimag(nf,matep), &e);

*pte = -e;

if (e>=0) return not_given(av,flun,PRECI);

if (e >= 0) return not_given(av,fl,PRECI);

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

if (!gcmp1(idealnorm(nf, (GEN)y[j]))) break;

if (j < RU) { *pte = 0; return not_given(av,flun,PRECI); }

if (j < RU) { *pte = 0; return not_given(av,fl,PRECI); }

xarch = gmul(xarch,u);

A = gmul(A,u);

/* y[i] are unit generators. Normalize: smallest L2 norm + lead coeff > 0 */

y = gmul((GEN)nf[7], y);

vec = cgetg(RU+1,t_COL); p2 = mppi(prec);

vec = cgetg(RU+1,t_COL);

p1 = pureimag(p2);

p1 = PiI2n(0,prec); for (i=1; i<=R1; i++) vec[i] = (long)p1;

p2 = pureimag(gmul2n(p2,1));

p2 = PiI2n(1,prec); for ( ; i<=RU; i++) vec[i] = (long)p2;

for (i=1; i<=R1; i++) vec[i]=(long)p1;

for ( ; i<=RU; i++) vec[i]=(long)p2;

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

{

p1 = (GEN)y[j]; p2 = ginvmod(p1, (GEN)nf[1]);

p1 = (GEN)y[j]; p2 = QX_invmod(p1, (GEN)nf[1]);

if (gcmp(QuickNormL2(p2,DEFAULTPREC),

QuickNormL2(p1,DEFAULTPREC)) < 0)

{

xarch[j] = lneg((GEN)xarch[j]);

A[j] = lneg((GEN)A[j]);

p1 = p2;

}

if (gsigne(leading_term(p1)) < 0)

{

xarch[j] = ladd((GEN)xarch[j], vec);

A[j] = ladd((GEN)A[j], vec);

p1 = gneg(p1);

}

y[j] = (long)p1;

}

if (DEBUGLEVEL) msgtimer("getfu");

*ptxarch=xarch; gptr[0]=ptxarch; gptr[1]=&y;

gerepileall(av,2, &A, &y);

gerepilemany(av,gptr,2); return y;

*ptA = A; return y;

}

#undef RELAT

#undef LARGE

#undef PRECI

GEN

buchfu(GEN bnf)

{

long av = avma, c;

long c;

GEN nf,xarch,reg,res, y = cgetg(3,t_VEC);

gpmem_t av = avma;

GEN nf,A,reg,res, y = cgetg(3,t_VEC);

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

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

res = (GEN)bnf[8]; reg = (GEN)res[2];

if (lg(res)==7 && lg(res[5])==lg(nf[6])-1)

{

y[1] = lcopy((GEN)res[5]);

y[2] = lcopy((GEN)res[6]); return y;

}

y[1] = (long)getfu(nf,&xarch,reg,2,&c,0);

y[1] = (long)getfu(nf, &A, reg, nf_UNITS, &c, 0);

y[2] = lstoi(c); return gerepilecopy(av, y);

}

Line 665 get_norm_fact(GEN gen, GEN ex, GEN *pd)

Line 695 get_norm_fact(GEN gen, GEN ex, GEN *pd)

*pd = d; return N;

}

/* try to split ideal / (some integer) on FB */

static GEN

static long

get_pr_lists(GEN FB, long N, int list_pr)

factorgensimple(GEN nf,GEN ideal)

{

long N,i,v,av1 = avma,lo, L = lg(vectbase);

GEN pr, L;

GEN x;

long i, l = lg(FB), p, pmax;

if (typ(ideal) != t_MAT) ideal = (GEN)ideal[1]; /* idele */

x = dethnf_i(ideal);

N = lg(ideal)-1;

if (gcmp1(x)) { avma=av1; primfact[0]=0; return 1; }

for (lo=0, i=1; i<L; i++)

{

GEN q,y, P=(GEN)vectbase[i], p=(GEN)P[1];

long vx0, vx, sum_ef, e,f,lo0,i0;

vx = pvaluation(x,p,&y);

pmax = 0;

if (!vx) continue;

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

/* p | x = Nideal */

{

vx0 = vx; sum_ef = 0; lo0 =lo; i0 = i;

pr = (GEN)FB[i]; p = itos((GEN)pr[1]);

for(;;)

if (p > pmax) pmax = p;

}

L = cgetg(pmax+1, t_VEC);

for (p=1; p<=pmax; p++) L[p] = 0;

if (list_pr)

{

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

{

e = itos((GEN)P[3]);

pr = (GEN)FB[i]; p = itos((GEN)pr[1]);

f = itos((GEN)P[4]); sum_ef += e * f;

if (!L[p]) L[p] = (long)cget1(N+1, t_VEC);

v = idealval(nf,ideal,P);

appendL((GEN)L[p], pr);

if (v)

{

lo++; primfact[lo]=i; expoprimfact[lo]=v;

vx -= v * f; if (!vx) break;

}

i++; if (i == L) break;

P = (GEN)vectbase[i]; q = (GEN)P[1];

if (!egalii(p,q)) break;

}

if (vx)

for (p=1; p<=pmax; p++)

{ /* divisible by P | p not in FB */

if (L[p]) L[p] = (long)gen_sort((GEN)L[p],0,cmp_prime_over_p);

long k,l,n;

}

k = N - sum_ef;

else

if (vx % k) break;

{

k = vx / k; /* ideal / p^k may factor on FB */

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

for (l = lo0+1; l <= lo; l++)

{

pr = (GEN)FB[i]; p = itos((GEN)pr[1]);

P = (GEN)vectbase[primfact[l]];

if (!L[p]) L[p] = (long)cget1(N+1, t_VECSMALL);

expoprimfact[l] -= k * itos((GEN)P[3]); /* may become 0 */

appendL((GEN)L[p], (GEN)i);

}

n = lo0+1;

for (l = i0; l < i; l++) /* update exponents for other P | p */

{

if (n <= lo && primfact[n] == l) { n++; continue; }

lo++; primfact[lo] = l; P = (GEN)vectbase[l];

expoprimfact[lo] = - k * itos((GEN)P[3]);

}

/* check that ideal / p^k / (tentative factorisation) is integral */

{

GEN z = ideal;

for (l = lo0+1; l <= lo; l++)

{

GEN n = stoi(- expoprimfact[l]);

z = idealmulpowprime(nf, z, (GEN)vectbase[primfact[l]], n);

}

z = gdiv(z, gpowgs(p, k));

if (!gcmp1(denom(z))) { avma=av1; return 0; }

ideal = z;

}

if (gcmp1(y)) { avma=av1; primfact[0]=lo; return 1; }

x = y;

}

avma=av1; return 0;

return L;

}

static void

/* recover FB, LV, iLP, KCZ from Vbase */

add_to_fact(long l, long v, long e)

static GEN

recover_partFB(FB_t *F, GEN Vbase, long N)

{

long i,lo;

GEN FB, LV, iLP, L = get_pr_lists(Vbase, N, 0);

if (!e) return;

long l = lg(L), p, ip, i;

for (i=1; i<=l && primfact[i] < v; i++)/*empty*/;

if (i <= l && primfact[i] == v)

i = ip = 0;

expoprimfact[i] += e;

FB = cgetg(l, t_VECSMALL);

else

iLP= cgetg(l, t_VECSMALL);

LV = cgetg(l, t_VEC);

#if 1 /* for debugging */

for (p=1;p<l;p++) FB[p]=iLP[p]=LV[p]=0;

#endif

for (p = 2; p < l; p++)

{

lo = ++primfact[0];

if (!L[p]) continue;

primfact[lo] = v;

FB[++i] = p;

expoprimfact[lo] = e;

LV[p] = (long)vecextract_p(Vbase, (GEN)L[p]);

iLP[p]= ip; ip += lg(L[p])-1;

}

F->KCZ = i;

F->FB = FB; setlg(FB, i+1);

F->LV = (GEN*)LV;

F->iLP= iLP; return L;

}

static GEN

init_famat(GEN x)

{

GEN y = cgetg(3, t_VEC);

y[1] = (long)x;

y[2] = lgetg(1,t_MAT); return y;

}

/* add v^e to factorization */

static void

init_sub(long l, GEN perm, GEN *v, GEN *ex)

add_to_fact(long v, long e)

{

long i;

long i, l = primfact[0];

*v = cgetg(l,t_VECSMALL);

for (i=1; i<=l && primfact[i] < v; i++)/*empty*/;

*ex= cgetg(l,t_VECSMALL);

if (i <= l && primfact[i] == v) exprimfact[i] += e; else store(v, e);

for (i=1; i<l; i++) (*v)[i] = itos((GEN)perm[i]);

}

/* factor x = x0[1] on vectbase (modulo principal ideals).

/* L (small) list of primes above the same p including pr. Return pr index */

* We may have x0[2] = NULL (principal part not wanted), in which case,

static int

* don't compute archimedean parts */

pr_index(GEN L, GEN pr)

{

long j, l = lg(L);

GEN al = (GEN)pr[2];

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

if (gegal(al, gmael(L,j,2))) return j;

err(bugparier,"codeprime");

return 0; /* not reached */

}

static long

Vbase_to_FB(FB_t *F, GEN pr)

{

long p = itos((GEN)pr[1]);

return F->iLP[p] + pr_index(F->LV[p], pr);

}

/* return famat y (principal ideal) such that x / y is smooth [wrt Vbase] */

static GEN

split_ideal(GEN nf, GEN x0, long prec, GEN vperm)

SPLIT(FB_t *F, GEN nf, GEN x, GEN Vbase)

{

GEN p1,vdir,id,z,v,ex,y, x = (GEN)x0[1];

GEN vdir, id, z, ex, y, x0;

long nbtest_lim,nbtest,bou,i,ru, lgsub;

long nbtest_lim, nbtest, bou, i, ru, lgsub;

int flag = (gexpo(gcoeff(x,1,1)) < 100);

y = x0;

/* try without reduction if x is small */

if (flag && factorgensimple(nf,y)) return y;

if (flag && can_factor(F, nf, x, NULL, dethnf_i(x))) return NULL;

y = ideallllred(nf,x0,NULL,prec);

/* if reduction fails (y scalar), do not retry can_factor */

if (flag && ((!x0[2] && gegal((GEN)y[1], (GEN)x[1]))

y = idealred_elt(nf,x);

|| (x0[2] && gcmp0((GEN)y[2])))) flag = 0; /* y == x0 */

if ((!flag || !isnfscalar(y)) && factorgen(F, nf, x, y)) return y;

if (flag && factorgensimple(nf,y)) return y;

z = init_idele(nf); ru = lg(z[2]);

/* reduce in various directions */

if (!x0[2]) { z[2] = 0; x0 = x; } /* stop cheating */

ru = lg(nf[6]);

vdir = cgetg(ru,t_VEC);

vdir = cgetg(ru,t_VECSMALL);

for (i=2; i<ru; i++) vdir[i]=zero;

for (i=2; i<ru; i++) vdir[i]=0;

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

{

vdir[i]=lstoi(10);

vdir[i] = 10;

y = ideallllred(nf,x0,vdir,prec);

y = ideallllred_elt(nf,x,vdir);

if (factorgensimple(nf,y)) return y;

if (factorgen(F, nf, x, y)) return y;

vdir[i]=zero;

vdir[i] = 0;

}

nbtest = 0; nbtest_lim = (ru-1) << 2; lgsub = 3;

init_sub(lgsub, vperm, &v, &ex);

/* tough case, multiply by random products */

lgsub = 3;

ex = cgetg(lgsub, t_VECSMALL);

z = init_famat(NULL);

x0 = init_famat(x);

nbtest = 1; nbtest_lim = 4;

for(;;)

{

int non0 = 0;

gpmem_t av = avma;

if (DEBUGLEVEL>2) fprintferr("# ideals tried = %ld\n",nbtest);

id = x0;

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

{

ex[i] = mymyrand() >> randshift;

if (ex[i])

{ /* don't let id become too large as lgsub gets bigger: avoid

{ /* avoid prec pb: don't let id become too large as lgsub increases */

prec problems */

if (id != x0) id = ideallllred(nf,id,NULL,0);

if (non0) id = ideallllred(nf,id,NULL,prec);

z[1] = Vbase[i];

non0++;

id = idealmulh(nf, id, idealpowred(nf,z,stoi(ex[i]),0));

z[1]=vectbase[v[i]]; p1=idealpowred(nf,z,stoi(ex[i]),prec);

id = idealmulh(nf,id,p1);

}

if (id == x0) continue;

for (i=1; i<ru; i++) vdir[i] = lstoi(mymyrand() >> randshift);

for (i=1; i<ru; i++) vdir[i] = mymyrand() >> randshift;

for (bou=1; bou<ru; bou++)

{

if (bou>1)

y = ideallllred_elt(nf, (GEN)id[1], vdir);

if (factorgen(F, nf, (GEN)id[1], y))

{

for (i=1; i<ru; i++) vdir[i]=zero;

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

vdir[bou]=lstoi(10);

if (ex[i]) add_to_fact(Vbase_to_FB(F,(GEN)Vbase[i]), -ex[i]);

return famat_mul((GEN)id[2], y);

}

nbtest++;

for (i=1; i<ru; i++) vdir[i] = 0;

y = ideallllred(nf,id,vdir,prec);

vdir[bou] = 10;

if (DEBUGLEVEL>2)

fprintferr("nbtest = %ld, ideal = %Z\n",nbtest,y[1]);

if (factorgensimple(nf,y))

{

long l = primfact[0];

for (i=1; i<lgsub; i++) add_to_fact(l,v[i],-ex[i]);

return y;

}

if (nbtest > nbtest_lim)

avma = av;

if (++nbtest > nbtest_lim)

{

nbtest = 0;

if (lgsub < 7)

if (++lgsub < 7)

{

lgsub++; nbtest_lim <<= 2;

nbtest_lim <<= 1;

init_sub(lgsub, vperm, &v, &ex);

ex = cgetg(lgsub, t_VECSMALL);

}

else nbtest_lim = VERYBIGINT; /* don't increase further */

if (DEBUGLEVEL)

if (DEBUGLEVEL) fprintferr("SPLIT: increasing factor base [%ld]\n",lgsub);

fprintferr("split_ideal: increasing factor base [%ld]\n",lgsub);

}

static void

static GEN

get_split_expo(GEN xalpha, GEN yalpha, GEN vperm)

split_ideal(GEN nf, GEN x, GEN Vbase)

{

long i, colW = lg(xalpha)-1;

FB_t F;

GEN vinvperm = new_chunk(lg(vectbase));

GEN L = recover_partFB(&F, Vbase, lg(x)-1);

for (i=1; i<lg(vectbase); i++) vinvperm[itos((GEN)vperm[i])]=i;

GEN y = SPLIT(&F, nf, x, Vbase);

long p,j, i, l = lg(F.FB);

p = j = 0; /* -Wall */

for (i=1; i<=primfact[0]; i++)

{ /* decode index C = ip+j --> (p,j) */

long q,k,t, C = primfact[i];

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

{

q = F.FB[t];

k = C - F.iLP[q];

if (k <= 0) break;

p = q;

j = k;

}

primfact[i] = mael(L, p, j);

}

return y;

}

/* return sorted vectbase [sorted in bnf since version 2.2.4] */

static GEN

get_Vbase(GEN bnf)

{

GEN vectbase = (GEN)bnf[5], perm = (GEN)bnf[6], Vbase;

long i, l, tx = typ(perm);

if (tx == t_INT) return vectbase;

/* old format */

l = lg(vectbase); Vbase = cgetg(l,t_VEC);

for (i=1; i<l; i++) Vbase[i] = vectbase[itos((GEN)perm[i])];

return Vbase;

}

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

void

testprimes(GEN bnf, long bound)

{

gpmem_t av0 = avma, av;

long p,i,nbideal,k,pmax;

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

byteptr d = diffptr;

FB_t F;

if (DEBUGLEVEL>1)

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

dK= (GEN)nf[3];

f = (GEN)nf[4];

if (!gcmp1(f))

{

long k = vinvperm[primfact[i]];

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

long l = k - colW;

if (DEBUGLEVEL>1) fprintferr("**** Testing Different = %Z\n",D);

if (l <= 0) xalpha[k]=lstoi(expoprimfact[i]);

p1 = isprincipalall(bnf, D, nf_FORCE);

else yalpha[l]=lstoi(expoprimfact[i]);

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

}

/* sort factorbase for tablesearch */

fb = gen_sort((GEN)bnf[5], 0, &cmp_prime_ideal);

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

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

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

Vbase = get_Vbase(bnf);

(void)recover_partFB(&F, Vbase, degpol(nf[1]));

av = avma;

for (p = 0; p < bound; )

{

NEXT_PRIME_VIADIFF(p, d);

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

vP = primedec(bnf, stoi(p));

nbideal = lg(vP)-1;

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

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

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

{

GEN P = (GEN)vP[i];

if (DEBUGLEVEL>1) fprintferr(" Testing P = %Z\n",P);

if (cmpis(idealnorm(bnf,P), bound) >= 0)

{

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

continue;

}

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

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

else if (DEBUGLEVEL>1)

fprintferr(" is %Z\n", isprincipal(bnf,P));

else /* faster: don't compute result */

SPLIT(&F, nf, prime_to_ideal(nf,P), Vbase);

}

avma = av;

}

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

avma = av0;

}

GEN

Line 897 red_mod_units(GEN col, GEN z, long prec)

Line 1010 red_mod_units(GEN col, GEN z, long prec)

RU = lg(mat); x = cgetg(RU+1,t_COL);

for (i=1; i<RU; i++) x[i]=lreal((GEN)col[i]);

x[RU] = (long)N2;

x = lllintern(concatsp(mat,x), 1,prec);

x = lllintern(concatsp(mat,x),100, 1,prec);

if (!x) return NULL;

x = (GEN)x[RU];

if (signe(x[RU]) < 0) x = gneg_i(x);

Line 939 get_Garch(GEN nf, GEN gen, GEN clg2, long prec)

Line 1052 get_Garch(GEN nf, GEN gen, GEN clg2, long prec)

static GEN

act_arch(GEN A, GEN x)

{

GEN z, a;

GEN a;

long i,l = lg(A);

long i,l = lg(A), tA = typ(A);

if (typ(A) == t_MAT)

if (tA == t_MAT)

{

/* assume lg(x) >= l */

z = cgetg(l, t_VEC);

a = cgetg(l, t_VEC);

for (i=1; i<l; i++) z[i] = (long)act_arch((GEN)A[i], x);

for (i=1; i<l; i++) a[i] = (long)act_arch((GEN)A[i], x);

return z;

return a;

}

/* A a t_COL of t_INT. Assume lg(A)==lg(x) */

if (l==1) return cgetg(1, t_VEC);

l = lg(A); z = cgetg(l, t_VEC);

if (tA == t_VECSMALL)

if (l==1) return z;

{

a = gmul((GEN)A[1], (GEN)x[1]);

a = gmulsg(A[1], (GEN)x[1]);

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

if (signe(A[i])) a = gadd(a, gmul((GEN)A[i], (GEN)x[i]));

if (A[i]) a = gadd(a, gmulsg(A[i], (GEN)x[i]));

}

else

{ /* A a t_COL of t_INT. Assume lg(A)==lg(x) */

a = gmul((GEN)A[1], (GEN)x[1]);

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

if (signe(A[i])) a = gadd(a, gmul((GEN)A[i], (GEN)x[i]));

}

settyp(a, t_VEC); return a;

}

Line 981 isprincipalarch(GEN bnf, GEN col, GEN kNx, GEN e, GEN

Line 1101 isprincipalarch(GEN bnf, GEN col, GEN kNx, GEN e, GEN

N = degpol(nf[1]);

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

RU = (N + R1)>>1;

col = cleancol(col,N,prec); settyp(col, t_COL);

col = cleanarch(col,N,prec); settyp(col, t_COL);

if (RU > 1)

{ /* reduce mod units */

GEN u, z = init_red_mod_units(bnf,prec);

Line 1001 isprincipalarch(GEN bnf, GEN col, GEN kNx, GEN e, GEN

Line 1121 isprincipalarch(GEN bnf, GEN col, GEN kNx, GEN e, GEN

static int

fact_ok(GEN nf, GEN y, GEN C, GEN g, GEN e)

{

gpmem_t av = avma;

long i, c = lg(e);

GEN z = C? C: gun;

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

if (signe(e[i])) z = idealmul(nf, z, idealpow(nf, (GEN)g[i], (GEN)e[i]));

if (typ(z) != t_MAT) z = idealhermite(nf,z);

if (typ(y) != t_MAT) y = idealhermite(nf,y);

return gegal(y, z);

i = gegal(y, z); avma = av; return i;

}

/* assume x in HNF. cf class_group_gen for notations */

static GEN

isprincipalall0(GEN bnf, GEN x, long *ptprec, long flag)

_isprincipal(GEN bnf, GEN x, long *ptprec, long flag)

{

long i,lW,lB,e,c, prec = *ptprec;

GEN Q,xar,Wex,Bex,U,y,p1,gen,cyc,xc,ex,d,col,A;

GEN W = (GEN)bnf[1];

GEN B = (GEN)bnf[2];

GEN WB_C = (GEN)bnf[4];

GEN vperm = (GEN)bnf[6];

GEN nf = (GEN)bnf[7];

GEN clg2 = (GEN)bnf[9];

int old_format = (typ(clg2[2]) == t_MAT);

U = (GEN)clg2[1];

U = (GEN)clg2[1]; if (old_format) U = ginv(U);

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

gen = gmael3(bnf,8,1,3);

ex = cgetg(c+1,t_COL);

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

vectbase = (GEN)bnf[5]; /* needed by factorgensimple */

/* factor x */

xc = content(x); if (!gcmp1(xc)) x = gdiv(x,xc);

x = Q_primitive_part(x, &xc);

xar = split_ideal(nf,x,get_Vbase(bnf));

lW = lg(W)-1; Wex = vecsmall_const(lW, 0);

lB = lg(B)-1; Bex = vecsmall_const(lB, 0);

for (i=1; i<=primfact[0]; i++)

{

long k = primfact[i];

long l = k - lW;

if (l <= 0) Wex[k] = exprimfact[i];

else Bex[l] = exprimfact[i];

}

p1 = init_idele(nf); p1[1] = (long)x;

if (!(flag & nf_GEN)) p1[2] = 0; /* don't compute archimedean part */

xar = split_ideal(nf,p1,prec,vperm);

lW = lg(W)-1; Wex = zerocol(lW);

lB = lg(B)-1; Bex = zerocol(lB); get_split_expo(Wex,Bex,vperm);

/* x = g_W Wex + g_B Bex - [xar]

* = g_W A - [xar] + [C_B]Bex since g_W B + g_B = [C_B] */

A = gsub(Wex, gmul(B,Bex));

A = gsub(small_to_col(Wex), gmul_mati_smallvec(B,Bex));

if (old_format) U = ginv(U);

Q = gmul(U, A);

ex = cgetg(c+1,t_COL);

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

Q[i] = (long)truedvmdii((GEN)Q[i],(GEN)cyc[i],(GEN*)(ex+i));

Q[i] = (long)truedvmdii((GEN)Q[i], (GEN)cyc[i], (GEN*)(ex+i));

if (!(flag & nf_GEN)) return gcopy(ex);

if ((flag & nf_GEN_IF_PRINCIPAL))

{ if (!gcmp0(ex)) return gzero; }

else if (!(flag & (nf_GEN|nf_GENMAT)))

return gcopy(ex);

/* compute arch component of the missing principal ideal */

if (old_format)

{

GEN Garch, V = (GEN)clg2[2];

Bex = small_to_col(Bex);

p1 = c? concatsp(gmul(V,Q), Bex): Bex;

col = act_arch(p1, WB_C);

if (c)

Line 1066 isprincipalall0(GEN bnf, GEN x, long *ptprec, long fla

Line 1191 isprincipalall0(GEN bnf, GEN x, long *ptprec, long fla

{ /* g A = G Ur A + [ga]A, Ur A = D Q + R as above (R = ex)

= G R + [GD]Q + [ga]A */

GEN ga = (GEN)clg2[2], GD = (GEN)clg2[3];

col = act_arch(Bex, WB_C + lW);

if (lB) col = act_arch(Bex, WB_C + lW); else col = zerovec(1); /* nf=Q ! */

if (lW) col = gadd(col, act_arch(A, ga));

if (c) col = gadd(col, act_arch(Q, GD));

}

col = gsub(col, (GEN)xar[2]);

if (xar) col = gadd(col, famat_to_arch(nf, xar, prec));

/* find coords on Zk; Q = N (x / \prod gj^ej) = N(alpha), denom(alpha) | d */

Q = gdiv(dethnf_i(x), get_norm_fact(gen, ex, &d));

col = isprincipalarch(bnf, col, Q, gun, d, &e);

if (col && !fact_ok(nf,x, col,gen,ex)) col = NULL;

if (!col && !gcmp0(ex))

{

p1 = isprincipalfact(bnf, gen, gneg(ex), x, flag);

col = (GEN)p1[2];

e = itos((GEN)p1[3]);

}

if (!col)

{

*ptprec = prec + (e >> TWOPOTBITS_IN_LONG) + (MEDDEFAULTPREC-2);

Line 1088 isprincipalall0(GEN bnf, GEN x, long *ptprec, long fla

Line 1219 isprincipalall0(GEN bnf, GEN x, long *ptprec, long fla

e = 0;

}

y = cgetg(4,t_VEC);

if (!xc) xc = gun;

col = e? gmul(xc,col): cgetg(1,t_COL);

if (flag & nf_GEN_IF_PRINCIPAL) return col;

y[1] = lcopy(ex);

y[2] = e? lmul(xc,col): lgetg(1,t_COL);

y[2] = (long)col;

y[3] = lstoi(-e); return y;

}

Line 1097 static GEN

Line 1231 static GEN

triv_gen(GEN nf, GEN x, long c, long flag)

{

GEN y;

if (!(flag & nf_GEN)) return zerocol(c);

if (flag & nf_GEN_IF_PRINCIPAL)

return (typ(x) == t_COL)? gcopy(x): algtobasis(nf,x);

if (!(flag & (nf_GEN|nf_GENMAT))) return zerocol(c);

y = cgetg(4,t_VEC);

y[1] = (long)zerocol(c);

x = (typ(x) == t_COL)? gcopy(x): algtobasis(nf,x);

y[2] = (long)((typ(x) == t_COL)? gcopy(x): algtobasis(nf,x));

y[2] = (long)x;

y[3] = lstoi(BIGINT); return y;

}

GEN

isprincipalall(GEN bnf,GEN x,long flag)

{

long av = avma,c,pr, tx = typ(x);

long c, pr, tx = typ(x);

gpmem_t av = avma;

GEN nf;

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

Line 1118 isprincipalall(GEN bnf,GEN x,long flag)

Line 1254 isprincipalall(GEN bnf,GEN x,long flag)

err(talker,"not the same number field in isprincipal");

x = (GEN)x[2]; tx = t_POL;

}

if (tx == t_POL || tx == t_COL)

if (tx == t_POL || tx == t_COL || tx == t_INT || tx == t_FRAC)

{

if (gcmp0(x)) err(talker,"zero ideal in isprincipal");

return triv_gen(nf, x, lg(mael3(bnf,8,1,2))-1, flag);

}

x = idealhermite(nf,x);

if (lg(x)==1) err(talker,"zero ideal in isprincipal");

if (lgef(nf[1])==4)

if (degpol(nf[1]) == 1)

return gerepileupto(av, triv_gen(nf, gcoeff(x,1,1), 0, flag));

pr = prec_arch(bnf); /* precision of unit matrix */

c = getrand();

for (;;)

{

long av1 = avma;

gpmem_t av1 = avma;

GEN y = isprincipalall0(bnf,x,&pr,flag);

GEN y = _isprincipal(bnf,x,&pr,flag);

if (y) return gerepileupto(av,y);

if (DEBUGLEVEL) err(warnprec,"isprincipalall0",pr);

if (DEBUGLEVEL) err(warnprec,"isprincipal",pr);

avma = av1; bnf = bnfnewprec(bnf,pr); setrand(c);

avma = av1; bnf = bnfnewprec(bnf,pr); (void)setrand(c);

}

Line 1145 isprincipalall(GEN bnf,GEN x,long flag)

Line 1281 isprincipalall(GEN bnf,GEN x,long flag)

GEN

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

{

long av = avma, l = lg(e), i,prec,c;

long l = lg(e), i, prec, c;

gpmem_t av = avma;

GEN id,id2, nf = checknf(bnf), z = NULL; /* gcc -Wall */

int gen = flag & (nf_GEN | nf_GENMAT);

int gen = flag & (nf_GEN|nf_GENMAT);

prec = prec_arch(bnf);

if (gen)

Line 1163 isprincipalfact(GEN bnf,GEN P, GEN e, GEN C, long flag

Line 1300 isprincipalfact(GEN bnf,GEN P, GEN e, GEN C, long flag

id2 = idealpowred(bnf,z, (GEN)e[i],prec);

id = id? idealmulred(nf,id,id2,prec): id2;

}

if (id == C)

if (id == C) /* e = 0 */

{

if (!C) id = gun;

if (!C) return isprincipalall(bnf, gun, flag);

return isprincipalall(bnf, id, flag);

C = idealhermite(nf,C); id = z;

if (gen) id[1] = (long)C; else id = C;

}

c = getrand();

for (;;)

{

long av1 = avma;

gpmem_t av1 = avma;

GEN y = isprincipalall0(bnf, gen? (GEN)id[1]: id,&prec,flag);

GEN y = _isprincipal(bnf, gen? (GEN)id[1]: id,&prec,flag);

if (y)

{

if (gen && typ(y) == t_VEC)

GEN u = (GEN)y[2];

if (!gen || typ(y) != t_VEC) return gerepileupto(av,y);

if (flag & nf_GENMAT)

y[2] = (isnfscalar(u) && gcmp1((GEN)u[1]))? id[2]

:(long)arch_mul((GEN)id[2],u);

else

{

GEN u = lift_intern(basistoalg(nf, (GEN)y[2]));

u = lift_intern(basistoalg(nf, u));

if (flag & nf_GENMAT)

y[2] = (long)algtobasis(nf, gmul((GEN)id[2], u));

y[2] = (gcmp1(u)&&lg(id[2])>1)? id[2]: (long)arch_mul((GEN)id[2], u);

else

y[2] = (long)algtobasis(nf, gmul((GEN)id[2], u));

y = gcopy(y);

}

return gerepileupto(av,y);

return gerepilecopy(av, y);

}

if (flag & nf_GIVEPREC)

Line 1193 isprincipalfact(GEN bnf,GEN P, GEN e, GEN C, long flag

Line 1332 isprincipalfact(GEN bnf,GEN P, GEN e, GEN C, long flag

err(warner,"insufficient precision for generators, not given");

avma = av; return stoi(prec);

}

if (DEBUGLEVEL) err(warnprec,"isprincipalall0",prec);

if (DEBUGLEVEL) err(warnprec,"isprincipal",prec);

avma = av1; bnf = bnfnewprec(bnf,prec); setrand(c);

avma = av1; bnf = bnfnewprec(bnf,prec); (void)setrand(c);

}

GEN

isprincipal(GEN bnf,GEN x)

{

return isprincipalall(bnf,x,nf_REGULAR);

return isprincipalall(bnf,x,0);

}

GEN

Line 1222 isprincipalgenforce(GEN bnf,GEN x)

Line 1361 isprincipalgenforce(GEN bnf,GEN x)

return isprincipalall(bnf,x,nf_GEN | nf_FORCE);

}

static GEN

rational_unit(GEN x, long n, long RU)

{

GEN y;

if (!gcmp1(x) && !gcmp_1(x)) return cgetg(1,t_COL);

y = zerocol(RU);

y[RU] = (long)gmodulss((gsigne(x)>0)? 0: n>>1, n);

return y;

}

/* if x a famat, assume it is an algebraic integer (very costly to check) */

GEN

isunit(GEN bnf,GEN x)

{

long av=avma,tetpil,tx = typ(x),i,R1,RU,n;

long tx = typ(x), i, R1, RU, n, prec;

GEN p1,logunit,y,ex,nf,z,pi2_sur_w,gn,emb;

gpmem_t av = avma;

GEN p1, v, rlog, logunit, ex, nf, z, pi2_sur_w, gn, emb;

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

logunit=(GEN)bnf[3]; RU=lg(logunit);

logunit = (GEN)bnf[3]; RU = lg(logunit);

p1 = gmael(bnf,8,4); /* roots of 1 */

gn = (GEN)p1[1]; n = itos(gn);

z = (GEN)p1[2];

z = algtobasis(nf, (GEN)p1[2]);

switch(tx)

{

case t_INT: case t_FRAC: case t_FRACN:

if (!gcmp1(x) && !gcmp_1(x)) return cgetg(1,t_COL);

return rational_unit(x, n, RU);

y = zerocol(RU); i = (gsigne(x) > 0)? 0: n>>1;

y[RU] = (long)gmodulss(i, n); return y;

case t_POLMOD:

case t_MAT: /* famat */

if (!gegal((GEN)nf[1],(GEN)x[1]))

if (lg(x) != 3 || lg(x[1]) != lg(x[2]))

err(talker,"not the same number field in isunit");

err(talker, "not a factorization matrix in isunit");

x = (GEN)x[2]; /* fall through */

break;

case t_POL:

p1 = x; x = algtobasis(bnf,x); break;

case t_COL:

if (lgef(nf[1])-2 == lg(x)) { p1 = basistoalg(nf,x); break; }

if (degpol(nf[1]) != lg(x)-1)

err(talker,"not an algebraic number in isunit");

break;

default:

default: x = algtobasis(nf, x);

err(talker,"not an algebraic number in isunit");

break;

}

if (!gcmp1(denom(x))) { avma = av; return cgetg(1,t_COL); }

/* assume a famat is integral */

if (typ(p1) != t_POLMOD) p1 = gmodulcp(p1,(GEN)nf[1]);

if (tx != t_MAT && !gcmp1(denom(x))) { avma = av; return cgetg(1,t_COL); }

p1 = gnorm(p1);

if (isnfscalar(x)) return gerepileupto(av, rational_unit((GEN)x[1],n,RU));

if (!is_pm1(p1)) { avma = av; return cgetg(1,t_COL); }

R1 = itos(gmael(nf,2,1)); p1 = cgetg(RU+1,t_COL);

R1 = nf_get_r1(nf); v = cgetg(RU+1,t_COL);

for (i=1; i<=R1; i++) p1[i] = un;

for (i=1; i<=R1; i++) v[i] = un;

for ( ; i<=RU; i++) p1[i] = deux;

for ( ; i<=RU; i++) v[i] = deux;

logunit = concatsp(logunit,p1);

logunit = concatsp(logunit, v);

/* ex = fundamental units exponents */

rlog = greal(logunit);

prec = nfgetprec(nf);

for (i=1;;)

{

GEN rx, rlog = greal(logunit);

GEN logN, rx = get_arch_real(nf,x,&emb, MEDDEFAULTPREC);

long e, prec = nfgetprec(nf);

long e;

for (i=1;;)

if (rx)

{

rx = get_arch_real(nf,x,&emb, MEDDEFAULTPREC);

logN = sum(rx, 1, RU); /* log(Nx), should be ~ 0 */

if (rx)

if (gexpo(logN) > -20)

{

ex = grndtoi(gauss(rlog, rx), &e);

long p = 2 + max(1, (nfgetprec(nf)-2) / 2);

if (gcmp0((GEN)ex[RU]) && e < -4) break;

if (typ(logN) != t_REAL || gprecision(rx) > p)

{ avma = av; return cgetg(1,t_COL); } /* not a precision problem */

rx = NULL;

}

if (++i > 4) err(precer,"isunit");

if (rx)

{

ex = grndtoi(gauss(rlog, rx), &e);

if (gcmp0((GEN)ex[RU]) && e < -4) break;

}

if (i == 1)

prec = MEDDEFAULTPREC + (gexpo(x) >> TWOPOTBITS_IN_LONG);

else

{

if (i > 4) err(precer,"isunit");

prec = (prec-1)<<1;

if (DEBUGLEVEL) err(warnprec,"isunit",prec);

nf = nfnewprec(nf, prec);

}

i++;

if (DEBUGLEVEL) err(warnprec,"isunit",prec);

nf = nfnewprec(nf, prec);

}

setlg(ex, RU);

setlg(p1, RU); settyp(p1, t_VEC);

p1 = row_i(logunit,1, 1,RU-1);

for (i=1; i<RU; i++) p1[i] = coeff(logunit, 1, i);

p1 = gneg(gimag(gmul(p1,ex))); if (!R1) p1 = gmul2n(p1, -1);

p1 = gadd(garg((GEN)emb[1],DEFAULTPREC), p1);

p1 = gadd(garg((GEN)emb[1],prec), p1);

/* p1 = arg(the missing root of 1) */

pi2_sur_w = divrs(mppi(DEFAULTPREC), n>>1);

pi2_sur_w = divrs(mppi(prec), n>>1); /* 2pi / n */

p1 = ground(gdiv(p1, pi2_sur_w));

if (n > 2)

{

GEN ro = gmael(nf,6,1);

GEN ro = gmul(row(gmael(nf,5,1), 1), z);

GEN p2 = ground(gdiv(garg(poleval(z,ro), DEFAULTPREC), pi2_sur_w));

GEN p2 = ground(gdiv(garg(ro, prec), pi2_sur_w));

p1 = mulii(p1, mpinvmod(p2, gn));

}

tetpil = avma; y = cgetg(RU+1,t_COL);

ex[RU] = lmodulcp(p1, gn);

for (i=1; i<RU; i++) y[i] = lcopy((GEN)ex[i]);

setlg(ex, RU+1); return gerepilecopy(av, ex);

y[RU] = lmodulcp(p1, gn); return gerepile(av,tetpil,y);

}

GEN

signunits(GEN bnf)

{

long av,i,j,R1,RU,mun;

long i, j, R1, RU, mun;

gpmem_t av;

GEN matunit,y,p1,p2,nf,pi;

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

Line 1326 signunits(GEN bnf)

Line 1488 signunits(GEN bnf)

return y;

}

/* LLL-reduce ideal and return T2 | ideal */

/* LLL-reduce ideal and return Cholesky for T2 | ideal */

static GEN

red_ideal(GEN *ideal,GEN T2vec,GEN prvec)

red_ideal(GEN *ideal,GEN Gvec,GEN prvec)

{

long i;

for (i=1; i<lg(T2vec); i++)

for (i=1; i<lg(Gvec); i++)

{

GEN p1,base, T2 = (GEN)T2vec[i];

GEN p1,u, G = (GEN)Gvec[i];

long prec = prvec[i];

p1 = qf_base_change(T2, *ideal, 1);

p1 = gmul(G, *ideal);

base = lllgramintern(p1,100,1,prec);

u = lllintern(p1,100,1,prec);

if (base)

if (u)

{

p1 = sqred1intern(qf_base_change(p1,base,1),1);

p1 = gmul(p1, u);

if (p1) { *ideal = gmul(*ideal,base); return p1; }

p1 = sqred1_from_QR(p1, prec);

if (p1) { *ideal = gmul(*ideal,u); return p1; }

}

if (DEBUGLEVEL) err(warner, "prec too low in red_ideal[%ld]: %ld",i,prec);

}

return NULL;

}

static double

get_minkovski(long N, long R1, GEN D, GEN gborne)

{

const long prec = DEFAULTPREC;

GEN p1,p2, pi = mppi(prec);

double bound;

p1 = gsqrt(gmulsg(N,gmul2n(pi,1)),prec);

p1 = gdiv(p1, gexp(stoi(N),prec));

p1 = gmulsg(N, gpui(p1, dbltor(2./(double)N),prec));

p2 = gpui(gdivsg(4,pi), gdivgs(stoi(N-R1),N),prec);

p1 = gmul(p1,p2);

bound = gtodouble(gmul(p1, gpui(absi(D), dbltor(1./(double)N),prec)));

bound *= gtodouble(gborne);

if (DEBUGLEVEL) { fprintferr("Bound for norms = %.0f\n",bound); flusherr(); }

return bound;

}

static void

wr_rel(long *col)

{

long i;

fprintferr("\nrel = ");

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

if (col[i]) fprintferr("%ld^%ld ",i,col[i]);

fprintferr("\n");

}

static void

set_log_embed(GEN col, GEN x, long R1, long prec)

{

long i, l = lg(x);

Line 1385 set_log_embed(GEN col, GEN x, long R1, long prec)

Line 1520 set_log_embed(GEN col, GEN x, long R1, long prec)

}

static void

set_fact(GEN c)

set_fact(GEN c, GEN first_nz, long cglob)

{

long i;

c[0] = primfact[1]; /* for already_found_relation */

first_nz[cglob] = primfact[1];

for (i=1; i<=primfact[0]; i++) c[primfact[i]] = expoprimfact[i];

for (i=1; i<=primfact[0]; i++) c[primfact[i]] = exprimfact[i];

}

static void

Line 1399 unset_fact(GEN c)

Line 1534 unset_fact(GEN c)

for (i=1; i<=primfact[0]; i++) c[primfact[i]] = 0;

}

#define MAXTRY 20

/* as small_to_mat, with explicit dimension d */

static GEN

small_to_mat_i(GEN z, long d)

{

long i,j, l = lg(z);

GEN x = cgetg(l,t_MAT);

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

{

GEN c = cgetg(d+1, t_COL), cz = (GEN)z[j];

x[j] = (long)c;

for (i=1; i<=d; i++) c[i] = lstoi(cz[i]);

}

return x;

}

/* return -1 in case of precision problems. t = current number of relations */

/* return -1 in case of precision problems. t = current # of relations */

static long

small_norm_for_buchall(long cglob,GEN *mat,GEN matarch,long LIMC, long PRECREG,

small_norm_for_buchall(long cglob,GEN *mat,GEN first_nz, GEN matarch,

GEN nf,GEN gborne,long nbrelpid,GEN invp,GEN L,GEN LLLnf)

long LIMC, long PRECREG, FB_t *F,

GEN nf,long nbrelpid,GEN invp,GEN L)

{

const int maxtry_DEP = 20;

const int maxtry_FACT = 500;

const double eps = 0.000001;

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

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

ulong av = avma, av1,av2,limpile;

gpmem_t av = avma, av1, av2, limpile;

long i,j,k,noideal, nbrel = lg(mat)-1;

long j,k,noideal, nbrel = lg(mat)-1;

long alldep = 0, nbsmallnorm,nbfact,R1, N = degpol(nf[1]);

long nbsmallnorm,nbfact,R1, N = degpol(nf[1]);

GEN x,xembed,M,T2,r,cbase,invcbase,T2vec,prvec;

GEN x,xembed,M,G,r,Gvec,prvec;

if (gsigne(gborne)<=0) return cglob;

if (DEBUGLEVEL)

fprintferr("\n#### Looking for %ld relations (small norms)\n",nbrel);

xembed = NULL; /* gcc -Wall */

nbsmallnorm = nbfact = 0;

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

R1 = nf_get_r1(nf);

T2 = (GEN)LLLnf[1];

M = gmael(nf,5,1);

M = (GEN)LLLnf[2];

G = gmael(nf,5,2);

cbase =(GEN)LLLnf[3];

invcbase=(GEN)LLLnf[4];

prvec = cgetg(3,t_VECSMALL);

prvec = cgetg(3,t_VECSMALL); Gvec = cgetg(3,t_VEC);

prvec[1] = MEDDEFAULTPREC;

prvec[1] = MEDDEFAULTPREC; Gvec[1] = (long)gprec_w(G,prvec[1]);

prvec[2] = PRECREG;

prvec[2] = PRECREG; Gvec[2] = (long)G;

T2vec = cgetg(3,t_VEC);

T2vec[1] = (long)gprec_w(T2,prvec[1]);

T2vec[2] = (long)T2;

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

av1 = avma;

MINKOVSKI_BOUND = get_minkovski(N,R1,(GEN)nf[3],gborne);

for (noideal=F->KC; noideal; noideal--)

for (noideal=KC; noideal; noideal--)

{

long nbrelideal=0, dependent = 0;

gpmem_t av0 = avma;

GEN IDEAL, ideal = (GEN)vectbase[noideal];

long nbrelideal=0, dependent = 0, try_factor = 0, oldcglob = cglob;

GEN normideal = idealnorm(nf,ideal);

GEN IDEAL, ideal = (GEN)F->LP[noideal];

if (alldep > 2*N) break;

if (DEBUGLEVEL>1) fprintferr("\n*** Ideal no %ld: %Z\n", noideal, ideal);

ideal = prime_to_ideal(nf,ideal);

IDEAL = invcbase? gmul(invcbase, ideal): ideal;

IDEAL = lllint_ip(ideal,4); /* should be almost T2-reduced */

IDEAL = gmul(IDEAL, lllint(IDEAL)); /* should be almost T2-reduced */

r = red_ideal(&IDEAL,Gvec,prvec);

r = red_ideal(&IDEAL,T2vec,prvec);

if (!r) return -1; /* precision problem */

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

{

v[k]=gtodouble(gcoeff(r,k,k));

for (j=1; j<k; j++) q[j][k]=gtodouble(gcoeff(r,j,k));

if (DEBUGLEVEL>3) fprintferr("v[%ld]=%.0f ",k,v[k]);

if (DEBUGLEVEL>3) fprintferr("v[%ld]=%.4g ",k,v[k]);

}

BOUND = MINKOVSKI_BOUND * pow(gtodouble(normideal),2./(double)N);

BOUND = v[2] + v[1] * q[1][2] * q[1][2];

if (BOUND < v[1]) BOUND = v[1];

BOUND *= 2; /* at most twice larger than smallest known vector */

if (DEBUGLEVEL>1)

{

if (DEBUGLEVEL>3) fprintferr("\n");

fprintferr("BOUND = %.0f\n",BOUND); flusherr();

fprintferr("BOUND = %.4g\n",BOUND); flusherr();

}

BOUND += eps; av2=avma; limpile = stack_lim(av2,1);

BOUND *= 1 + eps; av2=avma; limpile = stack_lim(av2,1);

k = N; y[N]=z[N]=0; x[N]= (long) sqrt(BOUND/v[N]);

for(;; x[1]--)

{

ulong av3 = avma;

gpmem_t av3 = avma;

double p;

GEN col;

Line 1489 small_norm_for_buchall(long cglob,GEN *mat,GEN matarch

Line 1632 small_norm_for_buchall(long cglob,GEN *mat,GEN matarch

}

if (k==1) /* element complete */

{

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

if (y[1]<=eps) goto ENDIDEAL; /* skip all scalars: [*,0...0] */

if (ccontent(x)==1) /* primitive */

{

GEN Nx, gx = gmul_mati_smallvec(IDEAL,x);

long av4;

gpmem_t av4;

if (!isnfscalar(gx))

{

xembed = gmul(M,gx); av4 = avma; nbsmallnorm++;

if (++try_factor > maxtry_FACT) goto ENDIDEAL;

Nx = ground(norm_by_embed(R1,xembed)); setsigne(Nx, 1);

if (factorelt(nf,cbase,gx,Nx,KCZ,LIMC)) { avma = av4; break; }

if (can_factor(F, nf, NULL, gx, Nx)) { avma = av4; break; }

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

}

avma = av3;

Line 1507 small_norm_for_buchall(long cglob,GEN *mat,GEN matarch

Line 1651 small_norm_for_buchall(long cglob,GEN *mat,GEN matarch

}

cglob++; col = mat[cglob];

set_fact(col);

set_fact(col, first_nz, cglob);

if (cglob > 1 && cglob <= KC && ! addcolumntomatrix(col,invp,L))

/* make sure we get maximal rank first, then allow all relations */

if (cglob > 1 && cglob <= F->KC && ! addcolumntomatrix(col,invp,L))

{ /* Q-dependent from previous ones: forget it */

cglob--; unset_fact(col);

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

if (++dependent > MAXTRY) { alldep++; break; }

if (++dependent > maxtry_DEP) break;

avma = av3; continue;

}

/* record archimedean part */

set_log_embed((GEN)matarch[cglob], xembed, R1, PRECREG);

alldep = dependent = 0;

dependent = 0;

if (DEBUGLEVEL)

if (DEBUGLEVEL) { nbfact++; dbg_rel(cglob, mat[cglob]); }

{

if (DEBUGLEVEL==1) fprintferr("%4ld",cglob);

else { fprintferr("cglob = %ld. ",cglob); wr_rel(mat[cglob]); }

flusherr(); nbfact++;

}

if (cglob >= nbrel) goto END; /* we have enough */

if (++nbrelideal == nbrelpid) break;

Line 1535 small_norm_for_buchall(long cglob,GEN *mat,GEN matarch

Line 1675 small_norm_for_buchall(long cglob,GEN *mat,GEN matarch

}

ENDIDEAL:

invp = gerepilecopy(av1, invp);

if (cglob == oldcglob) avma = av0; else invp = gerepilecopy(av1, invp);

if (DEBUGLEVEL>1) msgtimer("for this ideal");

}

END:

if (DEBUGLEVEL)

{

fprintferr("\n");

fprintferr("\n"); msgtimer("small norm relations");

msgtimer("small norm relations");

fprintferr(" small norms gave %ld relations, rank = %ld.\n",

if (DEBUGLEVEL>1)

cglob, rank(small_to_mat_i((GEN)mat, F->KC)));

{

if (nbsmallnorm)

GEN p1,tmp=cgetg(cglob+1,t_MAT);

fprintferr(" nb. fact./nb. small norm = %ld/%ld = %.3f\n",

fprintferr("Elements of small norm gave %ld relations.\n",cglob);

fprintferr("Computing rank: "); flusherr();

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

{

p1=cgetg(KC+1,t_COL); tmp[j]=(long)p1;

for(i=1;i<=KC;i++) p1[i]=lstoi(mat[j][i]);

}

tmp = (GEN)sindexrank(tmp)[2]; k=lg(tmp)-1;

fprintferr("%ld; independent columns: %Z\n",k,tmp);

}

if(nbsmallnorm)

fprintferr("\nnb. fact./nb. small norm = %ld/%ld = %.3f\n",

nbfact,nbsmallnorm,((double)nbfact)/nbsmallnorm);

else

fprintferr("\nnb. small norm = 0\n");

}

avma = av; return cglob;

}

#undef MAXTRY

/* I assumed to be integral HNF, T2 a weighted T2 matrix */

/* I assumed to be integral HNF, G the Cholesky form of a weighted T2 matrix.

* Return an irrational m in I with T2(m) small */

static GEN

ideallllredpart1(GEN I, GEN T2, long prec)

pseudomin(GEN I, GEN G)

{

GEN y,m,idealpro;

GEN m, y = lllintern(gmul(G, I), 100,1, 0);

y = lllgramintern(qf_base_change(T2,I,1),100,1,prec+1);

if (!y) return NULL;

/* I, m, y integral */

m = gmul(I,(GEN)y[1]);

if (isnfscalar(m)) m = gmul(I,(GEN)y[2]);

if (DEBUGLEVEL>5) fprintferr("\nm = %Z\n",m);

idealpro = cgetg(3,t_VEC);

return m;

idealpro[1] = (long)I;

idealpro[2] = (long)m; /* irrational element of small T2 norm in I */

if (DEBUGLEVEL>5) fprintferr("\nidealpro = %Z\n",idealpro);

return idealpro;

}

static void

Line 1605 dbg_newrel(long jideal, long jdir, long phase, long cg

Line 1723 dbg_newrel(long jideal, long jdir, long phase, long cg

static void

dbg_cancelrel(long jideal,long jdir,long phase, long *col)

{

fprintferr("rel. cancelled. phase %ld: %ld ",phase);

fprintferr("rel. cancelled. phase %ld: ",phase);

if (DEBUGLEVEL>3) fprintferr("(jideal=%ld,jdir=%ld)",jideal,jdir);

wr_rel(col); flusherr();

}

static void

dbg_outrel(long phase,long cglob, GEN vperm,GEN *mat,GEN maarch)

dbg_outrel(long cglob, GEN *mat,GEN maarch)

{

long i,j;

gpmem_t av = avma;

GEN p1,p2;

GEN p1;

long j;

if (phase == 0)

p1 = cgetg(cglob+1, t_MAT);

for (j=1; j<=cglob; j++) p1[j] = (long)small_to_col(mat[j]);

fprintferr("\nRank = %ld\n", rank(p1));

if (DEBUGLEVEL>3)

{

ulong av = avma; p2=cgetg(cglob+1,t_MAT);

fprintferr("relations = \n");

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

for (j=1; j <= cglob; j++) wr_rel(mat[j]);

{

fprintferr("\nmatarch = %Z\n",maarch);

p1=cgetg(KC+1,t_COL); p2[j]=(long)p1;

for (i=1; i<=KC; i++) p1[i]=lstoi(mat[j][i]);

}

fprintferr("\nRank = %ld\n", rank(p2));

if (DEBUGLEVEL>3)

{

fprintferr("relations = \n");

for (j=1; j <= cglob; j++) wr_rel(mat[j]);

fprintferr("\nmatarch = %Z\n",maarch);

}

avma = av;

}

else if (DEBUGLEVEL>6)

avma = av; flusherr();

{

fprintferr("before hnfadd:\nvectbase[vperm[]] = \n");

fprintferr("[");

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

{

bruterr((GEN)vectbase[vperm[i]],'g',-1);

if (i<KC) fprintferr(",");

}

fprintferr("]~\n");

}

flusherr();

}

/* Check if we already have a column mat[l] equal to mat[s]

/* Check if we already have a column mat[i] equal to mat[s]

* General check for colinearity useless since exceedingly rare */

static long

already_found_relation(GEN *mat,long s)

already_found_relation(GEN *mat, long s, GEN first_nz)

{

GEN coll, cols = mat[s];

GEN cols = mat[s];

long l,bs;

long i, bs, l = lg(cols);

bs = 1; while (bs<=KC && !cols[bs]) bs++;

bs = 1; while (bs < l && !cols[bs]) bs++;

if (bs > KC) return s; /* zero relation */

if (bs == l) return s; /* zero relation */

for (l=s-1; l; l--)

for (i=s-1; i; i--)

{

coll = mat[l];

if (bs == first_nz[i]) /* = index of first non zero elt in cols */

if (bs == coll[0]) /* = index of first non zero elt in coll */

{

GEN coll = mat[i];

long b = bs;

do b++; while (b<=KC && cols[b] == coll[b]);

do b++; while (b < l && cols[b] == coll[b]);

if (b > KC) return l;

if (b == l) return i;

}

cols[0] = bs; return 0;

first_nz[s] = bs; return 0;

}

/* I integral ideal in HNF form */

Line 1676 static GEN

Line 1776 static GEN

remove_content(GEN I)

{

long N = lg(I)-1;

if (!gcmp1(gcoeff(I,N,N))) { GEN y=content(I); if (!gcmp1(y)) I=gdiv(I,y); }

if (!gcmp1(gcoeff(I,N,N))) I = Q_primpart(I);

return I;

}

/* if phase != 1 re-initialize static variables. If <0 return immediately */

static long

random_relation(long phase,long cglob,long LIMC,long PRECLLL,long PRECREG,

random_relation(long phase,long cglob,long LIMC,long PRECREG,long MAXRELSUP,

GEN nf,GEN subFB,GEN vecT2,GEN *mat,GEN matarch,GEN list_jideal)

GEN nf,GEN vecG,GEN *mat,GEN first_nz,GEN matarch,

GEN list_jideal, FB_t *F)

{

static long jideal, jdir;

long lim,i,av,av1,cptzer,nbT2,lgsub,r1, jlist = 1;

long i, maxcglob, cptlist, cptzer, nbG, lgsub, r1, jlist = 1;

GEN arch,col,colarch,ideal,idealpro,P,ex;

gpmem_t av, av1;

GEN arch,col,colarch,ideal,m,P,ex;

if (phase != 1) { jideal=jdir=1; if (phase<0) return 0; }

r1 = nf_get_r1(nf);

lim = lg(mat)-1; /* requested number of relations */

maxcglob = lg(mat)-1; /* requested # of relations */

nbT2 = lg(vecT2)-1;

nbG = lg(vecG)-1;

lgsub = lg(subFB); ex = cgetg(lgsub, t_VECSMALL);

lgsub = lg(F->subFB); ex = cgetg(lgsub, t_VECSMALL);

cptzer = 0;

cptzer = cptlist = 0;

if (DEBUGLEVEL && list_jideal)

fprintferr("looking hard for %Z\n",list_jideal);

P = NULL; /* gcc -Wall */

for (av = avma;;)

{

if (list_jideal && jlist < lg(list_jideal) && jdir <= nbT2)

if (list_jideal && jlist < lg(list_jideal) && jdir <= nbG)

jideal = list_jideal[jlist++];

if (!list_jideal || jdir <= nbT2)

{

jideal = list_jideal[jlist++]; cptlist = 0;

}

if (!list_jideal || jdir <= nbG)

{

avma = av;

P = prime_to_ideal(nf, (GEN)vectbase[jideal]);

P = prime_to_ideal(nf, (GEN)F->LP[jideal]);

}

else

{

if (++cptlist > 300) return -1;

}

ideal = P;

do {

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

{

ex[i] = mymyrand()>>randshift;

if (ex[i]) ideal = idealmulh(nf,ideal, gmael3(powsubFB,i,ex[i],1));

if (ex[i]) ideal = idealmulh(nf,ideal, gmael3(F->powsubFB,i,ex[i],1));

}

while (ideal == P); /* If ex = 0, try another */

ideal = remove_content(ideal);

if (phase != 1) jdir = 1; else phase = 2;

for (av1 = avma; jdir <= nbT2; jdir++, avma = av1)

for (av1 = avma; jdir <= nbG; jdir++, avma = av1)

{ /* reduce along various directions */

if (DEBUGLEVEL>2)

fprintferr("phase=%ld,jideal=%ld,jdir=%ld,rand=%ld\n",

phase,jideal,jdir,getrand());

idealpro = ideallllredpart1(ideal,(GEN)vecT2[jdir],PRECLLL);

m = pseudomin(ideal,(GEN)vecG[jdir]);

if (!idealpro) return -2;

if (!m) return -2;

if (!factorgen(nf,idealpro,KCZ,LIMC))

if (!factorgen(F,nf,ideal,m))

{

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

continue;

}

/* can factor ideal, record relation */

cglob++; col = mat[cglob];

set_fact(col); col[jideal]--;

set_fact(col, first_nz, cglob); col[jideal]--;

for (i=1; i<lgsub; i++) col[subFB[i]] -= ex[i];

for (i=1; i<lgsub; i++) col[ F->subFB[i] ] -= ex[i];

if (already_found_relation(mat, cglob))

if (already_found_relation(mat, cglob, first_nz))

{ /* Already known: forget it */

if (DEBUGLEVEL>1) dbg_cancelrel(jideal,jdir,phase,col);

cglob--; unset_fact(col); col[jideal] = 0;

for (i=1; i<lgsub; i++) col[subFB[i]] = 0;

for (i=1; i<lgsub; i++) col[ F->subFB[i] ] = 0;

if (++cptzer > MAXRELSUP)

{

Line 1755 random_relation(long phase,long cglob,long LIMC,long P

Line 1863 random_relation(long phase,long cglob,long LIMC,long P

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

if (ex[i])

{

GEN p1 = gmael3(powsubFB,i,ex[i],2);

GEN p1 = gmael3(F->powsubFB,i,ex[i],2);

arch = arch? vecmul(arch, p1): p1;

}

colarch = (GEN)matarch[cglob];

/* arch = archimedean component (MULTIPLICATIVE form) of ideal */

arch = vecdiv(arch, gmul(gmael(nf,5,1), (GEN)idealpro[2]));

arch = vecdiv(arch, gmul(gmael(nf,5,1), m));

set_log_embed(colarch, arch, r1, PRECREG);

if (DEBUGLEVEL) dbg_newrel(jideal,jdir,phase,cglob,col,colarch,lim);

if (DEBUGLEVEL) dbg_newrel(jideal,jdir,phase,cglob,col,colarch,maxcglob);

/* Need more, try next P */

if (cglob < lim) break;

if (cglob < maxcglob) break;

/* We have found enough. Return */

if (phase)

{

jdir = 1;

if (jideal == KC) jideal=1; else jideal++;

if (jideal == F->KC) jideal=1; else jideal++;

}

if (DEBUGLEVEL>2)

fprintferr("Upon exit: jideal=%ld,jdir=%ld\n",jideal,jdir);

Line 1779 random_relation(long phase,long cglob,long LIMC,long P

Line 1887 random_relation(long phase,long cglob,long LIMC,long P

}

if (!list_jideal)

{

if (jideal == KC) jideal=1; else jideal++;

if (jideal == F->KC) jideal=1; else jideal++;

}

static long

/* remark: F->KCZ changes if be_honest() fails */

be_honest(GEN nf,GEN subFB,long PRECLLL)

static int

be_honest(FB_t *F, GEN nf, long PRECLLL)

{

ulong av;

long ex, i, j, J, k, iz, nbtest, ru, lgsub = lg(F->subFB), KCZ0 = F->KCZ;

GEN MC = gmael(nf,5,2), M = gmael(nf,5,1), D = (GEN)nf[3];

GEN G, M, P, ideal, m, vdir;

long ex,i,j,J,k,iz,nbtest, lgsub = lg(subFB), ru = lg(MC);

gpmem_t av;

GEN P,ideal,idealpro, exu = cgetg(ru, t_VECSMALL), MCtw= cgetg(ru, t_MAT);

if (F->KCZ2 <= F->KCZ) return 1;

if (DEBUGLEVEL)

{

fprintferr("Be honest for primes from %ld to %ld\n", FB[KCZ+1],FB[KCZ2]);

fprintferr("Be honest for primes from %ld to %ld\n",

F->FB[ F->KCZ+1 ], F->FB[ F->KCZ2 ]);

flusherr();

}

if (!F->powsubFB) powsubFBgen(F, nf, CBUCHG+1, 0);

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

G = gprec_w(gmael(nf,5,2), PRECLLL);

ru = lg(nf[6]);

vdir = cgetg(ru, t_VECSMALL);

av = avma;

for (iz=KCZ+1; iz<=KCZ2; iz++, avma = av)

for (iz=F->KCZ+1; iz<=F->KCZ2; iz++, avma = av)

{

if (DEBUGLEVEL>1) fprintferr("%ld ", FB[iz]);

long p = F->FB[iz];

P = idealbase[numFB[FB[iz]]]; J = lg(P);

if (DEBUGLEVEL>1) fprintferr("%ld ", p);

/* if unramified, check all but 1 */

P = F->LV[p]; J = lg(P);

if (J > 1 && !divise(D, gmael(P,1,1))) J--;

/* all P|p in FB + last is unramified --> check all but last */

if (isclone(P) && itou(gmael(P,J-1,3)) == 1UL /* e = 1 */) J--;

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

{

GEN ideal0 = prime_to_ideal(nf,(GEN)P[j]);

ulong av2 = avma;

gpmem_t av2 = avma;

for(nbtest=0;;)

{

ideal = ideal0;

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

{

ex = mymyrand()>>randshift;

if (ex) ideal = idealmulh(nf,ideal,gmael3(powsubFB,i,ex,1));

if (ex) ideal = idealmulh(nf,ideal,gmael3(F->powsubFB,i,ex,1));

}

ideal = remove_content(ideal);

for (i=1; i<ru; i++) vdir[i] = mymyrand()>>randshift;

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

{

if (k==1)

m = pseudomin(ideal, computeGtwist(nf, vdir));

for (i=1; i<ru; i++) exu[i] = mymyrand()>>randshift;

if (m && factorgen(F,nf,ideal,m)) break;

else

for (i=1; i<ru; i++) vdir[i] = 0;

{

vdir[k] = 10;

for (i=1; i<ru; i++) exu[i] = 0;

exu[k] = 10;

}

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

MCtw[i] = exu[i]? lmul2n((GEN)MC[i],exu[i]<<1): MC[i];

idealpro = ideallllredpart1(ideal, mulmat_real(MCtw,M), PRECLLL);

if (idealpro && factorgen(nf,idealpro,iz-1,FB[iz-1])) break;

nbtest++; if (nbtest==200) return 0;

}

avma = av2; if (k < ru) break;

if (++nbtest > 50)

{

err(warner,"be_honest() failure on prime %Z\n", P[j]);

return 0;

}

F->KCZ++; /* SUCCESS, "enlarge" factorbase */

}

if (DEBUGLEVEL)

Line 1841 be_honest(GEN nf,GEN subFB,long PRECLLL)

Line 1957 be_honest(GEN nf,GEN subFB,long PRECLLL)

if (DEBUGLEVEL>1) fprintferr("\n");

msgtimer("be honest");

}

F->KCZ = KCZ0;

avma = av; return 1;

}

Line 1851 trunc_error(GEN x)

Line 1968 trunc_error(GEN x)

&& (expo(x)>>TWOPOTBITS_IN_LONG) + 3 > lg(x);

}

/* xarch = complex logarithmic embeddings of units (u_j) found so far */

/* A = complex logarithmic embeddings of units (u_j) found so far */

static GEN

compute_multiple_of_R(GEN xarch,long RU,long N,GEN *ptsublambda)

compute_multiple_of_R(GEN A,long RU,long N,GEN *ptlambda)

{

GEN v,mdet,mdet_t,Im_mdet,kR,sublambda,lambda,xreal;

GEN T,v,mdet,mdet_t,Im_mdet,kR,xreal,lambda;

GEN *gptr[2];

long i, zc = lg(A)-1, R1 = 2*RU - N;

long av = avma, i,j, sreg = lg(xarch)-1, R1 = 2*RU - N;

gpmem_t av = avma;

if (DEBUGLEVEL) { fprintferr("\n#### Computing regulator\n"); flusherr(); }

if (DEBUGLEVEL) fprintferr("\n#### Computing regulator multiple\n");

xreal=greal(xarch); v=cgetg(RU+1,t_COL); /* xreal = (log |sigma_i(u_j)|) */

xreal = greal(A); /* = (log |sigma_i(u_j)|) */

for (i=1; i<=R1; i++) v[i]=un;

T = cgetg(RU+1,t_COL);

for ( ; i<=RU; i++) v[i]=deux;

for (i=1; i<=R1; i++) T[i] = un;

mdet=cgetg(sreg+2,t_MAT); mdet[1]=(long)v;

for ( ; i<=RU; i++) T[i] = deux;

for (j=2; j<=sreg+1; j++) mdet[j]=xreal[j-1]; /* det(Span(mdet)) = N * R */

mdet = concatsp(xreal,T); /* det(Span(mdet)) = N * R */

i = gprecision(mdet); /* truncate to avoid "near dependant" vectors */

mdet_t = (i <= 4)? mdet: gprec_w(mdet,i-1);

Line 1874 compute_multiple_of_R(GEN xarch,long RU,long N,GEN *pt

Line 1991 compute_multiple_of_R(GEN xarch,long RU,long N,GEN *pt

Im_mdet = vecextract_p(mdet,v);

/* integral multiple of R: the cols we picked form a Q-basis, they have an

* index in the full lattice */

* index in the full lattice. Last column is T */

kR = gdivgs(det2(Im_mdet), N);

/* R > 0.2 uniformly */

if (gexpo(kR) < -3) { avma=av; return NULL; }

kR = mpabs(kR);

sublambda = cgetg(sreg+1,t_MAT);

lambda = gauss(Im_mdet,xreal); /* approximate rational entries */

lambda = gauss(Im_mdet,xreal); /* rational entries */

for (i=1; i<=zc; i++) setlg(lambda[i], RU);

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

gerepileall(av,2, &lambda, &kR);

{

*ptlambda = lambda; return kR;

GEN p1 = cgetg(RU,t_COL), p2 = (GEN)lambda[i];

sublambda[i] = (long)p1;

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

{

p1[j] = p2[j+1];

if (trunc_error((GEN)p1[j])) { *ptsublambda = NULL; return gzero; }

}

*ptsublambda = sublambda;

gptr[0]=ptsublambda; gptr[1]=&kR;

gerepilemany(av,gptr,2); return kR;

}

/* Assuming enough relations, c = Rz is close to an even integer, according

* to Dirichlet's formula. Otherwise, close to a multiple.

* Compute a tentative regulator (not a multiple this time) */

static GEN

compute_check(GEN sublambda, GEN z, GEN *parch, GEN *reg)

bestappr_noer(GEN x, GEN k)

{

long av = avma, av2, tetpil;

GEN y;

GEN p1,c,den, R = *reg; /* multiple of regulator */

CATCH(precer) { y = NULL; }

TRY { y = bestappr(x,k); } ENDCATCH;

return y;

}

/* Input:

* lambda = approximate rational entries: coords of units found so far on a

* sublattice of maximal rank (sublambda)

* *ptkR = regulator of sublambda = multiple of regulator of lambda

* Compute R = true regulator of lambda.

* If c := Rz ~ 1, by Dirichlet's formula, then lambda is the full group of

* units AND the full set of relations for the class group has been computed.

* In fact z is a very rough approximation and we only expect 0.75 < Rz < 1.5

* Output: *ptkR = R, *ptU = basis of fundamental units (in terms lambda) */

static int

compute_R(GEN lambda, GEN z, GEN *ptL, GEN *ptkR)

{

gpmem_t av = avma;

long r;

GEN L,H,D,den,R;

double c;

if (DEBUGLEVEL) { fprintferr("\n#### Computing check\n"); flusherr(); }

c = gmul(R,z);

D = gmul2n(gmul(*ptkR,z), 1); /* bound for denom(lambda) */

sublambda = bestappr(sublambda,c); den = denom(sublambda);

lambda = bestappr_noer(lambda,D);

if (gcmp(den,c) > 0)

if (!lambda)

{

if (DEBUGLEVEL) fprintferr("c = %Z\nden = %Z\n",c,den);

if (DEBUGLEVEL) fprintferr("truncation error in bestappr\n");

avma=av; return NULL;

return PRECI;

}

den = denom(lambda);

if (gcmp(den,D) > 0)

{

if (DEBUGLEVEL) fprintferr("D = %Z\nden = %Z\n",D,den);

return PRECI;

}

L = Q_muli_to_int(lambda, den);

H = hnfall_i(L, NULL, 1); r = lg(H)-1;

p1 = gmul(sublambda,den); tetpil=avma;

/* tentative regulator */

*parch = lllint(p1);

R = mpabs( gmul(*ptkR, gdiv(dethnf_i(H), gpowgs(den, r))) );

c = gtodouble(gmul(R,z)); /* should be n (= 1 if we are done) */

av2=avma; p1 = det2(gmul(sublambda,*parch));

if (DEBUGLEVEL)

affrr(mpabs(gmul(R,p1)), R); avma=av2;

{

msgtimer("bestappr/regulator");

if (DEBUGLEVEL) msgtimer("bestappr/regulator");

fprintferr("\n ***** check = %f\n",c);

*parch = gerepile(av,tetpil,*parch); return gmul(R,z);

}

if (c < 0.75 || c > 1.5) { avma = av; return RELAT; }

*ptkR = R; *ptL = L; return LARGE;

}

/* find the smallest (wrt norm) among I, I^-1 and red(I^-1) */

Line 1930 static GEN

Line 2066 static GEN

inverse_if_smaller(GEN nf, GEN I, long prec)

{

GEN J, d, dmin, I1;

int inv;

J = (GEN)I[1];

dmin = dethnf_i(J); I1 = idealinv(nf,I);

J = (GEN)I1[1]; J = gmul(J,denom(J)); I1[1] = (long)J;

d = dethnf_i(J);

d = dethnf_i(J); if (cmpii(d,dmin) < 0) {I=I1; dmin=d;}

if (cmpii(d,dmin) < 0) {inv=1; I=I1; dmin=d;}

else inv=0;

/* try reducing (often _increases_ the norm) */

I1 = ideallllred(nf,I1,NULL,prec);

J = (GEN)I1[1];

d = dethnf_i(J);

d = dethnf_i(J); if (cmpii(d,dmin) < 0) I=I1;

if (cmpii(d,dmin) < 0) {inv=1; I=I1;}

return I;

}

Line 1962 setlg_col(GEN U, long l)

Line 2095 setlg_col(GEN U, long l)

/* compute class group (clg1) + data for isprincipal (clg2) */

static void

class_group_gen(GEN nf,GEN W,GEN C,GEN vperm,GEN *ptclg1,GEN *ptclg2,long prec)

class_group_gen(GEN nf,GEN W,GEN C,GEN Vbase,long prec, GEN nf0,

GEN *ptclg1,GEN *ptclg2)

{

GEN z,G,Ga,ga,GD,cyc,X,Y,D,U,V,Ur,Ui,Uir,I,J,Vbase;

GEN z,G,Ga,ga,GD,cyc,X,Y,D,U,V,Ur,Ui,Uir,I,J;

long i,j,s,lo,lo0;

long i,j,lo,lo0;

if (DEBUGLEVEL)

{ fprintferr("\n#### Computing class group generators\n"); timer2(); }

{ fprintferr("\n#### Computing class group generators\n"); (void)timer2(); }

z = smith2(W); /* U W V = D, D diagonal, G = g Ui (G=new gens, g=old) */

D = smithall(W,&U,&V); /* UWV = D, D diagonal, G = g Ui (G=new gens, g=old) */

U = (GEN)z[1]; Ui = ginv(U);

Ui = ginv(U);

V = (GEN)z[2];

D = (GEN)z[3];

lo0 = lo = lg(D);

/* we could set lo = lg(cyc) and truncate all matrices below

* setlg_col(D && U && Y, lo) + setlg(D && V && X && Ui, lo)

Line 1987 class_group_gen(GEN nf,GEN W,GEN C,GEN vperm,GEN *ptcl

Line 2119 class_group_gen(GEN nf,GEN W,GEN C,GEN vperm,GEN *ptcl

* P invertible diagonal matrix (\pm 1) which is only implicitly defined

* G = g Uir P + [Ga], Uir = Ui + WX

* g = G P Ur + [ga], Ur = U + DY */

Vbase = cgetg(lo0,t_VEC);

if (typ(vperm) == t_VECSMALL)

for (i=1; i<lo0; i++) Vbase[i] = vectbase[vperm[i]];

else

for (i=1; i<lo0; i++) Vbase[i] = vectbase[itos((GEN)vperm[i])];

G = cgetg(lo,t_VEC);

Ga= cgetg(lo,t_VEC);

z = init_idele(nf);

z = init_famat(NULL);

if (!nf0) nf0 = nf;

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

{

GEN p1 = gcoeff(Uir,1,j);

z[1]=Vbase[1]; I = idealpowred(nf,z,p1,prec);

z[1]=Vbase[1]; I = idealpowred(nf0,z,p1,prec);

if (signe(p1)<0) I[1] = lmul((GEN)I[1],denom((GEN)I[1]));

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

{

p1 = gcoeff(Uir,i,j); s = signe(p1);

p1 = gcoeff(Uir,i,j);

if (s)

if (signe(p1))

{

z[1]=Vbase[i]; J = idealpowred(nf,z,p1,prec);

z[1]=Vbase[i]; J = idealpowred(nf0,z,p1,prec);

if (s<0) J[1] = lmul((GEN)J[1],denom((GEN)J[1]));

I = idealmulh(nf0,I,J);

I = idealmulh(nf,I,J);

I = ideallllred(nf0,I,NULL,prec);

I = ideallllred(nf,I,NULL,prec);

}

J = inverse_if_smaller(nf, I, prec);

J = inverse_if_smaller(nf0, I, prec);

if (J != I)

{ /* update wrt P */

neg_row(Y ,j); V[j] = lneg((GEN)V[j]);

neg_row(Ur,j); X[j] = lneg((GEN)X[j]);

}

G[j] = (long)J[1]; /* generator, order cyc[j] */

Ga[j]= (long)J[2];

Ga[j]= lneg(famat_to_arch(nf, (GEN)J[2], prec));

}

/* at this point Y = PY, Ur = PUr, V = VP, X = XP */

Line 2056 class_group_gen(GEN nf,GEN W,GEN C,GEN vperm,GEN *ptcl

Line 2181 class_group_gen(GEN nf,GEN W,GEN C,GEN vperm,GEN *ptcl

if (DEBUGLEVEL) msgtimer("classgroup generators");

}

/* real(MC * M), where columns a and b of MC have been multiplied by 2^20+1 */

static void

static GEN

shift_embed(GEN G, GEN Gtw, long a, long r1, long r2)

shift_t2(GEN T2, GEN M, GEN MC, long a, long b)

{

long i,j,N = lg(T2);

long j, k, l = lg(G);

GEN z, t2 = cgetg(N,t_MAT);

if (a <= r1)

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

for (j=1; j<l; j++) coeff(G,a,j) = coeff(Gtw,a,j);

else

{

t2[j] = lgetg(N,t_COL);

k = (a<<1) - r1;

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

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

{

z = mul_real(gcoeff(MC,i,a), gcoeff(M,a,j));

coeff(G,k-1,j) = coeff(Gtw,k-1,j);

if (a!=b) z = gadd(z, mul_real(gcoeff(MC,i,b), gcoeff(M,b,j)));

coeff(G,k ,j) = coeff(Gtw,k, j);

coeff(t2,j,i) = coeff(t2,i,j) = ladd(gcoeff(T2,i,j), gmul2n(z,20));

}

return t2;

}

/* G where embeddings a and b are multiplied by 2^10 */

static GEN

compute_vecT2(long RU,GEN nf)

shift_G(GEN G, GEN Gtw, long a, long b, long r1, long r2)

{

GEN vecT2, M = gmael(nf,5,1), MC = gmael(nf,5,2), T2 = gmael(nf,5,3);

GEN g = dummycopy(G);

long i,j,n = min(RU,9), ind = 1;

shift_embed(g,Gtw,a,r1,r2);

shift_embed(g,Gtw,b,r1,r2); return g;

}

vecT2=cgetg(1 + n*(n+1)/2,t_VEC);

static GEN

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

compute_vecG(GEN nf,long prec)

for (i=1; i<=j; i++) vecT2[ind++] = (long)shift_t2(T2,M,MC,i,j);

{

if (DEBUGLEVEL) msgtimer("weighted T2 matrices");

GEN vecG, Gtw, M = gmael(nf,5,1), G = gmael(nf,5,2);

return vecT2;

long r1,r2,i,j,ind, n = min(lg(M[1])-1, 9);

nf_get_sign(nf,&r1,&r2);

vecG=cgetg(1 + n*(n+1)/2,t_VEC);

if (nfgetprec(nf) > prec)

{

M = gprec_w(M,prec);

G = gprec_w(G,prec);

}

Gtw = gmul2n(G, 10);

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

for (i=1; i<=j; i++) vecG[ind++] = (long)shift_G(G,Gtw,i,j,r1,r2);

if (DEBUGLEVEL) msgtimer("weighted G matrices");

return vecG;

}

/* cf. relationrank()

Line 2129 addcolumntomatrix(GEN V, GEN invp, GEN L)

Line 2268 addcolumntomatrix(GEN V, GEN invp, GEN L)

* Starting from the identity (canonical basis of Q^n), we obtain a base

* change matrix P by taking the independant columns of A and vectors from

* the canonical basis. Update invp = 1/P, and L in {0,1}^n, with L[i] = 0

* of P[i] = e_i, and 1 if P[i] = A_i (i-th column of A)

* if P[i] = e_i, and 1 if P[i] = A_i (i-th column of A)

static GEN

relationrank(GEN *A, long r, GEN L)

{

long i, n = lg(L)-1, av = avma, lim = stack_lim(av,1);

long i, n = lg(L)-1;

gpmem_t av = avma, lim = stack_lim(av, 1);

GEN invp = idmat(n);

if (!r) return invp;

Line 2152 relationrank(GEN *A, long r, GEN L)

Line 2292 relationrank(GEN *A, long r, GEN L)

return gerepilecopy(av, invp);

}

/* SMALLBUCHINIT */

static GEN

codeprime(GEN bnf, GEN pr)

decode_pr_lists(GEN nf, GEN pfc)

{

long j,av=avma,tetpil;

long i, p, pmax, n = degpol(nf[1]), l = lg(pfc);

GEN p,al,fa,p1;

GEN t, L;

p=(GEN)pr[1]; al=(GEN)pr[2]; fa=primedec(bnf,p);

pmax = 0;

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

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

if (gegal(al,gmael(fa,j,2)))

{

t = (GEN)pfc[i]; p = itos(t) / n;

p1=mulsi(lg(al)-1,p); tetpil=avma;

if (p > pmax) pmax = p;

return gerepile(av,tetpil,addsi(j-1,p1));

}

L = cgetg(pmax+1, t_VEC);

err(talker,"bug in codeprime/smallbuchinit");

for (p=1; p<=pmax; p++) L[p] = 0;

return NULL; /* not reached */

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

{

t = (GEN)pfc[i]; p = itos(t) / n;

if (!L[p]) L[p] = (long)primedec(nf, stoi(p));

}

return L;

}

static GEN

decodeprime(GEN nf, GEN co)

decodeprime(GEN T, GEN L, long n)

{

long n,indi,av=avma;

long t = itos(T);

GEN p,rem,p1;

return gmael(L, t/n, t%n + 1);

}

static GEN

codeprime(GEN L, long N, GEN pr)

{

long p = itos((GEN)pr[1]);

return utoi( N*p + pr_index((GEN)L[p], pr)-1 );

}

n=lg(nf[7])-1; p=dvmdis(co,n,&rem); indi=itos(rem)+1;

static GEN

p1=primedec(nf,p);

codeprimes(GEN Vbase, long N)

return gerepilecopy(av,(GEN)p1[indi]);

{

GEN v, L = get_pr_lists(Vbase, N, 1);

long i, l = lg(Vbase);

v = cgetg(l, t_VEC);

for (i=1; i<l; i++) v[i] = (long)codeprime(L, N, (GEN)Vbase[i]);

return v;

}

/* v = bnf[10] */

Line 2226 makecycgen(GEN bnf)

Line 2385 makecycgen(GEN bnf)

if (y && !fact_ok(nf,y,NULL,gen,(GEN)D[i])) y = NULL;

if (y) { h[i] = (long)to_famat_all(y,gun); continue; }

}

y = isprincipalfact(bnf, gen, (GEN)D[i], NULL,

y = isprincipalfact(bnf, gen, (GEN)D[i], NULL, nf_GENMAT|nf_FORCE);

nf_GEN|nf_GENMAT|nf_FORCE|nf_GIVEPREC);

h[i] = y[2];

if (typ(y) != t_INT)

h[i] = y[2];

else

{

y = idealpow(nf, (GEN)gen[i], (GEN)cyc[i]);

h[i] = isprincipalgenforce(bnf,y)[2];

}

return h;

}

Line 2243 makecycgen(GEN bnf)

Line 2395 makecycgen(GEN bnf)

static GEN

makematal(GEN bnf)

{

GEN W,B,pFB,vp,nf,ma, WB_C;

GEN W,B,pFB,nf,ma, WB_C;

long lW,lma,j,prec;

ma = get_matal((GEN)bnf[10]);

Line 2252 makematal(GEN bnf)

Line 2404 makematal(GEN bnf)

W = (GEN)bnf[1];

B = (GEN)bnf[2];

WB_C= (GEN)bnf[4];

vp = (GEN)bnf[6];

nf = (GEN)bnf[7];

lW=lg(W)-1; lma=lW+lg(B);

pFB=cgetg(lma,t_VEC); ma=(GEN)bnf[5]; /* reorder factor base */

pFB = get_Vbase(bnf);

for (j=1; j<lma; j++) pFB[j] = ma[itos((GEN)vp[j])];

ma = cgetg(lma,t_MAT);

prec = prec_arch(bnf);

Line 2274 makematal(GEN bnf)

Line 2424 makematal(GEN bnf)

ma[j] = (long)y; continue;

}

if (!y) y = isprincipalfact(bnf,pFB,ex,C, nf_GEN|nf_FORCE|nf_GIVEPREC);

if (!y) y = isprincipalfact(bnf,pFB,ex,C, nf_GENMAT|nf_FORCE|nf_GIVEPREC);

if (typ(y) != t_INT)

{

if (DEBUGLEVEL>1) fprintferr("%ld ",j);

Line 2284 makematal(GEN bnf)

Line 2434 makematal(GEN bnf)

prec = itos(y); j--;

if (DEBUGLEVEL) err(warnprec,"makematal",prec);

nf = nfnewprec(nf,prec);

bnf = bnfinit0(nf,1,NULL,prec); setrand(c);

bnf = bnfinit0(nf,1,NULL,prec); (void)setrand(c);

}

if (DEBUGLEVEL>1) fprintferr("\n");

return ma;

Line 2316 check_and_build_cycgen(GEN bnf)

Line 2466 check_and_build_cycgen(GEN bnf)

GEN cycgen = get_cycgen((GEN)bnf[10]);

if (!cycgen)

{

long av = avma;

gpmem_t av = avma;

if (DEBUGLEVEL) err(warner,"completing bnf (building cycgen)");

bnfinsert(bnf, makecycgen(bnf), 2); avma = av;

cycgen = get_cycgen((GEN)bnf[10]);

Line 2330 check_and_build_matal(GEN bnf)

Line 2480 check_and_build_matal(GEN bnf)

GEN matal = get_matal((GEN)bnf[10]);

if (!matal)

{

long av = avma;

gpmem_t av = avma;

if (DEBUGLEVEL) err(warner,"completing bnf (building matal)");

bnfinsert(bnf, makematal(bnf), 1); avma = av;

matal = get_matal((GEN)bnf[10]);

}

return matal;

}

GEN

smallbuchinit(GEN pol,GEN gcbach,GEN gcbach2,GEN gRELSUP,GEN gborne,long nbrelpid,long minsFB,long prec)

{

long av=avma,av1,k;

GEN y, bnf, nf, res, p1;

GEN y,bnf,pFB,vp,nf,mas,res,uni,v1,v2,v3;

gpmem_t av = avma;

if (typ(pol)==t_VEC) bnf = checkbnf(pol);

else

{

bnf=buchall(pol,gcbach,gcbach2,gRELSUP,gborne,nbrelpid,minsFB,-3,prec);

const long fl = nf_INIT | nf_UNITS | nf_FORCE;

bnf = buchall(pol,gcbach,gcbach2,gRELSUP,gborne,nbrelpid,minsFB,fl,prec);

bnf = checkbnf_discard(bnf);

}

pFB=(GEN)bnf[5]; vp=(GEN)bnf[6]; nf=(GEN)bnf[7];

nf = (GEN)bnf[7];

mas=(GEN)nf[5]; res=(GEN)bnf[8]; uni=(GEN)res[5];

res = (GEN)bnf[8];

y=cgetg(13,t_VEC); y[1]=lcopy((GEN)nf[1]); y[2]=lcopy(gmael(nf,2,1));

y = cgetg(13,t_VEC);

y[3]=lcopy((GEN)nf[3]); y[4]=lcopy((GEN)nf[7]);

y[1] = nf[1];

y[5]=lcopy((GEN)nf[6]); y[6]=lcopy((GEN)mas[5]);

y[2] = mael(nf,2,1);

y[7]=lcopy((GEN)bnf[1]); y[8]=lcopy((GEN)bnf[2]);

y[3] = nf[3];

v1=cgetg(lg(pFB),t_VEC); y[9]=(long)v1;

y[4] = nf[7];

for (k=1; k<lg(pFB); k++)

y[5] = nf[6];

v1[k]=(long)codeprime(bnf,(GEN)pFB[itos((GEN)vp[k])]);

y[6] = mael(nf,5,5);

v2=cgetg(3,t_VEC); y[10]=(long)v2;

y[7] = bnf[1];

v2[1]=lcopy(gmael(res,4,1));

y[8] = bnf[2];

v2[2]=(long)algtobasis(bnf,gmael(res,4,2));

y[9] = (long)codeprimes((GEN)bnf[5], degpol(nf[1]));

v3=cgetg(lg(uni),t_VEC); y[11]=(long)v3;

for (k=1; k<lg(uni); k++)

p1 = cgetg(3, t_VEC);

v3[k] = (long)algtobasis(bnf,(GEN)uni[k]);

p1[1] = mael(res,4,1);

av1 = avma;

p1[2] = (long)algtobasis(bnf,gmael(res,4,2));

y[12] = gcmp0((GEN)bnf[10])? lpilecopy(av1, makematal(bnf))

y[10] = (long)p1;

: lcopy((GEN)bnf[10]);

return gerepileupto(av,y);

y[11] = (long)algtobasis(bnf, (GEN)res[5]);

y[12] = gcmp0((GEN)bnf[10])? (long)makematal(bnf): bnf[10];

return gerepilecopy(av, y);

}

static GEN

get_regulator(GEN mun,long prec)

get_regulator(GEN mun)

{

long av,tetpil;

gpmem_t av = avma;

GEN p1;

GEN A;

if (lg(mun)==1) return gun;

av=avma; p1 = gtrans(greal(mun));

A = gtrans( greal(mun) );

setlg(p1,lg(p1)-1); p1 = det(p1);

setlg(A, lg(A)-1);

tetpil=avma; return gerepile(av,tetpil,gabs(p1,prec));

return gerepileupto(av, gabs(det(A), 0));

}

/* return corrected archimedian components for elts of x (vector)

* (= log(sigma_i(x)) - log(|Nx|) / [K:Q]) */

static GEN

log_poleval(GEN P, GEN *ro, long i, GEN nf, long prec0)

get_archclean(GEN nf, GEN x, long prec, int units)

{

GEN x = poleval(P, (GEN)(*ro)[i]);

long k,N, la = lg(x);

long prec = prec0, k = 0;

GEN M = cgetg(la,t_MAT);

while (gcmp0(x) || ((k = gprecision(x)) && k < DEFAULTPREC))

{

prec = (prec-2)<<1;

if (DEBUGLEVEL) err(warnprec,"log_poleval",prec);

*ro = get_roots((GEN)nf[1], itos(gmael(nf,2,1)),lg(*ro)-1,prec);

x = poleval(P, (GEN)(*ro)[i]);

}

if (k > prec0) x = gprec_w(x,prec0);

return glog(x, prec0);

}

/* return corrected archimedian components

if (la == 1) return M;

* (= log(sigma_i(x)) - log(|Nx|)/[K:Q]) for a (matrix of elements) */

N = degpol(nf[1]);

static GEN

get_arch2_i(GEN nf, GEN a, long prec, int units)

{

GEN M,N, ro = dummycopy((GEN)nf[6]);

long j,k, la = lg(a), ru = lg(ro), r1 = nf_get_r1(nf);

M = cgetg(la,t_MAT); if (la == 1) return M;

if (typ(a[1]) == t_COL) a = gmul((GEN)nf[7], a);

if (units) N = NULL; /* no correction necessary */

else

{

GEN pol = (GEN)nf[1];

N = cgetg(la,t_VEC);

for (k=1; k<la; k++) N[k] = (long)gabs(subres(pol,(GEN)a[k]),0);

N = gdivgs(glog(N,prec), - (degpol(pol))); /* - log(|norm|) / [K:Q] */

}

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

{

GEN p, c = cgetg(ru,t_COL); M[k] = (long)c;

GEN c = get_arch(nf, (GEN)x[k], prec);

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

if (!units) c = cleanarch(c, N, prec);

{

M[k] = (long)c;

p = log_poleval((GEN)a[k],&ro,j,nf,prec);

if (N) p = gadd(p, (GEN)N[k]);

c[j]=(j<=r1)? (long) p: lmul2n(p,1);

}

return M;

}

static GEN

get_arch2(GEN nf, GEN a, long prec, int units)

{

long av = avma;

return gerepilecopy(av, get_arch2_i(nf,a,prec,units));

}

static void

my_class_group_gen(GEN bnf, GEN *ptcl, GEN *ptcl2, long prec)

my_class_group_gen(GEN bnf, long prec, GEN nf0, GEN *ptcl, GEN *ptcl2)

{

GEN W=(GEN)bnf[1], C=(GEN)bnf[4], vperm=(GEN)bnf[6], nf=(GEN)bnf[7], *gptr[2];

GEN W=(GEN)bnf[1], C=(GEN)bnf[4], nf=(GEN)bnf[7];

long av = avma;

class_group_gen(nf,W,C,get_Vbase(bnf),prec,nf0, ptcl,ptcl2);

vectbase = (GEN)bnf[5]; /* needed by class_group_gen */

class_group_gen(nf,W,C,vperm,ptcl,ptcl2, prec);

gptr[0]=ptcl; gptr[1]=ptcl2; gerepilemany(av,gptr,2);

}

GEN

bnfnewprec(GEN bnf, long prec)

{

GEN nf,ro,res,p1,funits,mun,matal,clgp,clgp2,y;

GEN nf0 = (GEN)bnf[7], nf, res, funits, mun, matal, clgp, clgp2, y;

long r1,r2,ru,pl1,pl2,prec1;

gpmem_t av = avma;

long r1, r2, prec1;

bnf = checkbnf(bnf);

if (prec <= 0) return nfnewprec(checknf(bnf),prec);

y = cgetg(11,t_VEC);

funits = check_units(bnf,"bnfnewprec");

nf = (GEN)bnf[7];

ro = (GEN)nf[6]; ru = lg(ro)-1;

nf_get_sign(nf, &r1, &r2);

r1 = nf_get_r1(nf);

funits = algtobasis(nf, check_units(bnf,"bnfnewprec"));

r2 = (r1 + ru) >> 1;

pl1 = (ru == 1 && r1 == 0)? 0: gexpo(funits);

pl2 = gexpo(ro);

prec1 = prec;

prec += ((ru + r2 - 1) * (pl1 + (ru + r2) * pl2)) >> TWOPOTBITS_IN_LONG;

if (r2 > 1 || r1 != 0)

nf = nfnewprec((GEN)bnf[7],prec);

prec += 1 + (gexpo(funits) >> TWOPOTBITS_IN_LONG);

res = cgetg(7,t_VEC);

nf = nfnewprec(nf0,prec);

ro = (GEN)nf[6];

mun = get_archclean(nf,funits,prec,1);

mun = get_arch2(nf,funits,prec,1);

if (prec != prec1) { mun = gprec_w(mun,prec1); prec = prec1; }

res[2]=(long)get_regulator(mun,prec);

p1 = (GEN)bnf[8];

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

res[4]=lcopy((GEN)p1[4]);

res[5]=lcopy((GEN)p1[5]);

res[6]=lcopy((GEN)p1[6]);

y[1]=lcopy((GEN)bnf[1]);

y[2]=lcopy((GEN)bnf[2]);

y[3]=(long)mun;

matal = check_and_build_matal(bnf);

y[4]=(long)get_arch2(nf,matal,prec,0);

y[5]=lcopy((GEN)bnf[5]);

y = dummycopy(bnf);

y[6]=lcopy((GEN)bnf[6]);

y[3] = (long)mun;

y[7]=(long)nf;

y[4] = (long)get_archclean(nf,matal,prec,0);

y[8]=(long)res;

y[7] = (long)nf;

my_class_group_gen(y,&clgp,&clgp2,prec);

my_class_group_gen(y,prec,nf0, &clgp,&clgp2);

res[1]=(long)clgp;

res = dummycopy((GEN)bnf[8]);

y[9]=(long)clgp2;

res[1] = (long)clgp;

y[10]=lcopy((GEN)bnf[10]); return y;

res[2] = (long)get_regulator(mun);

y[8] = (long)res;

y[9] = (long)clgp2; return gerepilecopy(av, y);

}

GEN

Line 2505 bnrnewprec(GEN bnr, long prec)

Line 2607 bnrnewprec(GEN bnr, long prec)

return y;

}

static void

nfbasic_from_sbnf(GEN sbnf, nfbasic_t *T)

{

T->x = (GEN)sbnf[1];

T->dK = (GEN)sbnf[3];

T->bas = (GEN)sbnf[4];

T->index= get_nfindex(T->bas);

T->r1 = itos((GEN)sbnf[2]);

T->dx = NULL;

T->lead = NULL;

T->basden = NULL;

}

static GEN

get_clfu(GEN clgp, GEN reg, GEN c1, GEN zu, GEN fu, long k)

{

GEN z = cgetg(7, t_VEC);

z[1]=(long)clgp; z[2]=(long)reg; z[3]=(long)c1;

z[4]=(long)zu; z[5]=(long)fu; z[6]=lstoi(k); return z;

}

GEN

bnfmake(GEN sbnf, long prec)

{

long av = avma, j,k,n,r1,r2,ru,lpf;

long j, k, l, n;

GEN p1,x,bas,ro,nf,mun,funits,index;

gpmem_t av = avma;

GEN pfc,vp,mc,clgp,clgp2,res,y,W,racu,reg,matal;

GEN p1, bas, ro, nf, mun, fu, L;

GEN pfc, mc, clgp, clgp2, res, y, W, zu, reg, matal, Vbase;

nfbasic_t T;

if (typ(sbnf)!=t_VEC || lg(sbnf)!=13)

if (typ(sbnf) != t_VEC || lg(sbnf) != 13) err(typeer,"bnfmake");

err(talker,"incorrect sbnf in bnfmake");

if (prec < DEFAULTPREC) prec = DEFAULTPREC;

x=(GEN)sbnf[1]; bas=(GEN)sbnf[4]; n=lg(bas)-1;

r1=itos((GEN)sbnf[2]); r2=(n-r1)>>1; ru=r1+r2;

ro=(GEN)sbnf[5];

if (prec > gprecision(ro)) ro=get_roots(x,r1,ru,prec);

index = gun;

for (j=2; j<=n; j++) index = mulii(index, denom(leading_term(bas[j])));

nf=cgetg(10,t_VEC);

nfbasic_from_sbnf(sbnf, &T);

nf[1]=sbnf[1]; p1=cgetg(3,t_VEC); p1[1]=lstoi(r1); p1[2]=lstoi(r2);

ro = (GEN)sbnf[5];

nf[2]=(long)p1;

if (prec > gprecision(ro)) ro = get_roots(T.x,T.r1,prec);

nf[3]=sbnf[3];

nf = nfbasic_to_nf(&T, ro, prec);

nf[4]=(long)index;

bas = (GEN)nf[7];

nf[6]=(long)ro;

nf[7]=(long)bas;

get_nf_matrices(nf, 0);

funits=cgetg(ru,t_VEC); p1 = (GEN)sbnf[11];

p1 = (GEN)sbnf[11]; l = lg(p1); fu = cgetg(l, t_VEC);

for (k=1; k < lg(p1); k++)

for (k=1; k < l; k++) fu[k] = lmul(bas, (GEN)p1[k]);

funits[k] = lmul(bas,(GEN)p1[k]);

mun = get_archclean(nf,fu,prec,1);

mun = get_arch2_i(nf,funits,prec,1);

prec=gprecision(ro); if (prec<DEFAULTPREC) prec=DEFAULTPREC;

prec = gprecision(ro);

matal = get_matal((GEN)sbnf[12]);

if (!matal) matal = (GEN)sbnf[12];

mc = get_arch2_i(nf,matal,prec,0);

mc = get_archclean(nf,matal,prec,0);

pfc=(GEN)sbnf[9]; lpf=lg(pfc);

pfc = (GEN)sbnf[9];

vectbase=cgetg(lpf,t_COL); vp=cgetg(lpf,t_COL);

l = lg(pfc);

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

Vbase = cgetg(l,t_COL);

{

L = decode_pr_lists(nf, pfc);

vp[j]=lstoi(j);

n = degpol(nf[1]);

vectbase[j]=(long)decodeprime(nf,(GEN)pfc[j]);

for (j=1; j<l; j++) Vbase[j] = (long)decodeprime((GEN)pfc[j], L, n);

}

W = (GEN)sbnf[7];

class_group_gen(nf,W,mc,vp,&clgp,&clgp2, prec); /* uses vectbase */

class_group_gen(nf,W,mc,Vbase,prec,NULL, &clgp,&clgp2);

reg = get_regulator(mun,prec);

reg = get_regulator(mun);

p1=cgetg(3,t_VEC); racu=(GEN)sbnf[10];

p1 = cgetg(3,t_VEC); zu = (GEN)sbnf[10];

p1[1]=racu[1]; p1[2]=lmul(bas,(GEN)racu[2]);

p1[1] = zu[1];

racu=p1;

p1[2] = lmul(bas,(GEN)zu[2]); zu = p1;

res=cgetg(7,t_VEC);

res = get_clfu(clgp,reg,dbltor(1.0),zu,fu,1000);

res[1]=(long)clgp; res[2]=(long)reg; res[3]=(long)dbltor(1.0);

res[4]=(long)racu; res[5]=(long)funits; res[6]=lstoi(1000);

y=cgetg(11,t_VEC);

y[1]=(long)W; y[2]=sbnf[8]; y[3]=(long)mun;

y[4]=(long)mc; y[5]=(long)vectbase; y[6]=(long)vp;

y[4]=(long)mc; y[5]=(long)Vbase; y[6]=zero;

y[7]=(long)nf; y[8]=(long)res; y[9]=(long)clgp2; y[10]=sbnf[12];

y[7]=(long)nf; y[8]=(long)res; y[9]=(long)clgp2;

return gerepilecopy(av,y);

y[10] = sbnf[12]; return gerepilecopy(av,y);

}

static GEN

classgroupall(GEN P, GEN data, long flag, long prec)

{

long court[3],doubl[4];

long av=avma,flun,lx, minsFB=3,nbrelpid=4;

long fl, lx, minsFB=3, nbrelpid=4;

gpmem_t av=avma;

GEN bach=doubl,bach2=doubl,RELSUP=court,borne=gun;

if (!data) lx=1;

Line 2594 classgroupall(GEN P, GEN data, long flag, long prec)

Line 2706 classgroupall(GEN P, GEN data, long flag, long prec)

}

switch(flag)

{

case 0: flun=-2; break;

case 0: fl = nf_INIT | nf_UNITS; break;

case 1: flun=-3; break;

case 1: fl = nf_INIT | nf_UNITS | nf_FORCE; break;

case 2: flun=-1; break;

case 2: fl = nf_INIT | nf_ROOT1; break;

case 3: return smallbuchinit(P,bach,bach2,RELSUP,borne,nbrelpid,minsFB,prec);

case 4: flun=2; break;

case 4: fl = nf_UNITS; break;

case 5: flun=3; break;

case 5: fl = nf_UNITS | nf_FORCE; break;

case 6: flun=0; break;

case 6: fl = 0; break;

default: err(flagerr,"classgroupall");

return NULL; /* not reached */

}

return buchall(P,bach,bach2,RELSUP,borne,nbrelpid,minsFB,flun,prec);

return buchall(P,bach,bach2,RELSUP,borne,nbrelpid,minsFB,fl,prec);

}

GEN

Line 2633 bnfinit0(GEN P, long flag, GEN data, long prec)

Line 2745 bnfinit0(GEN P, long flag, GEN data, long prec)

GEN

classgrouponly(GEN P, GEN data, long prec)

{

ulong av = avma;

gpmem_t av = avma;

GEN z;

if (typ(P)==t_INT)

Line 2648 classgrouponly(GEN P, GEN data, long prec)

Line 2760 classgrouponly(GEN P, GEN data, long prec)

GEN

regulator(GEN P, GEN data, long prec)

{

ulong av = avma;

gpmem_t av = avma;

GEN z;

if (typ(P)==t_INT)

Line 2673 col_0(long n)

Line 2785 col_0(long n)

GEN c = (GEN) gpmalloc(sizeof(long)*(n+1));

long i;

for (i=1; i<=n; i++) c[i]=0;

c[0] = evaltyp(t_VECSMALL) | evallg(n);

c[0] = evaltyp(t_VECSMALL) | evallg(n+1);

return c;

}

Line 2687 cgetc_col(long n, long prec)

Line 2799 cgetc_col(long n, long prec)

}

static GEN

buchall_end(GEN nf,GEN CHANGE,long fl,long k, GEN fu, GEN clg1, GEN clg2,

buchall_end(GEN nf,GEN CHANGE,long fl,GEN res, GEN clg2, GEN W, GEN B,

GEN reg, GEN c_1, GEN zu, GEN W, GEN B,

GEN A, GEN matarch, GEN Vbase)

GEN xarch, GEN matarch, GEN vectbase, GEN vperm)

{

long i, l = labs(fl)>1? 11: fl? 9: 8;

long l = (fl & nf_UNITS)? 7

GEN p1,z, RES = cgetg(11,t_COL);

: (fl & nf_ROOT1)? 5: 4;

GEN p1, z;

setlg(RES,l);

setlg(res, l);

RES[5]=(long)clg1;

if (! (fl & nf_INIT))

RES[6]=(long)reg;

RES[7]=(long)c_1;

RES[8]=(long)zu;

RES[9]=(long)fu;

RES[10]=lstoi(k);

if (fl>=0)

{

RES[1]=nf[1];

GEN x = cgetg(5, t_VEC);

RES[2]=nf[2]; p1=cgetg(3,t_VEC); p1[1]=nf[3]; p1[2]=nf[4];

x[1]=nf[1];

RES[3]=(long)p1;

x[2]=nf[2]; p1=cgetg(3,t_VEC); p1[1]=nf[3]; p1[2]=nf[4];

RES[4]=nf[7];

x[3]=(long)p1;

z=cgetg(2,t_MAT); z[1]=lcopy(RES); return z;

x[4]=nf[7];

z=cgetg(2,t_MAT); z[1]=(long)concatsp(x, res); return z;

}

z=cgetg(11,t_VEC);

z[1]=(long)W;

z[2]=(long)B;

z[3]=(long)xarch;

z[3]=(long)A;

z[4]=(long)matarch;

z[5]=(long)vectbase;

z[5]=(long)Vbase;

for (i=lg(vperm)-1; i>0; i--) vperm[i]=lstoi(vperm[i]);

z[6]=zero;

settyp(vperm, t_VEC);

z[7]=(long)nf;

z[6]=(long)vperm;

z[8]=(long)res;

z[7]=(long)nf; RES+=4; RES[0]=evaltyp(t_VEC) | evallg(l-4);

z[8]=(long)RES;

z[9]=(long)clg2;

z[10]=zero; /* dummy: we MUST have lg(bnf) != lg(nf) */

if (CHANGE) { p1=cgetg(3,t_VEC); p1[1]=(long)z; p1[2]=(long)CHANGE; z=p1; }

return gcopy(z);

return z;

}

static GEN

buchall_for_degree_one_pol(GEN nf, GEN CHANGE, long flun)

{

long av = avma, k = EXP220;

gpmem_t av = avma;

GEN W,B,xarch,matarch,vectbase,vperm;

GEN W,B,A,matarch,Vbase,res;

GEN fu=cgetg(1,t_VEC), reg=gun, c_1=gun, zu=cgetg(3,t_VEC);

GEN fu=cgetg(1,t_VEC), R=gun, c1=gun, zu=cgetg(3,t_VEC);

GEN clg1=cgetg(4,t_VEC), clg2=cgetg(4,t_VEC);

clg1[1]=un; clg1[2]=clg1[3]=clg2[3]=lgetg(1,t_VEC);

clg1[1]=un; clg1[2]=clg1[3]=clg2[2]=clg2[3]=lgetg(1,t_VEC);

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

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

zu[1]=deux; zu[2]=lnegi(gun);

W=B=xarch=matarch=cgetg(1,t_MAT);

W=B=A=matarch=cgetg(1,t_MAT);

vectbase=cgetg(1,t_COL); vperm=cgetg(1,t_VEC);

Vbase=cgetg(1,t_COL);

return gerepileupto(av, buchall_end(nf,CHANGE,flun,k,fu,clg1,clg2,reg,c_1,zu,W,B,xarch,matarch,vectbase,vperm));

res = get_clfu(clg1, R, c1, zu, fu, EXP220);

return gerepilecopy(av, buchall_end(nf,CHANGE,flun,res,clg2,W,B,A,matarch,Vbase));

}

static GEN

/* return (small set of) indices of columns generating the same lattice as x.

get_LLLnf(GEN nf, long prec)

* Assume HNF(x) is inexpensive (few rows, many columns).

* Dichotomy approach since interesting columns may be at the very end */

GEN

extract_full_lattice(GEN x)

{

GEN M = gmael(nf,5,1);

long dj, j, k, l = lg(x);

GEN T2 = gmael(nf,5,3);

GEN h, h2, H, v;

GEN cbase = lllgramintern(T2,100,1,prec);

GEN v = cgetg(5,t_VEC);

if (l < 200) return NULL; /* not worth it */

if (!cbase) return NULL;

if (gegal(cbase, idmat(lg(T2)-1))) cbase = NULL;

v = cgetg(l, t_VECSMALL);

v[1] = (long) (cbase? qf_base_change(T2, cbase, 1): T2);

setlg(v, 1);

v[2] = (long) (cbase? gmul(M, cbase): M);

H = hnfall_i(x, NULL, 1);

v[3] = (long) cbase;

h = cgetg(1, t_MAT);

v[4] = (long) (cbase? ZM_inv(cbase,gun): NULL); return v;

dj = 1;

for (j = 1; j < l; )

{

gpmem_t av = avma;

long lv = lg(v);

for (k = 0; k < dj; k++) v[lv+k] = j+k;

setlg(v, lv + dj);

h2 = hnfall_i(vecextract_p(x, v), NULL, 1);

if (gegal(h, h2))

{ /* these dj columns can be eliminated */

avma = av; setlg(v, lv);

j += dj;

if (j >= l) break; /* shouldn't occur */

dj <<= 1;

if (j + dj >= l) dj = (l - j) >> 1;

}

else if (dj > 1)

{ /* at least one interesting column, try with first half of this set */

avma = av; setlg(v, lv);

dj >>= 1;

}

else

{ /* this column should be kept */

if (gegal(h2, H)) break;

h = h2; j++;

}

return v;

}

GEN

buchall(GEN P,GEN gcbach,GEN gcbach2,GEN gRELSUP,GEN gborne,long nbrelpid,

long minsFB,long flun,long prec)

{

ulong av = avma,av0,av1,limpile;

gpmem_t av = avma, av0, av1, limpile;

long N,R1,R2,RU,PRECREG,PRECLLL,KCCO,RELSUP,LIMC,LIMC2,lim;

long N,R1,R2,RU,PRECREG,PRECLLL,PRECLLLadd,KCCO,RELSUP,LIMC,LIMC2,lim;

long nlze,sreg,nrelsup,nreldep,phase,matmax,i,j,k,ss,cglob;

long nlze,zc,nrelsup,nreldep,phase,matmax,i,j,k,seed,MAXRELSUP;

long first=1, sfb_increase=0, sfb_trials=0, precdouble=0, precadd=0;

long sfb_increase, sfb_trials, precdouble=0, precadd=0;

double cbach,cbach2,drc,LOGD2;

long cglob; /* # of relations found so far */

GEN p1,vecT2,fu,zu,nf,LLLnf,D,xarch,W,reg,lfun,z,clh,vperm,subFB;

double cbach, cbach2, drc, LOGD2;

GEN B,C,c1,sublambda,pdep,parch,liste,invp,clg1,clg2, *mat;

GEN vecG,fu,zu,nf,D,A,W,R,Res,z,h,list_jideal;

GEN CHANGE=NULL, extramat=NULL, extraC=NULL, list_jideal=NULL;

GEN res,L,resc,B,C,lambda,pdep,liste,invp,clg1,clg2,Vbase;

GEN *mat; /* raw relations found (as VECSMALL, not HNF-reduced) */

GEN first_nz; /* first_nz[i] = index of 1st non-0 entry in mat[i] */

GEN CHANGE=NULL, extramat=NULL, extraC=NULL;

char *precpb = NULL;

FB_t F;

if (DEBUGLEVEL) timer2();

if (DEBUGLEVEL) (void)timer2();

if (typ(P)==t_POL) nf = NULL; else { nf = checknf(P); P = (GEN)nf[1]; }

P = get_nfpol(P, &nf);

if (typ(gRELSUP) != t_INT) gRELSUP = gtrunc(gRELSUP);

RELSUP = itos(gRELSUP);

if (RELSUP<=0) err(talker,"not enough relations in bnfxxx");

/* Initializations */

fu = NULL; /* gcc -Wall */

N = degpol(P); PRECREG = max(BIGDEFAULTPREC,prec);

N = degpol(P);

PRECREG = max(BIGDEFAULTPREC,prec);

PRECLLLadd = DEFAULTPREC;

if (!nf)

{

nf = initalgall0(P, nf_REGULAR, PRECREG);

nf = initalg(P, PRECREG);

/* P was non-monic and nfinit changed it ? */

if (lg(nf)==3) { CHANGE = (GEN)nf[2]; nf = (GEN)nf[1]; }

}

if (N<=1) return buchall_for_degree_one_pol(nf,CHANGE,flun);

if (N <= 1) return buchall_for_degree_one_pol(nf,CHANGE,flun);

zu = rootsof1(nf);

zu[2] = lmul((GEN)nf[7],(GEN)zu[2]);

if (DEBUGLEVEL) msgtimer("initalg & rootsof1");

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

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

D = (GEN)nf[3]; drc = fabs(gtodouble(D));

LOGD2 = log(drc); LOGD2 = LOGD2*LOGD2;

lim = (long) (exp(-(double)N) * sqrt(2*PI*N*drc) * pow(4/PI,(double)R2));

if (lim < 3) lim = 3;

cbach = min(12., gtodouble(gcbach)); cbach /= 2;

if (cbach == 0.) err(talker,"Bach constant = 0 in bnfxxx");

if (nbrelpid <= 0) gborne = gzero;

cbach2 = gtodouble(gcbach2);

/* resc ~ sqrt(D) w / 2^r1 (2pi)^r2 = hR / Res(zeta_K, s=1) */

resc = gdiv(mulri(gsqrt(absi(D),DEFAULTPREC), (GEN)zu[1]),

gmul2n(gpowgs(Pi2n(1,DEFAULTPREC), R2), R1));

if (DEBUGLEVEL)

{

fprintferr("N = %ld, R1 = %ld, R2 = %ld\nD = %Z\n",N,R1,R2,D);

fprintferr("N = %ld, R1 = %ld, R2 = %ld, RU = %ld\n",N,R1,R2,RU);

av0 = avma; mat = NULL; first_nz = NULL;

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

}

av0 = avma; mat = NULL; FB = NULL;

START:

avma = av0;

seed = getrand();

if (first) first = 0; else desallocate(&mat);

avma = av0; desallocate(&mat, &first_nz);

if (precpb)

{

precdouble++;

Line 2829 START:

Line 2975 START:

cbach = check_bach(cbach,12.);

precadd = 0;

/* T2-LLL-reduce nf.zk */

LIMC = (long)(cbach*LOGD2);

LLLnf = get_LLLnf(nf, PRECREG);

if (LIMC < 20) { LIMC = 20; cbach = LIMC / LOGD2; }

if (!LLLnf) { precpb = "LLLnf"; goto START; }

LIMC = (long)(cbach*LOGD2); if (LIMC < 20) LIMC = 20;

LIMC2= max(3 * N, (long)(cbach2*LOGD2));

if (LIMC2 < LIMC) LIMC2 = LIMC;

if (DEBUGLEVEL) { fprintferr("LIMC = %ld, LIMC2 = %ld\n",LIMC,LIMC2); }

/* initialize FB, [sub]vperm */

Res = FBgen(&F, nf, LIMC2, LIMC);

lfun = FBgen(nf,LIMC2,LIMC);

if (!Res || !subFBgen(&F, nf, min(lim,LIMC2) + 0.5, minsFB)) goto START;

if (!lfun) goto START;

vperm = cgetg(lg(vectbase), t_VECSMALL);

if (DEBUGLEVEL)

{

if (DEBUGLEVEL>3)

{

fprintferr("\n***** IDEALS IN FACTORBASE *****\n\n");

for (i=1; i<=F.KC; i++) fprintferr("no %ld = %Z\n",i,F.LP[i]);

fprintferr("\n***** IDEALS IN SUB FACTORBASE *****\n\n");

outerr(vecextract_p(F.LP,F.subFB));

fprintferr("\n***** INITIAL PERMUTATION *****\n\n");

fprintferr("perm = %Z\n\n",F.perm);

}

msgtimer("sub factorbase (%ld elements)",lg(F.subFB)-1);

}

sfb_trials = sfb_increase = 0;

subFB = subFBgen(N,min(lim,LIMC2), minsFB, vperm, &ss);

if (!subFB) goto START;

nreldep = nrelsup = 0;

/* PRECLLL = prec for LLL-reductions (idealred)

* PRECREG = prec for archimedean components */

PRECLLL = DEFAULTPREC

PRECLLL = PRECLLLadd

+ ((expi(D)*(lg(subFB)-2)+((N*N)>>2))>>TWOPOTBITS_IN_LONG);

+ ((expi(D)*(lg(F.subFB)-2) + ((N*N)>>2)) >> TWOPOTBITS_IN_LONG);

if (!precdouble) PRECREG = prec+1;

if (PRECREG < PRECLLL) PRECREG = PRECLLL;

if (PRECREG < PRECLLL)

KCCO = KC+RU-1 + max(ss,RELSUP); /* expected number of needed relations */

{ /* very rare */

PRECREG = PRECLLL;

nf = nfnewprec(nf,PRECREG); av0 = avma;

}

KCCO = F.KC+RU-1 + RELSUP; /* expected # of needed relations */

if (DEBUGLEVEL)

fprintferr("relsup = %ld, ss = %ld, KCZ = %ld, KC = %ld, KCCO = %ld\n",

fprintferr("relsup = %ld, KCZ = %ld, KC = %ld, KCCO = %ld\n",

RELSUP,ss,KCZ,KC,KCCO);

RELSUP, F.KCZ, F.KC, KCCO);

MAXRELSUP = min(50, 4*F.KC);

matmax = 10*KCCO + MAXRELSUP; /* make room for lots of relations */

mat = (GEN*)gpmalloc(sizeof(GEN)*(matmax + 1));

reallocate(matmax, (GEN*)&mat, &first_nz);

setlg(mat, KCCO+1);

C = cgetg(KCCO+1,t_MAT);

/* trivial relations (p) = prod P^e */

cglob = 0; z = zerocol(RU);

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

for (i=1; i<=F.KCZ; i++)

{

GEN P = idealbase[i];

long p = F.FB[i];

if (isclone(P))

GEN c, P = F.LV[p];

{ /* all prime divisors in FB */

if (!isclone(P)) continue;

unsetisclone(P); cglob++;

C[cglob] = (long)z; /* 0 */

/* all prime divisors in FB */

mat[cglob] = p1 = col_0(KC);

cglob++;

k = numideal[FB[i]];

C[cglob] = (long)z; /* 0 */

p1[0] = k+1; /* for already_found_relation */

mat[cglob] = c = col_0(F.KC);

p1 += k;

k = F.iLP[p];

for (j=lg(P)-1; j; j--) p1[j] = itos(gmael(P,j,3));

first_nz[cglob] = k+1;

}

c += k;

for (j=lg(P)-1; j; j--) c[j] = itos(gmael(P,j,3));

}

if (DEBUGLEVEL) fprintferr("After trivial relations, cglob = %ld\n",cglob);

/* initialize for other relations */

for (i=cglob+1; i<=KCCO; i++)

{

mat[i] = col_0(KC);

mat[i] = col_0(F.KC);

C[i] = (long)cgetc_col(RU,PRECREG);

}

av1 = avma; liste = cgetg(KC+1,t_VECSMALL);

av1 = avma; liste = vecsmall_const(F.KC, 0);

for (i=1; i<=KC; i++) liste[i] = 0;

invp = relationrank(mat,cglob,liste);

/* relations through elements of small norm */

cglob = small_norm_for_buchall(cglob,mat,C,(long)LIMC,PRECREG,

if (gsigne(gborne) > 0)

nf,gborne,nbrelpid,invp,liste,LLLnf);

cglob = small_norm_for_buchall(cglob,mat,first_nz,C,(long)LIMC,PRECREG,&F,

nf,nbrelpid,invp,liste);

if (cglob < 0) { precpb = "small_norm"; goto START; }

avma = av1; limpile = stack_lim(av1,1);

phase = 0;

nlze = matmax = 0; /* for lint */

nlze = 0; /* for lint */

vecT2 = NULL;

vecG = NULL;

list_jideal = NULL;

/* random relations */

if (cglob == KCCO) /* enough rels, but init random_relations just in case */

Line 2904 START:

Line 3064 START:

else

{

GEN matarch;

long ss;

if (DEBUGLEVEL) fprintferr("\n#### Looking for random relations\n");

Line 2911 MORE:

Line 3072 MORE:

{

if (DEBUGLEVEL) fprintferr("*** Increasing sub factor base\n");

sfb_increase = 0;

if (++sfb_trials >= SFB_MAX) goto START;

if (++sfb_trials > SFB_MAX || !subFBgen_increase(&F, nf, SFB_STEP))

subFB = subFBgen(N,min(lim,LIMC2), lg(subFB)-1+SFB_STEP,NULL,&ss);

goto START;

if (!subFB) goto START;

nreldep = nrelsup = 0;

}

if (phase == 0) matarch = C;

else

{ /* reduced the relation matrix at least once */

long lgex = max(nlze, MIN_EXTRA); /* number of new relations sought */

long lgex = max(nlze, MIN_EXTRA); /* # of new relations sought */

long slim; /* total number of relations (including lgex new ones) */

long slim; /* total # of relations (including lgex new ones) */

setlg(extraC, lgex+1);

setlg(extramat,lgex+1); /* were allocated after hnfspec */

slim = cglob+lgex;

if (slim > matmax)

{

mat = (long**)gprealloc(mat, (2*slim+1) * sizeof(long*),

(matmax+1) * sizeof(long*));

matmax = 2 * slim;

reallocate(matmax, (GEN*)&mat, &first_nz);

}

setlg(mat, slim+1);

if (DEBUGLEVEL)

fprintferr("\n(need %ld more relation%s)\n", lgex, (lgex==1)?"":"s");

for (j=cglob+1; j<=slim; j++) mat[j] = col_0(KC);

for (j=cglob+1; j<=slim; j++) mat[j] = col_0(F.KC);

matarch = extraC - cglob; /* start at 0, the others at cglob */

}

if (!vecT2)

if (!vecG)

{

vecT2 = compute_vecT2(RU,nf);

vecG = compute_vecG(nf,PRECLLL);

av1 = avma;

}

if (!powsubFB)

if (!F.powsubFB)

{

powsubFBgen(nf,subFB,CBUCHG+1,PRECREG);

powsubFBgen(&F, nf, CBUCHG+1, PRECREG);

av1 = avma;

}

ss = random_relation(phase,cglob,(long)LIMC,PRECLLL,PRECREG,

ss = random_relation(phase,cglob,(long)LIMC,PRECREG,MAXRELSUP,

nf,subFB,vecT2,mat,matarch,list_jideal);

nf,vecG,mat,first_nz,matarch,list_jideal,&F);

if (ss < 0)

{ /* could not find relations */

if (ss != -1) precpb = "random_relation"; /* precision pb */

if (ss != -1)

{

precpb = "random_relation"; /* precision pb */

PRECLLLadd = (PRECLLLadd<<1) - 2;

}

goto START;

}

if (DEBUGLEVEL > 2) dbg_outrel(phase,cglob,vperm,mat,matarch);

if (DEBUGLEVEL > 2 && phase == 0) dbg_outrel(cglob,mat,matarch);

if (phase)

for (j=lg(extramat)-1; j>0; j--)

{

GEN c = mat[cglob+j], *g = (GEN*) extramat[j];

for (k=1; k<=KC; k++) g[k] = stoi(c[vperm[k]]);

for (k=1; k<=F.KC; k++) g[k] = stoi(c[F.perm[k]]);

}

cglob = ss;

}

Line 2968 MORE:

Line 3131 MORE:

if (phase == 0)

{ /* never reduced before */

long lgex;

W = hnfspec(mat,vperm,&pdep,&B,&C, lg(subFB)-1);

W = hnfspec(mat,F.perm,&pdep,&B,&C, lg(F.subFB)-1);

phase = 1;

nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1;

if (nlze)

if (nlze) list_jideal = vecextract_i(F.perm, 1, nlze);

{

list_jideal = cgetg(nlze+1, t_VECSMALL);

for (i=1; i<=nlze; i++) list_jideal[i] = vperm[i];

}

lgex = max(nlze, MIN_EXTRA); /* set lgex > 0, in case regulator is 0 */

/* allocate now, reduce dimension later (setlg) when lgex goes down */

extramat= cgetg(lgex+1,t_MAT);

extraC = cgetg(lgex+1,t_MAT);

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

{

extramat[j]=lgetg(KC+1,t_COL);

extramat[j]=lgetg(F.KC+1,t_COL);

extraC[j] = (long)cgetc_col(RU,PRECREG);

}

Line 2990 MORE:

Line 3149 MORE:

{

if (low_stack(limpile, stack_lim(av1,1)))

{

GEN *gptr[6];

if(DEBUGMEM>1) err(warnmem,"buchall");

gptr[0]=&W; gptr[1]=&C; gptr[2]=&B; gptr[3]=&pdep;

gerepileall(av1,6, &W,&C,&B,&pdep,&extramat,&extraC);

gptr[4]=&extramat; gptr[5]=&extraC; gerepilemany(av1,gptr,6);

}

list_jideal = NULL;

W = hnfadd(W,vperm,&pdep,&B,&C, extramat,extraC);

W = hnfadd(W,F.perm,&pdep,&B,&C, extramat,extraC);

nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1;

if (nlze && ++nreldep > MAXRELSUP) { sfb_increase=1; goto MORE; }

}

if (nlze) goto MORE; /* dependent rows */

/* first attempt at regulator for the check */

sreg = (lg(mat)-1) - (lg(B)-1) - (lg(W)-1); /* = zc (cf hnffinal) */

zc = (lg(mat)-1) - (lg(B)-1) - (lg(W)-1);

xarch = cgetg(sreg+1,t_MAT); /* cols corresponding to units */

A = vecextract_i(C, 1, zc); /* cols corresponding to units */

for (j=1; j<=sreg; j++) xarch[j] = C[j];

R = compute_multiple_of_R(A,RU,N,&lambda);

reg = compute_multiple_of_R(xarch,RU,N,&sublambda);

if (!R)

if (!reg)

{ /* not full rank for units */

if (DEBUGLEVEL) fprintferr("regulator is zero.\n");

if (++nrelsup > MAXRELSUP) goto START;

nlze = MIN_EXTRA; goto MORE;

}

/* anticipate precision problems */

if (!sublambda) { precpb = "bestappr"; goto START; }

if (!lambda) { precpb = "bestappr"; goto START; }

/* class number */

h = dethnf_i(W);

clh = dethnf_i(W);

if (DEBUGLEVEL) fprintferr("\n#### Tentative class number: %Z\n", h);

if (DEBUGLEVEL) fprintferr("\n#### Tentative class number: %Z\n", clh);

/* check */

z = mulrr(Res, resc); /* ~ hR if enough relations, a multiple otherwise */

z = mulrr(lfun,gmul(gmul2n(gpuigs(shiftr(mppi(DEFAULTPREC),1),-R2),-R1),

switch (compute_R(lambda, divir(h,z), &L, &R))

gsqrt(absi(D),DEFAULTPREC)));

z = mulri(z,(GEN)zu[1]);

/* z = Res (zeta_K, s = 1) * w D^(1/2) / [ 2^r1 (2pi)^r2 ] = h R */

p1 = gmul2n(divir(clh,z), 1);

/* c1 should be close to 2, and not much smaller */

c1 = compute_check(sublambda,p1,&parch,&reg);

/* precision problems? */

if (!c1 || gcmpgs(gmul2n(c1,1),3) < 0)

{ /* has to be a precision problem unless we cheat on Bach constant */

if (!precdouble) precpb = "compute_check";

goto START;

}

if (gcmpgs(c1,3) > 0)

{

if (++nrelsup <= MAXRELSUP)

case PRECI: /* precision problem unless we cheat on Bach constant */

{

if (!precdouble) precpb = "compute_R";

if (DEBUGLEVEL) fprintferr("\n ***** check = %f\n",gtodouble(c1)/2);

goto START;

nlze = MIN_EXTRA; goto MORE;

case RELAT: /* not enough relations */

}

if (++nrelsup <= MAXRELSUP) nlze = MIN_EXTRA; else sfb_increase = 1;

if (cbach < 11.99) { sfb_increase = 1; goto MORE; }

goto MORE;

err(warner,"suspicious check. Try to increase extra relations");

}

/* DONE */

if (KCZ2 > KCZ)

if (!be_honest(&F, nf, PRECLLL)) goto START;

{ /* "be honest" */

if (!powsubFB) powsubFBgen(nf,subFB,CBUCHG+1,PRECREG);

if (!be_honest(nf,subFB,PRECLLL)) goto START;

}

/* fundamental units */

if (flun < 0 || flun > 1)

if (flun & nf_INIT || flun & nf_UNITS)

{

xarch = cleancol(gmul(xarch,parch),N,PRECREG);

GEN v = extract_full_lattice(L); /* L may be very large */

if (DEBUGLEVEL) msgtimer("cleancol");

if (v)

{

A = vecextract_p(A, v);

L = vecextract_p(L, v);

}

/* arch. components of fund. units */

A = cleanarch(gmul(A,lllint(L)), N, PRECREG);

if (DEBUGLEVEL) msgtimer("cleanarch");

}

if (labs(flun) > 1)

if (flun & nf_UNITS)

{

fu = getfu(nf,&xarch,reg,flun,&k,PRECREG);

fu = getfu(nf,&A,R,flun,&k,PRECREG);

if (k <= 0 && labs(flun) > 2)

if (k <= 0 && flun & nf_FORCE)

{

if (k < 0)

if (k < 0) precadd = (DEFAULTPREC-2) + ((-k) >> TWOPOTBITS_IN_LONG);

{

(void)setrand(seed);

precadd = (DEFAULTPREC-2) + ((-k) >> TWOPOTBITS_IN_LONG);

if (precadd <= 0) precadd = (DEFAULTPREC-2);

}

precpb = "getfu"; goto START;

}

desallocate(&mat, &first_nz);

/* class group generators */

i = lg(C)-sreg; C += sreg; C[0] = evaltyp(t_MAT)|evallg(i);

i = lg(C)-zc; C += zc; C[0] = evaltyp(t_MAT)|evallg(i);

C = cleancol(C,N,PRECREG);

C = cleanarch(C,N,PRECREG);

class_group_gen(nf,W,C,vperm, &clg1, &clg2, PRECREG);

Vbase = vecextract_p(F.LP, F.perm);

class_group_gen(nf,W,C,Vbase,PRECREG,NULL, &clg1, &clg2);

c1 = gdiv(gmul(reg,clh), z);

res = get_clfu(clg1, R, gdiv(mpmul(R,h), z), zu, fu, k);

desallocate(&mat);

return gerepilecopy(av, buchall_end(nf,CHANGE,flun,res,clg2,W,B,A,C,Vbase));

return gerepileupto(av, buchall_end(nf,CHANGE,flun,k,fu,clg1,clg2,reg,c1,zu,W,B,xarch,C,vectbase,vperm));

}

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