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