Annotation of OpenXM_contrib/pari-2.2/src/basemath/bibli1.c, Revision 1.1.1.1
1.1 noro 1: /* $Id: bibli1.c,v 1.76 2001/10/01 12:11:29 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: /********************************************************************/
17: /** **/
18: /** LLL Algorithm and close friends **/
19: /** **/
20: /********************************************************************/
21: #include "pari.h"
22: #include "parinf.h"
23: extern GEN ZV_lincomb(GEN u, GEN v, GEN X, GEN Y);
24: extern GEN roots_to_pol_r1r2(GEN a, long r1, long v);
25: extern GEN makebasis(GEN nf,GEN pol);
26: extern GEN caractducos(GEN p, GEN x, int v);
27:
28: /* scalar product x.x */
29: GEN
30: sqscal(GEN x)
31: {
32: long i,av,lx;
33: GEN z;
34: lx = lg(x);
35: if (lx == 1) return gzero;
36: av = avma;
37: z = gsqr((GEN)x[1]);
38: for (i=2; i<lx; i++)
39: z = gadd(z, gsqr((GEN)x[i]));
40: return gerepileupto(av,z);
41: }
42:
43: /* scalar product x.y */
44: GEN
45: gscal(GEN x,GEN y)
46: {
47: long i,av,lx;
48: GEN z;
49: if (x == y) return sqscal(x);
50: lx = lg(x);
51: if (lx == 1) return gzero;
52: av = avma;
53: z = gmul((GEN)x[1],(GEN)y[1]);
54: for (i=2; i<lx; i++)
55: z = gadd(z, gmul((GEN)x[i],(GEN)y[i]));
56: return gerepileupto(av,z);
57: }
58:
59: static GEN
60: sqscali(GEN x)
61: {
62: long i,av,lx;
63: GEN z;
64: lx = lg(x);
65: if (lx == 1) return gzero;
66: av = avma;
67: z = sqri((GEN)x[1]);
68: for (i=2; i<lx; i++)
69: z = addii(z, sqri((GEN)x[i]));
70: return gerepileuptoint(av,z);
71: }
72:
73: static GEN
74: gscali(GEN x,GEN y)
75: {
76: long i,av,lx;
77: GEN z;
78: if (x == y) return sqscali(x);
79: lx = lg(x);
80: if (lx == 1) return gzero;
81: av = avma;
82: z = mulii((GEN)x[1],(GEN)y[1]);
83: for (i=2; i<lx; i++)
84: z = addii(z, mulii((GEN)x[i],(GEN)y[i]));
85: return gerepileuptoint(av,z);
86: }
87:
88: static GEN
89: lllall_trivial(GEN x, long n, long flag)
90: {
91: GEN y;
92: if (!n)
93: {
94: if (flag != lll_ALL) return cgetg(1,t_MAT);
95: y=cgetg(3,t_VEC);
96: y[1]=lgetg(1,t_MAT);
97: y[2]=lgetg(1,t_MAT); return y;
98: }
99: /* here n = 1 */
100: if (gcmp0((GEN)x[1]))
101: {
102: switch(flag ^ lll_GRAM)
103: {
104: case lll_KER: return idmat(1);
105: case lll_IM : return cgetg(1,t_MAT);
106: default: y=cgetg(3,t_VEC);
107: y[1]=(long)idmat(1);
108: y[2]=lgetg(1,t_MAT); return y;
109: }
110: }
111: if (flag & lll_GRAM) flag ^= lll_GRAM; else x = NULL;
112: switch (flag)
113: {
114: case lll_KER: return cgetg(1,t_MAT);
115: case lll_IM : return idmat(1);
116: default: y=cgetg(3,t_VEC);
117: y[1]=lgetg(1,t_MAT);
118: y[2]=x? lcopy(x): (long)idmat(1); return y;
119: }
120: }
121:
122: static GEN
123: lllgramall_finish(GEN h,GEN fl,long flag,long n)
124: {
125: long k;
126: GEN y;
127:
128: k=1; while (k<=n && !fl[k]) k++;
129: switch(flag)
130: {
131: case lll_KER: setlg(h,k);
132: y = gcopy(h); break;
133:
134: case lll_IM: h += k-1; h[0] = evaltyp(t_MAT)| evallg(n-k+2);
135: y = gcopy(h); break;
136:
137: default: setlg(h,k); y=cgetg(3,t_VEC);
138: y[1] = lcopy(h);
139: h += k-1; h[0] = evaltyp(t_MAT)| evallg(n-k+2);
140: y[2] = lcopy(h);
141: break;
142: }
143: return y;
144: }
145:
146: /********************************************************************/
147: /** **/
148: /** LLL with content **/
149: /** **/
150: /********************************************************************/
151:
152: /* real Gram matrix has coeffs X[i,j] = x[i,j]*veccon[i]*veccon[j] */
153: static GEN
154: lllgramintwithcontent(GEN x, GEN veccon, long flag)
155: {
156: long av0=avma,av,tetpil,lx=lg(x),i,j,k,l,n,lim,kmax;
157: GEN u,u2,B,lam,q,r,h,la,bb,p1,p2,p3,p4,fl,corr,corr2,newcon;
158: GEN *gptr[8];
159:
160: if (typ(x) != t_MAT) err(typeer,"lllgramintwithcontent");
161: n=lx-1; if (n<=1) return lllall_trivial(x,n,flag);
162: if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramintwithcontent");
163: fl = new_chunk(lx);
164:
165: av=avma; lim=stack_lim(av,1);
166: x=dummycopy(x); veccon=dummycopy(veccon);
167: B=cgetg(lx+1,t_COL); B[1]=un;
168: /* B[i+1]=B_i; le vrai B_i est B_i*prod(1,j=1,i,veccon[j]^2) */
169:
170: for (i=1; i<lx; i++) { B[i+1]=zero; fl[i]=0; }
171: lam=cgetg(lx,t_MAT);
172: for (j=1; j<lx; j++)
173: { p2=cgetg(lx,t_COL); lam[j]=(long)p2; for (i=1; i<lx; i++) p2[i]=zero; }
174: /* le vrai lam[i,j] est
175: lam[i,j]*veccon[i]*veccon[j]*prod(1,l=1,j-1,veccon[l]^2) */
176: k=2; h=idmat(n); kmax=1;
177: u=gcoeff(x,1,1); if (typ(u)!=t_INT) err(lllger4);
178: if (signe(u)) { B[2]=(long)u; coeff(lam,1,1)=un; fl[1]=1; }
179: else { B[2]=un; fl[1]=0; }
180: if (DEBUGLEVEL>5) { fprintferr("k = %ld",k); flusherr(); }
181: for(;;)
182: {
183: if (k>kmax)
184: {
185: kmax=k;
186: for (j=1; j<=k; j++)
187: {
188: if (j==k || fl[j])
189: {
190: u=gcoeff(x,k,j); if (typ(u)!=t_INT) err(lllger4);
191: for (i=1; i<j; i++)
192: if (fl[i])
193: u=divii(subii(mulii((GEN)B[i+1],u),mulii(gcoeff(lam,k,i),gcoeff(lam,j,i))),(GEN)B[i]);
194: if (j<k) coeff(lam,k,j)=(long)u;
195: else
196: {
197: if (signe(u)) { B[k+1]=(long)u; coeff(lam,k,k)=un; fl[k]=1; }
198: else { B[k+1]=B[k]; fl[k]=0; }
199: }
200: }
201: }
202: if (low_stack(lim, stack_lim(av,1)))
203: {
204: if(DEBUGMEM>1) err(warnmem,"[1]: lllgramintwithcontent");
205: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
206: gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
207: }
208: }
209: u=shifti(mulii(gcoeff(lam,k,k-1),(GEN)veccon[k]),1);
210: u2=mulii((GEN)B[k],(GEN)veccon[k-1]);
211: if (cmpii(absi(u),u2)>0)
212: {
213: q=dvmdii(addii(u,u2),shifti(u2,1),&r);
214: if (signe(r)<0) q=addsi(-1,q);
215: h[k]=lsub((GEN)h[k],gmul(q,(GEN)h[k-1]));
216: newcon=mppgcd((GEN)veccon[k],(GEN)veccon[k-1]);
217: corr=divii((GEN)veccon[k],newcon); veccon[k]=(long)newcon;
218: if(!gcmp1(corr))
219: {
220: corr2=sqri(corr);
221: for (j=1; j<=n; j++)
222: coeff(x,j,k)=coeff(x,k,j)=lmulii(corr,gcoeff(x,k,j));
223: coeff(x,k,k)=lmulii(corr,gcoeff(x,k,k));
224: for (j=k; j<=kmax; j++) B[j+1]=lmulii(corr2,(GEN)B[j+1]);
225: for (i=1; i<k; i++) coeff(lam,k,i)=lmulii(corr,gcoeff(lam,k,i));
226: for (i=k+1; i<=kmax; i++)
227: {
228: coeff(lam,i,k)=lmulii(corr,gcoeff(lam,i,k));
229: for (j=k+1; j<i; j++)
230: coeff(lam,i,j)=lmulii(corr2,gcoeff(lam,i,j));
231: }
232: }
233: r=negi(mulii(q,divii((GEN)veccon[k-1],(GEN)veccon[k])));
234: p1=gcoeff(x,k,k-1);
235: for (j=1; j<=n; j++)
236: coeff(x,j,k)=coeff(x,k,j)=laddii(gcoeff(x,j,k),mulii(r,(j==k)?p1:gcoeff(x,j,k-1)));
237: coeff(x,k,k)=laddii(gcoeff(x,k,k),mulii(r,gcoeff(x,k-1,k)));
238: coeff(lam,k,k-1)=laddii(gcoeff(lam,k,k-1),mulii(r,(GEN)B[k]));
239: for (i=1; i<k-1; i++)
240: coeff(lam,k,i)=laddii(gcoeff(lam,k,i),mulii(r,gcoeff(lam,k-1,i)));
241: }
242: if (low_stack(lim, stack_lim(av,1)))
243: {
244: if(DEBUGMEM>1) err(warnmem,"[2]: lllgramintwithcontent");
245: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
246: gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
247: }
248: p3=mulii((GEN)B[k-1],(GEN)B[k+1]);la=gcoeff(lam,k,k-1);p4=mulii(la,la);
249: p2=mulsi(100,mulii(mulii((GEN)veccon[k],(GEN)veccon[k]),addii(p3,p4)));
250: p1=mulii((GEN)veccon[k-1],(GEN)B[k]);p1=mulsi(99,mulii(p1,p1));
251: if (fl[k-1] && (cmpii(p1,p2)>0 || !fl[k]))
252: {
253: if (DEBUGLEVEL>=4 && k==n)
254: { fprintferr(" (%ld)", expi(p1)-expi(p2)); flusherr(); }
255: p1=(GEN)h[k-1]; h[k-1]=h[k]; h[k]=(long)p1;
256: p1=(GEN)x[k-1]; x[k-1]=x[k]; x[k]=(long)p1;
257: for (j=1; j<=n; j++)
258: { p1=gcoeff(x,k-1,j); coeff(x,k-1,j)=coeff(x,k,j); coeff(x,k,j)=(long)p1; }
259: p1=(GEN)veccon[k-1]; veccon[k-1]=veccon[k]; veccon[k]=(long)p1;
260: for (j=1; j<=k-2; j++)
261: { p1=gcoeff(lam,k-1,j); coeff(lam,k-1,j)=coeff(lam,k,j); coeff(lam,k,j)=(long)p1; }
262: if (fl[k])
263: {
264: for (i=k+1; i<=kmax; i++)
265: {
266: bb=gcoeff(lam,i,k);
267: coeff(lam,i,k)=ldivii(subii(mulii((GEN)B[k+1],gcoeff(lam,i,k-1)),mulii(la,bb)),(GEN)B[k]);
268: coeff(lam,i,k-1)=ldivii(addii(mulii(la,gcoeff(lam,i,k-1)),mulii((GEN)B[k-1],bb)),(GEN)B[k]);
269: if (low_stack(lim, stack_lim(av,1)))
270: {
271: if(DEBUGMEM>1) err(warnmem,"[3]: lllgramintwithcontent");
272: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
273: gptr[3]=&x; gptr[4]=&la; gptr[5]=&veccon; gptr[6]=&p3;
274: gptr[7]=&p4; gerepilemany(av,gptr,8);
275: }
276: }
277: B[k]=ldivii(addii(p3,p4),(GEN)B[k]);
278: }
279: else
280: {
281: if (signe(la))
282: {
283: p2=(GEN)B[k]; p1=divii(p4,p2);
284: for (i=k+1; i<=kmax; i++)
285: coeff(lam,i,k-1)=ldivii(mulii(la,gcoeff(lam,i,k-1)),p2);
286: for (j=k+1; j<kmax; j++)
287: {
288: for (i=j+1; i<=kmax; i++)
289: coeff(lam,i,j)=ldivii(mulii(p1,gcoeff(lam,i,j)),p2);
290: if (low_stack(lim, stack_lim(av,1)))
291: {
292: if(DEBUGMEM>1) err(warnmem,"[4]: lllgramintwithcontent");
293: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
294: gptr[3]=&x; gptr[4]=&la; gptr[5]=&veccon; gptr[6]=&p1;
295: gptr[7]=&p2; gerepilemany(av,gptr,8);
296: }
297: }
298: B[k+1]=B[k]=(long)p1;
299: for (i=k+2; i<=lx; i++)
300: B[i]=ldivii(mulii(p1,(GEN)B[i]),p2);
301: }
302: else
303: {
304: coeff(lam,k,k-1)=zero;
305: for (i=k+1; i<=kmax; i++)
306: {
307: coeff(lam,i,k)=coeff(lam,i,k-1);
308: coeff(lam,i,k-1)=zero;
309: }
310: B[k]=B[k-1]; fl[k]=1; fl[k-1]=0;
311: }
312:
313: if (low_stack(lim, stack_lim(av,1)))
314: {
315: if(DEBUGMEM>1) err(warnmem,"[5]: lllgramintwithcontent");
316: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
317: gptr[3]=&x; gptr[4]=&veccon;
318: gerepilemany(av,gptr,5);
319: }
320: }
321: if (k>2) k--;
322: if (DEBUGLEVEL>5) { fprintferr(" %ld",k); flusherr(); }
323: }
324: else
325: {
326: for (l=k-2; l>=1; l--)
327: {
328: u=shifti(mulii(gcoeff(lam,k,l),(GEN)veccon[k]),1);
329: u2=mulii((GEN)B[l+1],(GEN)veccon[l]);
330: if (cmpii(absi(u),u2)>0)
331: {
332: q=dvmdii(addii(u,u2),shifti(u2,1),&r);
333: if (signe(r)<0) q=addsi(-1,q);
334: h[k]=lsub((GEN)h[k],gmul(q,(GEN)h[l]));
335: newcon=mppgcd((GEN)veccon[k],(GEN)veccon[l]);
336: corr=divii((GEN)veccon[k],newcon); veccon[k]=(long)newcon;
337: if(!gcmp1(corr))
338: {
339: corr2=sqri(corr);
340: for (j=1; j<=n; j++)
341: coeff(x,j,k)=coeff(x,k,j)=lmulii(corr,gcoeff(x,k,j));
342: coeff(x,k,k)=lmulii(corr,gcoeff(x,k,k));
343: for (j=k; j<=kmax; j++) B[j+1]=lmulii(corr2,(GEN)B[j+1]);
344: for (i=1; i<k; i++) coeff(lam,k,i)=lmulii(corr,gcoeff(lam,k,i));
345: for (i=k+1; i<=kmax; i++)
346: {
347: coeff(lam,i,k)=lmulii(corr,gcoeff(lam,i,k));
348: for (j=k+1; j<i; j++)
349: coeff(lam,i,j)=lmulii(corr2,gcoeff(lam,i,j));
350: }
351: }
352: r=negi(mulii(q,divii((GEN)veccon[l],(GEN)veccon[k])));
353: p1=gcoeff(x,k,l);
354: for (j=1; j<=n; j++)
355: coeff(x,j,k)=coeff(x,k,j)=laddii(gcoeff(x,j,k),mulii(r,(j==k)?p1:gcoeff(x,j,l)));
356: coeff(x,k,k)=laddii(gcoeff(x,k,k),mulii(r,gcoeff(x,l,k)));
357: coeff(lam,k,l)=laddii(gcoeff(lam,k,l),mulii(r,(GEN)B[l+1]));
358: for (i=1; i<l; i++)
359: coeff(lam,k,i)=laddii(gcoeff(lam,k,i),mulii(r,gcoeff(lam,l,i)));
360: }
361: if (low_stack(lim, stack_lim(av,1)))
362: {
363: if(DEBUGMEM>1) err(warnmem,"[6]: lllgramintwithcontent");
364: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
365: gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
366: }
367: }
368: k++; if (DEBUGLEVEL>5) { fprintferr(" %ld",k); flusherr(); }
369: if (k>n)
370: {
371: if (DEBUGLEVEL>5) { fprintferr("\n"); flusherr(); }
372: tetpil=avma;
373: return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
374: }
375: }
376: if (low_stack(lim, stack_lim(av,1)))
377: {
378: if(DEBUGMEM>1) err(warnmem,"[7]: lllgramintwithcontent");
379: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
380: gptr[3]=&x; gptr[4]=&veccon; gerepilemany(av,gptr,5);
381: }
382: }
383: }
384:
385: static GEN
386: lllintwithcontent(GEN x)
387: {
388: long lx=lg(x),i,j,av,tetpil;
389: GEN veccon,con,xred,g;
390:
391: if (typ(x) != t_MAT) err(typeer,"lllintwithcontent");
392: if (lx==1) return cgetg(1,t_MAT);
393: av=avma; veccon=cgetg(lx,t_VEC); g=cgetg(lx,t_MAT); xred=cgetg(lx,t_MAT);
394: for (j=1; j<lx; j++)
395: {
396: g[j]=lgetg(lx,t_COL); con=content((GEN)x[j]);
397: xred[j]=ldiv((GEN)x[j],con); veccon[j]=(long)con;
398: }
399: for (i=1; i<lx; i++)
400: for (j=1; j<=i; j++)
401: coeff(g,i,j)=coeff(g,j,i)=(long)gscal((GEN)xred[i],(GEN)xred[j]);
402: tetpil=avma;
403: return gerepile(av,tetpil,lllgramintwithcontent(g,veccon,2));
404: }
405:
406: /********************************************************************/
407: /** **/
408: /** LLL ALGORITHM **/
409: /** **/
410: /********************************************************************/
411: #define swap(x,y) { long _t=x; x=y; y=_t; }
412: #define gswap(x,y) { GEN _t=x; x=y; y=_t; }
413:
414: static void
415: REDI(GEN x, GEN h, GEN L, GEN B, long K, long k, long l)
416: {
417: long i,lx;
418: GEN q = truedvmdii(addii(B,shifti(gcoeff(L,k,l),1)), shifti(B,1), NULL);
419: GEN *hk,*hl,xk,xl;
420: if (!signe(q)) return;
421: q = negi(q); lx = lg(x);
422: xl = (GEN)x[l]; hl = (GEN*)h[l];
423: xk = (GEN)x[k]; hk = (GEN*)h[k];
424: if (is_pm1(q))
425: {
426: if (signe(q) > 0)
427: {
428: for (i=1;i<=K;i++) hk[i]=addii(hk[i],hl[i]);
429: xk[k] = laddii((GEN)xk[k], (GEN)xl[k]);
430: for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=laddii((GEN)xk[i], (GEN)xl[i]);
431: for (i=1;i<l; i++) coeff(L,k,i)=laddii(gcoeff(L,k,i),gcoeff(L,l,i));
432: q = B;
433: } else {
434: for (i=1;i<=K;i++) hk[i]=subii(hk[i],hl[i]);
435: xk[k] = lsubii((GEN)xk[k], (GEN)xl[k]);
436: for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=lsubii((GEN)xk[i], (GEN)xl[i]);
437: for (i=1;i<l; i++) coeff(L,k,i)=lsubii(gcoeff(L,k,i),gcoeff(L,l,i));
438: q = negi(B);
439: }
440: } else {
441: for(i=1;i<=K;i++) hk[i]=addii(hk[i],mulii(q,hl[i]));
442: xk[k] = laddii((GEN)xk[k], mulii(q,(GEN)xl[k]));
443: for(i=1;i<lx;i++) coeff(x,k,i)=xk[i]=laddii((GEN)xk[i],mulii(q,(GEN)xl[i]));
444: for(i=1;i<l;i++) coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(q,gcoeff(L,l,i)));
445: q = mulii(q,B);
446: }
447: coeff(L,k,l) = laddii(gcoeff(L,k,l), q);
448: }
449:
450: /* b[k] <-- b[k] - round(L[k,l]) b[l], only b[1] ... b[K] modified so far
451: * assume l < k and update x = Gram(b), L = Gram Schmidt coeffs. */
452: static void
453: RED(GEN x, GEN h, GEN L, long K, long k, long l)
454: {
455: long e,i,lx;
456: GEN q = grndtoi(gcoeff(L,k,l),&e);
457: GEN *hk,*hl,xk,xl;
458: if (DEBUGLEVEL>8)
459: fprintferr("error bits when rounding in lllgram: %ld\n",e);
460: if (!signe(q)) return;
461: q = negi(q); lx = lg(x);
462: xl = (GEN)x[l]; hl = (GEN*)h[l];
463: xk = (GEN)x[k]; hk = (GEN*)h[k];
464: if (is_pm1(q))
465: {
466: if (signe(q) > 0)
467: {
468: for (i=1;i<=K;i++) hk[i]=addii(hk[i],hl[i]);
469: xk[k] = ladd((GEN)xk[k], (GEN)xl[k]);
470: for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=ladd((GEN)xk[i], (GEN)xl[i]);
471: for (i=1;i<l; i++) coeff(L,k,i)=ladd(gcoeff(L,k,i),gcoeff(L,l,i));
472: } else {
473: for (i=1;i<=K;i++) hk[i]=subii(hk[i],hl[i]);
474: xk[k] = lsub((GEN)xk[k], (GEN)xl[k]);
475: for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=lsub((GEN)xk[i], (GEN)xl[i]);
476: for (i=1;i<l; i++) coeff(L,k,i)=lsub(gcoeff(L,k,i),gcoeff(L,l,i));
477: }
478: } else {
479: for (i=1;i<=K;i++) hk[i]=addii(hk[i],mulii(q,hl[i]));
480: xk[k] = ladd((GEN)xk[k], gmul(q,(GEN)xl[k]));
481: for (i=1;i<lx;i++) coeff(x,k,i)=xk[i]=ladd((GEN)xk[i], gmul(q,(GEN)xl[i]));
482: for (i=1;i<l; i++) coeff(L,k,i)=ladd(gcoeff(L,k,i),gmul(q,gcoeff(L,l,i)));
483: }
484: coeff(L,k,l) = ladd(gcoeff(L,k,l),q);
485: }
486:
487: static int
488: do_SWAPI(GEN x, GEN h, GEN L, GEN B, long kmax, long k, long alpha, GEN fl)
489: {
490: GEN la,la2,p1,p2,Bk;
491: long av,i,j,lx;
492:
493: if (!fl[k-1]) return 0;
494: lx = lg(x); av = avma;
495: la = gcoeff(L,k,k-1); la2 = sqri(la);
496: p2 = addii(la2, mulii((GEN)B[k-1],(GEN)B[k+1]));
497: Bk = (GEN)B[k];
498: if (fl[k] && cmpii(mulsi(alpha-1,sqri(Bk)), mulsi(alpha,p2)) <= 0)
499: { avma = av; return 0; }
500:
501: /* SWAPI(k-1,k) */
502: if (DEBUGLEVEL>3 && k==kmax) {
503: fprintferr(" (%ld)", expi(mulsi(alpha-1,sqri(Bk)))
504: - expi(mulsi(alpha,p2))); flusherr();
505: }
506: swap(h[k-1], h[k]);
507: swap(x[k-1], x[k]);
508: for (j=1; j< lx; j++) swap(coeff(x,k-1,j), coeff(x,k,j));
509: for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j))
510: if (fl[k])
511: {
512: av = avma;
513: for (i=k+1; i<=kmax; i++)
514: {
515: GEN t = gcoeff(L,i,k);
516: p1 = subii(mulii((GEN)B[k+1],gcoeff(L,i,k-1)), mulii(la,t));
517: p1 = divii(p1, Bk);
518: av = avma = coeff(L,i,k) = (long)icopy_av(p1,(GEN)av);
519:
520: p1 = addii(mulii(la,gcoeff(L,i,k-1)), mulii((GEN)B[k-1],t));
521: p1 = divii(p1, Bk);
522: av = avma = coeff(L,i,k-1) = (long)icopy_av(p1,(GEN)av);
523: }
524: B[k] = ldivii(p2,Bk);
525: }
526: else
527: {
528: if (signe(la))
529: {
530: p1 = divii(la2, Bk);
531: B[k+1] = B[k] = (long)p1;
532: for (i=k+2; i<=lx; i++) B[i] = ldivii(mulii(p1,(GEN)B[i]), Bk);
533: for (i=k+1; i<=kmax; i++)
534: coeff(L,i,k-1) = ldivii(mulii(la,gcoeff(L,i,k-1)), Bk);
535: for (j=k+1; j<kmax; j++)
536: for (i=j+1; i<=kmax; i++)
537: coeff(L,i,j) = ldivii(mulii(p1,gcoeff(L,i,j)), Bk);
538: }
539: else
540: {
541: for (i=k+1; i<=kmax; i++)
542: {
543: coeff(L,i,k) = coeff(L,i,k-1);
544: coeff(L,i,k-1) = zero;
545: }
546: B[k] = B[k-1]; fl[k] = 1; fl[k-1] = 0;
547: }
548: }
549: return 1;
550: }
551:
552: static int
553: do_SWAP(GEN x, GEN h, GEN L, GEN B, long kmax, long k, GEN QR)
554: {
555: GEN la,la2, BK,BB,q;
556: long av,i,j,lx;
557:
558: lx = lg(x); av = avma;
559: la = gcoeff(L,k,k-1); la2 = gsqr(la);
560: q = gmul((GEN)B[k-1], gsub(QR,la2));
561: if (gcmp(q, (GEN)B[k]) <= 0) { avma = av; return 0; }
562:
563: /* SWAP(k-1,k) */
564: if (DEBUGLEVEL>3 && k==kmax) {
565: fprintferr(" (%ld)", gexpo(q) - gexpo((GEN)B[k])); flusherr();
566: }
567: BB = gadd((GEN)B[k], gmul((GEN)B[k-1],la2));
568: if (gcmp0(BB)) { B[k] = 0; return 1; } /* precision pb */
569:
570: coeff(L,k,k-1) = ldiv(gmul(la,(GEN)B[k-1]), BB);
571: BK = gdiv((GEN)B[k], BB);
572: B[k] = lmul((GEN)B[k-1], BK);
573: B[k-1] = (long)BB;
574:
575: swap(h[k-1],h[k]);
576: swap(x[k-1],x[k]);
577: for (j=1; j<lx ; j++) swap(coeff(x,k-1,j), coeff(x,k,j));
578: for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j))
579: for (i=k+1; i<=kmax; i++)
580: {
581: GEN t = gcoeff(L,i,k);
582: coeff(L,i,k) = lsub(gcoeff(L,i,k-1), gmul(la,t));
583: coeff(L,i,k-1) = ladd(gmul(gcoeff(L,k,k-1),gcoeff(L,i,k-1)), gmul(BK,t));
584: /* = ladd(t, gmul(gcoeff(L,k,k-1), gcoeff(L,i,k)));
585: * would save one multiplication, but loses precision faster... */
586: }
587: return 1;
588: }
589:
590: /* x integer matrix */
591: GEN
592: lllgramall(GEN x, long alpha, long flag)
593: {
594: long av0=avma,av,tetpil,lim,lx=lg(x),i,j,k,l,n,s,kmax;
595: GEN u,B,L,h,fl, *gptr[4];
596:
597: if (typ(x) != t_MAT) err(typeer,"lllgramall");
598: n=lx-1; if (n<=1) return lllall_trivial(x,n,flag);
599: if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramall");
600: fl = cgetg(lx,t_VECSMALL);
601:
602: av=avma; lim=stack_lim(av,1); x=dummycopy(x);
603: B=gscalcol(gun, lx);
604: L=cgetg(lx,t_MAT);
605: for (j=1; j<lx; j++)
606: {
607: for (i=1; i<lx; i++)
608: if (typ(gcoeff(x,i,j))!=t_INT) err(lllger4);
609: fl[j] = 0; L[j] = (long)zerocol(n);
610: }
611: k=2; h=idmat(n); kmax=1;
612: u=gcoeff(x,1,1); s= signe(u);
613: if (s == 0) B[2]=un;
614: else
615: {
616: if (s < 0) err(lllger3);
617: B[2]=(long)u; coeff(L,1,1)=un; fl[1]=1;
618: }
619: if (DEBUGLEVEL>5) fprintferr("k =");
620: for(;;)
621: {
622: if (k > kmax)
623: {
624: if (DEBUGLEVEL>3) {fprintferr(" K%ld",k);flusherr();}
625: kmax = k;
626: for (j=1; j<=k; j++)
627: if (j==k || fl[j])
628: {
629: long av1 = avma;
630: u = gcoeff(x,k,j);
631: for (i=1; i<j; i++) if (fl[i])
632: u = divii(subii(mulii((GEN)B[i+1],u),
633: mulii(gcoeff(L,k,i),gcoeff(L,j,i))),
634: (GEN)B[i]);
635: u = gerepileuptoint(av1, u);
636: if (j<k) coeff(L,k,j)=(long)u;
637: else
638: {
639: s = signe(u);
640: if (s < 0) err(lllger3);
641: if (s)
642: { B[k+1]=(long)u; coeff(L,k,k)=un; fl[k]=1; }
643: else
644: { B[k+1]=B[k]; fl[k]=0; }
645: }
646: }
647: } else if (DEBUGLEVEL>5) {fprintferr(" %ld",k); flusherr();}
648: REDI(x,h,L,(GEN)B[k],kmax,k,k-1);
649: if (do_SWAPI(x,h,L,B,kmax,k,alpha,fl))
650: {
651: if (k>2) k--;
652: }
653: else
654: {
655: for (l=k-2; l; l--)
656: {
657: REDI(x,h,L,(GEN)B[l+1],kmax,k,l);
658: if (low_stack(lim, stack_lim(av,1)))
659: {
660: if(DEBUGMEM>1) err(warnmem,"lllgramall[1]");
661: gptr[0]=&B; gptr[1]=&L; gptr[2]=&h; gptr[3]=&x;
662: gerepilemany(av,gptr,4);
663: }
664: }
665: if (++k > n) break;
666: }
667: if (low_stack(lim, stack_lim(av,1)))
668: {
669: if(DEBUGMEM>1) err(warnmem,"lllgramall[2]");
670: gptr[0]=&B; gptr[1]=&L; gptr[2]=&h; gptr[3]=&x;
671: gerepilemany(av,gptr,4);
672: }
673: }
674: if (DEBUGLEVEL>3) fprintferr("\n");
675: tetpil=avma; return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
676: }
677:
678: static GEN
679: lllgramall0(GEN x, long flag) { return lllgramall(x,100,flag); }
680:
681: static int
682: pslg(GEN x)
683: {
684: long tx;
685: if (gcmp0(x)) return 2;
686: tx=typ(x); return is_scalar_t(tx)? 3: lgef(x);
687: }
688:
689: static GEN
690: lllgramallgen(GEN x, long flag)
691: {
692: long av0=avma,av,tetpil,lx=lg(x),tu,i,j,k,l,n,lim;
693: GEN u,B,lam,q,cq,h,la,bb,p1,p2,p3,p4,fl;
694: int ps1,ps2,flc;
695:
696: if (typ(x) != t_MAT) err(typeer,"lllgramallgen");
697: n=lx-1; if (n<=1) return lllall_trivial(x,n,flag);
698: if (lg((GEN)x[1])!=lx) err(mattype1,"lllgramallgen");
699:
700: fl = new_chunk(lx);
701:
702: av=avma; lim=stack_lim(av,1);
703: B=cgetg(lx+1,t_COL);
704: B[1]=un; lam=cgetg(lx,t_MAT);
705: for (j=1; j<lx; j++) lam[j]=lgetg(lx,t_COL);
706: for (i=1; i<lx; i++)
707: for (j=1; j<=i; j++)
708: {
709: if (j<i && !fl[j]) coeff(lam,i,j)=coeff(lam,j,i)=zero;
710: else
711: {
712: u=gcoeff(x,i,j); tu=typ(u);
713: if (! is_scalar_t(tu) && tu != t_POL) err(lllger4);
714: for (k=1; k<j; k++)
715: if (fl[k])
716: u=gdiv(gsub( gmul((GEN)B[k+1],u),
717: gmul(gcoeff(lam,i,k),gcoeff(lam,j,k)) ),
718: (GEN)B[k]);
719: if (j<i) { coeff(lam,i,j)=(long)u; coeff(lam,j,i)=zero; }
720: else
721: {
722: if (!gcmp0(u)) { B[i+1]=(long)u; coeff(lam,i,i)=un; fl[i]=1; }
723: else { B[i+1]=B[i]; coeff(lam,i,i)=zero; fl[i]=0; }
724: }
725: }
726: }
727: k=2; h=idmat(n); flc=0;
728: for(;;)
729: {
730: u=gcoeff(lam,k,k-1);
731: if (pslg(u) >= pslg((GEN)B[k]))
732: {
733: q=gdeuc(u,(GEN)B[k]); cq=gdivsg(1,content(q)); q=gmul(q,cq); flc=1;
734: h[k]=lsub(gmul(cq,(GEN)h[k]),gmul(q,(GEN)h[k-1]));
735: coeff(lam,k,k-1)=lsub(gmul(cq,gcoeff(lam,k,k-1)),gmul(q,(GEN)B[k]));
736: for (i=1; i<k-1; i++)
737: coeff(lam,k,i)=lsub(gmul(cq,gcoeff(lam,k,i)),gmul(q,gcoeff(lam,k-1,i)));
738: }
739: ps1 = pslg(gsqr((GEN)B[k]));
740: p3 = gmul((GEN)B[k-1],(GEN)B[k+1]);
741: la=gcoeff(lam,k,k-1); p4 = gmul(la,gcoeff(lam,k,k-1));
742: ps2=pslg(gadd(p3,p4));
743: if (fl[k-1] && (ps1>ps2 || (ps1==ps2 && flc) || !fl[k]))
744: {
745: flc = (ps1!=ps2);
746: swap(h[k-1],h[k]);
747: for (j=1; j<=k-2; j++) swap(coeff(lam,k-1,j), coeff(lam,k,j));
748: if (fl[k])
749: {
750: for (i=k+1; i<=n; i++)
751: {
752: bb=gcoeff(lam,i,k);
753: coeff(lam,i,k)=ldiv(gsub(gmul((GEN)B[k+1],gcoeff(lam,i,k-1)),gmul(la,bb)),(GEN)B[k]);
754: coeff(lam,i,k-1)=ldiv(gadd(gmul(la,gcoeff(lam,i,k-1)),gmul((GEN)B[k-1],bb)),(GEN)B[k]);
755: }
756: B[k]=ldiv(gadd(p3,p4),(GEN)B[k]);
757: }
758: else
759: {
760: if (!gcmp0(la))
761: {
762: p2=(GEN)B[k]; p1=gdiv(p4,p2);
763: for (i=k+1; i<lx; i++)
764: coeff(lam,i,k-1)=ldiv(gmul(la,gcoeff(lam,i,k-1)),p2);
765: for (j=k+1; j<lx-1; j++)
766: for (i=j+1; i<lx; i++)
767: coeff(lam,i,j)=ldiv(gmul(p1,gcoeff(lam,i,j)),p2);
768: B[k+1]=B[k]=(long)p1;
769: for (i=k+2; i<=lx; i++)
770: B[i]=ldiv(gmul(p1,(GEN)B[i]),p2);
771: }
772: else
773: {
774: coeff(lam,k,k-1)=zero;
775: for (i=k+1; i<lx; i++)
776: { coeff(lam,i,k)=coeff(lam,i,k-1); coeff(lam,i,k-1)=zero; }
777: B[k]=B[k-1]; fl[k]=1; fl[k-1]=0;
778: }
779: }
780: if (k>2) k--;
781: }
782: else
783: {
784: for (l=k-2; l>=1; l--)
785: {
786: u=gcoeff(lam,k,l);
787: if (pslg(u)>=pslg((GEN)B[l+1]))
788: {
789: q=gdeuc(u,(GEN)B[l+1]); cq=gdivsg(1,content(q));
790: q=gmul(q,cq); flc=1;
791: h[k]=lsub(gmul(cq,(GEN)h[k]),gmul(q,(GEN)h[l]));
792: coeff(lam,k,l)=lsub(gmul(cq,gcoeff(lam,k,l)),gmul(q,(GEN)B[l+1]));
793: for (i=1; i<l; i++)
794: coeff(lam,k,i)=lsub(gmul(cq,gcoeff(lam,k,i)),gmul(q,gcoeff(lam,l,i)));
795: }
796: }
797: if (++k > n) break;
798: }
799: if (low_stack(lim, stack_lim(av,1)))
800: {
801: GEN *gptr[4];
802: if(DEBUGMEM>1) err(warnmem,"lllgramallgen");
803: gptr[0]=&B; gptr[1]=&lam; gptr[2]=&h;
804: gerepilemany(av,gptr,3);
805: }
806: }
807: tetpil=avma;
808: return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
809: }
810:
811: /* compute B[k], update mu(k,1..k-1) */
812: static int
813: get_Gram_Schmidt(GEN x, GEN mu, GEN B, long k)
814: {
815: GEN s,A = cgetg(k+1, t_COL); /* scratch space */
816: long av,i,j;
817: A[1] = coeff(x,k,1);
818: for(j=1;j<k;)
819: {
820: coeff(mu,k,j) = ldiv((GEN)A[j],(GEN)B[j]);
821: j++; av = avma;
822: /* A[j] <-- x[k,j] - sum_{i<j} mu[j,i] A[i] */
823: s = gmul(gcoeff(mu,j,1),(GEN)A[1]);
824: for (i=2; i<j; i++) s = gadd(s, gmul(gcoeff(mu,j,i),(GEN)A[i]));
825: s = gneg(s); A[j] = lpileupto(av, gadd(gcoeff(x,k,j),s));
826: }
827: B[k] = A[k]; return (gsigne((GEN)B[k]) > 0);
828: }
829:
830: /* x = Gram(b_i). If precision problems return NULL if flag=1, error otherwise.
831: * Quality ratio = Q = (alpha-1)/alpha. Suggested value: alpha = 100
832: */
833: GEN
834: lllgramintern(GEN x, long alpha, long flag, long prec)
835: {
836: GEN xinit,L,h,B,L1,QR;
837: long retry = 2, av = avma,lim,l,i,j,k,k1,lx=lg(x),n,kmax,KMAX;
838: long last_prec;
839:
840: if (typ(x) != t_MAT) err(typeer,"lllgram");
841: n=lx-1; if (n && lg((GEN)x[1])!=lx) err(mattype1,"lllgram");
842: if (n<=1) return idmat(n);
843: xinit = x; lim = stack_lim(av,1);
844: for (k=2,j=1; j<lx; j++)
845: {
846: GEN p1=(GEN)x[j];
847: for (i=1; i<=j; i++)
848: if (typ(p1[i]) == t_REAL) { l = lg((GEN)p1[i]); if (l>k) k=l; }
849: }
850: if (k == 2)
851: {
852: if (!prec) return lllgramint(x);
853: x = gmul(x, realun(prec+1));
854: }
855: else
856: {
857: if (prec < k) prec = k;
858: x = gprec_w(x, prec+1);
859: }
860: /* kmax = maximum column index attained during this run
861: * KMAX = same over all runs (after PRECPB) */
862: kmax = KMAX = 1;
863:
864: PRECPB:
865: switch(retry--)
866: {
867: case 2: h = idmat(n); break; /* entry */
868: case 1:
869: if (kmax > 2)
870: { /* some progress but precision loss, try again */
871: prec = (prec<<1)-2; kmax = 1;
872: if (DEBUGLEVEL > 3) fprintferr("\n");
873: if (DEBUGLEVEL) err(warnprec,"lllgramintern",prec);
874: x = qf_base_change(gprec_w(xinit,prec),h,1);
875: {
876: GEN *gsav[2]; gsav[0]=&h; gsav[1]=&x;
877: gerepilemany(av, gsav, 2);
878: }
879: if (DEBUGLEVEL) err(warner,"lllgramintern starting over");
880: break;
881: } /* fall through */
882: case 0: /* give up */
883: if (DEBUGLEVEL > 3) fprintferr("\n");
884: if (DEBUGLEVEL) err(warner,"lllgramintern giving up");
885: if (flag) { avma=av; return NULL; }
886: if (DEBUGLEVEL) outerr(xinit);
887: err(lllger3);
888: }
889: last_prec = prec;
890: QR = cgetr(prec+1); affsr(alpha-1,QR);
891: QR = divrs(QR,alpha);
892:
893: L=cgetg(lx,t_MAT);
894: B=cgetg(lx,t_COL);
895: for (j=1; j<lx; j++) { L[j] = (long)zerocol(n); B[j] = zero; }
896: k=2; B[1]=coeff(x,1,1);
897: if (gcmp0((GEN)B[1]))
898: {
899: if (flag) return NULL;
900: err(lllger3);
901: }
902: if (DEBUGLEVEL>5) fprintferr("k =");
903: for(;;)
904: {
905: if (k>kmax)
906: {
907: kmax = k; if (KMAX < kmax) KMAX = kmax;
908: if (DEBUGLEVEL>3) {fprintferr(" K%ld",k);flusherr();}
909: if (!get_Gram_Schmidt(x,L,B,k)) goto PRECPB;
910: }
911: else if (DEBUGLEVEL>5) fprintferr(" %ld",k);
912: L1 = gcoeff(L,k,k-1);
913: if (typ(L1) == t_REAL &&
914: (2*gexpo(L1) > bit_accuracy(lg(L1)) || 2*lg(L1) < last_prec))
915: {
916: last_prec = lg(L1);
917: if (DEBUGLEVEL>3)
918: fprintferr("\nRecomputing Gram-Schmidt, kmax = %ld, prec was %ld\n",
919: kmax,last_prec);
920: for (k1=1; k1<=kmax; k1++)
921: if (!get_Gram_Schmidt(x,L,B,k1)) goto PRECPB;
922: }
923: RED(x,h,L,KMAX,k,k-1);
924: if (do_SWAP(x,h,L,B,kmax,k,QR))
925: {
926: if (!B[k]) goto PRECPB;
927: if (k>2) k--;
928: }
929: else
930: {
931: for (l=k-2; l; l--) RED(x,h,L,KMAX,k,l);
932: if (++k > n) break;
933: }
934: if (low_stack(lim, stack_lim(av,1)))
935: {
936: GEN *gptr[5]; gptr[0]=&B; gptr[1]=&L; gptr[2]=&h; gptr[3]=&x; gptr[4]=&QR;
937: if(DEBUGMEM>1)
938: {
939: if (DEBUGLEVEL > 3) fprintferr("\n");
940: err(warnmem,"lllgram");
941: }
942: gerepilemany(av,gptr,5);
943: }
944: }
945: if (DEBUGLEVEL>3) fprintferr("\n");
946: return gerepilecopy(av, h);
947: }
948:
949: static GEN
950: lllgram_noerr(GEN x,long prec) { return lllgramintern(x,100,1,prec); }
951:
952: GEN
953: lllgram(GEN x,long prec) { return lllgramintern(x,100,0,prec); }
954:
955: GEN
956: qflll0(GEN x, long flag, long prec)
957: {
958: switch(flag)
959: {
960: case 0: return lll(x,prec);
961: case 1: return lllint(x);
962: case 2: return lllintpartial(x);
963: case 3: return lllrat(x);
964: case 4: return lllkerim(x);
965: case 5: return lllkerimgen(x);
966: case 7: return lll1(x,prec);
967: case 8: return lllgen(x);
968: case 9: return lllintwithcontent(x);
969: default: err(flagerr,"qflll");
970: }
971: return NULL; /* not reached */
972: }
973:
974: GEN
975: qflllgram0(GEN x, long flag, long prec)
976: {
977: switch(flag)
978: {
979: case 0: return lllgram(x,prec);
980: case 1: return lllgramint(x);
981: case 4: return lllgramkerim(x);
982: case 5: return lllgramkerimgen(x);
983: case 7: return lllgram1(x,prec);
984: case 8: return lllgramgen(x);
985: default: err(flagerr,"qflllgram");
986: }
987: return NULL; /* not reached */
988: }
989:
990: /* x est la matrice d'une base b_i; retourne la matrice u (entiere
991: * unimodulaire) d'une base LLL-reduite c_i en fonction des b_i (la base
992: * reduite est c=b*u).
993: */
994: static GEN
995: lll_proto(GEN x, GEN f(GEN, long), long prec)
996: {
997: long lx=lg(x),i,j,av,av1;
998: GEN g;
999:
1000: if (typ(x) != t_MAT) err(typeer,"lll_proto");
1001: if (lx==1) return cgetg(1,t_MAT);
1002: av=avma; g=cgetg(lx,t_MAT);
1003: for (j=1; j<lx; j++) g[j]=lgetg(lx,t_COL);
1004: for (i=1; i<lx; i++)
1005: for (j=1; j<=i; j++)
1006: coeff(g,i,j)=coeff(g,j,i)=(long)gscal((GEN)x[i],(GEN)x[j]);
1007: av1=avma; x = f(g,prec);
1008: if (!x) { avma=av; return NULL; }
1009: return gerepile(av,av1,x);
1010: }
1011:
1012: GEN
1013: lllintern(GEN x,long flag,long prec)
1014: {
1015: return lll_proto(x,flag? lllgram_noerr: lllgram,prec);
1016: }
1017:
1018: GEN
1019: lll(GEN x,long prec) { return lll_proto(x,lllgram,prec); }
1020:
1021: GEN
1022: lll1(GEN x,long prec) { return lll_proto(x,lllgram1,prec); }
1023:
1024: GEN
1025: lllrat(GEN x) { return lll_proto(x,lllgram,lll_ALL); }
1026:
1027: GEN
1028: lllint(GEN x) { return lll_proto(x,lllgramall0,lll_IM); }
1029:
1030: GEN
1031: lllgen(GEN x) { return lll_proto(x,lllgramallgen,lll_IM); }
1032:
1033: GEN
1034: lllgramgen(GEN x) { return lllgramallgen(x,lll_IM); }
1035:
1036: GEN
1037: lllgramkerimgen(GEN x) { return lllgramallgen(x,lll_ALL); }
1038:
1039: static GEN
1040: lllkerim_proto(GEN x, GEN f(GEN,long))
1041: {
1042: long lx=lg(x), i,j,av,av1;
1043: GEN g;
1044:
1045: if (typ(x) != t_MAT) err(typeer,"lllkerim_proto");
1046: if (lx==1)
1047: {
1048: g=cgetg(3,t_VEC);
1049: g[1]=lgetg(1,t_MAT);
1050: g[2]=lgetg(1,t_MAT); return g;
1051: }
1052: if (lg((GEN)x[1])==1)
1053: {
1054: g=cgetg(3,t_VEC);
1055: g[1]=(long)idmat(lx-1);
1056: g[2]=lgetg(1,t_MAT); return g;
1057: }
1058: av=avma; g=cgetg(lx,t_MAT);
1059: for (j=1; j<lx; j++) g[j]=lgetg(lx,t_COL);
1060: for (i=1; i<lx; i++)
1061: for (j=1; j<=i; j++)
1062: coeff(g,i,j)=coeff(g,j,i)=(long)gscal((GEN)x[i],(GEN)x[j]);
1063: av1=avma; return gerepile(av,av1,f(g,lll_ALL));
1064: }
1065:
1066: GEN
1067: lllkerim(GEN x) { return lllkerim_proto(x,lllgramall0); }
1068:
1069: GEN
1070: lllkerimgen(GEN x) { return lllkerim_proto(x,lllgramallgen); }
1071:
1072: /* x est ici la matrice de GRAM des bi. */
1073: GEN
1074: lllgram1(GEN x, long prec)
1075: {
1076: GEN mu,u,B,BB,BK,p,q,r,cst,unreel,sv,mu1,mu2;
1077: long av,tetpil,lim,l,i,j,k,lx=lg(x),n,e;
1078:
1079: if (typ(x) != t_MAT) err(typeer,"lllgram1");
1080: if (lg((GEN)x[1])!=lx) err(mattype1,"lllgram1"); n=lx-1;
1081: if (n<=1) return idmat(n);
1082: cst=gdivgs(stoi(99),100); /* LLL proposent 0.75 */
1083: if (prec)
1084: {
1085: unreel = realun(prec+1);
1086: x = gmul(x,unreel);
1087: cst = gmul(cst,unreel);
1088: }
1089: av=avma; lim=stack_lim(av,1);
1090: mu=gtrans(sqred(x)); B=cgetg(lx,t_COL);
1091: for (i=1,l=0; i<=n; i++)
1092: {
1093: if (gsigne((GEN)(B[i]=coeff(mu,i,i)))>0) l++;
1094: coeff(mu,i,i)=un;
1095: }
1096: if (l<n) err(lllger3);
1097:
1098: u=idmat(n); k=2;
1099: do
1100: {
1101: if (!gcmp0(r=grndtoi(gcoeff(mu,k,k-1),&e)))
1102: {
1103: u[k]=lsub((GEN)u[k],gmul(r,(GEN)u[k-1]));
1104: for (j=1; j<k-1; j++)
1105: coeff(mu,k,j)=lsub(gcoeff(mu,k,j),gmul(r,gcoeff(mu,k-1,j)));
1106: mu1=(GEN)(coeff(mu,k,k-1)=lsub(gcoeff(mu,k,k-1),r));
1107: }
1108: else mu1=gcoeff(mu,k,k-1);
1109: q=gmul((GEN)B[k-1],gsub(cst,mu2=gsqr(mu1)));
1110: if (gcmp(q,(GEN)B[k])>0)
1111: {
1112: BB=gadd((GEN)B[k],gmul((GEN)B[k-1],mu2));
1113: coeff(mu,k,k-1)=ldiv(gmul(mu1,(GEN)B[k-1]),BB);
1114: B[k]=lmul((GEN)B[k-1],BK=gdiv((GEN)B[k],BB));
1115: B[k-1]=(long)BB;
1116: swap(u[k-1],u[k]);
1117: for (j=1; j<=k-2; j++) swap(coeff(mu,k-1,j), coeff(mu,k,j));
1118: for (i=k+1; i<=n; i++)
1119: {
1120: p=gcoeff(mu,i,k);
1121: coeff(mu,i,k)=lsub(gcoeff(mu,i,k-1),gmul(mu1,p));
1122: coeff(mu,i,k-1)=ladd(gmul(BK,p),gmul(gcoeff(mu,k,k-1),gcoeff(mu,i,k-1)));
1123: }
1124: if (k>2) k--;
1125: }
1126: else
1127: {
1128: for (l=k-2; l; l--)
1129: {
1130: if (!gcmp0(r=grndtoi(gcoeff(mu,k,l),&e)))
1131: {
1132: u[k]=lsub((GEN)u[k],gmul(r,(GEN)u[l]));
1133: for (j=1; j<l; j++)
1134: coeff(mu,k,j)=lsub(gcoeff(mu,k,j),gmul(r,gcoeff(mu,l,j)));
1135: coeff(mu,k,l)=lsub(gcoeff(mu,k,l),r);
1136: }
1137: }
1138: k++;
1139: }
1140: if (low_stack(lim, stack_lim(av,1)))
1141: {
1142: if(DEBUGMEM>1) err(warnmem,"lllgram1");
1143: tetpil=avma;
1144: sv=cgetg(4,t_VEC);
1145: sv[1]=lcopy(B); sv[2]=lcopy(u); sv[3]=lcopy(mu);
1146: sv=gerepile(av,tetpil,sv);
1147: B=(GEN)sv[1]; u=(GEN)sv[2]; mu=(GEN)sv[3];
1148: }
1149: }
1150: while (k<=n);
1151: return gerepilecopy(av,u);
1152: }
1153:
1154: GEN
1155: lllgramint(GEN x)
1156: {
1157: return lllgramall0(x,lll_IM);
1158: }
1159:
1160: GEN
1161: lllgramkerim(GEN x)
1162: {
1163: return lllgramall0(x,lll_ALL);
1164: }
1165:
1166: /* This routine is functionally similar to lllint when all = 0.
1167: *
1168: * The input is an integer matrix mat (not necessarily square) whose
1169: * columns are linearly independent. It outputs another matrix T such that
1170: * mat * T is partially reduced. If all = 0, the size of mat * T is the
1171: * same as the size of mat. If all = 1 the number of columns of mat * T
1172: * is at most equal to its number of rows. A matrix M is said to be
1173: * -partially reduced- if | m1 +- m2 | >= |m1| for any two distinct
1174: * columns m1, m2, in M.
1175: *
1176: * This routine is designed to quickly reduce lattices in which one row
1177: * is huge compared to the other rows. For example, when searching for a
1178: * polynomial of degree 3 with root a mod p, the four input vectors might
1179: * be the coefficients of
1180: * X^3 - (a^3 mod p), X^2 - (a^2 mod p), X - (a mod p), p.
1181: * All four constant coefficients are O(p) and the rest are O(1). By the
1182: * pigeon-hole principle, the coefficients of the smallest vector in the
1183: * lattice are O(p^(1/4)), Hence significant reduction of vector lengths
1184: * can be anticipated.
1185: *
1186: * Peter Montgomery (July, 1994)
1187: *
1188: * If flag = 1 complete the reduction using lllint, otherwise return
1189: * partially reduced basis.
1190: */
1191: GEN
1192: lllintpartialall(GEN m, long flag)
1193: {
1194: const long ncol = lg(m)-1;
1195: const long ltop1 = avma;
1196: long ncolnz;
1197: GEN tm1, tm2, mid, *gptr[4];
1198:
1199: if (typ(m) != t_MAT) err(typeer,"lllintpartialall");
1200: if (ncol <= 1) return idmat(ncol);
1201:
1202: {
1203: GEN dot11 = sqscali((GEN)m[1]);
1204: GEN dot22 = sqscali((GEN)m[2]);
1205: GEN dot12 = gscali((GEN)m[1], (GEN)m[2]);
1206: GEN tm = idmat(2); /* For first two columns only */
1207:
1208: int progress = 0;
1209: long npass2 = 0;
1210:
1211: /* Try to row reduce the first two columns of m.
1212: * Our best result so far is (first two columns of m)*tm.
1213: *
1214: * Initially tm = 2 x 2 identity matrix.
1215: * The inner products of the reduced matrix are in
1216: * dot11, dot12, dot22.
1217: */
1218: while (npass2 < 2 || progress)
1219: {
1220: GEN dot12new,q = gdivround(dot12, dot22);
1221:
1222: npass2++; progress = signe(q);
1223: if (progress)
1224: {
1225: /* Conceptually replace (v1, v2) by (v1 - q*v2, v2),
1226: * where v1 and v2 represent the reduced basis for the
1227: * first two columns of the matrix.
1228: *
1229: * We do this by updating tm and the inner products.
1230: *
1231: * An improved algorithm would look only at the leading
1232: * digits of dot11, dot12, dot22. It would use
1233: * single-precision calculations as much as possible.
1234: */
1235: q = negi(q);
1236: dot12new = addii(dot12, mulii(q, dot22));
1237: dot11 = addii(dot11, mulii(q, addii(dot12, dot12new)));
1238: dot12 = dot12new;
1239: tm[1] = (long)ZV_lincomb(gun,q, (GEN)tm[1],(GEN)tm[2]);
1240: }
1241:
1242: /* Interchange the output vectors v1 and v2. */
1243: gswap(dot11,dot22); swap(tm[1],tm[2]);
1244:
1245: /* Occasionally (including final pass) do garbage collection. */
1246: if (npass2 % 8 == 0 || !progress)
1247: {
1248: gptr[0] = &dot11; gptr[1] = &dot12;
1249: gptr[2] = &dot22; gptr[3] = &tm;
1250: gerepilemany(ltop1, gptr, 4);
1251: }
1252: } /* while npass2 < 2 || progress */
1253:
1254: {
1255: long icol;
1256: GEN det12 = subii(mulii(dot11, dot22), mulii(dot12, dot12));
1257:
1258: tm1 = idmat(ncol);
1259: mid = cgetg(ncol+1, t_MAT);
1260: for (icol = 1; icol <= 2; icol++)
1261: {
1262: coeff(tm1,1,icol) = coeff(tm,1,icol);
1263: coeff(tm1,2,icol) = coeff(tm,2,icol);
1264: mid[icol] = (long)ZV_lincomb(
1265: gcoeff(tm,1,icol),gcoeff(tm,2,icol), (GEN)m[1],(GEN)m[2]);
1266: }
1267: for (icol = 3; icol <= ncol; icol++)
1268: {
1269: GEN curcol = (GEN)m[icol];
1270: GEN dot1i = gscali((GEN)mid[1], curcol);
1271: GEN dot2i = gscali((GEN)mid[2], curcol);
1272: /* Try to solve
1273: *
1274: * ( dot11 dot12 ) (q1) ( dot1i )
1275: * ( dot12 dot22 ) (q2) = ( dot2i )
1276: *
1277: * Round -q1 and -q2 to the nearest integer.
1278: * Then compute curcol - q1*mid[1] - q2*mid[2].
1279: * This will be approximately orthogonal to the
1280: * first two vectors in the new basis.
1281: */
1282: GEN q1neg = subii(mulii(dot12, dot2i), mulii(dot22, dot1i));
1283: GEN q2neg = subii(mulii(dot12, dot1i), mulii(dot11, dot2i));
1284:
1285: q1neg = gdivround(q1neg, det12);
1286: q2neg = gdivround(q2neg, det12);
1287: coeff(tm1, 1, icol) = ladd(gmul(q1neg, gcoeff(tm,1,1)),
1288: gmul(q2neg, gcoeff(tm,1,2)));
1289: coeff(tm1, 2, icol) = ladd(gmul(q1neg, gcoeff(tm,2,1)),
1290: gmul(q2neg, gcoeff(tm,2,2)));
1291: mid[icol] = ladd(curcol,
1292: ZV_lincomb(q1neg,q2neg, (GEN)mid[1],(GEN)mid[2]));
1293: } /* for icol */
1294: } /* local block */
1295: }
1296: if (DEBUGLEVEL>6)
1297: {
1298: fprintferr("tm1 = %Z",tm1);
1299: fprintferr("mid = %Z",mid);
1300: }
1301: gptr[0] = &tm1; gptr[1] = ∣
1302: gerepilemany(ltop1, gptr, 2);
1303: {
1304: /* For each pair of column vectors v and w in mid * tm2,
1305: * try to replace (v, w) by (v, v - q*w) for some q.
1306: * We compute all inner products and check them repeatedly.
1307: */
1308: const long ltop3 = avma; /* Excludes region with tm1 and mid */
1309: long icol, lim, reductions, npass = 0;
1310: GEN dotprd = cgetg(ncol+1, t_MAT);
1311:
1312: tm2 = idmat(ncol);
1313: for (icol=1; icol <= ncol; icol++)
1314: {
1315: long jcol;
1316:
1317: dotprd[icol] = lgetg(ncol+1,t_COL);
1318: for (jcol=1; jcol <= icol; jcol++)
1319: coeff(dotprd,jcol,icol) = coeff(dotprd,icol,jcol) =
1320: (long)gscali((GEN)mid[icol], (GEN)mid[jcol]);
1321: } /* for icol */
1322: lim = stack_lim(ltop3,1);
1323: for(;;)
1324: {
1325: reductions = 0;
1326: for (icol=1; icol <= ncol; icol++)
1327: {
1328: long ijdif, jcol, k1, k2;
1329: GEN codi, q;
1330:
1331: for (ijdif=1; ijdif < ncol; ijdif++)
1332: {
1333: const long previous_avma = avma;
1334:
1335: jcol = (icol + ijdif - 1) % ncol; jcol++; /* Hack for NeXTgcc 2.5.8 */
1336: k1 = (cmpii(gcoeff(dotprd,icol,icol),
1337: gcoeff(dotprd,jcol,jcol) ) > 0)
1338: ? icol : jcol; /* index of larger column */
1339: k2 = icol + jcol - k1; /* index of smaller column */
1340: codi = gcoeff(dotprd,k2,k2);
1341: q = gcmp0(codi)? gzero: gdivround(gcoeff(dotprd,k1,k2), codi);
1342:
1343: /* Try to subtract a multiple of column k2 from column k1. */
1344: if (gcmp0(q)) avma = previous_avma;
1345: else
1346: {
1347: long dcol;
1348:
1349: reductions++; q = negi(q);
1350: tm2[k1]=(long)
1351: ZV_lincomb(gun,q, (GEN)tm2[k1], (GEN)tm2[k2]);
1352: dotprd[k1]=(long)
1353: ZV_lincomb(gun,q, (GEN)dotprd[k1], (GEN)dotprd[k2]);
1354: coeff(dotprd, k1, k1) = laddii(gcoeff(dotprd,k1,k1),
1355: mulii(q, gcoeff(dotprd,k2,k1)));
1356: for (dcol = 1; dcol <= ncol; dcol++)
1357: coeff(dotprd,k1,dcol) = coeff(dotprd,dcol,k1);
1358: } /* if q != 0 */
1359: } /* for ijdif */
1360: if (low_stack(lim, stack_lim(ltop3,1)))
1361: {
1362: if(DEBUGMEM>1) err(warnmem,"lllintpartialall");
1363: gptr[0] = &dotprd; gptr[1] = &tm2;
1364: gerepilemany(ltop3, gptr, 2);
1365: }
1366: } /* for icol */
1367: if (!reductions) break;
1368: if (DEBUGLEVEL>6)
1369: {
1370: GEN diag_prod = dbltor(1.0);
1371: for (icol = 1; icol <= ncol; icol++)
1372: diag_prod = gmul(diag_prod, gcoeff(dotprd,icol,icol));
1373: npass++;
1374: fprintferr("npass = %ld, red. last time = %ld, diag_prod = %Z\n\n",
1375: npass, reductions, diag_prod);
1376: }
1377: } /* for(;;)*/
1378:
1379: /* Sort columns so smallest comes first in m * tm1 * tm2.
1380: * Use insertion sort.
1381: */
1382: for (icol = 1; icol < ncol; icol++)
1383: {
1384: long jcol, s = icol;
1385:
1386: for (jcol = icol+1; jcol <= ncol; jcol++)
1387: if (cmpii(gcoeff(dotprd,s,s),gcoeff(dotprd,jcol,jcol)) > 0) s = jcol;
1388: if (icol != s)
1389: { /* Exchange with proper column */
1390: /* Only diagonal of dotprd is updated */
1391: swap(tm2[icol], tm2[s]);
1392: swap(coeff(dotprd,icol,icol), coeff(dotprd,s,s));
1393: }
1394: } /* for icol */
1395: icol=1;
1396: while (icol <= ncol && !signe(gcoeff(dotprd,icol,icol))) icol++;
1397: ncolnz = ncol - icol + 1;
1398: } /* local block */
1399:
1400: if (flag)
1401: {
1402: if (ncolnz == lg((GEN)m[1])-1)
1403: {
1404: tm2 += (ncol-ncolnz);
1405: tm2[0]=evaltyp(t_MAT)|evallg(ncolnz+1);
1406: }
1407: else
1408: {
1409: tm1 = gmul(tm1, tm2);
1410: tm2 = lllint(gmul(m, tm1));
1411: }
1412: }
1413: if (DEBUGLEVEL>6)
1414: fprintferr("lllintpartial output = %Z", gmul(tm1, tm2));
1415: return gerepileupto(ltop1, gmul(tm1, tm2));
1416: }
1417:
1418: GEN
1419: lllintpartial(GEN mat)
1420: {
1421: return lllintpartialall(mat,1);
1422: }
1423:
1424: /********************************************************************/
1425: /** **/
1426: /** LINEAR & ALGEBRAIC DEPENDANCE **/
1427: /** **/
1428: /********************************************************************/
1429:
1430: GEN
1431: lindep0(GEN x,long bit,long prec)
1432: {
1433: if (!bit) return lindep(x,prec);
1434: if (bit>0) return lindep2(x,bit);
1435: return deplin(x);
1436: }
1437:
1438: static int
1439: real_indep(GEN re, GEN im, long bitprec)
1440: {
1441: GEN p1 = gsub(gmul((GEN)re[1],(GEN)im[2]),
1442: gmul((GEN)re[2],(GEN)im[1]));
1443: return (!gcmp0(p1) && gexpo(p1) > - bitprec);
1444: }
1445:
1446: GEN
1447: lindep2(GEN x, long bit)
1448: {
1449: long tx=typ(x),lx=lg(x),ly,i,j,e, av = avma;
1450: GEN re,im,p1,p2;
1451:
1452: if (! is_vec_t(tx)) err(typeer,"lindep2");
1453: if (lx<=2) return cgetg(1,t_VEC);
1454: re = greal(x);
1455: im = gimag(x); bit = (long) (bit/L2SL10);
1456: /* independant over R ? */
1457: if (lx == 3 && real_indep(re,im,bit))
1458: { avma = av; return cgetg(1, t_VEC); }
1459:
1460: if (gcmp0(im)) im = NULL;
1461: ly = im? lx+2: lx+1;
1462: p2=cgetg(lx,t_MAT);
1463: for (i=1; i<lx; i++)
1464: {
1465: p1 = cgetg(ly,t_COL); p2[i] = (long)p1;
1466: for (j=1; j<lx; j++) p1[j] = (i==j)? un: zero;
1467: p1[lx] = lcvtoi(gshift((GEN)re[i],bit),&e);
1468: if (im) p1[lx+1] = lcvtoi(gshift((GEN)im[i],bit),&e);
1469: }
1470: p1 = (GEN)gmul(p2,lllint(p2))[1];
1471: p1[0] = evaltyp(t_VEC) | evallg(lx);
1472: return gerepilecopy(av, p1);
1473: }
1474:
1475: #define quazero(x) (gcmp0(x) || (typ(x)==t_REAL && expo(x) < EXP))
1476: GEN
1477: lindep(GEN x, long prec)
1478: {
1479: GEN *b,*be,*bn,**m,qzer;
1480: GEN c1,c2,c3,px,py,pxy,re,im,p1,p2,r,f,em;
1481: long i,j,fl,i1, lx = lg(x), tx = typ(x), n = lx-1;
1482: long av = avma, lim = stack_lim(av,1), av0,av1,tetpil;
1483: const long EXP = - bit_accuracy(prec) + 2*n;
1484:
1485: if (! is_vec_t(tx)) err(typeer,"lindep");
1486: if (lx<=2) return cgetg(1,t_VEC);
1487: x = gmul(x, realun(prec));
1488: re=greal(x); im=gimag(x);
1489: /* independant over R ? */
1490: if (lx == 3 && real_indep(re,im,bit_accuracy(prec)))
1491: { avma = av; return cgetg(1, t_VEC); }
1492:
1493: qzer = new_chunk(lx);
1494: b = (GEN*)idmat(n);
1495: be= (GEN*)new_chunk(lx);
1496: bn= (GEN*)new_chunk(lx);
1497: m = (GEN**)new_chunk(lx);
1498: for (i=1; i<=n; i++)
1499: {
1500: bn[i]=cgetr(prec+1);
1501: be[i]=cgetg(lx,t_COL);
1502: m[i] = (GEN*)new_chunk(lx);
1503: for (j=1; j<i ; j++) m[i][j]=cgetr(prec+1);
1504: for (j=1; j<=n; j++) be[i][j]=lgetr(prec+1);
1505: }
1506: px=sqscal(re);
1507: py=sqscal(im); pxy=gscal(re,im);
1508: p1=mpsub(mpmul(px,py),gsqr(pxy));
1509: if (quazero(re)) { re=im; px=py; fl=1; } else fl=quazero(p1);
1510: av0 = av1 = avma;
1511: for (i=1; i<=n; i++)
1512: {
1513: p2 = gscal(b[i],re);
1514: if (fl) p2=gmul(gdiv(p2,px),re);
1515: else
1516: {
1517: GEN p5,p6,p7;
1518: p5=gscal(b[i],im);
1519: p6=gdiv(gsub(gmul(py,p2),gmul(pxy,p5)),p1);
1520: p7=gdiv(gsub(gmul(px,p5),gmul(pxy,p2)),p1);
1521: p2=gadd(gmul(p6,re),gmul(p7,im));
1522: }
1523: if (tx!=t_COL) p2=gtrans(p2);
1524: p2=gsub(b[i],p2);
1525: for (j=1; j<i; j++)
1526: if (qzer[j]) affrr(bn[j],m[i][j]);
1527: else
1528: {
1529: gdivz(gscal(b[i],be[j]),bn[j],m[i][j]);
1530: p2=gsub(p2,gmul(m[i][j],be[j]));
1531: }
1532: gaffect(p2,be[i]); affrr(sqscal(be[i]),bn[i]);
1533: qzer[i]=quazero(bn[i]); avma=av1;
1534: }
1535: while (qzer[n])
1536: {
1537: long e;
1538: if (DEBUGLEVEL>9)
1539: {
1540: fprintferr("qzer[%ld]=%ld\n",n,qzer[n]);
1541: for (i1=1; i1<=n; i1++)
1542: for (i=1; i<i1; i++) output(m[i1][i]);
1543: }
1544: em=bn[1]; j=1;
1545: for (i=2; i<n; i++)
1546: {
1547: p1=shiftr(bn[i],i);
1548: if (cmprr(p1,em)>0) { em=p1; j=i; }
1549: }
1550: i=j; i1=i+1;
1551: avma = av1; r = grndtoi(m[i1][i], &e);
1552: if (e >= 0) err(precer,"lindep");
1553: r = negi(r);
1554: p1 = ZV_lincomb(gun,r, b[i1],b[i]);
1555: av1 = avma;
1556: b[i1]=b[i]; b[i]=p1; f=addir(r,m[i1][i]);
1557: for (j=1; j<i; j++)
1558: if (!qzer[j])
1559: {
1560: p1=mpadd(m[i1][j],mulir(r,m[i][j]));
1561: affrr(m[i][j],m[i1][j]); mpaff(p1,m[i][j]);
1562: }
1563: c1=addrr(bn[i1],mulrr(gsqr(f),bn[i]));
1564: if (!quazero(c1))
1565: {
1566: c2=divrr(mulrr(bn[i],f),c1); affrr(c2,m[i1][i]);
1567: c3=divrr(bn[i1],c1); mulrrz(c3,bn[i],bn[i1]);
1568: affrr(c1,bn[i]); qzer[i1]=quazero(bn[i1]); qzer[i]=0;
1569: for (j=i+2; j<=n; j++)
1570: {
1571: p1=addrr(mulrr(m[j][i1],c3),mulrr(m[j][i],c2));
1572: subrrz(m[j][i],mulrr(f,m[j][i1]), m[j][i1]);
1573: affrr(p1,m[j][i]);
1574: }
1575: }
1576: else
1577: {
1578: qzer[i1]=qzer[i]; qzer[i]=1;
1579: affrr(bn[i],bn[i1]); affrr(c1,bn[i]);
1580: for (j=i+2; j<=n; j++) affrr(m[j][i],m[j][i1]);
1581: }
1582: if (low_stack(lim, stack_lim(av,1)))
1583: {
1584: if(DEBUGMEM>1) err(warnmem,"lindep");
1585: b = (GEN*)gerepilecopy(av0, (GEN)b);
1586: av1 = avma;
1587: }
1588: }
1589: p1=cgetg(lx,t_COL); p1[n]=un; for (i=1; i<n; i++) p1[i]=zero;
1590: p2 = (GEN)b; p2[0] = evaltyp(t_MAT) | evallg(lx);
1591: p2=gauss(gtrans(p2),p1); tetpil=avma;
1592: return gerepile(av,tetpil,gtrans(p2));
1593: }
1594:
1595: /* x is a vector of elts of a p-adic field */
1596: GEN
1597: plindep(GEN x)
1598: {
1599: long av = avma,i,j, prec = VERYBIGINT, lx = lg(x)-1, ly,v;
1600: GEN p = NULL, pn,p1,m,a;
1601:
1602: if (lx < 2) return cgetg(1,t_VEC);
1603: for (i=1; i<=lx; i++)
1604: {
1605: p1 = (GEN)x[i];
1606: if (typ(p1) != t_PADIC) continue;
1607:
1608: j = precp(p1); if (j < prec) prec = j;
1609: if (!p) p = (GEN)p1[2];
1610: else if (!egalii(p, (GEN)p1[2]))
1611: err(talker,"inconsistent primes in plindep");
1612: }
1613: if (!p) err(talker,"not a p-adic vector in plindep");
1614: v = ggval(x,p); pn = gpowgs(p,prec);
1615: if (v != 0) x = gmul(x, gpowgs(p, -v));
1616: x = lift_intern(gmul(x, gmodulcp(gun, pn)));
1617:
1618: ly = 2*lx - 1;
1619: m = cgetg(ly+1,t_MAT);
1620: for (j=1; j<=ly; j++)
1621: {
1622: p1 = cgetg(lx+1, t_COL); m[j] = (long)p1;
1623: for (i=1; i<=lx; i++) p1[i] = zero;
1624: }
1625: a = negi((GEN)x[1]);
1626: for (i=1; i<lx; i++)
1627: {
1628: coeff(m,1+i,i) = (long)a;
1629: coeff(m,1 ,i) = x[i+1];
1630: }
1631: for (i=1; i<=lx; i++) coeff(m,i,lx-1+i) = (long)pn;
1632: p1 = lllint(m);
1633: return gerepileupto(av, gmul(m, (GEN)p1[1]));
1634: }
1635:
1636: GEN
1637: algdep0(GEN x, long n, long bit, long prec)
1638: {
1639: long tx=typ(x),av,i,k;
1640: GEN y,p1;
1641:
1642: if (! is_scalar_t(tx)) err(typeer,"algdep0");
1643: if (tx==t_POLMOD) { y=forcecopy((GEN)x[1]); setvarn(y,0); return y; }
1644: if (gcmp0(x)) return gzero;
1645: if (!n) return gun;
1646:
1647: av=avma; p1=cgetg(n+2,t_COL);
1648: p1[1] = un;
1649: p1[2] = (long)x; /* n >= 1 */
1650: for (i=3; i<=n+1; i++) p1[i]=lmul((GEN)p1[i-1],x);
1651: if (typ(x) == t_PADIC)
1652: p1 = plindep(p1);
1653: else
1654: p1 = bit? lindep2(p1,bit): lindep(p1,prec);
1655: if (lg(p1) < 2)
1656: err(talker,"higher degree than expected in algdep");
1657:
1658: y=cgetg(n+3,t_POL);
1659: y[1] = evalsigne(1) | evalvarn(0);
1660: k=1; while (gcmp0((GEN)p1[k])) k++;
1661: for (i=0; i<=n+1-k; i++) y[i+2] = p1[k+i];
1662: normalizepol_i(y, n+4-k);
1663: y = (gsigne(leading_term(y)) > 0)? gcopy(y): gneg(y);
1664: return gerepileupto(av,y);
1665: }
1666:
1667: GEN
1668: algdep2(GEN x, long n, long bit)
1669: {
1670: return algdep0(x,n,bit,0);
1671: }
1672:
1673: GEN
1674: algdep(GEN x, long n, long prec)
1675: {
1676: return algdep0(x,n,0,prec);
1677: }
1678:
1679: /********************************************************************/
1680: /** **/
1681: /** INTEGRAL KERNEL (LLL REDUCED) **/
1682: /** **/
1683: /********************************************************************/
1684:
1685: GEN
1686: matkerint0(GEN x, long flag)
1687: {
1688: switch(flag)
1689: {
1690: case 0: return kerint(x);
1691: case 1: return kerint1(x);
1692: case 2: return kerint2(x);
1693: default: err(flagerr,"matkerint");
1694: }
1695: return NULL; /* not reached */
1696: }
1697:
1698: GEN
1699: kerint1(GEN x)
1700: {
1701: long av=avma,tetpil;
1702: GEN p1,p2;
1703:
1704: p1=matrixqz3(ker(x)); p2=lllint(p1); tetpil=avma;
1705: return gerepile(av,tetpil,gmul(p1,p2));
1706: }
1707:
1708: GEN
1709: kerint2(GEN x)
1710: {
1711: long lx=lg(x), i,j,av,av1;
1712: GEN g,p1;
1713:
1714: if (typ(x) != t_MAT) err(typeer,"kerint2");
1715: av=avma; g=cgetg(lx,t_MAT);
1716: for (j=1; j<lx; j++) g[j]=lgetg(lx,t_COL);
1717: for (i=1; i<lx; i++)
1718: for (j=1; j<=i; j++)
1719: coeff(g,i,j) = coeff(g,j,i) = (long)gscal((GEN)x[i],(GEN)x[j]);
1720: g=lllgramall0(g,lll_KER); p1=lllint(g);
1721: av1=avma; return gerepile(av,av1,gmul(g,p1));
1722: }
1723:
1724: static GEN
1725: lllall0(GEN x, long flag)
1726: {
1727: long av0=avma,av,tetpil,lx=lg(x),i,j,k,l,n,lim,kmax;
1728: GEN u,B,L,q,r,h,la,p1,p2,p4,fl;
1729:
1730: if (typ(x) != t_MAT) err(typeer,"lllall0");
1731: n=lx-1; if (n<=1) return lllall_trivial(x,n, flag | lll_GRAM);
1732: fl = new_chunk(lx);
1733:
1734: av=avma; lim=stack_lim(av,1); x=dummycopy(x);
1735: B=gscalcol(gun, lx);
1736: L=cgetg(lx,t_MAT);
1737: for (k=lg(x[1]),j=1; j<lx; j++)
1738: {
1739: for (i=1; i<k; i++)
1740: if (typ(gcoeff(x,i,j))!=t_INT) err(lllger4);
1741: fl[j] = 0; L[j] = (long)zerocol(n);
1742: }
1743: k=2; h=idmat(n); kmax=1;
1744: u=sqscali((GEN)x[1]);
1745: if (signe(u)) { B[2]=(long)u; coeff(L,1,1)=un; fl[1]=1; } else B[2]=un;
1746: for(;;)
1747: {
1748: if (k > kmax)
1749: {
1750: kmax = k;
1751: for (j=1; j<=k; j++)
1752: {
1753: if (j==k || fl[j])
1754: {
1755: long av1 = avma;
1756: u=gscali((GEN)x[k],(GEN)x[j]);
1757: for (i=1; i<j; i++)
1758: if (fl[i])
1759: u = divii(subii(mulii((GEN)B[i+1],u),
1760: mulii(gcoeff(L,k,i),gcoeff(L,j,i))),
1761: (GEN)B[i]);
1762: u = gerepileuptoint(av1, u);
1763: if (j<k) coeff(L,k,j)=(long)u;
1764: else
1765: {
1766: if (signe(u)) { B[k+1]=(long)u; coeff(L,k,k)=un; fl[k]=1; }
1767: else { B[k+1]=B[k]; fl[k]=0; }
1768: }
1769: }
1770: }
1771: }
1772: if (fl[k-1] && !fl[k])
1773: {
1774: u = shifti(gcoeff(L,k,k-1),1);
1775: if (absi_cmp(u, (GEN)B[k]) > 0)
1776: {
1777: q = truedvmdii(addii(u,(GEN)B[k]),shifti((GEN)B[k],1), NULL);
1778: r = negi(q);
1779: h[k] = (long)ZV_lincomb(gun,r, (GEN)h[k],(GEN)h[k-1]);
1780: x[k] = (long)ZV_lincomb(gun,r, (GEN)x[k],(GEN)x[k-1]);
1781: coeff(L,k,k-1)=laddii(gcoeff(L,k,k-1),mulii(r,(GEN)B[k]));
1782: for (i=1; i<k-1; i++)
1783: coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(r,gcoeff(L,k-1,i)));
1784: }
1785: la=gcoeff(L,k,k-1); p4=sqri(la);
1786: swap(h[k-1], h[k]);
1787: swap(x[k-1], x[k]);
1788: for (j=1; j<k-1; j++) swap(coeff(L,k-1,j), coeff(L,k,j));
1789: if (signe(la))
1790: {
1791: p2=(GEN)B[k]; p1=divii(p4,p2);
1792: B[k+1]=B[k]=(long)p1;
1793: for (i=k+1; i<=kmax; i++)
1794: coeff(L,i,k-1)=ldivii(mulii(la,gcoeff(L,i,k-1)),p2);
1795: for (j=k+1; j<kmax; j++)
1796: for (i=j+1; i<=kmax; i++)
1797: coeff(L,i,j)=ldivii(mulii((GEN)p1,gcoeff(L,i,j)),p2);
1798: for (i=k+2; i<=kmax+1; i++)
1799: B[i]=ldivii(mulii((GEN)p1,(GEN)B[i]),p2);
1800: }
1801: else
1802: {
1803: for (i=k+1; i<=kmax; i++)
1804: { coeff(L,i,k)=coeff(L,i,k-1); coeff(L,i,k-1)=zero; }
1805: B[k]=B[k-1]; fl[k]=1; fl[k-1]=0;
1806: }
1807: if (k>2) k--;
1808: }
1809: else
1810: {
1811: for (l=k-1; l>=1; l--)
1812: {
1813: u = shifti(gcoeff(L,k,l),1);
1814: if (absi_cmp(u,(GEN)B[l+1]) > 0)
1815: {
1816: q = truedvmdii(addii(u,(GEN)B[l+1]),shifti((GEN)B[l+1],1), NULL);
1817: r = negi(q);
1818: h[k] = (long)ZV_lincomb(gun,r,(GEN)h[k],(GEN)h[l]);
1819: x[k] = (long)ZV_lincomb(gun,r,(GEN)x[k],(GEN)x[l]);
1820: coeff(L,k,l)=laddii(gcoeff(L,k,l),mulii(r,(GEN)B[l+1]));
1821: for (i=1; i<l; i++)
1822: coeff(L,k,i)=laddii(gcoeff(L,k,i),mulii(r,gcoeff(L,l,i)));
1823: }
1824: }
1825: if (++k > n) break;
1826: }
1827: if (low_stack(lim, stack_lim(av,1)))
1828: {
1829: GEN *gptr[4];
1830: if(DEBUGMEM>1) err(warnmem,"lllall0");
1831: gptr[0]=&B; gptr[1]=&L; gptr[2]=&h;
1832: gptr[3]=&x; gerepilemany(av,gptr,4);
1833: }
1834: }
1835: tetpil=avma;
1836: return gerepile(av0,tetpil,lllgramall_finish(h,fl,flag,n));
1837: }
1838:
1839: GEN
1840: kerint(GEN x)
1841: {
1842: long av=avma,av1;
1843: GEN g,p1;
1844:
1845: g=lllall0(x,lll_KER); if (lg(g)==1) return g;
1846: p1=lllint(g); av1=avma;
1847: return gerepile(av,av1,gmul(g,p1));
1848: }
1849:
1850: /********************************************************************/
1851: /** **/
1852: /** POLRED & CO. **/
1853: /** **/
1854: /********************************************************************/
1855: /* remove duplicate polynomials in y, updating a (same indexes), in place */
1856: static long
1857: remove_duplicates(GEN y, GEN a)
1858: {
1859: long k,i, nv = lg(y), av = avma;
1860: GEN z;
1861:
1862: if (nv < 2) return nv;
1863: z = new_chunk(3);
1864: z[1] = (long)y;
1865: z[2] = (long)a; (void)sort_factor(z, gcmp);
1866: for (k=1, i=2; i<nv; i++)
1867: if (!gegal((GEN)y[i], (GEN)y[k]))
1868: {
1869: k++;
1870: a[k] = a[i];
1871: y[k] = y[i];
1872: }
1873: nv = k+1; setlg(a,nv); setlg(y,nv);
1874: avma = av; return nv;
1875: }
1876:
1877: /* in place; make sure second non-zero coeff is negative (choose x or -x) */
1878: int
1879: canon_pol(GEN z)
1880: {
1881: long i,s;
1882:
1883: for (i=lgef(z)-2; i>=2; i-=2)
1884: {
1885: s = signe(z[i]);
1886: if (s > 0)
1887: {
1888: for (; i>=2; i-=2) z[i]=lnegi((GEN)z[i]);
1889: return -1;
1890: }
1891: if (s) return 1;
1892: }
1893: return 0;
1894: }
1895:
1896: static GEN
1897: pols_for_polred(GEN x, GEN base, GEN LLLbase, GEN *pta,
1898: int (*check)(GEN, GEN), GEN arg)
1899: {
1900: long i, v = varn(x), n = lg(base);
1901: GEN p1,p2,p3,y, a = cgetg(n,t_VEC);
1902:
1903: for (i=1; i<n; i++) a[i] = lmul(base,(GEN)LLLbase[i]);
1904: y=cgetg(n,t_VEC);
1905: for (i=1; i<n; i++)
1906: {
1907: if (DEBUGLEVEL > 2) { fprintferr("i = %ld\n",i); flusherr(); }
1908: p1=(GEN)a[i];
1909: p1 = primitive_part(p1, &p3);
1910: p1 = caractducos(x,p1,v);
1911: if (p3) p1 = ZX_rescale_pol(p1,p3);
1912: p2 = modulargcd(derivpol(p1),p1);
1913: p3 = leading_term(p2); if (!gcmp1(p3)) p2=gdiv(p2,p3);
1914: p1 = gdiv(p1,p2);
1915: if (canon_pol(p1) < 0 && pta) a[i] = (long)gneg_i((GEN)a[i]);
1916: y[i] = (long)p1; if (DEBUGLEVEL > 3) outerr(p1);
1917: if (check && check(arg, p1)) return p1;
1918: }
1919: if (check) return NULL; /* no suitable polynomial found */
1920: (void)remove_duplicates(y,a);
1921: if (pta) *pta = a;
1922: return y;
1923: }
1924:
1925: GEN
1926: nf_get_T2(GEN base, GEN polr)
1927: {
1928: long i,j, n = lg(base);
1929: GEN p1,p2=cgetg(n,t_MAT);
1930:
1931: for (i=1; i<n; i++)
1932: {
1933: p1=cgetg(n,t_COL); p2[i]=(long)p1;
1934: for (j=1; j<n; j++)
1935: p1[j] = (long) poleval((GEN)base[i],(GEN)polr[j]);
1936: }
1937: return mulmat_real(gconj(gtrans(p2)),p2);
1938: }
1939:
1940: /* compute Tr(w_i w_j) */
1941: static GEN
1942: nf_get_T(GEN x, GEN w)
1943: {
1944: long i,j,k, n = degpol(x);
1945: GEN p1,p2,p3;
1946: GEN ptrace = cgetg(n+2,t_VEC);
1947: GEN den = cgetg(n+1,t_VEC);
1948: GEN T = cgetg(n+1,t_MAT);
1949:
1950: ptrace[2]=lstoi(n);
1951: for (k=2; k<=n; k++)
1952: { /* cf polsym */
1953: GEN y = x + (n-k+1);
1954: p1 = mulsi(k-1,(GEN)y[2]);
1955: for (i=3; i<=k; i++)
1956: p1 = addii(p1,mulii((GEN)y[i],(GEN)ptrace[i]));
1957: ptrace[i] = lnegi(p1);
1958: }
1959: w = dummycopy(w);
1960: for (i=1; i<=n; i++)
1961: {
1962: den[i] = (long)denom(content((GEN)w[i]));
1963: w[i] = lmul((GEN)w[i],(GEN)den[i]);
1964: }
1965:
1966: for (i=1; i<=n; i++)
1967: {
1968: p1=cgetg(n+1,t_COL); T[i]=(long)p1;
1969: for (j=1; j<i ; j++) p1[j] = coeff(T,i,j);
1970: for ( ; j<=n; j++)
1971: { /* cf quicktrace */
1972: p2 = gres(gmul((GEN)w[i],(GEN)w[j]),x);
1973: p3 = gzero;
1974: for (k=lgef(p2)-1; k>1; k--)
1975: p3 = addii(p3, mulii((GEN)p2[k],(GEN)ptrace[k]));
1976: p1[j]=(long)divii(p3, mulii((GEN)den[i],(GEN)den[j]));
1977: }
1978: }
1979: return T;
1980: }
1981:
1982: /* Return the base change matrix giving the an LLL-reduced basis for the
1983: * maximal order of the nf given by x. Expressed in terms of the standard
1984: * HNF basis (as polynomials) given in base (ignored if x is an nf)
1985: */
1986: GEN
1987: LLL_nfbasis(GEN *ptx, GEN polr, GEN base, long prec)
1988: {
1989: GEN T2,p1, x = *ptx;
1990: int totally_real,n,i;
1991:
1992: if (typ(x) != t_POL)
1993: {
1994: p1=checknf(x); *ptx=x=(GEN)p1[1];
1995: base=(GEN)p1[7]; n=degpol(x);
1996: totally_real = !signe(gmael(p1,2,2));
1997: T2=gmael(p1,5,3); if (totally_real) T2 = ground(T2);
1998: }
1999: else
2000: {
2001: n=degpol(x); totally_real = (!prec || sturm(x)==n);
2002: if (typ(base) != t_VEC || lg(base)-1 != n)
2003: err(talker,"incorrect Zk basis in LLL_nfbasis");
2004: if (!totally_real)
2005: T2 = nf_get_T2(base,polr? polr: roots(x,prec));
2006: else
2007: T2 = nf_get_T(x, base);
2008: }
2009: if (totally_real) return lllgramint(T2);
2010: for (i=1; ; i++)
2011: {
2012: if ((p1 = lllgramintern(T2,100,1,prec))) return p1;
2013: if (i == MAXITERPOL) err(accurer,"allpolred");
2014: prec=(prec<<1)-2;
2015: if (DEBUGLEVEL) err(warnprec,"allpolred",prec);
2016: T2=nf_get_T2(base,roots(x,prec));
2017: }
2018: }
2019:
2020: /* x can be a polynomial, but also an nf or a bnf */
2021: /* if check != NULL, return the first polynomial pol
2022: found such that check(arg, pol) != 0; return NULL
2023: if there are no such pol */
2024: static GEN
2025: allpolred0(GEN x, GEN *pta, long code, long prec,
2026: int (*check)(GEN,GEN), GEN arg)
2027: {
2028: GEN y,p1, base = NULL, polr = NULL;
2029: long av = avma;
2030:
2031: if (typ(x) == t_POL)
2032: {
2033: if (!signe(x)) return gcopy(x);
2034: check_pol_int(x,"polred");
2035: if (!gcmp1(leading_term(x)))
2036: err(impl,"allpolred for nonmonic polynomials");
2037: base = allbase4(x,code,&p1,NULL); /* p1 is junk */
2038: }
2039: else
2040: {
2041: long i = lg(x);
2042: if (typ(x) == t_VEC && i<=4 && i>=3 && typ(x[1])==t_POL)
2043: { /* polynomial + integer basis */
2044: base=(GEN)x[2]; x=(GEN)x[1];
2045: }
2046: else
2047: {
2048: x = checknf(x);
2049: base=(GEN)x[7]; x=(GEN)x[1];
2050: }
2051: }
2052: p1 = LLL_nfbasis(&x,polr,base,prec);
2053: y = pols_for_polred(x,base,p1,pta,check,arg);
2054: if (check)
2055: {
2056: if (y) return gerepileupto(av, y);
2057: avma = av; return NULL;
2058: }
2059:
2060: if (pta)
2061: {
2062: GEN *gptr[2]; gptr[0]=&y; gptr[1]=pta;
2063: gerepilemany(av,gptr,2); return y;
2064: }
2065: return gerepileupto(av,y);
2066: }
2067:
2068: static GEN
2069: allpolred(GEN x, GEN *pta, long code, long prec)
2070: {
2071: return allpolred0(x,pta,code,prec,NULL,NULL);
2072: }
2073:
2074: GEN
2075: polred0(GEN x,long flag, GEN p, long prec)
2076: {
2077: GEN y;
2078: long smll = (flag & 1);
2079:
2080: if (p && !gcmp0(p)) smll=(long) p; /* factored polred */
2081: if (flag & 2) /* polred2 */
2082: {
2083: y=cgetg(3,t_MAT);
2084: y[2]=(long)allpolred(x,(GEN*)(y+1),smll,prec);
2085: return y;
2086: }
2087: return allpolred(x,NULL,smll,prec);
2088: }
2089:
2090: GEN
2091: ordred(GEN x, long prec)
2092: {
2093: GEN base,y;
2094: long n=degpol(x),i,av=avma,v = varn(x);
2095:
2096: if (typ(x) != t_POL) err(typeer,"ordred");
2097: if (!signe(x)) return gcopy(x);
2098: if (!gcmp1((GEN)x[n+2])) err(impl,"ordred for nonmonic polynomials");
2099:
2100: base = cgetg(n+1,t_VEC); /* power basis */
2101: for (i=1; i<=n; i++)
2102: base[i] = lpowgs(polx[v],i-1);
2103: y = cgetg(3,t_VEC);
2104: y[1] = (long)x;
2105: y[2] = (long)base;
2106: return gerepileupto(av, allpolred(y,NULL,0,prec));
2107: }
2108:
2109: GEN
2110: sum(GEN v, long a, long b)
2111: {
2112: GEN p;
2113: long i;
2114: if (a > b) return gzero;
2115: p = (GEN)v[a];
2116: for (i=a+1; i<=b; i++) p = gadd(p, (GEN)v[i]);
2117: return p;
2118: }
2119:
2120: GEN
2121: T2_from_embed_norm(GEN x, long r1)
2122: {
2123: GEN p = sum(x, 1, r1);
2124: GEN q = sum(x, r1+1, lg(x)-1);
2125: if (q != gzero) p = gadd(p, gmul2n(q,1));
2126: return p;
2127: }
2128:
2129: GEN
2130: T2_from_embed(GEN x, long r1)
2131: {
2132: return T2_from_embed_norm(gnorm(x), r1);
2133: }
2134:
2135: /* return T2 norm of the polynomial defining nf */
2136: static GEN
2137: get_Bnf(GEN nf)
2138: {
2139: return T2_from_embed((GEN)nf[6], nf_get_r1(nf));
2140: }
2141:
2142: /* characteristic pol of x */
2143: static GEN
2144: get_polchar(GEN data, GEN x)
2145: {
2146: GEN basis_embed = (GEN)data[1];
2147: GEN g = gmul(basis_embed,x);
2148: return ground(roots_to_pol_r1r2(g, data[0], 0));
2149: }
2150:
2151: /* return a defining polynomial for Q(x) */
2152: static GEN
2153: get_polmin(GEN data, GEN x)
2154: {
2155: GEN g = get_polchar(data,x);
2156: GEN h = modulargcd(g,derivpol(g));
2157: if (lgef(h) > 3) g = gdivexact(g,h);
2158: return g;
2159: }
2160:
2161: /* does x generate the correct field ? */
2162: static GEN
2163: chk_gen(GEN data, GEN x)
2164: {
2165: long av = avma;
2166: GEN g = get_polchar(data,x);
2167: if (lgef(modulargcd(g,derivpol(g))) > 3) { avma=av; return NULL; }
2168: if (DEBUGLEVEL>3) fprintferr(" generator: %Z\n",g);
2169: return g;
2170: }
2171:
2172: /* mat = base change matrix, gram = LLL-reduced T2 matrix */
2173: static GEN
2174: chk_gen_init(FP_chk_fun *chk, GEN nf, GEN gram, GEN mat, long *ptprec)
2175: {
2176: GEN P,bound,prev,x,B,data, M = gmael(nf,5,1);
2177: long N = lg(nf[7]), n = N-1,i,prec,prec2;
2178: int skipfirst = 1; /* [1,0...0] --> x rational */
2179:
2180: /* data[0] = r1
2181: * data[1] = embeddings of LLL-reduced Zk basis
2182: * data[2] = LLL reduced Zk basis (in M_n(Z)) */
2183: data = new_chunk(3);
2184: data[0] = itos(gmael(nf,2,1));
2185: data[1] = lmul(M, mat);
2186: data[2] = lmul((GEN)nf[7],mat);
2187: chk->data = data;
2188:
2189: x = cgetg(N,t_COL);
2190: bound = get_Bnf(nf); prev = NULL;
2191: for (i=1; i<N; i++) x[i]=zero;
2192: for (i=2; i<N; i++)
2193: {
2194: x[i] = un;
2195: P = get_polmin(data,x);
2196: x[i] = zero;
2197: if (degpol(P) == n)
2198: {
2199: B = gcoeff(gram,i,i);
2200: if (gcmp(B,bound) < 0) bound = B ;
2201: }
2202: else
2203: {
2204: if (DEBUGLEVEL>2) fprintferr("chk_gen_init: subfield %Z\n",P);
2205: if (skipfirst == i-1)
2206: {
2207: if (prev && !gegal(prev,P))
2208: {
2209: if (degpol(prev) * degpol(P) > 32) continue; /* too expensive */
2210: P = (GEN)compositum(prev,P)[1];
2211: if (degpol(P) == n) continue;
2212: if (DEBUGLEVEL>2 && lgef(P)>lgef(prev))
2213: fprintferr("chk_gen_init: subfield %Z\n",P);
2214: }
2215: skipfirst++; prev = P;
2216: }
2217: }
2218: }
2219: /* x_1,...,x_skipfirst generate a strict subfield [unless n=skipfirst=1] */
2220: chk->skipfirst = skipfirst;
2221: if (DEBUGLEVEL>2) fprintferr("chk_gen_init: skipfirst = %ld\n",skipfirst);
2222:
2223: /* should be gexpo( [max_k C^n_k (bound/k) ^ k] ^ (1/2) ) */
2224: prec2 = (1 + (((gexpo(bound)*n)/2) >> TWOPOTBITS_IN_LONG));
2225: if (prec2 < 0) prec2 = 0;
2226: prec = 3 + prec2;
2227: prec2= nfgetprec(nf);
2228: if (DEBUGLEVEL)
2229: fprintferr("chk_gen_init: estimated prec = %ld (initially %ld)\n",
2230: prec, prec2);
2231: if (prec > prec2) return NULL;
2232: if (prec < prec2) data[1] = (long)gprec_w((GEN)data[1], prec);
2233: *ptprec = prec; return bound;
2234: }
2235:
2236: static GEN
2237: chk_gen_post(GEN data, GEN res)
2238: {
2239: GEN basis = (GEN)data[2];
2240: GEN p1 = (GEN)res[2];
2241: long i, l = lg(p1);
2242: for (i=1; i<l; i++) p1[i] = lmul(basis, (GEN)p1[i]);
2243: return res;
2244: }
2245:
2246: /* no garbage collecting, done in polredabs0 */
2247: static GEN
2248: findmindisc(GEN nf, GEN y, GEN a, GEN phimax, long flun)
2249: {
2250: long i,k, c = lg(y);
2251: GEN v,dmin,z,beta,discs = cgetg(c,t_VEC);
2252:
2253: for (i=1; i<c; i++) discs[i] = labsi(discsr((GEN)y[i]));
2254: v = sindexsort(discs);
2255: k = v[1]; dmin = (GEN)discs[k]; z = (GEN)y[k]; beta = (GEN)a[k];
2256: for (i=2; i<c; i++)
2257: {
2258: k = v[i];
2259: if (!egalii((GEN)discs[k],dmin)) break;
2260: if (gpolcomp((GEN)y[k],z) < 0) { z = (GEN)y[k]; beta = (GEN)a[k]; }
2261: }
2262: if (flun & nf_RAW)
2263: {
2264: y=cgetg(3,t_VEC);
2265: y[1]=lcopy(z);
2266: y[2]=lcopy(beta);
2267: }
2268: else if (phimax)
2269: {
2270: y=cgetg(3,t_VEC);
2271: y[1]=lcopy(z);
2272: beta=polymodrecip(gmodulcp(beta,(GEN)nf[1]));
2273: y[2]=(long)poleval(phimax,beta);
2274: }
2275: else y = gcopy(z);
2276: return y;
2277: }
2278:
2279: /* no garbage collecting, done in polredabs0 */
2280: static GEN
2281: storeallpols(GEN nf, GEN z, GEN a, GEN phimax, long flun)
2282: {
2283: GEN p1,y,beta;
2284:
2285: if (flun & nf_RAW)
2286: {
2287: long i, c = lg(z);
2288: y=cgetg(c,t_VEC);
2289: for (i=1; i<c; i++)
2290: {
2291: p1=cgetg(3,t_VEC); y[i]=(long)p1;
2292: p1[1]=lcopy((GEN)z[i]);
2293: p1[2]=lcopy((GEN)a[i]);
2294: }
2295: }
2296: else if (phimax)
2297: {
2298: long i, c = lg(z);
2299: beta = new_chunk(c);
2300: for (i=1; i<c; i++)
2301: beta[i] = (long)polymodrecip(gmodulcp((GEN)a[i],(GEN)nf[1]));
2302:
2303: y=cgetg(c,t_VEC);
2304: for (i=1; i<c; i++)
2305: {
2306: p1=cgetg(3,t_VEC); y[i]=(long)p1;
2307: p1[1]=lcopy((GEN)z[i]);
2308: p1[2]=(long)poleval(phimax,(GEN)beta[i]);
2309: }
2310: }
2311: else y = gcopy(z);
2312: return y;
2313: }
2314:
2315: GEN
2316: polredabs0(GEN x, long flun, long prec)
2317: {
2318: long i,nv, av = avma;
2319: GEN nf,v,y,a,phimax;
2320: GEN (*storepols)(GEN, GEN, GEN, GEN, long);
2321: FP_chk_fun *chk;
2322:
2323: chk = (FP_chk_fun*)new_chunk(sizeof(FP_chk_fun));
2324: chk->f = &chk_gen;
2325: chk->f_init = &chk_gen_init;
2326: chk->f_post = &chk_gen_post;
2327:
2328: if ((ulong)flun >= 16) err(flagerr,"polredabs");
2329: nf = initalgall0(x,nf_SMALL|nf_REGULAR,prec);
2330: if (lg(nf) == 3)
2331: {
2332: phimax = lift_to_pol((GEN)nf[2]);
2333: nf = (GEN)nf[1];
2334: }
2335: else
2336: phimax = (flun & nf_ORIG)? polx[0]: (GEN)NULL;
2337: prec = nfgetprec(nf);
2338: x = (GEN)nf[1];
2339:
2340: if (lgef(x) == 4)
2341: {
2342: y = _vec(polx[varn(x)]);
2343: a = _vec(gsub((GEN)y[1], (GEN)x[2]));
2344: }
2345: else
2346: {
2347: for (i=1; ; i++)
2348: {
2349: v = fincke_pohst(nf,NULL,5000,3,prec, chk);
2350: if (v) break;
2351: if (i==MAXITERPOL) err(accurer,"polredabs0");
2352: prec = (prec<<1)-2; nf = nfnewprec(nf,prec);
2353: if (DEBUGLEVEL) err(warnprec,"polredabs0",prec);
2354: }
2355: a = (GEN)v[2];
2356: y = (GEN)v[1];
2357: }
2358: nv = lg(a);
2359: for (i=1; i<nv; i++)
2360: if (canon_pol((GEN)y[i]) < 0 && phimax)
2361: a[i] = (long) gneg_i((GEN)a[i]);
2362: nv = remove_duplicates(y,a);
2363:
2364: if (DEBUGLEVEL)
2365: { fprintferr("%ld minimal vectors found.\n",nv-1); flusherr(); }
2366: if (nv >= 10000) flun &= (~nf_ALL); /* should not happen */
2367: storepols = (flun & nf_ALL)? storeallpols: findmindisc;
2368:
2369: if (DEBUGLEVEL) fprintferr("\n");
2370: if (nv==1)
2371: {
2372: y = _vec(x);
2373: a = _vec(polx[varn(x)]);
2374: }
2375: if (varn(y[1]) != varn(x))
2376: {
2377: long vx = varn(x);
2378: for (i=1; i<nv; i++) setvarn(y[i], vx);
2379: }
2380: return gerepileupto(av, storepols(nf,y,a,phimax,flun));
2381: }
2382:
2383: GEN
2384: polredabsall(GEN x, long flun, long prec)
2385: {
2386: return polredabs0(x, flun | nf_ALL, prec);
2387: }
2388:
2389: GEN
2390: polredabs(GEN x, long prec)
2391: {
2392: return polredabs0(x,nf_REGULAR,prec);
2393: }
2394:
2395: GEN
2396: polredabs2(GEN x, long prec)
2397: {
2398: return polredabs0(x,nf_ORIG,prec);
2399: }
2400:
2401: GEN
2402: polredabsnored(GEN x, long prec)
2403: {
2404: return polredabs0(x,nf_NORED,prec);
2405: }
2406:
2407: GEN
2408: polred(GEN x, long prec)
2409: {
2410: return allpolred(x,(GEN*)0,0,prec);
2411: }
2412:
2413: GEN
2414: polredfirstpol(GEN x, long prec, int (*check)(GEN,GEN), GEN arg)
2415: {
2416: return allpolred0(x,(GEN*)0,0,prec,check,arg);
2417: }
2418:
2419: GEN
2420: smallpolred(GEN x, long prec)
2421: {
2422: return allpolred(x,(GEN*)0,1,prec);
2423: }
2424:
2425: GEN
2426: factoredpolred(GEN x, GEN p, long prec)
2427: {
2428: return allpolred(x,(GEN*)0,(long)p,prec);
2429: }
2430:
2431: GEN
2432: polred2(GEN x, long prec)
2433: {
2434: GEN y=cgetg(3,t_MAT);
2435:
2436: y[2]= (long) allpolred(x,(GEN*)(y+1),0,prec);
2437: return y;
2438: }
2439:
2440: GEN
2441: smallpolred2(GEN x, long prec)
2442: {
2443: GEN y=cgetg(3,t_MAT);
2444:
2445: y[2]= (long) allpolred(x,(GEN*)(y+1),1,prec);
2446: return y;
2447: }
2448:
2449: GEN
2450: factoredpolred2(GEN x, GEN p, long prec)
2451: {
2452: GEN y=cgetg(3,t_MAT);
2453:
2454: y[2]= (long) allpolred(x,(GEN*)(y+1),(long)p,prec);
2455: return y;
2456: }
2457:
2458: /* relative polredabs. Returns
2459: * flag = 0: relative polynomial
2460: * flag = 1: relative polynomial + element
2461: * flag = 2: absolute polynomial */
2462: GEN
2463: rnfpolredabs(GEN nf, GEN relpol, long flag, long prec)
2464: {
2465: GEN p1,bpol,rnf,elt,pol;
2466: long va, av = avma;
2467:
2468: if (typ(relpol)!=t_POL) err(typeer,"rnfpolredabs");
2469: nf=checknf(nf); va = varn(relpol);
2470: if (DEBUGLEVEL>1) timer2();
2471: p1 = makebasis(nf, unifpol(nf,relpol,1));
2472: rnf = (GEN)p1[3];
2473: if (DEBUGLEVEL>1)
2474: {
2475: msgtimer("absolute basis");
2476: fprintferr("original absolute generator: %Z\n",p1[1]);
2477: }
2478: p1 = polredabs0(p1, nf_RAW | nf_NORED, prec);
2479: bpol = (GEN)p1[1];
2480: if (DEBUGLEVEL>1) fprintferr("reduced absolute generator: %Z\n",bpol);
2481: if (flag==2) return gerepileupto(av,bpol);
2482:
2483: elt = rnfelementabstorel(rnf,(GEN)p1[2]);
2484: p1=cgetg(3,t_VEC);
2485: pol = rnfcharpoly(nf,relpol,elt,va);
2486: if (!flag) p1 = pol;
2487: else
2488: {
2489: p1[1]=(long)pol;
2490: p1[2]=(long)polymodrecip(elt);
2491: }
2492: return gerepileupto(av,p1);
2493: }
2494:
2495: /********************************************************************/
2496: /** **/
2497: /** MINIM **/
2498: /** **/
2499: /********************************************************************/
2500: int addcolumntomatrix(GEN V,GEN INVP,GEN L);
2501:
2502: /* x is a non-empty matrix, y is a compatible VECSMALL (dimension > 0). */
2503: GEN
2504: gmul_mat_smallvec(GEN x, GEN y)
2505: {
2506: long c = lg(x), l = lg(x[1]), av,i,j;
2507: GEN z=cgetg(l,t_COL), s;
2508:
2509: for (i=1; i<l; i++)
2510: {
2511: av = avma; s = gmulgs(gcoeff(x,i,1),y[1]);
2512: for (j=2; j<c; j++)
2513: if (y[j]) s = gadd(s, gmulgs(gcoeff(x,i,j),y[j]));
2514: z[i] = lpileupto(av,s);
2515: }
2516: return z;
2517: }
2518:
2519: /* same, x integral */
2520: GEN
2521: gmul_mati_smallvec(GEN x, GEN y)
2522: {
2523: long c = lg(x), l = lg(x[1]), av,i,j;
2524: GEN z=cgetg(l,t_COL), s;
2525:
2526: for (i=1; i<l; i++)
2527: {
2528: av = avma; s = mulis(gcoeff(x,i,1),y[1]);
2529: for (j=2; j<c; j++)
2530: if (y[j]) s = addii(s, mulis(gcoeff(x,i,j),y[j]));
2531: z[i]=(long)gerepileuptoint(av,s);
2532: }
2533: return z;
2534: }
2535:
2536: void
2537: minim_alloc(long n, double ***q, GEN *x, double **y, double **z, double **v)
2538: {
2539: long i, s;
2540: double **Q;
2541:
2542: *x = cgetg(n, t_VECSMALL);
2543: *q = (double**) new_chunk(n);
2544:
2545: /* correct alignment for the following */
2546: s = avma % sizeof(double); avma -= s;
2547: if (avma<bot) err(errpile);
2548:
2549: s = (n * sizeof(double))/sizeof(long);
2550: *y = (double*) new_chunk(s);
2551: *z = (double*) new_chunk(s);
2552: *v = (double*) new_chunk(s); Q = *q;
2553: for (i=1; i<n; i++) Q[i] = (double*) new_chunk(s);
2554: }
2555:
2556: /* Minimal vectors for the integral definite quadratic form: a.
2557: * Result u:
2558: * u[1]= Number of vectors of square norm <= BORNE
2559: * u[2]= maximum norm found
2560: * u[3]= list of vectors found (at most STOCKMAX)
2561: *
2562: * If BORNE = gzero: Minimal non-zero vectors.
2563: * flag = min_ALL, as above
2564: * flag = min_FIRST, exits when first suitable vector is found.
2565: * flag = min_PERF, only compute rank of the family of v.v~ (v min.)
2566: */
2567: static GEN
2568: minim00(GEN a, GEN BORNE, GEN STOCKMAX, long flag)
2569: {
2570: GEN x,res,p1,u,r,liste,gnorme,invp,V, *gptr[2];
2571: long n = lg(a), av0 = avma, av1,av,tetpil,lim, i,j,k,s,maxrank;
2572: double p,maxnorm,BOUND,*v,*y,*z,**q, eps = 0.000001;
2573:
2574: maxrank = 0; res = V = invp = NULL; /* gcc -Wall */
2575: switch(flag)
2576: {
2577: case min_FIRST:
2578: if (gcmp0(BORNE)) err(talker,"bound = 0 in minim2");
2579: res = cgetg(3,t_VEC); break;
2580: case min_ALL: res = cgetg(4,t_VEC); break;
2581: case min_PERF: break;
2582: default: err(talker, "incorrect flag in minim00");
2583: }
2584: if (n == 1)
2585: {
2586: switch(flag)
2587: {
2588: case min_FIRST: avma=av0; return cgetg(1,t_VEC);
2589: case min_PERF: avma=av0; return gzero;
2590: }
2591: res[1]=res[2]=zero;
2592: res[3]=lgetg(1,t_MAT); return res;
2593: }
2594:
2595: av = avma;
2596: minim_alloc(n, &q, &x, &y, &z, &v);
2597: av1 = avma;
2598:
2599: u = lllgramint(a);
2600: if (lg(u) != n)
2601: err(talker,"not a definite form in minim00");
2602: a = qf_base_change(a,u,1);
2603:
2604: n--;
2605: a = gmul(a, realun(DEFAULTPREC)); r = sqred1(a);
2606: if (DEBUGLEVEL>4) { fprintferr("minim: r = "); outerr(r); }
2607: for (j=1; j<=n; j++)
2608: {
2609: v[j] = rtodbl(gcoeff(r,j,j));
2610: for (i=1; i<j; i++) q[i][j] = rtodbl(gcoeff(r,i,j));
2611: }
2612:
2613: if (flag==min_PERF || gcmp0(BORNE))
2614: {
2615: double c, b = rtodbl(gcoeff(a,1,1));
2616:
2617: for (i=2; i<=n; i++)
2618: { c=rtodbl(gcoeff(a,i,i)); if (c<b) b=c; }
2619: BOUND = b+eps;
2620: BORNE = ground(dbltor(BOUND));
2621: maxnorm = -1.; /* don't update maxnorm */
2622: }
2623: else
2624: {
2625: BORNE = gfloor(BORNE);
2626: BOUND = gtodouble(BORNE)+eps;
2627: maxnorm = 0.;
2628: }
2629:
2630: switch(flag)
2631: {
2632: case min_ALL:
2633: maxrank=itos(STOCKMAX);
2634: liste = new_chunk(1+maxrank);
2635: break;
2636: case min_PERF:
2637: BORNE = gerepileupto(av1,BORNE);
2638: maxrank = (n*(n+1))>>1;
2639: liste = cgetg(1+maxrank, t_VECSMALL);
2640: V = cgetg(1+maxrank, t_VECSMALL);
2641: for (i=1; i<=maxrank; i++) liste[i]=0;
2642: }
2643:
2644: s=0; av1=avma; lim = stack_lim(av1,1);
2645: k = n; y[n] = z[n] = 0;
2646: x[n] = (long) sqrt(BOUND/v[n]);
2647: if (flag == min_PERF) invp = idmat(maxrank);
2648: for(;;x[1]--)
2649: {
2650: do
2651: {
2652: if (k>1)
2653: {
2654: long l = k-1;
2655: z[l]=0;
2656: for (j=k; j<=n; j++) z[l] += q[l][j]*x[j];
2657: p = x[k]+z[k];
2658: y[l] = y[k] + p*p*v[k];
2659: x[l] = (long) floor(sqrt((BOUND-y[l])/v[l])-z[l]);
2660: k = l;
2661: }
2662: for(;;)
2663: {
2664: p = x[k]+z[k];
2665: if (y[k] + p*p*v[k] <= BOUND) break;
2666: k++; x[k]--;
2667: }
2668: }
2669: while (k>1);
2670: if (! x[1] && y[1]<=eps) break;
2671: p = x[1]+z[1]; p = y[1] + p*p*v[1]; /* norm(x) */
2672: if (maxnorm >= 0)
2673: {
2674: if (flag == min_FIRST)
2675: {
2676: gnorme = dbltor(p);
2677: tetpil=avma; gnorme = ground(gnorme); r = gmul_mati_smallvec(u,x);
2678: gptr[0]=&gnorme; gptr[1]=&r; gerepilemanysp(av,tetpil,gptr,2);
2679: res[1]=(long)gnorme;
2680: res[2]=(long)r; return res;
2681: }
2682: if (p > maxnorm) maxnorm = p;
2683: }
2684: else
2685: {
2686: long av2 = avma;
2687: gnorme = ground(dbltor(p));
2688: if (gcmp(gnorme,BORNE) >= 0) avma = av2;
2689: else
2690: {
2691: BOUND=gtodouble(gnorme)+eps; s=0;
2692: affii(gnorme,BORNE); avma=av1;
2693: if (flag == min_PERF) invp = idmat(maxrank);
2694: }
2695: }
2696: s++;
2697:
2698: switch(flag)
2699: {
2700: case min_ALL:
2701: if (s<=maxrank)
2702: {
2703: p1 = new_chunk(n+1); liste[s] = (long) p1;
2704: for (i=1; i<=n; i++) p1[i] = x[i];
2705: }
2706: break;
2707:
2708: case min_PERF:
2709: {
2710: long av2=avma, I=1;
2711:
2712: for (i=1; i<=n; i++)
2713: for (j=i; j<=n; j++,I++) V[I] = x[i]*x[j];
2714: if (! addcolumntomatrix(V,invp,liste))
2715: {
2716: if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }
2717: s--; avma=av2; continue;
2718: }
2719:
2720: if (DEBUGLEVEL>1) { fprintferr("*"); flusherr(); }
2721: if (s == maxrank)
2722: {
2723: if (DEBUGLEVEL>1) { fprintferr("\n"); flusherr(); }
2724: avma=av0; return stoi(s);
2725: }
2726:
2727: if (low_stack(lim, stack_lim(av1,1)))
2728: {
2729: if(DEBUGMEM>1) err(warnmem,"minim00, rank>=%ld",s);
2730: invp = gerepilecopy(av1, invp);
2731: }
2732: }
2733: }
2734: }
2735: switch(flag)
2736: {
2737: case min_FIRST:
2738: avma=av0; return cgetg(1,t_VEC);
2739: case min_PERF:
2740: if (DEBUGLEVEL>1) { fprintferr("\n"); flusherr(); }
2741: avma=av0; return stoi(s);
2742: }
2743: k = min(s,maxrank);
2744:
2745: tetpil = avma; p1=cgetg(k+1,t_MAT);
2746: for (j=1; j<=k; j++)
2747: p1[j] = (long) gmul_mati_smallvec(u,(GEN)liste[j]);
2748: liste = p1;
2749: r = (maxnorm >= 0) ? ground(dbltor(maxnorm)): BORNE;
2750:
2751: r=icopy(r); gptr[0]=&r; gptr[1]=&liste;
2752: gerepilemanysp(av,tetpil,gptr,2);
2753: res[1]=lstoi(s<<1);
2754: res[2]=(long)r;
2755: res[3]=(long)liste; return res;
2756: }
2757:
2758: GEN
2759: qfminim0(GEN a, GEN borne, GEN stockmax, long flag, long prec)
2760: {
2761: switch(flag)
2762: {
2763: case 0: return minim00(a,borne,stockmax,min_ALL);
2764: case 1: return minim00(a,borne,gzero ,min_FIRST);
2765: case 2: return fincke_pohst(a,borne,itos(stockmax),0,prec,NULL);
2766: default: err(flagerr,"qfminim");
2767: }
2768: return NULL; /* not reached */
2769: }
2770:
2771: GEN
2772: minim(GEN a, GEN borne, GEN stockmax)
2773: {
2774: return minim00(a,borne,stockmax,min_ALL);
2775: }
2776:
2777: GEN
2778: minim2(GEN a, GEN borne, GEN stockmax)
2779: {
2780: return minim00(a,borne,stockmax,min_FIRST);
2781: }
2782:
2783: GEN
2784: perf(GEN a)
2785: {
2786: return minim00(a,gzero,gzero,min_PERF);
2787: }
2788:
2789: /* general program for positive definit quadratic forms (real coeffs).
2790: * One needs BORNE != 0; LLL reduction done in fincke_pohst().
2791: * If (flag & 2) stop as soon as stockmax is reached.
2792: * If (flag & 1) return NULL on precision problems (no error).
2793: * If (check != NULL consider only vectors passing the check [ assumes
2794: * stockmax > 0 and we only want the smallest possible vectors ] */
2795: static GEN
2796: smallvectors(GEN a, GEN BORNE, long stockmax, long flag, FP_chk_fun *CHECK)
2797: {
2798: long av,av1,lim,N,n,i,j,k,s,epsbit,prec, checkcnt = 1;
2799: GEN u,S,x,y,z,v,q,norme1,normax1,borne1,borne2,eps,p1,alpha,norms,dummy;
2800: GEN (*check)(GEN,GEN) = CHECK? CHECK->f: NULL;
2801: GEN data = CHECK? CHECK->data: NULL;
2802: int skipfirst = CHECK? CHECK->skipfirst: 0;
2803:
2804: if (DEBUGLEVEL)
2805: fprintferr("smallvectors looking for norm <= %Z\n",gprec_w(BORNE,3));
2806:
2807: q = sqred1intern(a,flag&1);
2808: if (q == NULL) return NULL;
2809: if (DEBUGLEVEL>5) fprintferr("q = %Z",q);
2810: prec = gprecision(q);
2811: epsbit = bit_accuracy(prec) >> 1;
2812: eps = realun(prec); setexpo(eps,-epsbit);
2813: alpha = dbltor(0.95);
2814: normax1 = gzero;
2815: borne1= gadd(BORNE,eps);
2816: borne2 = mpmul(borne1,alpha);
2817: N = lg(q); n = N-1;
2818: v = cgetg(N,t_VEC);
2819: dummy = cgetg(1,t_STR);
2820:
2821: av = avma; lim = stack_lim(av,1);
2822: if (check) norms = cgetg(stockmax+1,t_VEC);
2823: S = cgetg(stockmax+1,t_VEC);
2824: x = cgetg(N,t_COL);
2825: y = cgetg(N,t_COL);
2826: z = cgetg(N,t_COL);
2827: for (i=1; i<N; i++) { v[i] = coeff(q,i,i); x[i]=y[i]=z[i] = zero; }
2828:
2829: x[n] = lmpent(mpsqrt(gdiv(borne1,(GEN)v[n])));
2830: if (DEBUGLEVEL>3) { fprintferr("\nx[%ld] = %Z\n",n,x[n]); flusherr(); }
2831:
2832: s = 0; k = n;
2833: for(;; x[k] = laddis((GEN)x[k],-1)) /* main */
2834: {
2835: do
2836: {
2837: int fl = 0;
2838: if (k > 1)
2839: {
2840: av1=avma; k--; p1 = mpmul(gcoeff(q,k,k+1),(GEN)x[k+1]);
2841: for (j=k+2; j<N; j++)
2842: p1 = mpadd(p1, mpmul(gcoeff(q,k,j),(GEN)x[j]));
2843: z[k] = (long)gerepileuptoleaf(av1,p1);
2844:
2845: av1=avma; p1 = gsqr(mpadd((GEN)x[k+1],(GEN)z[k+1]));
2846: p1 = mpadd((GEN)y[k+1], mpmul(p1,(GEN)v[k+1]));
2847: y[k] = (long)gerepileuptoleaf(av1, p1);
2848:
2849: av1=avma; p1 = mpsub(borne1, (GEN)y[k]);
2850: if (signe(p1) < 0) { avma=av1; fl = 1; }
2851: else
2852: {
2853: p1 = mpadd(eps,mpsub(mpsqrt(gdiv(p1,(GEN)v[k])), (GEN)z[k]));
2854: x[k] = (long)gerepileuptoleaf(av1,mpent(p1));
2855: }
2856: /* reject the [x_1,...,x_skipfirst,0,...,0] */
2857: if (k <= skipfirst && !signe(y[skipfirst])) goto END;
2858: }
2859: for(;; x[k] = laddis((GEN)x[k],-1))
2860: {
2861: if (!fl)
2862: {
2863: av1 = avma; /* p1 >= 0 */
2864: p1 = mpmul((GEN)v[k], gsqr(mpadd((GEN)x[k],(GEN)z[k])));
2865: i = mpcmp(mpsub(mpadd(p1,(GEN)y[k]), borne1), gmul2n(p1,-epsbit));
2866: avma = av1; if (i <= 0) break;
2867: }
2868: k++; fl=0;
2869: }
2870:
2871: if (low_stack(lim, stack_lim(av,1)))
2872: {
2873: GEN *gptr[7];
2874: int cnt = 4;
2875: if(DEBUGMEM>1) err(warnmem,"smallvectors");
2876: gptr[0]=&x; gptr[1]=&y; gptr[2]=&z; gptr[3]=&normax1;
2877: if (stockmax)
2878: { /* initialize to dummy values */
2879: GEN T = S;
2880: for (i=s+1; i<=stockmax; i++) S[i]=(long)dummy;
2881: S = gclone(S); if (isclone(T)) gunclone(T);
2882: }
2883: if (check)
2884: {
2885: cnt+=3; gptr[4]=&borne1; gptr[5]=&borne2; gptr[6]=&norms;
2886: for (i=s+1; i<=stockmax; i++) norms[i]=(long)dummy;
2887: }
2888: gerepilemany(av,gptr,cnt);
2889: }
2890: if (DEBUGLEVEL>3)
2891: {
2892: if (DEBUGLEVEL>5) fprintferr("%ld ",k);
2893: if (k==n) fprintferr("\nx[%ld] = %Z\n",n,x[n]);
2894: flusherr();
2895: }
2896: }
2897: while (k > 1);
2898:
2899: /* x = 0: we're done */
2900: if (!signe(x[1]) && !signe(y[1])) goto END;
2901:
2902: av1=avma; p1 = gsqr(mpadd((GEN)x[1],(GEN)z[1]));
2903: norme1 = mpadd((GEN)y[1], mpmul(p1, (GEN)v[1]));
2904: if (mpcmp(norme1,borne1) > 0) { avma=av1; continue; /* main */ }
2905:
2906: norme1 = gerepileupto(av1,norme1);
2907: if (check)
2908: {
2909: if (checkcnt < 5 && mpcmp(norme1, borne2) < 0)
2910: {
2911: if (check(data,x))
2912: {
2913: borne1 = mpadd(norme1,eps);
2914: borne2 = mpmul(borne1,alpha);
2915: s = 0; checkcnt = 0;
2916: }
2917: else { checkcnt++ ; continue; /* main */}
2918: }
2919: }
2920: else if (mpcmp(norme1,normax1) > 0) normax1 = norme1;
2921: if (++s <= stockmax)
2922: {
2923: if (check) norms[s] = (long)norme1;
2924: S[s] = (long)dummycopy(x);
2925: if (s == stockmax && (flag&2) && check)
2926: {
2927: long av1 = avma;
2928: GEN per = sindexsort(norms);
2929: if (DEBUGLEVEL) fprintferr("sorting...\n");
2930: for (i=1; i<=s; i++)
2931: {
2932: long k = per[i];
2933: if (check(data,(GEN)S[k]))
2934: {
2935: S[1] = S[k]; avma = av1;
2936: borne1 = mpadd(norme1,eps);
2937: borne2 = mpmul(borne1,alpha);
2938: s = 1; checkcnt = 0; break;
2939: }
2940: }
2941: if (checkcnt) s = 0;
2942: }
2943: }
2944: }
2945: END:
2946: if (s < stockmax) stockmax = s;
2947: if (check)
2948: {
2949: GEN per, alph,pols,p;
2950: if (DEBUGLEVEL) fprintferr("final sort & check...\n");
2951: setlg(norms,s+1); per = sindexsort(norms);
2952: alph = cgetg(s+1,t_VEC);
2953: pols = cgetg(s+1,t_VEC);
2954: for (j=0,i=1; i<=s; i++)
2955: {
2956: long k = per[i];
2957: if (j && mpcmp((GEN)norms[k], borne1) > 0) break;
2958: if ((p = check(data,(GEN)S[k])))
2959: {
2960: if (!j) borne1 = gadd((GEN)norms[k],eps);
2961: j++; pols[j]=(long)p; alph[j]=S[k];
2962: }
2963: }
2964: u = cgetg(3,t_VEC);
2965: setlg(pols,j+1); u[1] = (long)pols;
2966: setlg(alph,j+1); u[2] = (long)alph;
2967: if (isclone(S)) { u[2] = (long)forcecopy(alph); gunclone(S); }
2968: return u;
2969: }
2970: u = cgetg(4,t_VEC);
2971: u[1] = lstoi(s<<1);
2972: u[2] = (long)normax1;
2973: if (stockmax)
2974: {
2975: setlg(S,stockmax+1);
2976: settyp(S,t_MAT);
2977: if (isclone(S)) { p1 = S; S = forcecopy(S); gunclone(p1); }
2978: }
2979: else
2980: S = cgetg(1,t_MAT);
2981: u[3] = (long)S; return u;
2982: }
2983:
2984: /* solve x~.a.x <= bound, a > 0. If check is non-NULL keep x only if check(x).
2985: * flag & 1, return NULL in case of precision problems (sqred1 or lllgram)
2986: * raise an error otherwise.
2987: * flag & 2, return as soon as stockmax vectors found.
2988: * If a is a number field, use its T2 matrix
2989: */
2990: GEN
2991: fincke_pohst(GEN a,GEN B0,long stockmax,long flag, long PREC, FP_chk_fun *CHECK)
2992: {
2993: VOLATILE long pr,av=avma,i,j,n;
2994: VOLATILE GEN B,nf,r,rinvtrans,u,v,res,z,vnorm,sperm,perm,uperm,gram;
2995: VOLATILE GEN bound = B0;
2996: void *catcherr = NULL;
2997: long prec = PREC;
2998:
2999: if (DEBUGLEVEL>2) { fprintferr("entering fincke_pohst\n"); flusherr(); }
3000: if (typ(a) == t_VEC) { nf = checknf(a); a = gmael(nf,5,3); } else nf = NULL;
3001: pr = gprecision(a);
3002: if (pr) prec = pr; else a = gmul(a,realun(prec));
3003: if (DEBUGLEVEL>2) fprintferr("first LLL: prec = %ld\n", prec);
3004: if (nf && !signe(gmael(nf,2,2))) /* totally real */
3005: {
3006: GEN T = nf_get_T((GEN)nf[1], (GEN)nf[7]);
3007: u = lllgramint(T);
3008: prec += 2 * ((gexpo(u) + gexpo((GEN)nf[1])) >> TWOPOTBITS_IN_LONG);
3009: nf = nfnewprec(nf, prec);
3010: r = qf_base_change(T,u,1);
3011: i = PREC + (gexpo(r) >> TWOPOTBITS_IN_LONG);
3012: if (i < prec) prec = i;
3013: r = gmul(r, realun(prec));
3014: }
3015: else
3016: {
3017: u = lllgramintern(a,4,flag&1, (prec<<1)-2);
3018: if (!u) goto PRECPB;
3019: r = qf_base_change(a,u,1);
3020: }
3021: r = sqred1intern(r,flag&1);
3022: if (!r) goto PRECPB;
3023:
3024: n = lg(a);
3025: if (n == 1)
3026: {
3027: if (CHECK) err(talker, "dimension 0 in fincke_pohst");
3028: avma = av; z=cgetg(4,t_VEC);
3029: z[1]=z[2]=zero;
3030: z[3]=lgetg(1,t_MAT); return z;
3031: }
3032: for (i=1; i<n; i++)
3033: {
3034: GEN p1 = gsqrt(gcoeff(r,i,i), prec);
3035: coeff(r,i,i)=(long)p1;
3036: for (j=i+1; j<n; j++)
3037: coeff(r,i,j) = lmul(p1, gcoeff(r,i,j));
3038: }
3039: /* now r~ * r = a in approximate LLL basis */
3040: rinvtrans = gtrans(invmat(r));
3041:
3042: v = NULL;
3043: for (i=1; i<6; i++) /* try to get close to a genuine LLL basis */
3044: {
3045: GEN p1;
3046: if (DEBUGLEVEL>2)
3047: fprintferr("final LLLs: prec = %ld, precision(rinvtrans) = %ld\n",
3048: prec,gprecision(rinvtrans));
3049: p1 = lllintern(rinvtrans,flag&1, (gprecision(rinvtrans)<<1)-2);
3050: if (!p1) goto PRECPB;
3051: if (ishnfall(p1)) break; /* upper triangular */
3052: if (v) v = gmul(v,p1); else v = p1;
3053: rinvtrans = gmul(rinvtrans,p1);
3054: }
3055: if (i == 6) goto PRECPB; /* diverges... */
3056:
3057: if (v)
3058: {
3059: GEN u2 = ZM_inv(gtrans(v),gun);
3060: r = gmul(r,u2); /* r = LLL basis now */
3061: u = gmul(u,u2);
3062: }
3063:
3064: vnorm = cgetg(n,t_VEC);
3065: for (j=1; j<n; j++) vnorm[j] = lnorml2((GEN)rinvtrans[j]);
3066: sperm = cgetg(n,t_MAT);
3067: uperm = cgetg(n,t_MAT); perm = sindexsort(vnorm);
3068: for (i=1; i<n; i++) { uperm[n-i] = u[perm[i]]; sperm[n-i] = r[perm[i]]; }
3069:
3070: gram = gram_matrix(sperm);
3071: B = gcoeff(gram,n-1,n-1);
3072: if (gexpo(B) >= bit_accuracy(lg(B)-2)) goto PRECPB;
3073:
3074: { /* catch precision problems (truncation error) */
3075: jmp_buf env;
3076: if (setjmp(env)) goto PRECPB;
3077: catcherr = err_catch(precer, env, NULL);
3078: }
3079: if (CHECK && CHECK->f_init)
3080: {
3081: bound = CHECK->f_init(CHECK,nf,gram,uperm, &prec);
3082: if (!bound) goto PRECPB;
3083: }
3084: if (!bound)
3085: { /* set bound */
3086: GEN x = cgetg(n,t_COL);
3087: if (nf) bound = get_Bnf(nf);
3088: for (i=2; i<n; i++) x[i]=zero;
3089: for (i=1; i<n; i++)
3090: {
3091: x[i] = un; B = gcoeff(gram,i,i);
3092: if (!bound || mpcmp(B,bound) < 0) bound = B;
3093: x[i] = zero;
3094: }
3095: }
3096:
3097: if (DEBUGLEVEL>2) {fprintferr("entering smallvectors\n"); flusherr();}
3098: for (i=1; i<n; i++)
3099: {
3100: res = smallvectors(gram, bound? bound: gcoeff(gram,i,i),
3101: stockmax,flag,CHECK);
3102: if (!res) goto PRECPB;
3103: if (!CHECK || bound || lg(res[2]) > 1) break;
3104: if (DEBUGLEVEL>2) fprintferr(" i = %ld failed\n",i);
3105: }
3106: err_leave(&catcherr); catcherr = NULL;
3107: if (CHECK)
3108: {
3109: if (CHECK->f_post) res = CHECK->f_post(CHECK->data, res);
3110: return res;
3111: }
3112:
3113: if (DEBUGLEVEL>2) {fprintferr("leaving fincke_pohst\n"); flusherr();}
3114: z=cgetg(4,t_VEC);
3115: z[1]=lcopy((GEN)res[1]);
3116: z[2]=pr? lcopy((GEN)res[2]) : lround((GEN)res[2]);
3117: z[3]=lmul(uperm, (GEN)res[3]); return gerepileupto(av,z);
3118: PRECPB:
3119: if (catcherr) err_leave(&catcherr);
3120: if (!(flag & 1))
3121: err(talker,"not a positive definite form in fincke_pohst");
3122: avma = av; return NULL;
3123: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to bibli1.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 / basemath