version 1.1, 2001/10/02 11:17:02 |
version 1.2, 2002/09/11 07:26:49 |
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 ZV_lincomb(GEN u, GEN v, GEN X, GEN Y); |
extern GEN ZV_lincomb(GEN u, GEN v, GEN X, GEN Y); |
extern GEN roots_to_pol_r1r2(GEN a, long r1, long v); |
extern int addcolumntomatrix(GEN V,GEN INVP,GEN L); |
extern GEN makebasis(GEN nf,GEN pol); |
extern long expodb(double x); |
extern GEN caractducos(GEN p, GEN x, int v); |
|
|
|
|
/* default quality ratio for LLL: 99/100 */ |
|
#define LLLDFT 100 |
|
|
/* scalar product x.x */ |
/* scalar product x.x */ |
GEN |
GEN |
sqscal(GEN x) |
sqscal(GEN x) |
{ |
{ |
long i,av,lx; |
long i, lx; |
|
gpmem_t av; |
GEN z; |
GEN z; |
lx = lg(x); |
lx = lg(x); |
if (lx == 1) return gzero; |
if (lx == 1) return gzero; |
|
|
GEN |
GEN |
gscal(GEN x,GEN y) |
gscal(GEN x,GEN y) |
{ |
{ |
long i,av,lx; |
long i, lx; |
|
gpmem_t av; |
GEN z; |
GEN z; |
if (x == y) return sqscal(x); |
if (x == y) return sqscal(x); |
lx = lg(x); |
lx = lg(x); |
Line 59 gscal(GEN x,GEN y) |
|
Line 63 gscal(GEN x,GEN y) |
|
static GEN |
static GEN |
sqscali(GEN x) |
sqscali(GEN x) |
{ |
{ |
long i,av,lx; |
long i, lx; |
|
gpmem_t av; |
GEN z; |
GEN z; |
lx = lg(x); |
lx = lg(x); |
if (lx == 1) return gzero; |
if (lx == 1) return gzero; |
|
|
static GEN |
static GEN |
gscali(GEN x,GEN y) |
gscali(GEN x,GEN y) |
{ |
{ |
long i,av,lx; |
long i, lx; |
|
gpmem_t av; |
GEN z; |
GEN z; |
if (x == y) return sqscali(x); |
if (x == y) return sqscali(x); |
lx = lg(x); |
lx = lg(x); |
Line 85 gscali(GEN x,GEN y) |
|
Line 91 gscali(GEN x,GEN y) |
|
return gerepileuptoint(av,z); |
return gerepileuptoint(av,z); |
} |
} |
|
|
|
/********************************************************************/ |
|
/** QR Factorization via Householder matrices **/ |
|
/********************************************************************/ |
|
|
|
/* zero x[1..k-1], fill mu */ |
|
static int |
|
FindApplyQ(GEN x, GEN mu, GEN B, long k, GEN Q, long prec) |
|
{ |
|
long i, lx = lg(x)-1, lv = lx - (k-1); |
|
GEN v, p1, beta, Nx, x2, x1, xd = x + (k-1); |
|
|
|
x1 = (GEN)xd[1]; |
|
x2 = gsqr(x1); |
|
if (k < lx) |
|
{ |
|
for (i=2; i<=lv; i++) x2 = mpadd(x2, gsqr((GEN)xd[i])); |
|
Nx = gsqrt(x2, prec); |
|
if (signe(x1) < 0) setsigne(Nx, -1); |
|
v = cgetg(lv+1, t_VEC); |
|
v[1] = lmpadd(x1, Nx); |
|
for (i=2; i<=lv; i++) v[i] = xd[i]; |
|
if (gcmp0(x2)) return 0; |
|
|
|
if (gcmp0(x1)) |
|
beta = mpmul(x2, realun(prec)); /* make sure typ(beta) != t_INT */ |
|
else |
|
beta = mpadd(x2, mpmul(Nx,x1)); |
|
beta = ginv(beta); |
|
p1 = cgetg(3,t_VEC); Q[k] = (long)p1; |
|
p1[1] = (long)beta; |
|
p1[2] = (long)v; |
|
|
|
coeff(mu,k,k) = lmpneg(Nx); |
|
} |
|
else |
|
coeff(mu,k,k) = x[k]; |
|
if (B) |
|
{ |
|
B[k] = (long)x2; |
|
for (i=1; i<k; i++) coeff(mu,k,i) = x[i]; |
|
} |
|
else |
|
for (i=1; i<k; i++) coeff(mu,i,k) = x[i]; |
|
return (typ(x2) != t_REAL || lg(x2) > 3 || expo(x2) < 10); |
|
} |
|
|
|
static void |
|
ApplyQ(GEN Q, GEN r) |
|
{ |
|
GEN s, rd, beta = (GEN)Q[1], v = (GEN)Q[2]; |
|
long i, l = lg(v), lr = lg(r); |
|
|
|
rd = r + (lr - l); |
|
s = mpmul((GEN)v[1], (GEN)rd[1]); |
|
for (i=2; i<l; i++) s = mpadd(s, mpmul((GEN)v[i] ,(GEN)rd[i])); |
|
s = mpneg(mpmul(beta, s)); |
|
for (i=1; i<l; i++) rd[i] = lmpadd((GEN)rd[i], mpmul(s, (GEN)v[i])); |
|
} |
|
|
static GEN |
static GEN |
lllall_trivial(GEN x, long n, long flag) |
ApplyAllQ(GEN Q, GEN r0, long k) |
{ |
{ |
GEN y; |
gpmem_t av = avma; |
if (!n) |
GEN r = dummycopy(r0); |
|
long j; |
|
for (j=1; j<k; j++) ApplyQ((GEN)Q[j], r); |
|
return gerepilecopy(av, r); |
|
} |
|
|
|
/* compute B[k] = | x[k] |^2, update mu(k, 1..k-1) using Householder matrices |
|
* (Q = Householder(x[1..k-1]) in factored form) */ |
|
static int |
|
incrementalQ(GEN x, GEN L, GEN B, GEN Q, long k, long prec) |
|
{ |
|
GEN r = ApplyAllQ(Q, (GEN)x[k], k); |
|
return FindApplyQ(r, L, B, k, Q, prec); |
|
} |
|
|
|
/* Q vector of Householder matrices orthogonalizing x[1..j0]. |
|
* Q[i] = 0 means not computed yet */ |
|
static int |
|
Householder_get_mu(GEN x, GEN L, GEN B, long k, GEN Q, long prec) |
|
{ |
|
GEN Nx, invNx, m; |
|
long i, j, j0; |
|
if (!Q) Q = zerovec(k); |
|
for (j=1; j<=k; j++) |
|
if (typ(Q[j]) == t_INT) break; |
|
j0 = j; |
|
for ( ; j<=k; j++) |
|
if (! incrementalQ(x, L, B, Q, j, prec)) return 0; |
|
for (j=1; j<=k; j++) |
{ |
{ |
if (flag != lll_ALL) return cgetg(1,t_MAT); |
m = (GEN)L[j]; Nx = (GEN)m[j]; /* should set m[j] = un; but need it later */ |
|
if (j == k) break; |
|
invNx = ginv(Nx); |
|
for (i=max(j0, j+1); i<=k; i++) m[i] = lmpmul(invNx, (GEN)m[i]); |
|
} |
|
return 1; |
|
} |
|
|
|
GEN |
|
sqred1_from_QR(GEN x, long prec) |
|
{ |
|
long j, k = lg(x)-1; |
|
GEN L, B = zerovec(k); |
|
L = cgetg(k+1, t_MAT); |
|
for (j=1; j<=k; j++) L[j] = (long)zerocol(k); |
|
if (!Householder_get_mu(x, L, B, k, NULL, prec)) return NULL; |
|
for (j=1; j<=k; j++) coeff(L,j,j) = B[j]; |
|
return gtrans_i(L); |
|
} |
|
|
|
GEN |
|
R_from_QR(GEN x, long prec) |
|
{ |
|
long j, k = lg(x)-1; |
|
GEN L, B = zerovec(k), Q = cgetg(k+1, t_VEC); |
|
L = cgetg(k+1, t_MAT); |
|
for (j=1; j<=k; j++) L[j] = (long)zerocol(k); |
|
for (j=1; j<=k; j++) |
|
if (!incrementalQ(x, L, B, Q, j, prec)) return NULL; |
|
return gtrans_i(L); |
|
} |
|
|
|
/********************************************************************/ |
|
/** QR Factorization via Gram-Schmidt **/ |
|
/********************************************************************/ |
|
|
|
/* compute B[k] = | x[k] |^2, update mu(k, 1..k-1). |
|
* Classical Gram-Schmidt (unstable!) */ |
|
static int |
|
incrementalGS(GEN x, GEN mu, GEN B, long k) |
|
{ |
|
GEN s,A = cgetg(k+1, t_COL); /* scratch space */ |
|
long i, j; |
|
gpmem_t av; |
|
A[1] = coeff(x,k,1); |
|
for(j=1;j<k;) |
|
{ |
|
coeff(mu,k,j) = lmpdiv((GEN)A[j], (GEN)B[j]); |
|
j++; av = avma; |
|
/* A[j] <-- x[k,j] - sum_{i<j} mu[j,i] A[i] */ |
|
s = mpmul(gcoeff(mu,j,1),(GEN)A[1]); |
|
for (i=2; i<j; i++) s = mpadd(s, mpmul(gcoeff(mu,j,i),(GEN)A[i])); |
|
s = mpneg(s); A[j] = (long)gerepileuptoleaf(av, mpadd(gcoeff(x,k,j), s)); |
|
} |
|
B[k] = A[k]; return (signe((GEN)B[k]) > 0); |
|
} |
|
|
|
#if 0 |
|
/* return Gram-Schmidt orthogonal basis (f) associated to (e), B is the |
|
* vector of the (f_i . f_i)*/ |
|
GEN |
|
gram_schmidt(GEN e, GEN *ptB) |
|
{ |
|
long i,j,lx = lg(e); |
|
GEN f = dummycopy(e), B, iB; |
|
|
|
B = cgetg(lx, t_VEC); |
|
iB= cgetg(lx, t_VEC); |
|
|
|
for (i=1;i<lx;i++) |
|
{ |
|
GEN p1 = NULL; |
|
gpmem_t av; |
|
B[i] = (long)sqscal((GEN)f[i]); |
|
iB[i]= linv((GEN)B[i]); av = avma; |
|
for (j=1; j<i; j++) |
|
{ |
|
GEN mu = gmul(gscal((GEN)e[i],(GEN)f[j]), (GEN)iB[j]); |
|
GEN p2 = gmul(mu, (GEN)f[j]); |
|
p1 = p1? gadd(p1,p2): p2; |
|
} |
|
p1 = p1? gerepileupto(av, gsub((GEN)e[i], p1)): (GEN)e[i]; |
|
f[i] = (long)p1; |
|
} |
|
*ptB = B; return f; |
|
} |
|
#endif |
|
|
|
/********************************************************************/ |
|
/** **/ |
|
/** LLL ALGORITHM **/ |
|
/** **/ |
|
/********************************************************************/ |
|
#define swap(x,y) { long _t=x; x=y; y=_t; } |
|
#define gswap(x,y) { GEN _t=x; x=y; y=_t; } |
|
|
|
static GEN |
|
lll_trivial(GEN x, long flag) |
|
{ |
|
GEN y; |
|
if (lg(x) == 1) |
|
{ /* dim x = 0 */ |
|
if (! (flag & lll_ALL)) return cgetg(1,t_MAT); |
y=cgetg(3,t_VEC); |
y=cgetg(3,t_VEC); |
y[1]=lgetg(1,t_MAT); |
y[1]=lgetg(1,t_MAT); |
y[2]=lgetg(1,t_MAT); return y; |
y[2]=lgetg(1,t_MAT); return y; |
} |
} |
/* here n = 1 */ |
/* here dim = 1 */ |
if (gcmp0((GEN)x[1])) |
if (gcmp0((GEN)x[1])) |
{ |
{ |
switch(flag ^ lll_GRAM) |
switch(flag & (~lll_GRAM)) |
{ |
{ |
case lll_KER: return idmat(1); |
case lll_KER: return idmat(1); |
case lll_IM : return cgetg(1,t_MAT); |
case lll_IM : return cgetg(1,t_MAT); |
Line 120 lllall_trivial(GEN x, long n, long flag) |
|
Line 315 lllall_trivial(GEN x, long n, long flag) |
|
} |
} |
|
|
static GEN |
static GEN |
lllgramall_finish(GEN h,GEN fl,long flag,long n) |
lll_finish(GEN h,GEN fl,long flag) |
{ |
{ |
long k; |
long i, k, l = lg(fl); |
GEN y; |
GEN y, g; |
|
|
k=1; while (k<=n && !fl[k]) k++; |
k=1; while (k<l && !fl[k]) k++; |
switch(flag) |
switch(flag & (~lll_GRAM)) |
{ |
{ |
case lll_KER: setlg(h,k); |
case lll_KER: setlg(h,k); |
y = gcopy(h); break; |
return h; |
|
|
case lll_IM: h += k-1; h[0] = evaltyp(t_MAT)| evallg(n-k+2); |
case lll_IM: h += k-1; h[0] = evaltyp(t_MAT) | evallg(l-k+1); |
y = gcopy(h); break; |
return h; |
|
|
default: setlg(h,k); y=cgetg(3,t_VEC); |
|
y[1] = lcopy(h); |
|
h += k-1; h[0] = evaltyp(t_MAT)| evallg(n-k+2); |
|
y[2] = lcopy(h); |
|
break; |
|
} |
} |
return y; |
y = cgetg(3,t_VEC); |
|
g = cgetg(k, t_MAT); for (i=1; i<k; i++) g[i] = h[i]; |
|
y[1] = (long)g; |
|
h += k-1; h[0] = evaltyp(t_MAT) | evallg(l-k+1); |
|
y[2] = (long)h; return y; |
} |
} |
|
|
/********************************************************************/ |
/* h[,k] += q * h[,l] */ |
/** **/ |
static void |
/** LLL with content **/ |
Zupdate_col(long k, long l, GEN q, long K, GEN h) |
/** **/ |
|
/********************************************************************/ |
|
|
|
/* real Gram matrix has coeffs X[i,j] = x[i,j]*veccon[i]*veccon[j] */ |
|
static GEN |
|
lllgramintwithcontent(GEN x, GEN veccon, long flag) |
|
{ |
{ |
long av0=avma,av,tetpil,lx=lg(x),i,j,k,l,n,lim,kmax; |
GEN *hl, *hk; |
GEN u,u2,B,lam,q,r,h,la,bb,p1,p2,p3,p4,fl,corr,corr2,newcon; |
long i; |
GEN *gptr[8]; |
|
|
|
if (typ(x) != t_MAT) err(typeer,"lllgramintwithcontent"); |
if (!h) return; |
n=lx-1; if (n<=1) return lllall_trivial(x,n,flag); |
hl = (GEN*)h[l]; hk = (GEN*)h[k]; |
if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramintwithcontent"); |
if (is_pm1(q)) |
fl = new_chunk(lx); |
{ |
|
if (signe(q) > 0) |
|
for (i=1;i<=K;i++) { if (signe(hl[i])) hk[i] = addii(hk[i],hl[i]); } |
|
else |
|
for (i=1;i<=K;i++) { if (signe(hl[i])) hk[i] = subii(hk[i],hl[i]); } |
|
} else |
|
for (i=1;i<=K;i++) if (signe(hl[i])) hk[i] = addii(hk[i],mulii(q,hl[i])); |
|
} |
|
|
av=avma; lim=stack_lim(av,1); |
/* L[k,] += q * L[l,], l < k */ |
x=dummycopy(x); veccon=dummycopy(veccon); |
static void |
B=cgetg(lx+1,t_COL); B[1]=un; |
Zupdate_row(long k, long l, GEN q, GEN L, GEN B) |
/* B[i+1]=B_i; le vrai B_i est B_i*prod(1,j=1,i,veccon[j]^2) */ |
{ |
|
long i; |
for (i=1; i<lx; i++) { B[i+1]=zero; fl[i]=0; } |
if (is_pm1(q)) |
lam=cgetg(lx,t_MAT); |
|
for (j=1; j<lx; j++) |
|
{ p2=cgetg(lx,t_COL); lam[j]=(long)p2; for (i=1; i<lx; i++) p2[i]=zero; } |
|
/* le vrai lam[i,j] est |
|
lam[i,j]*veccon[i]*veccon[j]*prod(1,l=1,j-1,veccon[l]^2) */ |
|
k=2; h=idmat(n); kmax=1; |
|
u=gcoeff(x,1,1); if (typ(u)!=t_INT) err(lllger4); |
|
if (signe(u)) { B[2]=(long)u; coeff(lam,1,1)=un; fl[1]=1; } |
|
else { B[2]=un; fl[1]=0; } |
|
if (DEBUGLEVEL>5) { fprintferr("k = %ld",k); flusherr(); } |
|
for(;;) |
|
{ |
{ |
if (k>kmax) |
if (signe(q) > 0) |
{ |
{ |
kmax=k; |
for (i=1;i<l; i++) coeff(L,k,i) = laddii(gcoeff(L,k,i),gcoeff(L,l,i)); |
for (j=1; j<=k; j++) |
coeff(L,k,l) = laddii(gcoeff(L,k,l), B); |
{ |
} else { |
if (j==k || fl[j]) |
for (i=1;i<l; i++) coeff(L,k,i) = lsubii(gcoeff(L,k,i),gcoeff(L,l,i)); |
{ |
coeff(L,k,l) = laddii(gcoeff(L,k,l), negi(B)); |
u=gcoeff(x,k,j); if (typ(u)!=t_INT) err(lllger4); |
|
for (i=1; i<j; i++) |
|
if (fl[i]) |
|
u=divii(subii(mulii((GEN)B[i+1],u),mulii(gcoeff(lam,k,i),gcoeff(lam,j,i))),(GEN)B[i]); |
|
if (j<k) coeff(lam,k,j)=(long)u; |
|
else |
|
{ |
|
if (signe(u)) { B[k+1]=(long)u; coeff(lam,k,k)=un; fl[k]=1; } |
|
else { B[k+1]=B[k]; fl[k]=0; } |
|
} |
|
} |
|
} |
|
if (low_stack(lim, stack_lim(av,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"[1]: lllgramintwithcontent"); |
|
gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h; |
|
gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5); |
|
} |
|
} |
} |
u=shifti(mulii(gcoeff(lam,k,k-1),(GEN)veccon[k]),1); |
} else { |
u2=mulii((GEN)B[k],(GEN)veccon[k-1]); |
for(i=1;i<l;i++) coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i))); |
if (cmpii(absi(u),u2)>0) |
coeff(L,k,l) = laddii(gcoeff(L,k,l), mulii(q,B)); |
{ |
|
q=dvmdii(addii(u,u2),shifti(u2,1),&r); |
|
if (signe(r)<0) q=addsi(-1,q); |
|
h[k]=lsub((GEN)h[k],gmul(q,(GEN)h[k-1])); |
|
newcon=mppgcd((GEN)veccon[k],(GEN)veccon[k-1]); |
|
corr=divii((GEN)veccon[k],newcon); veccon[k]=(long)newcon; |
|
if(!gcmp1(corr)) |
|
{ |
|
corr2=sqri(corr); |
|
for (j=1; j<=n; j++) |
|
coeff(x,j,k)=coeff(x,k,j)=lmulii(corr,gcoeff(x,k,j)); |
|
coeff(x,k,k)=lmulii(corr,gcoeff(x,k,k)); |
|
for (j=k; j<=kmax; j++) B[j+1]=lmulii(corr2,(GEN)B[j+1]); |
|
for (i=1; i<k; i++) coeff(lam,k,i)=lmulii(corr,gcoeff(lam,k,i)); |
|
for (i=k+1; i<=kmax; i++) |
|
{ |
|
coeff(lam,i,k)=lmulii(corr,gcoeff(lam,i,k)); |
|
for (j=k+1; j<i; j++) |
|
coeff(lam,i,j)=lmulii(corr2,gcoeff(lam,i,j)); |
|
} |
|
} |
|
r=negi(mulii(q,divii((GEN)veccon[k-1],(GEN)veccon[k]))); |
|
p1=gcoeff(x,k,k-1); |
|
for (j=1; j<=n; j++) |
|
coeff(x,j,k)=coeff(x,k,j)=laddii(gcoeff(x,j,k),mulii(r,(j==k)?p1:gcoeff(x,j,k-1))); |
|
coeff(x,k,k)=laddii(gcoeff(x,k,k),mulii(r,gcoeff(x,k-1,k))); |
|
coeff(lam,k,k-1)=laddii(gcoeff(lam,k,k-1),mulii(r,(GEN)B[k])); |
|
for (i=1; i<k-1; i++) |
|
coeff(lam,k,i)=laddii(gcoeff(lam,k,i),mulii(r,gcoeff(lam,k-1,i))); |
|
} |
|
if (low_stack(lim, stack_lim(av,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"[2]: lllgramintwithcontent"); |
|
gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h; |
|
gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5); |
|
} |
|
p3=mulii((GEN)B[k-1],(GEN)B[k+1]);la=gcoeff(lam,k,k-1);p4=mulii(la,la); |
|
p2=mulsi(100,mulii(mulii((GEN)veccon[k],(GEN)veccon[k]),addii(p3,p4))); |
|
p1=mulii((GEN)veccon[k-1],(GEN)B[k]);p1=mulsi(99,mulii(p1,p1)); |
|
if (fl[k-1] && (cmpii(p1,p2)>0 || !fl[k])) |
|
{ |
|
if (DEBUGLEVEL>=4 && k==n) |
|
{ fprintferr(" (%ld)", expi(p1)-expi(p2)); flusherr(); } |
|
p1=(GEN)h[k-1]; h[k-1]=h[k]; h[k]=(long)p1; |
|
p1=(GEN)x[k-1]; x[k-1]=x[k]; x[k]=(long)p1; |
|
for (j=1; j<=n; j++) |
|
{ p1=gcoeff(x,k-1,j); coeff(x,k-1,j)=coeff(x,k,j); coeff(x,k,j)=(long)p1; } |
|
p1=(GEN)veccon[k-1]; veccon[k-1]=veccon[k]; veccon[k]=(long)p1; |
|
for (j=1; j<=k-2; j++) |
|
{ p1=gcoeff(lam,k-1,j); coeff(lam,k-1,j)=coeff(lam,k,j); coeff(lam,k,j)=(long)p1; } |
|
if (fl[k]) |
|
{ |
|
for (i=k+1; i<=kmax; i++) |
|
{ |
|
bb=gcoeff(lam,i,k); |
|
coeff(lam,i,k)=ldivii(subii(mulii((GEN)B[k+1],gcoeff(lam,i,k-1)),mulii(la,bb)),(GEN)B[k]); |
|
coeff(lam,i,k-1)=ldivii(addii(mulii(la,gcoeff(lam,i,k-1)),mulii((GEN)B[k-1],bb)),(GEN)B[k]); |
|
if (low_stack(lim, stack_lim(av,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"[3]: lllgramintwithcontent"); |
|
gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h; |
|
gptr[3]=&x; gptr[4]=&la; gptr[5]=&veccon; gptr[6]=&p3; |
|
gptr[7]=&p4; gerepilemany(av,gptr,8); |
|
} |
|
} |
|
B[k]=ldivii(addii(p3,p4),(GEN)B[k]); |
|
} |
|
else |
|
{ |
|
if (signe(la)) |
|
{ |
|
p2=(GEN)B[k]; p1=divii(p4,p2); |
|
for (i=k+1; i<=kmax; i++) |
|
coeff(lam,i,k-1)=ldivii(mulii(la,gcoeff(lam,i,k-1)),p2); |
|
for (j=k+1; j<kmax; j++) |
|
{ |
|
for (i=j+1; i<=kmax; i++) |
|
coeff(lam,i,j)=ldivii(mulii(p1,gcoeff(lam,i,j)),p2); |
|
if (low_stack(lim, stack_lim(av,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"[4]: lllgramintwithcontent"); |
|
gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h; |
|
gptr[3]=&x; gptr[4]=&la; gptr[5]=&veccon; gptr[6]=&p1; |
|
gptr[7]=&p2; gerepilemany(av,gptr,8); |
|
} |
|
} |
|
B[k+1]=B[k]=(long)p1; |
|
for (i=k+2; i<=lx; i++) |
|
B[i]=ldivii(mulii(p1,(GEN)B[i]),p2); |
|
} |
|
else |
|
{ |
|
coeff(lam,k,k-1)=zero; |
|
for (i=k+1; i<=kmax; i++) |
|
{ |
|
coeff(lam,i,k)=coeff(lam,i,k-1); |
|
coeff(lam,i,k-1)=zero; |
|
} |
|
B[k]=B[k-1]; fl[k]=1; fl[k-1]=0; |
|
} |
|
|
|
if (low_stack(lim, stack_lim(av,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"[5]: lllgramintwithcontent"); |
|
gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h; |
|
gptr[3]=&x; gptr[4]=&veccon; |
|
gerepilemany(av,gptr,5); |
|
} |
|
} |
|
if (k>2) k--; |
|
if (DEBUGLEVEL>5) { fprintferr(" %ld",k); flusherr(); } |
|
} |
|
else |
|
{ |
|
for (l=k-2; l>=1; l--) |
|
{ |
|
u=shifti(mulii(gcoeff(lam,k,l),(GEN)veccon[k]),1); |
|
u2=mulii((GEN)B[l+1],(GEN)veccon[l]); |
|
if (cmpii(absi(u),u2)>0) |
|
{ |
|
q=dvmdii(addii(u,u2),shifti(u2,1),&r); |
|
if (signe(r)<0) q=addsi(-1,q); |
|
h[k]=lsub((GEN)h[k],gmul(q,(GEN)h[l])); |
|
newcon=mppgcd((GEN)veccon[k],(GEN)veccon[l]); |
|
corr=divii((GEN)veccon[k],newcon); veccon[k]=(long)newcon; |
|
if(!gcmp1(corr)) |
|
{ |
|
corr2=sqri(corr); |
|
for (j=1; j<=n; j++) |
|
coeff(x,j,k)=coeff(x,k,j)=lmulii(corr,gcoeff(x,k,j)); |
|
coeff(x,k,k)=lmulii(corr,gcoeff(x,k,k)); |
|
for (j=k; j<=kmax; j++) B[j+1]=lmulii(corr2,(GEN)B[j+1]); |
|
for (i=1; i<k; i++) coeff(lam,k,i)=lmulii(corr,gcoeff(lam,k,i)); |
|
for (i=k+1; i<=kmax; i++) |
|
{ |
|
coeff(lam,i,k)=lmulii(corr,gcoeff(lam,i,k)); |
|
for (j=k+1; j<i; j++) |
|
coeff(lam,i,j)=lmulii(corr2,gcoeff(lam,i,j)); |
|
} |
|
} |
|
r=negi(mulii(q,divii((GEN)veccon[l],(GEN)veccon[k]))); |
|
p1=gcoeff(x,k,l); |
|
for (j=1; j<=n; j++) |
|
coeff(x,j,k)=coeff(x,k,j)=laddii(gcoeff(x,j,k),mulii(r,(j==k)?p1:gcoeff(x,j,l))); |
|
coeff(x,k,k)=laddii(gcoeff(x,k,k),mulii(r,gcoeff(x,l,k))); |
|
coeff(lam,k,l)=laddii(gcoeff(lam,k,l),mulii(r,(GEN)B[l+1])); |
|
for (i=1; i<l; i++) |
|
coeff(lam,k,i)=laddii(gcoeff(lam,k,i),mulii(r,gcoeff(lam,l,i))); |
|
} |
|
if (low_stack(lim, stack_lim(av,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"[6]: lllgramintwithcontent"); |
|
gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h; |
|
gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5); |
|
} |
|
} |
|
k++; if (DEBUGLEVEL>5) { fprintferr(" %ld",k); flusherr(); } |
|
if (k>n) |
|
{ |
|
if (DEBUGLEVEL>5) { fprintferr("\n"); flusherr(); } |
|
tetpil=avma; |
|
return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n)); |
|
} |
|
} |
|
if (low_stack(lim, stack_lim(av,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"[7]: lllgramintwithcontent"); |
|
gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h; |
|
gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5); |
|
} |
|
} |
} |
} |
} |
|
|
static GEN |
static void |
lllintwithcontent(GEN x) |
update_row(long k, long l, GEN q, GEN L) |
{ |
{ |
long lx=lg(x),i,j,av,tetpil; |
long i; |
GEN veccon,con,xred,g; |
if (is_pm1(q)) |
|
|
if (typ(x) != t_MAT) err(typeer,"lllintwithcontent"); |
|
if (lx==1) return cgetg(1,t_MAT); |
|
av=avma; veccon=cgetg(lx,t_VEC); g=cgetg(lx,t_MAT); xred=cgetg(lx,t_MAT); |
|
for (j=1; j<lx; j++) |
|
{ |
{ |
g[j]=lgetg(lx,t_COL); con=content((GEN)x[j]); |
if (signe(q) > 0) |
xred[j]=ldiv((GEN)x[j],con); veccon[j]=(long)con; |
{ |
|
for (i=1;i<l; i++) coeff(L,k,i) = lmpadd(gcoeff(L,k,i),gcoeff(L,l,i)); |
|
} else { |
|
for (i=1;i<l; i++) coeff(L,k,i) = lmpsub(gcoeff(L,k,i),gcoeff(L,l,i)); |
|
} |
|
} else { |
|
for(i=1;i<l;i++) coeff(L,k,i)=lmpadd(gcoeff(L,k,i),mpmul(q,gcoeff(L,l,i))); |
} |
} |
for (i=1; i<lx; i++) |
coeff(L,k,l) = lmpadd(gcoeff(L,k,l), q); |
for (j=1; j<=i; j++) |
|
coeff(g,i,j)=coeff(g,j,i)=(long)gscal((GEN)xred[i],(GEN)xred[j]); |
|
tetpil=avma; |
|
return gerepile(av,tetpil,lllgramintwithcontent(g,veccon,2)); |
|
} |
} |
|
|
/********************************************************************/ |
|
/** **/ |
|
/** LLL ALGORITHM **/ |
|
/** **/ |
|
/********************************************************************/ |
|
#define swap(x,y) { long _t=x; x=y; y=_t; } |
|
#define gswap(x,y) { GEN _t=x; x=y; y=_t; } |
|
|
|
static void |
static void |
REDI(GEN x, GEN h, GEN L, GEN B, long K, long k, long l) |
ZRED_gram(long k, long l, GEN x, GEN h, GEN L, GEN B, long K) |
{ |
{ |
long i,lx; |
long i,lx; |
GEN q = truedvmdii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1), NULL); |
GEN q = truedvmdii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1), NULL); |
GEN *hk,*hl,xk,xl; |
GEN xk,xl; |
if (!signe(q)) return; |
if (!signe(q)) return; |
q = negi(q); lx = lg(x); |
q = negi(q); |
xl = (GEN)x[l]; hl = (GEN*)h[l]; |
xl = (GEN) x[l]; xk = (GEN) x[k]; |
xk = (GEN)x[k]; hk = (GEN*)h[k]; |
lx = lg(xl); |
if (is_pm1(q)) |
if (is_pm1(q)) |
{ |
{ |
if (signe(q) > 0) |
if (signe(q) > 0) |
{ |
{ |
for (i=1;i<=K;i++) hk[i]=addii(hk[i],hl[i]); |
|
xk[k] = laddii((GEN)xk[k], (GEN)xl[k]); |
xk[k] = laddii((GEN)xk[k], (GEN)xl[k]); |
for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=laddii((GEN)xk[i], (GEN)xl[i]); |
for (i=1;i<lx;i++) coeff(x,k,i) = xk[i] = laddii((GEN)xk[i], (GEN)xl[i]); |
for (i=1;i<l; i++) coeff(L,k,i)=laddii(gcoeff(L,k,i),gcoeff(L,l,i)); |
|
q = B; |
|
} else { |
} else { |
for (i=1;i<=K;i++) hk[i]=subii(hk[i],hl[i]); |
|
xk[k] = lsubii((GEN)xk[k], (GEN)xl[k]); |
xk[k] = lsubii((GEN)xk[k], (GEN)xl[k]); |
for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=lsubii((GEN)xk[i], (GEN)xl[i]); |
for (i=1;i<lx;i++) coeff(x,k,i) = xk[i] = lsubii((GEN)xk[i], (GEN)xl[i]); |
for (i=1;i<l; i++) coeff(L,k,i)=lsubii(gcoeff(L,k,i),gcoeff(L,l,i)); |
|
q = negi(B); |
|
} |
} |
} else { |
} else { /* h[,k] += q* h[,l]. x[,k] += q * x[,l]. L[k,] += q* L[l,] */ |
for(i=1;i<=K;i++) hk[i]=addii(hk[i],mulii(q,hl[i])); |
|
xk[k] = laddii((GEN)xk[k], mulii(q,(GEN)xl[k])); |
xk[k] = laddii((GEN)xk[k], mulii(q,(GEN)xl[k])); |
for(i=1;i<lx;i++) coeff(x,k,i)=xk[i]=laddii((GEN)xk[i],mulii(q,(GEN)xl[i])); |
for(i=1;i<lx;i++) coeff(x,k,i)=xk[i]=laddii((GEN)xk[i],mulii(q,(GEN)xl[i])); |
for(i=1;i<l;i++) coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i))); |
|
q = mulii(q,B); |
|
} |
} |
coeff(L,k,l) = laddii(gcoeff(L,k,l), q); |
Zupdate_row(k,l,q,L,B); |
|
Zupdate_col (k,l,q,K,h); |
} |
} |
|
|
/* b[k] <-- b[k] - round(L[k,l]) b[l], only b[1] ... b[K] modified so far |
|
* assume l < k and update x = Gram(b), L = Gram Schmidt coeffs. */ |
|
static void |
static void |
RED(GEN x, GEN h, GEN L, long K, long k, long l) |
ZRED(long k, long l, GEN x, GEN h, GEN L, GEN B, long K) |
{ |
{ |
long e,i,lx; |
GEN q = truedvmdii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1), NULL); |
GEN q = grndtoi(gcoeff(L,k,l),&e); |
|
GEN *hk,*hl,xk,xl; |
|
if (DEBUGLEVEL>8) |
|
fprintferr("error bits when rounding in lllgram: %ld\n",e); |
|
if (!signe(q)) return; |
if (!signe(q)) return; |
|
q = negi(q); |
|
Zupdate_row(k,l,q,L,B); |
|
Zupdate_col (k,l,q,K,h); |
|
x[k] = (long)ZV_lincomb(gun, q, (GEN)x[k], (GEN)x[l]); |
|
} |
|
|
|
static GEN |
|
round_safe(GEN q) |
|
{ |
|
long e; |
|
if (typ(q) == t_INT) return q; |
|
if (expo(q) > 30) |
|
{ |
|
q = grndtoi(q, &e); |
|
if (e > 0) return NULL; |
|
} else |
|
q = ground(q); |
|
return q; |
|
} |
|
|
|
/* b[k] <-- b[k] - round(L[k,l]) b[l], only b[1] ... b[K] modified so far |
|
* assume l < k and update x = Gram(b), L = Gram Schmidt coeffs. */ |
|
static int |
|
RED_gram(long k, long l, GEN x, GEN h, GEN L, long K) |
|
{ |
|
long i,lx; |
|
GEN q = round_safe(gcoeff(L,k,l)); |
|
GEN xk, xl; |
|
|
|
if (!q) return 0; |
|
if (!signe(q)) return 1; |
q = negi(q); lx = lg(x); |
q = negi(q); lx = lg(x); |
xl = (GEN)x[l]; hl = (GEN*)h[l]; |
xk = (GEN)x[k]; xl = (GEN)x[l]; |
xk = (GEN)x[k]; hk = (GEN*)h[k]; |
|
if (is_pm1(q)) |
if (is_pm1(q)) |
{ |
{ |
if (signe(q) > 0) |
if (signe(q) > 0) |
{ |
{ |
for (i=1;i<=K;i++) hk[i]=addii(hk[i],hl[i]); |
xk[k] = lmpadd((GEN)xk[k], (GEN)xl[k]); |
xk[k] = ladd((GEN)xk[k], (GEN)xl[k]); |
for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=lmpadd((GEN)xk[i], (GEN)xl[i]); |
for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=ladd((GEN)xk[i], (GEN)xl[i]); |
|
for (i=1;i<l; i++) coeff(L,k,i)=ladd(gcoeff(L,k,i),gcoeff(L,l,i)); |
|
} else { |
} else { |
for (i=1;i<=K;i++) hk[i]=subii(hk[i],hl[i]); |
xk[k] = lmpsub((GEN)xk[k], (GEN)xl[k]); |
xk[k] = lsub((GEN)xk[k], (GEN)xl[k]); |
for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=lmpsub((GEN)xk[i], (GEN)xl[i]); |
for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=lsub((GEN)xk[i], (GEN)xl[i]); |
|
for (i=1;i<l; i++) coeff(L,k,i)=lsub(gcoeff(L,k,i),gcoeff(L,l,i)); |
|
} |
} |
} else { |
} else { |
for (i=1;i<=K;i++) hk[i]=addii(hk[i],mulii(q,hl[i])); |
xk[k] = lmpadd((GEN)xk[k], mpmul(q,(GEN)xl[k])); |
xk[k] = ladd((GEN)xk[k], gmul(q,(GEN)xl[k])); |
for (i=1;i<lx;i++) |
for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=ladd((GEN)xk[i], gmul(q,(GEN)xl[i])); |
coeff(x,k,i)=xk[i]=lmpadd((GEN)xk[i], mpmul(q,(GEN)xl[i])); |
for (i=1;i<l; i++) coeff(L,k,i)=ladd(gcoeff(L,k,i),gmul(q,gcoeff(L,l,i))); |
|
} |
} |
coeff(L,k,l) = ladd(gcoeff(L,k,l),q); |
update_row(k,l,q,L); |
|
Zupdate_col(k,l,q,K,h); return 1; |
} |
} |
|
|
|
static void |
|
update_col(long k, long l, GEN q, GEN x) |
|
{ |
|
long i,lx; |
|
GEN xk = (GEN)x[k], xl = (GEN)x[l]; |
|
lx = lg(xk); |
|
if (is_pm1(q)) |
|
{ |
|
if (signe(q) > 0) |
|
{ |
|
for (i=1;i<lx;i++) xk[i] = lmpadd((GEN)xk[i], (GEN)xl[i]); |
|
} else { |
|
for (i=1;i<lx;i++) xk[i] = lmpsub((GEN)xk[i], (GEN)xl[i]); |
|
} |
|
} else { |
|
for (i=1;i<lx;i++) xk[i] = lmpadd((GEN)xk[i], mpmul(q,(GEN)xl[i])); |
|
} |
|
} |
|
|
static int |
static int |
do_SWAPI(GEN x, GEN h, GEN L, GEN B, long kmax, long k, long alpha, GEN fl) |
RED(long k, long l, GEN x, GEN h, GEN L, long K) |
{ |
{ |
GEN la,la2,p1,p2,Bk; |
GEN q = round_safe(gcoeff(L,k,l)); |
long av,i,j,lx; |
if (!q) return 0; |
|
if (!signe(q)) return 1; |
|
q = negi(q); |
|
update_col(k,l,q,x); |
|
update_row(k,l,q,L); |
|
Zupdate_col(k,l,q,K,h); return 1; |
|
} |
|
|
|
static int |
|
do_ZSWAP(GEN x, GEN h, GEN L, GEN B, long kmax, long k, long D, GEN fl, |
|
int gram) |
|
{ |
|
GEN la,la2,p1,Bk; |
|
long i, j, lx; |
|
gpmem_t av; |
|
|
if (!fl[k-1]) return 0; |
if (!fl[k-1]) return 0; |
lx = lg(x); av = avma; |
av = avma; |
la = gcoeff(L,k,k-1); la2 = sqri(la); |
la = gcoeff(L,k,k-1); la2 = sqri(la); |
p2 = addii(la2, mulii((GEN)B[k-1],(GEN)B[k+1])); |
|
Bk = (GEN)B[k]; |
Bk = (GEN)B[k]; |
if (fl[k] && cmpii(mulsi(alpha-1,sqri(Bk)), mulsi(alpha,p2)) <= 0) |
if (fl[k]) |
{ avma = av; return 0; } |
{ |
|
GEN q; |
/* SWAPI(k-1,k) */ |
if (!D) return 0; /* only swap non-kernel + kernel vector */ |
if (DEBUGLEVEL>3 && k==kmax) { |
q = addii(la2, mulii((GEN)B[k-1],(GEN)B[k+1])); |
fprintferr(" (%ld)", expi(mulsi(alpha-1,sqri(Bk))) |
if (cmpii(mulsi(D-1,sqri(Bk)), mulsi(D,q)) <= 0) { |
- expi(mulsi(alpha,p2))); flusherr(); |
avma = av; return 0; |
|
} |
|
B[k] = (long)diviiexact(q, Bk); |
} |
} |
swap(h[k-1], h[k]); |
/* ZSWAP(k-1,k) */ |
swap(x[k-1], x[k]); |
if (DEBUGLEVEL>3 && k==kmax) |
for (j=1; j< lx; j++) swap(coeff(x,k-1,j), coeff(x,k,j)); |
{ /* output diagnostics associated to re-normalized rational quantities */ |
|
gpmem_t av1 = avma; |
|
GEN d = mulii((GEN)B[k-1],(GEN)B[k+1]); |
|
p1 = subii(mulsi(D-1, sqri(Bk)), mulsi(D, la2)); |
|
fprintferr(" (%ld)", expi(p1) - expi(mulsi(D, d))); |
|
avma = av1; |
|
} |
|
if (h) swap(h[k-1], h[k]); |
|
swap(x[k-1], x[k]); lx = lg(x); |
|
if (gram) |
|
for (j=1; j < lx; j++) swap(coeff(x,k-1,j), coeff(x,k,j)); |
for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j)) |
for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j)) |
if (fl[k]) |
if (fl[k]) |
{ |
{ |
Line 514 do_SWAPI(GEN x, GEN h, GEN L, GEN B, long kmax, long k |
|
Line 553 do_SWAPI(GEN x, GEN h, GEN L, GEN B, long kmax, long k |
|
{ |
{ |
GEN t = gcoeff(L,i,k); |
GEN t = gcoeff(L,i,k); |
p1 = subii(mulii((GEN)B[k+1],gcoeff(L,i,k-1)), mulii(la,t)); |
p1 = subii(mulii((GEN)B[k+1],gcoeff(L,i,k-1)), mulii(la,t)); |
p1 = divii(p1, Bk); |
p1 = diviiexact(p1, Bk); |
av = avma = coeff(L,i,k) = (long)icopy_av(p1,(GEN)av); |
coeff(L,i,k) = (long)icopy_av(p1,(GEN)av); |
|
av = avma = (gpmem_t)coeff(L,i,k); |
|
|
p1 = addii(mulii(la,gcoeff(L,i,k-1)), mulii((GEN)B[k-1],t)); |
p1 = addii(mulii(la,gcoeff(L,i,k-1)), mulii((GEN)B[k-1],t)); |
p1 = divii(p1, Bk); |
p1 = diviiexact(p1, Bk); |
av = avma = coeff(L,i,k-1) = (long)icopy_av(p1,(GEN)av); |
coeff(L,i,k-1) = (long)icopy_av(p1,(GEN)av); |
|
av = avma = (gpmem_t)coeff(L,i,k-1); |
} |
} |
B[k] = ldivii(p2,Bk); |
|
} |
} |
|
else if (signe(la)) |
|
{ |
|
p1 = diviiexact(la2, Bk); |
|
B[k+1] = B[k] = (long)p1; |
|
for (i=k+2; i<=lx; i++) B[i] = (long)diviiexact(mulii(p1,(GEN)B[i]), Bk); |
|
for (i=k+1; i<=kmax; i++) |
|
coeff(L,i,k-1) = (long)diviiexact(mulii(la,gcoeff(L,i,k-1)), Bk); |
|
for (j=k+1; j<kmax; j++) |
|
for (i=j+1; i<=kmax; i++) |
|
coeff(L,i,j) = (long)diviiexact(mulii(p1,gcoeff(L,i,j)), Bk); |
|
} |
else |
else |
{ |
{ |
if (signe(la)) |
for (i=k+1; i<=kmax; i++) |
{ |
{ |
p1 = divii(la2, Bk); |
coeff(L,i,k) = coeff(L,i,k-1); |
B[k+1] = B[k] = (long)p1; |
coeff(L,i,k-1) = zero; |
for (i=k+2; i<=lx; i++) B[i] = ldivii(mulii(p1,(GEN)B[i]), Bk); |
|
for (i=k+1; i<=kmax; i++) |
|
coeff(L,i,k-1) = ldivii(mulii(la,gcoeff(L,i,k-1)), Bk); |
|
for (j=k+1; j<kmax; j++) |
|
for (i=j+1; i<=kmax; i++) |
|
coeff(L,i,j) = ldivii(mulii(p1,gcoeff(L,i,j)), Bk); |
|
} |
} |
else |
B[k] = B[k-1]; fl[k] = 1; fl[k-1] = 0; |
{ |
|
for (i=k+1; i<=kmax; i++) |
|
{ |
|
coeff(L,i,k) = coeff(L,i,k-1); |
|
coeff(L,i,k-1) = zero; |
|
} |
|
B[k] = B[k-1]; fl[k] = 1; fl[k-1] = 0; |
|
} |
|
} |
} |
return 1; |
return 1; |
} |
} |
|
|
static int |
static int |
do_SWAP(GEN x, GEN h, GEN L, GEN B, long kmax, long k, GEN QR) |
do_SWAP(GEN x, GEN h, GEN L, GEN B, long kmax, long k, GEN delta, int gram) |
{ |
{ |
GEN la,la2, BK,BB,q; |
GEN la,la2, BK,BB,q; |
long av,i,j,lx; |
long i, j, lx; |
|
gpmem_t av; |
|
|
lx = lg(x); av = avma; |
av = avma; |
la = gcoeff(L,k,k-1); la2 = gsqr(la); |
la = gcoeff(L,k,k-1); la2 = gsqr(la); |
q = gmul((GEN)B[k-1], gsub(QR,la2)); |
q = mpmul((GEN)B[k-1], mpsub(delta,la2)); |
if (gcmp(q, (GEN)B[k]) <= 0) { avma = av; return 0; } |
if (mpcmp(q, (GEN)B[k]) <= 0) { avma = av; return 0; } |
|
|
/* SWAP(k-1,k) */ |
/* SWAP(k-1,k) */ |
if (DEBUGLEVEL>3 && k==kmax) { |
if (DEBUGLEVEL>3 && k==kmax) { |
fprintferr(" (%ld)", gexpo(q) - gexpo((GEN)B[k])); flusherr(); |
fprintferr(" (%ld)", gexpo(q) - gexpo((GEN)B[k])); flusherr(); |
} |
} |
BB = gadd((GEN)B[k], gmul((GEN)B[k-1],la2)); |
BB = mpadd((GEN)B[k], mpmul((GEN)B[k-1],la2)); |
if (gcmp0(BB)) { B[k] = 0; return 1; } /* precision pb */ |
if (!signe(BB)) { B[k] = 0; return 1; } /* precision pb */ |
|
|
coeff(L,k,k-1) = ldiv(gmul(la,(GEN)B[k-1]), BB); |
coeff(L,k,k-1) = lmpdiv(mpmul(la,(GEN)B[k-1]), BB); |
BK = gdiv((GEN)B[k], BB); |
BK = mpdiv((GEN)B[k], BB); |
B[k] = lmul((GEN)B[k-1], BK); |
B[k] = lmpmul((GEN)B[k-1], BK); |
B[k-1] = (long)BB; |
B[k-1] = (long)BB; |
|
|
swap(h[k-1],h[k]); |
swap(h[k-1],h[k]); |
swap(x[k-1],x[k]); |
swap(x[k-1],x[k]); lx = lg(x); |
for (j=1; j<lx ; j++) swap(coeff(x,k-1,j), coeff(x,k,j)); |
if (gram) |
|
for (j=1; j < lx; j++) swap(coeff(x,k-1,j), coeff(x,k,j)); |
for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j)) |
for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j)) |
for (i=k+1; i<=kmax; i++) |
for (i=k+1; i<=kmax; i++) |
{ |
{ |
GEN t = gcoeff(L,i,k); |
GEN t = gcoeff(L,i,k); |
coeff(L,i,k) = lsub(gcoeff(L,i,k-1), gmul(la,t)); |
coeff(L,i,k) = lmpsub(gcoeff(L,i,k-1), mpmul(la,t)); |
coeff(L,i,k-1) = ladd(gmul(gcoeff(L,k,k-1),gcoeff(L,i,k-1)), gmul(BK,t)); |
coeff(L,i,k-1) = lmpadd(t, mpmul(gcoeff(L,k,k-1), gcoeff(L,i,k))); |
/* = ladd(t, gmul(gcoeff(L,k,k-1), gcoeff(L,i,k))); |
|
* would save one multiplication, but loses precision faster... */ |
|
} |
} |
return 1; |
return 1; |
} |
} |
|
static void |
|
ZincrementalGS(GEN x, GEN L, GEN B, long k, GEN fl, int gram) |
|
{ |
|
GEN u = NULL; /* gcc -Wall */ |
|
long i, j, s; |
|
for (j=1; j<=k; j++) |
|
if (j==k || fl[j]) |
|
{ |
|
gpmem_t av = avma; |
|
u = gram? gcoeff(x,k,j): gscali((GEN)x[k], (GEN)x[j]); |
|
for (i=1; i<j; i++) |
|
if (fl[i]) |
|
{ |
|
u = subii(mulii((GEN)B[i+1], u), mulii(gcoeff(L,k,i), gcoeff(L,j,i))); |
|
u = diviiexact(u, (GEN)B[i]); |
|
} |
|
coeff(L,k,j) = lpileuptoint(av, u); |
|
} |
|
s = signe(u); |
|
if (s == 0) B[k+1] = B[k]; |
|
else |
|
{ |
|
if (s < 0) err(lllger3); |
|
B[k+1] = coeff(L,k,k); coeff(L,k,k) = un; fl[k] = 1; |
|
} |
|
} |
|
|
/* x integer matrix */ |
/* x integer matrix. Beware: this function can return NULL */ |
GEN |
GEN |
lllgramall(GEN x, long alpha, long flag) |
lllint_marked(long MARKED, GEN x, long D, int gram, |
|
GEN *pth, GEN *ptfl, GEN *ptB) |
{ |
{ |
long av0=avma,av,tetpil,lim,lx=lg(x),i,j,k,l,n,s,kmax; |
long lx = lg(x), hx, i, j, k, l, n, kmax; |
GEN u,B,L,h,fl, *gptr[4]; |
gpmem_t av, lim; |
|
GEN B,L,h,fl; |
|
|
if (typ(x) != t_MAT) err(typeer,"lllgramall"); |
if (typ(x) != t_MAT) err(typeer,"lllint"); |
n=lx-1; if (n<=1) return lllall_trivial(x,n,flag); |
n = lx-1; if (n <= 1) return NULL; |
if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramall"); |
hx = lg(x[1]); |
|
if (gram && hx != lx) err(mattype1,"lllint"); |
fl = cgetg(lx,t_VECSMALL); |
fl = cgetg(lx,t_VECSMALL); |
|
|
av=avma; lim=stack_lim(av,1); x=dummycopy(x); |
av = avma; lim = stack_lim(av,1); x = dummycopy(x); |
B=gscalcol(gun, lx); |
B = gscalcol(gun, lx); |
L=cgetg(lx,t_MAT); |
L = cgetg(lx,t_MAT); |
for (j=1; j<lx; j++) |
for (j=1; j<lx; j++) |
{ |
{ |
for (i=1; i<lx; i++) |
for (i=1; i<hx; i++) |
if (typ(gcoeff(x,i,j))!=t_INT) err(lllger4); |
if (typ(gcoeff(x,i,j)) != t_INT) err(lllger4); |
fl[j] = 0; L[j] = (long)zerocol(n); |
fl[j] = 0; L[j] = (long)zerocol(n); |
} |
} |
k=2; h=idmat(n); kmax=1; |
h = pth? idmat(n): NULL; |
u=gcoeff(x,1,1); s= signe(u); |
ZincrementalGS(x, L, B, 1, fl, gram); |
if (s == 0) B[2]=un; |
kmax = 1; |
else |
if (DEBUGLEVEL>5) fprintferr("k = "); |
|
for (k=2;;) |
{ |
{ |
if (s < 0) err(lllger3); |
|
B[2]=(long)u; coeff(L,1,1)=un; fl[1]=1; |
|
} |
|
if (DEBUGLEVEL>5) fprintferr("k ="); |
|
for(;;) |
|
{ |
|
if (k > kmax) |
if (k > kmax) |
{ |
{ |
if (DEBUGLEVEL>3) {fprintferr(" K%ld",k);flusherr();} |
if (DEBUGLEVEL>3) fprintferr("K%ld ",k); |
|
ZincrementalGS(x, L, B, k, fl, gram); |
kmax = k; |
kmax = k; |
for (j=1; j<=k; j++) |
} |
if (j==k || fl[j]) |
if (k != MARKED) |
{ |
|
long av1 = avma; |
|
u = gcoeff(x,k,j); |
|
for (i=1; i<j; i++) if (fl[i]) |
|
u = divii(subii(mulii((GEN)B[i+1],u), |
|
mulii(gcoeff(L,k,i),gcoeff(L,j,i))), |
|
(GEN)B[i]); |
|
u = gerepileuptoint(av1, u); |
|
if (j<k) coeff(L,k,j)=(long)u; |
|
else |
|
{ |
|
s = signe(u); |
|
if (s < 0) err(lllger3); |
|
if (s) |
|
{ B[k+1]=(long)u; coeff(L,k,k)=un; fl[k]=1; } |
|
else |
|
{ B[k+1]=B[k]; fl[k]=0; } |
|
} |
|
} |
|
} else if (DEBUGLEVEL>5) {fprintferr(" %ld",k); flusherr();} |
|
REDI(x,h,L,(GEN)B[k],kmax,k,k-1); |
|
if (do_SWAPI(x,h,L,B,kmax,k,alpha,fl)) |
|
{ |
{ |
if (k>2) k--; |
if (!gram) ZRED(k,k-1, x,h,L,(GEN)B[k],kmax); |
|
else ZRED_gram(k,k-1, x,h,L,(GEN)B[k],kmax); |
} |
} |
|
if (do_ZSWAP(x,h,L,B,kmax,k,D,fl,gram)) |
|
{ |
|
if (MARKED == k) MARKED = k-1; |
|
else if (MARKED == k-1) MARKED = k; |
|
if (k > 2) k--; |
|
} |
else |
else |
{ |
{ |
for (l=k-2; l; l--) |
if (k != MARKED) |
{ |
for (l=k-2; l; l--) |
REDI(x,h,L,(GEN)B[l+1],kmax,k,l); |
|
if (low_stack(lim, stack_lim(av,1))) |
|
{ |
{ |
if(DEBUGMEM>1) err(warnmem,"lllgramall[1]"); |
if (!gram) ZRED(k,l, x,h,L,(GEN)B[l+1],kmax); |
gptr[0]=&B; gptr[1]=&L; gptr[2]=&h; gptr[3]=&x; |
else ZRED_gram(k,l, x,h,L,(GEN)B[l+1],kmax); |
gerepilemany(av,gptr,4); |
if (low_stack(lim, stack_lim(av,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"lllint[1], kmax = %ld", kmax); |
|
gerepileall(av,h?4:3,&B,&L,&x,&h); |
|
} |
} |
} |
} |
|
if (++k > n) break; |
if (++k > n) break; |
} |
} |
if (low_stack(lim, stack_lim(av,1))) |
if (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
if(DEBUGMEM>1) err(warnmem,"lllgramall[2]"); |
if(DEBUGMEM>1) err(warnmem,"lllint[2], kmax = %ld", kmax); |
gptr[0]=&B; gptr[1]=&L; gptr[2]=&h; gptr[3]=&x; |
gerepileall(av,h?4:3,&B,&L,&x,&h); |
gerepilemany(av,gptr,4); |
|
} |
} |
} |
} |
if (DEBUGLEVEL>3) fprintferr("\n"); |
if (DEBUGLEVEL>3) fprintferr("\n"); |
tetpil=avma; return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n)); |
if (ptB) *ptB = B; |
|
if (ptfl) *ptfl = fl; |
|
if (pth) *pth = h; |
|
if (MARKED && MARKED != 1) |
|
{ |
|
if (h) { swap(h[1], h[MARKED]); } else swap(x[1], x[MARKED]); |
|
} |
|
return h? h: x; |
} |
} |
|
|
static GEN |
/* Beware: this function can return NULL (dim x <= 1) */ |
lllgramall0(GEN x, long flag) { return lllgramall(x,100,flag); } |
GEN |
|
lllint_i(GEN x, long D, int gram, GEN *pth, GEN *ptfl, GEN *ptB) |
|
{ |
|
return lllint_marked(0, x,D,gram,pth,ptfl,ptB); |
|
} |
|
|
|
/* return x * lllint(x). No garbage collection */ |
|
GEN |
|
lllint_ip(GEN x, long D) |
|
{ |
|
GEN h = lllint_i(x, D, 0, NULL, NULL, NULL); |
|
if (!h) h = x; |
|
return h; |
|
} |
|
|
|
GEN |
|
lllall(GEN x, long D, int gram, long flag) |
|
{ |
|
gpmem_t av = avma; |
|
GEN fl, junk, h = lllint_i(x, D, gram, &junk, &fl, NULL); |
|
if (!h) return lll_trivial(x,flag); |
|
return gerepilecopy(av, lll_finish(h,fl,flag)); |
|
} |
|
|
|
GEN |
|
lllint(GEN x) { return lllall(x,LLLDFT,0, lll_IM); } |
|
|
|
GEN |
|
lllkerim(GEN x) { return lllall(x,LLLDFT,0, lll_ALL); } |
|
|
|
GEN |
|
lllgramint(GEN x) { return lllall(x,LLLDFT,1, lll_IM | lll_GRAM); } |
|
|
|
GEN |
|
lllgramkerim(GEN x) { return lllall(x,LLLDFT,1, lll_ALL | lll_GRAM); } |
|
|
static int |
static int |
pslg(GEN x) |
pslg(GEN x) |
{ |
{ |
long tx; |
long tx; |
if (gcmp0(x)) return 2; |
if (gcmp0(x)) return 2; |
tx=typ(x); return is_scalar_t(tx)? 3: lgef(x); |
tx = typ(x); return is_scalar_t(tx)? 3: lgef(x); |
} |
} |
|
|
static GEN |
static GEN |
lllgramallgen(GEN x, long flag) |
to_MP(GEN x, long prec) |
{ |
{ |
long av0=avma,av,tetpil,lx=lg(x),tu,i,j,k,l,n,lim; |
GEN y = cgetr(prec); gaffect(x, y); return y; |
GEN u,B,lam,q,cq,h,la,bb,p1,p2,p3,p4,fl; |
} |
int ps1,ps2,flc; |
static GEN |
|
col_to_MP(GEN x, long prec) |
|
{ |
|
long j, l = lg(x); |
|
GEN y = cgetg(l, t_COL); |
|
for (j=1; j<l; j++) y[j] = (long)to_MP((GEN)x[j], prec); |
|
return y; |
|
} |
|
|
if (typ(x) != t_MAT) err(typeer,"lllgramallgen"); |
static GEN |
n=lx-1; if (n<=1) return lllall_trivial(x,n,flag); |
to_mp(GEN x, long prec) |
if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramallgen"); |
{ |
|
GEN y; |
|
if (typ(x) == t_INT) return x; |
|
y = cgetr(prec); gaffect(x, y); return y; |
|
} |
|
static GEN |
|
col_to_mp(GEN x, long prec) |
|
{ |
|
long j, l = lg(x); |
|
GEN y = cgetg(l, t_COL); |
|
for (j=1; j<l; j++) y[j] = (long)to_mp((GEN)x[j], prec); |
|
return y; |
|
} |
|
static GEN |
|
mat_to_mp(GEN x, long prec) |
|
{ |
|
long j, l = lg(x); |
|
GEN y = cgetg(l, t_MAT); |
|
for (j=1; j<l; j++) y[j] = (long)col_to_mp((GEN)x[j], prec); |
|
return y; |
|
} |
|
|
fl = new_chunk(lx); |
static int |
|
REDgen(long k, long l, GEN h, GEN L, GEN B) |
|
{ |
|
GEN u, q, cq; |
|
long i; |
|
|
av=avma; lim=stack_lim(av,1); |
u = gcoeff(L,k,l); |
B=cgetg(lx+1,t_COL); |
if (pslg(u) < pslg((GEN)B[l+1])) return 0; |
B[1]=un; lam=cgetg(lx,t_MAT); |
|
for (j=1; j<lx; j++) lam[j]=lgetg(lx,t_COL); |
q = gdeuc(u,(GEN)B[l+1]); |
for (i=1; i<lx; i++) |
cq = ginv(content(q)); |
for (j=1; j<=i; j++) |
q = gneg(gmul(q,cq)); |
|
h[k] = ladd(gmul(cq,(GEN)h[k]), gmul(q,(GEN)h[l])); |
|
coeff(L,k,l) = ladd(gmul(cq,gcoeff(L,k,l)), gmul(q,(GEN)B[l+1])); |
|
for (i=1; i<l; i++) |
|
coeff(L,k,i) = ladd(gmul(cq,gcoeff(L,k,i)), gmul(q,gcoeff(L,l,i))); |
|
return 1; |
|
} |
|
|
|
static int |
|
do_SWAPgen(GEN h, GEN L, GEN B, long k, GEN fl, int *flc) |
|
{ |
|
GEN p1, la, la2, Bk; |
|
long ps1, ps2, i, j, lx; |
|
|
|
if (!fl[k-1]) return 0; |
|
|
|
la = gcoeff(L,k,k-1); la2 = gsqr(la); |
|
Bk = (GEN)B[k]; |
|
if (fl[k]) |
|
{ |
|
GEN q = gadd(la2, gmul((GEN)B[k-1],(GEN)B[k+1])); |
|
ps1 = pslg(gsqr(Bk)); |
|
ps2 = pslg(q); |
|
if (ps1 <= ps2 && (ps1 < ps2 || !*flc)) return 0; |
|
*flc = (ps1 != ps2); |
|
B[k] = ldiv(q, Bk); |
|
} |
|
|
|
swap(h[k-1], h[k]); lx = lg(L); |
|
for (j=1; j<=k-2; j++) swap(coeff(L,k-1,j), coeff(L,k,j)); |
|
if (fl[k]) |
|
{ |
|
for (i=k+1; i<lx; i++) |
{ |
{ |
if (j<i && !fl[j]) coeff(lam,i,j)=coeff(lam,j,i)=zero; |
GEN t = gcoeff(L,i,k); |
else |
p1 = gsub(gmul((GEN)B[k+1],gcoeff(L,i,k-1)), gmul(la,t)); |
{ |
coeff(L,i,k) = ldiv(p1, Bk); |
u=gcoeff(x,i,j); tu=typ(u); |
p1 = gadd(gmul(la,gcoeff(L,i,k-1)), gmul((GEN)B[k-1],t)); |
if (! is_scalar_t(tu) && tu != t_POL) err(lllger4); |
coeff(L,i,k-1) = ldiv(p1, Bk); |
for (k=1; k<j; k++) |
|
if (fl[k]) |
|
u=gdiv(gsub( gmul((GEN)B[k+1],u), |
|
gmul(gcoeff(lam,i,k),gcoeff(lam,j,k)) ), |
|
(GEN)B[k]); |
|
if (j<i) { coeff(lam,i,j)=(long)u; coeff(lam,j,i)=zero; } |
|
else |
|
{ |
|
if (!gcmp0(u)) { B[i+1]=(long)u; coeff(lam,i,i)=un; fl[i]=1; } |
|
else { B[i+1]=B[i]; coeff(lam,i,i)=zero; fl[i]=0; } |
|
} |
|
} |
|
} |
} |
k=2; h=idmat(n); flc=0; |
} |
for(;;) |
else if (!gcmp0(la)) |
{ |
{ |
u=gcoeff(lam,k,k-1); |
p1 = gdiv(la2, Bk); |
if (pslg(u) >= pslg((GEN)B[k])) |
B[k+1] = B[k] = (long)p1; |
|
for (i=k+2; i<=lx; i++) B[i] = ldiv(gmul(p1,(GEN)B[i]),Bk); |
|
for (i=k+1; i<lx; i++) |
|
coeff(L,i,k-1) = ldiv(gmul(la,gcoeff(L,i,k-1)), Bk); |
|
for (j=k+1; j<lx-1; j++) |
|
for (i=j+1; i<lx; i++) |
|
coeff(L,i,j) = ldiv(gmul(p1,gcoeff(L,i,j)), Bk); |
|
} |
|
else |
|
{ |
|
coeff(L,k,k-1) = zero; |
|
for (i=k+1; i<lx; i++) |
{ |
{ |
q=gdeuc(u,(GEN)B[k]); cq=gdivsg(1,content(q)); q=gmul(q,cq); flc=1; |
coeff(L,i,k) = coeff(L,i,k-1); |
h[k]=lsub(gmul(cq,(GEN)h[k]),gmul(q,(GEN)h[k-1])); |
coeff(L,i,k-1) = zero; |
coeff(lam,k,k-1)=lsub(gmul(cq,gcoeff(lam,k,k-1)),gmul(q,(GEN)B[k])); |
|
for (i=1; i<k-1; i++) |
|
coeff(lam,k,i)=lsub(gmul(cq,gcoeff(lam,k,i)),gmul(q,gcoeff(lam,k-1,i))); |
|
} |
} |
ps1 = pslg(gsqr((GEN)B[k])); |
B[k] = B[k-1]; fl[k] = 1; fl[k-1] = 0; |
p3 = gmul((GEN)B[k-1],(GEN)B[k+1]); |
} |
la=gcoeff(lam,k,k-1); p4 = gmul(la,gcoeff(lam,k,k-1)); |
return 1; |
ps2=pslg(gadd(p3,p4)); |
} |
if (fl[k-1] && (ps1>ps2 || (ps1==ps2 && flc) || !fl[k])) |
|
|
static void |
|
incrementalGSgen(GEN x, GEN L, GEN B, long k, GEN fl) |
|
{ |
|
GEN u = NULL; /* gcc -Wall */ |
|
long i, j, tu; |
|
for (j=1; j<=k; j++) |
|
if (j==k || fl[j]) |
{ |
{ |
flc = (ps1!=ps2); |
u = gcoeff(x,k,j); tu = typ(u); |
swap(h[k-1],h[k]); |
if (! is_scalar_t(tu) && tu != t_POL) err(lllger4); |
for (j=1; j<=k-2; j++) swap(coeff(lam,k-1,j), coeff(lam,k,j)); |
for (i=1; i<j; i++) |
if (fl[k]) |
if (fl[i]) |
{ |
{ |
for (i=k+1; i<=n; i++) |
u = gsub(gmul((GEN)B[i+1],u), gmul(gcoeff(L,k,i),gcoeff(L,j,i))); |
{ |
u = gdiv(u, (GEN)B[i]); |
bb=gcoeff(lam,i,k); |
} |
coeff(lam,i,k)=ldiv(gsub(gmul((GEN)B[k+1],gcoeff(lam,i,k-1)),gmul(la,bb)),(GEN)B[k]); |
coeff(L,k,j) = (long)u; |
coeff(lam,i,k-1)=ldiv(gadd(gmul(la,gcoeff(lam,i,k-1)),gmul((GEN)B[k-1],bb)),(GEN)B[k]); |
|
} |
|
B[k]=ldiv(gadd(p3,p4),(GEN)B[k]); |
|
} |
|
else |
|
{ |
|
if (!gcmp0(la)) |
|
{ |
|
p2=(GEN)B[k]; p1=gdiv(p4,p2); |
|
for (i=k+1; i<lx; i++) |
|
coeff(lam,i,k-1)=ldiv(gmul(la,gcoeff(lam,i,k-1)),p2); |
|
for (j=k+1; j<lx-1; j++) |
|
for (i=j+1; i<lx; i++) |
|
coeff(lam,i,j)=ldiv(gmul(p1,gcoeff(lam,i,j)),p2); |
|
B[k+1]=B[k]=(long)p1; |
|
for (i=k+2; i<=lx; i++) |
|
B[i]=ldiv(gmul(p1,(GEN)B[i]),p2); |
|
} |
|
else |
|
{ |
|
coeff(lam,k,k-1)=zero; |
|
for (i=k+1; i<lx; i++) |
|
{ coeff(lam,i,k)=coeff(lam,i,k-1); coeff(lam,i,k-1)=zero; } |
|
B[k]=B[k-1]; fl[k]=1; fl[k-1]=0; |
|
} |
|
} |
|
if (k>2) k--; |
|
} |
} |
|
if (gcmp0(u)) B[k+1] = B[k]; |
|
else |
|
{ |
|
B[k+1] = coeff(L,k,k); coeff(L,k,k) = un; fl[k] = 1; |
|
} |
|
} |
|
|
|
static GEN |
|
lllgramallgen(GEN x, long flag) |
|
{ |
|
long lx = lg(x), i, j, k, l, n; |
|
gpmem_t av0 = avma, av, lim; |
|
GEN B, L, h, fl; |
|
int flc; |
|
|
|
if (typ(x) != t_MAT) err(typeer,"lllgramallgen"); |
|
n = lx-1; if (n<=1) return lll_trivial(x,flag); |
|
if (lg(x[1]) != lx) err(mattype1,"lllgramallgen"); |
|
|
|
fl = cgetg(lx, t_VECSMALL); |
|
|
|
av = avma; lim = stack_lim(av,1); |
|
B = gscalcol(gun, lx); |
|
L = cgetg(lx,t_MAT); |
|
for (j=1; j<lx; j++) { L[j] = (long)zerocol(n); fl[j] = 0; } |
|
|
|
h = idmat(n); |
|
for (i=1; i<lx; i++) |
|
incrementalGSgen(x, L, B, i, fl); |
|
flc = 0; |
|
for(k=2;;) |
|
{ |
|
if (REDgen(k, k-1, h, L, B)) flc = 1; |
|
if (do_SWAPgen(h, L, B, k, fl, &flc)) { if (k > 2) k--; } |
else |
else |
{ |
{ |
for (l=k-2; l>=1; l--) |
for (l=k-2; l>=1; l--) |
{ |
if (REDgen(k, l, h, L, B)) flc = 1; |
u=gcoeff(lam,k,l); |
|
if (pslg(u)>=pslg((GEN)B[l+1])) |
|
{ |
|
q=gdeuc(u,(GEN)B[l+1]); cq=gdivsg(1,content(q)); |
|
q=gmul(q,cq); flc=1; |
|
h[k]=lsub(gmul(cq,(GEN)h[k]),gmul(q,(GEN)h[l])); |
|
coeff(lam,k,l)=lsub(gmul(cq,gcoeff(lam,k,l)),gmul(q,(GEN)B[l+1])); |
|
for (i=1; i<l; i++) |
|
coeff(lam,k,i)=lsub(gmul(cq,gcoeff(lam,k,i)),gmul(q,gcoeff(lam,l,i))); |
|
} |
|
} |
|
if (++k > n) break; |
if (++k > n) break; |
} |
} |
if (low_stack(lim, stack_lim(av,1))) |
if (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
GEN *gptr[4]; |
|
if(DEBUGMEM>1) err(warnmem,"lllgramallgen"); |
if(DEBUGMEM>1) err(warnmem,"lllgramallgen"); |
gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h; |
gerepileall(av,3,&B,&L,&h); |
gerepilemany(av,gptr,3); |
|
} |
} |
} |
} |
tetpil=avma; |
return gerepilecopy(av0, lll_finish(h,fl,flag)); |
return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n)); |
|
} |
} |
|
|
/* compute B[k], update mu(k,1..k-1) */ |
static GEN |
|
_mul2n(GEN x, long e) { return e? gmul2n(x, e): x; } |
|
|
|
/* E = scaling exponents |
|
* F = backward scaling exponents |
|
* X = lattice basis, Xs = associated scaled basis |
|
* Q = vector of Householder transformations |
|
* h = base change matrix (from X0) |
|
* R = from QR factorization of Xs[1..k-1] */ |
static int |
static int |
get_Gram_Schmidt(GEN x, GEN mu, GEN B, long k) |
HRS(int MARKED, long k, int prim, long kmax, GEN X, GEN Xs, GEN h, GEN R, |
|
GEN Q, GEN E, GEN F) |
{ |
{ |
GEN s,A = cgetg(k+1, t_COL); /* scratch space */ |
long e, i, N = lg(X[k]); |
long av,i,j; |
const long prec = MEDDEFAULTPREC; /* 128 bits */ |
A[1] = coeff(x,k,1); |
GEN q, tau2, rk; |
for(j=1;j<k;) |
int rounds = 0, overf; |
|
|
|
E[k] = prim? E[k-1]: 0; |
|
F[k] = 0; |
|
Xs[k] = E[k]? lmul2n((GEN)X[k], E[k]): (long)dummycopy((GEN)X[k]); |
|
rounds = 0; |
|
if (k == MARKED) goto DONE; /* no size reduction/scaling */ |
|
|
|
UP: |
|
rk = ApplyAllQ(Q, col_to_MP((GEN)Xs[k], prec), k); |
|
overf = 0; |
|
for (i = k-1; i > 0; i--) |
|
{ /* size seduction of Xs[k] */ |
|
GEN q2, mu, muint; |
|
long ex; |
|
gpmem_t av = avma; |
|
|
|
mu = mpdiv((GEN)rk[i], gcoeff(R,i,i)); |
|
if (!signe(mu)) { avma = av; continue; } |
|
mu = _mul2n(mu, -F[i]); /*backward scaling*/ |
|
ex = gexpo(mu); |
|
if (ex > 10) { |
|
overf = 1; /* large size reduction */ |
|
muint = ex > 30? ceil_safe(mu): ground(mu); |
|
} |
|
else if (fabs(gtodouble(mu)) > 0.51) |
|
muint = ground(mu); |
|
else { avma = av; continue; } |
|
|
|
q = gerepileuptoint(av, negi(muint)); |
|
q2= _mul2n(q, F[i]); /*backward scaling*/ |
|
Zupdate_col(k, i, q2, N-1, Xs); |
|
rk = gadd(rk, gmul(q2, (GEN)R[i])); |
|
/* if in HRS', propagate the last SR on X (for SWAP test) */ |
|
if (prim && (i == k-1)) { |
|
Zupdate_col(k, i, q, kmax, h); |
|
update_col(k, i, q, X); |
|
} |
|
} |
|
/* If large size-reduction performed, restart from exact data */ |
|
if (overf) |
{ |
{ |
coeff(mu,k,j) = ldiv((GEN)A[j],(GEN)B[j]); |
if (++rounds > 100) return 0; /* detect infinite loop */ |
j++; av = avma; |
goto UP; |
/* A[j] <-- x[k,j] - sum_{i<j} mu[j,i] A[i] */ |
|
s = gmul(gcoeff(mu,j,1),(GEN)A[1]); |
|
for (i=2; i<j; i++) s = gadd(s, gmul(gcoeff(mu,j,i),(GEN)A[i])); |
|
s = gneg(s); A[j] = lpileupto(av, gadd(gcoeff(x,k,j),s)); |
|
} |
} |
B[k] = A[k]; return (gsigne((GEN)B[k]) > 0); |
|
|
if (prim || k == 1) goto DONE; /* DONE */ |
|
|
|
rounds = 0; |
|
/* rescale Xs[k] ? */ |
|
tau2 = signe(rk[k])? gsqr((GEN)rk[k]): gzero; |
|
for (i = k+1; i < N; i++) |
|
if (signe(rk[i])) tau2 = mpadd(tau2, gsqr((GEN)rk[i])); |
|
q = mpdiv(gsqr(gcoeff(R,1,1)), tau2); |
|
e = gexpo(q)/2 + F[1]; /* ~ log_2 (||Xs[1]|| / tau), backward scaling*/ |
|
if (e > 0) |
|
{ /* rescale */ |
|
if (e > 30) e = 30; |
|
Xs[k] = lmul2n((GEN)Xs[k], e); |
|
E[k] += e; goto UP; /* one more round */ |
|
} |
|
else if (e < 0) |
|
{ /* enable backward scaling */ |
|
for (i = 1; i < k; i++) F[i] -= e; |
|
} |
|
|
|
DONE: |
|
rk = ApplyAllQ(Q, col_to_MP((GEN)Xs[k], prec), k); |
|
(void)FindApplyQ(rk, R, NULL, k, Q, prec); |
|
return 1; |
} |
} |
|
|
/* x = Gram(b_i). If precision problems return NULL if flag=1, error otherwise. |
static GEN |
* Quality ratio = Q = (alpha-1)/alpha. Suggested value: alpha = 100 |
rescale_to_int(GEN x) |
*/ |
{ |
|
long e, prec = gprecision(x); |
|
GEN y; |
|
if (!prec) return Q_primpart(x); |
|
|
|
y = gmul2n(x, bit_accuracy(prec) - gexpo(x)); |
|
return gcvtoi(y, &e); |
|
} |
|
|
GEN |
GEN |
lllgramintern(GEN x, long alpha, long flag, long prec) |
lll_scaled(int MARKED, GEN X0, long D) |
{ |
{ |
GEN xinit,L,h,B,L1,QR; |
GEN delta, X, Xs, h, R, Q, E, F; |
long retry = 2, av = avma,lim,l,i,j,k,k1,lx=lg(x),n,kmax,KMAX; |
long j, kmax = 1, k, N = lg(X0); |
long last_prec; |
gpmem_t lim, av, av0 = avma; |
|
int retry = 0; |
|
|
if (typ(x) != t_MAT) err(typeer,"lllgram"); |
delta = stor(D-1, DEFAULTPREC); |
n=lx-1; if (n && lg((GEN)x[1])!=lx) err(mattype1,"lllgram"); |
delta = divrs(delta,D); |
if (n<=1) return idmat(n); |
E = vecsmall_const(N-1, 0); |
xinit = x; lim = stack_lim(av,1); |
F = vecsmall_const(N-1, 0); |
for (k=2,j=1; j<lx; j++) |
|
|
av = avma; lim = stack_lim(av, 1); |
|
h = idmat(N-1); |
|
X = rescale_to_int(X0); |
|
|
|
PRECPB: |
|
k = 1; |
|
if (retry++) return _vec(h); |
|
Q = zerovec(N-1); |
|
Xs = zeromat(N-1, N-1); |
|
R = cgetg(N, t_MAT); |
|
for (j=1; j<N; j++) R[j] = (long)zerocol(N-1); |
|
|
|
while (k < N) |
{ |
{ |
GEN p1=(GEN)x[j]; |
gpmem_t av1; |
for (i=1; i<=j; i++) |
if (k == 1) { |
if (typ(p1[i]) == t_REAL) { l = lg((GEN)p1[i]); if (l>k) k=l; } |
HRS(MARKED, 1, 0, kmax, X, Xs, h, R, Q, E, F); |
|
k = 2; |
|
} |
|
if (k > kmax) { |
|
kmax = k; |
|
if (DEBUGLEVEL>3) {fprintferr(" K%ld",k);flusherr();} |
|
} |
|
if (!HRS(MARKED, k, 1, kmax, X, Xs, h, R, Q, E, F)) goto PRECPB; |
|
|
|
av1 = avma; |
|
if (mpcmp( mpmul(delta, gsqr(gcoeff(R, k-1,k-1))), |
|
mpadd(gsqr(gcoeff(R,k-1,k)), gsqr(gcoeff(R, k,k)))) > 0) |
|
{ |
|
if (DEBUGLEVEL>3 && k == kmax) |
|
{ |
|
GEN q = mpsub( mpmul(delta, gsqr(gcoeff(R, k-1,k-1))), |
|
gsqr(gcoeff(R,k-1,k))); |
|
fprintferr(" (%ld)", gexpo(q) - gexpo(gsqr(gcoeff(R, k,k)))); |
|
} |
|
swap(X[k], X[k-1]); |
|
swap(h[k], h[k-1]); |
|
if (MARKED == k) MARKED = k-1; |
|
else if (MARKED == k-1) MARKED = k; |
|
avma = av1; k--; |
|
} |
|
else |
|
{ |
|
avma = av1; |
|
if (!HRS(MARKED, k, 0, kmax, X, Xs, h, R, Q, E, F)) goto PRECPB; |
|
k++; |
|
} |
|
if (low_stack(lim, stack_lim(av,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"lllfp[1]"); |
|
gerepileall(av,5,&X,&Xs, &R,&h,&Q); |
|
} |
} |
} |
|
return gerepilecopy(av0, h); |
|
} |
|
|
|
/* x = Gram(b_i). If precision problems return NULL if flag=1, error otherwise. |
|
* Quality ratio = delta = (D-1)/D. Suggested value: D = 100 |
|
* |
|
* If MARKED != 0 make sure e[MARKED] is the first vector of the output basis |
|
* (which may then not be LLL-reduced) */ |
|
GEN |
|
lllfp_marked(int MARKED, GEN x, long D, long flag, long prec, int gram) |
|
{ |
|
GEN xinit,L,h,B,L1,delta, Q; |
|
long retry = 2, lx = lg(x), hx, l, i, j, k, k1, n, kmax, KMAX; |
|
gpmem_t av0 = avma, av, lim; |
|
|
|
if (typ(x) != t_MAT) err(typeer,"lllfp"); |
|
n = lx-1; if (n <= 1) return idmat(n); |
|
#if 0 /* doesn't work yet */ |
|
return lll_scaled(MARKED, x, D); |
|
#endif |
|
|
|
hx = lg(x[1]); |
|
if (gram && hx != lx) err(mattype1,"lllfp"); |
|
delta = stor(D-1, DEFAULTPREC); |
|
delta = divrs(delta,D); |
|
xinit = x; |
|
av = avma; lim = stack_lim(av,1); |
|
if (gram) { |
|
for (k=2,j=1; j<lx; j++) |
|
{ |
|
GEN p1=(GEN)x[j]; |
|
for (i=1; i<=j; i++) |
|
if (typ(p1[i]) == t_REAL) { l = lg(p1[i]); if (l>k) k=l; } |
|
} |
|
} else { |
|
for (k=2,j=1; j<lx; j++) |
|
{ |
|
GEN p1=(GEN)x[j]; |
|
for (i=1; i<hx; i++) |
|
if (typ(p1[i]) == t_REAL) { l = lg(p1[i]); if (l>k) k=l; } |
|
} |
|
} |
if (k == 2) |
if (k == 2) |
{ |
{ |
if (!prec) return lllgramint(x); |
if (!prec) return gram? lllgramint(x): lllint(x); |
x = gmul(x, realun(prec+1)); |
x = gmul(x, realun(prec+1)); |
} |
} |
else |
else |
{ |
{ |
if (prec < k) prec = k; |
if (prec < k) prec = k; |
x = gprec_w(x, prec+1); |
x = mat_to_mp(x, prec+1); |
} |
} |
/* kmax = maximum column index attained during this run |
/* kmax = maximum column index attained during this run |
* KMAX = same over all runs (after PRECPB) */ |
* KMAX = same over all runs (after PRECPB) */ |
kmax = KMAX = 1; |
kmax = KMAX = 1; |
|
|
|
Q = zerovec(n); |
|
|
PRECPB: |
PRECPB: |
switch(retry--) |
switch(retry--) |
{ |
{ |
case 2: h = idmat(n); break; /* entry */ |
case 2: h = idmat(n); break; /* entry */ |
case 1: |
case 1: |
if (kmax > 2) |
if (DEBUGLEVEL > 3) fprintferr("\n"); |
|
if (flag == 2) return _vec(h); |
|
if (gram && kmax > 2) |
{ /* some progress but precision loss, try again */ |
{ /* some progress but precision loss, try again */ |
prec = (prec<<1)-2; kmax = 1; |
prec = (prec<<1)-2; kmax = 1; |
if (DEBUGLEVEL > 3) fprintferr("\n"); |
if (DEBUGLEVEL) err(warnprec,"lllfp",prec); |
if (DEBUGLEVEL) err(warnprec,"lllgramintern",prec); |
x = gprec_w(xinit,prec); |
x = qf_base_change(gprec_w(xinit,prec),h,1); |
if (gram) x = qf_base_change(x,h,1); else x = gmul(x,h); |
{ |
gerepileall(av,2,&h,&x); |
GEN *gsav[2]; gsav[0]=&h; gsav[1]=&x; |
if (DEBUGLEVEL) err(warner,"lllfp starting over"); |
gerepilemany(av, gsav, 2); |
|
} |
|
if (DEBUGLEVEL) err(warner,"lllgramintern starting over"); |
|
break; |
break; |
} /* fall through */ |
} /* fall through */ |
case 0: /* give up */ |
case 0: /* give up */ |
if (DEBUGLEVEL > 3) fprintferr("\n"); |
if (DEBUGLEVEL > 3) fprintferr("\n"); |
if (DEBUGLEVEL) err(warner,"lllgramintern giving up"); |
if (DEBUGLEVEL) err(warner,"lllfp giving up"); |
if (flag) { avma=av; return NULL; } |
if (flag) { avma=av; return NULL; } |
if (DEBUGLEVEL) outerr(xinit); |
|
err(lllger3); |
err(lllger3); |
} |
} |
last_prec = prec; |
|
QR = cgetr(prec+1); affsr(alpha-1,QR); |
|
QR = divrs(QR,alpha); |
|
|
|
L=cgetg(lx,t_MAT); |
L=cgetg(lx,t_MAT); |
B=cgetg(lx,t_COL); |
B=cgetg(lx,t_COL); |
for (j=1; j<lx; j++) { L[j] = (long)zerocol(n); B[j] = zero; } |
for (j=1; j<lx; j++) { L[j] = (long)zerocol(n); B[j] = zero; } |
k=2; B[1]=coeff(x,1,1); |
if (gram && !incrementalGS(x, L, B, 1)) |
if (gcmp0((GEN)B[1])) |
|
{ |
{ |
if (flag) return NULL; |
if (!flag) return NULL; |
err(lllger3); |
err(lllger3); |
} |
} |
if (DEBUGLEVEL>5) fprintferr("k ="); |
if (DEBUGLEVEL>5) fprintferr("k ="); |
for(;;) |
for(k=2;;) |
{ |
{ |
if (k>kmax) |
if (k > kmax) |
{ |
{ |
kmax = k; if (KMAX < kmax) KMAX = kmax; |
kmax = k; if (KMAX < kmax) KMAX = kmax; |
if (DEBUGLEVEL>3) {fprintferr(" K%ld",k);flusherr();} |
if (DEBUGLEVEL>3) {fprintferr(" K%ld",k);flusherr();} |
if (!get_Gram_Schmidt(x,L,B,k)) goto PRECPB; |
if (gram) j = incrementalGS(x, L, B, k); |
|
else j = Householder_get_mu(x, L, B, k, Q, prec); |
|
if (!j) goto PRECPB; |
} |
} |
else if (DEBUGLEVEL>5) fprintferr(" %ld",k); |
else if (DEBUGLEVEL>5) fprintferr(" %ld",k); |
L1 = gcoeff(L,k,k-1); |
L1 = gcoeff(L,k,k-1); |
if (typ(L1) == t_REAL && |
if (typ(L1) == t_REAL && expo(L1) + 20 > bit_accuracy(lg(L1))) |
(2*gexpo(L1) > bit_accuracy(lg(L1)) || 2*lg(L1) < last_prec)) |
|
{ |
{ |
last_prec = lg(L1); |
|
if (DEBUGLEVEL>3) |
if (DEBUGLEVEL>3) |
fprintferr("\nRecomputing Gram-Schmidt, kmax = %ld, prec was %ld\n", |
fprintferr("\nRecomputing Gram-Schmidt, kmax = %ld\n", kmax); |
kmax,last_prec); |
if (!gram) goto PRECPB; |
for (k1=1; k1<=kmax; k1++) |
for (k1=1; k1<=kmax; k1++) |
if (!get_Gram_Schmidt(x,L,B,k1)) goto PRECPB; |
if (!incrementalGS(x, L, B, k1)) goto PRECPB; |
} |
} |
RED(x,h,L,KMAX,k,k-1); |
if (k != MARKED) |
if (do_SWAP(x,h,L,B,kmax,k,QR)) |
|
{ |
{ |
|
if (!gram) j = RED(k,k-1, x,h,L,KMAX); |
|
else j = RED_gram(k,k-1, x,h,L,KMAX); |
|
if (!j) goto PRECPB; |
|
} |
|
if (do_SWAP(x,h,L,B,kmax,k,delta,gram)) |
|
{ |
|
if (MARKED == k) MARKED = k-1; |
|
else if (MARKED == k-1) MARKED = k; |
if (!B[k]) goto PRECPB; |
if (!B[k]) goto PRECPB; |
|
Q[k] = Q[k-1] = zero; |
if (k>2) k--; |
if (k>2) k--; |
} |
} |
else |
else |
{ |
{ |
for (l=k-2; l; l--) RED(x,h,L,KMAX,k,l); |
if (k != MARKED) |
if (++k > n) break; |
for (l=k-2; l; l--) |
|
{ |
|
if (!gram) j = RED(k,l, x,h,L,KMAX); |
|
else j = RED_gram(k,l, x,h,L,KMAX); |
|
if (!j) goto PRECPB; |
|
if (low_stack(lim, stack_lim(av,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"lllfp[1], kmax = %ld", kmax); |
|
gerepileall(av,5,&B,&L,&h,&x,&Q); |
|
} |
|
} |
|
if (++k > n) |
|
{ |
|
if (!gram && Q[n-1] == zero) |
|
{ |
|
if (DEBUGLEVEL>3) fprintferr("\nChecking LLL basis\n"); |
|
j = Householder_get_mu(gmul(xinit,h), L, B, n, Q, prec); |
|
if (!j) goto PRECPB; |
|
k = 2; continue; |
|
} |
|
break; |
|
} |
} |
} |
if (low_stack(lim, stack_lim(av,1))) |
if (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
GEN *gptr[5]; gptr[0]=&B; gptr[1]=&L; gptr[2]=&h; gptr[3]=&x; gptr[4]=&QR; |
if(DEBUGMEM>1) err(warnmem,"lllfp[2], kmax = %ld", kmax); |
if(DEBUGMEM>1) |
gerepileall(av,5,&B,&L,&h,&x,&Q); |
{ |
|
if (DEBUGLEVEL > 3) fprintferr("\n"); |
|
err(warnmem,"lllgram"); |
|
} |
|
gerepilemany(av,gptr,5); |
|
} |
} |
} |
} |
if (DEBUGLEVEL>3) fprintferr("\n"); |
if (DEBUGLEVEL>3) fprintferr("\n"); |
return gerepilecopy(av, h); |
if (MARKED && MARKED != 1) swap(h[1], h[MARKED]); |
|
return gerepilecopy(av0, h); |
} |
} |
|
|
static GEN |
GEN |
lllgram_noerr(GEN x,long prec) { return lllgramintern(x,100,1,prec); } |
lllgramintern(GEN x, long D, long flag, long prec) |
|
{ |
|
return lllfp_marked(0, x,D,flag,prec,1); |
|
} |
|
|
GEN |
GEN |
lllgram(GEN x,long prec) { return lllgramintern(x,100,0,prec); } |
lllintern(GEN x, long D, long flag, long prec) |
|
{ |
|
return lllfp_marked(0, x,D,flag,prec,0); |
|
} |
|
|
GEN |
GEN |
|
lllgram(GEN x,long prec) { return lllgramintern(x,LLLDFT,0,prec); } |
|
|
|
GEN |
|
lll(GEN x,long prec) { return lllintern(x,LLLDFT,0,prec); } |
|
|
|
GEN |
qflll0(GEN x, long flag, long prec) |
qflll0(GEN x, long flag, long prec) |
{ |
{ |
switch(flag) |
switch(flag) |
Line 960 qflll0(GEN x, long flag, long prec) |
|
Line 1308 qflll0(GEN x, long flag, long prec) |
|
case 0: return lll(x,prec); |
case 0: return lll(x,prec); |
case 1: return lllint(x); |
case 1: return lllint(x); |
case 2: return lllintpartial(x); |
case 2: return lllintpartial(x); |
case 3: return lllrat(x); |
|
case 4: return lllkerim(x); |
case 4: return lllkerim(x); |
case 5: return lllkerimgen(x); |
case 5: return lllkerimgen(x); |
case 7: return lll1(x,prec); |
|
case 8: return lllgen(x); |
case 8: return lllgen(x); |
case 9: return lllintwithcontent(x); |
|
default: err(flagerr,"qflll"); |
default: err(flagerr,"qflll"); |
} |
} |
return NULL; /* not reached */ |
return NULL; /* not reached */ |
Line 980 qflllgram0(GEN x, long flag, long prec) |
|
Line 1325 qflllgram0(GEN x, long flag, long prec) |
|
case 1: return lllgramint(x); |
case 1: return lllgramint(x); |
case 4: return lllgramkerim(x); |
case 4: return lllgramkerim(x); |
case 5: return lllgramkerimgen(x); |
case 5: return lllgramkerimgen(x); |
case 7: return lllgram1(x,prec); |
|
case 8: return lllgramgen(x); |
case 8: return lllgramgen(x); |
default: err(flagerr,"qflllgram"); |
default: err(flagerr,"qflllgram"); |
} |
} |
return NULL; /* not reached */ |
return NULL; /* not reached */ |
} |
} |
|
|
/* x est la matrice d'une base b_i; retourne la matrice u (entiere |
GEN |
* unimodulaire) d'une base LLL-reduite c_i en fonction des b_i (la base |
gram_matrix(GEN b) |
* reduite est c=b*u). |
|
*/ |
|
static GEN |
|
lll_proto(GEN x, GEN f(GEN, long), long prec) |
|
{ |
{ |
long lx=lg(x),i,j,av,av1; |
long i,j, lx = lg(b); |
GEN g; |
GEN g; |
|
if (typ(b) != t_MAT) err(typeer,"gram"); |
if (typ(x) != t_MAT) err(typeer,"lll_proto"); |
g = cgetg(lx,t_MAT); |
if (lx==1) return cgetg(1,t_MAT); |
|
av=avma; g=cgetg(lx,t_MAT); |
|
for (j=1; j<lx; j++) g[j]=lgetg(lx,t_COL); |
|
for (i=1; i<lx; i++) |
for (i=1; i<lx; i++) |
|
{ |
|
g[i] = lgetg(lx,t_COL); |
for (j=1; j<=i; j++) |
for (j=1; j<=i; j++) |
coeff(g,i,j)=coeff(g,j,i)=(long)gscal((GEN)x[i],(GEN)x[j]); |
coeff(g,i,j) = coeff(g,j,i) = (long)gscal((GEN)b[i],(GEN)b[j]); |
av1=avma; x = f(g,prec); |
} |
if (!x) { avma=av; return NULL; } |
return g; |
return gerepile(av,av1,x); |
|
} |
} |
|
|
GEN |
GEN |
lllintern(GEN x,long flag,long prec) |
lllgen(GEN x) { |
{ |
gpmem_t av = avma; |
return lll_proto(x,flag? lllgram_noerr: lllgram,prec); |
return gerepileupto(av, lllgramgen(gram_matrix(x))); |
} |
} |
|
|
GEN |
GEN |
lll(GEN x,long prec) { return lll_proto(x,lllgram,prec); } |
lllkerimgen(GEN x) { |
|
gpmem_t av = avma; |
|
return gerepileupto(av, lllgramallgen(gram_matrix(x),lll_ALL)); |
|
} |
|
|
GEN |
GEN |
lll1(GEN x,long prec) { return lll_proto(x,lllgram1,prec); } |
|
|
|
GEN |
|
lllrat(GEN x) { return lll_proto(x,lllgram,lll_ALL); } |
|
|
|
GEN |
|
lllint(GEN x) { return lll_proto(x,lllgramall0,lll_IM); } |
|
|
|
GEN |
|
lllgen(GEN x) { return lll_proto(x,lllgramallgen,lll_IM); } |
|
|
|
GEN |
|
lllgramgen(GEN x) { return lllgramallgen(x,lll_IM); } |
lllgramgen(GEN x) { return lllgramallgen(x,lll_IM); } |
|
|
GEN |
GEN |
lllgramkerimgen(GEN x) { return lllgramallgen(x,lll_ALL); } |
lllgramkerimgen(GEN x) { return lllgramallgen(x,lll_ALL); } |
|
|
static GEN |
|
lllkerim_proto(GEN x, GEN f(GEN,long)) |
|
{ |
|
long lx=lg(x), i,j,av,av1; |
|
GEN g; |
|
|
|
if (typ(x) != t_MAT) err(typeer,"lllkerim_proto"); |
|
if (lx==1) |
|
{ |
|
g=cgetg(3,t_VEC); |
|
g[1]=lgetg(1,t_MAT); |
|
g[2]=lgetg(1,t_MAT); return g; |
|
} |
|
if (lg((GEN)x[1])==1) |
|
{ |
|
g=cgetg(3,t_VEC); |
|
g[1]=(long)idmat(lx-1); |
|
g[2]=lgetg(1,t_MAT); return g; |
|
} |
|
av=avma; g=cgetg(lx,t_MAT); |
|
for (j=1; j<lx; j++) g[j]=lgetg(lx,t_COL); |
|
for (i=1; i<lx; i++) |
|
for (j=1; j<=i; j++) |
|
coeff(g,i,j)=coeff(g,j,i)=(long)gscal((GEN)x[i],(GEN)x[j]); |
|
av1=avma; return gerepile(av,av1,f(g,lll_ALL)); |
|
} |
|
|
|
GEN |
|
lllkerim(GEN x) { return lllkerim_proto(x,lllgramall0); } |
|
|
|
GEN |
|
lllkerimgen(GEN x) { return lllkerim_proto(x,lllgramallgen); } |
|
|
|
/* x est ici la matrice de GRAM des bi. */ |
|
GEN |
|
lllgram1(GEN x, long prec) |
|
{ |
|
GEN mu,u,B,BB,BK,p,q,r,cst,unreel,sv,mu1,mu2; |
|
long av,tetpil,lim,l,i,j,k,lx=lg(x),n,e; |
|
|
|
if (typ(x) != t_MAT) err(typeer,"lllgram1"); |
|
if (lg((GEN)x[1])!=lx) err(mattype1,"lllgram1"); n=lx-1; |
|
if (n<=1) return idmat(n); |
|
cst=gdivgs(stoi(99),100); /* LLL proposent 0.75 */ |
|
if (prec) |
|
{ |
|
unreel = realun(prec+1); |
|
x = gmul(x,unreel); |
|
cst = gmul(cst,unreel); |
|
} |
|
av=avma; lim=stack_lim(av,1); |
|
mu=gtrans(sqred(x)); B=cgetg(lx,t_COL); |
|
for (i=1,l=0; i<=n; i++) |
|
{ |
|
if (gsigne((GEN)(B[i]=coeff(mu,i,i)))>0) l++; |
|
coeff(mu,i,i)=un; |
|
} |
|
if (l<n) err(lllger3); |
|
|
|
u=idmat(n); k=2; |
|
do |
|
{ |
|
if (!gcmp0(r=grndtoi(gcoeff(mu,k,k-1),&e))) |
|
{ |
|
u[k]=lsub((GEN)u[k],gmul(r,(GEN)u[k-1])); |
|
for (j=1; j<k-1; j++) |
|
coeff(mu,k,j)=lsub(gcoeff(mu,k,j),gmul(r,gcoeff(mu,k-1,j))); |
|
mu1=(GEN)(coeff(mu,k,k-1)=lsub(gcoeff(mu,k,k-1),r)); |
|
} |
|
else mu1=gcoeff(mu,k,k-1); |
|
q=gmul((GEN)B[k-1],gsub(cst,mu2=gsqr(mu1))); |
|
if (gcmp(q,(GEN)B[k])>0) |
|
{ |
|
BB=gadd((GEN)B[k],gmul((GEN)B[k-1],mu2)); |
|
coeff(mu,k,k-1)=ldiv(gmul(mu1,(GEN)B[k-1]),BB); |
|
B[k]=lmul((GEN)B[k-1],BK=gdiv((GEN)B[k],BB)); |
|
B[k-1]=(long)BB; |
|
swap(u[k-1],u[k]); |
|
for (j=1; j<=k-2; j++) swap(coeff(mu,k-1,j), coeff(mu,k,j)); |
|
for (i=k+1; i<=n; i++) |
|
{ |
|
p=gcoeff(mu,i,k); |
|
coeff(mu,i,k)=lsub(gcoeff(mu,i,k-1),gmul(mu1,p)); |
|
coeff(mu,i,k-1)=ladd(gmul(BK,p),gmul(gcoeff(mu,k,k-1),gcoeff(mu,i,k-1))); |
|
} |
|
if (k>2) k--; |
|
} |
|
else |
|
{ |
|
for (l=k-2; l; l--) |
|
{ |
|
if (!gcmp0(r=grndtoi(gcoeff(mu,k,l),&e))) |
|
{ |
|
u[k]=lsub((GEN)u[k],gmul(r,(GEN)u[l])); |
|
for (j=1; j<l; j++) |
|
coeff(mu,k,j)=lsub(gcoeff(mu,k,j),gmul(r,gcoeff(mu,l,j))); |
|
coeff(mu,k,l)=lsub(gcoeff(mu,k,l),r); |
|
} |
|
} |
|
k++; |
|
} |
|
if (low_stack(lim, stack_lim(av,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"lllgram1"); |
|
tetpil=avma; |
|
sv=cgetg(4,t_VEC); |
|
sv[1]=lcopy(B); sv[2]=lcopy(u); sv[3]=lcopy(mu); |
|
sv=gerepile(av,tetpil,sv); |
|
B=(GEN)sv[1]; u=(GEN)sv[2]; mu=(GEN)sv[3]; |
|
} |
|
} |
|
while (k<=n); |
|
return gerepilecopy(av,u); |
|
} |
|
|
|
GEN |
|
lllgramint(GEN x) |
|
{ |
|
return lllgramall0(x,lll_IM); |
|
} |
|
|
|
GEN |
|
lllgramkerim(GEN x) |
|
{ |
|
return lllgramall0(x,lll_ALL); |
|
} |
|
|
|
/* This routine is functionally similar to lllint when all = 0. |
/* This routine is functionally similar to lllint when all = 0. |
* |
* |
* The input is an integer matrix mat (not necessarily square) whose |
* The input is an integer matrix mat (not necessarily square) whose |
|
|
lllintpartialall(GEN m, long flag) |
lllintpartialall(GEN m, long flag) |
{ |
{ |
const long ncol = lg(m)-1; |
const long ncol = lg(m)-1; |
const long ltop1 = avma; |
const gpmem_t ltop1 = avma; |
long ncolnz; |
long ncolnz; |
GEN tm1, tm2, mid, *gptr[4]; |
GEN tm1, tm2, mid; |
|
|
if (typ(m) != t_MAT) err(typeer,"lllintpartialall"); |
if (typ(m) != t_MAT) err(typeer,"lllintpartialall"); |
if (ncol <= 1) return idmat(ncol); |
if (ncol <= 1) return idmat(ncol); |
Line 1217 lllintpartialall(GEN m, long flag) |
|
Line 1419 lllintpartialall(GEN m, long flag) |
|
*/ |
*/ |
while (npass2 < 2 || progress) |
while (npass2 < 2 || progress) |
{ |
{ |
GEN dot12new,q = gdivround(dot12, dot22); |
GEN dot12new,q = diviiround(dot12, dot22); |
|
|
npass2++; progress = signe(q); |
npass2++; progress = signe(q); |
if (progress) |
if (progress) |
Line 1244 lllintpartialall(GEN m, long flag) |
|
Line 1446 lllintpartialall(GEN m, long flag) |
|
|
|
/* Occasionally (including final pass) do garbage collection. */ |
/* Occasionally (including final pass) do garbage collection. */ |
if (npass2 % 8 == 0 || !progress) |
if (npass2 % 8 == 0 || !progress) |
{ |
gerepileall(ltop1, 4, &dot11,&dot12,&dot22,&tm); |
gptr[0] = &dot11; gptr[1] = &dot12; |
|
gptr[2] = &dot22; gptr[3] = &tm; |
|
gerepilemany(ltop1, gptr, 4); |
|
} |
|
} /* while npass2 < 2 || progress */ |
} /* while npass2 < 2 || progress */ |
|
|
{ |
{ |
Line 1282 lllintpartialall(GEN m, long flag) |
|
Line 1480 lllintpartialall(GEN m, long flag) |
|
GEN q1neg = subii(mulii(dot12, dot2i), mulii(dot22, dot1i)); |
GEN q1neg = subii(mulii(dot12, dot2i), mulii(dot22, dot1i)); |
GEN q2neg = subii(mulii(dot12, dot1i), mulii(dot11, dot2i)); |
GEN q2neg = subii(mulii(dot12, dot1i), mulii(dot11, dot2i)); |
|
|
q1neg = gdivround(q1neg, det12); |
q1neg = diviiround(q1neg, det12); |
q2neg = gdivround(q2neg, det12); |
q2neg = diviiround(q2neg, det12); |
coeff(tm1, 1, icol) = ladd(gmul(q1neg, gcoeff(tm,1,1)), |
coeff(tm1, 1, icol) = ladd(gmul(q1neg, gcoeff(tm,1,1)), |
gmul(q2neg, gcoeff(tm,1,2))); |
gmul(q2neg, gcoeff(tm,1,2))); |
coeff(tm1, 2, icol) = ladd(gmul(q1neg, gcoeff(tm,2,1)), |
coeff(tm1, 2, icol) = ladd(gmul(q1neg, gcoeff(tm,2,1)), |
Line 1298 lllintpartialall(GEN m, long flag) |
|
Line 1496 lllintpartialall(GEN m, long flag) |
|
fprintferr("tm1 = %Z",tm1); |
fprintferr("tm1 = %Z",tm1); |
fprintferr("mid = %Z",mid); |
fprintferr("mid = %Z",mid); |
} |
} |
gptr[0] = &tm1; gptr[1] = ∣ |
gerepileall(ltop1,2, &tm1, &mid); |
gerepilemany(ltop1, gptr, 2); |
|
{ |
{ |
/* For each pair of column vectors v and w in mid * tm2, |
/* For each pair of column vectors v and w in mid * tm2, |
* try to replace (v, w) by (v, v - q*w) for some q. |
* try to replace (v, w) by (v, v - q*w) for some q. |
* We compute all inner products and check them repeatedly. |
* We compute all inner products and check them repeatedly. |
*/ |
*/ |
const long ltop3 = avma; /* Excludes region with tm1 and mid */ |
const gpmem_t ltop3 = avma; /* Excludes region with tm1 and mid */ |
long icol, lim, reductions, npass = 0; |
const gpmem_t lim = stack_lim(ltop3,1); |
|
long icol, reductions, npass = 0; |
GEN dotprd = cgetg(ncol+1, t_MAT); |
GEN dotprd = cgetg(ncol+1, t_MAT); |
|
|
tm2 = idmat(ncol); |
tm2 = idmat(ncol); |
Line 1319 lllintpartialall(GEN m, long flag) |
|
Line 1517 lllintpartialall(GEN m, long flag) |
|
coeff(dotprd,jcol,icol) = coeff(dotprd,icol,jcol) = |
coeff(dotprd,jcol,icol) = coeff(dotprd,icol,jcol) = |
(long)gscali((GEN)mid[icol], (GEN)mid[jcol]); |
(long)gscali((GEN)mid[icol], (GEN)mid[jcol]); |
} /* for icol */ |
} /* for icol */ |
lim = stack_lim(ltop3,1); |
|
for(;;) |
for(;;) |
{ |
{ |
reductions = 0; |
reductions = 0; |
Line 1330 lllintpartialall(GEN m, long flag) |
|
Line 1527 lllintpartialall(GEN m, long flag) |
|
|
|
for (ijdif=1; ijdif < ncol; ijdif++) |
for (ijdif=1; ijdif < ncol; ijdif++) |
{ |
{ |
const long previous_avma = avma; |
const gpmem_t previous_avma = avma; |
|
|
jcol = (icol + ijdif - 1) % ncol; jcol++; /* Hack for NeXTgcc 2.5.8 */ |
jcol = (icol + ijdif - 1) % ncol; jcol++; /* Hack for NeXTgcc 2.5.8 */ |
k1 = (cmpii(gcoeff(dotprd,icol,icol), |
k1 = (cmpii(gcoeff(dotprd,icol,icol), |
Line 1338 lllintpartialall(GEN m, long flag) |
|
Line 1535 lllintpartialall(GEN m, long flag) |
|
? icol : jcol; /* index of larger column */ |
? icol : jcol; /* index of larger column */ |
k2 = icol + jcol - k1; /* index of smaller column */ |
k2 = icol + jcol - k1; /* index of smaller column */ |
codi = gcoeff(dotprd,k2,k2); |
codi = gcoeff(dotprd,k2,k2); |
q = gcmp0(codi)? gzero: gdivround(gcoeff(dotprd,k1,k2), codi); |
q = gcmp0(codi)? gzero: diviiround(gcoeff(dotprd,k1,k2), codi); |
|
|
/* Try to subtract a multiple of column k2 from column k1. */ |
/* Try to subtract a multiple of column k2 from column k1. */ |
if (gcmp0(q)) avma = previous_avma; |
if (gcmp0(q)) avma = previous_avma; |
Line 1360 lllintpartialall(GEN m, long flag) |
|
Line 1557 lllintpartialall(GEN m, long flag) |
|
if (low_stack(lim, stack_lim(ltop3,1))) |
if (low_stack(lim, stack_lim(ltop3,1))) |
{ |
{ |
if(DEBUGMEM>1) err(warnmem,"lllintpartialall"); |
if(DEBUGMEM>1) err(warnmem,"lllintpartialall"); |
gptr[0] = &dotprd; gptr[1] = &tm2; |
gerepileall(ltop3, 2, &dotprd,&tm2); |
gerepilemany(ltop3, gptr, 2); |
|
} |
} |
} /* for icol */ |
} /* for icol */ |
if (!reductions) break; |
if (!reductions) break; |
Line 1423 lllintpartial(GEN mat) |
|
Line 1619 lllintpartial(GEN mat) |
|
|
|
/********************************************************************/ |
/********************************************************************/ |
/** **/ |
/** **/ |
/** LINEAR & ALGEBRAIC DEPENDANCE **/ |
/** LINEAR & ALGEBRAIC DEPENDENCE **/ |
/** **/ |
/** **/ |
/********************************************************************/ |
/********************************************************************/ |
|
|
|
GEN pslq(GEN x, long prec); |
|
GEN pslqL2(GEN x, long prec); |
|
|
GEN |
GEN |
lindep0(GEN x,long bit,long prec) |
lindep0(GEN x,long bit,long prec) |
{ |
{ |
if (!bit) return lindep(x,prec); |
switch (bit) |
if (bit>0) return lindep2(x,bit); |
{ |
return deplin(x); |
case 0: return pslq(x,prec); |
|
case -1: return lindep(x,prec); |
|
case -2: return deplin(x); |
|
case -3: return pslqL2(x,prec); |
|
default: return lindep2(x,labs(bit)); |
|
} |
} |
} |
|
|
static int |
static int |
Line 1446 real_indep(GEN re, GEN im, long bitprec) |
|
Line 1650 real_indep(GEN re, GEN im, long bitprec) |
|
GEN |
GEN |
lindep2(GEN x, long bit) |
lindep2(GEN x, long bit) |
{ |
{ |
long tx=typ(x),lx=lg(x),ly,i,j,e, av = avma; |
long tx=typ(x), lx=lg(x), ly, i, j, e; |
|
gpmem_t av = avma; |
GEN re,im,p1,p2; |
GEN re,im,p1,p2; |
|
|
if (! is_vec_t(tx)) err(typeer,"lindep2"); |
if (! is_vec_t(tx)) err(typeer,"lindep2"); |
Line 1467 lindep2(GEN x, long bit) |
|
Line 1672 lindep2(GEN x, long bit) |
|
p1[lx] = lcvtoi(gshift((GEN)re[i],bit),&e); |
p1[lx] = lcvtoi(gshift((GEN)re[i],bit),&e); |
if (im) p1[lx+1] = lcvtoi(gshift((GEN)im[i],bit),&e); |
if (im) p1[lx+1] = lcvtoi(gshift((GEN)im[i],bit),&e); |
} |
} |
p1 = (GEN)gmul(p2,lllint(p2))[1]; |
p1 = (GEN)lllint_ip(p2,100)[1]; |
p1[0] = evaltyp(t_VEC) | evallg(lx); |
p1[0] = evaltyp(t_VEC) | evallg(lx); |
return gerepilecopy(av, p1); |
return gerepilecopy(av, p1); |
} |
} |
Line 1478 lindep(GEN x, long prec) |
|
Line 1683 lindep(GEN x, long prec) |
|
{ |
{ |
GEN *b,*be,*bn,**m,qzer; |
GEN *b,*be,*bn,**m,qzer; |
GEN c1,c2,c3,px,py,pxy,re,im,p1,p2,r,f,em; |
GEN c1,c2,c3,px,py,pxy,re,im,p1,p2,r,f,em; |
long i,j,fl,i1, lx = lg(x), tx = typ(x), n = lx-1; |
long i,j,fl,k, lx = lg(x), tx = typ(x), n = lx-1; |
long av = avma, lim = stack_lim(av,1), av0,av1,tetpil; |
gpmem_t av = avma, lim = stack_lim(av,1), av0, av1; |
const long EXP = - bit_accuracy(prec) + 2*n; |
const long EXP = - bit_accuracy(prec) + 2*n; |
|
|
if (! is_vec_t(tx)) err(typeer,"lindep"); |
if (! is_vec_t(tx)) err(typeer,"lindep"); |
if (lx<=2) return cgetg(1,t_VEC); |
if (n <= 1) return cgetg(1,t_VEC); |
x = gmul(x, realun(prec)); |
x = gmul(x, realun(prec)); |
re=greal(x); im=gimag(x); |
re = greal(x); |
|
im = gimag(x); |
/* independant over R ? */ |
/* independant over R ? */ |
if (lx == 3 && real_indep(re,im,bit_accuracy(prec))) |
if (n == 2 && real_indep(re,im,bit_accuracy(prec))) |
{ avma = av; return cgetg(1, t_VEC); } |
{ avma = av; return cgetg(1, t_VEC); } |
|
if (EXP > -10) err(precer,"lindep"); |
|
|
qzer = new_chunk(lx); |
qzer = cgetg(lx, t_VECSMALL); |
b = (GEN*)idmat(n); |
b = (GEN*)idmat(n); |
be= (GEN*)new_chunk(lx); |
be= (GEN*)new_chunk(lx); |
bn= (GEN*)new_chunk(lx); |
bn= (GEN*)new_chunk(lx); |
m = (GEN**)new_chunk(lx); |
m = (GEN**)new_chunk(lx); |
for (i=1; i<=n; i++) |
for (i=1; i<=n; i++) |
{ |
{ |
bn[i]=cgetr(prec+1); |
bn[i] = cgetr(prec+1); |
be[i]=cgetg(lx,t_COL); |
be[i] = cgetg(lx,t_COL); |
m[i] = (GEN*)new_chunk(lx); |
m[i] = (GEN*)new_chunk(lx); |
for (j=1; j<i ; j++) m[i][j]=cgetr(prec+1); |
for (j=1; j<i ; j++) m[i][j] = cgetr(prec+1); |
for (j=1; j<=n; j++) be[i][j]=lgetr(prec+1); |
for (j=1; j<=n; j++) be[i][j]= lgetr(prec+1); |
} |
} |
px=sqscal(re); |
px = sqscal(re); |
py=sqscal(im); pxy=gscal(re,im); |
py = sqscal(im); pxy = gscal(re,im); |
p1=mpsub(mpmul(px,py),gsqr(pxy)); |
p1 = mpsub(mpmul(px,py), gsqr(pxy)); |
if (quazero(re)) { re=im; px=py; fl=1; } else fl=quazero(p1); |
if (quazero(px)) { re = im; px = py; fl = 1; } else fl = quazero(p1); |
av0 = av1 = avma; |
av0 = av1 = avma; |
for (i=1; i<=n; i++) |
for (i=1; i<=n; i++) |
{ |
{ |
p2 = gscal(b[i],re); |
p2 = gscal(b[i],re); |
if (fl) p2=gmul(gdiv(p2,px),re); |
if (fl) p2 = gmul(gdiv(p2,px),re); |
else |
else |
{ |
{ |
GEN p5,p6,p7; |
GEN p5,p6,p7; |
p5=gscal(b[i],im); |
p5 = gscal(b[i],im); |
p6=gdiv(gsub(gmul(py,p2),gmul(pxy,p5)),p1); |
p6 = gdiv(gsub(gmul(py,p2),gmul(pxy,p5)), p1); |
p7=gdiv(gsub(gmul(px,p5),gmul(pxy,p2)),p1); |
p7 = gdiv(gsub(gmul(px,p5),gmul(pxy,p2)), p1); |
p2=gadd(gmul(p6,re),gmul(p7,im)); |
p2 = gadd(gmul(p6,re), gmul(p7,im)); |
} |
} |
if (tx!=t_COL) p2=gtrans(p2); |
if (tx != t_COL) p2 = gtrans(p2); |
p2=gsub(b[i],p2); |
p2 = gsub(b[i],p2); |
for (j=1; j<i; j++) |
for (j=1; j<i; j++) |
if (qzer[j]) affrr(bn[j],m[i][j]); |
if (qzer[j]) affrr(bn[j], m[i][j]); |
else |
else |
{ |
{ |
gdivz(gscal(b[i],be[j]),bn[j],m[i][j]); |
gdivz(gscal(b[i],be[j]),bn[j], m[i][j]); |
p2=gsub(p2,gmul(m[i][j],be[j])); |
p2 = gsub(p2, gmul(m[i][j],be[j])); |
} |
} |
gaffect(p2,be[i]); affrr(sqscal(be[i]),bn[i]); |
gaffect(p2, be[i]); |
qzer[i]=quazero(bn[i]); avma=av1; |
affrr(sqscal(be[i]), bn[i]); |
|
qzer[i] = quazero(bn[i]); avma = av1; |
} |
} |
while (qzer[n]) |
while (qzer[n]) |
{ |
{ |
long e; |
long e; |
if (DEBUGLEVEL>9) |
if (DEBUGLEVEL>8) |
{ |
{ |
fprintferr("qzer[%ld]=%ld\n",n,qzer[n]); |
fprintferr("qzer[%ld]=%ld\n",n,qzer[n]); |
for (i1=1; i1<=n; i1++) |
for (k=1; k<=n; k++) |
for (i=1; i<i1; i++) output(m[i1][i]); |
for (i=1; i<k; i++) output(m[k][i]); |
} |
} |
em=bn[1]; j=1; |
em=bn[1]; j=1; |
for (i=2; i<n; i++) |
for (i=2; i<n; i++) |
Line 1547 lindep(GEN x, long prec) |
|
Line 1755 lindep(GEN x, long prec) |
|
p1=shiftr(bn[i],i); |
p1=shiftr(bn[i],i); |
if (cmprr(p1,em)>0) { em=p1; j=i; } |
if (cmprr(p1,em)>0) { em=p1; j=i; } |
} |
} |
i=j; i1=i+1; |
i=j; k=i+1; |
avma = av1; r = grndtoi(m[i1][i], &e); |
avma = av1; r = grndtoi(m[k][i], &e); |
if (e >= 0) err(precer,"lindep"); |
if (e >= 0) err(precer,"lindep"); |
r = negi(r); |
r = negi(r); |
p1 = ZV_lincomb(gun,r, b[i1],b[i]); |
p1 = ZV_lincomb(gun,r, b[k],b[i]); |
|
b[k] = b[i]; |
|
b[i] = p1; |
av1 = avma; |
av1 = avma; |
b[i1]=b[i]; b[i]=p1; f=addir(r,m[i1][i]); |
f = addir(r,m[k][i]); |
for (j=1; j<i; j++) |
for (j=1; j<i; j++) |
if (!qzer[j]) |
if (!qzer[j]) |
{ |
{ |
p1=mpadd(m[i1][j],mulir(r,m[i][j])); |
p1 = mpadd(m[k][j], mulir(r,m[i][j])); |
affrr(m[i][j],m[i1][j]); mpaff(p1,m[i][j]); |
affrr(m[i][j],m[k][j]); |
|
mpaff(p1,m[i][j]); |
} |
} |
c1=addrr(bn[i1],mulrr(gsqr(f),bn[i])); |
c1 = addrr(bn[k], mulrr(gsqr(f),bn[i])); |
if (!quazero(c1)) |
if (!quazero(c1)) |
{ |
{ |
c2=divrr(mulrr(bn[i],f),c1); affrr(c2,m[i1][i]); |
c2 = divrr(mulrr(bn[i],f),c1); affrr(c2, m[k][i]); |
c3=divrr(bn[i1],c1); mulrrz(c3,bn[i],bn[i1]); |
c3 = divrr(bn[k],c1); |
affrr(c1,bn[i]); qzer[i1]=quazero(bn[i1]); qzer[i]=0; |
mulrrz(c3,bn[i], bn[k]); qzer[k] = quazero(bn[k]); |
|
affrr(c1, bn[i]); qzer[i] = 0; |
for (j=i+2; j<=n; j++) |
for (j=i+2; j<=n; j++) |
{ |
{ |
p1=addrr(mulrr(m[j][i1],c3),mulrr(m[j][i],c2)); |
p1 = addrr(mulrr(m[j][k],c3), mulrr(m[j][i],c2)); |
subrrz(m[j][i],mulrr(f,m[j][i1]), m[j][i1]); |
subrrz(m[j][i],mulrr(f,m[j][k]), m[j][k]); |
affrr(p1,m[j][i]); |
affrr(p1, m[j][i]); |
} |
} |
} |
} |
else |
else |
{ |
{ |
qzer[i1]=qzer[i]; qzer[i]=1; |
affrr(bn[i], bn[k]); qzer[k] = qzer[i]; |
affrr(bn[i],bn[i1]); affrr(c1,bn[i]); |
affrr(c1, bn[i]); qzer[i] = 1; |
for (j=i+2; j<=n; j++) affrr(m[j][i],m[j][i1]); |
for (j=i+2; j<=n; j++) affrr(m[j][i], m[j][k]); |
} |
} |
if (low_stack(lim, stack_lim(av,1))) |
if (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
Line 1586 lindep(GEN x, long prec) |
|
Line 1798 lindep(GEN x, long prec) |
|
av1 = avma; |
av1 = avma; |
} |
} |
} |
} |
p1=cgetg(lx,t_COL); p1[n]=un; for (i=1; i<n; i++) p1[i]=zero; |
p1 = cgetg(lx,t_COL); p1[n] = un; for (i=1; i<n; i++) p1[i] = zero; |
p2 = (GEN)b; p2[0] = evaltyp(t_MAT) | evallg(lx); |
p2 = (GEN)b; p2[0] = evaltyp(t_MAT) | evallg(lx); |
p2=gauss(gtrans(p2),p1); tetpil=avma; |
p2 = gauss(gtrans_i(p2),p1); |
return gerepile(av,tetpil,gtrans(p2)); |
return gerepileupto(av, gtrans(p2)); |
} |
} |
|
|
/* x is a vector of elts of a p-adic field */ |
/* PSLQ Programs */ |
GEN |
|
plindep(GEN x) |
typedef struct { |
|
long vmind, t12, t1234, reda, fin; |
|
long ct; |
|
} pslq_timer; |
|
|
|
/* WARNING: for double ** matrices, A[i] is the i-th ROW of A */ |
|
typedef struct { |
|
double *y, **H, **A, **B; |
|
double *W; /* scratch space */ |
|
long n; |
|
pslq_timer *T; |
|
} pslqL2_M; |
|
|
|
typedef struct { |
|
GEN y, H, A, B; |
|
long n, EXP; |
|
int flreal; |
|
pslq_timer *T; |
|
} pslq_M; |
|
|
|
void |
|
init_dalloc() |
|
{ /* correct alignment for dalloc */ |
|
avma -= avma % sizeof(double); |
|
if (avma < bot) err(errpile); |
|
} |
|
|
|
double * |
|
dalloc(size_t n) |
{ |
{ |
long av = avma,i,j, prec = VERYBIGINT, lx = lg(x)-1, ly,v; |
if (avma - bot < n) err(errpile); |
GEN p = NULL, pn,p1,m,a; |
avma -= n; return (double*)avma; |
|
} |
|
|
if (lx < 2) return cgetg(1,t_VEC); |
static double |
for (i=1; i<=lx; i++) |
conjd(double x) { return x; } |
|
|
|
static double |
|
sqrd(double a) { return a*a; } |
|
|
|
static void |
|
redall(pslq_M *M, long i, long jsup) |
|
{ |
|
long j, k, n = M->n; |
|
GEN t,b; |
|
GEN A = M->A, B = M->B, H = M->H, y = M->y; |
|
const GEN p1 = (GEN)B[i]; |
|
|
|
for (j=jsup; j>=1; j--) |
{ |
{ |
p1 = (GEN)x[i]; |
/* FIXME: use grndtoi */ |
if (typ(p1) != t_PADIC) continue; |
t = ground(gdiv(gcoeff(H,i,j), gcoeff(H,j,j))); |
|
if (gcmp0(t)) continue; |
|
|
j = precp(p1); if (j < prec) prec = j; |
b = (GEN)B[j]; |
if (!p) p = (GEN)p1[2]; |
y[j] = ladd((GEN)y[j], gmul(t,(GEN)y[i])); |
else if (!egalii(p, (GEN)p1[2])) |
for (k=1; k<=j; k++) |
err(talker,"inconsistent primes in plindep"); |
coeff(H,i,k) = lsub(gcoeff(H,i,k), gmul(t,gcoeff(H,j,k))); |
|
for (k=1; k<=n; k++) |
|
{ |
|
coeff(A,i,k) = lsub(gcoeff(A,i,k), gmul(t,gcoeff(A,j,k))); |
|
b[k] = ladd((GEN)b[k], gmul(t,(GEN)p1[k])); |
|
} |
} |
} |
if (!p) err(talker,"not a p-adic vector in plindep"); |
} |
v = ggval(x,p); pn = gpowgs(p,prec); |
|
if (v != 0) x = gmul(x, gpowgs(p, -v)); |
|
x = lift_intern(gmul(x, gmodulcp(gun, pn))); |
|
|
|
ly = 2*lx - 1; |
static void |
m = cgetg(ly+1,t_MAT); |
redallbar(pslqL2_M *Mbar, long i, long jsup) |
for (j=1; j<=ly; j++) |
{ |
|
long j, k, n = Mbar->n; |
|
double t; |
|
double *hi = Mbar->H[i], *ai = Mbar->A[i], *hj, *aj; |
|
|
|
#if 0 |
|
fprintferr("%ld:\n==\n",i); |
|
#endif |
|
for (j=jsup; j>=1; j--) |
{ |
{ |
p1 = cgetg(lx+1, t_COL); m[j] = (long)p1; |
hj = Mbar->H[j]; |
for (i=1; i<=lx; i++) p1[i] = zero; |
t = floor(0.5 + hi[j] / hj[j]); |
|
if (!t) continue; |
|
#if 0 |
|
fprintferr("%15.15e ",t); |
|
#endif |
|
aj = Mbar->A[j]; |
|
|
|
Mbar->y[j] += t * Mbar->y[i]; |
|
for (k=1; k<=j; k++) hi[k] -= t * hj[k]; |
|
for (k=1; k<=n; k++) { |
|
ai[k] -= t * aj[k]; |
|
Mbar->B[k][j] += t * Mbar->B[k][i]; |
|
} |
|
#if 0 |
|
fprintferr(" %ld:\n",j); dprintmat(Mbar->H,n,n-1); |
|
#endif |
} |
} |
a = negi((GEN)x[1]); |
} |
for (i=1; i<lx; i++) |
|
|
static long |
|
vecabsminind(GEN v) |
|
{ |
|
long l = lg(v), m = 1, i; |
|
GEN t, la = mpabs((GEN)v[1]); |
|
for (i=2; i<l; i++) |
{ |
{ |
coeff(m,1+i,i) = (long)a; |
t = mpabs((GEN)v[i]); |
coeff(m,1 ,i) = x[i+1]; |
if (mpcmp(t,la) < 0) { la = t; m = i; } |
} |
} |
for (i=1; i<=lx; i++) coeff(m,i,lx-1+i) = (long)pn; |
return m; |
p1 = lllint(m); |
|
return gerepileupto(av, gmul(m, (GEN)p1[1])); |
|
} |
} |
|
|
GEN |
static long |
algdep0(GEN x, long n, long bit, long prec) |
vecmaxind(GEN v) |
{ |
{ |
long tx=typ(x),av,i,k; |
long l = lg(v), m = 1, i; |
GEN y,p1; |
GEN t, la = (GEN)v[1]; |
|
for (i=2; i<l; i++) |
|
{ |
|
t = (GEN)v[i]; |
|
if (mpcmp(t,la) > 0) { la = t; m = i; } |
|
} |
|
return m; |
|
} |
|
|
if (! is_scalar_t(tx)) err(typeer,"algdep0"); |
static long |
if (tx==t_POLMOD) { y=forcecopy((GEN)x[1]); setvarn(y,0); return y; } |
vecmaxindbar(double *v, long n) |
if (gcmp0(x)) return gzero; |
{ |
if (!n) return gun; |
long m = 1, i; |
|
double la = v[1]; |
|
for (i=2; i<=n; i++) |
|
if (v[i] > la) { la = v[i]; m = i; } |
|
return m; |
|
} |
|
|
av=avma; p1=cgetg(n+2,t_COL); |
static GEN |
p1[1] = un; |
maxnorml2(pslq_M *M) |
p1[2] = (long)x; /* n >= 1 */ |
{ |
for (i=3; i<=n+1; i++) p1[i]=lmul((GEN)p1[i-1],x); |
long n = M->n, i, j; |
if (typ(x) == t_PADIC) |
GEN ma = gzero, s; |
p1 = plindep(p1); |
|
else |
|
p1 = bit? lindep2(p1,bit): lindep(p1,prec); |
|
if (lg(p1) < 2) |
|
err(talker,"higher degree than expected in algdep"); |
|
|
|
y=cgetg(n+3,t_POL); |
for (i=1; i<=n; i++) |
y[1] = evalsigne(1) | evalvarn(0); |
{ |
k=1; while (gcmp0((GEN)p1[k])) k++; |
s = gzero; |
for (i=0; i<=n+1-k; i++) y[i+2] = p1[k+i]; |
for (j=1; j<n; j++) s = gadd(s, gnorm(gcoeff(M->H,i,j))); |
normalizepol_i(y, n+4-k); |
ma = gmax(ma, s); |
y = (gsigne(leading_term(y)) > 0)? gcopy(y): gneg(y); |
} |
return gerepileupto(av,y); |
return gsqrt(ma, DEFAULTPREC); |
} |
} |
|
|
GEN |
static void |
algdep2(GEN x, long n, long bit) |
init_timer(pslq_timer *T) |
{ |
{ |
return algdep0(x,n,bit,0); |
T->vmind = T->t12 = T->t1234 = T->reda = T->fin = T->ct = 0; |
} |
} |
|
|
GEN |
static int |
algdep(GEN x, long n, long prec) |
init_pslq(pslq_M *M, GEN x, long prec) |
{ |
{ |
return algdep0(x,n,0,prec); |
long lx = lg(x), n = lx-1, i, j, k; |
} |
GEN s1, s, sinv; |
|
|
/********************************************************************/ |
M->EXP = - bit_accuracy(prec) + 2*n; |
/** **/ |
M->flreal = (gexpo(gimag(x)) < M->EXP); |
/** INTEGRAL KERNEL (LLL REDUCED) **/ |
if (!M->flreal) |
/** **/ |
return 0; /* FIXME */ |
/********************************************************************/ |
else |
|
x = greal(x); |
|
|
GEN |
if (DEBUGLEVEL>=3) { (void)timer(); init_timer(M->T); } |
matkerint0(GEN x, long flag) |
x = gmul(x, realun(prec)); settyp(x,t_VEC); |
{ |
M->n = n; |
switch(flag) |
M->A = idmat(n); |
|
M->B = idmat(n); |
|
s1 = cgetg(lx,t_VEC); s1[n] = lnorm((GEN)x[n]); |
|
s = cgetg(lx,t_VEC); s[n] = (long)gabs((GEN)x[n],prec); |
|
for (k=n-1; k>=1; k--) |
{ |
{ |
case 0: return kerint(x); |
s1[k] = ladd((GEN)s1[k+1], gnorm((GEN)x[k])); |
case 1: return kerint1(x); |
s[k] = (long)gsqrt((GEN)s1[k], prec); |
case 2: return kerint2(x); |
|
default: err(flagerr,"matkerint"); |
|
} |
} |
return NULL; /* not reached */ |
sinv = ginv((GEN)s[1]); |
|
s = gmul(sinv,s); |
|
M->y = gmul(sinv, x); |
|
M->H = cgetg(n,t_MAT); |
|
for (j=1; j<n; j++) |
|
{ |
|
GEN d, c = cgetg(lx,t_COL); |
|
|
|
M->H[j] = (long)c; |
|
for (i=1; i<j; i++) c[i] = zero; |
|
c[j] = ldiv((GEN)s[j+1],(GEN)s[j]); |
|
d = gneg( gdiv((GEN)M->y[j], gmul((GEN)s[j],(GEN)s[j+1]) )); |
|
for (i=j+1; i<=n; i++) c[i] = lmul(gconj((GEN)M->y[i]), d); |
|
} |
|
for (i=2; i<=n; i++) redall(M, i, i-1); |
|
return 1; |
} |
} |
|
|
GEN |
static void |
kerint1(GEN x) |
SWAP(pslq_M *M, long m) |
{ |
{ |
long av=avma,tetpil; |
long j, n = M->n; |
GEN p1,p2; |
swap(M->y[m], M->y[m+1]); |
|
swap(M->B[m], M->B[m+1]); |
p1=matrixqz3(ker(x)); p2=lllint(p1); tetpil=avma; |
for (j=1; j<=n; j++) swap(coeff(M->A,m,j), coeff(M->A,m+1,j)); |
return gerepile(av,tetpil,gmul(p1,p2)); |
for (j=1; j<n; j++) swap(coeff(M->H,m,j), coeff(M->H,m+1,j)); |
} |
} |
|
|
GEN |
#define dswap(x,y) { double _t=x; x=y; y=_t; } |
kerint2(GEN x) |
#define ddswap(x,y) { double* _t=x; x=y; y=_t; } |
{ |
|
long lx=lg(x), i,j,av,av1; |
|
GEN g,p1; |
|
|
|
if (typ(x) != t_MAT) err(typeer,"kerint2"); |
static void |
av=avma; g=cgetg(lx,t_MAT); |
SWAPbar(pslqL2_M *M, long m) |
for (j=1; j<lx; j++) g[j]=lgetg(lx,t_COL); |
{ |
for (i=1; i<lx; i++) |
long j, n = M->n; |
for (j=1; j<=i; j++) |
dswap(M->y[m], M->y[m+1]); |
coeff(g,i,j) = coeff(g,j,i) = (long)gscal((GEN)x[i],(GEN)x[j]); |
ddswap(M->A[m], M->A[m+1]); |
g=lllgramall0(g,lll_KER); p1=lllint(g); |
ddswap(M->H[m], M->H[m+1]); |
av1=avma; return gerepile(av,av1,gmul(g,p1)); |
for (j=1; j<=n; j++) dswap(M->B[j][m], M->B[j][m+1]); |
} |
} |
|
|
static GEN |
static GEN |
lllall0(GEN x, long flag) |
one_step_gen(pslq_M *M, GEN tabga, long prec) |
{ |
{ |
long av0=avma,av,tetpil,lx=lg(x),i,j,k,l,n,lim,kmax; |
GEN H = M->H, p1, t0, tinv, t1,t2,t3,t4; |
GEN u,B,L,q,r,h,la,p1,p2,p4,fl; |
long n = M->n, i, m; |
|
|
if (typ(x) != t_MAT) err(typeer,"lllall0"); |
p1 = cgetg(n,t_VEC); |
n=lx-1; if (n<=1) return lllall_trivial(x,n, flag | lll_GRAM); |
for (i=1; i<n; i++) p1[i] = lmul((GEN)tabga[i], gabs(gcoeff(H,i,i),prec)); |
fl = new_chunk(lx); |
m = vecmaxind(p1); |
|
if (DEBUGLEVEL>=4) M->T->vmind += timer(); |
av=avma; lim=stack_lim(av,1); x=dummycopy(x); |
SWAP(M, m); |
B=gscalcol(gun, lx); |
if (m <= n-2) |
L=cgetg(lx,t_MAT); |
|
for (k=lg(x[1]),j=1; j<lx; j++) |
|
{ |
{ |
for (i=1; i<k; i++) |
t0 = gadd(gnorm(gcoeff(H,m,m)), gnorm(gcoeff(H,m,m+1))); |
if (typ(gcoeff(x,i,j))!=t_INT) err(lllger4); |
tinv = ginv(gsqrt(t0, prec)); |
fl[j] = 0; L[j] = (long)zerocol(n); |
t1 = gmul(tinv, gcoeff(H,m,m)); |
} |
t2 = gmul(tinv, gcoeff(H,m,m+1)); |
k=2; h=idmat(n); kmax=1; |
if (DEBUGLEVEL>=4) M->T->t12 += timer(); |
u=sqscali((GEN)x[1]); |
for (i=m; i<=n; i++) |
if (signe(u)) { B[2]=(long)u; coeff(L,1,1)=un; fl[1]=1; } else B[2]=un; |
|
for(;;) |
|
{ |
|
if (k > kmax) |
|
{ |
{ |
kmax = k; |
t3 = gcoeff(H,i,m); |
for (j=1; j<=k; j++) |
t4 = gcoeff(H,i,m+1); |
{ |
if (M->flreal) |
if (j==k || fl[j]) |
coeff(H,i,m) = ladd(gmul(t1,t3), gmul(t2,t4)); |
{ |
|
long av1 = avma; |
|
u=gscali((GEN)x[k],(GEN)x[j]); |
|
for (i=1; i<j; i++) |
|
if (fl[i]) |
|
u = divii(subii(mulii((GEN)B[i+1],u), |
|
mulii(gcoeff(L,k,i),gcoeff(L,j,i))), |
|
(GEN)B[i]); |
|
u = gerepileuptoint(av1, u); |
|
if (j<k) coeff(L,k,j)=(long)u; |
|
else |
|
{ |
|
if (signe(u)) { B[k+1]=(long)u; coeff(L,k,k)=un; fl[k]=1; } |
|
else { B[k+1]=B[k]; fl[k]=0; } |
|
} |
|
} |
|
} |
|
} |
|
if (fl[k-1] && !fl[k]) |
|
{ |
|
u = shifti(gcoeff(L,k,k-1),1); |
|
if (absi_cmp(u, (GEN)B[k]) > 0) |
|
{ |
|
q = truedvmdii(addii(u,(GEN)B[k]),shifti((GEN)B[k],1), NULL); |
|
r = negi(q); |
|
h[k] = (long)ZV_lincomb(gun,r, (GEN)h[k],(GEN)h[k-1]); |
|
x[k] = (long)ZV_lincomb(gun,r, (GEN)x[k],(GEN)x[k-1]); |
|
coeff(L,k,k-1)=laddii(gcoeff(L,k,k-1),mulii(r,(GEN)B[k])); |
|
for (i=1; i<k-1; i++) |
|
coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(r,gcoeff(L,k-1,i))); |
|
} |
|
la=gcoeff(L,k,k-1); p4=sqri(la); |
|
swap(h[k-1], h[k]); |
|
swap(x[k-1], x[k]); |
|
for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j)); |
|
if (signe(la)) |
|
{ |
|
p2=(GEN)B[k]; p1=divii(p4,p2); |
|
B[k+1]=B[k]=(long)p1; |
|
for (i=k+1; i<=kmax; i++) |
|
coeff(L,i,k-1)=ldivii(mulii(la,gcoeff(L,i,k-1)),p2); |
|
for (j=k+1; j<kmax; j++) |
|
for (i=j+1; i<=kmax; i++) |
|
coeff(L,i,j)=ldivii(mulii((GEN)p1,gcoeff(L,i,j)),p2); |
|
for (i=k+2; i<=kmax+1; i++) |
|
B[i]=ldivii(mulii((GEN)p1,(GEN)B[i]),p2); |
|
} |
|
else |
else |
{ |
coeff(H,i,m) = ladd(gmul(gconj(t1),t3), gmul(gconj(t2),t4)); |
for (i=k+1; i<=kmax; i++) |
coeff(H,i,m+1) = lsub(gmul(t1,t4), gmul(t2,t3)); |
{ coeff(L,i,k)=coeff(L,i,k-1); coeff(L,i,k-1)=zero; } |
|
B[k]=B[k-1]; fl[k]=1; fl[k-1]=0; |
|
} |
|
if (k>2) k--; |
|
} |
} |
else |
if (DEBUGLEVEL>=4) M->T->t1234 += timer(); |
{ |
} |
for (l=k-1; l>=1; l--) |
for (i=1; i<=n-1; i++) |
{ |
if (gexpo(gcoeff(H,i,i)) <= M->EXP) { |
u = shifti(gcoeff(L,k,l),1); |
m = vecabsminind(M->y); return (GEN)M->B[m]; |
if (absi_cmp(u,(GEN)B[l+1]) > 0) |
|
{ |
|
q = truedvmdii(addii(u,(GEN)B[l+1]),shifti((GEN)B[l+1],1), NULL); |
|
r = negi(q); |
|
h[k] = (long)ZV_lincomb(gun,r,(GEN)h[k],(GEN)h[l]); |
|
x[k] = (long)ZV_lincomb(gun,r,(GEN)x[k],(GEN)x[l]); |
|
coeff(L,k,l)=laddii(gcoeff(L,k,l),mulii(r,(GEN)B[l+1])); |
|
for (i=1; i<l; i++) |
|
coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(r,gcoeff(L,l,i))); |
|
} |
|
} |
|
if (++k > n) break; |
|
} |
} |
if (low_stack(lim, stack_lim(av,1))) |
for (i=m+1; i<=n; i++) redall(M, i, min(i-1,m+1)); |
|
|
|
if (DEBUGLEVEL>=4) M->T->reda += timer(); |
|
if (gexpo(M->A) >= -M->EXP) return ginv(maxnorml2(M)); |
|
m = vecabsminind(M->y); |
|
if (gexpo((GEN)M->y[m]) <= M->EXP) return (GEN)M->B[m]; |
|
|
|
if (DEBUGLEVEL>=3) |
|
{ |
|
if (DEBUGLEVEL>=4) M->T->fin += timer(); |
|
M->T->ct++; |
|
if ((M->T->ct&0xff) == 0) |
{ |
{ |
GEN *gptr[4]; |
if (DEBUGLEVEL == 3) |
if(DEBUGMEM>1) err(warnmem,"lllall0"); |
fprintferr("time for ct = %ld : %ld\n",M->T->ct,timer()); |
gptr[0]=&B; gptr[1]=&L; gptr[2]=&h; |
else |
gptr[3]=&x; gerepilemany(av,gptr,4); |
fprintferr("time [max,t12,loop,reds,fin] = [%ld, %ld, %ld, %ld, %ld]\n", |
|
M->T->vmind, M->T->t12, M->T->t1234, M->T->reda, M->T->fin); |
} |
} |
} |
} |
tetpil=avma; |
return NULL; /* nothing interesting */ |
return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n)); |
|
} |
} |
|
|
GEN |
static GEN |
kerint(GEN x) |
get_tabga(int flreal, long n, long prec) |
{ |
{ |
long av=avma,av1; |
GEN ga = mpsqrt( flreal? divrs(stor(4, prec), 3): stor(2, prec) ); |
GEN g,p1; |
GEN tabga = cgetg(n,t_VEC); |
|
long i; |
g=lllall0(x,lll_KER); if (lg(g)==1) return g; |
tabga[1] = (long)ga; |
p1=lllint(g); av1=avma; |
for (i = 2; i < n; i++) tabga[i] = lmul((GEN)tabga[i-1],ga); |
return gerepile(av,av1,gmul(g,p1)); |
return tabga; |
} |
} |
|
|
/********************************************************************/ |
GEN |
/** **/ |
pslq(GEN x, long prec) |
/** POLRED & CO. **/ |
|
/** **/ |
|
/********************************************************************/ |
|
/* remove duplicate polynomials in y, updating a (same indexes), in place */ |
|
static long |
|
remove_duplicates(GEN y, GEN a) |
|
{ |
{ |
long k,i, nv = lg(y), av = avma; |
GEN tabga, p1; |
GEN z; |
long tx = typ(x); |
|
gpmem_t av0 = avma, lim = stack_lim(av0,1), av; |
|
pslq_M M; |
|
pslq_timer T; |
|
|
if (nv < 2) return nv; |
if (! is_vec_t(tx)) err(typeer,"pslq"); |
z = new_chunk(3); |
if (lg(x) <= 2) return cgetg(1,t_VEC); |
z[1] = (long)y; |
|
z[2] = (long)a; (void)sort_factor(z, gcmp); |
M.T = &T; init_pslq(&M, x, prec); |
for (k=1, i=2; i<nv; i++) |
if (!M.flreal) return lindep(x, prec); /* FIXME */ |
if (!gegal((GEN)y[i], (GEN)y[k])) |
|
|
tabga = get_tabga(M.flreal, M.n, prec); |
|
av = avma; |
|
if (DEBUGLEVEL>=3) printf("Initialization time = %ld\n",timer()); |
|
for (;;) |
|
{ |
|
if ((p1 = one_step_gen(&M, tabga, prec))) |
|
return gerepilecopy(av0, p1); |
|
|
|
if (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
k++; |
if(DEBUGMEM>1) err(warnmem,"pslq"); |
a[k] = a[i]; |
gerepileall(av,4,&M.y,&M.H,&M.A,&M.B); |
y[k] = y[i]; |
|
} |
} |
nv = k+1; setlg(a,nv); setlg(y,nv); |
} |
avma = av; return nv; |
|
} |
} |
|
|
/* in place; make sure second non-zero coeff is negative (choose x or -x) */ |
/* W de longueur n-1 */ |
int |
|
canon_pol(GEN z) |
static double |
|
dnorml2(double *W, long n, long row) |
{ |
{ |
long i,s; |
long i; |
|
double s = 0.; |
|
|
for (i=lgef(z)-2; i>=2; i-=2) |
for (i=row; i<n; i++) s += W[i]*W[i]; |
|
return s; |
|
} |
|
|
|
/* Hbar *= Pbar */ |
|
static void |
|
dmatmul(pslqL2_M *Mbar, double **Pbar, long row) |
|
{ |
|
const long n = Mbar->n; /* > row */ |
|
long i, j, k; |
|
double s, *W = Mbar->W, **H = Mbar->H; |
|
|
|
for (i = row; i <= n; i++) |
{ |
{ |
s = signe(z[i]); |
for (j = row; j < n; j++) |
if (s > 0) |
|
{ |
{ |
for (; i>=2; i-=2) z[i]=lnegi((GEN)z[i]); |
k = row; s = H[i][k] * Pbar[k][j]; |
return -1; |
for ( ; k < n; k++) s += H[i][k] * Pbar[k][j]; |
|
W[j] = s; |
} |
} |
if (s) return 1; |
for (j = row; j < n; j++) H[i][j] = W[j]; |
} |
} |
return 0; |
|
} |
} |
|
|
static GEN |
/* compute n-1 x n-1 matrix Pbar */ |
pols_for_polred(GEN x, GEN base, GEN LLLbase, GEN *pta, |
static void |
int (*check)(GEN, GEN), GEN arg) |
dmakep(pslqL2_M *Mbar, double **Pbar, long row) |
{ |
{ |
long i, v = varn(x), n = lg(base); |
long i, j, n = Mbar->n; |
GEN p1,p2,p3,y, a = cgetg(n,t_VEC); |
double pro, nc, *C = Mbar->H[row], *W = Mbar->W; |
|
|
for (i=1; i<n; i++) a[i] = lmul(base,(GEN)LLLbase[i]); |
nc = sqrt(dnorml2(C,n,row)); |
y=cgetg(n,t_VEC); |
W[row] = (C[row] < 0) ? C[row] - nc : C[row] + nc; |
for (i=1; i<n; i++) |
for (i=row; i<n; i++) W[i] = C[i]; |
|
pro = -2.0 / dnorml2(W, n, row); |
|
/* must have dnorml2(W,n,row) = 2*nc*(nc+fabs(C[1])) */ |
|
for (j=row; j<n; j++) |
{ |
{ |
if (DEBUGLEVEL > 2) { fprintferr("i = %ld\n",i); flusherr(); } |
for (i=j+1; i<n; i++) |
p1=(GEN)a[i]; |
Pbar[j][i] = Pbar[i][j] = pro * W[i] * W[j]; |
p1 = primitive_part(p1, &p3); |
Pbar[j][j] = 1.0 + pro * W[j] * W[j]; |
p1 = caractducos(x,p1,v); |
|
if (p3) p1 = ZX_rescale_pol(p1,p3); |
|
p2 = modulargcd(derivpol(p1),p1); |
|
p3 = leading_term(p2); if (!gcmp1(p3)) p2=gdiv(p2,p3); |
|
p1 = gdiv(p1,p2); |
|
if (canon_pol(p1) < 0 && pta) a[i] = (long)gneg_i((GEN)a[i]); |
|
y[i] = (long)p1; if (DEBUGLEVEL > 3) outerr(p1); |
|
if (check && check(arg, p1)) return p1; |
|
} |
} |
if (check) return NULL; /* no suitable polynomial found */ |
|
(void)remove_duplicates(y,a); |
|
if (pta) *pta = a; |
|
return y; |
|
} |
} |
|
|
GEN |
static void |
nf_get_T2(GEN base, GEN polr) |
dLQdec(pslqL2_M *Mbar, double **Pbar) |
{ |
{ |
long i,j, n = lg(base); |
long row, j, n = Mbar->n; |
GEN p1,p2=cgetg(n,t_MAT); |
|
|
|
for (i=1; i<n; i++) |
for (row=1; row<n; row++) |
{ |
{ |
p1=cgetg(n,t_COL); p2[i]=(long)p1; |
dmakep(Mbar, Pbar, row); |
for (j=1; j<n; j++) |
dmatmul(Mbar, Pbar, row); |
p1[j] = (long) poleval((GEN)base[i],(GEN)polr[j]); |
for (j=row+1; j<n; j++) Mbar->H[row][j] = 0.; |
} |
} |
return mulmat_real(gconj(gtrans(p2)),p2); |
|
} |
} |
|
|
/* compute Tr(w_i w_j) */ |
void |
static GEN |
dprintvec(double *V, long m) |
nf_get_T(GEN x, GEN w) |
|
{ |
{ |
long i,j,k, n = degpol(x); |
long i; |
GEN p1,p2,p3; |
fprintferr("["); |
GEN ptrace = cgetg(n+2,t_VEC); |
for (i=1; i<m; i++) fprintferr("%15.15e, ",V[i]); |
GEN den = cgetg(n+1,t_VEC); |
fprintferr("%15.15e]\n",V[m]); pariflush(); |
GEN T = cgetg(n+1,t_MAT); |
} |
|
|
ptrace[2]=lstoi(n); |
void |
for (k=2; k<=n; k++) |
dprintmat(double **M, long r, long c) |
{ /* cf polsym */ |
{ |
GEN y = x + (n-k+1); |
long i, j; |
p1 = mulsi(k-1,(GEN)y[2]); |
fprintferr("["); |
for (i=3; i<=k; i++) |
for (i=1; i<r; i++) |
p1 = addii(p1,mulii((GEN)y[i],(GEN)ptrace[i])); |
|
ptrace[i] = lnegi(p1); |
|
} |
|
w = dummycopy(w); |
|
for (i=1; i<=n; i++) |
|
{ |
{ |
den[i] = (long)denom(content((GEN)w[i])); |
for (j=1; j<c; j++) fprintferr("%15.15e, ",M[i][j]); |
w[i] = lmul((GEN)w[i],(GEN)den[i]); |
fprintferr("%15.15e;\n ",M[i][c]); |
} |
} |
|
for (j=1; j<c; j++) fprintferr("%15.15e, ",M[r][j]); |
for (i=1; i<=n; i++) |
fprintferr("%15.15e]\n",M[r][c]); pariflush(); |
{ |
|
p1=cgetg(n+1,t_COL); T[i]=(long)p1; |
|
for (j=1; j<i ; j++) p1[j] = coeff(T,i,j); |
|
for ( ; j<=n; j++) |
|
{ /* cf quicktrace */ |
|
p2 = gres(gmul((GEN)w[i],(GEN)w[j]),x); |
|
p3 = gzero; |
|
for (k=lgef(p2)-1; k>1; k--) |
|
p3 = addii(p3, mulii((GEN)p2[k],(GEN)ptrace[k])); |
|
p1[j]=(long)divii(p3, mulii((GEN)den[i],(GEN)den[j])); |
|
} |
|
} |
|
return T; |
|
} |
} |
|
|
/* Return the base change matrix giving the an LLL-reduced basis for the |
static long |
* maximal order of the nf given by x. Expressed in terms of the standard |
initializedoubles(pslqL2_M *Mbar, pslq_M *M, long prec) |
* HNF basis (as polynomials) given in base (ignored if x is an nf) |
|
*/ |
|
GEN |
|
LLL_nfbasis(GEN *ptx, GEN polr, GEN base, long prec) |
|
{ |
{ |
GEN T2,p1, x = *ptx; |
long i, j, n = Mbar->n; |
int totally_real,n,i; |
GEN ypro; |
|
gpmem_t av = avma; |
|
|
if (typ(x) != t_POL) |
ypro = gdiv(M->y, vecmax(gabs(M->y,prec))); |
|
for (i=1; i<=n; i++) |
{ |
{ |
p1=checknf(x); *ptx=x=(GEN)p1[1]; |
if (gexpo((GEN)ypro[i]) < -0x3ff) return 0; |
base=(GEN)p1[7]; n=degpol(x); |
Mbar->y[i] = rtodbl((GEN)ypro[i]); |
totally_real = !signe(gmael(p1,2,2)); |
|
T2=gmael(p1,5,3); if (totally_real) T2 = ground(T2); |
|
} |
} |
else |
avma = av; |
{ |
for (j=1; j<=n; j++) |
n=degpol(x); totally_real = (!prec || sturm(x)==n); |
for (i=1; i<=n; i++) |
if (typ(base) != t_VEC || lg(base)-1 != n) |
{ |
err(talker,"incorrect Zk basis in LLL_nfbasis"); |
if (i==j) Mbar->A[i][j] = Mbar->B[i][j] = 1.; |
if (!totally_real) |
else Mbar->A[i][j] = Mbar->B[i][j] = 0.; |
T2 = nf_get_T2(base,polr? polr: roots(x,prec)); |
if (j < n) |
else |
{ |
T2 = nf_get_T(x, base); |
GEN h = gcoeff(M->H,i,j); |
} |
if (!gcmp0(h) && labs(gexpo(h)) > 0x3ff) return 0; |
if (totally_real) return lllgramint(T2); |
Mbar->H[i][j] = rtodbl(h); |
for (i=1; ; i++) |
} |
{ |
} |
if ((p1 = lllgramintern(T2,100,1,prec))) return p1; |
return 1; |
if (i == MAXITERPOL) err(accurer,"allpolred"); |
|
prec=(prec<<1)-2; |
|
if (DEBUGLEVEL) err(warnprec,"allpolred",prec); |
|
T2=nf_get_T2(base,roots(x,prec)); |
|
} |
|
} |
} |
|
|
/* x can be a polynomial, but also an nf or a bnf */ |
/* T(arget) := S(ource) */ |
/* if check != NULL, return the first polynomial pol |
static void |
found such that check(arg, pol) != 0; return NULL |
storeprecdoubles(pslqL2_M *T, pslqL2_M *S) |
if there are no such pol */ |
|
static GEN |
|
allpolred0(GEN x, GEN *pta, long code, long prec, |
|
int (*check)(GEN,GEN), GEN arg) |
|
{ |
{ |
GEN y,p1, base = NULL, polr = NULL; |
long n = T->n, i, j; |
long av = avma; |
|
|
|
if (typ(x) == t_POL) |
for (i=1; i<=n; i++) |
{ |
{ |
if (!signe(x)) return gcopy(x); |
for (j=1; j<n; j++) |
check_pol_int(x,"polred"); |
|
if (!gcmp1(leading_term(x))) |
|
err(impl,"allpolred for nonmonic polynomials"); |
|
base = allbase4(x,code,&p1,NULL); /* p1 is junk */ |
|
} |
|
else |
|
{ |
|
long i = lg(x); |
|
if (typ(x) == t_VEC && i<=4 && i>=3 && typ(x[1])==t_POL) |
|
{ /* polynomial + integer basis */ |
|
base=(GEN)x[2]; x=(GEN)x[1]; |
|
} |
|
else |
|
{ |
{ |
x = checknf(x); |
T->H[i][j] = S->H[i][j]; |
base=(GEN)x[7]; x=(GEN)x[1]; |
T->A[i][j] = S->A[i][j]; |
|
T->B[i][j] = S->B[i][j]; |
} |
} |
|
T->A[i][n] = S->A[i][n]; |
|
T->B[i][n] = S->B[i][n]; |
|
T->y[i] = S->y[i]; |
} |
} |
p1 = LLL_nfbasis(&x,polr,base,prec); |
|
y = pols_for_polred(x,base,p1,pta,check,arg); |
|
if (check) |
|
{ |
|
if (y) return gerepileupto(av, y); |
|
avma = av; return NULL; |
|
} |
|
|
|
if (pta) |
|
{ |
|
GEN *gptr[2]; gptr[0]=&y; gptr[1]=pta; |
|
gerepilemany(av,gptr,2); return y; |
|
} |
|
return gerepileupto(av,y); |
|
} |
} |
|
|
static GEN |
static long |
allpolred(GEN x, GEN *pta, long code, long prec) |
checkentries(pslqL2_M *Mbar) |
{ |
{ |
return allpolred0(x,pta,code,prec,NULL,NULL); |
long n = Mbar->n, i, j; |
} |
double *p1, *p2; |
|
|
GEN |
for (i=1; i<=n; i++) |
polred0(GEN x,long flag, GEN p, long prec) |
|
{ |
|
GEN y; |
|
long smll = (flag & 1); |
|
|
|
if (p && !gcmp0(p)) smll=(long) p; /* factored polred */ |
|
if (flag & 2) /* polred2 */ |
|
{ |
{ |
y=cgetg(3,t_MAT); |
if (expodb(Mbar->y[i]) < -46) return 0; |
y[2]=(long)allpolred(x,(GEN*)(y+1),smll,prec); |
p1 = Mbar->A[i]; |
return y; |
p2 = Mbar->B[i]; |
|
for (j=1; j<=n; j++) |
|
if (expodb(p1[j]) > 43 || expodb(p2[j]) > 43) return 0; |
} |
} |
return allpolred(x,NULL,smll,prec); |
return 1; |
} |
} |
|
|
GEN |
static long |
ordred(GEN x, long prec) |
applybar(pslq_M *M, pslqL2_M *Mbar, GEN Abargen, GEN Bbargen) |
{ |
{ |
GEN base,y; |
long n = Mbar->n, i, j; |
long n=degpol(x),i,av=avma,v = varn(x); |
double *p1, *p2; |
|
|
if (typ(x) != t_POL) err(typeer,"ordred"); |
|
if (!signe(x)) return gcopy(x); |
|
if (!gcmp1((GEN)x[n+2])) err(impl,"ordred for nonmonic polynomials"); |
|
|
|
base = cgetg(n+1,t_VEC); /* power basis */ |
|
for (i=1; i<=n; i++) |
for (i=1; i<=n; i++) |
base[i] = lpowgs(polx[v],i-1); |
{ |
y = cgetg(3,t_VEC); |
p1 = Mbar->A[i]; |
y[1] = (long)x; |
p2 = Mbar->B[i]; |
y[2] = (long)base; |
for (j=1; j<=n; j++) |
return gerepileupto(av, allpolred(y,NULL,0,prec)); |
{ |
|
if (expodb(p1[j]) >= 52 || expodb(p2[j]) >= 52) return 0; |
|
coeff(Abargen,i,j) = (long)mpent(dbltor(p1[j])); |
|
coeff(Bbargen,i,j) = (long)mpent(dbltor(p2[j])); |
|
} |
|
} |
|
M->y = gmul(M->y, Bbargen); |
|
M->B = gmul(M->B, Bbargen); |
|
M->A = gmul(Abargen, M->A); |
|
M->H = gmul(Abargen, M->H); return 1; |
} |
} |
|
|
GEN |
static GEN |
sum(GEN v, long a, long b) |
checkend(pslq_M *M, long prec) |
{ |
{ |
GEN p; |
long i, m, n = M->n; |
long i; |
|
if (a > b) return gzero; |
|
p = (GEN)v[a]; |
|
for (i=a+1; i<=b; i++) p = gadd(p, (GEN)v[i]); |
|
return p; |
|
} |
|
|
|
GEN |
for (i=1; i<=n-1; i++) |
T2_from_embed_norm(GEN x, long r1) |
if (gexpo(gcoeff(M->H,i,i)) <= M->EXP) |
{ |
{ |
GEN p = sum(x, 1, r1); |
m = vecabsminind(M->y); |
GEN q = sum(x, r1+1, lg(x)-1); |
return (GEN)M->B[m]; |
if (q != gzero) p = gadd(p, gmul2n(q,1)); |
} |
return p; |
if (gexpo(M->A) >= -M->EXP) |
|
return ginv( maxnorml2(M) ); |
|
m = vecabsminind(M->y); |
|
if (gexpo((GEN)M->y[m]) <= M->EXP) return (GEN)M->B[m]; |
|
return NULL; |
} |
} |
|
|
GEN |
GEN |
T2_from_embed(GEN x, long r1) |
pslqL2(GEN x, long prec) |
{ |
{ |
return T2_from_embed_norm(gnorm(x), r1); |
GEN Abargen, Bbargen, tabga, p1; |
} |
long lx = lg(x), tx = typ(x), n = lx-1, i, m, ctpro, flreal, flit; |
|
gpmem_t av0 = avma, lim = stack_lim(av0,1), av; |
|
double *tabgabar, gabar, tinvbar, t1bar, t2bar, t3bar, t4bar; |
|
double **Pbar, **Abar, **Bbar, **Hbar, *ybar; |
|
pslqL2_M Mbar, Mbarst; |
|
pslq_M M; |
|
pslq_timer T; |
|
|
/* return T2 norm of the polynomial defining nf */ |
if (! is_vec_t(tx)) err(typeer,"pslq"); |
static GEN |
if (n <= 1) return cgetg(1,t_COL); |
get_Bnf(GEN nf) |
M.T = &T; init_pslq(&M, x, prec); |
{ |
if (!M.flreal) return lindep(x, prec); |
return T2_from_embed((GEN)nf[6], nf_get_r1(nf)); |
|
} |
|
|
|
/* characteristic pol of x */ |
flreal = M.flreal; |
static GEN |
tabga = get_tabga(flreal, n, prec); |
get_polchar(GEN data, GEN x) |
Abargen = idmat(n); |
{ |
Bbargen = idmat(n); |
GEN basis_embed = (GEN)data[1]; |
|
GEN g = gmul(basis_embed,x); |
|
return ground(roots_to_pol_r1r2(g, data[0], 0)); |
|
} |
|
|
|
/* return a defining polynomial for Q(x) */ |
Mbarst.n = Mbar.n = n; |
static GEN |
Mbar.A = Abar = (double**)new_chunk(n+1); |
get_polmin(GEN data, GEN x) |
Mbar.B = Bbar = (double**)new_chunk(n+1); |
{ |
Mbar.H = Hbar = (double**)new_chunk(n+1); |
GEN g = get_polchar(data,x); |
Mbarst.A = (double**)new_chunk(n+1); |
GEN h = modulargcd(g,derivpol(g)); |
Mbarst.B = (double**)new_chunk(n+1); |
if (lgef(h) > 3) g = gdivexact(g,h); |
Mbarst.H = (double**)new_chunk(n+1); |
return g; |
Pbar = (double**)new_chunk(n); |
} |
|
|
|
/* does x generate the correct field ? */ |
tabgabar = dalloc((n+1)*sizeof(double)); |
static GEN |
Mbar.y = ybar = dalloc((n+1)*sizeof(double)); |
chk_gen(GEN data, GEN x) |
Mbarst.y = dalloc((n+1)*sizeof(double)); |
{ |
|
long av = avma; |
|
GEN g = get_polchar(data,x); |
|
if (lgef(modulargcd(g,derivpol(g))) > 3) { avma=av; return NULL; } |
|
if (DEBUGLEVEL>3) fprintferr(" generator: %Z\n",g); |
|
return g; |
|
} |
|
|
|
/* mat = base change matrix, gram = LLL-reduced T2 matrix */ |
Mbar.W = dalloc((n+1)*sizeof(double)); |
static GEN |
for (i=1; i< n; i++) Pbar[i] = dalloc((n+1)*sizeof(double)); |
chk_gen_init(FP_chk_fun *chk, GEN nf, GEN gram, GEN mat, long *ptprec) |
for (i=1; i<=n; i++) Abar[i] = dalloc((n+1)*sizeof(double)); |
{ |
for (i=1; i<=n; i++) Bbar[i] = dalloc((n+1)*sizeof(double)); |
GEN P,bound,prev,x,B,data, M = gmael(nf,5,1); |
for (i=1; i<=n; i++) Hbar[i] = dalloc(n*sizeof(double)); |
long N = lg(nf[7]), n = N-1,i,prec,prec2; |
for (i=1; i<=n; i++) Mbarst.A[i] = dalloc((n+1)*sizeof(double)); |
int skipfirst = 1; /* [1,0...0] --> x rational */ |
for (i=1; i<=n; i++) Mbarst.B[i] = dalloc((n+1)*sizeof(double)); |
|
for (i=1; i<=n; i++) Mbarst.H[i] = dalloc(n*sizeof(double)); |
|
|
/* data[0] = r1 |
gabar = gtodouble((GEN)tabga[1]); tabgabar[1] = gabar; |
* data[1] = embeddings of LLL-reduced Zk basis |
for (i=2; i<n; i++) tabgabar[i] = tabgabar[i-1]*gabar; |
* data[2] = LLL reduced Zk basis (in M_n(Z)) */ |
|
data = new_chunk(3); |
|
data[0] = itos(gmael(nf,2,1)); |
|
data[1] = lmul(M, mat); |
|
data[2] = lmul((GEN)nf[7],mat); |
|
chk->data = data; |
|
|
|
x = cgetg(N,t_COL); |
av = avma; |
bound = get_Bnf(nf); prev = NULL; |
if (DEBUGLEVEL>=3) printf("Initialization time = %ld\n",timer()); |
for (i=1; i<N; i++) x[i]=zero; |
RESTART: |
for (i=2; i<N; i++) |
flit = initializedoubles(&Mbar, &M, prec); |
|
storeprecdoubles(&Mbarst, &Mbar); |
|
if (flit) dLQdec(&Mbar, Pbar); |
|
ctpro = 0; |
|
for (;;) |
{ |
{ |
x[i] = un; |
if (low_stack(lim, stack_lim(av,1))) |
P = get_polmin(data,x); |
|
x[i] = zero; |
|
if (degpol(P) == n) |
|
{ |
{ |
B = gcoeff(gram,i,i); |
if(DEBUGMEM>1) err(warnmem,"pslq"); |
if (gcmp(B,bound) < 0) bound = B ; |
gerepileall(av,4,&M.y,&M.H,&M.A,&M.B); |
} |
} |
else |
if (flit) |
{ |
{ |
if (DEBUGLEVEL>2) fprintferr("chk_gen_init: subfield %Z\n",P); |
ctpro++; |
if (skipfirst == i-1) |
for (i=1; i<n; i++) Mbar.W[i] = tabgabar[i]*fabs(Hbar[i][i]); |
|
m = vecmaxindbar(Mbar.W, n-1); |
|
SWAPbar(&Mbar, m); |
|
if (m <= n-2) |
{ |
{ |
if (prev && !gegal(prev,P)) |
tinvbar = 1.0 / sqrt(sqrd(Hbar[m][m]) + sqrd(Hbar[m][m+1])); |
|
t1bar = tinvbar*Hbar[m][m]; |
|
t2bar = tinvbar*Hbar[m][m+1]; |
|
if (DEBUGLEVEL>=4) T.t12 += timer(); |
|
for (i=m; i<=n; i++) |
|
{ |
|
t3bar = Hbar[i][m]; |
|
t4bar = Hbar[i][m+1]; |
|
if (flreal) |
|
Hbar[i][m] = t1bar*t3bar + t2bar*t4bar; |
|
else |
|
Hbar[i][m] = conjd(t1bar)*t3bar + conjd(t2bar)*t4bar; |
|
Hbar[i][m+1] = t1bar*t4bar - t2bar*t3bar; |
|
} |
|
if (DEBUGLEVEL>=4) T.t1234 += timer(); |
|
} |
|
|
|
flit = checkentries(&Mbar); |
|
if (flit) |
|
{ |
|
storeprecdoubles(&Mbarst, &Mbar); |
|
for (i=m+1; i<=n; i++) redallbar(&Mbar, i, min(i-1,m+1)); |
|
} |
|
else |
|
{ |
|
if (applybar(&M, &Mbar, Abargen,Bbargen)) |
|
{ |
|
if ( (p1 = checkend(&M,prec)) ) return gerepilecopy(av0, p1); |
|
goto RESTART; |
|
} |
|
else |
{ |
{ |
if (degpol(prev) * degpol(P) > 32) continue; /* too expensive */ |
if (ctpro == 1) goto DOGEN; |
P = (GEN)compositum(prev,P)[1]; |
storeprecdoubles(&Mbar, &Mbarst); /* restore */ |
if (degpol(P) == n) continue; |
if (! applybar(&M, &Mbar, Abargen,Bbargen)) err(bugparier,"pslqL2"); |
if (DEBUGLEVEL>2 && lgef(P)>lgef(prev)) |
if ( (p1 = checkend(&M, prec)) ) return gerepilecopy(av0, p1); |
fprintferr("chk_gen_init: subfield %Z\n",P); |
goto RESTART; |
} |
} |
skipfirst++; prev = P; |
|
} |
} |
} |
} |
|
else |
|
{ |
|
DOGEN: |
|
if ((p1 = one_step_gen(&M, tabga, prec))) |
|
return gerepilecopy(av, p1); |
|
} |
} |
} |
/* x_1,...,x_skipfirst generate a strict subfield [unless n=skipfirst=1] */ |
|
chk->skipfirst = skipfirst; |
|
if (DEBUGLEVEL>2) fprintferr("chk_gen_init: skipfirst = %ld\n",skipfirst); |
|
|
|
/* should be gexpo( [max_k C^n_k (bound/k) ^ k] ^ (1/2) ) */ |
|
prec2 = (1 + (((gexpo(bound)*n)/2) >> TWOPOTBITS_IN_LONG)); |
|
if (prec2 < 0) prec2 = 0; |
|
prec = 3 + prec2; |
|
prec2= nfgetprec(nf); |
|
if (DEBUGLEVEL) |
|
fprintferr("chk_gen_init: estimated prec = %ld (initially %ld)\n", |
|
prec, prec2); |
|
if (prec > prec2) return NULL; |
|
if (prec < prec2) data[1] = (long)gprec_w((GEN)data[1], prec); |
|
*ptprec = prec; return bound; |
|
} |
} |
|
|
static GEN |
/* x is a vector of elts of a p-adic field */ |
chk_gen_post(GEN data, GEN res) |
GEN |
|
plindep(GEN x) |
{ |
{ |
GEN basis = (GEN)data[2]; |
long i, j, prec = VERYBIGINT, lx = lg(x)-1, ly, v; |
GEN p1 = (GEN)res[2]; |
gpmem_t av = avma; |
long i, l = lg(p1); |
GEN p = NULL, pn,p1,m,a; |
for (i=1; i<l; i++) p1[i] = lmul(basis, (GEN)p1[i]); |
|
return res; |
|
} |
|
|
|
/* no garbage collecting, done in polredabs0 */ |
if (lx < 2) return cgetg(1,t_VEC); |
static GEN |
for (i=1; i<=lx; i++) |
findmindisc(GEN nf, GEN y, GEN a, GEN phimax, long flun) |
|
{ |
|
long i,k, c = lg(y); |
|
GEN v,dmin,z,beta,discs = cgetg(c,t_VEC); |
|
|
|
for (i=1; i<c; i++) discs[i] = labsi(discsr((GEN)y[i])); |
|
v = sindexsort(discs); |
|
k = v[1]; dmin = (GEN)discs[k]; z = (GEN)y[k]; beta = (GEN)a[k]; |
|
for (i=2; i<c; i++) |
|
{ |
{ |
k = v[i]; |
p1 = (GEN)x[i]; |
if (!egalii((GEN)discs[k],dmin)) break; |
if (typ(p1) != t_PADIC) continue; |
if (gpolcomp((GEN)y[k],z) < 0) { z = (GEN)y[k]; beta = (GEN)a[k]; } |
|
|
j = precp(p1); if (j < prec) prec = j; |
|
if (!p) p = (GEN)p1[2]; |
|
else if (!egalii(p, (GEN)p1[2])) |
|
err(talker,"inconsistent primes in plindep"); |
} |
} |
if (flun & nf_RAW) |
if (!p) err(talker,"not a p-adic vector in plindep"); |
{ |
v = ggval(x,p); pn = gpowgs(p,prec); |
y=cgetg(3,t_VEC); |
if (v != 0) x = gmul(x, gpowgs(p, -v)); |
y[1]=lcopy(z); |
x = lift_intern(gmul(x, gmodulcp(gun, pn))); |
y[2]=lcopy(beta); |
|
} |
|
else if (phimax) |
|
{ |
|
y=cgetg(3,t_VEC); |
|
y[1]=lcopy(z); |
|
beta=polymodrecip(gmodulcp(beta,(GEN)nf[1])); |
|
y[2]=(long)poleval(phimax,beta); |
|
} |
|
else y = gcopy(z); |
|
return y; |
|
} |
|
|
|
/* no garbage collecting, done in polredabs0 */ |
ly = 2*lx - 1; |
static GEN |
m = cgetg(ly+1,t_MAT); |
storeallpols(GEN nf, GEN z, GEN a, GEN phimax, long flun) |
for (j=1; j<=ly; j++) |
{ |
|
GEN p1,y,beta; |
|
|
|
if (flun & nf_RAW) |
|
{ |
{ |
long i, c = lg(z); |
p1 = cgetg(lx+1, t_COL); m[j] = (long)p1; |
y=cgetg(c,t_VEC); |
for (i=1; i<=lx; i++) p1[i] = zero; |
for (i=1; i<c; i++) |
|
{ |
|
p1=cgetg(3,t_VEC); y[i]=(long)p1; |
|
p1[1]=lcopy((GEN)z[i]); |
|
p1[2]=lcopy((GEN)a[i]); |
|
} |
|
} |
} |
else if (phimax) |
a = negi((GEN)x[1]); |
|
for (i=1; i<lx; i++) |
{ |
{ |
long i, c = lg(z); |
coeff(m,1+i,i) = (long)a; |
beta = new_chunk(c); |
coeff(m,1 ,i) = x[i+1]; |
for (i=1; i<c; i++) |
|
beta[i] = (long)polymodrecip(gmodulcp((GEN)a[i],(GEN)nf[1])); |
|
|
|
y=cgetg(c,t_VEC); |
|
for (i=1; i<c; i++) |
|
{ |
|
p1=cgetg(3,t_VEC); y[i]=(long)p1; |
|
p1[1]=lcopy((GEN)z[i]); |
|
p1[2]=(long)poleval(phimax,(GEN)beta[i]); |
|
} |
|
} |
} |
else y = gcopy(z); |
for (i=1; i<=lx; i++) coeff(m,i,lx-1+i) = (long)pn; |
return y; |
p1 = lllint(m); |
|
return gerepileupto(av, gmul(m, (GEN)p1[1])); |
} |
} |
|
|
GEN |
GEN |
polredabs0(GEN x, long flun, long prec) |
algdep0(GEN x, long n, long bit, long prec) |
{ |
{ |
long i,nv, av = avma; |
long tx=typ(x), i, k; |
GEN nf,v,y,a,phimax; |
gpmem_t av; |
GEN (*storepols)(GEN, GEN, GEN, GEN, long); |
GEN y,p1; |
FP_chk_fun *chk; |
|
|
|
chk = (FP_chk_fun*)new_chunk(sizeof(FP_chk_fun)); |
if (! is_scalar_t(tx)) err(typeer,"algdep0"); |
chk->f = &chk_gen; |
if (tx==t_POLMOD) { y=forcecopy((GEN)x[1]); setvarn(y,0); return y; } |
chk->f_init = &chk_gen_init; |
if (gcmp0(x)) return gzero; |
chk->f_post = &chk_gen_post; |
if (!n) return gun; |
|
|
if ((ulong)flun >= 16) err(flagerr,"polredabs"); |
av=avma; p1=cgetg(n+2,t_COL); |
nf = initalgall0(x,nf_SMALL|nf_REGULAR,prec); |
p1[1] = un; |
if (lg(nf) == 3) |
p1[2] = (long)x; /* n >= 1 */ |
{ |
for (i=3; i<=n+1; i++) p1[i]=lmul((GEN)p1[i-1],x); |
phimax = lift_to_pol((GEN)nf[2]); |
if (typ(x) == t_PADIC) |
nf = (GEN)nf[1]; |
p1 = plindep(p1); |
} |
|
else |
else |
phimax = (flun & nf_ORIG)? polx[0]: (GEN)NULL; |
|
prec = nfgetprec(nf); |
|
x = (GEN)nf[1]; |
|
|
|
if (lgef(x) == 4) |
|
{ |
{ |
y = _vec(polx[varn(x)]); |
switch(bit) |
a = _vec(gsub((GEN)y[1], (GEN)x[2])); |
|
} |
|
else |
|
{ |
|
for (i=1; ; i++) |
|
{ |
{ |
v = fincke_pohst(nf,NULL,5000,3,prec, chk); |
case 0: p1 = pslq(p1,prec); break; |
if (v) break; |
case -1: p1 = lindep(p1,prec); break; |
if (i==MAXITERPOL) err(accurer,"polredabs0"); |
case -2: p1 = deplin(p1); break; |
prec = (prec<<1)-2; nf = nfnewprec(nf,prec); |
default: p1 = lindep2(p1,bit); |
if (DEBUGLEVEL) err(warnprec,"polredabs0",prec); |
|
} |
} |
a = (GEN)v[2]; |
|
y = (GEN)v[1]; |
|
} |
} |
nv = lg(a); |
if ((!bit) && (typ(p1) == t_REAL)) |
for (i=1; i<nv; i++) |
|
if (canon_pol((GEN)y[i]) < 0 && phimax) |
|
a[i] = (long) gneg_i((GEN)a[i]); |
|
nv = remove_duplicates(y,a); |
|
|
|
if (DEBUGLEVEL) |
|
{ fprintferr("%ld minimal vectors found.\n",nv-1); flusherr(); } |
|
if (nv >= 10000) flun &= (~nf_ALL); /* should not happen */ |
|
storepols = (flun & nf_ALL)? storeallpols: findmindisc; |
|
|
|
if (DEBUGLEVEL) fprintferr("\n"); |
|
if (nv==1) |
|
{ |
{ |
y = _vec(x); |
y = gcopy(p1); return gerepileupto(av,y); |
a = _vec(polx[varn(x)]); |
|
} |
} |
if (varn(y[1]) != varn(x)) |
if (lg(p1) < 2) |
{ |
err(talker,"higher degree than expected in algdep"); |
long vx = varn(x); |
y=cgetg(n+3,t_POL); |
for (i=1; i<nv; i++) setvarn(y[i], vx); |
y[1] = evalsigne(1) | evalvarn(0); |
} |
k=1; while (gcmp0((GEN)p1[k])) k++; |
return gerepileupto(av, storepols(nf,y,a,phimax,flun)); |
for (i=0; i<=n+1-k; i++) y[i+2] = p1[k+i]; |
|
(void)normalizepol_i(y, n+4-k); |
|
y = (gsigne(leading_term(y)) > 0)? gcopy(y): gneg(y); |
|
return gerepileupto(av,y); |
} |
} |
|
|
GEN |
GEN |
polredabsall(GEN x, long flun, long prec) |
algdep2(GEN x, long n, long bit) |
{ |
{ |
return polredabs0(x, flun | nf_ALL, prec); |
return algdep0(x,n,bit,0); |
} |
} |
|
|
GEN |
GEN |
polredabs(GEN x, long prec) |
algdep(GEN x, long n, long prec) |
{ |
{ |
return polredabs0(x,nf_REGULAR,prec); |
return algdep0(x,n,0,prec); |
} |
} |
|
|
GEN |
/********************************************************************/ |
polredabs2(GEN x, long prec) |
/** **/ |
{ |
/** INTEGRAL KERNEL (LLL REDUCED) **/ |
return polredabs0(x,nf_ORIG,prec); |
/** **/ |
} |
/********************************************************************/ |
|
|
GEN |
GEN |
polredabsnored(GEN x, long prec) |
matkerint0(GEN x, long flag) |
{ |
{ |
return polredabs0(x,nf_NORED,prec); |
switch(flag) |
|
{ |
|
case 0: return kerint(x); |
|
case 1: return kerint1(x); |
|
default: err(flagerr,"matkerint"); |
|
} |
|
return NULL; /* not reached */ |
} |
} |
|
|
GEN |
GEN |
polred(GEN x, long prec) |
kerint1(GEN x) |
{ |
{ |
return allpolred(x,(GEN*)0,0,prec); |
gpmem_t av = avma; |
|
return gerepilecopy(av, lllint_ip(matrixqz3(ker(x)), 100)); |
} |
} |
|
|
GEN |
|
polredfirstpol(GEN x, long prec, int (*check)(GEN,GEN), GEN arg) |
|
{ |
|
return allpolred0(x,(GEN*)0,0,prec,check,arg); |
|
} |
|
|
|
GEN |
GEN |
smallpolred(GEN x, long prec) |
kerint(GEN x) |
{ |
{ |
return allpolred(x,(GEN*)0,1,prec); |
gpmem_t av = avma; |
|
GEN fl, junk, h = lllint_i(x, 0, 0, &junk, &fl, NULL); |
|
if (h) h = lll_finish(h,fl, lll_KER); |
|
else h = lll_trivial(x, lll_KER); |
|
if (lg(h)==1) { avma = av; return cgetg(1, t_MAT); } |
|
return gerepilecopy(av, lllint_ip(h, 100)); |
} |
} |
|
|
GEN |
|
factoredpolred(GEN x, GEN p, long prec) |
|
{ |
|
return allpolred(x,(GEN*)0,(long)p,prec); |
|
} |
|
|
|
GEN |
|
polred2(GEN x, long prec) |
|
{ |
|
GEN y=cgetg(3,t_MAT); |
|
|
|
y[2]= (long) allpolred(x,(GEN*)(y+1),0,prec); |
|
return y; |
|
} |
|
|
|
GEN |
|
smallpolred2(GEN x, long prec) |
|
{ |
|
GEN y=cgetg(3,t_MAT); |
|
|
|
y[2]= (long) allpolred(x,(GEN*)(y+1),1,prec); |
|
return y; |
|
} |
|
|
|
GEN |
|
factoredpolred2(GEN x, GEN p, long prec) |
|
{ |
|
GEN y=cgetg(3,t_MAT); |
|
|
|
y[2]= (long) allpolred(x,(GEN*)(y+1),(long)p,prec); |
|
return y; |
|
} |
|
|
|
/* relative polredabs. Returns |
|
* flag = 0: relative polynomial |
|
* flag = 1: relative polynomial + element |
|
* flag = 2: absolute polynomial */ |
|
GEN |
|
rnfpolredabs(GEN nf, GEN relpol, long flag, long prec) |
|
{ |
|
GEN p1,bpol,rnf,elt,pol; |
|
long va, av = avma; |
|
|
|
if (typ(relpol)!=t_POL) err(typeer,"rnfpolredabs"); |
|
nf=checknf(nf); va = varn(relpol); |
|
if (DEBUGLEVEL>1) timer2(); |
|
p1 = makebasis(nf, unifpol(nf,relpol,1)); |
|
rnf = (GEN)p1[3]; |
|
if (DEBUGLEVEL>1) |
|
{ |
|
msgtimer("absolute basis"); |
|
fprintferr("original absolute generator: %Z\n",p1[1]); |
|
} |
|
p1 = polredabs0(p1, nf_RAW | nf_NORED, prec); |
|
bpol = (GEN)p1[1]; |
|
if (DEBUGLEVEL>1) fprintferr("reduced absolute generator: %Z\n",bpol); |
|
if (flag==2) return gerepileupto(av,bpol); |
|
|
|
elt = rnfelementabstorel(rnf,(GEN)p1[2]); |
|
p1=cgetg(3,t_VEC); |
|
pol = rnfcharpoly(nf,relpol,elt,va); |
|
if (!flag) p1 = pol; |
|
else |
|
{ |
|
p1[1]=(long)pol; |
|
p1[2]=(long)polymodrecip(elt); |
|
} |
|
return gerepileupto(av,p1); |
|
} |
|
|
|
/********************************************************************/ |
/********************************************************************/ |
/** **/ |
/** **/ |
/** MINIM **/ |
/** MINIM **/ |
/** **/ |
/** **/ |
/********************************************************************/ |
/********************************************************************/ |
int addcolumntomatrix(GEN V,GEN INVP,GEN L); |
|
|
|
/* x is a non-empty matrix, y is a compatible VECSMALL (dimension > 0). */ |
/* x is a non-empty matrix, y is a compatible VECSMALL (dimension > 0). */ |
GEN |
GEN |
gmul_mat_smallvec(GEN x, GEN y) |
gmul_mat_smallvec(GEN x, GEN y) |
{ |
{ |
long c = lg(x), l = lg(x[1]), av,i,j; |
long c = lg(x), l = lg(x[1]), i, j; |
|
gpmem_t av; |
GEN z=cgetg(l,t_COL), s; |
GEN z=cgetg(l,t_COL), s; |
|
|
for (i=1; i<l; i++) |
for (i=1; i<l; i++) |
Line 2520 gmul_mat_smallvec(GEN x, GEN y) |
|
Line 2599 gmul_mat_smallvec(GEN x, GEN y) |
|
GEN |
GEN |
gmul_mati_smallvec(GEN x, GEN y) |
gmul_mati_smallvec(GEN x, GEN y) |
{ |
{ |
long c = lg(x), l = lg(x[1]), av,i,j; |
long c = lg(x), l = lg(x[1]), i, j; |
|
gpmem_t av; |
GEN z=cgetg(l,t_COL), s; |
GEN z=cgetg(l,t_COL), s; |
|
|
for (i=1; i<l; i++) |
for (i=1; i<l; i++) |
Line 2528 gmul_mati_smallvec(GEN x, GEN y) |
|
Line 2608 gmul_mati_smallvec(GEN x, GEN y) |
|
av = avma; s = mulis(gcoeff(x,i,1),y[1]); |
av = avma; s = mulis(gcoeff(x,i,1),y[1]); |
for (j=2; j<c; j++) |
for (j=2; j<c; j++) |
if (y[j]) s = addii(s, mulis(gcoeff(x,i,j),y[j])); |
if (y[j]) s = addii(s, mulis(gcoeff(x,i,j),y[j])); |
z[i]=(long)gerepileuptoint(av,s); |
z[i] = lpileuptoint(av,s); |
} |
} |
return z; |
return z; |
} |
} |
|
|
minim_alloc(long n, double ***q, GEN *x, double **y, double **z, double **v) |
minim_alloc(long n, double ***q, GEN *x, double **y, double **z, double **v) |
{ |
{ |
long i, s; |
long i, s; |
double **Q; |
|
|
|
*x = cgetg(n, t_VECSMALL); |
*x = cgetg(n, t_VECSMALL); |
*q = (double**) new_chunk(n); |
*q = (double**) new_chunk(n); |
|
s = n * sizeof(double); |
/* correct alignment for the following */ |
init_dalloc(); |
s = avma % sizeof(double); avma -= s; |
*y = dalloc(s); |
if (avma<bot) err(errpile); |
*z = dalloc(s); |
|
*v = dalloc(s); |
s = (n * sizeof(double))/sizeof(long); |
for (i=1; i<n; i++) (*q)[i] = dalloc(s); |
*y = (double*) new_chunk(s); |
|
*z = (double*) new_chunk(s); |
|
*v = (double*) new_chunk(s); Q = *q; |
|
for (i=1; i<n; i++) Q[i] = (double*) new_chunk(s); |
|
} |
} |
|
|
/* Minimal vectors for the integral definite quadratic form: a. |
/* Minimal vectors for the integral definite quadratic form: a. |
|
|
minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag) |
minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag) |
{ |
{ |
GEN x,res,p1,u,r,liste,gnorme,invp,V, *gptr[2]; |
GEN x,res,p1,u,r,liste,gnorme,invp,V, *gptr[2]; |
long n = lg(a), av0 = avma, av1,av,tetpil,lim, i,j,k,s,maxrank; |
long n = lg(a), i, j, k, s, maxrank; |
|
gpmem_t av0 = avma, av1, av, tetpil, lim; |
double p,maxnorm,BOUND,*v,*y,*z,**q, eps = 0.000001; |
double p,maxnorm,BOUND,*v,*y,*z,**q, eps = 0.000001; |
|
|
maxrank = 0; res = V = invp = NULL; /* gcc -Wall */ |
maxrank = 0; res = V = invp = NULL; /* gcc -Wall */ |
Line 2683 minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag) |
|
Line 2759 minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag) |
|
} |
} |
else |
else |
{ |
{ |
long av2 = avma; |
gpmem_t av2 = avma; |
gnorme = ground(dbltor(p)); |
gnorme = ground(dbltor(p)); |
if (gcmp(gnorme,BORNE) >= 0) avma = av2; |
if (gcmp(gnorme,BORNE) >= 0) avma = av2; |
else |
else |
Line 2707 minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag) |
|
Line 2783 minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag) |
|
|
|
case min_PERF: |
case min_PERF: |
{ |
{ |
long av2=avma, I=1; |
long I=1; |
|
gpmem_t av2=avma; |
|
|
for (i=1; i<=n; i++) |
for (i=1; i<=n; i++) |
for (j=i; j<=n; j++,I++) V[I] = x[i]*x[j]; |
for (j=i; j<=n; j++,I++) V[I] = x[i]*x[j]; |
|
|
|
|
/* general program for positive definit quadratic forms (real coeffs). |
/* general program for positive definit quadratic forms (real coeffs). |
* One needs BORNE != 0; LLL reduction done in fincke_pohst(). |
* One needs BORNE != 0; LLL reduction done in fincke_pohst(). |
* If (flag & 2) stop as soon as stockmax is reached. |
* If (noer) return NULL on precision problems (no error). |
* If (flag & 1) return NULL on precision problems (no error). |
|
* If (check != NULL consider only vectors passing the check [ assumes |
* If (check != NULL consider only vectors passing the check [ assumes |
* stockmax > 0 and we only want the smallest possible vectors ] */ |
* stockmax > 0 and we only want the smallest possible vectors ] */ |
static GEN |
static GEN |
smallvectors(GEN a, GEN BORNE, long stockmax, long flag, FP_chk_fun *CHECK) |
smallvectors(GEN a, GEN BORNE, long stockmax, long noer, FP_chk_fun *CHECK) |
{ |
{ |
long av,av1,lim,N,n,i,j,k,s,epsbit,prec, checkcnt = 1; |
long N, n, i, j, k, s, epsbit, prec, checkcnt = 1; |
|
gpmem_t av, av1, lim; |
GEN u,S,x,y,z,v,q,norme1,normax1,borne1,borne2,eps,p1,alpha,norms,dummy; |
GEN u,S,x,y,z,v,q,norme1,normax1,borne1,borne2,eps,p1,alpha,norms,dummy; |
GEN (*check)(GEN,GEN) = CHECK? CHECK->f: NULL; |
GEN (*check)(void *,GEN) = CHECK? CHECK->f: NULL; |
GEN data = CHECK? CHECK->data: NULL; |
void *data = CHECK? CHECK->data: NULL; |
int skipfirst = CHECK? CHECK->skipfirst: 0; |
int skipfirst = CHECK? CHECK->skipfirst: 0; |
|
|
if (DEBUGLEVEL) |
if (DEBUGLEVEL) |
fprintferr("smallvectors looking for norm <= %Z\n",gprec_w(BORNE,3)); |
fprintferr("smallvectors looking for norm <= %Z\n",gprec_w(BORNE,3)); |
|
|
q = sqred1intern(a,flag&1); |
q = sqred1intern(a, noer); |
if (q == NULL) return NULL; |
if (q == NULL) return NULL; |
if (DEBUGLEVEL>5) fprintferr("q = %Z",q); |
if (DEBUGLEVEL>5) fprintferr("q = %Z",q); |
prec = gprecision(q); |
prec = gprecision(q); |
epsbit = bit_accuracy(prec) >> 1; |
epsbit = bit_accuracy(prec) >> 1; |
eps = realun(prec); setexpo(eps,-epsbit); |
eps = real2n(-epsbit, prec); |
alpha = dbltor(0.95); |
alpha = dbltor(0.95); |
normax1 = gzero; |
normax1 = gzero; |
borne1= gadd(BORNE,eps); |
borne1= gadd(BORNE,eps); |
Line 2870 smallvectors(GEN a, GEN BORNE, long stockmax, long fla |
|
Line 2947 smallvectors(GEN a, GEN BORNE, long stockmax, long fla |
|
|
|
if (low_stack(lim, stack_lim(av,1))) |
if (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
GEN *gptr[7]; |
|
int cnt = 4; |
int cnt = 4; |
if(DEBUGMEM>1) err(warnmem,"smallvectors"); |
if(DEBUGMEM>1) err(warnmem,"smallvectors"); |
gptr[0]=&x; gptr[1]=&y; gptr[2]=&z; gptr[3]=&normax1; |
|
if (stockmax) |
if (stockmax) |
{ /* initialize to dummy values */ |
{ /* initialize to dummy values */ |
GEN T = S; |
GEN T = S; |
Line 2882 smallvectors(GEN a, GEN BORNE, long stockmax, long fla |
|
Line 2957 smallvectors(GEN a, GEN BORNE, long stockmax, long fla |
|
} |
} |
if (check) |
if (check) |
{ |
{ |
cnt+=3; gptr[4]=&borne1; gptr[5]=&borne2; gptr[6]=&norms; |
cnt+=3; |
for (i=s+1; i<=stockmax; i++) norms[i]=(long)dummy; |
for (i=s+1; i<=stockmax; i++) norms[i]=(long)dummy; |
} |
} |
gerepilemany(av,gptr,cnt); |
gerepileall(av,cnt,&x,&y,&z,&normax1,&borne1,&borne2,&norms); |
} |
} |
if (DEBUGLEVEL>3) |
if (DEBUGLEVEL>3) |
{ |
{ |
Line 2922 smallvectors(GEN a, GEN BORNE, long stockmax, long fla |
|
Line 2997 smallvectors(GEN a, GEN BORNE, long stockmax, long fla |
|
{ |
{ |
if (check) norms[s] = (long)norme1; |
if (check) norms[s] = (long)norme1; |
S[s] = (long)dummycopy(x); |
S[s] = (long)dummycopy(x); |
if (s == stockmax && (flag&2) && check) |
if (s == stockmax) |
{ |
{ |
long av1 = avma; |
gpmem_t av2 = avma; |
GEN per = sindexsort(norms); |
GEN per; |
|
|
|
if (!check) goto END; |
|
per = sindexsort(norms); |
if (DEBUGLEVEL) fprintferr("sorting...\n"); |
if (DEBUGLEVEL) fprintferr("sorting...\n"); |
for (i=1; i<=s; i++) |
for (i=1; i<=s; i++) |
{ |
{ |
long k = per[i]; |
long k = per[i]; |
if (check(data,(GEN)S[k])) |
if (check(data,(GEN)S[k])) |
{ |
{ |
S[1] = S[k]; avma = av1; |
S[1] = S[k]; avma = av2; |
borne1 = mpadd(norme1,eps); |
borne1 = mpadd(norme1,eps); |
borne2 = mpmul(borne1,alpha); |
borne2 = mpmul(borne1,alpha); |
s = 1; checkcnt = 0; break; |
s = 1; checkcnt = 0; break; |
|
|
u[3] = (long)S; return u; |
u[3] = (long)S; return u; |
} |
} |
|
|
/* solve x~.a.x <= bound, a > 0. If check is non-NULL keep x only if check(x). |
/* solve q(x) = x~.a.x <= bound, a > 0. |
* flag & 1, return NULL in case of precision problems (sqred1 or lllgram) |
* If check is non-NULL keep x only if check(x). |
* raise an error otherwise. |
* If noer, return NULL in case of precision problems, raise an error otherwise. |
* flag & 2, return as soon as stockmax vectors found. |
* If a is a vector, assume a[1] is the LLL-reduced Cholesky form of q */ |
* If a is a number field, use its T2 matrix |
static GEN |
*/ |
_fincke_pohst(GEN a,GEN B0,long stockmax,long noer, long PREC, FP_chk_fun *CHECK) |
GEN |
|
fincke_pohst(GEN a,GEN B0,long stockmax,long flag, long PREC, FP_chk_fun *CHECK) |
|
{ |
{ |
VOLATILE long pr,av=avma,i,j,n; |
gpmem_t av = avma; |
VOLATILE GEN B,nf,r,rinvtrans,u,v,res,z,vnorm,sperm,perm,uperm,gram; |
long i,j,l, round = 0; |
VOLATILE GEN bound = B0; |
GEN B,r,rinvtrans,u,v,res,z,vnorm,sperm,perm,uperm,gram, bound = B0; |
void *catcherr = NULL; |
|
long prec = PREC; |
|
|
|
if (DEBUGLEVEL>2) { fprintferr("entering fincke_pohst\n"); flusherr(); } |
if (DEBUGLEVEL>2) fprintferr("entering fincke_pohst\n"); |
if (typ(a) == t_VEC) { nf = checknf(a); a = gmael(nf,5,3); } else nf = NULL; |
if (typ(a) == t_VEC) |
pr = gprecision(a); |
|
if (pr) prec = pr; else a = gmul(a,realun(prec)); |
|
if (DEBUGLEVEL>2) fprintferr("first LLL: prec = %ld\n", prec); |
|
if (nf && !signe(gmael(nf,2,2))) /* totally real */ |
|
{ |
{ |
GEN T = nf_get_T((GEN)nf[1], (GEN)nf[7]); |
r = (GEN)a[1]; |
u = lllgramint(T); |
u = NULL; |
prec += 2 * ((gexpo(u) + gexpo((GEN)nf[1])) >> TWOPOTBITS_IN_LONG); |
|
nf = nfnewprec(nf, prec); |
|
r = qf_base_change(T,u,1); |
|
i = PREC + (gexpo(r) >> TWOPOTBITS_IN_LONG); |
|
if (i < prec) prec = i; |
|
r = gmul(r, realun(prec)); |
|
} |
} |
else |
else |
{ |
{ |
u = lllgramintern(a,4,flag&1, (prec<<1)-2); |
long prec = PREC; |
if (!u) goto PRECPB; |
l = lg(a); |
|
if (l == 1) |
|
{ |
|
if (CHECK) err(talker, "dimension 0 in fincke_pohst"); |
|
z = cgetg(4,t_VEC); |
|
z[1] = z[2] = zero; |
|
z[3] = lgetg(1,t_MAT); return z; |
|
} |
|
i = gprecision(a); |
|
if (i) prec = i; else { a = gmul(a,realun(prec)); round = 1; } |
|
if (DEBUGLEVEL>2) fprintferr("first LLL: prec = %ld\n", prec); |
|
u = lllgramintern(a,4,noer, (prec<<1)-2); |
|
if (!u) return NULL; |
r = qf_base_change(a,u,1); |
r = qf_base_change(a,u,1); |
|
r = sqred1intern(r,noer); |
|
if (!r) return NULL; |
|
for (i=1; i<l; i++) |
|
{ |
|
GEN s = gsqrt(gcoeff(r,i,i), prec); |
|
coeff(r,i,i) = (long)s; |
|
for (j=i+1; j<l; j++) coeff(r,i,j) = lmul(s, gcoeff(r,i,j)); |
|
} |
} |
} |
r = sqred1intern(r,flag&1); |
/* now r~ * r = a in LLL basis */ |
if (!r) goto PRECPB; |
|
|
|
n = lg(a); |
|
if (n == 1) |
|
{ |
|
if (CHECK) err(talker, "dimension 0 in fincke_pohst"); |
|
avma = av; z=cgetg(4,t_VEC); |
|
z[1]=z[2]=zero; |
|
z[3]=lgetg(1,t_MAT); return z; |
|
} |
|
for (i=1; i<n; i++) |
|
{ |
|
GEN p1 = gsqrt(gcoeff(r,i,i), prec); |
|
coeff(r,i,i)=(long)p1; |
|
for (j=i+1; j<n; j++) |
|
coeff(r,i,j) = lmul(p1, gcoeff(r,i,j)); |
|
} |
|
/* now r~ * r = a in approximate LLL basis */ |
|
rinvtrans = gtrans(invmat(r)); |
rinvtrans = gtrans(invmat(r)); |
|
if (DEBUGLEVEL>2) |
|
fprintferr("final LLL: prec = %ld\n", gprecision(rinvtrans)); |
|
v = lllintern(rinvtrans, 100, noer, 0); |
|
if (!v) return NULL; |
|
rinvtrans = gmul(rinvtrans, v); |
|
|
v = NULL; |
v = ZM_inv(gtrans_i(v),gun); |
for (i=1; i<6; i++) /* try to get close to a genuine LLL basis */ |
r = gmul(r,v); |
{ |
u = u? gmul(u,v): v; |
GEN p1; |
|
if (DEBUGLEVEL>2) |
|
fprintferr("final LLLs: prec = %ld, precision(rinvtrans) = %ld\n", |
|
prec,gprecision(rinvtrans)); |
|
p1 = lllintern(rinvtrans,flag&1, (gprecision(rinvtrans)<<1)-2); |
|
if (!p1) goto PRECPB; |
|
if (ishnfall(p1)) break; /* upper triangular */ |
|
if (v) v = gmul(v,p1); else v = p1; |
|
rinvtrans = gmul(rinvtrans,p1); |
|
} |
|
if (i == 6) goto PRECPB; /* diverges... */ |
|
|
|
if (v) |
l = lg(r); |
{ |
vnorm = cgetg(l,t_VEC); |
GEN u2 = ZM_inv(gtrans(v),gun); |
for (j=1; j<l; j++) vnorm[j] = lnorml2((GEN)rinvtrans[j]); |
r = gmul(r,u2); /* r = LLL basis now */ |
sperm = cgetg(l,t_MAT); |
u = gmul(u,u2); |
uperm = cgetg(l,t_MAT); perm = sindexsort(vnorm); |
} |
for (i=1; i<l; i++) { uperm[l-i] = u[perm[i]]; sperm[l-i] = r[perm[i]]; } |
|
|
vnorm = cgetg(n,t_VEC); |
|
for (j=1; j<n; j++) vnorm[j] = lnorml2((GEN)rinvtrans[j]); |
|
sperm = cgetg(n,t_MAT); |
|
uperm = cgetg(n,t_MAT); perm = sindexsort(vnorm); |
|
for (i=1; i<n; i++) { uperm[n-i] = u[perm[i]]; sperm[n-i] = r[perm[i]]; } |
|
|
|
gram = gram_matrix(sperm); |
gram = gram_matrix(sperm); |
B = gcoeff(gram,n-1,n-1); |
B = gcoeff(gram,l-1,l-1); |
if (gexpo(B) >= bit_accuracy(lg(B)-2)) goto PRECPB; |
if (gexpo(B) >= bit_accuracy(lg(B)-2)) return NULL; |
|
|
{ /* catch precision problems (truncation error) */ |
res = NULL; |
jmp_buf env; |
CATCH(precer) { } |
if (setjmp(env)) goto PRECPB; |
TRY { |
catcherr = err_catch(precer, env, NULL); |
if (CHECK && CHECK->f_init) |
} |
|
if (CHECK && CHECK->f_init) |
|
{ |
|
bound = CHECK->f_init(CHECK,nf,gram,uperm, &prec); |
|
if (!bound) goto PRECPB; |
|
} |
|
if (!bound) |
|
{ /* set bound */ |
|
GEN x = cgetg(n,t_COL); |
|
if (nf) bound = get_Bnf(nf); |
|
for (i=2; i<n; i++) x[i]=zero; |
|
for (i=1; i<n; i++) |
|
{ |
{ |
x[i] = un; B = gcoeff(gram,i,i); |
bound = CHECK->f_init(CHECK,gram,uperm); |
if (!bound || mpcmp(B,bound) < 0) bound = B; |
if (!bound) err(precer,"fincke_pohst"); |
x[i] = zero; |
|
} |
} |
} |
if (!bound) bound = gcoeff(gram,1,1); |
|
|
if (DEBUGLEVEL>2) {fprintferr("entering smallvectors\n"); flusherr();} |
if (DEBUGLEVEL>2) fprintferr("entering smallvectors\n"); |
for (i=1; i<n; i++) |
for (i=1; i<l; i++) |
|
{ |
|
res = smallvectors(gram, bound, stockmax,noer,CHECK); |
|
if (!res) err(precer,"fincke_pohst"); |
|
if (!CHECK || bound || lg(res[2]) > 1) break; |
|
if (DEBUGLEVEL>2) fprintferr(" i = %ld failed\n",i); |
|
} |
|
} ENDCATCH; |
|
if (DEBUGLEVEL>2) fprintferr("leaving fincke_pohst\n"); |
|
if (CHECK || !res) return res; |
|
|
|
z = cgetg(4,t_VEC); |
|
z[1] = lcopy((GEN)res[1]); |
|
z[2] = round? lround((GEN)res[2]): lcopy((GEN)res[2]); |
|
z[3] = lmul(uperm, (GEN)res[3]); return gerepileupto(av,z); |
|
} |
|
|
|
GEN |
|
fincke_pohst(GEN a,GEN B0,long stockmax,long flag, long PREC, FP_chk_fun *CHECK) |
|
{ |
|
gpmem_t av = avma; |
|
GEN z = _fincke_pohst(a,B0,stockmax,flag & 1,PREC,CHECK); |
|
if (!z) |
{ |
{ |
res = smallvectors(gram, bound? bound: gcoeff(gram,i,i), |
if (!(flag & 1)) err(precer,"fincke_pohst"); |
stockmax,flag,CHECK); |
avma = av; |
if (!res) goto PRECPB; |
|
if (!CHECK || bound || lg(res[2]) > 1) break; |
|
if (DEBUGLEVEL>2) fprintferr(" i = %ld failed\n",i); |
|
} |
} |
err_leave(&catcherr); catcherr = NULL; |
return z; |
if (CHECK) |
|
{ |
|
if (CHECK->f_post) res = CHECK->f_post(CHECK->data, res); |
|
return res; |
|
} |
|
|
|
if (DEBUGLEVEL>2) {fprintferr("leaving fincke_pohst\n"); flusherr();} |
|
z=cgetg(4,t_VEC); |
|
z[1]=lcopy((GEN)res[1]); |
|
z[2]=pr? lcopy((GEN)res[2]) : lround((GEN)res[2]); |
|
z[3]=lmul(uperm, (GEN)res[3]); return gerepileupto(av,z); |
|
PRECPB: |
|
if (catcherr) err_leave(&catcherr); |
|
if (!(flag & 1)) |
|
err(talker,"not a positive definite form in fincke_pohst"); |
|
avma = av; return NULL; |
|
} |
} |