File: [local] / OpenXM_contrib / pari / src / modules / Attic / thue.c (download)
Revision 1.1.1.1 (vendor branch), Sun Jan 9 17:35:33 2000 UTC (24 years, 5 months ago) by maekawa
Branch: PARI_GP
CVS Tags: maekawa-ipv6, VERSION_2_0_17_BETA, RELEASE_20000124, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, RELEASE_1_1_3, RELEASE_1_1_2
Changes since 1.1: +0 -0
lines
Import PARI/GP 2.0.17 beta.
|
/* thue.c -- Thue equation solver. In all the forthcoming remarks, "paper"
* designs the paper "Thue Equations of High Degree", by Yu. Bilu and
* G. Hanrot, J. Number Theory (1996). Note that numbering of the constants
* is that of Hanrot's thesis rather than that of the paper.
* The last part of the program (bnfisintnorm) was written by K. Belabas.
*/
/* $Id: thue.c,v 1.1.1.1 1999/09/16 13:48:20 karim Exp $ */
#include "pari.h"
static int curne,r,s,t,deg,Prec,ConstPrec,numroot;
static GEN c1,c2,c3,c4,c5,c6,c7,c8,c10,c11,c12,c13,c14,c15,B0,x1,x2,x3,x0;
static GEN halpha,eps3,gdeg,A,MatNE,roo,MatFU,Delta,Lambda,delta,lambda;
static GEN Vect2,SOL,uftot;
GEN bnfisintnorm(GEN, GEN);
/* Check whether tnf is a valid structure */
static int
checktnf(GEN tnf)
{
if (typ(tnf)!=t_VEC || (lg(tnf)!=8 && lg(tnf)!=3)) return 0;
if (typ(tnf[1])!=t_POL) return 0;
if (lg(tnf) != 8) return 1; /* s=0 */
deg=lgef(tnf[1])-3;
if (deg<=2) err(talker,"invalid polynomial in thue (need deg>2)");
s=sturm((GEN)tnf[1]); t=(deg-s)>>1; r=s+t-1;
(void)checkbnf((GEN)tnf[2]);
if (typ(tnf[3]) != t_COL || lg(tnf[3]) != deg+1) return 0;
if (typ(tnf[4]) != t_COL || lg(tnf[4]) != r+1) return 0;
if (typ(tnf[5]) != t_MAT || lg(tnf[5]) != r+1
|| lg(gmael(tnf,5,1)) != deg+1) return 0;
if (typ(tnf[6]) != t_MAT || lg(tnf[6])!=r+1
|| lg(gmael(tnf,6,1)) != r+1) return 0;
if (typ(tnf[7]) != t_VEC || lg(tnf[7]) != 7) return 0;
return 1;
}
static GEN distoZ(GEN z)
{
GEN p1=gfrac(z);
return gmin(p1, gsub(gun,p1));
}
/* Check whether a solution has already been found */
static int
_thue_new(GEN zz)
{
int i;
for (i=1; i<lg(SOL); i++)
if (gegal(zz,(GEN)SOL[i])) return 0;
return 1;
}
/* Compensates rounding errors for computation/display of the constants */
static GEN
myround(GEN Cst, GEN upd)
{
return gmul(Cst, gadd(gun, gmul(upd,gpuigs(stoi(10),-10))));
}
/* Returns the index of the largest element in a vector */
static int
Vecmax(GEN Vec, int size)
{
int k, tno = 1;
GEN tmax = (GEN)Vec[1];
for (k=2; k<=size; k++)
if (gcmp((GEN)Vec[k],tmax)==1) { tmax=(GEN)Vec[k]; tno=k; }
return tno;
}
/* Performs basic computations concerning the equation: buchinitfu, c1, c2 */
static void
inithue(GEN poly, long flag)
{
GEN roo2, tmp, gpmin, de;
int k,j;
x0=gzero; x1=gzero; deg=lgef(poly)-3; gdeg=stoi(deg);
if (!uftot)
{
uftot=bnfinit0(poly,1,NULL,Prec);
if (flag) certifybuchall(uftot);
s=itos(gmael3(uftot,7,2,1));
t=itos(gmael3(uftot,7,2,2));
}
/* Switch the roots to get the usual order */
roo=roots(poly, Prec); roo2=cgetg(deg+1,t_COL);
for (k=1; k<=s; k++) roo2[k]=roo[k];
for (k=1; k<=t; k++)
{
roo2[k+s]=roo[2*k-1+s];
roo2[k+s+t]=lconj((GEN)roo2[k+s]);
}
roo=roo2;
r=s+t-1; halpha=gun;
for (k=1; k<=deg; k++)
halpha = gmul(halpha,gmax(gun,gabs((GEN)roo[k],Prec)));
halpha=gdiv(glog(halpha,Prec),gdeg);
de=derivpol(poly); c1=gabs(poleval(de,(GEN)roo[1]),Prec);
for (k=2; k<=s; k++)
{
tmp=gabs(poleval(de,(GEN)roo[k]),Prec);
if (gcmp(tmp,c1)==-1) c1=tmp;
}
c1=gdiv(gpui(gdeux,gsub(gdeg,gun),Prec),c1); c1=myround(c1,gun);
c2=gabs(gsub((GEN)roo[1],(GEN)roo[2]),Prec);
for (k=1; k<=deg; k++)
for (j=k+1; j<=deg; j++)
{
tmp=gabs(gsub((GEN)roo[j],(GEN)roo[k]),Prec);
if (gcmp(c2,tmp)==1) c2=tmp;
}
c2=myround(c2,stoi(-1));
if (t==0) x0=gun;
else
{
gpmin=gabs(poleval(de,(GEN)roo[s+1]),Prec);
for (k=2; k<=t; k++)
{
tmp=gabs(poleval(de,(GEN)roo[s+k]),Prec);
if (gcmp(tmp,gpmin)==-1) gpmin=tmp;
}
/* Compute x0. See paper, Prop. 2.2.1 */
x0=gpui(gdiv(gpui(gdeux,gsub(gdeg,gun),Prec),
gmul(gpmin,
gabs((GEN)gimag(roo)[Vecmax(gabs(gimag(roo),Prec),deg)],Prec))),
ginv(gdeg),Prec);
}
if (DEBUGLEVEL >=2) fprintferr("c1 = %Z\nc2 = %Z\n", c1, c2);
}
/* Computation of M, its inverse A and check of the precision (see paper) */
static void
T_A_Matrices()
{
GEN mask, eps1, eps2, nia, m1, IntM;
int i,j;
m1=glog(gabs(MatFU,Prec),Prec); mask=gsub(gpui(gdeux,stoi(r),Prec),gun);
m1=matextract(m1,mask,mask);
A=invmat(m1); IntM=gsub(gmul(A,m1),idmat(r));
eps2=gzero;
for (i=1; i<=r; i++)
for (j=1; j<=r; j++)
if (gcmp(eps2,gabs(gcoeff(IntM,i,j),20)) == -1)
eps2=gabs(gcoeff(IntM,i,j),20);
eps1=gpui(gdeux,stoi((Prec-2)*32),Prec); eps2=gadd(eps2,ginv(eps1));
nia=gzero;
for (i=1; i<=r; i++)
for (j=1; j<=r; j++)
if (gcmp(nia,gabs(gcoeff(A,i,j),20)) == -1)
nia = gabs(gcoeff(A,i,j),20);
/* Check for the precision in matrix inversion. See paper, Lemma 2.4.2. */
if (gcmp(gmul(gadd(gmul(stoi(r),gmul(nia,eps1)),eps2),
gmul(gdeux,stoi(r))),gun) == -1)
err(talker,"not enough precision in thue");
eps3=gmul(gdeux,gmul(gmul(gsqr(stoi(r)),nia),
gadd(gmul(stoi(r),gdiv(nia,eps1)),eps2)));
myround(eps3,gun);
if (DEBUGLEVEL>=2) fprintferr("epsilon_3 -> %Z\n",eps3);
}
/* Computation of the constants c5, c7, c10, c12, c15 */
static void
ComputeConstants()
{
int k;
GEN Vect;
Vect=cgetg(r+1,t_COL); for (k=1; k<=r; k++) Vect[k]=un;
if (numroot <= r) Vect[numroot]=lsub(gun,gdeg);
Delta=gmul(A,Vect);
c5=(GEN)(gabs(Delta,Prec)[Vecmax(gabs(Delta,Prec),r)]);
c5=myround(c5,gun); c7=mulsr(r,c5);
c10=divsr(deg,c7); c13=divsr(deg,c5);
if (DEBUGLEVEL>=2)
{
fprintferr("c5 = %Z\n",c5);
fprintferr("c7 = %Z\n",c7);
fprintferr("c10 = %Z\n",c10);
fprintferr("c13 = %Z\n",c13);
}
}
/* Computation of the constants c6, c8, c9 */
static void
ComputeConstants2(GEN poly, GEN rhs)
{
GEN Vect1, tmp;
int k;
Vect1=cgetg(r+1,t_COL);
for (k=1; k<=r; k++) Vect1[k]=un;
Vect1=gmul(gabs(A,ConstPrec),Vect1);
Vect2=cgetg(r+1,t_COL);
for (k=1; k<=r; k++)
{ if (k!=numroot)
{ Vect2[k]=llog(gabs(gdiv(gsub((GEN)roo[numroot],(GEN)roo[k]),
gcoeff(MatNE,k,curne)),Prec),Prec); }
else { Vect2[k]=llog(gabs(gdiv(rhs,gmul(poleval(derivpol(poly)
,(GEN)roo[numroot]),
gcoeff(MatNE,k,curne))),Prec),Prec);
}
}
Lambda=gmul(A,Vect2);
tmp=(GEN)Vect1[Vecmax(Vect1,r)];
x2=gmax(x1,gpui(mulsr(10,mulrr(c4,tmp)),ginv(gdeg),ConstPrec));
c14=mulrr(c4,tmp);
c6=gabs((GEN)Lambda[Vecmax(gabs(Lambda,ConstPrec),r)],ConstPrec);
c6=addrr(c6,dbltor(0.1)); c6=myround(c6,gun);
c8=addrr(dbltor(1.23),mulsr(r,c6));
c11=mulrr(mulsr(2,c3),gexp(divrr(mulsr(deg,c8),c7),ConstPrec));
tmp=gexp(divrr(mulsr(deg,c6),c5),ConstPrec);
c12=mulrr(mulsr(2,c3),tmp);
c15=mulsr(2,mulrr(c14,tmp));
if (DEBUGLEVEL>=2)
{
fprintferr("c6 = %Z\n",c6);
fprintferr("c8 = %Z\n",c8);
fprintferr("c11 = %Z\n",c11);
fprintferr("c12 = %Z\n",c12);
fprintferr("c14 = %Z\n",c14);
fprintferr("c15 = %Z\n",c15);
}
}
/* Computation of the logarithmic heights of the units */
static GEN
Logarithmic_Height(int s)
{
int i;
GEN LH=gun,mat;
mat=gabs(MatFU,Prec);
for (i=1; i<=deg; i++)
LH = gmul(LH,gmax(gun,gabs(gcoeff(mat,i,s),Prec)));
return gmul(gdeux,gdiv(glog(LH,Prec),gdeg));
}
/* Computation of the matrix ((\sigma_i(\eta_j)) used for M_1 and
the computation of logarithmic heights */
static void
Compute_Fund_Units(GEN uf)
{
int i,j;
MatFU=cgetg(r+1,t_MAT);
for (i=1; i<=r; i++)
MatFU[i]=lgetg(deg+1,t_COL);
for (i=1; i<=r; i++)
{
if (typ(uf[i])!=t_POL) err(talker,"incorrect system of units");
for (j=1; j<=deg; j++)
coeff(MatFU,j,i)=(long)poleval((GEN)uf[i],(GEN)roo[j]);
}
}
/* Computation of the conjugates of the solutions to the norm equation */
static void
Conj_Norm_Eq(GEN ne, GEN *Hmu)
{
int p,q,nesol,x;
nesol=lg(ne); MatNE=(GEN)cgetg(nesol,t_MAT); (*Hmu)=cgetg(nesol,t_COL);
for (p=1; p<nesol; p++) { MatNE[p]=lgetg(deg+1,t_COL); (*Hmu)[p]=un; }
for (q=1; q<nesol; q++)
{
x=typ(ne[q]);
if (x!=t_INT && x!=t_POL)
err(talker,"incorrect solutions of norm equation");
for (p=1; p<=deg; p++)
{
coeff(MatNE,p,q)=(long)poleval((GEN)ne[q],(GEN)roo[p]);
/* Logarithmic height of the solutions of the norm equation */
(*Hmu)[q]=lmul((GEN)(*Hmu)[q],gmax(gun,gabs(gcoeff(MatNE,p,q),Prec)));
}
}
for (q=1; q<nesol; q++)
(*Hmu)[q]=ldiv(glog((GEN)(*Hmu)[q],Prec),gdeg);
}
/* Compute Baker's bound c11, and B_0, the bound for the b_i's .*/
static void
Baker(GEN ALH, GEN Hmu)
{
GEN c9=gun, gbak, hb0=gzero;
int k,i1,i2;
gbak = gmul(gmul(gdeg,gsub(gdeg,gun)),gsub(gdeg,gdeux));
switch (numroot) {
case 1: i1=2; i2=3; break;
case 2: i1=1; i2=3; break;
case 3: i1=1; i2=2; break;
default: i1=1; i2=2; break;
}
/* For the symbols used for the computation of c11, see paper, Thm 2.3.1 */
/* Compute h_1....h_r */
for (k=1; k<=r; k++)
{
c9=gmul(c9,gmax((GEN)ALH[k],
gmax(ginv(gbak),
gdiv(gabs(glog(gdiv(gcoeff(MatFU,i1,k),
gcoeff(MatFU,i2,k)),
Prec),Prec),gbak))));
}
/* Compute a bound for the h_0 */
hb0=gadd(gmul(stoi(4),halpha),gadd(gmul(gdeux,(GEN)Hmu[curne]),
gmul(gdeux,glog(gdeux,Prec))));
hb0=gmax(hb0,gmax(ginv(gbak),
gdiv(gabs(glog(gdiv(gmul(gsub((GEN)roo[numroot],
(GEN)roo[i2]),
gcoeff(MatNE,i1,curne)),
gmul(gsub((GEN)roo[numroot],
(GEN)roo[i1]),
gcoeff(MatNE,i2,curne))),
Prec),Prec),gbak)));
c9=gmul(c9,hb0);
/* Multiply c9 by the "constant" factor */
c9=gmul(c9,gmul(gmul(gmul(stoi(18),mppi(Prec)),gpui(stoi(32),stoi(4+r),Prec)),
gmul(gmul(mpfact(r+3), gpowgs(gmul(gbak,stoi(r+2)), r+3)),
glog(gmul(gdeux,gmul(gbak,stoi(r+2))),Prec))));
c9=myround(c9,gun);
/* Compute B0 according to Lemma 2.3.3 */
B0=gmax(gexp(gun,Prec),
mulsr(2,divrr(addrr(mulrr(c9,glog(divrr(c9,c10),ConstPrec)),
glog(c11,ConstPrec)),c10)));
if (DEBUGLEVEL>=2) fprintferr("Baker -> %Z\nB0 -> %Z\n",c9,B0);
}
/* Compute delta and lambda */
static void
Create_CF_Values(int i1, int i2, GEN *errdelta)
{
GEN eps5;
if (r>1)
{
delta=divrr((GEN)Delta[i2],(GEN)Delta[i1]);
eps5=mulrs(eps3,r);
*errdelta=mulrr(addsr(1,delta),
divrr(eps5,subrr(gabs((GEN)Delta[i1],ConstPrec),eps5)));
lambda=gdiv(gsub(gmul((GEN)Delta[i2],(GEN)Lambda[i1]),
gmul((GEN)Delta[i1],(GEN)Lambda[i2])),
(GEN)Delta[i1]);
}
else
{
GEN Pi2 = gmul2n(mppi(Prec),1);
delta=gdiv(gcoeff(MatFU,2,1),gcoeff(MatFU,3,1));
delta=gdiv(garg(delta,Prec),Pi2);
*errdelta=gdiv(gpui(gdeux,stoi(-(Prec-2)*32+1),Prec),
gabs(gcoeff(MatFU,2,1),Prec));
lambda=gmul(gdiv(gsub((GEN)roo[numroot],(GEN)roo[2]),
gsub((GEN)roo[numroot],(GEN)roo[3])),
gdiv(gcoeff(MatNE,3,curne),gcoeff(MatNE,2,curne)));
lambda=gdiv(garg(lambda,Prec),Pi2);
}
if (DEBUGLEVEL>=2) outerr(*errdelta);
}
/* Try to reduce the bound through continued fractions; see paper. */
static int
CF_First_Pass(GEN kappa, GEN errdelta)
{
GEN q,ql,qd,l0;
if (gcmp(gmul(dbltor(0.1),gsqr(mulri(B0,kappa))),ginv(errdelta))==1)
{
return -1;
}
q=denom(bestappr(delta, mulir(kappa, B0))); ql=mulir(q,lambda);
qd=gmul(q,delta); ql=gabs(subri(ql, ground(ql)),Prec);
qd=gabs(mulrr(subri(qd,ground(qd)),B0),Prec);
l0=subrr(ql,addrr(qd,divri(dbltor(0.1),kappa)));
if (signe(l0) <= 0)
{
if (DEBUGLEVEL>=2)
fprintferr("CF_First_Pass failed. Trying again with larger kappa\n");
return 0;
}
if (r>1)
B0=divrr(glog(divrr(mulir(q,c15),l0),ConstPrec),c13);
else
B0=divrr(glog(divrr(mulir(q,c11),mulrr(l0,gmul2n(mppi(ConstPrec),1))),ConstPrec),c10);
if (DEBUGLEVEL>=2)
fprintferr("CF_First_Pass successful !!\nB0 -> %Z\n",B0);
return 1;
}
/* Check whether a "potential" solution is truly a solution. */
static void
Check_Solutions(GEN xmay1, GEN xmay2, GEN poly, GEN rhs)
{
GEN xx1, xx2, yy1, yy2, zz, u;
yy1=ground(greal(gdiv(gsub(xmay2,xmay1), gsub((GEN)roo[1],(GEN)roo[2]))));
yy2=gneg(yy1);
xx1=greal(gadd(xmay1,gmul((GEN)roo[1],yy1)));
xx2=gneg(xx1);
if (gcmp(distoZ(xx1),dbltor(0.0001))==-1)
{
xx1=ground(xx1);
if (!gcmp0(yy1))
{
if (gegal(gmul(poleval(poly,gdiv(xx1,yy1)),
gpui(yy1,gdeg,Prec)),(GEN)rhs))
{
zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
u[1]=(long)xx1; u[2]=(long)yy1; zz[1]=(long)u;
if (_thue_new(u)) SOL=concatsp(SOL,zz);
}
}
else if (gegal(gpui(xx1,gdeg,Prec),(GEN)rhs))
{
zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
u[1]=(long)xx1; u[2]=(long)gzero; zz[1]=(long)u;
if (_thue_new(u)) SOL=concatsp(SOL,zz);
}
xx2=ground(xx2);
if (!gcmp0(yy2))
{
if (gegal(gmul(poleval(poly,gdiv(xx2,yy2)),
gpui(yy2,gdeg,Prec)),(GEN)rhs))
{
zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
u[1]=(long)xx2; u[2]=(long)yy2; zz[1]=(long)u;
if (_thue_new(u)) SOL=concatsp(SOL,zz);
}
}
else if (gegal(gpui(xx2,gdeg,Prec),(GEN)rhs))
{
zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
u[1]=(long)xx2; u[2]=(long)gzero; zz[1]=(long)u;
if (_thue_new(u)) SOL=concatsp(SOL,zz);
}
}
}
static GEN
GuessQi(GEN *q1, GEN *q2, GEN *q3)
{
GEN C, Lat, eps;
C=gpui(stoi(10),stoi(10),Prec);
Lat=idmat(3); mael(Lat,1,3)=(long)C; mael(Lat,2,3)=lround(gmul(C,delta));
mael(Lat,3,3)=lround(gmul(C,lambda));
Lat=lllint(Lat);
*q1=gmael(Lat,1,1); *q2=gmael(Lat,1,2); *q3=gmael(Lat,1,3);
eps=gabs(gadd(gadd(*q1,gmul(*q2,delta)),gmul(*q3,lambda)),Prec);
return(eps);
}
static void
TotRat()
{
err(bugparier,"thue (totally rational case)");
}
/* Check for solutions under a small bound (see paper) */
static void
Check_Small(int bound, GEN poly, GEN rhs)
{
GEN interm, xx, zz, u, maxr, tmp, ypot, xxn, xxnm1, yy;
long av = avma, lim = stack_lim(av,1);
int x, j, bsupy, y;
double bndyx;
maxr=gabs((GEN)roo[1],Prec); for(j=1; j<=deg; j++)
{ if(gcmp(tmp=gabs((GEN)roo[j],Prec),maxr)==1) maxr=tmp; }
bndyx=gtodouble(gadd(gpui(gabs(rhs,Prec),ginv(gdeg),Prec),maxr));
for (x=-bound; x<=bound; x++)
{
if (low_stack(lim,stack_lim(av,1)))
{
if(DEBUGMEM>1) err(warnmem,"Check_small");
SOL = gerepileupto(av, gcopy(SOL));
}
if (x==0)
{
ypot=gmul(stoi(gsigne(rhs)),ground(gpui(gabs(rhs,0),ginv(gdeg),Prec)));
if (gegal(gpui(ypot,gdeg,0),rhs))
{
zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
u[1]=(long)ypot; u[2]=(long)gzero; zz[1]=(long)u;
if (_thue_new(u)) SOL=concatsp(SOL,zz);
}
if (gegal(gpui(gneg(ypot),gdeg,0),rhs))
{
zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
u[1]=(long)gneg(ypot); u[2]=(long)gzero; zz[1]=(long)u;
if (_thue_new(u)) SOL=concatsp(SOL,zz);
}
}
else
{
bsupy=(int)(x>0?bndyx*x:-bndyx*x);
xx=stoi(x); xxn=gpui(xx,gdeg,Prec);
interm=gsub((GEN)rhs,gmul(xxn,(GEN)poly[2]));
/* Verifier ... */
xxnm1=xxn; j=2;
while(gegal(interm,gzero))
{
xxnm1=gdiv(xxnm1,xx);
interm=gmul((GEN)poly[++j],xxnm1);
}
for(y=-bsupy; y<=bsupy; y++)
{
yy=stoi(y);
if(y==0) {
if (gegal(gmul((GEN)poly[2],xxn),rhs))
{
zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
u[1]=(long)gzero; u[2]=(long)xx; zz[1]=(long)u;
if (_thue_new(u)) SOL=concatsp(SOL,zz);
}
}
else if(gegal(gmod(interm,yy),gzero))
if(gegal(poleval(poly,gdiv(yy,xx)),gdiv(rhs,xxn)))
/* Remplacer par un eval *homogene* */
{
zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
u[1]=(long)yy; u[2]=(long)xx; zz[1]=(long)u;
if (_thue_new(u)) SOL=concatsp(SOL,zz);
}
}
}
}
}
/* Computes [x]! */
static double
fact(double x)
{
double ft=1.0;
x=(int)x; while (x>1) { ft*=x; x--; }
return ft ;
}
/* From a polynomial and optionally a system of fundamental units for the
* field defined by poly, computes all the relevant constants needed to solve
* the equation P(x,y)=a whenever the solutions of N_{ K/Q }(x)=a. Returns a
* "tnf" structure containing
* 1) the polynomial
* 2) the bnf (used to solve the norm equation)
* 3) roots, with presumably enough precision
* 4) The logarithmic heights of units
* 5) The matrix of conjugates of units
* 6) its inverse
* 7) a few technical constants
*/
GEN
thueinit(GEN poly, long flag, long prec)
{
GEN thueres,ALH,csts,c0;
long av,tetpil,k,st;
double d,dr;
av=avma;
uftot=0;
if (checktnf(poly)) { uftot=(GEN)poly[2]; poly=(GEN)poly[1]; }
else
if (typ(poly)!=t_POL) err(notpoler,"thueinit");
if (!gisirreducible(poly)) err(redpoler,"thueinit");
st=sturm(poly);
if (st)
{
dr=(double)((st+lgef(poly)-5)>>1);
d=(double)lgef(poly)-3; d=d*(d-1)*(d-2);
/* Try to guess the precision by approximating Baker's bound.
* Note that the guess is most of the time pretty generous,
* ie 10 to 30 decimal digits above what is *really* necessary.
* Note that the limiting step is the reduction. See paper.
*/
Prec=3 + (long)((5.83 + (dr+4)*5 + log(fact(dr+3)) + (dr+3)*log(dr+2) +
(dr+3)*log(d) + log(log(2*d*(dr+2))) + (dr+1)) / 10.);
ConstPrec=4;
if (Prec<prec) Prec = prec;
if (!checktnf(poly)) inithue(poly,flag);
thueres=cgetg(8,t_VEC);
thueres[1]=(long)poly;
thueres[2]=(long)uftot;
thueres[3]=(long)roo;
Compute_Fund_Units(gmael(uftot,8,5));
ALH=cgetg(r+1,t_COL);
for (k=1; k<=r; k++) ALH[k]=(long)Logarithmic_Height(k);
thueres[4]=(long)ALH;
thueres[5]=(long)MatFU;
T_A_Matrices();
thueres[6]=(long)A;
csts=cgetg(7,t_VEC);
csts[1]=(long)c1; csts[2]=(long)c2; csts[3]=(long)halpha;
csts[4]=(long)x0; csts[5]=(long)eps3;
csts[6]=(long)stoi(Prec);
thueres[7]=(long)csts; tetpil=avma;
return gerepile(av,tetpil,gcopy(thueres));
}
thueres=cgetg(3,t_VEC); c0=gun; Prec=4;
roo=roots(poly,Prec);
for (k=1; k<lg(roo); k++)
c0=gmul(c0, gimag((GEN)roo[k]));
c0=ginv(gabs(c0,Prec));
thueres[1]=(long)poly; thueres[2]=(long)c0;
tetpil=avma; return gerepile(av,tetpil,gcopy(thueres));
}
/* Given a tnf structure as returned by thueinit, and a right-hand-side and
* optionally the solutions to the norm equation, returns the solutions to
* the Thue equation F(x,y)=a
*/
GEN
thue(GEN thueres, GEN rhs, GEN ne)
{
GEN Kstart,Old_B0,ALH,errdelta,Hmu,c0,poly,csts,bd;
GEN xmay1,xmay2,b,zp1,tmp,q1,q2,q3,ep;
long Nb_It=0,upb,bi1,av,tetpil,i1,i2,i, flag,cf,fs;
if (!checktnf(thueres)) err(talker,"not a tnf in thue");
av=avma; SOL=cgetg(1,t_VEC);
if (lg(thueres)==8)
{
poly=(GEN)thueres[1]; deg=lgef(poly)-3; gdeg=stoi(deg);
uftot=(GEN)thueres[2];
roo=gcopy((GEN)thueres[3]);
ALH=(GEN)thueres[4];
MatFU=gcopy((GEN)thueres[5]);
A=gcopy((GEN)thueres[6]);
csts=(GEN)thueres[7];
if (!ne) ne = bnfisintnorm(uftot, rhs);
if (lg(ne)==1) { avma=av; return cgetg(1,t_VEC); }
c1=gmul(gabs(rhs,Prec), (GEN)csts[1]); c2=(GEN)csts[2];
halpha=(GEN)csts[3];
if (t)
x0 = gmul((GEN)csts[4],gpui(gabs(rhs,Prec), ginv(gdeg), Prec));
eps3=(GEN)csts[5];
Prec=gtolong((GEN)csts[6]);
b=cgetg(r+1,t_COL);
tmp=divrr(c1,c2);
x1=gmax(x0,gpui(mulsr(2,tmp),ginv(gdeg),ConstPrec));
if(DEBUGLEVEL >=2) fprintferr("x1 -> %Z\n",x1);
c3=mulrr(dbltor(1.39),tmp);
c4=mulsr(deg-1,c3);
for (numroot=1; numroot<=s; numroot++)
{
ComputeConstants();
Conj_Norm_Eq(ne,&Hmu);
for (curne=1; curne<lg(ne); curne++)
{
ComputeConstants2(poly,rhs);
Baker(ALH,Hmu);
i1=Vecmax(gabs(Delta,Prec),r);
if (i1!=1) i2=1; else i2=2;
do
{
fs=0;
Create_CF_Values(i1,i2,&errdelta);
if (DEBUGLEVEL>=2) fprintferr("Entering CF\n");
Old_B0=gmul(B0,gdeux);
/* Try to reduce B0 while
* 1) there was less than 10 reductions
* 2) the previous reduction improved the bound of more than 0.1.
*/
while (Nb_It<10 && gcmp(Old_B0,gadd(B0,dbltor(0.1))) == 1 && fs<2)
{
cf=0; Old_B0=B0; Nb_It++; Kstart=stoi(10);
while (!fs && CF_First_Pass(Kstart,errdelta) == 0 && cf < 8 )
{
cf++;
Kstart=mulis(Kstart,10);
}
if ( CF_First_Pass(Kstart,errdelta) == -1 )
{ ne = gerepileupto(av, gcopy(ne)); Prec+=10;
if(DEBUGLEVEL>=2) fprintferr("Increasing precision\n");
thueres=thueinit(thueres,0,Prec);
return(thue(thueres,rhs,ne)); }
if (cf==8 || fs) /* Semirational or totally rational case */
{
if (!fs)
{ ep=GuessQi(&q1,&q2,&q3); }
bd=gmul(denom(bestappr(delta,gadd(B0,gabs(q2,Prec))))
,delta);
bd=gabs(gsub(bd,ground(bd)),Prec);
if(gcmp(bd,ep)==1 && !gegal(q3, gzero))
{
fs=1;
B0=gdiv(glog(gdiv(gmul(q3,c15),gsub(bd,ep)), Prec),c13);
if (DEBUGLEVEL>=2)
fprintferr("Semirat. reduction ok. B0 -> %Z\n",B0);
}
else fs=2;
}
else fs=0;
Nb_It=0; B0=gmin(Old_B0,B0); upb=gtolong(gceil(B0));
}
i2++; if (i2==i1) i2++;
}
while (fs == 2 && i2 <= r);
if (fs == 2) TotRat();
if (DEBUGLEVEL >=2) fprintferr("x2 -> %Z\n",x2);
/* For each possible value of b_i1, compute the b_i's
* and 2 conjugates of x-\alpha y. Then check.
*/
zp1=dbltor(0.01);
x3=gmax(x2,gpui(mulsr(2,divrr(c14,zp1)),ginv(gdeg),ConstPrec));
for (bi1=-upb; bi1<=upb; bi1++)
{
flag=1;
xmay1=gun; xmay2=gun;
for (i=1; i<=r; i++)
{
b[i]=(long)gdiv(gsub(gmul((GEN)Delta[i],stoi(bi1)),
gsub(gmul((GEN)Delta[i],(GEN)Lambda[i1]),
gmul((GEN)Delta[i1],(GEN)Lambda[i]))),
(GEN)Delta[i1]);
if (gcmp(distoZ((GEN)b[i]),zp1)==1) { flag=0; break; }
}
if(flag)
{
for(i=1; i<=r; i++)
{
b[i]=lround((GEN)b[i]);
xmay1=gmul(xmay1,gpui(gcoeff(MatFU,1,i),(GEN)b[i],Prec));
xmay2=gmul(xmay2,gpui(gcoeff(MatFU,2,i),(GEN)b[i],Prec));
}
xmay1=gmul(xmay1,gcoeff(MatNE,1,curne));
xmay2=gmul(xmay2,gcoeff(MatNE,2,curne));
Check_Solutions(xmay1,xmay2,poly,rhs);
}
}
}
}
if (DEBUGLEVEL>=2) fprintferr("Checking for small solutions\n");
Check_Small((int)(gtodouble(x3)),poly,rhs);
tetpil=avma; return gerepile(av,tetpil,gcopy(SOL));
}
/* Case s=0. All the solutions are "small". */
Prec=DEFAULTPREC; poly=(GEN)thueres[1]; deg=lgef(poly)-3;
gdeg=stoi(deg); c0=(GEN)thueres[2];
roo=roots(poly,Prec);
Check_Small(gtolong(ground(gpui(gmul((GEN)poly[2],gdiv(gabs(rhs,Prec),c0)),
ginv(gdeg),Prec) )), poly, rhs);
tetpil=avma; return gerepile(av,tetpil,gcopy(SOL));
}
static GEN *Relations; /* primes above a, expressed on generators of Cl(K) */
static GEN *Partial; /* list of vvectors, Partial[i] = Relations[1..i] * u[1..i] */
static GEN *gen_ord; /* orders of generators of Cl(K) given in bnf */
static long *f; /* f[i] = f(Primes[i]/p), inertial degree */
static long *n; /* a = prod p^{ n_p }. n[i]=n_p if Primes[i] divides p */
static long *inext; /* index of first P above next p, 0 if p is last */
static long *S; /* S[i] = n[i] - sum_{ 1<=k<=i } f[k].u[k] */
static long *u; /* We want principal ideals I = prod Primes[i]^u[i] */
static long **normsol; /* lists of copies of the u[] which are solutions */
static long Nprimes; /* length(Relations) = #{max ideal above divisors of a} */
static long sindex, smax; /* current index in normsol; max. index */
/* u[1..i] has been filled. Norm(u) is correct.
* Check relations in class group then save it.
*/
static void
test_sol(long i)
{
long k,*sol;
if (Partial)
{
long av=avma;
for (k=1; k<lg(Partial[1]); k++)
if ( signe(modii( (GEN)Partial[i][k], gen_ord[k] )) )
{ avma=av; return; }
avma=av;
}
if (sindex == smax)
{
long new_smax = smax << 1;
long **new_normsol = (long **) new_chunk(new_smax+1);
for (k=1; k<=smax; k++) new_normsol[k] = normsol[k];
normsol = new_normsol; smax = new_smax;
}
sol = cgetg(Nprimes+1,t_VECSMALL);
normsol[++sindex] = sol;
for (k=1; k<=i; k++) sol[k]=u[k];
for ( ; k<=Nprimes; k++) sol[k]=0;
if (DEBUGLEVEL > 2)
{
fprintferr("sol = %Z\n",sol);
if (Partial) fprintferr("Partial = %Z\n",Partial);
flusherr();
}
}
static void
fix_Partial(long i)
{
long k, av = avma;
for (k=1; k<lg(Partial[1]); k++)
addiiz(
(GEN) Partial[i-1][k],
mulis((GEN) Relations[i][k], u[i]),
(GEN) Partial[i][k]
);
avma=av;
}
/* Recursive loop. Suppose u[1..i] has been filled
* Find possible solutions u such that, Norm(prod Prime[i]^u[i]) = a, taking
* into account:
* 1) the relations in the class group if need be.
* 2) the factorization of a.
*/
static void
isintnorm_loop(long i)
{
if (S[i] == 0) /* sum u[i].f[i] = n[i], do another prime */
{
int k;
if (inext[i] == 0) { test_sol(i); return; }
/* there are some primes left */
if (Partial) gaffect(Partial[i], Partial[inext[i]-1]);
for (k=i+1; k < inext[i]; k++) u[k]=0;
i=inext[i]-1;
}
else if (i == inext[i]-2 || i == Nprimes-1)
{
/* only one Prime left above prime; change prime, fix u[i+1] */
if (S[i] % f[i+1]) return;
i++; u[i] = S[i-1] / f[i];
if (Partial) fix_Partial(i);
if (inext[i]==0) { test_sol(i); return; }
}
i++; u[i]=0;
if (Partial) gaffect(Partial[i-1], Partial[i]);
if (i == inext[i-1])
{ /* change prime */
if (inext[i] == i+1 || i == Nprimes) /* only one Prime above p */
{
S[i]=0; u[i] = n[i] / f[i]; /* we already know this is exact */
if (Partial) fix_Partial(i);
}
else S[i] = n[i];
}
else S[i] = S[i-1]; /* same prime, different Prime */
for(;;)
{
isintnorm_loop(i); S[i]-=f[i]; if (S[i]<0) break;
if (Partial)
gaddz(Partial[i], Relations[i], Partial[i]);
u[i]++;
}
}
static void
get_sol_abs(GEN bnf, GEN a, GEN **ptPrimes)
{
GEN dec, fact, primes, *Primes, *Fact;
long *gcdlist, gcd,nprimes,Ngen,i,j;
if (gcmp1(a))
{
long *sol = new_chunk(Nprimes+1);
sindex = 1; normsol = (long**) new_chunk(2);
normsol[1] = sol; for (i=1; i<=Nprimes; i++) sol[i]=0;
return;
}
fact=factor(a); primes=(GEN)fact[1];
nprimes=lg(primes)-1; sindex = 0;
gen_ord = (GEN *) mael3(bnf,8,1,2); Ngen = lg(gen_ord)-1;
Fact = (GEN*) new_chunk(1+nprimes);
gcdlist = new_chunk(1+nprimes);
Nprimes=0;
for (i=1; i<=nprimes; i++)
{
long ldec;
dec = primedec(bnf,(GEN)primes[i]); ldec = lg(dec)-1;
gcd = itos(gmael(dec,1,4));
/* check that gcd_{P|p} f_P | n_p */
for (j=2; gcd>1 && j<=ldec; j++)
gcd = cgcd(gcd,itos(gmael(dec,j,4)));
gcdlist[i]=gcd;
if (gcd != 1 && smodis(gmael(fact,2,i),gcd))
{
if (DEBUGLEVEL>2)
{ fprintferr("gcd f_P does not divide n_p\n"); flusherr(); }
return;
}
Fact[i] = dec; Nprimes += ldec;
}
f = new_chunk(1+Nprimes); u = new_chunk(1+Nprimes);
n = new_chunk(1+Nprimes); S = new_chunk(1+Nprimes);
inext = new_chunk(1+Nprimes);
Primes = (GEN*) new_chunk(1+Nprimes);
*ptPrimes = Primes;
if (Ngen)
{
Partial = (GEN*) new_chunk(1+Nprimes);
Relations = (GEN*) new_chunk(1+Nprimes);
}
else /* trivial class group, no relations to check */
Relations = Partial = NULL;
j=0;
for (i=1; i<=nprimes; i++)
{
long k,lim,v, vn=itos(gmael(fact,2,i));
gcd = gcdlist[i];
dec = Fact[i]; lim = lg(dec);
v = (i==nprimes)? 0: j + lim;
for (k=1; k < lim; k++)
{
j++; Primes[j] = (GEN)dec[k];
f[j] = itos(gmael(dec,k,4)) / gcd;
n[j] = vn / gcd; inext[j] = v;
if (Partial)
Relations[j] = isprincipal(bnf,Primes[j]);
}
}
if (Partial)
{
for (i=1; i <= Nprimes; i++)
if (!gcmp0(Relations[i])) break;
if (i > Nprimes) Partial = NULL; /* all ideals dividing a are principal */
}
if (Partial)
for (i=0; i<=Nprimes; i++) /* Partial[0] needs to be initialized */
{
Partial[i]=cgetg(Ngen+1,t_COL);
for (j=1; j<=Ngen; j++)
{
Partial[i][j]=lgeti(4);
gaffect(gzero,(GEN) Partial[i][j]);
}
}
smax=511; normsol = (long**) new_chunk(smax+1);
S[0]=n[1]; inext[0]=1; isintnorm_loop(0);
}
static long
get_unit_1(GEN bnf, GEN pol, GEN *unit)
{
GEN fu;
long j;
if (DEBUGLEVEL > 2)
fprintferr("looking for a fundamental unit of norm -1\n");
*unit = gmael3(bnf,8,4,2); /* torsion unit */
if (signe( gnorm(gmodulcp(*unit,pol)) ) < 0) return 1;
fu = gmael(bnf,8,5); /* fundamental units */
for (j=1; j<lg(fu); j++)
{
*unit = (GEN)fu[j];
if (signe( gnorm(gmodulcp(*unit,pol)) ) < 0) return 1;
}
return 0;
}
GEN
bnfisintnorm(GEN bnf, GEN a)
{
GEN nf,pol,res,unit,x,id, *Primes;
long av = avma, tetpil,sa,i,j,norm_1;
bnf = checkbnf(bnf); nf = (GEN)bnf[7]; pol = (GEN)nf[1];
if (typ(a)!=t_INT)
err(talker,"expected an integer in bnfisintnorm");
sa = signe(a);
if (sa == 0) { res=cgetg(2,t_VEC); res[1]=zero; return res; }
if (gcmp1(a)) { res=cgetg(2,t_VEC); res[1]=un; return res; }
get_sol_abs(bnf, absi(a), &Primes);
res = cgetg(1,t_VEC); unit = NULL;
for (i=1; i<=sindex; i++)
{
id = gun; x = normsol[i];
for (j=1; j<=Nprimes; j++) /* compute prod Primes[i]^u[i] */
if (x[j])
{
GEN id2 = Primes[j];
if (x[j] != 1) id2 = idealpow(nf,id2, stoi(x[j]));
id = idealmul(nf,id,id2);
}
x = (GEN) isprincipalgen(bnf,id)[2];
x = gmul((GEN)nf[7],x); /* x possible solution */
if (signe(gnorm(gmodulcp(x,(GEN)nf[1]))) != sa)
{
if (! unit) norm_1 = get_unit_1(bnf,pol,&unit);
if (norm_1) x = gmul(unit,x);
else
{
if (DEBUGLEVEL > 2)
fprintferr("%Z eliminated because of sign\n",x);
x = NULL;
}
}
if (x) res = concatsp(res, gmod(x,pol));
}
tetpil=avma; return gerepile(av,tetpil,gcopy(res));
}
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to thue.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari / src / modules