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