[BACK]Return to thue.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / pari-2.2 / src / modules

Diff for /OpenXM_contrib/pari-2.2/src/modules/Attic/thue.c between version 1.1 and 1.2

version 1.1, 2001/10/02 11:17:12 version 1.2, 2002/09/11 07:27:06
Line 13  Check the License for details. You should have receive
Line 13  Check the License for details. You should have receive
 with the package; see the file 'COPYING'. If not, write to the Free Software  with the package; see the file 'COPYING'. If not, write to the Free Software
 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
   
   #include "pari.h"
   
 /* Thue equation solver. In all the forthcoming remarks, "paper"  /* Thue equation solver. In all the forthcoming remarks, "paper"
  * designs the paper "Thue Equations of High Degree", by Yu. Bilu and   * designs the paper "Thue Equations of High Degree", by Yu. Bilu and
  * G. Hanrot, J. Number Theory (1996). Note that numbering of the constants   * G. Hanrot, J. Number Theory (1996). Note that numbering of the constants
  * is that of Hanrot's thesis rather than that of the paper.   * is that of Hanrot's thesis rather than that of the paper.
  * The last part of the program (bnfisintnorm) was written by K. Belabas.   * The last part of the program (bnfisintnorm) was written by K. Belabas. */
  */  
 #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 */  /* Check whether tnf is a valid structure */
 static int  static int
 checktnf(GEN tnf)  checktnf(GEN tnf)
 {  {
     int n, R, S, T;
   if (typ(tnf)!=t_VEC || (lg(tnf)!=8 && lg(tnf)!=3)) return 0;    if (typ(tnf)!=t_VEC || (lg(tnf)!=8 && lg(tnf)!=3)) return 0;
   if (typ(tnf[1])!=t_POL) return 0;    if (typ(tnf[1]) != t_POL) return 0;
   if (lg(tnf) != 8) return 1; /* s=0 */    if (lg(tnf) != 8) return 1; /* S=0 */
   
   deg=degpol(tnf[1]);    n = degpol(tnf[1]);
   if (deg<=2) err(talker,"invalid polynomial in thue (need deg>2)");    if (n <= 2) err(talker,"invalid polynomial in thue (need n>2)");
   s=sturm((GEN)tnf[1]); t=(deg-s)>>1; r=s+t-1;    S = sturm((GEN)tnf[1]); T = (n-S)>>1; R = S+T-1;
   (void)checkbnf((GEN)tnf[2]);    (void)checkbnf((GEN)tnf[2]);
   if (typ(tnf[3]) != t_COL || lg(tnf[3]) != deg+1) return 0;    if (typ(tnf[3]) != t_COL || lg(tnf[3]) != n+1) return 0;
   if (typ(tnf[4]) != t_COL || lg(tnf[4]) != r+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    if (typ(tnf[5]) != t_MAT || lg(tnf[5])!=R+1 || lg(gmael(tnf,5,1)) != n+1) return 0;
       || 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[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;    if (typ(tnf[7]) != t_VEC || lg(tnf[7]) != 7) return 0;
   return 1;    return 1;
 }  }
   
 static GEN distoZ(GEN z)  static GEN
   distoZ(GEN z)
 {  {
   GEN p1=gfrac(z);    GEN p1 = gfrac(z);
   return gmin(p1, gsub(gun,p1));    return gmin(p1, gsub(gun,p1));
 }  }
   
 /* Check whether a solution has already been found */  /* Compensates rounding errors for computation/display of the constants.
 static int   * Round up if dir > 0, down otherwise */
 _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  static GEN
 myround(GEN Cst, GEN upd)  myround(GEN x, long dir)
 {  {
   return gmul(Cst, gadd(gun, gmul(upd,gpuigs(stoi(10),-10))));    GEN eps = gpowgs(stoi(dir > 0? 10: -10), -10);
     return gmul(x, gadd(gun, eps));
 }  }
   
 /* Returns the index of the largest element in a vector */  /* Returns the index of the largest element in a vector */
 static int  static GEN
 Vecmax(GEN Vec, int size)  _Vecmax(GEN Vec, int *ind)
 {  {
   int k, tno = 1;    int k, tno = 1, l = lg(Vec);
   GEN tmax = (GEN)Vec[1];    GEN tmax = (GEN)Vec[1];
   for (k=2; k<=size; k++)    for (k = 2; k < l; k++)
     if (gcmp((GEN)Vec[k],tmax)==1) { tmax=(GEN)Vec[k]; tno=k; }      if (gcmp((GEN)Vec[k],tmax) > 0) { tmax = (GEN)Vec[k]; tno = k; }
   return tno;    if (ind) *ind = tno;
     return tmax;
 }  }
   
 /* Performs basic computations concerning the equation: buchinitfu, c1, c2 */  static GEN
 static void  Vecmax(GEN v) { return _Vecmax(v, NULL); }
 inithue(GEN poly, long flag)  
 {  
   GEN roo2, tmp, gpmin, de;  
   int k,j;  
   
   x0=gzero; x1=gzero; deg=degpol(poly); gdeg=stoi(deg);  static int
   Vecmaxind(GEN v) { int i; (void)_Vecmax(v, &i); return i; }
   
   if (!uftot)  static GEN
   {  tnf_get_roots(GEN poly, long prec, int S, int T)
     uftot=bnfinit0(poly,1,NULL,Prec);  {
     if (uftot != checkbnf_discard(uftot))    GEN R0 = roots(poly, prec), R = cgetg(lg(R0), t_COL);
       err(talker,"non-monic polynomial in thue");    int k;
     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<=S; k++) R[k] = lreal((GEN)R0[k]);
   for (k=1; k<=deg; k++)    /* swap roots to get the usual order */
     halpha = gmul(halpha,gmax(gun,gabs((GEN)roo[k],Prec)));    for (k=1; k<=T; k++)
   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);      R[k+S]  = R0[2*k+S-1];
     if (gcmp(tmp,c1)==-1) c1=tmp;      R[k+S+T]= R0[2*k+S];
   }    }
     return R;
   }
   
   c1=gdiv(gpui(gdeux,gsub(gdeg,gun),Prec),c1); c1=myround(c1,gun);  /* Computation of the logarithmic height of x (given by embeddings) */
   c2=gabs(gsub((GEN)roo[1],(GEN)roo[2]),Prec);  static GEN
   LogHeight(GEN x, long prec)
   for (k=1; k<=deg; k++)  {
     for (j=k+1; j<=deg; j++)    int i, n = lg(x)-1;
     {    GEN LH = gun;
       tmp=gabs(gsub((GEN)roo[j],(GEN)roo[k]),Prec);    for (i=1; i<=n; i++) LH = gmul(LH, gmax(gun, gabs((GEN)x[i], prec)));
       if (gcmp(c2,tmp)==1) c2=tmp;    return gdivgs(glog(LH,prec), n);
     }  
   
   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) */  /* x^(1/n) */
 static void  static GEN
 T_A_Matrices(void)  mpsqrtn(GEN x, long n)
 {  {
   GEN mask, eps1, eps2, nia, m1, IntM;    if (typ(x) != t_REAL) err(typeer,"mpsqrtn");
   int i,j;    return mpexp(divrs(mplog(x), n));
   }
   
   m1=glog(gabs(MatFU,Prec),Prec); mask=gsub(gpui(gdeux,stoi(r),Prec),gun);  /* |x|^(1/n), x t_INT */
   m1=matextract(m1,mask,mask);  static GEN
   absisqrtn(GEN x, long n, long prec) {
   A=invmat(m1); IntM=gsub(gmul(A,m1),idmat(r));    GEN r = itor(x,prec);
     if (signe(r) < 0) r = negr(r);
   eps2=gzero;    return mpsqrtn(r, n);
   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(precer,"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 GEN
 static void  get_emb(GEN x, GEN r)
 ComputeConstants(void)  
 {  {
   int k;    int l = lg(r), i, tx;
     GEN e, y = cgetg(l, t_COL);
   
   GEN Vect;    tx = typ(x);
     if (tx != t_INT && tx != t_POL) err(typeer,"get_emb");
   Vect=cgetg(r+1,t_COL); for (k=1; k<=r; k++) Vect[k]=un;    for (i=1; i<l; i++)
   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);      e = poleval(x, (GEN)r[i]);
     fprintferr("c7 = %Z\n",c7);      if (gcmp0(e)) return NULL;
     fprintferr("c10 = %Z\n",c10);      y[i] = (long)e;
     fprintferr("c13 = %Z\n",c13);  
   }    }
     return y;
 }  }
   
 /* Computation of the constants c6, c8, c9 */  /* Computation of the conjugates (sigma_i(v_j)), and log. heights, of elts of v */
 static void  static GEN
 ComputeConstants2(GEN poly, GEN rhs)  Conj_LH(GEN v, GEN *H, GEN r, long prec)
 {  {
   GEN Vect1, tmp;    int j, l = lg(v);
   int k;    GEN e, M = (GEN)cgetg(l,t_MAT);
   
   Vect1=cgetg(r+1,t_COL);    (*H) = cgetg(l,t_COL);
   for (k=1; k<=r; k++) Vect1[k]=un;    for (j = 1; j < l; j++)
   Vect1=gmul(gabs(A,ConstPrec),Vect1);    {
       if (! (e = get_emb((GEN)v[j], r)) ) return NULL; /* FAIL */
       M[j] = (long)e;
       (*H)[j]= (long)LogHeight(e, prec);
     }
     return M;
   }
   
   static GEN abslog(GEN x, long prec) { return gabs(glog(x,prec), prec); }
   static GEN logabs(GEN x, long prec) { return glog(gabs(x,prec), prec); }
   
   Vect2=cgetg(r+1,t_COL);  /* Computation of M, its inverse A and precision check (see paper) */
   for (k=1; k<=r; k++)  static GEN
     { if (k!=numroot)  T_A_Matrices(GEN MatFU, int r, GEN *eps5, long prec)
         { Vect2[k]=llog(gabs(gdiv(gsub((GEN)roo[numroot],(GEN)roo[k]),  {
                                   gcoeff(MatNE,k,curne)),Prec),Prec); }    GEN A, p1, m1, IntM, nia, eps3, eps2, eps1 = shifti(gun, bit_accuracy(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)];    m1 = rowextract_i(vecextract_i(MatFU, 1,r), 1,r); /* minor order r */
   x2=gmax(x1,gpui(mulsr(10,mulrr(c4,tmp)),ginv(gdeg),ConstPrec));    m1 = logabs(m1,prec);
   c14=mulrr(c4,tmp);  
   
   c6=gabs((GEN)Lambda[Vecmax(gabs(Lambda,ConstPrec),r)],ConstPrec);    A = invmat(m1);
   c6=addrr(c6,dbltor(0.1)); c6=myround(c6,gun);    IntM = gsub(gmul(A,m1), idmat(r));
   
   c8=addrr(dbltor(1.23),mulsr(r,c6));    eps2 = gadd(vecmax(gabs(IntM,prec)), ginv(eps1));
   c11=mulrr(mulsr(2,c3),gexp(divrr(mulsr(deg,c8),c7),ConstPrec));    nia = vecmax(gabs(A,prec));
   
   tmp=gexp(divrr(mulsr(deg,c6),c5),ConstPrec);    /* Check for the precision in matrix inversion. See paper, Lemma 2.4.2. */
   c12=mulrr(mulsr(2,c3),tmp);    p1 = gadd(gmulsg(r, gmul(nia,eps1)), eps2);
   c15=mulsr(2,mulrr(c14,tmp));    if (gcmp(gmulgs(p1, 2*r), gun) < 0) err(precer,"thue");
   
   if (DEBUGLEVEL>=2)    p1 = gadd(gmulsg(r, gdiv(nia,eps1)), eps2);
   {    eps3 = gmul(gmulsg(2*r*r,nia), p1);
     fprintferr("c6 = %Z\n",c6);    eps3 = myround(eps3, 1);
     fprintferr("c8 = %Z\n",c8);  
     fprintferr("c11 = %Z\n",c11);    if (DEBUGLEVEL>1) fprintferr("epsilon_3 -> %Z\n",eps3);
     fprintferr("c12 = %Z\n",c12);    *eps5 = mulsr(r, eps3); return A;
     fprintferr("c14 = %Z\n",c14);  
     fprintferr("c15 = %Z\n",c15);  
   }  
 }  }
   
 /* Computation of the logarithmic heights of the units */  /* Performs basic computations concerning the equation.
    * 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 */
 static GEN  static GEN
 Logarithmic_Height(int s)  inithue(GEN P, GEN bnf, long flag, long prec)
 {  {
   int i;    GEN MatFU, x0, tnf, tmp, gpmin, dP, csts, ALH, eps5, ro, c1, c2;
   GEN LH=gun,mat;    int k,j, s,t, n = degpol(P);
   
   mat=gabs(MatFU,Prec);    if (!bnf)
   for (i=1; i<=deg; i++)    {
     LH = gmul(LH,gmax(gun,gabs(gcoeff(mat,i,s),Prec)));      bnf = bnfinit0(P,1,NULL,prec);
   return gmul(gdeux,gdiv(glog(LH,Prec),gdeg));      if (bnf != checkbnf_discard(bnf)) err(talker,"non-monic polynomial in thue");
 }      if (flag) (void)certifybuchall(bnf);
     }
     nf_get_sign(checknf(bnf), (long*)&s, (long*)&t);
     ro = tnf_get_roots(P, prec, s, t);
     MatFU = Conj_LH(gmael(bnf,8,5), &ALH, ro, prec);
     if (!MatFU) return NULL; /* FAIL */
   
 /* Computation of the matrix ((\sigma_i(\eta_j)) used for M_1 and    dP = derivpol(P);
    the computation of logarithmic heights */    c1 = NULL; /* min |P'(r_i)|, i <= s */
 static void    for (k=1; k<=s; k++)
 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");      tmp = gabs(poleval(dP,(GEN)ro[k]),prec);
     for (j=1; j<=deg; j++)      if (!c1 || gcmp(tmp,c1) < 0) c1 = tmp;
       coeff(MatFU,j,i)=(long)poleval((GEN)uf[i],(GEN)roo[j]);  
   }    }
 }    c1 = gdiv(shifti(gun,n-1), c1);
     c1 = gprec_w(myround(c1, 1), DEFAULTPREC);
   
 /* Computation of the conjugates of the solutions to the norm equation */    c2 = NULL; /* max |r_i - r_j|, i!=j */
 static void    for (k=1; k<=n; k++)
 Conj_Norm_Eq(GEN ne, GEN *Hmu)      for (j=k+1; j<=n; j++)
 {      {
   int p,q,nesol,x;        tmp = gabs(gsub((GEN)ro[j],(GEN)ro[k]), prec);
         if (!c2 || gcmp(c2,tmp) > 0) c2 = tmp;
       }
     c2 = gprec_w(myround(c2, -1), DEFAULTPREC);
   
   nesol=lg(ne); MatNE=(GEN)cgetg(nesol,t_MAT); (*Hmu)=cgetg(nesol,t_COL);    if (t==0) x0 = gun;
     else
   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]);      gpmin = NULL; /* min |P'(r_i)|, i > s */
     if (x!=t_INT && x!=t_POL)      for (k=1; k<=t; k++)
       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]);        tmp = gabs(poleval(dP,(GEN)ro[s+k]), prec);
       /* Logarithmic height of the solutions of the norm equation */        if (!gpmin || gcmp(tmp,gpmin) < 0) gpmin = tmp;
       (*Hmu)[q]=lmul((GEN)(*Hmu)[q],gmax(gun,gabs(gcoeff(MatNE,p,q),Prec)));  
     }      }
       gpmin = gprec_w(gpmin, DEFAULTPREC);
   
       /* Compute x0. See paper, Prop. 2.2.1 */
       x0 = gmul(gpmin, Vecmax(gabs(gimag(ro), prec)));
       x0 = mpsqrtn(gdiv(shifti(gun,n-1), x0), n);
   }    }
   for (q=1; q<nesol; q++)    if (DEBUGLEVEL>1) fprintferr("c1 = %Z\nc2 = %Z\n", c1, c2);
     (*Hmu)[q]=ldiv(glog((GEN)(*Hmu)[q],Prec),gdeg);  
     ALH = gmul2n(ALH, 1);
     tnf = cgetg(8,t_VEC); csts = cgetg(7,t_VEC);
     tnf[1] = (long)P;
     tnf[2] = (long)bnf;
     tnf[3] = (long)ro;
     tnf[4] = (long)ALH;
     tnf[5] = (long)MatFU;
     tnf[6] = (long)T_A_Matrices(MatFU, s+t-1, &eps5, prec);
     tnf[7] = (long)csts;
     csts[1] = (long)c1; csts[2] = (long)c2;   csts[3] = (long)LogHeight(ro, prec);
     csts[4] = (long)x0; csts[5] = (long)eps5; csts[6] = (long)stoi(prec);
     return tnf;
 }  }
   
 /* Compute Baker's bound c11, and B_0, the bound for the b_i's .*/  typedef struct {
 static void    GEN c10, c11, c13, c15, bak, NE, ALH, hal, MatFU, ro, Hmu;
 Baker(GEN ALH, GEN Hmu)    GEN delta, lambda, errdelta;
     int r, iroot;
   } baker_s;
   
   /* Compute Baker's bound c9 and B_0, the bound for the b_i's. See Thm 2.3.1 */
   static GEN
   Baker(baker_s *BS)
 {  {
   GEN c9=gun, gbak, hb0=gzero;    const long prec = DEFAULTPREC;
   int k,i1,i2;    GEN tmp, B0, hb0, c9 = gun, ro = BS->ro, ro0 = (GEN)ro[BS->iroot];
     int k, i1, i2, r = BS->r;
   
   gbak = gmul(gmul(gdeg,gsub(gdeg,gun)),gsub(gdeg,gdeux));    switch (BS->iroot) {
       case 1: i1=2; i2=3; break;
   switch (numroot) {      case 2: i1=1; i2=3; break;
   case 1: i1=2; i2=3; break;     default: i1=1; i2=2; 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 */    /* Compute h_1....h_r */
   for (k=1; k<=r; k++)    for (k=1; k<=r; k++)
   {    {
     c9=gmul(c9,gmax((GEN)ALH[k],      tmp = gdiv(gcoeff(BS->MatFU,i1,k), gcoeff(BS->MatFU,i2,k));
                       gmax(ginv(gbak),      tmp = gmax(gun, abslog(tmp,prec));
                            gdiv(gabs(glog(gdiv(gcoeff(MatFU,i1,k),      c9 = gmul(c9, gmax((GEN)BS->ALH[k], gdiv(tmp, BS->bak)));
                                                gcoeff(MatFU,i2,k)),  
                                           Prec),Prec),gbak))));  
   }    }
   
   /* Compute a bound for the h_0 */    /* Compute a bound for the h_0 */
   hb0=gadd(gmul(stoi(4),halpha),gadd(gmul(gdeux,(GEN)Hmu[curne]),    hb0 = gadd(gmul2n(BS->hal,2), gmul(gdeux, gadd(BS->Hmu,mplog2(prec))));
                                      gmul(gdeux,glog(gdeux,Prec))));    tmp = gdiv(gmul(gsub(ro0, (GEN)ro[i2]), (GEN)BS->NE[i1]),
   hb0=gmax(hb0,gmax(ginv(gbak),               gmul(gsub(ro0, (GEN)ro[i1]), (GEN)BS->NE[i2]));
                     gdiv(gabs(glog(gdiv(gmul(gsub((GEN)roo[numroot],    tmp = gmax(gun, abslog(tmp, prec));
                                                   (GEN)roo[i2]),    hb0 = gmax(hb0, gdiv(tmp, BS->bak));
                                              gcoeff(MatNE,i1,curne)),    c9 = gmul(c9,hb0);
                                         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 */    /* Multiply c9 by the "constant" factor */
   c9=gmul(c9,gmul(gmul(gmul(stoi(18),mppi(Prec)),gpui(stoi(32),stoi(4+r),Prec)),    c9 = gmul(c9, gmul(mulri(mulsr(18,mppi(prec)), gpowgs(stoi(32),4+r)),
                   gmul(gmul(mpfact(r+3), gpowgs(gmul(gbak,stoi(r+2)), r+3)),                       gmul(gmul(mpfact(r+3), gpowgs(mulis(BS->bak,r+2), r+3)),
                          glog(gmul(gdeux,gmul(gbak,stoi(r+2))),Prec))));                            glog(mulis(BS->bak,2*(r+2)),prec))));
   c9=myround(c9,gun);    c9 = gprec_w(myround(c9, 1), DEFAULTPREC);
   /* Compute B0 according to Lemma 2.3.3 */    /* Compute B0 according to Lemma 2.3.3 */
   B0=gmax(gexp(gun,Prec),    B0 = mulsr(2, divrr(addrr(mulrr(c9,mplog(divrr(c9,BS->c10))), mplog(BS->c11)),
           mulsr(2,divrr(addrr(mulrr(c9,glog(divrr(c9,c10),ConstPrec)),                        BS->c10));
                               glog(c11,ConstPrec)),c10)));    B0 = gmax(B0, dbltor(2.71828183));
   
   if (DEBUGLEVEL>=2) fprintferr("Baker -> %Z\nB0 -> %Z\n",c9,B0);    if (DEBUGLEVEL>1) {
       fprintferr("  B0  = %Z\n",B0);
       fprintferr("  Baker = %Z\n",c9);
     }
     return B0;
 }  }
   
 /* Compute delta and lambda */  /* || x d ||, x t_REAL, d t_INT */
 static void  static GEN
 Create_CF_Values(int i1, int i2, GEN *errdelta)  errnum(GEN x, GEN d)
 {  {
   GEN eps5;    GEN dx = mulir(d, x);
     return mpabs(subri(dx, ground(dx)));
   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. */  /* Try to reduce the bound through continued fractions; see paper. */
 static int  static int
 CF_First_Pass(GEN kappa, GEN errdelta)  CF_1stPass(GEN *B0, GEN kappa, baker_s *BS)
 {  {
   GEN q,ql,qd,l0;    GEN q, ql, qd, l0, denbound = mulri(*B0, kappa);
   
   if (gcmp(gmul(dbltor(0.1),gsqr(mulri(B0,kappa))),ginv(errdelta))==1)    if (gcmp(gmul(dbltor(0.1),gsqr(denbound)), ginv(BS->errdelta)) > 0) return -1;
   {  
     return -1;  
   }  
   
   q=denom(bestappr(delta, mulir(kappa, B0))); ql=mulir(q,lambda);    q = denom( bestappr(BS->delta, denbound) );
   qd=gmul(q,delta); ql=gabs(subri(ql, ground(ql)),Prec);    qd = errnum(BS->delta, q);
   qd=gabs(mulrr(subri(qd,ground(qd)),B0),Prec);    ql = errnum(BS->lambda,q);
   l0=subrr(ql,addrr(qd,divri(dbltor(0.1),kappa)));  
   if (signe(l0) <= 0)    l0 = subrr(ql, addrr(mulrr(qd, *B0), divri(dbltor(0.1),kappa)));
     if (signe(l0) <= 0) return 0;
   
     if (BS->r > 1)
       *B0 = divrr(mplog(divrr(mulir(q,BS->c15), l0)), BS->c13);
     else
   {    {
     if (DEBUGLEVEL>=2)      l0 = mulrr(l0,Pi2n(1, DEFAULTPREC));
       fprintferr("CF_First_Pass failed. Trying again with larger kappa\n");      *B0 = divrr(mplog(divrr(mulir(q,BS->c11), l0)), BS->c10);
     return 0;  
   }    }
     if (DEBUGLEVEL>1) fprintferr("    B0 -> %Z\n",*B0);
     return 1;
   }
   
     if (r>1)  /* Check whether a solution has already been found */
       B0=divrr(glog(divrr(mulir(q,c15),l0),ConstPrec),c13);  static int
     else  new_sol(GEN z, GEN S)
     B0=divrr(glog(divrr(mulir(q,c11),mulrr(l0,gmul2n(mppi(ConstPrec),1))),ConstPrec),c10);  {
     int i, l = lg(S);
     for (i=1; i<l; i++)
       if (gegal(z,(GEN)S[i])) return 0;
     return 1;
   }
   
     if (DEBUGLEVEL>=2)  static void
     fprintferr("CF_First_Pass successful !!\nB0 -> %Z\n",B0);  add_sol(GEN *pS, GEN x, GEN y)
     return 1;  {
   }    GEN u = cgetg(3,t_VEC); u[1] = (long)x; u[2] = (long)y;
     if (new_sol(u, *pS)) *pS = concatsp(*pS, _vec(u));
   }
   
 /* Check whether a "potential" solution is truly a solution. */  
 static void  static void
 Check_Solutions(GEN xmay1, GEN xmay2, GEN poly, GEN rhs)  check_sol(GEN x, GEN y, GEN P, GEN rhs, GEN *pS)
 {  {
   GEN xx1, xx2, yy1, yy2, zz, u;    if (gcmp0(y))
     { if (gegal(gpowgs(x,degpol(P)), rhs)) add_sol(pS, x, gzero); }
     else
     { if (gegal(poleval(rescale_pol(P,y),x), rhs)) add_sol(pS, x, y); }
   }
   
   yy1=ground(greal(gdiv(gsub(xmay2,xmay1), gsub((GEN)roo[1],(GEN)roo[2]))));  /* Check whether a potential solution is a true solution. Return 0 if
   yy2=gneg(yy1);   * truncation error (increase precision) */
   static int
   CheckSol(GEN *pS, GEN z1, GEN z2, GEN P, GEN rhs, GEN ro)
   {
     GEN x, y, ro1 = (GEN)ro[1], ro2 = (GEN)ro[2];
     long e;
   
   xx1=greal(gadd(xmay1,gmul((GEN)roo[1],yy1)));    y = grndtoi(greal(gdiv(gsub(z2,z1), gsub(ro1,ro2))), &e);
   xx2=gneg(xx1);    if (e > 0) return 0;
     x = gadd(z1, gmul(ro1, y));
   if (gcmp(distoZ(xx1),dbltor(0.0001))==-1)    x = grndtoi(greal(x), &e);
     if (e > 0) return 0;
     if (e <= -13)
   {    {
     xx1=ground(xx1);      check_sol(     x ,      y,  P,rhs,pS);
     if (!gcmp0(yy1))      check_sol(gneg(x), gneg(y), P,rhs,pS);
     {  
       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);  
     }  
   }    }
     return 1;
 }  }
   
   /* find q1,q2,q3 st q1 + b q2 + c q3 ~ 0 */
 static GEN  static GEN
 GuessQi(GEN *q1, GEN *q2, GEN *q3)  GuessQi(GEN b, GEN c, GEN *eps)
 {  {
   GEN C, Lat, eps;    GEN Q, Lat, C = gpowgs(stoi(10), 10);
   
   C=gpui(stoi(10),stoi(10),Prec);    Lat = idmat(3);
   Lat=idmat(3); mael(Lat,1,3)=(long)C; mael(Lat,2,3)=lround(gmul(C,delta));    mael(Lat,1,3) = (long)C;
   mael(Lat,3,3)=lround(gmul(C,lambda));    mael(Lat,2,3) = lround(gmul(C,b));
     mael(Lat,3,3) = lround(gmul(C,c));
   
   Lat=lllint(Lat);    Q = (GEN)lllint(Lat)[1];
   *q1=gmael(Lat,1,1); *q2=gmael(Lat,1,2); *q3=gmael(Lat,1,3);    if (gcmp0((GEN)Q[3])) return NULL; /* FAIL */
   
   eps=gabs(gadd(gadd(*q1,gmul(*q2,delta)),gmul(*q3,lambda)),Prec);    *eps = gadd(gadd((GEN)Q[1], gmul((GEN)Q[2],b)), gmul((GEN)Q[3],c));
   return(eps);    *eps = mpabs(*eps); return Q;
 }  }
   
 static void  
 TotRat(void)  
 {  
   err(bugparier,"thue (totally rational case)");  
 }  
   
 /* Check for solutions under a small bound (see paper) */  /* Check for solutions under a small bound (see paper) */
 static void  static GEN
 Check_Small(int bound, GEN poly, GEN rhs)  SmallSols(GEN S, int Bx, GEN poly, GEN rhs, GEN ro)
 {  {
   GEN interm, xx, zz, u, maxr, tmp, ypot, xxn, xxnm1, yy;    gpmem_t av = avma, lim = stack_lim(av, 1);
   long av = avma, lim = stack_lim(av,1);    const long prec = DEFAULTPREC;
   int x, j, bsupy, y;    GEN Hpoly, interm, X, Xn, Xnm1, Y, sqrtnRHS;
     int x, y, j, By, n = degpol(poly);
   double bndyx;    double bndyx;
   
   maxr=gabs((GEN)roo[1],Prec); for(j=1; j<=deg; j++)    if (DEBUGLEVEL>1) fprintferr("* Checking for small solutions\n");
     { if(gcmp(tmp=gabs((GEN)roo[j],Prec),maxr)==1) maxr=tmp; }  
     sqrtnRHS = absisqrtn(rhs, n, prec);
     bndyx = gtodouble(gadd(sqrtnRHS, Vecmax(gabs(ro,prec))));
   
   bndyx=gtodouble(gadd(gpui(gabs(rhs,Prec),ginv(gdeg),Prec),maxr));    /* x = 0 first */
     Y = ground(sqrtnRHS);
   for (x=-bound; x<=bound; x++)    if (gegal(gpowgs(Y,n), rhs)) add_sol(&S, Y, gzero);
     Y = negi(Y);
     if (gegal(gpowgs(Y,n), rhs)) add_sol(&S, Y, gzero);
     /* x != 0 */
     for (x = -Bx; x <= Bx; x++)
   {    {
     if (low_stack(lim,stack_lim(av,1)))      if (!x) continue;
   
       X = stoi(x); Xn = gpowgs(X,n);
       interm = subii(rhs, mulii(Xn, (GEN)poly[2]));
       Xnm1 = Xn; j = 2;
       while (gcmp0(interm))
     {      {
       if(DEBUGMEM>1) err(warnmem,"Check_small");        Xnm1 = divis(Xnm1, x);
       SOL = gerepilecopy(av, SOL);        interm = mulii((GEN)poly[++j], Xnm1);
     }      }
     if (x==0)      /* y | interm */
       {  
         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 ... */      Hpoly = rescale_pol(poly,X); /* homogenize */
         xxnm1=xxn; j=2;      By = (int)(x>0? bndyx*x: -bndyx*x);
         while(gcmp0(interm))      if (gegal(gmul((GEN)poly[2],Xn),rhs)) add_sol(&S, gzero, X); /* y = 0 */
           {      for(y = -By; y <= By; y++)
             xxnm1=gdiv(xxnm1,xx);      {
             interm=gmul((GEN)poly[++j],xxnm1);        if (!y || smodis(interm, y)) continue;
           }        Y = stoi(y);
         if (gegal(poleval(Hpoly, Y), rhs)) add_sol(&S, Y, X);
       }
   
         for(y=-bsupy; y<=bsupy; y++)      if (low_stack(lim,stack_lim(av,1)))
           {      {
             yy=stoi(y);        if(DEBUGMEM>1) err(warnmem,"Check_small");
             if(y==0) {        S = gerepilecopy(av, S);
               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 (gcmp0(gmod(interm,yy)))  
                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);  
                 }  
           }  
       }  
   }    }
     return S;
 }  }
   
 /* Computes [x]! */  /* Computes [x]! */
 static double  static double
 fact(double x)  fact(double x)
 {  {
   double ft=1.0;    double ft = 1.0;
   x=(int)x; while (x>1) { ft*=x; x--; }    x = (int)x; while (x>1) { ft *= x; x--; }
   return ft ;    return ft ;
 }  }
   
 /* From a polynomial and optionally a system of fundamental units for the  /* From a polynomial and optionally a system of fundamental units for the
  * field defined by poly, computes all the relevant constants needed to solve   * field defined by poly, computes all relevant constants needed to solve
  * the equation P(x,y)=a whenever the solutions of N_{ K/Q }(x)=a.  Returns a   * the equation P(x,y)=a given the solutions of N_{K/Q}(x)=a (see inithue). */
  * "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  GEN
 thueinit(GEN poly, long flag, long prec)  thueinit(GEN poly, long flag, long prec)
 {  {
   GEN thueres,ALH,csts,c0;    GEN tnf, bnf = NULL;
   ulong av = avma;    gpmem_t av = avma;
   long k,st;    long k, s;
   double d,dr;  
   
   uftot = NULL;    if (checktnf(poly)) { bnf = checkbnf((GEN)poly[2]); poly = (GEN)poly[1]; }
   if (checktnf(poly)) { uftot=(GEN)poly[2]; poly=(GEN)poly[1]; }    if (typ(poly)!=t_POL) err(notpoler,"thueinit");
   else    if (degpol(poly) <= 2) err(talker,"invalid polynomial in thue (need deg>2)");
     if (typ(poly)!=t_POL) err(notpoler,"thueinit");  
   if (degpol(poly)<=2) err(talker,"invalid polynomial in thue (need deg>2)");  
   
   if (!gisirreducible(poly)) err(redpoler,"thueinit");    if (!gisirreducible(poly)) err(redpoler,"thueinit");
   st=sturm(poly);    s = sturm(poly);
   if (st)    if (s)
   {    {
     dr=(double)((st+lgef(poly)-5)>>1);      long PREC, n = degpol(poly);
     d=(double)degpol(poly); d=d*(d-1)*(d-2);      double d, dr, dn = (double)n;
     /* Try to guess the precision by approximating Baker's bound.  
      * Note that the guess is most of the time pretty generous,      dr = (double)((s+n-2)>>1); /* s+t-1 */
      * ie 10 to 30 decimal digits above what is *really* necessary.      d = dn*(dn-1)*(dn-2);
      * Note that the limiting step is the reduction. See paper.      /* Guess precision by approximating Baker's bound. The guess is most of
      */       * the time not sharp, ie 10 to 30 decimal digits above what is _really_
     Prec=3 + (long)((5.83 + (dr+4)*5 + log(fact(dr+3)) + (dr+3)*log(dr+2) +       * 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.);                       (dr+3)*log(d) + log(log(2*d*(dr+2))) + (dr+1)) / 10.);
     ConstPrec=4;      if (PREC < prec) PREC = prec;
     if (Prec<prec) Prec = prec;      for (;;)
     if (!checktnf(poly)) inithue(poly,flag);      {
         if (( tnf = inithue(poly, bnf, flag, PREC) )) break;
         PREC = (PREC<<1)-2;
         if (DEBUGLEVEL>1) err(warnprec,"thueinit",PREC);
         bnf = NULL; avma = av;
       }
     }
     else
     {
       GEN c0 = gun, ro = roots(poly, DEFAULTPREC);
       for (k=1; k<lg(ro); k++) c0 = gmul(c0, gimag((GEN)ro[k]));
       c0 = ginv( mpabs(c0) );
       tnf = cgetg(3,t_VEC);
       tnf[1] = (long)poly;
       tnf[2] = (long)c0;
     }
     return gerepilecopy(av,tnf);
   }
   
     thueres=cgetg(8,t_VEC);  static GEN
     thueres[1]=(long)poly;  get_B0(int i1, GEN Delta, GEN Lambda, GEN eps5, long prec, baker_s *BS)
     thueres[2]=(long)uftot;  {
     thueres[3]=(long)roo;    GEN delta, lambda, errdelta, B0 = Baker(BS);
     Compute_Fund_Units(gmael(uftot,8,5));    int i2, r = BS->r;
     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);    i2 = (i1 == 1)? 2: 1;
     csts[1]=(long)c1; csts[2]=(long)c2; csts[3]=(long)halpha;    for(;;) /* i2 from 1 to r unless r = 1 [then i2 = 2] */
     csts[4]=(long)x0; csts[5]=(long)eps3;    {
     csts[6]=(long)stoi(Prec);      if (r > 1)
       {
         delta = divrr((GEN)Delta[i2],(GEN)Delta[i1]);
         lambda = gdiv(gsub(gmul((GEN)Delta[i2],(GEN)Lambda[i1]),
                            gmul((GEN)Delta[i1],(GEN)Lambda[i2])),
                       (GEN)Delta[i1]);
         errdelta = mulrr(addsr(1,delta),
                          divrr(eps5, subrr(mpabs((GEN)Delta[i1]),eps5)));
       }
       else
       { /* r == 1, single fundamental unit (i1 = s = t = 1) */
         GEN p1, Pi2 = Pi2n(1, prec);
         GEN fu = (GEN)BS->MatFU[1], ro = BS->ro;
   
     thueres[7]=(long)csts;        p1 = gdiv((GEN)fu[2], (GEN)fu[3]);
     return gerepilecopy(av,thueres);        delta = divrr(garg(p1,prec), Pi2);
   }  
   
   thueres=cgetg(3,t_VEC); c0=gun; Prec=4;        p1 = gmul(gdiv(gsub((GEN)ro[1], (GEN)ro[2]),
   roo=roots(poly,Prec);                       gsub((GEN)ro[1], (GEN)ro[3])),
   for (k=1; k<lg(roo); k++)                  gdiv((GEN)BS->NE[3], (GEN)BS->NE[2]));
     c0=gmul(c0, gimag((GEN)roo[k]));        lambda = divrr(garg(p1,prec), Pi2);
   c0=ginv(gabs(c0,Prec));  
   thueres[1]=(long)poly; thueres[2]=(long)c0;        errdelta = gdiv(gmul2n(gun, 1 - bit_accuracy(prec)),
   return gerepilecopy(av,thueres);                        gabs((GEN)fu[2],prec));
       }
       BS->delta = delta;
       BS->lambda= lambda; BS->errdelta = errdelta;
       if (DEBUGLEVEL>1) fprintferr("  errdelta = %Z\n",errdelta);
   
       if (DEBUGLEVEL>1) fprintferr("  Entering CF...\n");
       /* Reduce B0 as long as we make progress: newB0 < oldB0 - 0.1 */
       for (;;)
       {
         GEN oldB0 = B0, kappa = stoi(10);
         int cf;
   
         for (cf = 0; cf < 10; cf++, kappa = mulis(kappa,10))
         {
           int res = CF_1stPass(&B0, kappa, BS);
           if (res < 0) return NULL; /* prec problem */
           if (res) break;
           if (DEBUGLEVEL>1) fprintferr("CF failed. Increasing kappa\n");
         }
         if (cf == 10)
         { /* Semirational or totally rational case */
           GEN Q, ep, q, l0, denbound;
   
           if (! (Q = GuessQi(delta, lambda, &ep)) ) break;
   
           denbound = gadd(B0, absi((GEN)Q[2]));
           q = denom( bestappr(delta, denbound) );
           l0 = subrr(errnum(delta, q), ep);
           if (signe(l0) <= 0) break;
   
           B0 = divrr(mplog(divrr(mulir((GEN)Q[3], BS->c15), l0)),  BS->c13);
           if (DEBUGLEVEL>1) fprintferr("Semirat. reduction: B0 -> %Z\n",B0);
         }
         /* if no progress, stop */
         if (gcmp(oldB0, gadd(B0,dbltor(0.1))) <= 0) return gmin(oldB0, B0);
       }
       i2++; if (i2 == i1) i2++;
       if (i2 > r) break;
     }
     err(bugparier,"thue (totally rational case)");
     return NULL; /* not reached */
 }  }
   
 /* Given a tnf structure as returned by thueinit, and a right-hand-side and  static GEN
  * optionally the solutions to the norm equation, returns the solutions to  LargeSols(GEN tnf, GEN rhs, GEN ne, GEN *pro, GEN *pS)
  * 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 Vect, P, ro, bnf, MatFU, A, csts, dP, vecdP;
   GEN xmay1,xmay2,b,zp1,tmp,q1,q2,q3,ep;    GEN c1,c2,c3,c4,c10,c11,c13,c14,c15, x0, x1, x2, x3, b, zp1, tmp, eps5;
   long Nb_It=0,upb,bi1,i1,i2,i, flag,cf,fs;    int iroot, ine, n, i, r,s,t;
   ulong av = avma;    long upb, bi1, Prec, prec;
     baker_s BS;
     gpmem_t av = avma;
   
   if (!checktnf(thueres)) err(talker,"not a tnf in thue");    bnf  = (GEN)tnf[2];
     if (!ne) ne = bnfisintnorm(bnf, rhs);
     if (lg(ne)==1) return NULL;
     nf_get_sign(checknf(bnf), (long*)&s, (long*)&t);
     BS.r = r = s+t-1;
     P      = (GEN)tnf[1]; n = degpol(P);
     ro     = (GEN)tnf[3];
     BS.ALH = (GEN)tnf[4];
     MatFU  = (GEN)tnf[5];
     A      = (GEN)tnf[6];
     csts   = (GEN)tnf[7];
     c1     = (GEN)csts[1]; c1 = gmul(absi(rhs), c1);
     c2     = (GEN)csts[2];
     BS.hal = (GEN)csts[3];
     x0     = (GEN)csts[4];
     eps5   = (GEN)csts[5];
     Prec = gtolong((GEN)csts[6]);
     BS.MatFU = MatFU;
     BS.bak = mulss(n, (n-1)*(n-2)); /* safe */
     *pS = cgetg(1, t_VEC);
   
   SOL = cgetg(1,t_VEC);    if (t) x0 = gmul(x0, absisqrtn(rhs, n, Prec));
   upb = 0; ep = NULL; /* gcc -Wall */    tmp = divrr(c1,c2);
     c3 = mulrr(dbltor(1.39), tmp);
     c4 = mulsr(n-1, c3);
     x1 = gmax(x0, mpsqrtn(mulsr(2,tmp),n));
   
   if (lg(thueres)==8)    Vect = cgetg(r+1,t_COL); for (i=1; i<=r; i++) Vect[i]=un;
     Vect = gmul(gabs(A,DEFAULTPREC), Vect);
     c14 = mulrr(c4, Vecmax(Vect));
     x2 = gmax(x1, mpsqrtn(mulsr(10,c14), n));
     if (DEBUGLEVEL>1) {
       fprintferr("x1 -> %Z\n",x1);
       fprintferr("x2 -> %Z\n",x2);
       fprintferr("c14 = %Z\n",c14);
     }
   
     dP = derivpol(P);
     vecdP = cgetg(s+1, t_VEC);
     for (i=1; i<=s; i++) vecdP[i] = (long)poleval(dP, (GEN)ro[i]);
   
     zp1 = dbltor(0.01);
     x3 = gmax(x2, mpsqrtn(mulsr(2,divrr(c14,zp1)),n));
   
     b = cgetg(r+1,t_COL);
     for (iroot=1; iroot<=s; iroot++)
   {    {
     poly=(GEN)thueres[1]; deg=degpol(poly); gdeg=stoi(deg);      GEN Delta, MatNE, Hmu, c5, c7;
     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);      Vect = cgetg(r+1,t_COL); for (i=1; i<=r; i++) Vect[i] = un;
     if (lg(ne)==1) { avma=av; return cgetg(1,t_VEC); }      if (iroot <= r) Vect[iroot] = lstoi(1-n);
       Delta = gmul(A,Vect);
   
     c1=gmul(gabs(rhs,Prec), (GEN)csts[1]); c2=(GEN)csts[2];      c5 = Vecmax(gabs(Delta,Prec));
     halpha=(GEN)csts[3];      c5  = myround(gprec_w(c5,DEFAULTPREC), 1);
     if (t)      c7  = mulsr(r,c5);
       x0 = gmul((GEN)csts[4],gpui(gabs(rhs,Prec), ginv(gdeg), Prec));      c10 = divsr(n,c7); BS.c10 = c10;
     eps3=(GEN)csts[5];      c13 = divsr(n,c5); BS.c13 = c13;
     Prec=gtolong((GEN)csts[6]);      if (DEBUGLEVEL>1) {
     b=cgetg(r+1,t_COL);        fprintferr("* real root no %ld/%ld\n", iroot,s);
     tmp=divrr(c1,c2);        fprintferr("  c10 = %Z\n",c10);
     x1=gmax(x0,gpui(mulsr(2,tmp),ginv(gdeg),ConstPrec));        fprintferr("  c13 = %Z\n",c13);
     if(DEBUGLEVEL >=2) fprintferr("x1 -> %Z\n",x1);      }
   
     c3=mulrr(dbltor(1.39),tmp);      prec = Prec;
     c4=mulsr(deg-1,c3);      for (;;)
       {
         if (( MatNE = Conj_LH(ne, &Hmu, ro, prec) )) break;
         prec = (prec<<1)-2;
         if (DEBUGLEVEL>1) err(warnprec,"thue",prec);
         ro = tnf_get_roots(P, prec, s, t);
       }
       BS.ro    = ro;
       BS.iroot = iroot;
   
     for (numroot=1; numroot<=s; numroot++)      for (ine=1; ine<lg(ne); ine++)
     {      {
       ComputeConstants();        GEN Lambda, B0, c6, c8;
       Conj_Norm_Eq(ne,&Hmu);        GEN NE = (GEN)MatNE[ine], Vect2 = cgetg(r+1,t_COL);
       for (curne=1; curne<lg(ne); curne++)        int k, i1;
   
         if (DEBUGLEVEL>1) fprintferr("  - norm sol. no %ld/%ld\n",ine,lg(ne)-1);
         for (k=1; k<=r; k++)
       {        {
         ComputeConstants2(poly,rhs);          if (k == iroot)
         Baker(ALH,Hmu);            tmp = gdiv(rhs, gmul((GEN)vecdP[k], (GEN)NE[k]));
           else
         i1=Vecmax(gabs(Delta,Prec),r);            tmp = gdiv(gsub((GEN)ro[iroot],(GEN)ro[k]), (GEN)NE[k]);
         if (i1!=1) i2=1; else i2=2;          Vect2[k] = llog(gabs(tmp,prec), prec);
         do        }
           {        Lambda = gmul(A,Vect2);
             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        c6 = addrr(dbltor(0.1), Vecmax(gabs(Lambda,DEFAULTPREC)));
              * 1) there was less than 10 reductions        c6 = myround(c6, 1);
              * 2) the previous reduction improved the bound of more than 0.1.        c8 = addrr(dbltor(1.23), mulsr(r,c6));
              */        c11= mulrr(mulsr(2,c3) , mpexp(divrr(mulsr(n,c8),c7)));
             while (Nb_It<10 && gcmp(Old_B0,gadd(B0,dbltor(0.1))) == 1 && fs<2)        c15= mulrr(mulsr(2,c14), mpexp(divrr(mulsr(n,c6),c5)));
               {  
                 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 = gerepilecopy(av, 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 (DEBUGLEVEL>1) {
                   {          fprintferr("  c6  = %Z\n",c6);
                     if (!fs)          fprintferr("  c8  = %Z\n",c8);
                       { ep=GuessQi(&q1,&q2,&q3); }          fprintferr("  c11 = %Z\n",c11);
                     bd=gmul(denom(bestappr(delta,gadd(B0,gabs(q2,Prec))))          fprintferr("  c15 = %Z\n",c15);
                             ,delta);        }
                     bd=gabs(gsub(bd,ground(bd)),Prec);        BS.c11 = c11;
                     if(gcmp(bd,ep)==1 && !gcmp0(q3))        BS.c15 = c15;
                       {        BS.NE = NE;
                         fs=1;        BS.Hmu = (GEN)Hmu[ine];
                       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);        i1 = Vecmaxind(gabs(Delta,prec));
         if (! (B0 = get_B0(i1, Delta, Lambda, eps5, prec, &BS)) ) goto PRECPB;
   
       /* For each possible value of b_i1, compute the b_i's       /* For each possible value of b_i1, compute the b_i's
        * and 2 conjugates of x-\alpha y. Then check.        * and 2 conjugates of z = x - alpha y. Then check. */
        */        upb = gtolong(gceil(B0));
         zp1=dbltor(0.01);        for (bi1=-upb; bi1<=upb; bi1++)
         x3=gmax(x2,gpui(mulsr(2,divrr(c14,zp1)),ginv(gdeg),ConstPrec));        {
           GEN z1, z2;
           for (i=1; i<=r; i++)
           {
             b[i] = ldiv(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) > 0) break;
           }
           if (i <= r) continue;
   
         for (bi1=-upb; bi1<=upb; bi1++)          z1 = z2 = gun;
         {          for(i=1; i<=r; i++)
           flag=1;          {
           xmay1=gun; xmay2=gun;            GEN c = ground((GEN)b[i]);
           for (i=1; i<=r; i++)            z1 = gmul(z1, powgi(gcoeff(MatFU,1,i), c));
           {            z2 = gmul(z2, powgi(gcoeff(MatFU,2,i), c));
             b[i]=(long)gdiv(gsub(gmul((GEN)Delta[i],stoi(bi1)),          }
                                   gsub(gmul((GEN)Delta[i],(GEN)Lambda[i1]),          z1 = gmul(z1, (GEN)NE[1]);
                                        gmul((GEN)Delta[i1],(GEN)Lambda[i]))),          z2 = gmul(z2, (GEN)NE[2]);
                              (GEN)Delta[i1]);          if (!CheckSol(pS, z1,z2,P,rhs,ro)) goto PRECPB;
             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);  
     return gerepilecopy(av,SOL);  
   }    }
     *pro = ro; return x3;
   
   /* Case s=0. All the solutions are "small". */  PRECPB:
   Prec=DEFAULTPREC; poly=(GEN)thueres[1]; deg=degpol(poly);    ne = gerepilecopy(av, ne);
   gdeg=stoi(deg); c0=(GEN)thueres[2];    prec += 5 * (DEFAULTPREC-2);
   roo=roots(poly,Prec);    if (DEBUGLEVEL>1) err(warnprec,"thue",prec);
   Check_Small(gtolong(ground(gpui(gmul((GEN)poly[2],gdiv(gabs(rhs,Prec),c0)),    tnf = inithue(P, bnf, 0, prec);
                                   ginv(gdeg),Prec) )), poly, rhs);    return LargeSols(tnf, rhs, ne, pro, pS);
   return gerepilecopy(av,SOL);  
 }  }
   
   /* Given a tnf structure as returned by thueinit, a RHS and
    * optionally the solutions to the norm equation, returns the solutions to
    * the Thue equation F(x,y)=a
    */
   GEN
   thue(GEN tnf, GEN rhs, GEN ne)
   {
     gpmem_t av = avma;
     GEN P, ro, x3, S;
   
     if (!checktnf(tnf)) err(talker,"not a tnf in thue");
     if (typ(rhs) != t_INT) err(typeer,"thue");
   
     P = (GEN)tnf[1];
     if (lg(tnf) == 8)
     {
       x3 = LargeSols(tnf, rhs, ne, &ro, &S);
       if (!x3) { avma = av; return cgetg(1,t_VEC); }
     }
     else
     { /* Case s=0. All solutions are "small". */
       GEN c0   = (GEN)tnf[2]; /* t_REAL */
       S = cgetg(1,t_VEC);
       ro = roots(P, DEFAULTPREC);
       x3 = mpsqrtn(mulir(constant_term(P), divir(absi(rhs),c0)), degpol(P));
     }
     return gerepilecopy(av, SmallSols(S, gtolong(x3), P, rhs, ro));
   }
   
 static GEN *Relations; /* primes above a, expressed on generators of Cl(K) */  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 *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 GEN *gen_ord;   /* orders of generators of Cl(K) given in bnf */
Line 853  test_sol(long i)
Line 811  test_sol(long i)
   
   if (Partial)    if (Partial)
   {    {
     long av=avma;      gpmem_t av=avma;
     for (k=1; k<lg(Partial[1]); k++)      for (k=1; k<lg(Partial[1]); k++)
       if ( signe(modii( (GEN)Partial[i][k], gen_ord[k] )) )        if ( signe(modii( (GEN)Partial[i][k], gen_ord[k] )) )
         { avma=av; return; }          { avma=av; return; }
Line 872  test_sol(long i)
Line 830  test_sol(long i)
   
   for (k=1; k<=i; k++) sol[k]=u[k];    for (k=1; k<=i; k++) sol[k]=u[k];
   for (   ; k<=Nprimes; k++) sol[k]=0;    for (   ; k<=Nprimes; k++) sol[k]=0;
   if (DEBUGLEVEL > 2)    if (DEBUGLEVEL>2)
   {    {
     fprintferr("sol = %Z\n",sol);      fprintferr("sol = %Z\n",sol);
     if (Partial) fprintferr("Partial = %Z\n",Partial);      if (Partial) fprintferr("Partial = %Z\n",Partial);
Line 882  test_sol(long i)
Line 840  test_sol(long i)
 static void  static void
 fix_Partial(long i)  fix_Partial(long i)
 {  {
   long k, av = avma;    long k;
     gpmem_t av = avma;
   for (k=1; k<lg(Partial[1]); k++)    for (k=1; k<lg(Partial[1]); k++)
     addiiz(      addiiz(
       (GEN) Partial[i-1][k],        (GEN) Partial[i-1][k],
Line 1062  bnfisintnorm(GEN bnf, GEN a)
Line 1021  bnfisintnorm(GEN bnf, GEN a)
 {  {
   GEN nf,pol,res,unit,x,id, *Primes;    GEN nf,pol,res,unit,x,id, *Primes;
   long sa,i,j,norm_1;    long sa,i,j,norm_1;
   ulong av = avma;    gpmem_t av = avma;
   
   bnf = checkbnf(bnf); nf = (GEN)bnf[7]; pol = (GEN)nf[1];    bnf = checkbnf(bnf); nf = (GEN)bnf[7]; pol = (GEN)nf[1];
   if (typ(a)!=t_INT)    if (typ(a)!=t_INT)
     err(talker,"expected an integer in bnfisintnorm");      err(talker,"expected an integer in bnfisintnorm");
   sa = signe(a);    sa = signe(a);
   if (sa == 0)  { res=cgetg(2,t_VEC); res[1]=zero; return res; }    if (sa == 0)  return _vec(gzero);
   if (gcmp1(a)) { res=cgetg(2,t_VEC); res[1]=un;   return res; }    if (gcmp1(a)) return _vec(gun);
   
   get_sol_abs(bnf, absi(a), &Primes);    get_sol_abs(bnf, absi(a), &Primes);
   
Line 1077  bnfisintnorm(GEN bnf, GEN a)
Line 1036  bnfisintnorm(GEN bnf, GEN a)
   norm_1 = 0; /* gcc -Wall */    norm_1 = 0; /* gcc -Wall */
   for (i=1; i<=sindex; i++)    for (i=1; i<=sindex; i++)
   {    {
     id = gun; x = normsol[i];      x = normsol[i];
     for (j=1; j<=Nprimes; j++) /* compute prod Primes[i]^u[i] */      id = idealpow(nf,Primes[1], stoi(x[1]));
       if (x[j])      for (j=2; j<=Nprimes; j++) /* compute prod Primes[i]^u[i] */
       {        id = idealmulpowprime(nf,id, Primes[j], stoi(x[j]));
         GEN id2 = Primes[j];  
         if (x[j] != 1) id2 = idealpow(nf,id2, stoi(x[j]));  
         id = idealmul(nf,id,id2);  
       }  
     x = (GEN) isprincipalgenforce(bnf,id)[2];      x = (GEN) isprincipalgenforce(bnf,id)[2];
     x = gmul((GEN)nf[7],x); /* x possible solution */      x = gmul((GEN)nf[7],x); /* x possible solution */
     if (signe(gnorm(gmodulcp(x,(GEN)nf[1]))) != sa)      if (signe(gnorm(gmodulcp(x,(GEN)nf[1]))) != sa)
Line 1093  bnfisintnorm(GEN bnf, GEN a)
Line 1048  bnfisintnorm(GEN bnf, GEN a)
       if (norm_1) x = gmul(unit,x);        if (norm_1) x = gmul(unit,x);
       else        else
       {        {
         if (DEBUGLEVEL > 2)          if (DEBUGLEVEL > 2) fprintferr("%Z eliminated because of sign\n",x);
           fprintferr("%Z eliminated because of sign\n",x);  
         x = NULL;          x = NULL;
       }        }
     }      }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>