version 1.1, 2001/10/02 11:17:03 |
version 1.2, 2002/09/11 07:26:50 |
Line 23 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
Line 23 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
|
|
const int narrow = 0; /* should set narrow = flag in buchquad, but buggy */ |
const int narrow = 0; /* should set narrow = flag in buchquad, but buggy */ |
|
|
|
/* For largeprime() hashtable. Note that hashed pseudoprimes are odd (unless |
|
* 2 | index), hence the low order bit is not useful. So we hash |
|
* HASHBITS bits starting at bit 1, not bit 0 */ |
|
#define HASHBITS 10 |
|
const int HASHT = 1 << HASHBITS; |
|
|
|
static long |
|
hash(long q) { return (q & ((1 << (HASHBITS+1)) - 1)) >> 1; } |
|
#undef HASHBITS |
|
|
/* See buch2.c: |
/* See buch2.c: |
* precision en digits decimaux=2*(#digits decimaux de Disc)+50 |
* precision en digits decimaux=2*(#digits decimaux de Disc)+50 |
* on prendra les p decomposes tels que prod(p)>lim dans la subbase |
* on prendra les p decomposes tels que prod(p)>lim dans la subFB |
* LIMC=Max(c.(log(Disc))^2,exp((1/8).sqrt(log(Disc).loglog(Disc)))) |
* LIMC=Max(c.(log(Disc))^2,exp((1/8).sqrt(log(Disc).loglog(Disc)))) |
* LIMC2=Max(6.(log(Disc))^2,exp((1/8).sqrt(log(Disc).loglog(Disc)))) |
* LIMC2=Max(6.(log(Disc))^2,exp((1/8).sqrt(log(Disc).loglog(Disc)))) |
* subbase contient les p decomposes tels que prod(p)>sqrt(Disc) |
* subFB contient les indices des p decomposes tels que prod(p) > sqrt(Disc) |
* lgsub=subbase[0]=#subbase; |
* powsubFB est la table des puissances des formes dans subFB |
* subfactorbase est la table des form[p] pour p dans subbase |
|
* nbram est le nombre de p divisant Disc elimines dans subbase |
|
* powsubfactorbase est la table des puissances des formes dans subfactorbase |
|
*/ |
*/ |
#define HASHT 1024 /* power of 2 */ |
|
static const long CBUCH = 15; /* of the form 2^k-1 */ |
|
static const long randshift=BITS_IN_RANDOM-1 - 4; /* BITS_IN_RANDOM-1-k */ |
|
|
|
static long KC,KC2,lgsub,limhash,RELSUP,PRECREG; |
#define RANDOM_BITS 4 |
static long *primfact,*primfact1, *exprimfact,*exprimfact1, *badprim; |
static const long CBUCH = (1<<RANDOM_BITS)-1; |
static long *factorbase,*numfactorbase, *subbase, *vectbase, **hashtab; |
static const long randshift=BITS_IN_RANDOM-1 - RANDOM_BITS; |
static GEN **powsubfactorbase,vperm,subfactorbase,Disc,sqrtD,isqrtD; |
#undef RANDOM_BITS |
|
|
|
static long KC,KC2,limhash,RELSUP,PRECREG; |
|
static long *primfact,*exprimfact,*badprim; |
|
static long *FB,*numFB, **hashtab; |
|
static GEN powsubFB,vperm,subFB,Disc,sqrtD,isqrtD; |
|
|
GEN buchquad(GEN D, double c, double c2, long RELSUP0, long flag, long prec); |
GEN buchquad(GEN D, double c, double c2, long RELSUP0, long flag, long prec); |
extern GEN roots_to_pol_intern(GEN L, GEN a, long v, int plus); |
extern GEN roots_to_pol_intern(GEN L, GEN a, long v, int plus); |
extern GEN roots_to_pol(GEN L, long v); |
extern GEN colreducemodHNF(GEN x, GEN y, GEN *Q); |
extern GEN colreducemodmat(GEN x, GEN y, GEN *Q); |
extern GEN quadhilbertreal(GEN D, GEN flag, long prec); |
|
extern void comp_gen(GEN z,GEN x,GEN y); |
|
extern GEN codeform5(GEN x, long prec); |
|
extern GEN comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqrtD); |
|
extern GEN redrealform5(GEN x, GEN D, GEN sqrtD, GEN isqrtD); |
|
extern GEN rhoreal_aux(GEN x, GEN D, GEN sqrtD, GEN isqrtD); |
|
extern GEN cgetalloc(GEN x, size_t l, long t); |
|
|
GEN |
GEN |
quadclassunit0(GEN x, long flag, GEN data, long prec) |
quadclassunit0(GEN x, long flag, GEN data, long prec) |
Line 109 getallforms(GEN D, long *pth, GEN *ptz) |
|
Line 124 getallforms(GEN D, long *pth, GEN *ptz) |
|
|
|
#define MOD4(x) ((((GEN)x)[2])&3) |
#define MOD4(x) ((((GEN)x)[2])&3) |
#define MAXL 300 |
#define MAXL 300 |
/* find P and Q two non principal prime ideals (above p,q) such that |
/* find P and Q two non principal prime ideals (above p,q) such that |
* (pq, z) = 1 [coprime to all class group representatives] |
* (pq, z) = 1 [coprime to all class group representatives] |
* cl(P) = cl(Q) if P has order 2 in Cl(K) |
* cl(P) = cl(Q) if P has order 2 in Cl(K) |
* Try to have e = 24 / gcd(24, (p-1)(q-1)) as small as possible */ |
* Try to have e = 24 / gcd(24, (p-1)(q-1)) as small as possible */ |
Line 143 get_pq(GEN D, GEN z, GEN flag, GEN *ptp, GEN *ptq) |
|
Line 158 get_pq(GEN D, GEN z, GEN flag, GEN *ptp, GEN *ptq) |
|
ell = 3; |
ell = 3; |
while (l < 3 || ell<=MAXL) |
while (l < 3 || ell<=MAXL) |
{ |
{ |
ell += *diffell++; if (!*diffell) err(primer1); |
NEXT_PRIME_VIADIFF_CHECK(ell, diffell); |
if (smodis(z,ell) && kross(d,ell) > 0) |
if (smodis(z,ell) && kross(d,ell) > 0) |
{ |
{ |
court[2]=ell; form = redimag(primeform(D,court,0)); |
court[2]=ell; form = redimag(primeform(D,court,0)); |
Line 210 gpq(GEN form, GEN p, GEN q, long e, GEN sqd, GEN u, lo |
|
Line 225 gpq(GEN form, GEN p, GEN q, long e, GEN sqd, GEN u, lo |
|
static GEN |
static GEN |
quadhilbertimag(GEN D, GEN flag) |
quadhilbertimag(GEN D, GEN flag) |
{ |
{ |
long av=avma,h,i,e,prec; |
long h, i, e, prec; |
|
gpmem_t av=avma; |
GEN z,L,P,p,q,qfp,qfq,up,uq,u; |
GEN z,L,P,p,q,qfp,qfq,up,uq,u; |
int raw = ((typ(flag)==t_INT && signe(flag))); |
int raw = ((typ(flag)==t_INT && signe(flag))); |
|
|
if (DEBUGLEVEL>=2) timer2(); |
if (DEBUGLEVEL>=2) (void)timer2(); |
if (gcmpgs(D,-11) >= 0) return polx[0]; |
if (gcmpgs(D,-11) >= 0) return polx[0]; |
L = getallforms(D,&h,&z); |
L = getallforms(D,&h,&z); |
if (DEBUGLEVEL>=2) msgtimer("class number = %ld",h); |
if (DEBUGLEVEL>=2) msgtimer("class number = %ld",h); |
Line 237 quadhilbertimag(GEN D, GEN flag) |
|
Line 253 quadhilbertimag(GEN D, GEN flag) |
|
prec = raw? DEFAULTPREC: 3; |
prec = raw? DEFAULTPREC: 3; |
for(;;) |
for(;;) |
{ |
{ |
long av0 = avma, ex, exmax = 0; |
long ex, exmax = 0; |
|
gpmem_t av0 = avma; |
GEN lead, sqd = gsqrt(negi(D),prec); |
GEN lead, sqd = gsqrt(negi(D),prec); |
P = cgetg(h+1,t_VEC); |
P = cgetg(h+1,t_VEC); |
for (i=1; i<=h; i++) |
for (i=1; i<=h; i++) |
{ |
{ |
GEN v, s = gpq((GEN)L[i], p, q, e, sqd, u, prec); |
GEN v, s = gpq((GEN)L[i], p, q, e, sqd, u, prec); |
Line 250 quadhilbertimag(GEN D, GEN flag) |
|
Line 267 quadhilbertimag(GEN D, GEN flag) |
|
v[2] = (long)s; |
v[2] = (long)s; |
} |
} |
else P[i] = (long)s; |
else P[i] = (long)s; |
ex = gexpo(s); if (ex > 0) exmax += ex; |
ex = gexpo(s); if (ex > 0) exmax += ex; |
} |
} |
if (DEBUGLEVEL>=2) msgtimer("roots"); |
if (DEBUGLEVEL>=2) msgtimer("roots"); |
if (raw) { P = gcopy(P); break; } |
if (raw) { P = gcopy(P); break; } |
Line 271 quadhilbertimag(GEN D, GEN flag) |
|
Line 288 quadhilbertimag(GEN D, GEN flag) |
|
return gerepileupto(av,P); |
return gerepileupto(av,P); |
} |
} |
|
|
GEN quadhilbertreal(GEN D, GEN flag, long prec); |
|
|
|
GEN |
GEN |
quadhilbert(GEN D, GEN flag, long prec) |
quadhilbert(GEN D, GEN flag, long prec) |
{ |
{ |
Line 310 get_om(GEN nf, GEN a) |
|
Line 325 get_om(GEN nf, GEN a) |
|
* Set list[j + 1] = g1^e1...gk^ek where j is the integer |
* Set list[j + 1] = g1^e1...gk^ek where j is the integer |
* ek + ck [ e(k-1) + c(k-1) [... + c2 [e1]]...] */ |
* ek + ck [ e(k-1) + c(k-1) [... + c2 [e1]]...] */ |
static GEN |
static GEN |
getallelts(GEN nf, GEN G) |
getallelts(GEN bnr) |
{ |
{ |
GEN C,c,g, *list, **pows, *gk; |
GEN nf,G,C,c,g, *list, **pows, *gk; |
long lc,i,j,k,no; |
long lc,i,j,k,no; |
|
|
|
nf = checknf(bnr); |
|
G = (GEN)bnr[5]; |
|
|
no = itos((GEN)G[1]); |
no = itos((GEN)G[1]); |
c = (GEN)G[2]; |
c = (GEN)G[2]; |
g = (GEN)G[3]; lc = lg(c)-1; |
g = (GEN)G[3]; lc = lg(c)-1; |
Line 330 getallelts(GEN nf, GEN G) |
|
Line 348 getallelts(GEN nf, GEN G) |
|
{ |
{ |
c[i] = k = itos((GEN)c[i]); |
c[i] = k = itos((GEN)c[i]); |
gk = (GEN*)cgetg(k, t_VEC); gk[1] = (GEN)g[i]; |
gk = (GEN*)cgetg(k, t_VEC); gk[1] = (GEN)g[i]; |
for (j=2; j<k; j++) gk[j] = idealmul(nf, gk[j-1], gk[1]); |
for (j=2; j<k; j++) |
|
gk[j] = idealmodidele(bnr, idealmul(nf, gk[j-1], gk[1])); |
pows[i] = gk; /* powers of g[i] */ |
pows[i] = gk; /* powers of g[i] */ |
} |
} |
|
|
C = cgetg(lc+1, t_VECSMALL); C[1] = c[lc]; |
C = cgetg(lc+1, t_VECSMALL); C[1] = c[lc]; |
for (i=2; i<=lc; i++) C[i] = C[i-1] * c[lc-i+1]; |
for (i=2; i<=lc; i++) C[i] = C[i-1] * c[lc-i+1]; |
/* C[i] = c(k-i+1) * ... * ck */ |
/* C[i] = c(k-i+1) * ... * ck */ |
/* j < C[i+1] <==> j only involves g(k-i)...gk */ |
/* j < C[i+1] <==> j only involves g(k-i)...gk */ |
Line 345 getallelts(GEN nf, GEN G) |
|
Line 364 getallelts(GEN nf, GEN G) |
|
{ |
{ |
GEN p1,p2; |
GEN p1,p2; |
if (j == C[i+1]) i++; |
if (j == C[i+1]) i++; |
p2 = pows[lc-i][j/C[i]]; |
p2 = pows[lc-i][j/C[i]]; |
p1 = list[j%C[i] + 1]; if (p1) p2 = idealmul(nf,p2,p1); |
p1 = list[j%C[i] + 1]; |
|
if (p1) p2 = idealmodidele(bnr, idealmul(nf,p2,p1)); |
list[j + 1] = p2; |
list[j + 1] = p2; |
} |
} |
list[1] = idealhermite(nf,gun); |
list[1] = idealhermite(nf,gun); |
Line 419 get_prec(GEN P, long prec) |
|
Line 439 get_prec(GEN P, long prec) |
|
} |
} |
|
|
/* AUXILLIARY ROUTINES FOR QUADRAYSIGMA */ |
/* AUXILLIARY ROUTINES FOR QUADRAYSIGMA */ |
GEN |
|
PiI2(long prec) |
|
{ |
|
GEN z = cgetg(3,t_COMPLEX); |
|
GEN x = mppi(prec); setexpo(x,2); |
|
z[1] = zero; |
|
z[2] = (long)x; return z; /* = 2I Pi */ |
|
} |
|
|
|
/* Compute data for ellphist */ |
/* Compute data for ellphist */ |
static GEN |
static GEN |
Line 436 ellphistinit(GEN om, long prec) |
|
Line 448 ellphistinit(GEN om, long prec) |
|
|
|
if (gsigne(gimag(gdiv(om1,om2))) < 0) |
if (gsigne(gimag(gdiv(om1,om2))) < 0) |
{ |
{ |
p1=om1; om1=om2; om2=p1; |
p1 = om1; om1 = om2; om2 = p1; |
p1=cgetg(3,t_VEC); p1[1]=(long)om1; p1[2]=(long)om2; |
om = cgetg(3,t_VEC); |
om = p1; |
om[1] = (long)om1; |
|
om[2] = (long)om2; |
} |
} |
om1b = gconj(om1); |
om1b = gconj(om1); |
om2b = gconj(om2); res = cgetg(4,t_VEC); |
om2b = gconj(om2); res = cgetg(4,t_VEC); |
Line 475 computeth2(GEN om, GEN la, long prec) |
|
Line 488 computeth2(GEN om, GEN la, long prec) |
|
static GEN |
static GEN |
computeP2(GEN bnr, GEN la, int raw, long prec) |
computeP2(GEN bnr, GEN la, int raw, long prec) |
{ |
{ |
long av=avma, av2, clrayno,i, first = 1; |
long clrayno, i, first = 1; |
|
gpmem_t av=avma, av2; |
GEN listray,P0,P,f,lanum, nf = checknf(bnr); |
GEN listray,P0,P,f,lanum, nf = checknf(bnr); |
|
|
f = gmael3(bnr,2,1,1); |
f = gmael3(bnr,2,1,1); |
if (typ(la) != t_COL) la = algtobasis(nf,la); |
if (typ(la) != t_COL) la = algtobasis(nf,la); |
listray = getallelts(nf,(GEN)bnr[5]); |
listray = getallelts(bnr); |
clrayno = lg(listray)-1; av2 = avma; |
clrayno = lg(listray)-1; av2 = avma; |
PRECPB: |
PRECPB: |
if (!first) |
if (!first) |
Line 518 do_compo(GEN x, GEN y) |
|
Line 532 do_compo(GEN x, GEN y) |
|
{ |
{ |
long a, ph = degpol(y); |
long a, ph = degpol(y); |
GEN z; |
GEN z; |
y = gmul(gpuigs(polx[0],ph), gsubst(y,0,gdiv(polx[MAXVARN],polx[0]))); |
y = gmul(gpowgs(polx[0],ph), gsubst(y,0,gdiv(polx[MAXVARN],polx[0]))); |
for (a = 0;; a = nexta(a)) |
for (a = 0;; a = nexta(a)) |
{ |
{ |
if (a) x = gsubst(x, 0, gaddsg(a, polx[0])); |
if (a) x = gsubst(x, 0, gaddsg(a, polx[0])); |
Line 543 galoisapplypol(GEN nf, GEN s, GEN x) |
|
Line 557 galoisapplypol(GEN nf, GEN s, GEN x) |
|
static GEN |
static GEN |
findquad(GEN a, GEN x, GEN p) |
findquad(GEN a, GEN x, GEN p) |
{ |
{ |
long tu,tv, av = avma; |
long tu, tv; |
|
gpmem_t av = avma; |
GEN u,v; |
GEN u,v; |
if (typ(x) == t_POLMOD) x = (GEN)x[2]; |
if (typ(x) == t_POLMOD) x = (GEN)x[2]; |
if (typ(a) == t_POLMOD) a = (GEN)a[2]; |
if (typ(a) == t_POLMOD) a = (GEN)a[2]; |
Line 645 treatspecialsigma(GEN nf, GEN gf, int raw, long prec) |
|
Line 660 treatspecialsigma(GEN nf, GEN gf, int raw, long prec) |
|
return NULL; |
return NULL; |
} |
} |
|
|
p1 = gcoeff(gf,1,1); |
p1 = gcoeff(gf,1,1); /* integer > 0 */ |
p2 = gcoeff(gf,2,2); |
p2 = gcoeff(gf,2,2); |
if (gcmp1(p2)) { fl = 0; tryf = p1; } |
if (gcmp1(p2)) { fl = 0; tryf = p1; } |
else { |
else { |
if (Ds % 16 != 8 || !egalii(content(gf),gdeux)) return NULL; |
if (Ds % 16 != 8 || !egalii(content(gf),gdeux)) return NULL; |
fl = 1; tryf = shifti(p1,-1); |
fl = 1; tryf = shifti(p1,-1); |
} |
} |
if (cmpis(tryf, 3) <= 0 || !gcmp0(resii(D, tryf)) || !isprime(tryf)) |
/* tryf integer > 0 */ |
|
if (cmpis(tryf, 3) <= 0 || signe(resii(D, tryf)) || !isprime(tryf)) |
return NULL; |
return NULL; |
|
|
i = itos(tryf); if (fl) i *= 4; |
i = itos(tryf); if (fl) i *= 4; |
Line 701 quadrayimagsigma(GEN bnr, int raw, long prec) |
|
Line 717 quadrayimagsigma(GEN bnr, int raw, long prec) |
|
f2 = 2 * itos(gcoeff(f,1,1)); |
f2 = 2 * itos(gcoeff(f,1,1)); |
u = getallrootsof1(bnf); lu = lg(u); |
u = getallrootsof1(bnf); lu = lg(u); |
for (i=1; i<lu; i++) |
for (i=1; i<lu; i++) |
u[i] = (long)colreducemodmat((GEN)u[i], f, NULL); /* roots of 1, mod f */ |
u[i] = (long)colreducemodHNF((GEN)u[i], f, NULL); /* roots of 1, mod f */ |
if (DEBUGLEVEL>1) |
if (DEBUGLEVEL>1) |
fprintferr("quadray: looking for [a,b] != unit mod 2f\n[a,b] = "); |
fprintferr("quadray: looking for [a,b] != unit mod 2f\n[a,b] = "); |
for (a=0; a<f2; a++) |
for (a=0; a<f2; a++) |
Line 712 quadrayimagsigma(GEN bnr, int raw, long prec) |
|
Line 728 quadrayimagsigma(GEN bnr, int raw, long prec) |
|
if (DEBUGLEVEL>1) fprintferr("[%ld,%ld] ",a,b); |
if (DEBUGLEVEL>1) fprintferr("[%ld,%ld] ",a,b); |
|
|
labas = algtobasis(nf, la); |
labas = algtobasis(nf, la); |
lamodf = colreducemodmat(labas, f, NULL); |
lamodf = colreducemodHNF(labas, f, NULL); |
for (i=1; i<lu; i++) |
for (i=1; i<lu; i++) |
if (gegal(lamodf, (GEN)u[i])) break; |
if (gegal(lamodf, (GEN)u[i])) break; |
if (i < lu) continue; /* la = unit mod f */ |
if (i < lu) continue; /* la = unit mod f */ |
|
|
quadray(GEN D, GEN f, GEN flag, long prec) |
quadray(GEN D, GEN f, GEN flag, long prec) |
{ |
{ |
GEN bnr,y,p1,pol,bnf,lambda; |
GEN bnr,y,p1,pol,bnf,lambda; |
long av = avma, raw; |
long raw; |
|
gpmem_t av = avma; |
|
|
if (!flag) flag = gzero; |
if (!flag) flag = gzero; |
if (typ(flag)==t_INT) lambda = NULL; |
if (typ(flag)==t_INT) lambda = NULL; |
Line 781 quadray(GEN D, GEN f, GEN flag, long prec) |
|
Line 798 quadray(GEN D, GEN f, GEN flag, long prec) |
|
/* Routines related to binary quadratic forms (for internal use) */ |
/* Routines related to binary quadratic forms (for internal use) */ |
/* */ |
/* */ |
/*******************************************************************/ |
/*******************************************************************/ |
extern void comp_gen(GEN z,GEN x,GEN y); |
|
extern GEN codeform5(GEN x, long prec); |
|
extern GEN comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqrtD); |
|
extern GEN redrealform5(GEN x, GEN D, GEN sqrtD, GEN isqrtD); |
|
extern GEN rhoreal_aux(GEN x, GEN D, GEN sqrtD, GEN isqrtD); |
|
|
|
#define rhorealform(x) rhoreal_aux(x,Disc,sqrtD,isqrtD) |
#define rhorealform(x) rhoreal_aux(x,Disc,sqrtD,isqrtD) |
#define redrealform(x) gcopy(fix_signs(redrealform5(x,Disc,sqrtD,isqrtD))) |
#define redrealform(x) gcopy(fix_signs(redrealform5(x,Disc,sqrtD,isqrtD))) |
#define comprealform(x,y) fix_signs(comprealform5(x,y,Disc,sqrtD,isqrtD)) |
|
|
|
|
static GEN |
|
fix_signs(GEN x) |
|
{ |
|
GEN a = (GEN)x[1]; |
|
GEN c = (GEN)x[3]; |
|
if (signe(a) < 0) |
|
{ |
|
if (narrow || absi_equal(a,c)) return rhorealform(x); |
|
setsigne(a,1); setsigne(c,-1); |
|
} |
|
return x; |
|
} |
|
|
|
static GEN |
|
comprealform(GEN x,GEN y) { |
|
return fix_signs( comprealform5(x,y,Disc,sqrtD,isqrtD) ); |
|
} |
|
|
/* compute rho^n(x) */ |
/* compute rho^n(x) */ |
static GEN |
static GEN |
rhoreal_pow(GEN x, long n) |
rhoreal_pow(GEN x, long n) |
{ |
{ |
long i, av = avma, lim = stack_lim(av,1); |
long i; |
|
gpmem_t av = avma, lim = stack_lim(av, 1); |
for (i=1; i<=n; i++) |
for (i=1; i<=n; i++) |
{ |
{ |
x = rhorealform(x); |
x = rhorealform(x); |
if (low_stack(lim, stack_lim(av,1))) |
if (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
if(DEBUGMEM>1) err(warnmem,"rhoreal_pow"); |
if(DEBUGMEM>1) err(warnmem,"rhoreal_pow"); |
x = gerepilecopy(av, x); |
x = gerepilecopy(av, x); |
} |
} |
} |
} |
Line 809 rhoreal_pow(GEN x, long n) |
|
Line 838 rhoreal_pow(GEN x, long n) |
|
} |
} |
|
|
static GEN |
static GEN |
fix_signs(GEN x) |
realpf5(GEN D, long p) |
{ |
{ |
GEN a = (GEN)x[1]; |
gpmem_t av = avma; |
GEN c = (GEN)x[3]; |
GEN y = primeform(D,stoi(p),PRECREG); |
if (signe(a) < 0) |
|
{ |
|
if (narrow || absi_equal(a,c)) return rhorealform(x); |
|
setsigne(a,1); setsigne(c,-1); |
|
} |
|
return x; |
|
} |
|
|
|
static GEN |
|
redrealprimeform5(GEN Disc, long p) |
|
{ |
|
long av = avma; |
|
GEN y = primeform(Disc,stoi(p),PRECREG); |
|
y = codeform5(y,PRECREG); |
y = codeform5(y,PRECREG); |
return gerepileupto(av, redrealform(y)); |
return gerepileupto(av, redrealform(y)); |
} |
} |
|
|
static GEN |
static GEN |
redrealprimeform(GEN Disc, long p) |
realpf(GEN D, long p) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
GEN y = primeform(Disc,stoi(p),PRECREG); |
GEN y = primeform(D,stoi(p),PRECREG); |
return gerepileupto(av, redrealform(y)); |
return gerepileupto(av, redrealform(y)); |
} |
} |
|
|
static GEN |
static GEN |
|
imagpf(GEN D, long p) { return primeform(D,stoi(p),0); } |
|
|
|
static GEN |
comprealform3(GEN x, GEN y) |
comprealform3(GEN x, GEN y) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
GEN z = cgetg(4,t_VEC); comp_gen(z,x,y); |
GEN z = cgetg(4,t_VEC); comp_gen(z,x,y); |
return gerepileupto(av, redrealform(z)); |
return gerepileupto(av, redrealform(z)); |
} |
} |
|
|
|
/* Warning: ex[0] not set */ |
static GEN |
static GEN |
initrealform5(long *ex) |
init_form(long *ex, GEN (*comp)(GEN,GEN)) |
{ |
{ |
GEN form = powsubfactorbase[1][ex[1]]; |
long i, l = lg(powsubFB); |
long i; |
GEN F = NULL; |
|
for (i=1; i<l; i++) |
|
if (ex[i]) |
|
{ |
|
GEN t = gmael(powsubFB,i,ex[i]); |
|
F = F? comp(F,t): t; |
|
} |
|
return F; |
|
} |
|
static GEN |
|
initrealform5(long *ex) { return init_form(ex, &comprealform); } |
|
static GEN |
|
initimagform(long *ex) { return init_form(ex, &compimag); } |
|
|
for (i=2; i<=lgsub; i++) |
static GEN |
form = comprealform(form, powsubfactorbase[i][ex[i]]); |
random_form(GEN ex, GEN (*comp)(GEN,GEN)) |
return form; |
{ |
|
long i, l = lg(ex); |
|
gpmem_t av = avma; |
|
GEN F; |
|
for(;;) |
|
{ |
|
for (i=1; i<l; i++) ex[i] = mymyrand()>>randshift; |
|
if ((F = init_form(ex, comp))) return F; |
|
avma = av; |
|
} |
} |
} |
|
|
|
static GEN |
|
real_random_form(GEN ex) { return random_form(ex, &comprealform3); } |
|
static GEN |
|
imag_random_form(GEN ex) { return random_form(ex, &compimag); } |
|
|
/*******************************************************************/ |
/*******************************************************************/ |
/* */ |
/* */ |
/* Common subroutines */ |
/* Common subroutines */ |
Line 865 initrealform5(long *ex) |
|
Line 911 initrealform5(long *ex) |
|
static void |
static void |
buch_init(void) |
buch_init(void) |
{ |
{ |
if (DEBUGLEVEL) timer2(); |
if (DEBUGLEVEL) (void)timer2(); |
primfact = new_chunk(100); |
primfact = new_chunk(100); |
primfact1 = new_chunk(100); |
|
exprimfact = new_chunk(100); |
exprimfact = new_chunk(100); |
exprimfact1 = new_chunk(100); |
|
badprim = new_chunk(100); |
badprim = new_chunk(100); |
hashtab = (long**) new_chunk(HASHT); |
hashtab = (long**) new_chunk(HASHT); |
} |
} |
Line 885 check_bach(double cbach, double B) |
|
Line 929 check_bach(double cbach, double B) |
|
} |
} |
|
|
static long |
static long |
factorisequad(GEN f, long kcz, long limp) |
factorquad(GEN f, long kcz, long limp) |
{ |
{ |
long i,p,k,av,lo; |
long i, p, k, lo; |
|
gpmem_t av; |
GEN q,r, x = (GEN)f[1]; |
GEN q,r, x = (GEN)f[1]; |
|
|
if (is_pm1(x)) { primfact[0]=0; return 1; } |
if (is_pm1(x)) { primfact[0]=0; return 1; } |
Line 895 factorisequad(GEN f, long kcz, long limp) |
|
Line 940 factorisequad(GEN f, long kcz, long limp) |
|
if (signe(x) < 0) x = absi(x); |
if (signe(x) < 0) x = absi(x); |
for (i=1; ; i++) |
for (i=1; ; i++) |
{ |
{ |
p=factorbase[i]; q=dvmdis(x,p,&r); |
p=FB[i]; q=dvmdis(x,p,&r); |
if (!signe(r)) |
if (!signe(r)) |
{ |
{ |
k=0; while (!signe(r)) { x=q; k++; q=dvmdis(x,p,&r); } |
for (k=0; !signe(r); k++) { x=q; q=dvmdis(x,p,&r); } |
lo++; primfact[lo]=p; exprimfact[lo]=k; |
primfact[++lo]=p; exprimfact[lo]=k; |
} |
} |
if (cmpis(q,p)<=0) break; |
if (cmpis(q,p)<=0) break; |
if (i==kcz) { avma=av; return 0; } |
if (i==kcz) { avma=av; return 0; } |
} |
} |
p = x[2]; avma=av; |
avma=av; |
/* p = itos(x) if lgefint(x)=3 */ |
if (lgefint(x)!=3 || (p=x[2]) > limhash) return 0; |
if (lgefint(x)!=3 || p > limhash) return 0; |
|
|
|
if (p != 1 && p <= limp) |
if (p != 1 && p <= limp) |
{ |
{ |
for (i=1; i<=badprim[0]; i++) |
for (i=1; i<=badprim[0]; i++) |
if (p % badprim[i] == 0) return 0; |
if (p % badprim[i] == 0) return 0; |
lo++; primfact[lo]=p; exprimfact[lo]=1; |
primfact[++lo]=p; exprimfact[lo]=1; |
p = 1; |
p = 1; |
} |
} |
primfact[0]=lo; return p; |
primfact[0]=lo; return p; |
Line 922 factorisequad(GEN f, long kcz, long limp) |
|
Line 966 factorisequad(GEN f, long kcz, long limp) |
|
static long * |
static long * |
largeprime(long q, long *ex, long np, long nrho) |
largeprime(long q, long *ex, long np, long nrho) |
{ |
{ |
const long hashv = ((q & (2 * HASHT - 1)) - 1) >> 1; |
const long hashv = hash(q); |
long *pt, i; |
long *pt, i, l = lg(subFB); |
|
|
/* If q = 0 (2048), very slim chance of getting a relation. |
|
* And hashtab[-1] is undefined anyway */ |
|
if (hashv < 0) return NULL; |
|
for (pt = hashtab[hashv]; ; pt = (long*) pt[0]) |
for (pt = hashtab[hashv]; ; pt = (long*) pt[0]) |
{ |
{ |
if (!pt) |
if (!pt) |
{ |
{ |
pt = (long*) gpmalloc((lgsub+4)<<TWOPOTBYTES_IN_LONG); |
pt = (long*) gpmalloc((l+3)<<TWOPOTBYTES_IN_LONG); |
*pt++ = nrho; /* nrho = pt[-3] */ |
*pt++ = nrho; /* nrho = pt[-3] */ |
*pt++ = np; /* np = pt[-2] */ |
*pt++ = np; /* np = pt[-2] */ |
*pt++ = q; /* q = pt[-1] */ |
*pt++ = q; /* q = pt[-1] */ |
pt[0] = (long)hashtab[hashv]; |
pt[0] = (long)hashtab[hashv]; |
for (i=1; i<=lgsub; i++) pt[i]=ex[i]; |
for (i=1; i<l; i++) pt[i]=ex[i]; |
hashtab[hashv]=pt; return NULL; |
hashtab[hashv]=pt; return NULL; |
} |
} |
if (pt[-1] == q) break; |
if (pt[-1] == q) break; |
} |
} |
for(i=1; i<=lgsub; i++) |
for(i=1; i<l; i++) |
if (pt[i] != ex[i]) return pt; |
if (pt[i] != ex[i]) return pt; |
return (pt[-2]==np)? (GEN)NULL: pt; |
return (pt[-2]==np)? (GEN)NULL: pt; |
} |
} |
|
|
static long |
/* p | conductor of order of disc D ? */ |
badmod8(GEN x) |
static int |
|
is_bad(GEN D, ulong p) |
{ |
{ |
long r = mod8(x); |
gpmem_t av = avma; |
if (!r) return 1; |
int r; |
if (signe(Disc) < 0) r = 8-r; |
if (p == 2) |
return (r < 4); |
{ |
|
r = mod16(D) >> 1; |
|
if (r && signe(D) < 0) r = 8-r; |
|
return (r < 4); |
|
} |
|
r = (resii(D, muluu(p,p)) == gzero); /* p^2 | D ? */ |
|
avma = av; return r; |
} |
} |
|
|
/* cree factorbase, numfactorbase, vectbase; affecte badprim */ |
/* create FB, numFB; fill badprim. Return L(kro_D, 1) */ |
static void |
static GEN |
factorbasequad(GEN Disc, long n2, long n) |
FBquad(GEN Disc, long n2, long n) |
{ |
{ |
long i,p,bad, av = avma; |
GEN Res = realun(DEFAULTPREC); |
byteptr d=diffptr; |
long i, p, bad, s; |
|
gpmem_t av; |
|
byteptr d = diffptr; |
|
|
numfactorbase = (long*) gpmalloc(sizeof(long)*(n2+1)); |
numFB = cgetg(n2+1, t_VECSMALL); |
factorbase = (long*) gpmalloc(sizeof(long)*(n2+1)); |
FB = cgetg(n2+1, t_VECSMALL); |
KC=0; bad=0; i=0; p = *d++; |
av = avma; |
while (p<=n2) |
KC = 0; bad = 0; i = 0; |
|
if (maxprime() < n2) err(primer1); |
|
for (p = 0;;) /* p <= n2 */ |
{ |
{ |
switch(krogs(Disc,p)) |
NEXT_PRIME_VIADIFF(p, d); |
|
if (!KC && p > n) KC = i; |
|
if (p > n2) break; |
|
s = krogs(Disc,p); |
|
switch (s) |
{ |
{ |
case -1: break; /* inert */ |
case -1: break; /* inert */ |
case 0: /* ramified */ |
case 0: /* ramified */ |
{ |
if (is_bad(Disc, (ulong)p)) { badprim[++bad]=p; break; } |
GEN p1 = divis(Disc,p); |
/* fall through */ |
if (smodis(p1,p) == 0) |
|
if (p!=2 || badmod8(p1)) { badprim[++bad]=p; break; } |
|
i++; numfactorbase[p]=i; factorbase[i] = -p; break; |
|
} |
|
default: /* split */ |
default: /* split */ |
i++; numfactorbase[p]=i; factorbase[i] = p; |
i++; numFB[p]=i; FB[i] = p; break; |
} |
} |
p += *d++; if (!*d) err(primer1); |
Res = mulsr(p, divrs(Res, p - s)); |
if (KC == 0 && p>n) KC = i; |
|
} |
} |
if (!KC) { free(factorbase); free(numfactorbase); return; } |
if (!KC) return NULL; |
|
Res = gerepileupto(av, Res); |
KC2 = i; |
KC2 = i; |
vectbase = (long*) gpmalloc(sizeof(long)*(KC2+1)); |
setlg(FB, KC2+1); |
for (i=1; i<=KC2; i++) |
|
{ |
|
p = factorbase[i]; |
|
vectbase[i]=p; factorbase[i]=labs(p); |
|
} |
|
if (DEBUGLEVEL) |
if (DEBUGLEVEL) |
{ |
{ |
msgtimer("factor base"); |
msgtimer("factor base"); |
if (DEBUGLEVEL>7) |
if (DEBUGLEVEL>7) |
{ |
{ |
fprintferr("factorbase:\n"); |
fprintferr("FB:\n"); |
for (i=1; i<=KC; i++) fprintferr("%ld ",factorbase[i]); |
for (i=1; i<=KC; i++) fprintferr("%ld ", FB[i]); |
fprintferr("\n"); flusherr(); |
fprintferr("\n"); |
} |
} |
} |
} |
avma=av; badprim[0] = bad; |
badprim[0] = bad; return Res; |
} |
} |
|
|
/* cree vectbase and subfactorbase. Affecte lgsub */ |
/* create vperm, return subFB */ |
static long |
static GEN |
subfactorbasequad(double ll, long KC) |
subFBquad(GEN D, double PROD, long KC) |
{ |
{ |
long i,j,k,nbidp,p,pro[100], ss; |
long i, j, lgsub, ino, lv = KC+1; |
GEN y; |
double prod = 1.; |
double prod; |
gpmem_t av; |
|
GEN no; |
|
|
i=0; ss=0; prod=1; |
vperm = cgetg(lv, t_VECSMALL); |
for (j=1; j<=KC && prod<=ll; j++) |
av = avma; |
|
no = cgetg(lv, t_VECSMALL); ino = 1; |
|
for (i=j=1; j < lv; j++) |
{ |
{ |
p = vectbase[j]; |
long p = FB[j]; |
if (p>0) { pro[++i]=p; prod*=p; vperm[i]=j; } else ss++; |
if (smodis(D, p) == 0) no[ino++] = j; /* ramified */ |
|
else |
|
{ |
|
vperm[i] = j; i++; |
|
prod *= p; |
|
if (prod > PROD) break; |
|
} |
} |
} |
if (prod<=ll) return -1; |
if (j == lv) return NULL; |
nbidp=i; |
lgsub = i; |
for (k=1; k<j; k++) |
for (j = 1; j <ino; i++,j++) vperm[i] = no[j]; |
if (vectbase[k]<=0) vperm[++i]=k; |
for ( ; i < lv; i++) vperm[i] = i; |
|
avma = av; |
y=cgetg(nbidp+1,t_COL); |
if (DEBUGLEVEL) msgtimer("subFBquad (%ld elt.)", lgsub-1); |
if (PRECREG) /* real */ |
return vecextract_i(vperm, 1, lgsub-1); |
for (j=1; j<=nbidp; j++) |
|
y[j] = (long)redrealprimeform5(Disc, pro[j]); |
|
else |
|
for (j=1; j<=nbidp; j++) /* imaginary */ |
|
y[j] = (long)primeform(Disc,stoi(pro[j]),0); |
|
subbase = (long*) gpmalloc(sizeof(long)*(nbidp+1)); |
|
lgsub = nbidp; for (j=1; j<=lgsub; j++) subbase[j]=pro[j]; |
|
if (DEBUGLEVEL>7) |
|
{ |
|
fprintferr("subfactorbase: "); |
|
for (i=1; i<=lgsub; i++) |
|
{ fprintferr("%ld: ",subbase[i]); outerr((GEN)y[i]); } |
|
fprintferr("\n"); flusherr(); |
|
} |
|
subfactorbase = y; return ss; |
|
} |
} |
|
|
static void |
/* assume n >= 1, x[i][j] = subFB[i]^j, for j = 1..n */ |
powsubfact(long n, long a) |
static GEN |
|
powsubFBquad(long n) |
{ |
{ |
GEN unform, *y, **x = (GEN**) gpmalloc(sizeof(GEN*)*(n+1)); |
long i,j, l = lg(subFB); |
long i,j; |
GEN F, y, x = cgetg(l, t_VEC); |
|
|
for (i=1; i<=n; i++) |
|
x[i] = (GEN*) gpmalloc(sizeof(GEN)*(a+1)); |
|
if (PRECREG) /* real */ |
if (PRECREG) /* real */ |
{ |
{ |
unform=cgetg(6,t_VEC); |
for (i=1; i<l; i++) |
unform[1]=un; |
|
unform[2]=(mod2(Disc) == mod2(isqrtD))? (long)isqrtD: laddsi(-1,isqrtD); |
|
unform[3]=lshifti(subii(sqri((GEN)unform[2]),Disc),-2); |
|
unform[4]=zero; |
|
unform[5]=(long)realun(PRECREG); |
|
for (i=1; i<=n; i++) |
|
{ |
{ |
y = x[i]; y[0] = unform; |
F = realpf5(Disc, FB[subFB[i]]); |
for (j=1; j<=a; j++) |
y = cgetg(n+1, t_VEC); x[i] = (long)y; |
y[j] = comprealform(y[j-1],(GEN)subfactorbase[i]); |
y[1] = (long)F; |
|
for (j=2; j<=n; j++) y[j] = (long)comprealform((GEN)y[j-1], F); |
} |
} |
} |
} |
else /* imaginary */ |
else /* imaginary */ |
{ |
{ |
unform = primeform(Disc, gun, 0); |
for (i=1; i<l; i++) |
for (i=1; i<=n; i++) |
|
{ |
{ |
y = x[i]; y[0] = unform; |
F = imagpf(Disc, FB[subFB[i]]); |
for (j=1; j<=a; j++) |
y = cgetg(n+1, t_VEC); x[i] = (long)y; |
y[j] = compimag(y[j-1],(GEN)subfactorbase[i]); |
y[1] = (long)F; |
|
for (j=2; j<=n; j++) y[j] = (long)compimag((GEN)y[j-1], F); |
} |
} |
} |
} |
if (DEBUGLEVEL) msgtimer("powsubfact"); |
if (DEBUGLEVEL) msgtimer("powsubFBquad"); |
powsubfactorbase = x; |
return x; |
} |
} |
|
|
static void |
static void |
desalloc(long **mat) |
desalloc(long **mat) |
{ |
{ |
long i,*p,*q; |
long i,*p,*q; |
|
for (i=1; i<lg(mat); i++) free(mat[i]); |
|
free(mat); |
|
for (i=1; i<HASHT; i++) |
|
for (p = hashtab[i]; p; p = q) { q=(long*)p[0]; free(p-3); } |
|
} |
|
|
free(vectbase); free(factorbase); free(numfactorbase); |
static void |
if (mat) |
sub_fact(GEN col, GEN F) |
|
{ |
|
GEN b = (GEN)F[2]; |
|
long i, p, e; |
|
for (i=1; i<=primfact[0]; i++) |
{ |
{ |
free(subbase); |
p = primfact[i]; e = exprimfact[i]; |
for (i=1; i<lg(subfactorbase); i++) free(powsubfactorbase[i]); |
if (smodis(b, p<<1) > p) e = -e; |
for (i=1; i<lg(mat); i++) free(mat[i]); |
col[numFB[p]] -= e; |
free(mat); free(powsubfactorbase); |
|
for (i=1; i<HASHT; i++) |
|
for (p = hashtab[i]; p; p = q) { q=(long*)p[0]; free(p-3); } |
|
} |
} |
} |
} |
|
static void |
/* L-function */ |
add_fact(GEN col, GEN F) |
static GEN |
|
lfunc(GEN Disc) |
|
{ |
{ |
long av=avma, p; |
GEN b = (GEN)F[2]; |
GEN y=realun(DEFAULTPREC); |
long i, p, e; |
byteptr d=diffptr; |
for (i=1; i<=primfact[0]; i++) |
|
|
for(p = *d++; p<=30000; p += *d++) |
|
{ |
{ |
if (!*d) err(primer1); |
p = primfact[i]; e = exprimfact[i]; |
y = mulsr(p, divrs(y, p-krogs(Disc,p))); |
if (smodis(b, p<<1) > p) e = -e; |
|
col[numFB[p]] += e; |
} |
} |
return gerepileupto(av,y); |
|
} |
} |
|
|
#define comp(x,y) x? (PRECREG? compreal(x,y): compimag(x,y)): y |
#define comp(x,y) x? (PRECREG? compreal(x,y): compimag(x,y)): y |
static GEN |
static GEN |
get_clgp(GEN Disc, GEN W, GEN *ptmet, long prec) |
get_clgp(GEN Disc, GEN W, GEN *ptmet, long prec) |
{ |
{ |
GEN p1,p2,res,*init, u1u2=smith2(W), u1=(GEN)u1u2[1], met=(GEN)u1u2[3]; |
GEN res,*init, u1, met = smithrel(W,NULL,&u1); |
long c,i,j, l = lg(met); |
long c,i,j, l = lg(W); |
|
|
u1 = reducemodmatrix(ginv(u1), W); |
c = lg(met); |
for (c=1; c<l; c++) |
|
if (gcmp1(gcoeff(met,c,c))) break; |
|
if (DEBUGLEVEL) msgtimer("smith/class group"); |
if (DEBUGLEVEL) msgtimer("smith/class group"); |
res=cgetg(c,t_VEC); init = (GEN*)cgetg(l,t_VEC); |
res=cgetg(c,t_VEC); init = (GEN*)cgetg(l,t_VEC); |
for (i=1; i<l; i++) |
for (i=1; i<l; i++) |
init[i] = primeform(Disc,stoi(labs(vectbase[vperm[i]])),prec); |
init[i] = primeform(Disc,stoi(FB[vperm[i]]),prec); |
for (j=1; j<c; j++) |
for (j=1; j<c; j++) |
{ |
{ |
p1 = NULL; |
GEN p1 = NULL; |
for (i=1; i<l; i++) |
for (i=1; i<l; i++) |
{ |
p1 = comp(p1, powgi(init[i], gcoeff(u1,i,j))); |
p2 = gpui(init[i], gcoeff(u1,i,j), prec); |
|
p1 = comp(p1,p2); |
|
} |
|
res[j] = (long)p1; |
res[j] = (long)p1; |
} |
} |
if (DEBUGLEVEL) msgtimer("generators"); |
if (DEBUGLEVEL) msgtimer("generators"); |
Line 1142 get_clgp(GEN Disc, GEN W, GEN *ptmet, long prec) |
|
Line 1174 get_clgp(GEN Disc, GEN W, GEN *ptmet, long prec) |
|
} |
} |
|
|
static GEN |
static GEN |
extra_relations(long LIMC, long *ex, long nlze, GEN extramatc) |
extra_relations(long LIMC, long nlze, GEN *ptextraC) |
{ |
{ |
long av,fpc,p,ep,i,j,k,nlze2, *col, *colg, s = 0, extrarel = nlze+2; |
long fpc, i, k, nlze2; |
|
long s = 0, extrarel = nlze+2, lgsub = lg(subFB); |
|
gpmem_t av; |
long MAXRELSUP = min(50,4*KC); |
long MAXRELSUP = min(50,4*KC); |
GEN p1,form, extramat = cgetg(extrarel+1,t_MAT); |
GEN p1, form, ex, extramat, extraC, col; |
|
|
if (DEBUGLEVEL) |
extramat = cgetg(extrarel+1, t_MAT); |
{ |
extraC = cgetg(extrarel+1, t_VEC); |
fprintferr("looking for %ld extra relations\n",extrarel); |
for (i=1; i<=extrarel; i++) extramat[i] = lgetg(KC+1,t_VECSMALL); |
flusherr(); |
if (!PRECREG) |
} |
for (i=1; i<=extrarel; i++) extraC[i] = lgetg(1,t_COL); |
for (j=1; j<=extrarel; j++) extramat[j]=lgetg(KC+1,t_COL); |
|
nlze2 = PRECREG? max(nlze,lgsub): min(nlze+1,KC); |
if (DEBUGLEVEL) fprintferr("looking for %ld extra relations\n",extrarel); |
|
nlze2 = PRECREG? max(nlze,lgsub-1): min(nlze+1,KC); |
if (nlze2 < 3 && KC > 2) nlze2 = 3; |
if (nlze2 < 3 && KC > 2) nlze2 = 3; |
|
ex = cgetg(nlze2+1, t_VECSMALL); |
av = avma; |
av = avma; |
while (s<extrarel) |
while (s<extrarel) |
{ |
{ |
form = NULL; |
form = NULL; |
for (i=1; i<=nlze2; i++) |
for (i=1; i<=nlze2; i++) |
{ |
{ |
ex[i]=mymyrand()>>randshift; |
ex[i] = mymyrand()>>randshift; |
if (ex[i]) |
if (ex[i]) |
{ |
{ |
p1 = primeform(Disc,stoi(factorbase[vperm[i]]),PRECREG); |
p1 = primeform(Disc,stoi(FB[vperm[i]]),PRECREG); |
p1 = gpuigs(p1,ex[i]); form = comp(form,p1); |
p1 = gpowgs(p1,ex[i]); form = comp(form,p1); |
} |
} |
} |
} |
if (!form) continue; |
if (!form) continue; |
|
|
fpc = factorisequad(form,KC,LIMC); |
fpc = factorquad(form,KC,LIMC); |
if (fpc==1) |
if (fpc==1) |
{ |
{ |
s++; col = (GEN)extramat[s]; |
s++; col = (GEN)extramat[s]; |
for (i=1; i<=nlze2; i++) col[vperm[i]] = -ex[i]; |
for (i=1; i<=nlze2; i++) col[vperm[i]] = -ex[i]; |
for ( ; i<=KC; i++) col[vperm[i]]= 0; |
for ( ; i<=KC; i++) col[vperm[i]] = 0; |
for (j=1; j<=primfact[0]; j++) |
add_fact(col, form); |
{ |
|
p=primfact[j]; ep=exprimfact[j]; |
|
if (smodis((GEN)form[2], p<<1) > p) ep = -ep; |
|
col[numfactorbase[p]] += ep; |
|
} |
|
for (i=1; i<=KC; i++) |
for (i=1; i<=KC; i++) |
if (col[i]) break; |
if (col[i]) break; |
if (i>KC) |
if (i>KC) |
Line 1190 extra_relations(long LIMC, long *ex, long nlze, GEN ex |
|
Line 1221 extra_relations(long LIMC, long *ex, long nlze, GEN ex |
|
s--; avma = av; |
s--; avma = av; |
if (--MAXRELSUP == 0) return NULL; |
if (--MAXRELSUP == 0) return NULL; |
} |
} |
else { av = avma; if (PRECREG) coeff(extramatc,1,s) = form[4]; } |
else { av = avma; if (PRECREG) extraC[s] = form[4]; } |
} |
} |
else avma = av; |
else avma = av; |
if (DEBUGLEVEL) |
if (DEBUGLEVEL) |
{ |
{ |
if (fpc == 1) fprintferr(" %ld",s); |
if (fpc == 1) fprintferr(" %ld",s); |
else if (DEBUGLEVEL>1) fprintferr("."); |
else if (DEBUGLEVEL>1) fprintferr("."); |
flusherr(); |
|
} |
} |
} |
} |
for (j=1; j<=extrarel; j++) |
|
|
for (i=1; i<=extrarel; i++) |
{ |
{ |
colg = cgetg(KC+1,t_COL); |
GEN colg = cgetg(KC+1,t_COL); |
col = (GEN)extramat[j]; extramat[j] = (long) colg; |
col = (GEN)extramat[i]; |
for (k=1; k<=KC; k++) |
extramat[i] = (long)colg; |
colg[k] = lstoi(col[vperm[k]]); |
for (k=1; k<=KC; k++) colg[k] = lstoi(col[vperm[k]]); |
} |
} |
if (DEBUGLEVEL) |
if (DEBUGLEVEL) { fprintferr("\n"); msgtimer("extra relations"); } |
|
*ptextraC = extraC; return extramat; |
|
} |
|
#undef comp |
|
|
|
static long |
|
trivial_relations(long **mat, GEN C) |
|
{ |
|
long i, s = 0; |
|
GEN col; |
|
for (i=1; i<=KC; i++) |
{ |
{ |
fprintferr("\n"); |
if (smodis(Disc, FB[i])) continue; |
msgtimer("extra relations"); |
/* ramified prime ==> trivial relation */ |
|
col = mat[++s]; |
|
col[i] = 2; |
|
if (C) affsr(0, (GEN)C[s]); |
} |
} |
return extramat; |
return s; |
} |
} |
#undef comp |
|
|
|
|
static void |
|
dbg_all(long phase, long s, long n) |
|
{ |
|
if (DEBUGLEVEL>1) fprintferr("\n"); |
|
msgtimer("%s rel [#rel/#test = %ld/%ld]", phase? "random": "initial", s, n); |
|
} |
|
|
|
void |
|
wr_rel(GEN col) |
|
{ |
|
long i, l = lg(col); |
|
fprintferr("\nrel = "); |
|
for (i=1; i<l; i++) |
|
if (col[i]) fprintferr("%ld^%ld ",i,col[i]); |
|
fprintferr("\n"); |
|
} |
|
|
|
void |
|
dbg_rel(long s, GEN col) |
|
{ |
|
if (DEBUGLEVEL == 1) fprintferr("%4ld",s); |
|
else { fprintferr("cglob = %ld. ", s); wr_rel(col); } |
|
flusherr(); |
|
} |
/*******************************************************************/ |
/*******************************************************************/ |
/* */ |
/* */ |
/* Imaginary Quadratic fields */ |
/* Imaginary Quadratic fields */ |
/* */ |
/* */ |
/*******************************************************************/ |
/*******************************************************************/ |
|
|
|
/* Strategy 1 until lim relation found, then Strategy 2 until LIM. */ |
static GEN |
static GEN |
imag_random_form(long current, long *ex) |
imag_relations(long LIM, long lim, long LIMC, long **mat) |
{ |
{ |
long av = avma,i; |
long lgsub = lg(subFB), i, s, fpc, current, nbtest = 0; |
GEN form,pc; |
int first = 1; |
|
gpmem_t av; |
|
GEN C, col, form, ex = cgetg(lgsub, t_VECSMALL); |
|
GEN dummy = cgetg(1,t_COL); |
|
|
|
C = cgetg(LIM+1, t_VEC); |
|
for (i=1; i<=LIM; i++) C[i] = (long)dummy; |
|
av = avma; |
|
s = trivial_relations(mat, NULL); |
|
lim += s; if (lim > LIM) lim = LIM; |
for(;;) |
for(;;) |
{ |
{ |
form = pc = primeform(Disc,stoi(factorbase[current]),PRECREG); |
if (s >= lim) { |
for (i=1; i<=lgsub; i++) |
if (s >= LIM) break; |
{ |
lim = LIM; first = 0; |
ex[i] = mymyrand()>>randshift; |
if (DEBUGLEVEL) dbg_all(0, s, nbtest); |
if (ex[i]) |
|
form = compimag(form,powsubfactorbase[i][ex[i]]); |
|
} |
} |
if (form != pc) return form; |
avma = av; current = first? 1+(s%KC): 1+s-RELSUP; |
avma = av; /* ex = 0, try again */ |
form = imagpf(Disc, FB[current]); |
} |
form = compimag(form, imag_random_form(ex)); |
} |
nbtest++; fpc = factorquad(form,KC,LIMC); |
|
|
static void |
|
imag_relations(long lim, long s, long LIMC, long *ex, long **mat) |
|
{ |
|
static long nbtest; |
|
long av = avma, i,j,pp,fpc,b1,b2,ep,current, first = (s==0); |
|
long *col,*fpd,*oldfact,*oldexp; |
|
GEN pc,form,form1; |
|
|
|
if (first) nbtest = 0 ; |
|
while (s<lim) |
|
{ |
|
avma=av; nbtest++; current = first? 1+(s%KC): 1+s-RELSUP; |
|
form = imag_random_form(current,ex); |
|
fpc = factorisequad(form,KC,LIMC); |
|
if (!fpc) |
if (!fpc) |
{ |
{ |
if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); } |
if (DEBUGLEVEL>1) fprintferr("."); |
continue; |
continue; |
} |
} |
if (fpc > 1) |
if (fpc > 1) |
{ |
{ |
fpd = largeprime(fpc,ex,current,0); |
long *fpd = largeprime(fpc,ex,current,0); |
|
long b1, b2, p; |
|
GEN form2; |
if (!fpd) |
if (!fpd) |
{ |
{ |
if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); } |
if (DEBUGLEVEL>1) fprintferr("."); |
continue; |
continue; |
} |
} |
form1 = powsubfactorbase[1][fpd[1]]; |
form2 = initimagform(fpd); |
for (i=2; i<=lgsub; i++) |
form2 = compimag(form2, imagpf(Disc, FB[fpd[-2]])); |
form1 = compimag(form1,powsubfactorbase[i][fpd[i]]); |
p = fpc << 1; |
pc=primeform(Disc,stoi(factorbase[fpd[-2]]),0); |
b1 = smodis((GEN)form2[2], p); |
form1=compimag(form1,pc); |
b2 = smodis((GEN)form[2], p); |
pp = fpc << 1; |
if (b1 != b2 && b1+b2 != p) continue; |
b1=smodis((GEN)form1[2], pp); |
|
b2=smodis((GEN)form[2], pp); |
|
if (b1 != b2 && b1+b2 != pp) continue; |
|
|
|
s++; col = mat[s]; |
col = mat[++s]; |
if (DEBUGLEVEL) { fprintferr(" %ld",s); flusherr(); } |
add_fact(col, form); |
oldfact = primfact; oldexp = exprimfact; |
(void)factorquad(form2,KC,LIMC); |
primfact = primfact1; exprimfact = exprimfact1; |
|
factorisequad(form1,KC,LIMC); |
|
|
|
if (b1==b2) |
if (b1==b2) |
{ |
{ |
for (i=1; i<=lgsub; i++) |
for (i=1; i<lgsub; i++) col[subFB[i]] += fpd[i]-ex[i]; |
col[numfactorbase[subbase[i]]] = fpd[i]-ex[i]; |
sub_fact(col, form2); col[fpd[-2]]++; |
col[fpd[-2]]++; |
|
for (j=1; j<=primfact[0]; j++) |
|
{ |
|
pp=primfact[j]; ep=exprimfact[j]; |
|
if (smodis((GEN)form1[2], pp<<1) > pp) ep = -ep; |
|
col[numfactorbase[pp]] -= ep; |
|
} |
|
} |
} |
else |
else |
{ |
{ |
for (i=1; i<=lgsub; i++) |
for (i=1; i<lgsub; i++) col[subFB[i]] += -fpd[i]-ex[i]; |
col[numfactorbase[subbase[i]]] = -fpd[i]-ex[i]; |
add_fact(col, form2); col[fpd[-2]]--; |
col[fpd[-2]]--; |
|
for (j=1; j<=primfact[0]; j++) |
|
{ |
|
pp=primfact[j]; ep=exprimfact[j]; |
|
if (smodis((GEN)form1[2], pp<<1) > pp) ep = -ep; |
|
col[numfactorbase[pp]] += ep; |
|
} |
|
} |
} |
primfact = oldfact; exprimfact = oldexp; |
} |
} |
|
else |
else |
{ |
{ |
s++; col = mat[s]; |
col = mat[++s]; |
if (DEBUGLEVEL) { fprintferr(" %ld",s); flusherr(); } |
for (i=1; i<lgsub; i++) col[subFB[i]] = -ex[i]; |
for (i=1; i<=lgsub; i++) |
add_fact(col, form); |
col[numfactorbase[subbase[i]]] = -ex[i]; |
|
} |
} |
for (j=1; j<=primfact[0]; j++) |
|
{ |
|
pp=primfact[j]; ep=exprimfact[j]; |
|
if (smodis((GEN)form[2], pp<<1) > pp) ep = -ep; |
|
col[numfactorbase[pp]] += ep; |
|
} |
|
col[current]--; |
col[current]--; |
|
if (DEBUGLEVEL) dbg_rel(s, col); |
if (!first && fpc == 1 && col[current] == 0) |
if (!first && fpc == 1 && col[current] == 0) |
{ |
{ |
s--; for (i=1; i<=KC; i++) col[i]=0; |
s--; for (i=1; i<=KC; i++) col[i]=0; |
} |
} |
} |
} |
if (DEBUGLEVEL) |
if (DEBUGLEVEL) dbg_all(1, s, nbtest); |
{ |
return C; |
fprintferr("\nnbrelations/nbtest = %ld/%ld\n",s,nbtest); |
|
msgtimer("%s relations", first? "initial": "random"); |
|
} |
|
} |
} |
|
|
static int |
static int |
imag_be_honest(long *ex) |
imag_be_honest() |
{ |
{ |
long p,fpc, s = KC, nbtest = 0, av = avma; |
long p, fpc, s = KC, nbtest = 0; |
GEN form; |
gpmem_t av = avma; |
|
GEN F, ex = cgetg(lg(subFB), t_VECSMALL); |
|
|
while (s<KC2) |
while (s<KC2) |
{ |
{ |
p = factorbase[s+1]; |
p = FB[s+1]; if (DEBUGLEVEL) fprintferr(" %ld",p); |
if (DEBUGLEVEL) { fprintferr(" %ld",p); flusherr(); } |
F = compimag(imagpf(Disc, p), imag_random_form(ex)); |
form = imag_random_form(s+1,ex); |
fpc = factorquad(F,s,p-1); |
fpc = factorisequad(form,s,p-1); |
|
if (fpc == 1) { nbtest=0; s++; } |
if (fpc == 1) { nbtest=0; s++; } |
else |
else |
if (++nbtest > 20) return 0; |
if (++nbtest > 20) return 0; |
Line 1363 imag_be_honest(long *ex) |
|
Line 1394 imag_be_honest(long *ex) |
|
/* */ |
/* */ |
/*******************************************************************/ |
/*******************************************************************/ |
|
|
static GEN |
/* (1/2) log (d * 2^{e * EXP220}) */ |
real_random_form(long *ex) |
GEN |
|
get_dist(GEN e, GEN d, long prec) |
{ |
{ |
long av = avma, i; |
GEN t = mplog(absr(d)); |
GEN p1,form = NULL; |
if (signe(e)) t = addrr(t, mulir(mulsi(EXP220,e), mplog2(prec))); |
|
return shiftr(t, -1); |
for(;;) |
|
{ |
|
for (i=1; i<=lgsub; i++) |
|
{ |
|
ex[i] = mymyrand()>>randshift; |
|
/* if (ex[i]) KB: less efficient if I put this in. Why ??? */ |
|
{ |
|
p1 = powsubfactorbase[i][ex[i]]; |
|
form = form? comprealform3(form,p1): p1; |
|
} |
|
} |
|
if (form) return form; |
|
avma = av; |
|
} |
|
} |
} |
|
|
static void |
static GEN |
real_relations(long lim, long s, long LIMC, long *ex, long **mat, GEN glog2, |
real_relations(long LIM, long lim, long LIMC, long **mat) |
GEN vecexpo) |
|
{ |
{ |
static long nbtest; |
long lgsub = lg(subFB), i, s, fpc, current, nbtest = 0, endcycle, rhoacc, rho; |
long av = avma,av1,i,j,p,fpc,b1,b2,ep,current, first = (s==0); |
int first = 1; |
long *col,*fpd,*oldfact,*oldexp,limstack; |
gpmem_t av, av1, limstack; |
long findecycle,nbrhocumule,nbrho; |
GEN C, d, col, form, form0, form1, ex = cgetg(lgsub, t_VECSMALL); |
GEN p1,p2,form,form0,form1,form2; |
|
|
|
limstack=stack_lim(av,1); |
C = cgetg(LIM+1, t_VEC); |
if (first) nbtest = 0; |
for (i=1; i<=LIM; i++) C[i] = lgetr(PRECREG); |
|
|
current = 0; |
current = 0; |
p1 = NULL; /* gcc -Wall */ |
av = avma; limstack = stack_lim(av,1); |
while (s<lim) |
s = trivial_relations(mat, C); |
|
lim += s; if (lim > LIM) lim = LIM; |
|
NEW: |
|
for(;;) |
{ |
{ |
|
if (s >= lim) { |
|
if (lim == LIM) break; |
|
lim = LIM; first = 0; |
|
if (DEBUGLEVEL) dbg_all(0, s, nbtest); |
|
} |
|
avma = av; |
form = real_random_form(ex); |
form = real_random_form(ex); |
if (!first) |
if (!first) |
{ |
{ |
current = 1+s-RELSUP; |
current = 1+s-RELSUP; |
p1 = redrealprimeform(Disc, factorbase[current]); |
form = comprealform3(form, realpf(Disc, FB[current])); |
form = comprealform3(form,p1); |
|
} |
} |
|
av1 = avma; |
form0 = form; form1 = NULL; |
form0 = form; form1 = NULL; |
findecycle = nbrhocumule = 0; |
endcycle = rhoacc = 0; |
nbrho = -1; av1 = avma; |
rho = -1; |
while (s<lim) |
|
|
CYCLE: |
|
if (endcycle) goto NEW; |
|
if (low_stack(limstack, stack_lim(av,1))) |
{ |
{ |
if (low_stack(limstack, stack_lim(av,1))) |
if(DEBUGMEM>1) err(warnmem,"real_relations"); |
{ |
gerepileall(av1, form1? 2: 1, &form, &form1); |
GEN *gptr[2]; |
} |
int c = 1; |
if (rho < 0) rho = 0; /* first time in */ |
if(DEBUGMEM>1) err(warnmem,"real_relations"); |
else |
gptr[0]=&form; if (form1) gptr[c++]=&form1; |
{ |
gerepilemany(av1,gptr,c); |
form = rhorealform(form); rho++; |
} |
rhoacc++; |
if (nbrho < 0) nbrho = 0; /* first time in */ |
if (first) |
|
endcycle = (absi_equal((GEN)form[1],(GEN)form0[1]) |
|
&& egalii((GEN)form[2],(GEN)form0[2]) |
|
&& (!narrow || signe(form0[1])==signe(form[1]))); |
else |
else |
{ |
{ |
if (findecycle) break; |
if (narrow) |
form = rhorealform(form); |
{ form = rhorealform(form); rho++; } |
nbrho++; nbrhocumule++; |
else if (absi_equal((GEN)form[1], (GEN)form[3])) /* a = -c */ |
if (first) |
|
{ |
{ |
if (absi_equal((GEN)form[1],(GEN)form0[1]) |
if (absi_equal((GEN)form[1],(GEN)form0[1]) && |
&& egalii((GEN)form[2],(GEN)form0[2]) |
egalii((GEN)form[2],(GEN)form0[2])) goto NEW; |
&& (!narrow || signe(form0[1])==signe(form[1]))) findecycle=1; |
form = rhorealform(form); rho++; |
} |
} |
else |
else |
{ |
{ setsigne(form[1],1); setsigne(form[3],-1); } |
if (narrow) |
if (egalii((GEN)form[1],(GEN)form0[1]) && |
{ form=rhorealform(form); nbrho++; } |
egalii((GEN)form[2],(GEN)form0[2])) goto NEW; |
else if (absi_equal((GEN)form[1], (GEN)form[3])) /* a = -c */ |
|
{ |
|
if (absi_equal((GEN)form[1],(GEN)form0[1]) && |
|
egalii((GEN)form[2],(GEN)form0[2])) break; |
|
form=rhorealform(form); nbrho++; |
|
} |
|
else |
|
{ setsigne(form[1],1); setsigne(form[3],-1); } |
|
if (egalii((GEN)form[1],(GEN)form0[1]) && |
|
egalii((GEN)form[2],(GEN)form0[2])) break; |
|
} |
|
} |
} |
nbtest++; fpc = factorisequad(form,KC,LIMC); |
} |
if (!fpc) |
nbtest++; fpc = factorquad(form,KC,LIMC); |
|
if (!fpc) |
|
{ |
|
if (DEBUGLEVEL>1) fprintferr("."); |
|
goto CYCLE; |
|
} |
|
if (fpc > 1) |
|
{ /* look for Large Prime relation */ |
|
long *fpd = largeprime(fpc,ex,current,rhoacc); |
|
long b1, b2, p; |
|
GEN form2; |
|
if (!fpd) |
{ |
{ |
if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); } |
if (DEBUGLEVEL>1) fprintferr("."); |
continue; |
goto CYCLE; |
} |
} |
if (fpc > 1) |
if (!form1) form1 = initrealform5(ex); |
|
if (!first) form1 = comprealform(form1, realpf5(Disc, FB[current])); |
|
form1 = rhoreal_pow(form1, rho); |
|
rho = 0; |
|
|
|
form2 = initrealform5(fpd); |
|
if (fpd[-2]) form2 = comprealform(form2, realpf5(Disc, FB[fpd[-2]])); |
|
form2 = rhoreal_pow(form2, fpd[-3]); |
|
if (!narrow && !absi_equal((GEN)form2[1],(GEN)form2[3])) |
{ |
{ |
fpd = largeprime(fpc,ex,current,nbrhocumule); |
setsigne(form2[1],1); |
if (!fpd) |
setsigne(form2[3],-1); |
{ |
} |
if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); } |
p = fpc << 1; |
continue; |
b1 = smodis((GEN)form2[2], p); |
} |
b2 = smodis((GEN)form1[2], p); |
if (!form1) form1 = initrealform5(ex); |
if (b1 != b2 && b1+b2 != p) goto CYCLE; |
if (!first) |
|
{ |
|
p1 = redrealprimeform5(Disc, factorbase[current]); |
|
form1=comprealform(form1,p1); |
|
} |
|
form1 = rhoreal_pow(form1, nbrho); nbrho = 0; |
|
form2 = initrealform5(fpd); |
|
if (fpd[-2]) |
|
{ |
|
p1 = redrealprimeform5(Disc, factorbase[fpd[-2]]); |
|
form2=comprealform(form2,p1); |
|
} |
|
form2 = rhoreal_pow(form2, fpd[-3]); |
|
if (!narrow && !absi_equal((GEN)form2[1],(GEN)form2[3])) |
|
{ |
|
setsigne(form2[1],1); |
|
setsigne(form2[3],-1); |
|
} |
|
p = fpc << 1; |
|
b1=smodis((GEN)form2[2], p); |
|
b2=smodis((GEN)form1[2], p); |
|
if (b1 != b2 && b1+b2 != p) continue; |
|
|
|
s++; col = mat[s]; if (DEBUGLEVEL) fprintferr(" %ld",s); |
col = mat[++s]; |
oldfact = primfact; oldexp = exprimfact; |
add_fact(col, form1); |
primfact = primfact1; exprimfact = exprimfact1; |
(void)factorquad(form2,KC,LIMC); |
factorisequad(form2,KC,LIMC); |
if (b1==b2) |
if (b1==b2) |
{ |
{ |
for (i=1; i<lgsub; i++) col[subFB[i]] += fpd[i]-ex[i]; |
for (i=1; i<=lgsub; i++) |
sub_fact(col, form2); |
col[numfactorbase[subbase[i]]] = fpd[i]-ex[i]; |
if (fpd[-2]) col[fpd[-2]]++; /* implies !first */ |
for (j=1; j<=primfact[0]; j++) |
d = get_dist(subii((GEN)form1[4],(GEN)form2[4]), |
{ |
divrr((GEN)form1[5],(GEN)form2[5]), PRECREG); |
p=primfact[j]; ep=exprimfact[j]; |
|
if (smodis((GEN)form2[2], p<<1) > p) ep = -ep; |
|
col[numfactorbase[p]] -= ep; |
|
} |
|
if (fpd[-2]) col[fpd[-2]]++; /* implies !first */ |
|
p1 = subii((GEN)form1[4],(GEN)form2[4]); |
|
p2 = divrr((GEN)form1[5],(GEN)form2[5]); |
|
} |
|
else |
|
{ |
|
for (i=1; i<=lgsub; i++) |
|
col[numfactorbase[subbase[i]]] = -fpd[i]-ex[i]; |
|
for (j=1; j<=primfact[0]; j++) |
|
{ |
|
p=primfact[j]; ep=exprimfact[j]; |
|
if (smodis((GEN)form2[2], p<<1) > p) ep = -ep; |
|
col[numfactorbase[p]] += ep; |
|
} |
|
if (fpd[-2]) col[fpd[-2]]--; |
|
p1 = addii((GEN)form1[4],(GEN)form2[4]); |
|
p2 = mulrr((GEN)form1[5],(GEN)form2[5]); |
|
} |
|
primfact = oldfact; exprimfact = oldexp; |
|
} |
} |
else |
else |
{ |
{ |
if (!form1) form1 = initrealform5(ex); |
for (i=1; i<lgsub; i++) col[subFB[i]] += -fpd[i]-ex[i]; |
if (!first) |
add_fact(col, form2); |
{ |
if (fpd[-2]) col[fpd[-2]]--; |
p1 = redrealprimeform5(Disc, factorbase[current]); |
d = get_dist(addii((GEN)form1[4],(GEN)form2[4]), |
form1=comprealform(form1,p1); |
mulrr((GEN)form1[5],(GEN)form2[5]), PRECREG); |
} |
|
form1 = rhoreal_pow(form1,nbrho); nbrho = 0; |
|
|
|
s++; col = mat[s]; if (DEBUGLEVEL) fprintferr(" %ld",s); |
|
for (i=1; i<=lgsub; i++) |
|
col[numfactorbase[subbase[i]]] = -ex[i]; |
|
p1 = (GEN) form1[4]; |
|
p2 = (GEN) form1[5]; |
|
} |
} |
for (j=1; j<=primfact[0]; j++) |
} |
|
else |
|
{ /* standard relation */ |
|
if (!form1) form1 = initrealform5(ex); |
|
if (!first) form1 = comprealform(form1, realpf5(Disc, FB[current])); |
|
form1 = rhoreal_pow(form1,rho); |
|
rho = 0; |
|
|
|
col = mat[++s]; |
|
for (i=1; i<lgsub; i++) col[subFB[i]] = -ex[i]; |
|
add_fact(col, form1); |
|
d = get_dist((GEN)form1[4], (GEN)form1[5], PRECREG); |
|
} |
|
if (DEBUGLEVEL) fprintferr(" %ld",s); |
|
affrr(d, (GEN)C[s]); |
|
if (first) |
|
{ |
|
if (s >= lim) goto NEW; |
|
goto CYCLE; |
|
} |
|
else |
|
{ |
|
col[current]--; |
|
if (fpc == 1 && col[current] == 0) |
{ |
{ |
p=primfact[j]; ep=exprimfact[j]; |
s--; for (i=1; i<=KC; i++) col[i]=0; |
if (smodis((GEN)form1[2], p<<1) > p) ep = -ep; |
|
col[numfactorbase[p]] += ep; |
|
} |
} |
p1 = mpadd(mulir(mulsi(EXP220,p1), glog2), mplog(absr(p2))); |
|
affrr(shiftr(p1,-1), (GEN)vecexpo[s]); |
|
if (!first) |
|
{ |
|
col[current]--; |
|
if (fpc == 1 && col[current] == 0) |
|
{ s--; for (i=1; i<=KC; i++) col[i]=0; } |
|
break; |
|
} |
|
} |
} |
avma = av; |
|
} |
} |
if (DEBUGLEVEL) |
if (DEBUGLEVEL) dbg_all(1, s, nbtest); |
{ |
return C; |
fprintferr("\nnbrelations/nbtest = %ld/%ld\n",s,nbtest); |
|
msgtimer("%s relations", first? "initial": "random"); |
|
} |
|
} |
} |
|
|
static int |
static int |
real_be_honest(long *ex) |
real_be_honest() |
{ |
{ |
long p,fpc,s = KC,nbtest = 0,av = avma; |
long p, fpc, s = KC, nbtest = 0; |
GEN p1,form,form0; |
gpmem_t av = avma; |
|
GEN F,F0, ex = cgetg(lg(subFB), t_VECSMALL); |
|
|
while (s<KC2) |
while (s<KC2) |
{ |
{ |
p = factorbase[s+1]; |
p = FB[s+1]; if (DEBUGLEVEL) fprintferr(" %ld",p); |
if (DEBUGLEVEL) { fprintferr(" %ld",p); flusherr(); } |
F = comprealform3(realpf(Disc, p), real_random_form(ex)); |
form = real_random_form(ex); |
for (F0 = F;;) |
p1 = redrealprimeform(Disc, p); |
|
form0 = form = comprealform3(form,p1); |
|
for(;;) |
|
{ |
{ |
fpc = factorisequad(form,s,p-1); |
fpc = factorquad(F,s,p-1); |
if (fpc == 1) { nbtest=0; s++; break; } |
if (fpc == 1) { nbtest=0; s++; break; } |
form = rhorealform(form); |
|
if (++nbtest > 20) return 0; |
if (++nbtest > 20) return 0; |
if (narrow || absi_equal((GEN)form[1],(GEN)form[3])) |
F = fix_signs( rhorealform(F) ); |
form = rhorealform(form); |
if (egalii((GEN)F[1],(GEN)F0[1]) |
else |
&& egalii((GEN)F[2],(GEN)F0[2])) break; |
{ |
|
setsigne(form[1],1); |
|
setsigne(form[3],-1); |
|
} |
|
if (egalii((GEN)form[1],(GEN)form0[1]) |
|
&& egalii((GEN)form[2],(GEN)form0[2])) break; |
|
} |
} |
avma = av; |
avma = av; |
} |
} |
Line 1597 real_be_honest(long *ex) |
|
Line 1582 real_be_honest(long *ex) |
|
} |
} |
|
|
static GEN |
static GEN |
gcdrealnoer(GEN a,GEN b,long *pte) |
gcdreal(GEN a,GEN b) |
{ |
{ |
long e; |
long e; |
GEN k1,r; |
GEN k1,r; |
|
|
|
if (!signe(a)) return mpabs(b); |
|
if (!signe(b)) return mpabs(a); |
|
|
if (typ(a)==t_INT) |
if (typ(a)==t_INT) |
{ |
{ |
if (typ(b)==t_INT) return mppgcd(a,b); |
if (typ(b)==t_INT) return mppgcd(a,b); |
k1=cgetr(lg(b)); affir(a,k1); a=k1; |
a = itor(a, lg(b)); |
} |
} |
else if (typ(b)==t_INT) |
else if (typ(b)==t_INT) |
{ k1=cgetr(lg(a)); affir(b,k1); b=k1; } |
{ |
|
b = itor(b, lg(a)); |
|
} |
if (expo(a)<-5) return absr(b); |
if (expo(a)<-5) return absr(b); |
if (expo(b)<-5) return absr(a); |
if (expo(b)<-5) return absr(a); |
a=absr(a); b=absr(b); |
a=absr(a); b=absr(b); |
while (expo(b) >= -5 && signe(b)) |
while (expo(b) >= -5 && signe(b)) |
{ |
{ |
k1=gcvtoi(divrr(a,b),&e); |
k1 = gcvtoi(divrr(a,b),&e); |
if (e > 0) return NULL; |
if (e > 0) return NULL; |
r=subrr(a,mulir(k1,b)); a=b; b=r; |
r=subrr(a,mulir(k1,b)); a=b; b=r; |
} |
} |
*pte=expo(b); return absr(a); |
return absr(a); |
} |
} |
|
|
static GEN |
enum { RELAT, LARGE, PRECI }; |
get_reg(GEN matc, long sreg) |
|
|
static int |
|
get_R(GEN C, long sreg, GEN z, GEN *ptR) |
{ |
{ |
long i,e,maxe; |
GEN R = gun; |
GEN reg = mpabs(gcoeff(matc,1,1)); |
double c; |
|
long i; |
|
|
e = maxe = 0; |
if (PRECREG) |
for (i=2; i<=sreg; i++) |
|
{ |
{ |
reg = gcdrealnoer(gcoeff(matc,1,i),reg,&e); |
R = mpabs((GEN)C[1]); |
if (!reg) return NULL; |
for (i=2; i<=sreg; i++) |
maxe = maxe? max(maxe,e): e; |
{ |
|
R = gcdreal((GEN)C[i], R); |
|
if (!R) return PRECI; |
|
} |
|
if (DEBUGLEVEL) |
|
{ |
|
if (DEBUGLEVEL>7) fprintferr("R = %Z",R); |
|
msgtimer("regulator"); |
|
} |
|
if (gexpo(R) <= -3) |
|
{ |
|
if (DEBUGLEVEL) fprintferr("regulator is zero.\n"); |
|
return RELAT; |
|
} |
} |
} |
|
c = gtodouble(gmul(z, R)); |
|
if (c < 0.75 || c > 1.5) return RELAT; |
|
*ptR = R; return LARGE; |
|
} |
|
|
|
static int |
|
quad_be_honest() |
|
{ |
|
int r; |
|
if (KC2 <= KC) return 1; |
if (DEBUGLEVEL) |
if (DEBUGLEVEL) |
{ |
fprintferr("be honest for primes from %ld to %ld\n", FB[KC+1],FB[KC2]); |
if (DEBUGLEVEL>7) { fprintferr("reg = "); outerr(reg); } |
r = PRECREG? real_be_honest(): imag_be_honest(); |
msgtimer("regulator"); |
if (DEBUGLEVEL) { fprintferr("\n"); msgtimer("be honest"); } |
} |
return r; |
return reg; |
|
} |
} |
|
|
GEN |
GEN |
buchquad(GEN D, double cbach, double cbach2, long RELSUP0, long flag, long prec) |
buchquad(GEN D, double cbach, double cbach2, long RELSUP0, long flag, long prec) |
{ |
{ |
long av0 = avma,av,tetpil,KCCO,KCCOPRO,i,j,s, *ex,**mat; |
gpmem_t av0 = avma, av; |
long extrarel,nrelsup,nreldep,LIMC,LIMC2,cp,nbram,nlze; |
long KCCO, i, j, s, **mat; |
GEN p1,h,W,met,res,basecl,dr,c_1,pdep,C,B,extramat,extraC; |
long nrelsup, nreldep, LIMC, LIMC2, cp, nlze; |
GEN reg,vecexpo,glog2,cst; |
GEN h, W, cyc, res, gen, dep, C, B, extramat, extraC; |
double drc,lim,LOGD; |
GEN R, resc, Res, z; |
|
double drc, lim, LOGD, LOGD2; |
|
const long MAXRELSUP = 7; |
|
|
Disc = D; if (typ(Disc)!=t_INT) err(typeer,"buchquad"); |
Disc = D; if (typ(Disc)!=t_INT) err(typeer,"buchquad"); |
s=mod4(Disc); |
s = mod4(Disc); |
glog2 = vecexpo = NULL; /* gcc -Wall */ |
|
switch(signe(Disc)) |
switch(signe(Disc)) |
{ |
{ |
case -1: |
case -1: |
if (cmpis(Disc,-4) >= 0) |
if (cmpis(Disc,-4) >= 0) |
{ |
{ |
p1=cgetg(6,t_VEC); p1[1]=p1[4]=p1[5]=un; |
GEN p1=cgetg(6,t_VEC); p1[1]=p1[4]=p1[5]=un; |
p1[2]=p1[3]=lgetg(1,t_VEC); return p1; |
p1[2]=p1[3]=lgetg(1,t_VEC); return p1; |
} |
} |
if (s==2 || s==1) err(funder2,"buchquad"); |
if (s==2 || s==1) err(funder2,"buchquad"); |
Line 1677 buchquad(GEN D, double cbach, double cbach2, long RELS |
|
Line 1692 buchquad(GEN D, double cbach, double cbach2, long RELS |
|
if (!isfundamental(Disc)) |
if (!isfundamental(Disc)) |
err(warner,"not a fundamental discriminant in quadclassunit"); |
err(warner,"not a fundamental discriminant in quadclassunit"); |
buch_init(); RELSUP = RELSUP0; |
buch_init(); RELSUP = RELSUP0; |
dr=cgetr(3); affir(Disc,dr); drc=fabs(rtodbl(dr)); LOGD=log(drc); |
drc = fabs(gtodouble(Disc)); |
lim=sqrt(drc); cst = mulrr(lfunc(Disc), dbltor(lim)); |
LOGD = log(drc); |
|
LOGD2 = LOGD * LOGD; |
|
|
|
lim = sqrt(drc); |
|
/* resc = sqrt(D) w / 2^r1 (2pi)^r2 ~ hR / L(chi,1) */ |
|
if (PRECREG) resc = dbltor(lim / 2.); |
|
else resc = dbltor(lim / PI); |
if (!PRECREG) lim /= sqrt(3.); |
if (!PRECREG) lim /= sqrt(3.); |
|
R = gun; |
cp = (long)exp(sqrt(LOGD*log(LOGD)/8.0)); |
cp = (long)exp(sqrt(LOGD*log(LOGD)/8.0)); |
if (cp < 13) cp = 13; |
if (cp < 13) cp = 13; |
av = avma; cbach /= 2; |
av = avma; cbach /= 2; |
|
mat = NULL; |
|
|
INCREASE: avma = av; cbach = check_bach(cbach,6.); |
START: avma = av; cbach = check_bach(cbach,6.); |
|
if (mat) { desalloc(mat); mat = NULL; } |
nreldep = nrelsup = 0; |
nreldep = nrelsup = 0; |
LIMC = (long)(cbach*LOGD*LOGD); |
LIMC = (long)(cbach*LOGD2); |
if (LIMC < cp) LIMC=cp; |
if (LIMC < cp) { LIMC = cp; cbach = LIMC / LOGD2; } |
LIMC2 = max(20, (long)(max(cbach,cbach2)*LOGD*LOGD)); |
LIMC2 = (long)(max(cbach,cbach2)*LOGD2); |
if (LIMC2 < LIMC) LIMC2 = LIMC; |
if (LIMC2 < LIMC) LIMC2 = LIMC; |
if (PRECREG) |
if (PRECREG) |
{ |
{ |
PRECREG = max(prec+1, MEDDEFAULTPREC + 2*(expi(Disc)>>TWOPOTBITS_IN_LONG)); |
PRECREG = max(prec+1, MEDDEFAULTPREC + 2*(expi(Disc)>>TWOPOTBITS_IN_LONG)); |
glog2 = glog(gdeux,PRECREG); |
|
sqrtD = gsqrt(Disc,PRECREG); |
sqrtD = gsqrt(Disc,PRECREG); |
isqrtD = gfloor(sqrtD); |
isqrtD = gfloor(sqrtD); |
} |
} |
factorbasequad(Disc,LIMC2,LIMC); |
|
if (!KC) goto INCREASE; |
|
|
|
vperm = cgetg(KC+1,t_VECSMALL); for (i=1; i<=KC; i++) vperm[i]=i; |
Res = FBquad(Disc,LIMC2,LIMC); |
nbram = subfactorbasequad(lim,KC); |
if (!Res) goto START; |
if (nbram == -1) { desalloc(NULL); goto INCREASE; } |
subFB = subFBquad(Disc, lim + 0.5, KC); |
KCCO = KC + RELSUP; |
if (!subFB) goto START; |
if (DEBUGLEVEL) { fprintferr("KC = %ld, KCCO = %ld\n",KC,KCCO); flusherr(); } |
powsubFB = powsubFBquad(CBUCH+1); |
powsubfact(lgsub,CBUCH+7); |
|
|
|
mat = (long**) gpmalloc((KCCO+1)*sizeof(long*)); |
|
setlg(mat, KCCO+1); |
|
for (i=1; i<=KCCO; i++) |
|
{ |
|
mat[i] = (long*) gpmalloc((KC+1)*sizeof(long)); |
|
for (j=1; j<=KC; j++) mat[i][j]=0; |
|
} |
|
ex = new_chunk(lgsub+1); |
|
limhash = ((ulong)LIMC < (MAXHALFULONG>>1))? LIMC*LIMC: HIGHBIT>>1; |
limhash = ((ulong)LIMC < (MAXHALFULONG>>1))? LIMC*LIMC: HIGHBIT>>1; |
for (i=0; i<HASHT; i++) hashtab[i]=NULL; |
for (i=0; i<HASHT; i++) hashtab[i]=NULL; |
|
|
s = lgsub+nbram+RELSUP; |
KCCO = KC + RELSUP; |
if (PRECREG) |
if (DEBUGLEVEL) fprintferr("KC = %ld, KCCO = %ld\n",KC,KCCO); |
|
mat = (long**)cgetalloc(NULL, KCCO+1, t_VEC); |
|
for (i=1; i<=KCCO; i++) |
{ |
{ |
vecexpo=cgetg(KCCO+1,t_VEC); |
GEN t = cgetalloc(NULL, KC+1, t_VECSMALL); |
for (i=1; i<=KCCO; i++) vecexpo[i]=lgetr(PRECREG); |
for (j=1; j<=KC; j++) t[j]=0; |
real_relations(s,0,LIMC,ex,mat,glog2,vecexpo); |
mat[i] = t; |
real_relations(KCCO,s,LIMC,ex,mat,glog2,vecexpo); |
|
} |
} |
else |
|
{ |
|
imag_relations(s,0,LIMC,ex,mat); |
|
imag_relations(KCCO,s,LIMC,ex,mat); |
|
} |
|
if (KC2 > KC) |
|
{ |
|
if (DEBUGLEVEL) |
|
fprintferr("be honest for primes from %ld to %ld\n", |
|
factorbase[KC+1],factorbase[KC2]); |
|
s = PRECREG? real_be_honest(ex): imag_be_honest(ex); |
|
if (DEBUGLEVEL) |
|
{ |
|
fprintferr("\n"); |
|
msgtimer("be honest"); |
|
} |
|
if (!s) { desalloc(mat); goto INCREASE; } |
|
} |
|
C=cgetg(KCCO+1,t_MAT); |
|
if (PRECREG) |
|
{ |
|
for (i=1; i<=KCCO; i++) |
|
{ |
|
C[i]=lgetg(2,t_COL); coeff(C,1,i)=vecexpo[i]; |
|
} |
|
if (DEBUGLEVEL>7) { fprintferr("C: "); outerr(C); flusherr(); } |
|
} |
|
else |
|
for (i=1; i<=KCCO; i++) C[i]=lgetg(1,t_COL); |
|
W = hnfspec(mat,vperm,&pdep,&B,&C,lgsub); |
|
nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1; |
|
|
|
KCCOPRO=KCCO; |
s = lg(subFB)-1 + RELSUP; |
|
C = PRECREG? real_relations(KCCO,s,LIMC,mat) |
|
: imag_relations(KCCO,s,LIMC,mat); |
|
W = hnfspec(mat,vperm,&dep,&B,&C,lg(subFB)-1); |
|
nlze = lg(dep)>1? lg(dep[1])-1: lg(B[1])-1; |
|
|
if (nlze) |
if (nlze) |
{ |
{ |
EXTRAREL: |
MORE: |
s = PRECREG? 2: 1; extrarel=nlze+2; |
extramat = extra_relations(LIMC,nlze, &extraC); |
extraC=cgetg(extrarel+1,t_MAT); |
if (!extramat) { goto START; } |
for (i=1; i<=extrarel; i++) extraC[i]=lgetg(s,t_COL); |
W = hnfadd(W,vperm,&dep,&B,&C, extramat,extraC); |
extramat = extra_relations(LIMC,ex,nlze,extraC); |
nlze = lg(dep)>1? lg(dep[1])-1: lg(B[1])-1; |
if (!extramat) { desalloc(mat); goto INCREASE; } |
KCCO += lg(extramat)-1; |
W = hnfadd(W,vperm,&pdep,&B,&C, extramat,extraC); |
|
nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1; |
|
KCCOPRO += extrarel; |
|
if (nlze) |
if (nlze) |
{ |
{ |
if (++nreldep > 5) { desalloc(mat); goto INCREASE; } |
if (++nreldep > 5) goto START; |
goto EXTRAREL; |
goto MORE; |
} |
} |
} |
} |
/* tentative class number */ |
|
h=gun; for (i=1; i<lg(W); i++) h=mulii(h,gcoeff(W,i,i)); |
|
if (PRECREG) |
|
{ |
|
/* tentative regulator */ |
|
reg = get_reg(C, KCCOPRO - (lg(B)-1) - (lg(W)-1)); |
|
if (!reg) |
|
{ |
|
desalloc(mat); |
|
prec = (PRECREG<<1)-2; goto INCREASE; |
|
} |
|
if (gexpo(reg)<=-3) |
|
{ |
|
if (++nrelsup <= 7) |
|
{ |
|
if (DEBUGLEVEL) { fprintferr("regulateur nul\n"); flusherr(); } |
|
nlze=min(KC,nrelsup); goto EXTRAREL; |
|
} |
|
desalloc(mat); goto INCREASE; |
|
} |
|
c_1 = divrr(gmul2n(gmul(h,reg),1), cst); |
|
} |
|
else |
|
{ |
|
reg = gun; |
|
c_1 = divrr(gmul(h,mppi(DEFAULTPREC)), cst); |
|
} |
|
|
|
if (gcmpgs(gmul2n(c_1,2),3)<0) { c_1=stoi(10); nrelsup=7; } |
h = dethnf_i(W); |
if (gcmpgs(gmul2n(c_1,1),3)>0) |
if (DEBUGLEVEL) fprintferr("\n#### Tentative class number: %Z\n", h); |
|
|
|
z = mulrr(Res, resc); /* ~ hR if enough relations, a multiple otherwise */ |
|
switch(get_R(C, KCCO - (lg(B)-1) - (lg(W)-1), divir(h,z), &R)) |
{ |
{ |
nrelsup++; |
case PRECI: |
if (nrelsup<=7) |
prec = (PRECREG<<1)-2; |
{ |
goto START; |
if (DEBUGLEVEL) |
|
{ fprintferr("***** check = %f\n\n",gtodouble(c_1)); flusherr(); } |
case RELAT: |
nlze=min(KC,nrelsup); goto EXTRAREL; |
if (++nrelsup <= MAXRELSUP) { nlze = min(KC, nrelsup); goto MORE; } |
} |
goto START; |
if (cbach < 5.99) { desalloc(mat); goto INCREASE; } |
|
err(warner,"suspicious check. Suggest increasing extra relations."); |
|
} |
} |
basecl = get_clgp(Disc,W,&met,PRECREG); |
/* DONE */ |
s = lg(basecl); desalloc(mat); tetpil=avma; |
if (!quad_be_honest()) goto START; |
|
|
|
gen = get_clgp(Disc,W,&cyc,PRECREG); |
|
desalloc(mat); |
|
|
res=cgetg(6,t_VEC); |
res=cgetg(6,t_VEC); |
res[1]=lcopy(h); p1=cgetg(s,t_VEC); |
res[1]=(long)h; |
for (i=1; i<s; i++) p1[i] = (long)icopy(gcoeff(met,i,i)); |
res[2]=(long)cyc; |
res[2]=(long)p1; |
res[3]=(long)gen; |
res[3]=lcopy(basecl); |
res[4]=(long)R; |
res[4]=lcopy(reg); |
res[5]=ldiv(mpmul(R,h), z); return gerepilecopy(av0,res); |
res[5]=lcopy(c_1); return gerepile(av0,tetpil,res); |
|
} |
} |
|
|
GEN |
GEN |