Annotation of OpenXM_contrib/pari/src/modules/thue.c, Revision 1.1.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