Annotation of OpenXM_contrib/pari-2.2/src/modules/thue.c, Revision 1.1
1.1 ! noro 1: /* $Id: thue.c,v 1.12 2001/10/01 12:11:34 karim Exp $
! 2:
! 3: Copyright (C) 2000 The PARI group.
! 4:
! 5: This file is part of the PARI/GP package.
! 6:
! 7: PARI/GP is free software; you can redistribute it and/or modify it under the
! 8: terms of the GNU General Public License as published by the Free Software
! 9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
! 10: ANY WARRANTY WHATSOEVER.
! 11:
! 12: Check the License for details. You should have received a copy of it, along
! 13: with the package; see the file 'COPYING'. If not, write to the Free Software
! 14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
! 15:
! 16: /* Thue equation solver. In all the forthcoming remarks, "paper"
! 17: * designs the paper "Thue Equations of High Degree", by Yu. Bilu and
! 18: * G. Hanrot, J. Number Theory (1996). Note that numbering of the constants
! 19: * is that of Hanrot's thesis rather than that of the paper.
! 20: * The last part of the program (bnfisintnorm) was written by K. Belabas.
! 21: */
! 22: #include "pari.h"
! 23:
! 24: static int curne,r,s,t,deg,Prec,ConstPrec,numroot;
! 25: static GEN c1,c2,c3,c4,c5,c6,c7,c8,c10,c11,c12,c13,c14,c15,B0,x1,x2,x3,x0;
! 26: static GEN halpha,eps3,gdeg,A,MatNE,roo,MatFU,Delta,Lambda,delta,lambda;
! 27: static GEN Vect2,SOL,uftot;
! 28:
! 29: GEN bnfisintnorm(GEN, GEN);
! 30:
! 31: /* Check whether tnf is a valid structure */
! 32: static int
! 33: checktnf(GEN tnf)
! 34: {
! 35: if (typ(tnf)!=t_VEC || (lg(tnf)!=8 && lg(tnf)!=3)) return 0;
! 36: if (typ(tnf[1])!=t_POL) return 0;
! 37: if (lg(tnf) != 8) return 1; /* s=0 */
! 38:
! 39: deg=degpol(tnf[1]);
! 40: if (deg<=2) err(talker,"invalid polynomial in thue (need deg>2)");
! 41: s=sturm((GEN)tnf[1]); t=(deg-s)>>1; r=s+t-1;
! 42: (void)checkbnf((GEN)tnf[2]);
! 43: if (typ(tnf[3]) != t_COL || lg(tnf[3]) != deg+1) return 0;
! 44: if (typ(tnf[4]) != t_COL || lg(tnf[4]) != r+1) return 0;
! 45: if (typ(tnf[5]) != t_MAT || lg(tnf[5]) != r+1
! 46: || lg(gmael(tnf,5,1)) != deg+1) return 0;
! 47: if (typ(tnf[6]) != t_MAT || lg(tnf[6])!=r+1
! 48: || lg(gmael(tnf,6,1)) != r+1) return 0;
! 49: if (typ(tnf[7]) != t_VEC || lg(tnf[7]) != 7) return 0;
! 50: return 1;
! 51: }
! 52:
! 53: static GEN distoZ(GEN z)
! 54: {
! 55: GEN p1=gfrac(z);
! 56: return gmin(p1, gsub(gun,p1));
! 57: }
! 58:
! 59: /* Check whether a solution has already been found */
! 60: static int
! 61: _thue_new(GEN zz)
! 62: {
! 63: int i;
! 64: for (i=1; i<lg(SOL); i++)
! 65: if (gegal(zz,(GEN)SOL[i])) return 0;
! 66: return 1;
! 67: }
! 68:
! 69: /* Compensates rounding errors for computation/display of the constants */
! 70: static GEN
! 71: myround(GEN Cst, GEN upd)
! 72: {
! 73: return gmul(Cst, gadd(gun, gmul(upd,gpuigs(stoi(10),-10))));
! 74: }
! 75:
! 76: /* Returns the index of the largest element in a vector */
! 77: static int
! 78: Vecmax(GEN Vec, int size)
! 79: {
! 80: int k, tno = 1;
! 81: GEN tmax = (GEN)Vec[1];
! 82: for (k=2; k<=size; k++)
! 83: if (gcmp((GEN)Vec[k],tmax)==1) { tmax=(GEN)Vec[k]; tno=k; }
! 84: return tno;
! 85: }
! 86:
! 87: /* Performs basic computations concerning the equation: buchinitfu, c1, c2 */
! 88: static void
! 89: inithue(GEN poly, long flag)
! 90: {
! 91: GEN roo2, tmp, gpmin, de;
! 92: int k,j;
! 93:
! 94: x0=gzero; x1=gzero; deg=degpol(poly); gdeg=stoi(deg);
! 95:
! 96: if (!uftot)
! 97: {
! 98: uftot=bnfinit0(poly,1,NULL,Prec);
! 99: if (uftot != checkbnf_discard(uftot))
! 100: err(talker,"non-monic polynomial in thue");
! 101: if (flag) certifybuchall(uftot);
! 102: s=itos(gmael3(uftot,7,2,1));
! 103: t=itos(gmael3(uftot,7,2,2));
! 104: }
! 105: /* Switch the roots to get the usual order */
! 106: roo=roots(poly, Prec); roo2=cgetg(deg+1,t_COL);
! 107: for (k=1; k<=s; k++) roo2[k]=roo[k];
! 108: for (k=1; k<=t; k++)
! 109: {
! 110: roo2[k+s]=roo[2*k-1+s];
! 111: roo2[k+s+t]=lconj((GEN)roo2[k+s]);
! 112: }
! 113: roo=roo2;
! 114:
! 115: r=s+t-1; halpha=gun;
! 116: for (k=1; k<=deg; k++)
! 117: halpha = gmul(halpha,gmax(gun,gabs((GEN)roo[k],Prec)));
! 118: halpha=gdiv(glog(halpha,Prec),gdeg);
! 119:
! 120: de=derivpol(poly); c1=gabs(poleval(de,(GEN)roo[1]),Prec);
! 121: for (k=2; k<=s; k++)
! 122: {
! 123: tmp=gabs(poleval(de,(GEN)roo[k]),Prec);
! 124: if (gcmp(tmp,c1)==-1) c1=tmp;
! 125: }
! 126:
! 127: c1=gdiv(gpui(gdeux,gsub(gdeg,gun),Prec),c1); c1=myround(c1,gun);
! 128: c2=gabs(gsub((GEN)roo[1],(GEN)roo[2]),Prec);
! 129:
! 130: for (k=1; k<=deg; k++)
! 131: for (j=k+1; j<=deg; j++)
! 132: {
! 133: tmp=gabs(gsub((GEN)roo[j],(GEN)roo[k]),Prec);
! 134: if (gcmp(c2,tmp)==1) c2=tmp;
! 135: }
! 136:
! 137: c2=myround(c2,stoi(-1));
! 138: if (t==0) x0=gun;
! 139: else
! 140: {
! 141: gpmin=gabs(poleval(de,(GEN)roo[s+1]),Prec);
! 142: for (k=2; k<=t; k++)
! 143: {
! 144: tmp=gabs(poleval(de,(GEN)roo[s+k]),Prec);
! 145: if (gcmp(tmp,gpmin)==-1) gpmin=tmp;
! 146: }
! 147:
! 148: /* Compute x0. See paper, Prop. 2.2.1 */
! 149: x0=gpui(gdiv(gpui(gdeux,gsub(gdeg,gun),Prec),
! 150: gmul(gpmin,
! 151: gabs((GEN)gimag(roo)[Vecmax(gabs(gimag(roo),Prec),deg)],Prec))),
! 152: ginv(gdeg),Prec);
! 153: }
! 154:
! 155: if (DEBUGLEVEL >=2) fprintferr("c1 = %Z\nc2 = %Z\n", c1, c2);
! 156:
! 157: }
! 158:
! 159: /* Computation of M, its inverse A and check of the precision (see paper) */
! 160: static void
! 161: T_A_Matrices(void)
! 162: {
! 163: GEN mask, eps1, eps2, nia, m1, IntM;
! 164: int i,j;
! 165:
! 166: m1=glog(gabs(MatFU,Prec),Prec); mask=gsub(gpui(gdeux,stoi(r),Prec),gun);
! 167: m1=matextract(m1,mask,mask);
! 168:
! 169: A=invmat(m1); IntM=gsub(gmul(A,m1),idmat(r));
! 170:
! 171: eps2=gzero;
! 172: for (i=1; i<=r; i++)
! 173: for (j=1; j<=r; j++)
! 174: if (gcmp(eps2,gabs(gcoeff(IntM,i,j),20)) == -1)
! 175: eps2=gabs(gcoeff(IntM,i,j),20);
! 176:
! 177: eps1=gpui(gdeux,stoi((Prec-2)*32),Prec); eps2=gadd(eps2,ginv(eps1));
! 178:
! 179: nia=gzero;
! 180: for (i=1; i<=r; i++)
! 181: for (j=1; j<=r; j++)
! 182: if (gcmp(nia,gabs(gcoeff(A,i,j),20)) == -1)
! 183: nia = gabs(gcoeff(A,i,j),20);
! 184:
! 185: /* Check for the precision in matrix inversion. See paper, Lemma 2.4.2. */
! 186: if (gcmp(gmul(gadd(gmul(stoi(r),gmul(nia,eps1)),eps2),
! 187: gmul(gdeux,stoi(r))),gun) == -1)
! 188: err(precer,"thue");
! 189:
! 190: eps3=gmul(gdeux,gmul(gmul(gsqr(stoi(r)),nia),
! 191: gadd(gmul(stoi(r),gdiv(nia,eps1)),eps2)));
! 192: myround(eps3,gun);
! 193:
! 194: if (DEBUGLEVEL>=2) fprintferr("epsilon_3 -> %Z\n",eps3);
! 195: }
! 196:
! 197: /* Computation of the constants c5, c7, c10, c12, c15 */
! 198: static void
! 199: ComputeConstants(void)
! 200: {
! 201: int k;
! 202:
! 203: GEN Vect;
! 204:
! 205: Vect=cgetg(r+1,t_COL); for (k=1; k<=r; k++) Vect[k]=un;
! 206: if (numroot <= r) Vect[numroot]=lsub(gun,gdeg);
! 207:
! 208: Delta=gmul(A,Vect);
! 209:
! 210: c5=(GEN)(gabs(Delta,Prec)[Vecmax(gabs(Delta,Prec),r)]);
! 211: c5=myround(c5,gun); c7=mulsr(r,c5);
! 212: c10=divsr(deg,c7); c13=divsr(deg,c5);
! 213:
! 214:
! 215: if (DEBUGLEVEL>=2)
! 216: {
! 217: fprintferr("c5 = %Z\n",c5);
! 218: fprintferr("c7 = %Z\n",c7);
! 219: fprintferr("c10 = %Z\n",c10);
! 220: fprintferr("c13 = %Z\n",c13);
! 221: }
! 222: }
! 223:
! 224: /* Computation of the constants c6, c8, c9 */
! 225: static void
! 226: ComputeConstants2(GEN poly, GEN rhs)
! 227: {
! 228: GEN Vect1, tmp;
! 229: int k;
! 230:
! 231: Vect1=cgetg(r+1,t_COL);
! 232: for (k=1; k<=r; k++) Vect1[k]=un;
! 233: Vect1=gmul(gabs(A,ConstPrec),Vect1);
! 234:
! 235:
! 236: Vect2=cgetg(r+1,t_COL);
! 237: for (k=1; k<=r; k++)
! 238: { if (k!=numroot)
! 239: { Vect2[k]=llog(gabs(gdiv(gsub((GEN)roo[numroot],(GEN)roo[k]),
! 240: gcoeff(MatNE,k,curne)),Prec),Prec); }
! 241: else { Vect2[k]=llog(gabs(gdiv(rhs,gmul(poleval(derivpol(poly)
! 242: ,(GEN)roo[numroot]),
! 243: gcoeff(MatNE,k,curne))),Prec),Prec);
! 244: }
! 245: }
! 246: Lambda=gmul(A,Vect2);
! 247:
! 248: tmp=(GEN)Vect1[Vecmax(Vect1,r)];
! 249: x2=gmax(x1,gpui(mulsr(10,mulrr(c4,tmp)),ginv(gdeg),ConstPrec));
! 250: c14=mulrr(c4,tmp);
! 251:
! 252: c6=gabs((GEN)Lambda[Vecmax(gabs(Lambda,ConstPrec),r)],ConstPrec);
! 253: c6=addrr(c6,dbltor(0.1)); c6=myround(c6,gun);
! 254:
! 255: c8=addrr(dbltor(1.23),mulsr(r,c6));
! 256: c11=mulrr(mulsr(2,c3),gexp(divrr(mulsr(deg,c8),c7),ConstPrec));
! 257:
! 258: tmp=gexp(divrr(mulsr(deg,c6),c5),ConstPrec);
! 259: c12=mulrr(mulsr(2,c3),tmp);
! 260: c15=mulsr(2,mulrr(c14,tmp));
! 261:
! 262: if (DEBUGLEVEL>=2)
! 263: {
! 264: fprintferr("c6 = %Z\n",c6);
! 265: fprintferr("c8 = %Z\n",c8);
! 266: fprintferr("c11 = %Z\n",c11);
! 267: fprintferr("c12 = %Z\n",c12);
! 268: fprintferr("c14 = %Z\n",c14);
! 269: fprintferr("c15 = %Z\n",c15);
! 270: }
! 271: }
! 272:
! 273: /* Computation of the logarithmic heights of the units */
! 274: static GEN
! 275: Logarithmic_Height(int s)
! 276: {
! 277: int i;
! 278: GEN LH=gun,mat;
! 279:
! 280: mat=gabs(MatFU,Prec);
! 281: for (i=1; i<=deg; i++)
! 282: LH = gmul(LH,gmax(gun,gabs(gcoeff(mat,i,s),Prec)));
! 283: return gmul(gdeux,gdiv(glog(LH,Prec),gdeg));
! 284: }
! 285:
! 286: /* Computation of the matrix ((\sigma_i(\eta_j)) used for M_1 and
! 287: the computation of logarithmic heights */
! 288: static void
! 289: Compute_Fund_Units(GEN uf)
! 290: {
! 291: int i,j;
! 292: MatFU=cgetg(r+1,t_MAT);
! 293:
! 294: for (i=1; i<=r; i++)
! 295: MatFU[i]=lgetg(deg+1,t_COL);
! 296: for (i=1; i<=r; i++)
! 297: {
! 298: if (typ(uf[i])!=t_POL) err(talker,"incorrect system of units");
! 299: for (j=1; j<=deg; j++)
! 300: coeff(MatFU,j,i)=(long)poleval((GEN)uf[i],(GEN)roo[j]);
! 301: }
! 302: }
! 303:
! 304: /* Computation of the conjugates of the solutions to the norm equation */
! 305: static void
! 306: Conj_Norm_Eq(GEN ne, GEN *Hmu)
! 307: {
! 308: int p,q,nesol,x;
! 309:
! 310: nesol=lg(ne); MatNE=(GEN)cgetg(nesol,t_MAT); (*Hmu)=cgetg(nesol,t_COL);
! 311:
! 312: for (p=1; p<nesol; p++) { MatNE[p]=lgetg(deg+1,t_COL); (*Hmu)[p]=un; }
! 313: for (q=1; q<nesol; q++)
! 314: {
! 315: x=typ(ne[q]);
! 316: if (x!=t_INT && x!=t_POL)
! 317: err(talker,"incorrect solutions of norm equation");
! 318: for (p=1; p<=deg; p++)
! 319: {
! 320: coeff(MatNE,p,q)=(long)poleval((GEN)ne[q],(GEN)roo[p]);
! 321: /* Logarithmic height of the solutions of the norm equation */
! 322: (*Hmu)[q]=lmul((GEN)(*Hmu)[q],gmax(gun,gabs(gcoeff(MatNE,p,q),Prec)));
! 323: }
! 324: }
! 325: for (q=1; q<nesol; q++)
! 326: (*Hmu)[q]=ldiv(glog((GEN)(*Hmu)[q],Prec),gdeg);
! 327: }
! 328:
! 329: /* Compute Baker's bound c11, and B_0, the bound for the b_i's .*/
! 330: static void
! 331: Baker(GEN ALH, GEN Hmu)
! 332: {
! 333: GEN c9=gun, gbak, hb0=gzero;
! 334: int k,i1,i2;
! 335:
! 336: gbak = gmul(gmul(gdeg,gsub(gdeg,gun)),gsub(gdeg,gdeux));
! 337:
! 338: switch (numroot) {
! 339: case 1: i1=2; i2=3; break;
! 340: case 2: i1=1; i2=3; break;
! 341: case 3: i1=1; i2=2; break;
! 342: default: i1=1; i2=2; break;
! 343: }
! 344:
! 345: /* For the symbols used for the computation of c11, see paper, Thm 2.3.1 */
! 346: /* Compute h_1....h_r */
! 347: for (k=1; k<=r; k++)
! 348: {
! 349: c9=gmul(c9,gmax((GEN)ALH[k],
! 350: gmax(ginv(gbak),
! 351: gdiv(gabs(glog(gdiv(gcoeff(MatFU,i1,k),
! 352: gcoeff(MatFU,i2,k)),
! 353: Prec),Prec),gbak))));
! 354: }
! 355:
! 356: /* Compute a bound for the h_0 */
! 357: hb0=gadd(gmul(stoi(4),halpha),gadd(gmul(gdeux,(GEN)Hmu[curne]),
! 358: gmul(gdeux,glog(gdeux,Prec))));
! 359: hb0=gmax(hb0,gmax(ginv(gbak),
! 360: gdiv(gabs(glog(gdiv(gmul(gsub((GEN)roo[numroot],
! 361: (GEN)roo[i2]),
! 362: gcoeff(MatNE,i1,curne)),
! 363: gmul(gsub((GEN)roo[numroot],
! 364: (GEN)roo[i1]),
! 365: gcoeff(MatNE,i2,curne))),
! 366: Prec),Prec),gbak)));
! 367: c9=gmul(c9,hb0);
! 368: /* Multiply c9 by the "constant" factor */
! 369: c9=gmul(c9,gmul(gmul(gmul(stoi(18),mppi(Prec)),gpui(stoi(32),stoi(4+r),Prec)),
! 370: gmul(gmul(mpfact(r+3), gpowgs(gmul(gbak,stoi(r+2)), r+3)),
! 371: glog(gmul(gdeux,gmul(gbak,stoi(r+2))),Prec))));
! 372: c9=myround(c9,gun);
! 373: /* Compute B0 according to Lemma 2.3.3 */
! 374: B0=gmax(gexp(gun,Prec),
! 375: mulsr(2,divrr(addrr(mulrr(c9,glog(divrr(c9,c10),ConstPrec)),
! 376: glog(c11,ConstPrec)),c10)));
! 377:
! 378: if (DEBUGLEVEL>=2) fprintferr("Baker -> %Z\nB0 -> %Z\n",c9,B0);
! 379: }
! 380:
! 381: /* Compute delta and lambda */
! 382: static void
! 383: Create_CF_Values(int i1, int i2, GEN *errdelta)
! 384: {
! 385: GEN eps5;
! 386:
! 387: if (r>1)
! 388: {
! 389: delta=divrr((GEN)Delta[i2],(GEN)Delta[i1]);
! 390: eps5=mulrs(eps3,r);
! 391: *errdelta=mulrr(addsr(1,delta),
! 392: divrr(eps5,subrr(gabs((GEN)Delta[i1],ConstPrec),eps5)));
! 393:
! 394: lambda=gdiv(gsub(gmul((GEN)Delta[i2],(GEN)Lambda[i1]),
! 395: gmul((GEN)Delta[i1],(GEN)Lambda[i2])),
! 396: (GEN)Delta[i1]);
! 397: }
! 398: else
! 399: {
! 400: GEN Pi2 = gmul2n(mppi(Prec),1);
! 401: delta=gdiv(gcoeff(MatFU,2,1),gcoeff(MatFU,3,1));
! 402: delta=gdiv(garg(delta,Prec),Pi2);
! 403: *errdelta=gdiv(gpui(gdeux,stoi(-(Prec-2)*32+1),Prec),
! 404: gabs(gcoeff(MatFU,2,1),Prec));
! 405: lambda=gmul(gdiv(gsub((GEN)roo[numroot],(GEN)roo[2]),
! 406: gsub((GEN)roo[numroot],(GEN)roo[3])),
! 407: gdiv(gcoeff(MatNE,3,curne),gcoeff(MatNE,2,curne)));
! 408: lambda=gdiv(garg(lambda,Prec),Pi2);
! 409: }
! 410: if (DEBUGLEVEL>=2) outerr(*errdelta);
! 411: }
! 412:
! 413: /* Try to reduce the bound through continued fractions; see paper. */
! 414: static int
! 415: CF_First_Pass(GEN kappa, GEN errdelta)
! 416: {
! 417: GEN q,ql,qd,l0;
! 418:
! 419: if (gcmp(gmul(dbltor(0.1),gsqr(mulri(B0,kappa))),ginv(errdelta))==1)
! 420: {
! 421: return -1;
! 422: }
! 423:
! 424: q=denom(bestappr(delta, mulir(kappa, B0))); ql=mulir(q,lambda);
! 425: qd=gmul(q,delta); ql=gabs(subri(ql, ground(ql)),Prec);
! 426: qd=gabs(mulrr(subri(qd,ground(qd)),B0),Prec);
! 427: l0=subrr(ql,addrr(qd,divri(dbltor(0.1),kappa)));
! 428: if (signe(l0) <= 0)
! 429: {
! 430: if (DEBUGLEVEL>=2)
! 431: fprintferr("CF_First_Pass failed. Trying again with larger kappa\n");
! 432: return 0;
! 433: }
! 434:
! 435: if (r>1)
! 436: B0=divrr(glog(divrr(mulir(q,c15),l0),ConstPrec),c13);
! 437: else
! 438: B0=divrr(glog(divrr(mulir(q,c11),mulrr(l0,gmul2n(mppi(ConstPrec),1))),ConstPrec),c10);
! 439:
! 440: if (DEBUGLEVEL>=2)
! 441: fprintferr("CF_First_Pass successful !!\nB0 -> %Z\n",B0);
! 442: return 1;
! 443: }
! 444:
! 445: /* Check whether a "potential" solution is truly a solution. */
! 446: static void
! 447: Check_Solutions(GEN xmay1, GEN xmay2, GEN poly, GEN rhs)
! 448: {
! 449: GEN xx1, xx2, yy1, yy2, zz, u;
! 450:
! 451: yy1=ground(greal(gdiv(gsub(xmay2,xmay1), gsub((GEN)roo[1],(GEN)roo[2]))));
! 452: yy2=gneg(yy1);
! 453:
! 454: xx1=greal(gadd(xmay1,gmul((GEN)roo[1],yy1)));
! 455: xx2=gneg(xx1);
! 456:
! 457: if (gcmp(distoZ(xx1),dbltor(0.0001))==-1)
! 458: {
! 459: xx1=ground(xx1);
! 460: if (!gcmp0(yy1))
! 461: {
! 462: if (gegal(gmul(poleval(poly,gdiv(xx1,yy1)),
! 463: gpui(yy1,gdeg,Prec)),(GEN)rhs))
! 464: {
! 465: zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
! 466: u[1]=(long)xx1; u[2]=(long)yy1; zz[1]=(long)u;
! 467: if (_thue_new(u)) SOL=concatsp(SOL,zz);
! 468: }
! 469: }
! 470: else if (gegal(gpui(xx1,gdeg,Prec),(GEN)rhs))
! 471: {
! 472: zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
! 473: u[1]=(long)xx1; u[2]=(long)gzero; zz[1]=(long)u;
! 474: if (_thue_new(u)) SOL=concatsp(SOL,zz);
! 475: }
! 476:
! 477: xx2=ground(xx2);
! 478: if (!gcmp0(yy2))
! 479: {
! 480: if (gegal(gmul(poleval(poly,gdiv(xx2,yy2)),
! 481: gpui(yy2,gdeg,Prec)),(GEN)rhs))
! 482: {
! 483: zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
! 484: u[1]=(long)xx2; u[2]=(long)yy2; zz[1]=(long)u;
! 485: if (_thue_new(u)) SOL=concatsp(SOL,zz);
! 486: }
! 487: }
! 488: else if (gegal(gpui(xx2,gdeg,Prec),(GEN)rhs))
! 489: {
! 490: zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
! 491: u[1]=(long)xx2; u[2]=(long)gzero; zz[1]=(long)u;
! 492: if (_thue_new(u)) SOL=concatsp(SOL,zz);
! 493: }
! 494: }
! 495: }
! 496:
! 497: static GEN
! 498: GuessQi(GEN *q1, GEN *q2, GEN *q3)
! 499: {
! 500: GEN C, Lat, eps;
! 501:
! 502: C=gpui(stoi(10),stoi(10),Prec);
! 503: Lat=idmat(3); mael(Lat,1,3)=(long)C; mael(Lat,2,3)=lround(gmul(C,delta));
! 504: mael(Lat,3,3)=lround(gmul(C,lambda));
! 505:
! 506: Lat=lllint(Lat);
! 507: *q1=gmael(Lat,1,1); *q2=gmael(Lat,1,2); *q3=gmael(Lat,1,3);
! 508:
! 509: eps=gabs(gadd(gadd(*q1,gmul(*q2,delta)),gmul(*q3,lambda)),Prec);
! 510: return(eps);
! 511: }
! 512:
! 513: static void
! 514: TotRat(void)
! 515: {
! 516: err(bugparier,"thue (totally rational case)");
! 517: }
! 518:
! 519: /* Check for solutions under a small bound (see paper) */
! 520: static void
! 521: Check_Small(int bound, GEN poly, GEN rhs)
! 522: {
! 523: GEN interm, xx, zz, u, maxr, tmp, ypot, xxn, xxnm1, yy;
! 524: long av = avma, lim = stack_lim(av,1);
! 525: int x, j, bsupy, y;
! 526: double bndyx;
! 527:
! 528: maxr=gabs((GEN)roo[1],Prec); for(j=1; j<=deg; j++)
! 529: { if(gcmp(tmp=gabs((GEN)roo[j],Prec),maxr)==1) maxr=tmp; }
! 530:
! 531: bndyx=gtodouble(gadd(gpui(gabs(rhs,Prec),ginv(gdeg),Prec),maxr));
! 532:
! 533: for (x=-bound; x<=bound; x++)
! 534: {
! 535: if (low_stack(lim,stack_lim(av,1)))
! 536: {
! 537: if(DEBUGMEM>1) err(warnmem,"Check_small");
! 538: SOL = gerepilecopy(av, SOL);
! 539: }
! 540: if (x==0)
! 541: {
! 542: ypot=gmul(stoi(gsigne(rhs)),ground(gpui(gabs(rhs,0),ginv(gdeg),Prec)));
! 543: if (gegal(gpui(ypot,gdeg,0),rhs))
! 544: {
! 545: zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
! 546: u[1]=(long)ypot; u[2]=(long)gzero; zz[1]=(long)u;
! 547: if (_thue_new(u)) SOL=concatsp(SOL,zz);
! 548: }
! 549: if (gegal(gpui(gneg(ypot),gdeg,0),rhs))
! 550: {
! 551: zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
! 552: u[1]=(long)gneg(ypot); u[2]=(long)gzero; zz[1]=(long)u;
! 553: if (_thue_new(u)) SOL=concatsp(SOL,zz);
! 554: }
! 555: }
! 556: else
! 557: {
! 558: bsupy=(int)(x>0?bndyx*x:-bndyx*x);
! 559:
! 560: xx=stoi(x); xxn=gpui(xx,gdeg,Prec);
! 561: interm=gsub((GEN)rhs,gmul(xxn,(GEN)poly[2]));
! 562:
! 563: /* Verifier ... */
! 564: xxnm1=xxn; j=2;
! 565: while(gcmp0(interm))
! 566: {
! 567: xxnm1=gdiv(xxnm1,xx);
! 568: interm=gmul((GEN)poly[++j],xxnm1);
! 569: }
! 570:
! 571: for(y=-bsupy; y<=bsupy; y++)
! 572: {
! 573: yy=stoi(y);
! 574: if(y==0) {
! 575: if (gegal(gmul((GEN)poly[2],xxn),rhs))
! 576: {
! 577: zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
! 578: u[1]=(long)gzero; u[2]=(long)xx; zz[1]=(long)u;
! 579: if (_thue_new(u)) SOL=concatsp(SOL,zz);
! 580: }
! 581: }
! 582: else if (gcmp0(gmod(interm,yy)))
! 583: if(gegal(poleval(poly,gdiv(yy,xx)),gdiv(rhs,xxn)))
! 584: /* Remplacer par un eval *homogene* */
! 585: {
! 586: zz=cgetg(2,t_VEC); u=cgetg(3,t_VEC);
! 587: u[1]=(long)yy; u[2]=(long)xx; zz[1]=(long)u;
! 588: if (_thue_new(u)) SOL=concatsp(SOL,zz);
! 589: }
! 590: }
! 591: }
! 592: }
! 593: }
! 594:
! 595: /* Computes [x]! */
! 596: static double
! 597: fact(double x)
! 598: {
! 599: double ft=1.0;
! 600: x=(int)x; while (x>1) { ft*=x; x--; }
! 601: return ft ;
! 602: }
! 603:
! 604: /* From a polynomial and optionally a system of fundamental units for the
! 605: * field defined by poly, computes all the relevant constants needed to solve
! 606: * the equation P(x,y)=a whenever the solutions of N_{ K/Q }(x)=a. Returns a
! 607: * "tnf" structure containing
! 608: * 1) the polynomial
! 609: * 2) the bnf (used to solve the norm equation)
! 610: * 3) roots, with presumably enough precision
! 611: * 4) The logarithmic heights of units
! 612: * 5) The matrix of conjugates of units
! 613: * 6) its inverse
! 614: * 7) a few technical constants
! 615: */
! 616: GEN
! 617: thueinit(GEN poly, long flag, long prec)
! 618: {
! 619: GEN thueres,ALH,csts,c0;
! 620: ulong av = avma;
! 621: long k,st;
! 622: double d,dr;
! 623:
! 624: uftot = NULL;
! 625: if (checktnf(poly)) { uftot=(GEN)poly[2]; poly=(GEN)poly[1]; }
! 626: else
! 627: if (typ(poly)!=t_POL) err(notpoler,"thueinit");
! 628: if (degpol(poly)<=2) err(talker,"invalid polynomial in thue (need deg>2)");
! 629:
! 630: if (!gisirreducible(poly)) err(redpoler,"thueinit");
! 631: st=sturm(poly);
! 632: if (st)
! 633: {
! 634: dr=(double)((st+lgef(poly)-5)>>1);
! 635: d=(double)degpol(poly); d=d*(d-1)*(d-2);
! 636: /* Try to guess the precision by approximating Baker's bound.
! 637: * Note that the guess is most of the time pretty generous,
! 638: * ie 10 to 30 decimal digits above what is *really* necessary.
! 639: * Note that the limiting step is the reduction. See paper.
! 640: */
! 641: Prec=3 + (long)((5.83 + (dr+4)*5 + log(fact(dr+3)) + (dr+3)*log(dr+2) +
! 642: (dr+3)*log(d) + log(log(2*d*(dr+2))) + (dr+1)) / 10.);
! 643: ConstPrec=4;
! 644: if (Prec<prec) Prec = prec;
! 645: if (!checktnf(poly)) inithue(poly,flag);
! 646:
! 647: thueres=cgetg(8,t_VEC);
! 648: thueres[1]=(long)poly;
! 649: thueres[2]=(long)uftot;
! 650: thueres[3]=(long)roo;
! 651: Compute_Fund_Units(gmael(uftot,8,5));
! 652: ALH=cgetg(r+1,t_COL);
! 653: for (k=1; k<=r; k++) ALH[k]=(long)Logarithmic_Height(k);
! 654: thueres[4]=(long)ALH;
! 655: thueres[5]=(long)MatFU;
! 656: T_A_Matrices();
! 657: thueres[6]=(long)A;
! 658:
! 659: csts=cgetg(7,t_VEC);
! 660: csts[1]=(long)c1; csts[2]=(long)c2; csts[3]=(long)halpha;
! 661: csts[4]=(long)x0; csts[5]=(long)eps3;
! 662: csts[6]=(long)stoi(Prec);
! 663:
! 664: thueres[7]=(long)csts;
! 665: return gerepilecopy(av,thueres);
! 666: }
! 667:
! 668: thueres=cgetg(3,t_VEC); c0=gun; Prec=4;
! 669: roo=roots(poly,Prec);
! 670: for (k=1; k<lg(roo); k++)
! 671: c0=gmul(c0, gimag((GEN)roo[k]));
! 672: c0=ginv(gabs(c0,Prec));
! 673: thueres[1]=(long)poly; thueres[2]=(long)c0;
! 674: return gerepilecopy(av,thueres);
! 675: }
! 676:
! 677: /* Given a tnf structure as returned by thueinit, and a right-hand-side and
! 678: * optionally the solutions to the norm equation, returns the solutions to
! 679: * the Thue equation F(x,y)=a
! 680: */
! 681: GEN
! 682: thue(GEN thueres, GEN rhs, GEN ne)
! 683: {
! 684: GEN Kstart,Old_B0,ALH,errdelta,Hmu,c0,poly,csts,bd;
! 685: GEN xmay1,xmay2,b,zp1,tmp,q1,q2,q3,ep;
! 686: long Nb_It=0,upb,bi1,i1,i2,i, flag,cf,fs;
! 687: ulong av = avma;
! 688:
! 689: if (!checktnf(thueres)) err(talker,"not a tnf in thue");
! 690:
! 691: SOL = cgetg(1,t_VEC);
! 692: upb = 0; ep = NULL; /* gcc -Wall */
! 693:
! 694: if (lg(thueres)==8)
! 695: {
! 696: poly=(GEN)thueres[1]; deg=degpol(poly); gdeg=stoi(deg);
! 697: uftot=(GEN)thueres[2];
! 698: roo=gcopy((GEN)thueres[3]);
! 699: ALH=(GEN)thueres[4];
! 700: MatFU=gcopy((GEN)thueres[5]);
! 701: A=gcopy((GEN)thueres[6]);
! 702: csts=(GEN)thueres[7];
! 703:
! 704: if (!ne) ne = bnfisintnorm(uftot, rhs);
! 705: if (lg(ne)==1) { avma=av; return cgetg(1,t_VEC); }
! 706:
! 707: c1=gmul(gabs(rhs,Prec), (GEN)csts[1]); c2=(GEN)csts[2];
! 708: halpha=(GEN)csts[3];
! 709: if (t)
! 710: x0 = gmul((GEN)csts[4],gpui(gabs(rhs,Prec), ginv(gdeg), Prec));
! 711: eps3=(GEN)csts[5];
! 712: Prec=gtolong((GEN)csts[6]);
! 713: b=cgetg(r+1,t_COL);
! 714: tmp=divrr(c1,c2);
! 715: x1=gmax(x0,gpui(mulsr(2,tmp),ginv(gdeg),ConstPrec));
! 716: if(DEBUGLEVEL >=2) fprintferr("x1 -> %Z\n",x1);
! 717:
! 718: c3=mulrr(dbltor(1.39),tmp);
! 719: c4=mulsr(deg-1,c3);
! 720:
! 721: for (numroot=1; numroot<=s; numroot++)
! 722: {
! 723: ComputeConstants();
! 724: Conj_Norm_Eq(ne,&Hmu);
! 725: for (curne=1; curne<lg(ne); curne++)
! 726: {
! 727: ComputeConstants2(poly,rhs);
! 728: Baker(ALH,Hmu);
! 729:
! 730: i1=Vecmax(gabs(Delta,Prec),r);
! 731: if (i1!=1) i2=1; else i2=2;
! 732: do
! 733: {
! 734: fs=0;
! 735: Create_CF_Values(i1,i2,&errdelta);
! 736: if (DEBUGLEVEL>=2) fprintferr("Entering CF\n");
! 737: Old_B0=gmul(B0,gdeux);
! 738:
! 739: /* Try to reduce B0 while
! 740: * 1) there was less than 10 reductions
! 741: * 2) the previous reduction improved the bound of more than 0.1.
! 742: */
! 743: while (Nb_It<10 && gcmp(Old_B0,gadd(B0,dbltor(0.1))) == 1 && fs<2)
! 744: {
! 745: cf=0; Old_B0=B0; Nb_It++; Kstart=stoi(10);
! 746: while (!fs && CF_First_Pass(Kstart,errdelta) == 0 && cf < 8 )
! 747: {
! 748: cf++;
! 749: Kstart=mulis(Kstart,10);
! 750: }
! 751: if ( CF_First_Pass(Kstart,errdelta) == -1 )
! 752: { ne = gerepilecopy(av, ne); Prec+=10;
! 753: if(DEBUGLEVEL>=2) fprintferr("Increasing precision\n");
! 754: thueres=thueinit(thueres,0,Prec);
! 755: return(thue(thueres,rhs,ne)); }
! 756:
! 757: if (cf==8 || fs) /* Semirational or totally rational case */
! 758: {
! 759: if (!fs)
! 760: { ep=GuessQi(&q1,&q2,&q3); }
! 761: bd=gmul(denom(bestappr(delta,gadd(B0,gabs(q2,Prec))))
! 762: ,delta);
! 763: bd=gabs(gsub(bd,ground(bd)),Prec);
! 764: if(gcmp(bd,ep)==1 && !gcmp0(q3))
! 765: {
! 766: fs=1;
! 767: B0=gdiv(glog(gdiv(gmul(q3,c15),gsub(bd,ep)), Prec),c13);
! 768: if (DEBUGLEVEL>=2)
! 769: fprintferr("Semirat. reduction ok. B0 -> %Z\n",B0);
! 770: }
! 771: else fs=2;
! 772: }
! 773: else fs=0;
! 774: Nb_It=0; B0=gmin(Old_B0,B0); upb=gtolong(gceil(B0));
! 775: }
! 776: i2++; if (i2==i1) i2++;
! 777: }
! 778: while (fs == 2 && i2 <= r);
! 779:
! 780: if (fs == 2) TotRat();
! 781:
! 782: if (DEBUGLEVEL >=2) fprintferr("x2 -> %Z\n",x2);
! 783:
! 784: /* For each possible value of b_i1, compute the b_i's
! 785: * and 2 conjugates of x-\alpha y. Then check.
! 786: */
! 787: zp1=dbltor(0.01);
! 788: x3=gmax(x2,gpui(mulsr(2,divrr(c14,zp1)),ginv(gdeg),ConstPrec));
! 789:
! 790: for (bi1=-upb; bi1<=upb; bi1++)
! 791: {
! 792: flag=1;
! 793: xmay1=gun; xmay2=gun;
! 794: for (i=1; i<=r; i++)
! 795: {
! 796: b[i]=(long)gdiv(gsub(gmul((GEN)Delta[i],stoi(bi1)),
! 797: gsub(gmul((GEN)Delta[i],(GEN)Lambda[i1]),
! 798: gmul((GEN)Delta[i1],(GEN)Lambda[i]))),
! 799: (GEN)Delta[i1]);
! 800: if (gcmp(distoZ((GEN)b[i]),zp1)==1) { flag=0; break; }
! 801: }
! 802:
! 803: if(flag)
! 804: {
! 805: for(i=1; i<=r; i++)
! 806: {
! 807: b[i]=lround((GEN)b[i]);
! 808: xmay1=gmul(xmay1,gpui(gcoeff(MatFU,1,i),(GEN)b[i],Prec));
! 809: xmay2=gmul(xmay2,gpui(gcoeff(MatFU,2,i),(GEN)b[i],Prec));
! 810: }
! 811: xmay1=gmul(xmay1,gcoeff(MatNE,1,curne));
! 812: xmay2=gmul(xmay2,gcoeff(MatNE,2,curne));
! 813: Check_Solutions(xmay1,xmay2,poly,rhs);
! 814: }
! 815: }
! 816: }
! 817: }
! 818: if (DEBUGLEVEL>=2) fprintferr("Checking for small solutions\n");
! 819: Check_Small((int)(gtodouble(x3)),poly,rhs);
! 820: return gerepilecopy(av,SOL);
! 821: }
! 822:
! 823: /* Case s=0. All the solutions are "small". */
! 824: Prec=DEFAULTPREC; poly=(GEN)thueres[1]; deg=degpol(poly);
! 825: gdeg=stoi(deg); c0=(GEN)thueres[2];
! 826: roo=roots(poly,Prec);
! 827: Check_Small(gtolong(ground(gpui(gmul((GEN)poly[2],gdiv(gabs(rhs,Prec),c0)),
! 828: ginv(gdeg),Prec) )), poly, rhs);
! 829: return gerepilecopy(av,SOL);
! 830: }
! 831:
! 832: static GEN *Relations; /* primes above a, expressed on generators of Cl(K) */
! 833: static GEN *Partial; /* list of vvectors, Partial[i] = Relations[1..i] * u[1..i] */
! 834: static GEN *gen_ord; /* orders of generators of Cl(K) given in bnf */
! 835:
! 836: static long *f; /* f[i] = f(Primes[i]/p), inertial degree */
! 837: static long *n; /* a = prod p^{ n_p }. n[i]=n_p if Primes[i] divides p */
! 838: static long *inext; /* index of first P above next p, 0 if p is last */
! 839: static long *S; /* S[i] = n[i] - sum_{ 1<=k<=i } f[k].u[k] */
! 840: static long *u; /* We want principal ideals I = prod Primes[i]^u[i] */
! 841: static long **normsol; /* lists of copies of the u[] which are solutions */
! 842:
! 843: static long Nprimes; /* length(Relations) = #{max ideal above divisors of a} */
! 844: static long sindex, smax; /* current index in normsol; max. index */
! 845:
! 846: /* u[1..i] has been filled. Norm(u) is correct.
! 847: * Check relations in class group then save it.
! 848: */
! 849: static void
! 850: test_sol(long i)
! 851: {
! 852: long k,*sol;
! 853:
! 854: if (Partial)
! 855: {
! 856: long av=avma;
! 857: for (k=1; k<lg(Partial[1]); k++)
! 858: if ( signe(modii( (GEN)Partial[i][k], gen_ord[k] )) )
! 859: { avma=av; return; }
! 860: avma=av;
! 861: }
! 862: if (sindex == smax)
! 863: {
! 864: long new_smax = smax << 1;
! 865: long **new_normsol = (long **) new_chunk(new_smax+1);
! 866:
! 867: for (k=1; k<=smax; k++) new_normsol[k] = normsol[k];
! 868: normsol = new_normsol; smax = new_smax;
! 869: }
! 870: sol = cgetg(Nprimes+1,t_VECSMALL);
! 871: normsol[++sindex] = sol;
! 872:
! 873: for (k=1; k<=i; k++) sol[k]=u[k];
! 874: for ( ; k<=Nprimes; k++) sol[k]=0;
! 875: if (DEBUGLEVEL > 2)
! 876: {
! 877: fprintferr("sol = %Z\n",sol);
! 878: if (Partial) fprintferr("Partial = %Z\n",Partial);
! 879: flusherr();
! 880: }
! 881: }
! 882: static void
! 883: fix_Partial(long i)
! 884: {
! 885: long k, av = avma;
! 886: for (k=1; k<lg(Partial[1]); k++)
! 887: addiiz(
! 888: (GEN) Partial[i-1][k],
! 889: mulis((GEN) Relations[i][k], u[i]),
! 890: (GEN) Partial[i][k]
! 891: );
! 892: avma=av;
! 893: }
! 894:
! 895: /* Recursive loop. Suppose u[1..i] has been filled
! 896: * Find possible solutions u such that, Norm(prod Prime[i]^u[i]) = a, taking
! 897: * into account:
! 898: * 1) the relations in the class group if need be.
! 899: * 2) the factorization of a.
! 900: */
! 901: static void
! 902: isintnorm_loop(long i)
! 903: {
! 904: if (S[i] == 0) /* sum u[i].f[i] = n[i], do another prime */
! 905: {
! 906: int k;
! 907: if (inext[i] == 0) { test_sol(i); return; }
! 908:
! 909: /* there are some primes left */
! 910: if (Partial) gaffect(Partial[i], Partial[inext[i]-1]);
! 911: for (k=i+1; k < inext[i]; k++) u[k]=0;
! 912: i=inext[i]-1;
! 913: }
! 914: else if (i == inext[i]-2 || i == Nprimes-1)
! 915: {
! 916: /* only one Prime left above prime; change prime, fix u[i+1] */
! 917: if (S[i] % f[i+1]) return;
! 918: i++; u[i] = S[i-1] / f[i];
! 919: if (Partial) fix_Partial(i);
! 920: if (inext[i]==0) { test_sol(i); return; }
! 921: }
! 922:
! 923: i++; u[i]=0;
! 924: if (Partial) gaffect(Partial[i-1], Partial[i]);
! 925: if (i == inext[i-1])
! 926: { /* change prime */
! 927: if (inext[i] == i+1 || i == Nprimes) /* only one Prime above p */
! 928: {
! 929: S[i]=0; u[i] = n[i] / f[i]; /* we already know this is exact */
! 930: if (Partial) fix_Partial(i);
! 931: }
! 932: else S[i] = n[i];
! 933: }
! 934: else S[i] = S[i-1]; /* same prime, different Prime */
! 935: for(;;)
! 936: {
! 937: isintnorm_loop(i); S[i]-=f[i]; if (S[i]<0) break;
! 938: if (Partial)
! 939: gaddz(Partial[i], Relations[i], Partial[i]);
! 940: u[i]++;
! 941: }
! 942: }
! 943:
! 944: static void
! 945: get_sol_abs(GEN bnf, GEN a, GEN **ptPrimes)
! 946: {
! 947: GEN dec, fact, primes, *Primes, *Fact;
! 948: long *gcdlist, gcd,nprimes,Ngen,i,j;
! 949:
! 950: if (gcmp1(a))
! 951: {
! 952: long *sol = new_chunk(Nprimes+1);
! 953: sindex = 1; normsol = (long**) new_chunk(2);
! 954: normsol[1] = sol; for (i=1; i<=Nprimes; i++) sol[i]=0;
! 955: return;
! 956: }
! 957:
! 958: fact=factor(a); primes=(GEN)fact[1];
! 959: nprimes=lg(primes)-1; sindex = 0;
! 960: gen_ord = (GEN *) mael3(bnf,8,1,2); Ngen = lg(gen_ord)-1;
! 961:
! 962: Fact = (GEN*) new_chunk(1+nprimes);
! 963: gcdlist = new_chunk(1+nprimes);
! 964:
! 965: Nprimes=0;
! 966: for (i=1; i<=nprimes; i++)
! 967: {
! 968: long ldec;
! 969:
! 970: dec = primedec(bnf,(GEN)primes[i]); ldec = lg(dec)-1;
! 971: gcd = itos(gmael(dec,1,4));
! 972:
! 973: /* check that gcd_{P|p} f_P | n_p */
! 974: for (j=2; gcd>1 && j<=ldec; j++)
! 975: gcd = cgcd(gcd,itos(gmael(dec,j,4)));
! 976:
! 977: gcdlist[i]=gcd;
! 978:
! 979: if (gcd != 1 && smodis(gmael(fact,2,i),gcd))
! 980: {
! 981: if (DEBUGLEVEL>2)
! 982: { fprintferr("gcd f_P does not divide n_p\n"); flusherr(); }
! 983: return;
! 984: }
! 985: Fact[i] = dec; Nprimes += ldec;
! 986: }
! 987:
! 988: f = new_chunk(1+Nprimes); u = new_chunk(1+Nprimes);
! 989: n = new_chunk(1+Nprimes); S = new_chunk(1+Nprimes);
! 990: inext = new_chunk(1+Nprimes);
! 991: Primes = (GEN*) new_chunk(1+Nprimes);
! 992: *ptPrimes = Primes;
! 993:
! 994: if (Ngen)
! 995: {
! 996: Partial = (GEN*) new_chunk(1+Nprimes);
! 997: Relations = (GEN*) new_chunk(1+Nprimes);
! 998: }
! 999: else /* trivial class group, no relations to check */
! 1000: Relations = Partial = NULL;
! 1001:
! 1002: j=0;
! 1003: for (i=1; i<=nprimes; i++)
! 1004: {
! 1005: long k,lim,v, vn=itos(gmael(fact,2,i));
! 1006:
! 1007: gcd = gcdlist[i];
! 1008: dec = Fact[i]; lim = lg(dec);
! 1009: v = (i==nprimes)? 0: j + lim;
! 1010: for (k=1; k < lim; k++)
! 1011: {
! 1012: j++; Primes[j] = (GEN)dec[k];
! 1013: f[j] = itos(gmael(dec,k,4)) / gcd;
! 1014: n[j] = vn / gcd; inext[j] = v;
! 1015: if (Partial)
! 1016: Relations[j] = isprincipal(bnf,Primes[j]);
! 1017: }
! 1018: }
! 1019: if (Partial)
! 1020: {
! 1021: for (i=1; i <= Nprimes; i++)
! 1022: if (!gcmp0(Relations[i])) break;
! 1023: if (i > Nprimes) Partial = NULL; /* all ideals dividing a are principal */
! 1024: }
! 1025: if (Partial)
! 1026: for (i=0; i<=Nprimes; i++) /* Partial[0] needs to be initialized */
! 1027: {
! 1028: Partial[i]=cgetg(Ngen+1,t_COL);
! 1029: for (j=1; j<=Ngen; j++)
! 1030: {
! 1031: Partial[i][j]=lgeti(4);
! 1032: gaffect(gzero,(GEN) Partial[i][j]);
! 1033: }
! 1034: }
! 1035: smax=511; normsol = (long**) new_chunk(smax+1);
! 1036: S[0]=n[1]; inext[0]=1; isintnorm_loop(0);
! 1037: }
! 1038:
! 1039: static long
! 1040: get_unit_1(GEN bnf, GEN pol, GEN *unit)
! 1041: {
! 1042: GEN fu;
! 1043: long j;
! 1044:
! 1045: if (DEBUGLEVEL > 2)
! 1046: fprintferr("looking for a fundamental unit of norm -1\n");
! 1047:
! 1048: *unit = gmael3(bnf,8,4,2); /* torsion unit */
! 1049: if (signe( gnorm(gmodulcp(*unit,pol)) ) < 0) return 1;
! 1050:
! 1051: fu = gmael(bnf,8,5); /* fundamental units */
! 1052: for (j=1; j<lg(fu); j++)
! 1053: {
! 1054: *unit = (GEN)fu[j];
! 1055: if (signe( gnorm(gmodulcp(*unit,pol)) ) < 0) return 1;
! 1056: }
! 1057: return 0;
! 1058: }
! 1059:
! 1060: GEN
! 1061: bnfisintnorm(GEN bnf, GEN a)
! 1062: {
! 1063: GEN nf,pol,res,unit,x,id, *Primes;
! 1064: long sa,i,j,norm_1;
! 1065: ulong av = avma;
! 1066:
! 1067: bnf = checkbnf(bnf); nf = (GEN)bnf[7]; pol = (GEN)nf[1];
! 1068: if (typ(a)!=t_INT)
! 1069: err(talker,"expected an integer in bnfisintnorm");
! 1070: sa = signe(a);
! 1071: if (sa == 0) { res=cgetg(2,t_VEC); res[1]=zero; return res; }
! 1072: if (gcmp1(a)) { res=cgetg(2,t_VEC); res[1]=un; return res; }
! 1073:
! 1074: get_sol_abs(bnf, absi(a), &Primes);
! 1075:
! 1076: res = cgetg(1,t_VEC); unit = NULL;
! 1077: norm_1 = 0; /* gcc -Wall */
! 1078: for (i=1; i<=sindex; i++)
! 1079: {
! 1080: id = gun; x = normsol[i];
! 1081: for (j=1; j<=Nprimes; j++) /* compute prod Primes[i]^u[i] */
! 1082: if (x[j])
! 1083: {
! 1084: GEN id2 = Primes[j];
! 1085: if (x[j] != 1) id2 = idealpow(nf,id2, stoi(x[j]));
! 1086: id = idealmul(nf,id,id2);
! 1087: }
! 1088: x = (GEN) isprincipalgenforce(bnf,id)[2];
! 1089: x = gmul((GEN)nf[7],x); /* x possible solution */
! 1090: if (signe(gnorm(gmodulcp(x,(GEN)nf[1]))) != sa)
! 1091: {
! 1092: if (! unit) norm_1 = get_unit_1(bnf,pol,&unit);
! 1093: if (norm_1) x = gmul(unit,x);
! 1094: else
! 1095: {
! 1096: if (DEBUGLEVEL > 2)
! 1097: fprintferr("%Z eliminated because of sign\n",x);
! 1098: x = NULL;
! 1099: }
! 1100: }
! 1101: if (x) res = concatsp(res, gmod(x,pol));
! 1102: }
! 1103: return gerepilecopy(av,res);
! 1104: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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-2.2 / src / modules