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, |
|
Line 21 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
/*******************************************************************/ |
/*******************************************************************/ |
#include "pari.h" |
#include "pari.h" |
#include "parinf.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 GEN to_famat_all(GEN x, GEN y); |
extern int addcolumntomatrix(GEN V, GEN invp, GEN L); |
extern int addcolumntomatrix(GEN V, GEN invp, GEN L); |
extern double check_bach(double cbach, double B); |
extern double check_bach(double cbach, double B); |
extern GEN gmul_mat_smallvec(GEN x, GEN y); |
extern GEN gmul_mat_smallvec(GEN x, GEN y); |
extern GEN gmul_mati_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_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 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 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 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_MAX 2 |
#define SFB_STEP 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; |
static const int randshift = BITS_IN_RANDOM-1 - RANDOM_BITS; |
#undef 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 */ |
/* x a t_VECSMALL */ |
static long |
static long |
ccontent(GEN x) |
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; |
return s; |
} |
} |
|
|
static void |
static void |
desallocate(GEN **M) |
desallocate(GEN **M, GEN *first_nz) |
{ |
{ |
GEN *m = *M; |
GEN *m = *M; |
long i; |
long i; |
if (m) |
if (m) |
{ |
{ |
for (i=lg(m)-1; i; i--) free(m[i]); |
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 |
static GEN |
mulred(GEN nf,GEN x, GEN I, long prec) |
mulred(GEN nf,GEN x, GEN I, long prec) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
GEN y = cgetg(3,t_VEC); |
GEN y = cgetg(3,t_VEC); |
|
|
y[1] = (long)idealmulh(nf,I,(GEN)x[1]); |
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) |
/* 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 |
* 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 |
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); |
GEN *pow, arch0 = cgetg(RU,t_COL); |
for (i=1; i<RU; i++) arch0[i] = un; /* 0 in multiplicative notations */ |
for (i=1; i<RU; i++) arch0[i] = un; /* 0 in multiplicative notations */ |
|
|
if (DEBUGLEVEL) fprintferr("Computing powers for sub-factor base:\n"); |
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++) |
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]; |
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]=cgetg(3,t_VEC); |
pow[1][1] = (long)z; |
pow[1][1] = (long)z; |
pow[1][2] = (long)arch0; |
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"); |
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()) |
* 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)), |
* 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 |
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) |
{ |
{ |
if (DEBUGLEVEL>1) fprintferr("\n"); |
if (DEBUGLEVEL>1) fprintferr("\n"); |
if (DEBUGLEVEL>6) |
if (DEBUGLEVEL>6) |
{ |
{ |
fprintferr("########## FACTORBASE ##########\n\n"); |
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"); |
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; |
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++) |
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 |
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; |
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 |
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; |
char *s; |
switch(reason) |
switch(reason) |
{ |
{ |
case RELAT: s = "not enough relations for fundamental units"; break; |
|
case LARGE: s = "fundamental units too large"; break; |
case LARGE: s = "fundamental units too large"; break; |
case PRECI: s = "insufficient precision for fundamental units"; break; |
case PRECI: s = "insufficient precision for fundamental units"; break; |
default: s = "unknown problem with fundamental units"; |
default: s = "unknown problem with fundamental units"; |
} |
} |
err(warner,"%s, not given",s); |
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 */ |
/* check whether exp(x) will get too big */ |
static long |
static long |
expgexpo(GEN x) |
expgexpo(GEN x) |
{ |
{ |
long i,j,e, E = -HIGHEXPOBIT; |
long i,j,e, E = - (long)HIGHEXPOBIT; |
GEN p1; |
GEN p1; |
|
|
for (i=1; i<lg(x); i++) |
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); |
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"); |
if (DEBUGLEVEL) fprintferr("\n#### Computing fundamental units\n"); |
R1 = itos(gmael(nf,2,1)); RU = (N+R1)>>1; |
R1 = itos(gmael(nf,2,1)); RU = (N+R1)>>1; |
if (RU==1) { *pte=BIGINT; return cgetg(1,t_VEC); } |
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); |
matep = cgetg(RU,t_MAT); |
for (j=1; j<RU; j++) |
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); |
s = gdivgs(s, -N); |
p1=cgetg(RU+1,t_COL); matep[j]=(long)p1; |
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); |
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); |
matep = gexp(p1,prec); |
y = grndtoi(gauss_realimag(nf,matep), &e); |
y = grndtoi(gauss_realimag(nf,matep), &e); |
*pte = -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++) |
for (j=1; j<RU; j++) |
if (!gcmp1(idealnorm(nf, (GEN)y[j]))) break; |
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[i] are unit generators. Normalize: smallest L2 norm + lead coeff > 0 */ |
y = gmul((GEN)nf[7], y); |
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++) |
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), |
if (gcmp(QuickNormL2(p2,DEFAULTPREC), |
QuickNormL2(p1,DEFAULTPREC)) < 0) |
QuickNormL2(p1,DEFAULTPREC)) < 0) |
{ |
{ |
xarch[j] = lneg((GEN)xarch[j]); |
A[j] = lneg((GEN)A[j]); |
p1 = p2; |
p1 = p2; |
} |
} |
if (gsigne(leading_term(p1)) < 0) |
if (gsigne(leading_term(p1)) < 0) |
{ |
{ |
xarch[j] = ladd((GEN)xarch[j], vec); |
A[j] = ladd((GEN)A[j], vec); |
p1 = gneg(p1); |
p1 = gneg(p1); |
} |
} |
y[j] = (long)p1; |
y[j] = (long)p1; |
} |
} |
if (DEBUGLEVEL) msgtimer("getfu"); |
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 |
GEN |
buchfu(GEN bnf) |
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]; |
res = (GEN)bnf[8]; reg = (GEN)res[2]; |
if (lg(res)==7 && lg(res[5])==lg(nf[6])-1) |
if (lg(res)==7 && lg(res[5])==lg(nf[6])-1) |
{ |
{ |
y[1] = lcopy((GEN)res[5]); |
y[1] = lcopy((GEN)res[5]); |
y[2] = lcopy((GEN)res[6]); return y; |
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); |
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; |
*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 |
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 |
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); |
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++) |
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(;;) |
for(;;) |
{ |
{ |
int non0 = 0; |
gpmem_t av = avma; |
|
if (DEBUGLEVEL>2) fprintferr("# ideals tried = %ld\n",nbtest); |
id = x0; |
id = x0; |
for (i=1; i<lgsub; i++) |
for (i=1; i<lgsub; i++) |
{ |
{ |
ex[i] = mymyrand() >> randshift; |
ex[i] = mymyrand() >> randshift; |
if (ex[i]) |
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; |
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++) |
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; |
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 */ |
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++) |
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 |
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); |
RU = lg(mat); x = cgetg(RU+1,t_COL); |
for (i=1; i<RU; i++) x[i]=lreal((GEN)col[i]); |
for (i=1; i<RU; i++) x[i]=lreal((GEN)col[i]); |
x[RU] = (long)N2; |
x[RU] = (long)N2; |
x = lllintern(concatsp(mat,x), 1,prec); |
x = lllintern(concatsp(mat,x),100, 1,prec); |
if (!x) return NULL; |
if (!x) return NULL; |
x = (GEN)x[RU]; |
x = (GEN)x[RU]; |
if (signe(x[RU]) < 0) x = gneg_i(x); |
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 |
static GEN |
act_arch(GEN A, GEN x) |
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 */ |
/* 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++) |
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; |
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]); |
N = degpol(nf[1]); |
R1 = itos(gmael(nf,2,1)); |
R1 = itos(gmael(nf,2,1)); |
RU = (N + R1)>>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) |
if (RU > 1) |
{ /* reduce mod units */ |
{ /* reduce mod units */ |
GEN u, z = init_red_mod_units(bnf,prec); |
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 |
static int |
fact_ok(GEN nf, GEN y, GEN C, GEN g, GEN e) |
fact_ok(GEN nf, GEN y, GEN C, GEN g, GEN e) |
{ |
{ |
|
gpmem_t av = avma; |
long i, c = lg(e); |
long i, c = lg(e); |
GEN z = C? C: gun; |
GEN z = C? C: gun; |
for (i=1; i<c; i++) |
for (i=1; i<c; i++) |
if (signe(e[i])) z = idealmul(nf, z, idealpow(nf, (GEN)g[i], (GEN)e[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(z) != t_MAT) z = idealhermite(nf,z); |
if (typ(y) != t_MAT) y = idealhermite(nf,y); |
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 */ |
/* assume x in HNF. cf class_group_gen for notations */ |
static GEN |
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; |
long i,lW,lB,e,c, prec = *ptprec; |
GEN Q,xar,Wex,Bex,U,y,p1,gen,cyc,xc,ex,d,col,A; |
GEN Q,xar,Wex,Bex,U,y,p1,gen,cyc,xc,ex,d,col,A; |
GEN W = (GEN)bnf[1]; |
GEN W = (GEN)bnf[1]; |
GEN B = (GEN)bnf[2]; |
GEN B = (GEN)bnf[2]; |
GEN WB_C = (GEN)bnf[4]; |
GEN WB_C = (GEN)bnf[4]; |
GEN vperm = (GEN)bnf[6]; |
|
GEN nf = (GEN)bnf[7]; |
GEN nf = (GEN)bnf[7]; |
GEN clg2 = (GEN)bnf[9]; |
GEN clg2 = (GEN)bnf[9]; |
int old_format = (typ(clg2[2]) == t_MAT); |
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; |
cyc = gmael3(bnf,8,1,2); c = lg(cyc)-1; |
gen = gmael3(bnf,8,1,3); |
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 */ |
/* 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] |
/* x = g_W Wex + g_B Bex - [xar] |
* = g_W A - [xar] + [C_B]Bex since g_W B + g_B = [C_B] */ |
* = 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); |
Q = gmul(U, A); |
ex = cgetg(c+1,t_COL); |
|
for (i=1; i<=c; i++) |
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 */ |
/* compute arch component of the missing principal ideal */ |
if (old_format) |
if (old_format) |
{ |
{ |
GEN Garch, V = (GEN)clg2[2]; |
GEN Garch, V = (GEN)clg2[2]; |
|
Bex = small_to_col(Bex); |
p1 = c? concatsp(gmul(V,Q), Bex): Bex; |
p1 = c? concatsp(gmul(V,Q), Bex): Bex; |
col = act_arch(p1, WB_C); |
col = act_arch(p1, WB_C); |
if (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 A = G Ur A + [ga]A, Ur A = D Q + R as above (R = ex) |
= G R + [GD]Q + [ga]A */ |
= G R + [GD]Q + [ga]A */ |
GEN ga = (GEN)clg2[2], GD = (GEN)clg2[3]; |
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 (lW) col = gadd(col, act_arch(A, ga)); |
if (c) col = gadd(col, act_arch(Q, GD)); |
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 */ |
/* 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)); |
Q = gdiv(dethnf_i(x), get_norm_fact(gen, ex, &d)); |
col = isprincipalarch(bnf, col, Q, gun, d, &e); |
col = isprincipalarch(bnf, col, Q, gun, d, &e); |
if (col && !fact_ok(nf,x, col,gen,ex)) col = NULL; |
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) |
if (!col) |
{ |
{ |
*ptprec = prec + (e >> TWOPOTBITS_IN_LONG) + (MEDDEFAULTPREC-2); |
*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; |
e = 0; |
} |
} |
y = cgetg(4,t_VEC); |
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[1] = lcopy(ex); |
y[2] = e? lmul(xc,col): lgetg(1,t_COL); |
y[2] = (long)col; |
y[3] = lstoi(-e); return y; |
y[3] = lstoi(-e); return y; |
} |
} |
|
|
|
|
triv_gen(GEN nf, GEN x, long c, long flag) |
triv_gen(GEN nf, GEN x, long c, long flag) |
{ |
{ |
GEN y; |
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 = cgetg(4,t_VEC); |
y[1] = (long)zerocol(c); |
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; |
y[3] = lstoi(BIGINT); return y; |
} |
} |
|
|
GEN |
GEN |
isprincipalall(GEN bnf,GEN x,long flag) |
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; |
GEN nf; |
|
|
bnf = checkbnf(bnf); nf = (GEN)bnf[7]; |
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"); |
err(talker,"not the same number field in isprincipal"); |
x = (GEN)x[2]; tx = t_POL; |
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"); |
if (gcmp0(x)) err(talker,"zero ideal in isprincipal"); |
return triv_gen(nf, x, lg(mael3(bnf,8,1,2))-1, flag); |
return triv_gen(nf, x, lg(mael3(bnf,8,1,2))-1, flag); |
} |
} |
x = idealhermite(nf,x); |
x = idealhermite(nf,x); |
if (lg(x)==1) err(talker,"zero ideal in isprincipal"); |
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)); |
return gerepileupto(av, triv_gen(nf, gcoeff(x,1,1), 0, flag)); |
|
|
pr = prec_arch(bnf); /* precision of unit matrix */ |
pr = prec_arch(bnf); /* precision of unit matrix */ |
c = getrand(); |
c = getrand(); |
for (;;) |
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 (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 |
GEN |
isprincipalfact(GEN bnf,GEN P, GEN e, GEN C, long flag) |
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 */ |
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); |
prec = prec_arch(bnf); |
if (gen) |
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); |
id2 = idealpowred(bnf,z, (GEN)e[i],prec); |
id = id? idealmulred(nf,id,id2,prec): id2; |
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(); |
c = getrand(); |
for (;;) |
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 (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) |
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"); |
err(warner,"insufficient precision for generators, not given"); |
avma = av; return stoi(prec); |
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 |
GEN |
isprincipal(GEN bnf,GEN x) |
isprincipal(GEN bnf,GEN x) |
{ |
{ |
return isprincipalall(bnf,x,nf_REGULAR); |
return isprincipalall(bnf,x,0); |
} |
} |
|
|
GEN |
GEN |
Line 1222 isprincipalgenforce(GEN bnf,GEN x) |
|
Line 1361 isprincipalgenforce(GEN bnf,GEN x) |
|
return isprincipalall(bnf,x,nf_GEN | nf_FORCE); |
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 |
GEN |
isunit(GEN bnf,GEN x) |
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]; |
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 */ |
p1 = gmael(bnf,8,4); /* roots of 1 */ |
gn = (GEN)p1[1]; n = itos(gn); |
gn = (GEN)p1[1]; n = itos(gn); |
z = (GEN)p1[2]; |
z = algtobasis(nf, (GEN)p1[2]); |
switch(tx) |
switch(tx) |
{ |
{ |
case t_INT: case t_FRAC: case t_FRACN: |
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: |
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 */ |
/* 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; |
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(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 = 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) */ |
/* 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)); |
p1 = ground(gdiv(p1, pi2_sur_w)); |
if (n > 2) |
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)); |
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 |
GEN |
signunits(GEN bnf) |
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; |
GEN matunit,y,p1,p2,nf,pi; |
|
|
bnf = checkbnf(bnf); nf = (GEN)bnf[7]; |
bnf = checkbnf(bnf); nf = (GEN)bnf[7]; |
Line 1326 signunits(GEN bnf) |
|
Line 1488 signunits(GEN bnf) |
|
return y; |
return y; |
} |
} |
|
|
/* LLL-reduce ideal and return T2 | ideal */ |
/* LLL-reduce ideal and return Cholesky for T2 | ideal */ |
static GEN |
static GEN |
red_ideal(GEN *ideal,GEN T2vec,GEN prvec) |
red_ideal(GEN *ideal,GEN Gvec,GEN prvec) |
{ |
{ |
long i; |
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]; |
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); |
if (DEBUGLEVEL) err(warner, "prec too low in red_ideal[%ld]: %ld",i,prec); |
} |
} |
return NULL; |
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 |
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) |
set_log_embed(GEN col, GEN x, long R1, long prec) |
{ |
{ |
long i, l = lg(x); |
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 |
static void |
set_fact(GEN c) |
set_fact(GEN c, GEN first_nz, long cglob) |
{ |
{ |
long i; |
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 |
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; |
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 |
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; |
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) |
if (DEBUGLEVEL) |
fprintferr("\n#### Looking for %ld relations (small norms)\n",nbrel); |
fprintferr("\n#### Looking for %ld relations (small norms)\n",nbrel); |
xembed = NULL; /* gcc -Wall */ |
xembed = NULL; /* gcc -Wall */ |
nbsmallnorm = nbfact = 0; |
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); |
minim_alloc(N+1, &q, &x, &y, &z, &v); |
av1 = avma; |
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); |
if (DEBUGLEVEL>1) fprintferr("\n*** Ideal no %ld: %Z\n", noideal, ideal); |
ideal = prime_to_ideal(nf,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 */ |
if (!r) return -1; /* precision problem */ |
|
|
for (k=1; k<=N; k++) |
for (k=1; k<=N; k++) |
{ |
{ |
v[k]=gtodouble(gcoeff(r,k,k)); |
v[k]=gtodouble(gcoeff(r,k,k)); |
for (j=1; j<k; j++) q[j][k]=gtodouble(gcoeff(r,j,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>1) |
{ |
{ |
if (DEBUGLEVEL>3) fprintferr("\n"); |
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]); |
k = N; y[N]=z[N]=0; x[N]= (long) sqrt(BOUND/v[N]); |
for(;; x[1]--) |
for(;; x[1]--) |
{ |
{ |
ulong av3 = avma; |
gpmem_t av3 = avma; |
double p; |
double p; |
GEN col; |
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 (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 */ |
if (ccontent(x)==1) /* primitive */ |
{ |
{ |
GEN Nx, gx = gmul_mati_smallvec(IDEAL,x); |
GEN Nx, gx = gmul_mati_smallvec(IDEAL,x); |
long av4; |
gpmem_t av4; |
if (!isnfscalar(gx)) |
if (!isnfscalar(gx)) |
{ |
{ |
xembed = gmul(M,gx); av4 = avma; nbsmallnorm++; |
xembed = gmul(M,gx); av4 = avma; nbsmallnorm++; |
|
if (++try_factor > maxtry_FACT) goto ENDIDEAL; |
Nx = ground(norm_by_embed(R1,xembed)); setsigne(Nx, 1); |
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(); } |
if (DEBUGLEVEL > 1) { fprintferr("."); flusherr(); } |
} |
} |
avma = av3; |
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]; |
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 */ |
{ /* Q-dependent from previous ones: forget it */ |
cglob--; unset_fact(col); |
cglob--; unset_fact(col); |
if (DEBUGLEVEL>1) { fprintferr("*"); flusherr(); } |
if (DEBUGLEVEL>1) { fprintferr("*"); flusherr(); } |
if (++dependent > MAXTRY) { alldep++; break; } |
if (++dependent > maxtry_DEP) break; |
avma = av3; continue; |
avma = av3; continue; |
} |
} |
/* record archimedean part */ |
/* record archimedean part */ |
set_log_embed((GEN)matarch[cglob], xembed, R1, PRECREG); |
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 (cglob >= nbrel) goto END; /* we have enough */ |
if (++nbrelideal == nbrelpid) break; |
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: |
ENDIDEAL: |
invp = gerepilecopy(av1, invp); |
if (cglob == oldcglob) avma = av0; else invp = gerepilecopy(av1, invp); |
if (DEBUGLEVEL>1) msgtimer("for this ideal"); |
if (DEBUGLEVEL>1) msgtimer("for this ideal"); |
} |
} |
END: |
END: |
if (DEBUGLEVEL) |
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); |
nbfact,nbsmallnorm,((double)nbfact)/nbsmallnorm); |
else |
|
fprintferr("\nnb. small norm = 0\n"); |
|
} |
} |
avma = av; return cglob; |
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 |
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; |
if (!y) return NULL; |
|
|
/* I, m, y integral */ |
|
m = gmul(I,(GEN)y[1]); |
m = gmul(I,(GEN)y[1]); |
if (isnfscalar(m)) m = gmul(I,(GEN)y[2]); |
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 |
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 |
static void |
dbg_cancelrel(long jideal,long jdir,long phase, long *col) |
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); |
if (DEBUGLEVEL>3) fprintferr("(jideal=%ld,jdir=%ld)",jideal,jdir); |
wr_rel(col); flusherr(); |
wr_rel(col); flusherr(); |
} |
} |
|
|
static void |
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 */ |
* General check for colinearity useless since exceedingly rare */ |
static long |
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; |
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 */ |
/* I integral ideal in HNF form */ |
|
|
remove_content(GEN I) |
remove_content(GEN I) |
{ |
{ |
long N = lg(I)-1; |
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; |
return I; |
} |
} |
|
|
/* if phase != 1 re-initialize static variables. If <0 return immediately */ |
/* if phase != 1 re-initialize static variables. If <0 return immediately */ |
static long |
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; |
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; } |
if (phase != 1) { jideal=jdir=1; if (phase<0) return 0; } |
|
|
r1 = nf_get_r1(nf); |
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) |
if (DEBUGLEVEL && list_jideal) |
fprintferr("looking hard for %Z\n",list_jideal); |
fprintferr("looking hard for %Z\n",list_jideal); |
P = NULL; /* gcc -Wall */ |
P = NULL; /* gcc -Wall */ |
for (av = avma;;) |
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; |
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; |
ideal = P; |
do { |
do { |
for (i=1; i<lgsub; i++) |
for (i=1; i<lgsub; i++) |
{ |
{ |
ex[i] = mymyrand()>>randshift; |
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 */ |
while (ideal == P); /* If ex = 0, try another */ |
ideal = remove_content(ideal); |
ideal = remove_content(ideal); |
|
|
if (phase != 1) jdir = 1; else phase = 2; |
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 */ |
{ /* reduce along various directions */ |
if (DEBUGLEVEL>2) |
if (DEBUGLEVEL>2) |
fprintferr("phase=%ld,jideal=%ld,jdir=%ld,rand=%ld\n", |
fprintferr("phase=%ld,jideal=%ld,jdir=%ld,rand=%ld\n", |
phase,jideal,jdir,getrand()); |
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(); } |
if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); } |
continue; |
continue; |
} |
} |
/* can factor ideal, record relation */ |
/* can factor ideal, record relation */ |
cglob++; col = mat[cglob]; |
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 */ |
{ /* Already known: forget it */ |
if (DEBUGLEVEL>1) dbg_cancelrel(jideal,jdir,phase,col); |
if (DEBUGLEVEL>1) dbg_cancelrel(jideal,jdir,phase,col); |
cglob--; unset_fact(col); col[jideal] = 0; |
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) |
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++) |
for (i=1; i<lgsub; i++) |
if (ex[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; |
arch = arch? vecmul(arch, p1): p1; |
} |
} |
colarch = (GEN)matarch[cglob]; |
colarch = (GEN)matarch[cglob]; |
/* arch = archimedean component (MULTIPLICATIVE form) of ideal */ |
/* 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); |
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 */ |
/* Need more, try next P */ |
if (cglob < lim) break; |
if (cglob < maxcglob) break; |
|
|
/* We have found enough. Return */ |
/* We have found enough. Return */ |
if (phase) |
if (phase) |
{ |
{ |
jdir = 1; |
jdir = 1; |
if (jideal == KC) jideal=1; else jideal++; |
if (jideal == F->KC) jideal=1; else jideal++; |
} |
} |
if (DEBUGLEVEL>2) |
if (DEBUGLEVEL>2) |
fprintferr("Upon exit: jideal=%ld,jdir=%ld\n",jideal,jdir); |
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 (!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) |
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(); |
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; |
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++) |
for (j=1; j<J; j++) |
{ |
{ |
GEN ideal0 = prime_to_ideal(nf,(GEN)P[j]); |
GEN ideal0 = prime_to_ideal(nf,(GEN)P[j]); |
ulong av2 = avma; |
gpmem_t av2 = avma; |
for(nbtest=0;;) |
for(nbtest=0;;) |
{ |
{ |
ideal = ideal0; |
ideal = ideal0; |
for (i=1; i<lgsub; i++) |
for (i=1; i<lgsub; i++) |
{ |
{ |
ex = mymyrand()>>randshift; |
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); |
ideal = remove_content(ideal); |
|
for (i=1; i<ru; i++) vdir[i] = mymyrand()>>randshift; |
for (k=1; k<ru; k++) |
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; |
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) |
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"); |
if (DEBUGLEVEL>1) fprintferr("\n"); |
msgtimer("be honest"); |
msgtimer("be honest"); |
} |
} |
|
F->KCZ = KCZ0; |
avma = av; return 1; |
avma = av; return 1; |
} |
} |
|
|
Line 1851 trunc_error(GEN x) |
|
Line 1968 trunc_error(GEN x) |
|
&& (expo(x)>>TWOPOTBITS_IN_LONG) + 3 > lg(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 |
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 */ |
i = gprecision(mdet); /* truncate to avoid "near dependant" vectors */ |
mdet_t = (i <= 4)? mdet: gprec_w(mdet,i-1); |
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); |
Im_mdet = vecextract_p(mdet,v); |
/* integral multiple of R: the cols we picked form a Q-basis, they have an |
/* 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); |
kR = gdivgs(det2(Im_mdet), N); |
/* R > 0.2 uniformly */ |
/* R > 0.2 uniformly */ |
if (gexpo(kR) < -3) { avma=av; return NULL; } |
if (gexpo(kR) < -3) { avma=av; return NULL; } |
|
|
kR = mpabs(kR); |
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 |
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(); } |
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) */ |
/* find the smallest (wrt norm) among I, I^-1 and red(I^-1) */ |
|
|
inverse_if_smaller(GEN nf, GEN I, long prec) |
inverse_if_smaller(GEN nf, GEN I, long prec) |
{ |
{ |
GEN J, d, dmin, I1; |
GEN J, d, dmin, I1; |
int inv; |
|
J = (GEN)I[1]; |
J = (GEN)I[1]; |
dmin = dethnf_i(J); I1 = idealinv(nf,I); |
dmin = dethnf_i(J); I1 = idealinv(nf,I); |
J = (GEN)I1[1]; J = gmul(J,denom(J)); I1[1] = (long)J; |
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) */ |
/* try reducing (often _increases_ the norm) */ |
I1 = ideallllred(nf,I1,NULL,prec); |
I1 = ideallllred(nf,I1,NULL,prec); |
J = (GEN)I1[1]; |
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; |
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) */ |
/* compute class group (clg1) + data for isprincipal (clg2) */ |
static void |
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) |
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); |
lo0 = lo = lg(D); |
/* we could set lo = lg(cyc) and truncate all matrices below |
/* we could set lo = lg(cyc) and truncate all matrices below |
* setlg_col(D && U && Y, lo) + setlg(D && V && X && Ui, lo) |
* 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 |
* P invertible diagonal matrix (\pm 1) which is only implicitly defined |
* G = g Uir P + [Ga], Uir = Ui + WX |
* G = g Uir P + [Ga], Uir = Ui + WX |
* g = G P Ur + [ga], Ur = U + DY */ |
* 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); |
G = cgetg(lo,t_VEC); |
Ga= 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++) |
for (j=1; j<lo; j++) |
{ |
{ |
GEN p1 = gcoeff(Uir,1,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++) |
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) |
if (J != I) |
{ /* update wrt P */ |
{ /* update wrt P */ |
neg_row(Y ,j); V[j] = lneg((GEN)V[j]); |
neg_row(Y ,j); V[j] = lneg((GEN)V[j]); |
neg_row(Ur,j); X[j] = lneg((GEN)X[j]); |
neg_row(Ur,j); X[j] = lneg((GEN)X[j]); |
} |
} |
G[j] = (long)J[1]; /* generator, order cyc[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 */ |
/* 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"); |
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 |
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() |
/* 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 |
* 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 |
* 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 |
* 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 |
static GEN |
relationrank(GEN *A, long r, GEN L) |
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); |
GEN invp = idmat(n); |
|
|
if (!r) return invp; |
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); |
return gerepilecopy(av, invp); |
} |
} |
|
|
|
/* SMALLBUCHINIT */ |
|
|
static GEN |
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 |
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] */ |
/* 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 && !fact_ok(nf,y,NULL,gen,(GEN)D[i])) y = NULL; |
if (y) { h[i] = (long)to_famat_all(y,gun); continue; } |
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; |
return h; |
} |
} |
Line 2243 makecycgen(GEN bnf) |
|
Line 2395 makecycgen(GEN bnf) |
|
static GEN |
static GEN |
makematal(GEN bnf) |
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; |
long lW,lma,j,prec; |
|
|
ma = get_matal((GEN)bnf[10]); |
ma = get_matal((GEN)bnf[10]); |
Line 2252 makematal(GEN bnf) |
|
Line 2404 makematal(GEN bnf) |
|
W = (GEN)bnf[1]; |
W = (GEN)bnf[1]; |
B = (GEN)bnf[2]; |
B = (GEN)bnf[2]; |
WB_C= (GEN)bnf[4]; |
WB_C= (GEN)bnf[4]; |
vp = (GEN)bnf[6]; |
|
nf = (GEN)bnf[7]; |
nf = (GEN)bnf[7]; |
lW=lg(W)-1; lma=lW+lg(B); |
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); |
ma = cgetg(lma,t_MAT); |
|
|
prec = prec_arch(bnf); |
prec = prec_arch(bnf); |
Line 2274 makematal(GEN bnf) |
|
Line 2424 makematal(GEN bnf) |
|
ma[j] = (long)y; continue; |
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 (typ(y) != t_INT) |
{ |
{ |
if (DEBUGLEVEL>1) fprintferr("%ld ",j); |
if (DEBUGLEVEL>1) fprintferr("%ld ",j); |
Line 2284 makematal(GEN bnf) |
|
Line 2434 makematal(GEN bnf) |
|
prec = itos(y); j--; |
prec = itos(y); j--; |
if (DEBUGLEVEL) err(warnprec,"makematal",prec); |
if (DEBUGLEVEL) err(warnprec,"makematal",prec); |
nf = nfnewprec(nf,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"); |
if (DEBUGLEVEL>1) fprintferr("\n"); |
return ma; |
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]); |
GEN cycgen = get_cycgen((GEN)bnf[10]); |
if (!cycgen) |
if (!cycgen) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
if (DEBUGLEVEL) err(warner,"completing bnf (building cycgen)"); |
if (DEBUGLEVEL) err(warner,"completing bnf (building cycgen)"); |
bnfinsert(bnf, makecycgen(bnf), 2); avma = av; |
bnfinsert(bnf, makecycgen(bnf), 2); avma = av; |
cycgen = get_cycgen((GEN)bnf[10]); |
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]); |
GEN matal = get_matal((GEN)bnf[10]); |
if (!matal) |
if (!matal) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
if (DEBUGLEVEL) err(warner,"completing bnf (building matal)"); |
if (DEBUGLEVEL) err(warner,"completing bnf (building matal)"); |
bnfinsert(bnf, makematal(bnf), 1); avma = av; |
bnfinsert(bnf, makematal(bnf), 1); avma = av; |
matal = get_matal((GEN)bnf[10]); |
matal = get_matal((GEN)bnf[10]); |
} |
} |
return matal; |
return matal; |
} |
} |
|
|
GEN |
GEN |
smallbuchinit(GEN pol,GEN gcbach,GEN gcbach2,GEN gRELSUP,GEN gborne,long nbrelpid,long minsFB,long prec) |
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); |
if (typ(pol)==t_VEC) bnf = checkbnf(pol); |
else |
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); |
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 |
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; |
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 |
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++) |
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; |
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 |
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 |
GEN |
bnfnewprec(GEN bnf, long prec) |
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); |
bnf = checkbnf(bnf); |
if (prec <= 0) return nfnewprec(checknf(bnf),prec); |
if (prec <= 0) return nfnewprec(checknf(bnf),prec); |
y = cgetg(11,t_VEC); |
|
funits = check_units(bnf,"bnfnewprec"); |
|
nf = (GEN)bnf[7]; |
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; |
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; } |
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); |
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 |
GEN |
Line 2505 bnrnewprec(GEN bnr, long prec) |
|
Line 2607 bnrnewprec(GEN bnr, long prec) |
|
return y; |
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 |
GEN |
bnfmake(GEN sbnf, long prec) |
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]); |
matal = get_matal((GEN)sbnf[12]); |
if (!matal) 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]; |
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=cgetg(11,t_VEC); |
y[1]=(long)W; y[2]=sbnf[8]; y[3]=(long)mun; |
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 |
static GEN |
classgroupall(GEN P, GEN data, long flag, long prec) |
classgroupall(GEN P, GEN data, long flag, long prec) |
{ |
{ |
long court[3],doubl[4]; |
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; |
GEN bach=doubl,bach2=doubl,RELSUP=court,borne=gun; |
|
|
if (!data) lx=1; |
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) |
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 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"); |
default: err(flagerr,"classgroupall"); |
return NULL; /* not reached */ |
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 |
GEN |
Line 2633 bnfinit0(GEN P, long flag, GEN data, long prec) |
|
Line 2745 bnfinit0(GEN P, long flag, GEN data, long prec) |
|
GEN |
GEN |
classgrouponly(GEN P, GEN data, long prec) |
classgrouponly(GEN P, GEN data, long prec) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
GEN z; |
GEN z; |
|
|
if (typ(P)==t_INT) |
if (typ(P)==t_INT) |
Line 2648 classgrouponly(GEN P, GEN data, long prec) |
|
Line 2760 classgrouponly(GEN P, GEN data, long prec) |
|
GEN |
GEN |
regulator(GEN P, GEN data, long prec) |
regulator(GEN P, GEN data, long prec) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
GEN z; |
GEN z; |
|
|
if (typ(P)==t_INT) |
if (typ(P)==t_INT) |
|
|
GEN c = (GEN) gpmalloc(sizeof(long)*(n+1)); |
GEN c = (GEN) gpmalloc(sizeof(long)*(n+1)); |
long i; |
long i; |
for (i=1; i<=n; i++) c[i]=0; |
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; |
return c; |
} |
} |
|
|
Line 2687 cgetc_col(long n, long prec) |
|
Line 2799 cgetc_col(long n, long prec) |
|
} |
} |
|
|
static GEN |
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=cgetg(11,t_VEC); |
z[1]=(long)W; |
z[1]=(long)W; |
z[2]=(long)B; |
z[2]=(long)B; |
z[3]=(long)xarch; |
z[3]=(long)A; |
z[4]=(long)matarch; |
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[9]=(long)clg2; |
z[10]=zero; /* dummy: we MUST have lg(bnf) != lg(nf) */ |
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; } |
if (CHANGE) { p1=cgetg(3,t_VEC); p1[1]=(long)z; p1[2]=(long)CHANGE; z=p1; } |
return gcopy(z); |
return z; |
} |
} |
|
|
static GEN |
static GEN |
buchall_for_degree_one_pol(GEN nf, GEN CHANGE, long flun) |
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); |
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); |
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 |
GEN |
buchall(GEN P,GEN gcbach,GEN gcbach2,GEN gRELSUP,GEN gborne,long nbrelpid, |
buchall(GEN P,GEN gcbach,GEN gcbach2,GEN gRELSUP,GEN gborne,long nbrelpid, |
long minsFB,long flun,long prec) |
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; |
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); |
if (typ(gRELSUP) != t_INT) gRELSUP = gtrunc(gRELSUP); |
RELSUP = itos(gRELSUP); |
RELSUP = itos(gRELSUP); |
if (RELSUP<=0) err(talker,"not enough relations in bnfxxx"); |
if (RELSUP<=0) err(talker,"not enough relations in bnfxxx"); |
|
|
/* Initializations */ |
/* Initializations */ |
fu = NULL; /* gcc -Wall */ |
fu = NULL; /* gcc -Wall */ |
N = degpol(P); PRECREG = max(BIGDEFAULTPREC,prec); |
N = degpol(P); |
|
PRECREG = max(BIGDEFAULTPREC,prec); |
|
PRECLLLadd = DEFAULTPREC; |
if (!nf) |
if (!nf) |
{ |
{ |
nf = initalgall0(P, nf_REGULAR, PRECREG); |
nf = initalg(P, PRECREG); |
/* P was non-monic and nfinit changed it ? */ |
/* P was non-monic and nfinit changed it ? */ |
if (lg(nf)==3) { CHANGE = (GEN)nf[2]; nf = (GEN)nf[1]; } |
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 = rootsof1(nf); |
zu[2] = lmul((GEN)nf[7],(GEN)zu[2]); |
zu[2] = lmul((GEN)nf[7],(GEN)zu[2]); |
if (DEBUGLEVEL) msgtimer("initalg & rootsof1"); |
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)); |
D = (GEN)nf[3]; drc = fabs(gtodouble(D)); |
LOGD2 = log(drc); LOGD2 = LOGD2*LOGD2; |
LOGD2 = log(drc); LOGD2 = LOGD2*LOGD2; |
lim = (long) (exp(-(double)N) * sqrt(2*PI*N*drc) * pow(4/PI,(double)R2)); |
lim = (long) (exp(-(double)N) * sqrt(2*PI*N*drc) * pow(4/PI,(double)R2)); |
if (lim < 3) lim = 3; |
if (lim < 3) lim = 3; |
cbach = min(12., gtodouble(gcbach)); cbach /= 2; |
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); |
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) |
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: |
START: |
avma = av0; |
seed = getrand(); |
if (first) first = 0; else desallocate(&mat); |
avma = av0; desallocate(&mat, &first_nz); |
if (precpb) |
if (precpb) |
{ |
{ |
precdouble++; |
precdouble++; |
|
|
cbach = check_bach(cbach,12.); |
cbach = check_bach(cbach,12.); |
precadd = 0; |
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)); |
LIMC2= max(3 * N, (long)(cbach2*LOGD2)); |
if (LIMC2 < LIMC) LIMC2 = LIMC; |
if (LIMC2 < LIMC) LIMC2 = LIMC; |
if (DEBUGLEVEL) { fprintferr("LIMC = %ld, LIMC2 = %ld\n",LIMC,LIMC2); } |
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; |
sfb_trials = sfb_increase = 0; |
subFB = subFBgen(N,min(lim,LIMC2), minsFB, vperm, &ss); |
|
if (!subFB) goto START; |
|
nreldep = nrelsup = 0; |
nreldep = nrelsup = 0; |
|
|
/* PRECLLL = prec for LLL-reductions (idealred) |
/* PRECLLL = prec for LLL-reductions (idealred) |
* PRECREG = prec for archimedean components */ |
* 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 (!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) |
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 */ |
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); |
setlg(mat, KCCO+1); |
C = cgetg(KCCO+1,t_MAT); |
C = cgetg(KCCO+1,t_MAT); |
/* trivial relations (p) = prod P^e */ |
/* trivial relations (p) = prod P^e */ |
cglob = 0; z = zerocol(RU); |
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); |
if (DEBUGLEVEL) fprintferr("After trivial relations, cglob = %ld\n",cglob); |
/* initialize for other relations */ |
/* initialize for other relations */ |
for (i=cglob+1; i<=KCCO; i++) |
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); |
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); |
invp = relationrank(mat,cglob,liste); |
|
|
/* relations through elements of small norm */ |
/* 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; } |
if (cglob < 0) { precpb = "small_norm"; goto START; } |
avma = av1; limpile = stack_lim(av1,1); |
avma = av1; limpile = stack_lim(av1,1); |
|
|
phase = 0; |
phase = 0; |
nlze = matmax = 0; /* for lint */ |
nlze = 0; /* for lint */ |
vecT2 = NULL; |
vecG = NULL; |
|
list_jideal = NULL; |
|
|
/* random relations */ |
/* random relations */ |
if (cglob == KCCO) /* enough rels, but init random_relations just in case */ |
if (cglob == KCCO) /* enough rels, but init random_relations just in case */ |
|
|
else |
else |
{ |
{ |
GEN matarch; |
GEN matarch; |
|
long ss; |
|
|
if (DEBUGLEVEL) fprintferr("\n#### Looking for random relations\n"); |
if (DEBUGLEVEL) fprintferr("\n#### Looking for random relations\n"); |
MORE: |
MORE: |
|
|
{ |
{ |
if (DEBUGLEVEL) fprintferr("*** Increasing sub factor base\n"); |
if (DEBUGLEVEL) fprintferr("*** Increasing sub factor base\n"); |
sfb_increase = 0; |
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; |
nreldep = nrelsup = 0; |
} |
} |
|
|
if (phase == 0) matarch = C; |
if (phase == 0) matarch = C; |
else |
else |
{ /* reduced the relation matrix at least once */ |
{ /* 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(extraC, lgex+1); |
setlg(extramat,lgex+1); /* were allocated after hnfspec */ |
setlg(extramat,lgex+1); /* were allocated after hnfspec */ |
slim = cglob+lgex; |
slim = cglob+lgex; |
if (slim > matmax) |
if (slim > matmax) |
{ |
{ |
mat = (long**)gprealloc(mat, (2*slim+1) * sizeof(long*), |
|
(matmax+1) * sizeof(long*)); |
|
matmax = 2 * slim; |
matmax = 2 * slim; |
|
reallocate(matmax, (GEN*)&mat, &first_nz); |
} |
} |
setlg(mat, slim+1); |
setlg(mat, slim+1); |
if (DEBUGLEVEL) |
if (DEBUGLEVEL) |
fprintferr("\n(need %ld more relation%s)\n", lgex, (lgex==1)?"":"s"); |
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 */ |
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; |
av1 = avma; |
} |
} |
if (!powsubFB) |
if (!F.powsubFB) |
{ |
{ |
powsubFBgen(nf,subFB,CBUCHG+1,PRECREG); |
powsubFBgen(&F, nf, CBUCHG+1, PRECREG); |
av1 = avma; |
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) |
if (ss < 0) |
{ /* could not find relations */ |
{ /* 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; |
goto START; |
} |
} |
if (DEBUGLEVEL > 2) dbg_outrel(phase,cglob,vperm,mat,matarch); |
if (DEBUGLEVEL > 2 && phase == 0) dbg_outrel(cglob,mat,matarch); |
if (phase) |
if (phase) |
for (j=lg(extramat)-1; j>0; j--) |
for (j=lg(extramat)-1; j>0; j--) |
{ |
{ |
GEN c = mat[cglob+j], *g = (GEN*) extramat[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; |
cglob = ss; |
} |
} |
|
|
if (phase == 0) |
if (phase == 0) |
{ /* never reduced before */ |
{ /* never reduced before */ |
long lgex; |
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; |
phase = 1; |
nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-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 */ |
lgex = max(nlze, MIN_EXTRA); /* set lgex > 0, in case regulator is 0 */ |
/* allocate now, reduce dimension later (setlg) when lgex goes down */ |
/* allocate now, reduce dimension later (setlg) when lgex goes down */ |
extramat= cgetg(lgex+1,t_MAT); |
extramat= cgetg(lgex+1,t_MAT); |
extraC = cgetg(lgex+1,t_MAT); |
extraC = cgetg(lgex+1,t_MAT); |
for (j=1; j<=lgex; j++) |
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); |
extraC[j] = (long)cgetc_col(RU,PRECREG); |
} |
} |
} |
} |
|
|
{ |
{ |
if (low_stack(limpile, stack_lim(av1,1))) |
if (low_stack(limpile, stack_lim(av1,1))) |
{ |
{ |
GEN *gptr[6]; |
|
if(DEBUGMEM>1) err(warnmem,"buchall"); |
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; |
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; |
nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1; |
if (nlze && ++nreldep > MAXRELSUP) { sfb_increase=1; goto MORE; } |
if (nlze && ++nreldep > MAXRELSUP) { sfb_increase=1; goto MORE; } |
} |
} |
if (nlze) goto MORE; /* dependent rows */ |
if (nlze) goto MORE; /* dependent rows */ |
|
|
/* first attempt at regulator for the check */ |
/* 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 */ |
{ /* not full rank for units */ |
if (DEBUGLEVEL) fprintferr("regulator is zero.\n"); |
if (DEBUGLEVEL) fprintferr("regulator is zero.\n"); |
if (++nrelsup > MAXRELSUP) goto START; |
if (++nrelsup > MAXRELSUP) goto START; |
nlze = MIN_EXTRA; goto MORE; |
nlze = MIN_EXTRA; goto MORE; |
} |
} |
/* anticipate precision problems */ |
/* 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,®); |
|
/* 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 */ |
/* 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; |
precpb = "getfu"; goto START; |
} |
} |
} |
} |
|
desallocate(&mat, &first_nz); |
|
|
/* class group generators */ |
/* 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)); |
|
} |
} |