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