Annotation of OpenXM_contrib/pari-2.2/src/basemath/buch2.c, Revision 1.1.1.1
1.1 noro 1: /* $Id: buch2.c,v 1.90 2001/10/01 12:11:30 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: /* CLASS GROUP AND REGULATOR (McCURLEY, BUCHMANN) */
19: /* GENERAL NUMBER FIELDS */
20: /* */
21: /*******************************************************************/
22: #include "pari.h"
23: #include "parinf.h"
24: extern GEN to_famat_all(GEN x, GEN y);
25: extern int addcolumntomatrix(GEN V, GEN invp, GEN L);
26: extern double check_bach(double cbach, double B);
27: extern GEN gmul_mat_smallvec(GEN x, GEN y);
28: extern GEN gmul_mati_smallvec(GEN x, GEN y);
29: extern GEN get_arch_real(GEN nf,GEN x,GEN *emb,long prec);
30: extern GEN get_roots(GEN x,long r1,long ru,long prec);
31: extern void get_nf_matrices(GEN nf, long small);
32: extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *t, long v);
33: extern GEN init_idele(GEN nf);
34: extern GEN norm_by_embed(long r1, GEN x);
35: extern void minim_alloc(long n,double ***q,long **x,double **y,double **z,double **v);
36: extern GEN idealmulpowprime(GEN nf, GEN x, GEN vp, GEN n);
37: extern GEN arch_mul(GEN x, GEN y);
38: extern GEN vecdiv(GEN x, GEN y);
39: extern GEN vecmul(GEN x, GEN y);
40: extern GEN mul_real(GEN x, GEN y);
41:
42: #define SFB_MAX 2
43: #define SFB_STEP 2
44: #define MIN_EXTRA 1
45:
46: #define RANDOM_BITS 4
47: static const int CBUCHG = (1<<RANDOM_BITS) - 1;
48: static const int randshift = BITS_IN_RANDOM-1 - RANDOM_BITS;
49: #undef RANDOM_BITS
50:
51: static long KC,KCZ,KCZ2,MAXRELSUP;
52: static long primfact[500],expoprimfact[500];
53: static GEN *idealbase, vectbase, FB, numFB, powsubFB, numideal;
54:
55: /* FB[i] i-th rational prime used in factor base
56: * numFB[i] index k such that FB[k]=i (0 if i is not prime)
57: *
58: * vectbase vector of all ideals in FB
59: * vecbase o subFB part of FB used to build random relations
60: * powsubFB array lg(subFB) x (CBUCHG+1) all powers up to CBUCHG
61: *
62: * idealbase[i] prime ideals above i in FB
63: * numideal[i] index k such that idealbase[k] = i.
64: *
65: * mat all relations found (as long integers, not reduced)
66: * cglob lg(mat) = total number of relations found so far
67: *
68: * Use only non-inert primes, coprime to discriminant index F:
69: * KC = number of prime ideals in factor base (norm < Bach cst)
70: * KC2= number of prime ideals assumed to generate class group (>= KC)
71: *
72: * KCZ = number of rational primes under ideals counted by KC
73: * KCZ2= same for KC2
74: */
75:
76: /* x a t_VECSMALL */
77: static long
78: ccontent(GEN x)
79: {
80: long i, l = lg(x), s=labs(x[1]);
81: for (i=2; i<l && s>1; i++) s = cgcd(x[i],s);
82: return s;
83: }
84:
85: static void
86: desallocate(GEN **M)
87: {
88: GEN *m = *M;
89: long i;
90: if (m)
91: {
92: for (i=lg(m)-1; i; i--) free(m[i]);
93: free(m); *M = NULL;
94: }
95: }
96:
97: /* Return the list of indexes or the primes chosen for the subFB.
98: * Fill vperm (if !=0): primes ideals sorted by increasing norm (except the
99: * ones in subFB come first [dense rows come first for hnfspec])
100: * ss = number of rational primes whose divisors are all in FB
101: */
102: static GEN
103: subFBgen(long N,long m,long minsFB,GEN vperm, long *ptss)
104: {
105: long av = avma,i,j, lv=lg(vectbase),s=0,s1=0,n=0,ss=0,z=0;
106: GEN y1,y2,subFB,perm,perm1,P,Q;
107: double prod;
108:
109: (void)new_chunk(lv); /* room for subFB */
110: y1 = cgetg(lv,t_COL);
111: y2 = cgetg(lv,t_COL);
112: for (i=1,P=(GEN)vectbase[i];;P=Q)
113: { /* we'll sort ideals by norm (excluded ideals = "zero") */
114: long e = itos((GEN)P[3]);
115: long ef= e*itos((GEN)P[4]);
116:
117: s1 += ef;
118: y2[i] = (long)powgi((GEN)P[1],(GEN)P[4]);
119: /* take only unramified ideals */
120: if (e>1) { y1[i]=zero; s=0; z++; } else { y1[i]=y2[i]; s += ef; }
121:
122: i++; Q = (GEN)vectbase[i];
123: if (i == lv || !egalii((GEN)P[1], (GEN)Q[1]))
124: { /* don't take all P above a given p (delete the last one) */
125: if (s == N) { y1[i-1]=zero; z++; }
126: if (s1== N) ss++;
127: if (i == lv) break;
128: s = s1 = 0;
129: }
130: }
131: if (z+minsFB >= lv) return NULL;
132:
133: prod = 1.0;
134: perm = sindexsort(y1) + z; /* skip "zeroes" (excluded ideals) */
135: for(;;)
136: {
137: if (++n > minsFB && (z+n >= lv || prod > m + 0.5)) break;
138: prod *= gtodouble((GEN)y1[perm[n]]);
139: }
140: if (prod < m) return NULL;
141: n--;
142:
143: /* take the first (non excluded) n ideals (wrt norm), put them first, and
144: * sort the rest by increasing norm */
145: for (j=1; j<=n; j++) y2[perm[j]] = zero;
146: perm1 = sindexsort(y2); avma = av;
147:
148: subFB = cgetg(n+1, t_VECSMALL);
149: if (vperm)
150: {
151: for (j=1; j<=n; j++) vperm[j] = perm[j];
152: for ( ; j<lv; j++) vperm[j] = perm1[j];
153: }
154: for (j=1; j<=n; j++) subFB[j] = perm[j];
155:
156: if (DEBUGLEVEL)
157: {
158: if (DEBUGLEVEL>3)
159: {
160: fprintferr("\n***** IDEALS IN FACTORBASE *****\n\n");
161: for (i=1; i<=KC; i++) fprintferr("no %ld = %Z\n",i,vectbase[i]);
162: fprintferr("\n***** IDEALS IN SUB FACTORBASE *****\n\n");
163: outerr(vecextract_p(vectbase,subFB));
164: fprintferr("\n***** INITIAL PERMUTATION *****\n\n");
165: fprintferr("vperm = %Z\n\n",vperm);
166: }
167: msgtimer("sub factorbase (%ld elements)",n);
168: }
169: powsubFB = NULL;
170: *ptss = ss; return subFB;
171: }
172:
173: static GEN
174: mulred(GEN nf,GEN x, GEN I, long prec)
175: {
176: ulong av = avma;
177: GEN y = cgetg(3,t_VEC);
178:
179: y[1] = (long)idealmulh(nf,I,(GEN)x[1]);
180: y[2] = x[2];
181: y = ideallllred(nf,y,NULL,prec);
182: y[1] = (long)ideal_two_elt(nf,(GEN)y[1]);
183: return gerepilecopy(av, y);
184: }
185:
186: /* Compute powers of prime ideals (P^0,...,P^a) in subFB (a > 1)
187: * powsubFB[j][i] contains P_i^j in LLL form + archimedean part in
188: * MULTIPLICATIVE form (logs are expensive)
189: */
190: static void
191: powsubFBgen(GEN nf,GEN subFB,long a,long prec)
192: {
193: long i,j, n = lg(subFB), RU = lg(nf[6]);
194: GEN *pow, arch0 = cgetg(RU,t_COL);
195: for (i=1; i<RU; i++) arch0[i] = un; /* 0 in multiplicative notations */
196:
197: if (DEBUGLEVEL) fprintferr("Computing powers for sub-factor base:\n");
198: powsubFB = cgetg(n, t_VEC);
199: for (i=1; i<n; i++)
200: {
201: GEN vp = (GEN)vectbase[subFB[i]];
202: GEN z = cgetg(3,t_VEC); z[1]=vp[1]; z[2]=vp[2];
203: pow = (GEN*)cgetg(a+1,t_VEC); powsubFB[i] = (long)pow;
204: pow[1]=cgetg(3,t_VEC);
205: pow[1][1] = (long)z;
206: pow[1][2] = (long)arch0;
207: vp = prime_to_ideal(nf,vp);
208: for (j=2; j<=a; j++)
209: {
210: pow[j] = mulred(nf,pow[j-1],vp,prec);
211: if (DEBUGLEVEL>1) fprintferr(" %ld",j);
212: }
213: if (DEBUGLEVEL>1) { fprintferr("\n"); flusherr(); }
214: }
215: if (DEBUGLEVEL) msgtimer("powsubFBgen");
216: }
217:
218: /* Compute FB, numFB, idealbase, vectbase, numideal.
219: * n2: bound for norm of tested prime ideals (includes be_honest())
220: * n : bound for prime ideals used to build relations (Bach cst) ( <= n2 )
221:
222: * Return prod_{p<=n2} (1-1/p) / prod_{Norm(P)<=n2} (1-1/Norm(P)),
223: * close to residue of zeta_K at 1 = 2^r1 (2pi)^r2 h R / (w D)
224: */
225: static GEN
226: FBgen(GEN nf,long n2,long n)
227: {
228: byteptr delta=diffptr;
229: long KC2,i,j,k,p,lon,ip,nor, N = degpol(nf[1]);
230: GEN p2,p1,NormP,lfun;
231: long prim[] = { evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),0 };
232:
233: numFB = cgetg(n2+1,t_VECSMALL);
234: FB = cgetg(n2+1,t_VECSMALL);
235: numideal = cgetg(n2+1,t_VECSMALL);
236: idealbase= (GEN*)cgetg(n2+1,t_VEC);
237:
238: lfun=realun(DEFAULTPREC);
239: p=*delta++; i=0; ip=0; KC=0;
240: while (p<=n2)
241: {
242: long av = avma, av1;
243: if (DEBUGLEVEL>=2) { fprintferr(" %ld",p); flusherr(); }
244: prim[2] = p; p1 = primedec(nf,prim); lon=lg(p1);
245: av1 = avma;
246: divrsz(mulsr(p-1,lfun),p,lfun);
247: if (itos(gmael(p1,1,4)) == N) /* p inert */
248: {
249: NormP = gpowgs(prim,N);
250: if (!is_bigint(NormP) && (nor=NormP[2]) <= n2)
251: divrsz(mulsr(nor,lfun),nor-1, lfun);
252: avma = av1;
253: }
254: else
255: {
256: numideal[p]=ip;
257: i++; numFB[p]=i; FB[i]=p;
258: for (k=1; k<lon; k++,ip++)
259: {
260: NormP = powgi(prim,gmael(p1,k,4));
261: if (is_bigint(NormP) || (nor=NormP[2]) > n2) break;
262:
263: divrsz(mulsr(nor,lfun),nor-1, lfun);
264: }
265: /* keep all ideals with Norm <= n2 */
266: avma = av1;
267: if (k == lon)
268: setisclone(p1); /* flag it: all prime divisors in FB */
269: else
270: { setlg(p1,k); p1 = gerepilecopy(av,p1); }
271: idealbase[i] = p1;
272: }
273: if (!*delta) err(primer1);
274: p += *delta++;
275: if (KC == 0 && p>n) { KCZ=i; KC=ip; }
276: }
277: if (!KC) return NULL;
278: KCZ2=i; KC2=ip; MAXRELSUP = min(50,4*KC);
279:
280: setlg(FB,KCZ2);
281: setlg(numFB,KCZ2);
282: setlg(numideal,KCZ2);
283: setlg(idealbase,KCZ2);
284: vectbase=cgetg(KC+1,t_COL);
285: for (i=1; i<=KCZ; i++)
286: {
287: p1 = idealbase[i]; k=lg(p1);
288: p2 = vectbase + numideal[FB[i]];
289: for (j=1; j<k; j++) p2[j]=p1[j];
290: }
291: if (DEBUGLEVEL)
292: {
293: if (DEBUGLEVEL>1) fprintferr("\n");
294: if (DEBUGLEVEL>6)
295: {
296: fprintferr("########## FACTORBASE ##########\n\n");
297: fprintferr("KC2=%ld, KC=%ld, KCZ=%ld, KCZ2=%ld, MAXRELSUP=%ld\n",
298: KC2, KC, KCZ, KCZ2, MAXRELSUP);
299: for (i=1; i<=KCZ; i++)
300: fprintferr("++ idealbase[%ld] = %Z",i,idealbase[i]);
301: }
302: msgtimer("factor base");
303: }
304: return lfun;
305: }
306:
307: /* can we factor I / m ? (m pseudo minimum, computed in ideallllredpart1) */
308: static long
309: factorgen(GEN nf,GEN idealvec,long kcz,long limp)
310: {
311: long i,j,n1,ip,v,p,k,lo,ifinal;
312: GEN x,q,r,P,p1,listexpo;
313: GEN I = (GEN)idealvec[1];
314: GEN m = (GEN)idealvec[2];
315: GEN Nm= absi( subres(gmul((GEN)nf[7],m), (GEN)nf[1]) ); /* |Nm| */
316:
317: x = divii(Nm, dethnf_i(I)); /* m in I, so NI | Nm */
318: if (is_pm1(x)) { primfact[0]=0; return 1; }
319: listexpo = new_chunk(kcz+1);
320: for (i=1; ; i++)
321: {
322: p=FB[i]; q=dvmdis(x,p,&r);
323: for (k=0; !signe(r); k++) { x=q; q=dvmdis(x,p,&r); }
324: listexpo[i] = k;
325: if (cmpis(q,p)<=0) break;
326: if (i==kcz) return 0;
327: }
328: if (cmpis(x,limp) > 0) return 0;
329:
330: ifinal = i; lo = 0;
331: for (i=1; i<=ifinal; i++)
332: {
333: k = listexpo[i];
334: if (k)
335: {
336: p = FB[i]; p1 = idealbase[numFB[p]];
337: n1 = lg(p1); ip = numideal[p];
338: for (j=1; j<n1; j++)
339: {
340: P = (GEN)p1[j];
341: v = idealval(nf,I, P) - element_val2(nf,m,Nm, P);
342: if (v) /* hence < 0 */
343: {
344: primfact[++lo]=ip+j; expoprimfact[lo]=v;
345: k += v * itos((GEN)P[4]);
346: if (!k) break;
347: }
348: }
349: if (k) return 0;
350: }
351: }
352: if (is_pm1(x)) { primfact[0]=lo; return 1; }
353:
354: p = itos(x); p1 = idealbase[numFB[p]];
355: n1 = lg(p1); ip = numideal[p];
356: for (k=1,j=1; j<n1; j++)
357: {
358: P = (GEN)p1[j];
359: v = idealval(nf,I, P) - element_val2(nf,m,Nm, P);
360: if (v)
361: {
362: primfact[++lo]=ip+j; expoprimfact[lo]=v;
363: k += v*itos((GEN)P[4]);
364: if (!k) { primfact[0]=lo; return 1; }
365: }
366: }
367: return 0;
368: }
369:
370: /* can we factor x ? Nx = norm(x) */
371: static long
372: factorelt(GEN nf,GEN cbase,GEN x,GEN Nx,long kcz,long limp)
373: {
374: long i,j,n1,ip,v,p,k,lo,ifinal;
375: GEN q,r,P,p1,listexpo;
376:
377: if (is_pm1(Nx)) { primfact[0]=0; return 1; }
378: listexpo = new_chunk(kcz+1);
379: for (i=1; ; i++)
380: {
381: p=FB[i]; q=dvmdis(Nx,p,&r);
382: for (k=0; !signe(r); k++) { Nx=q; q=dvmdis(Nx,p,&r); }
383: listexpo[i] = k;
384: if (cmpis(q,p)<=0) break;
385: if (i==kcz) return 0;
386: }
387: if (cmpis(Nx,limp) > 0) return 0;
388:
389: if (cbase) x = gmul(cbase,x);
390: ifinal=i; lo = 0;
391: for (i=1; i<=ifinal; i++)
392: {
393: k = listexpo[i];
394: if (k)
395: {
396: p = FB[i]; p1 = idealbase[numFB[p]];
397: n1 = lg(p1); ip = numideal[p];
398: for (j=1; j<n1; j++)
399: {
400: P = (GEN)p1[j];
401: v = int_elt_val(nf,x,(GEN)P[1],(GEN)P[5], NULL, k);
402: if (v)
403: {
404: primfact[++lo]=ip+j; expoprimfact[lo]=v;
405: k -= v * itos((GEN)P[4]);
406: if (!k) break;
407: }
408: }
409: if (k) return 0;
410: }
411: }
412: if (is_pm1(Nx)) { primfact[0]=lo; return 1; }
413:
414: p = itos(Nx); p1 = idealbase[numFB[p]];
415: n1 = lg(p1); ip = numideal[p];
416: for (k=1,j=1; j<n1; j++)
417: {
418: P = (GEN)p1[j];
419: v = int_elt_val(nf,x,(GEN)P[1],(GEN)P[5], NULL, k);
420: if (v)
421: {
422: primfact[++lo]=ip+j; expoprimfact[lo]=v;
423: k -= v*itos((GEN)P[4]);
424: if (!k) { primfact[0]=lo; return 1; }
425: }
426: }
427: return 0;
428: }
429:
430: static GEN
431: cleancol(GEN x,long N,long PRECREG)
432: {
433: long i,j,av,tetpil,tx=typ(x),R1,RU;
434: GEN s,s2,re,pi4,im,y;
435:
436: if (tx==t_MAT)
437: {
438: y=cgetg(lg(x),tx);
439: for (j=1; j<lg(x); j++)
440: y[j]=(long)cleancol((GEN)x[j],N,PRECREG);
441: return y;
442: }
443: if (!is_vec_t(tx)) err(talker,"not a vector/matrix in cleancol");
444: av = avma; RU=lg(x)-1; R1 = (RU<<1)-N;
445: re = greal(x); s=(GEN)re[1]; for (i=2; i<=RU; i++) s=gadd(s,(GEN)re[i]);
446: im = gimag(x); s = gdivgs(s,-N);
447: s2 = (N>R1)? gmul2n(s,1): NULL;
448: pi4=gmul2n(mppi(PRECREG),2);
449: tetpil=avma; y=cgetg(RU+1,tx);
450: for (i=1; i<=RU; i++)
451: {
452: GEN p1=cgetg(3,t_COMPLEX); y[i]=(long)p1;
453: p1[1] = ladd((GEN)re[i], (i<=R1)?s:s2);
454: p1[2] = lmod((GEN)im[i], pi4);
455: }
456: return gerepile(av,tetpil,y);
457: }
458:
459: #define RELAT 0
460: #define LARGE 1
461: #define PRECI 2
462: static GEN
463: not_given(long av, long flun, long reason)
464: {
465: if (labs(flun)==2)
466: {
467: char *s;
468: switch(reason)
469: {
470: case RELAT: s = "not enough relations for fundamental units"; break;
471: case LARGE: s = "fundamental units too large"; break;
472: case PRECI: s = "insufficient precision for fundamental units"; break;
473: default: s = "unknown problem with fundamental units";
474: }
475: err(warner,"%s, not given",s);
476: }
477: avma=av; return cgetg(1,t_MAT);
478: }
479:
480: /* check whether exp(x) will get too big */
481: static long
482: expgexpo(GEN x)
483: {
484: long i,j,e, E = -HIGHEXPOBIT;
485: GEN p1;
486:
487: for (i=1; i<lg(x); i++)
488: for (j=1; j<lg(x[1]); j++)
489: {
490: p1 = gmael(x,i,j);
491: if (typ(p1)==t_COMPLEX) p1 = (GEN)p1[1];
492: e = gexpo(p1); if (e>E) E=e;
493: }
494: return E;
495: }
496:
497: static GEN
498: split_realimag_col(GEN z, long r1, long r2)
499: {
500: long i, ru = r1+r2;
501: GEN a, x = cgetg(ru+r2+1,t_COL), y = x + r2;
502: for (i=1; i<=r1; i++) { a = (GEN)z[i]; x[i] = lreal(a); }
503: for ( ; i<=ru; i++) { a = (GEN)z[i]; x[i] = lreal(a); y[i] = limag(a); }
504: return x;
505: }
506:
507: static GEN
508: split_realimag(GEN x, long r1, long r2)
509: {
510: long i,l; GEN y;
511: if (typ(x) == t_COL) return split_realimag_col(x,r1,r2);
512: l = lg(x); y = cgetg(l, t_MAT);
513: for (i=1; i<l; i++) y[i] = (long)split_realimag_col((GEN)x[i], r1, r2);
514: return y;
515: }
516:
517: /* assume x = (r1+r2) x (r1+2r2) matrix and y compatible vector
518: * r1 first lines of x,y are real. Solve the system obtained by splitting
519: * real and imaginary parts. If x is of nf type, use M instead.
520: */
521: static GEN
522: gauss_realimag(GEN x, GEN y)
523: {
524: GEN M = (typ(x)==t_VEC)? gmael(checknf(x),5,1): x;
525: long l = lg(M), r2 = l - lg(M[1]), r1 = l-1 - 2*r2;
526: M = split_realimag(M,r1,r2);
527: y = split_realimag(y,r1,r2); return gauss(M, y);
528: }
529:
530: static GEN
531: getfu(GEN nf,GEN *ptxarch,GEN reg,long flun,long *pte,long prec)
532: {
533: long av = avma,e,i,j,R1,RU,N=degpol(nf[1]);
534: GEN p1,p2,u,y,matep,s,xarch,vec;
535: GEN *gptr[2];
536:
537: if (DEBUGLEVEL) fprintferr("\n#### Computing fundamental units\n");
538: R1 = itos(gmael(nf,2,1)); RU = (N+R1)>>1;
539: if (RU==1) { *pte=BIGINT; return cgetg(1,t_VEC); }
540:
541: *pte = 0; xarch = *ptxarch;
542: if (gexpo(reg) < -8) return not_given(av,flun,RELAT);
543:
544: matep = cgetg(RU,t_MAT);
545: for (j=1; j<RU; j++)
546: {
547: s = gzero; for (i=1; i<=RU; i++) s = gadd(s,greal(gcoeff(xarch,i,j)));
548: s = gdivgs(s, -N);
549: p1=cgetg(RU+1,t_COL); matep[j]=(long)p1;
550: for (i=1; i<=R1; i++) p1[i] = ladd(s, gcoeff(xarch,i,j));
551: for ( ; i<=RU; i++) p1[i] = ladd(s, gmul2n(gcoeff(xarch,i,j),-1));
552: }
553: if (prec <= 0) prec = gprecision(xarch);
554: u = lllintern(greal(matep),1,prec);
555: if (!u) return not_given(av,flun,PRECI);
556:
557: p1 = gmul(matep,u);
558: if (expgexpo(p1) > 20) { *pte = BIGINT; return not_given(av,flun,LARGE); }
559: matep = gexp(p1,prec);
560: y = grndtoi(gauss_realimag(nf,matep), &e);
561: *pte = -e;
562: if (e>=0) return not_given(av,flun,PRECI);
563: for (j=1; j<RU; j++)
564: if (!gcmp1(idealnorm(nf, (GEN)y[j]))) break;
565: if (j < RU) { *pte = 0; return not_given(av,flun,PRECI); }
566: xarch = gmul(xarch,u);
567:
568: /* y[i] are unit generators. Normalize: smallest L2 norm + lead coeff > 0 */
569: y = gmul((GEN)nf[7], y);
570: vec = cgetg(RU+1,t_COL); p2 = mppi(prec);
571: p1 = pureimag(p2);
572: p2 = pureimag(gmul2n(p2,1));
573: for (i=1; i<=R1; i++) vec[i]=(long)p1;
574: for ( ; i<=RU; i++) vec[i]=(long)p2;
575: for (j=1; j<RU; j++)
576: {
577: p1 = (GEN)y[j]; p2 = ginvmod(p1, (GEN)nf[1]);
578: if (gcmp(QuickNormL2(p2,DEFAULTPREC),
579: QuickNormL2(p1,DEFAULTPREC)) < 0)
580: {
581: xarch[j] = lneg((GEN)xarch[j]);
582: p1 = p2;
583: }
584: if (gsigne(leading_term(p1)) < 0)
585: {
586: xarch[j] = ladd((GEN)xarch[j], vec);
587: p1 = gneg(p1);
588: }
589: y[j] = (long)p1;
590: }
591: if (DEBUGLEVEL) msgtimer("getfu");
592: *ptxarch=xarch; gptr[0]=ptxarch; gptr[1]=&y;
593: gerepilemany(av,gptr,2); return y;
594: }
595: #undef RELAT
596: #undef LARGE
597: #undef PRECI
598:
599: GEN
600: buchfu(GEN bnf)
601: {
602: long av = avma, c;
603: GEN nf,xarch,reg,res, y = cgetg(3,t_VEC);
604:
605: bnf = checkbnf(bnf); xarch = (GEN)bnf[3]; nf = (GEN)bnf[7];
606: res = (GEN)bnf[8]; reg = (GEN)res[2];
607: if (lg(res)==7 && lg(res[5])==lg(nf[6])-1)
608: {
609: y[1] = lcopy((GEN)res[5]);
610: y[2] = lcopy((GEN)res[6]); return y;
611: }
612: y[1] = (long)getfu(nf,&xarch,reg,2,&c,0);
613: y[2] = lstoi(c); return gerepilecopy(av, y);
614: }
615:
616: /*******************************************************************/
617: /* */
618: /* PRINCIPAL IDEAL ALGORITHM (DISCRETE LOG) */
619: /* */
620: /*******************************************************************/
621:
622: /* gen: prime ideals, return Norm (product), bound for denominator.
623: * C = possible extra prime (^1) or NULL */
624: static GEN
625: get_norm_fact_primes(GEN gen, GEN ex, GEN C, GEN *pd)
626: {
627: GEN N,d,P,p,e;
628: long i,s,c = lg(ex);
629: d = N = gun;
630: for (i=1; i<c; i++)
631: if ((s = signe(ex[i])))
632: {
633: P = (GEN)gen[i]; e = (GEN)ex[i]; p = (GEN)P[1];
634: N = gmul(N, powgi(p, mulii((GEN)P[4], e)));
635: if (s < 0)
636: {
637: e = gceil(gdiv(negi(e), (GEN)P[3]));
638: d = mulii(d, powgi(p, e));
639: }
640: }
641: if (C)
642: {
643: P = C; p = (GEN)P[1];
644: N = gmul(N, powgi(p, (GEN)P[4]));
645: }
646: *pd = d; return N;
647: }
648:
649: /* gen: HNF ideals */
650: static GEN
651: get_norm_fact(GEN gen, GEN ex, GEN *pd)
652: {
653: long i, c = lg(ex);
654: GEN d,N,I,e,n,ne,de;
655: d = N = gun;
656: for (i=1; i<c; i++)
657: if (signe(ex[i]))
658: {
659: I = (GEN)gen[i]; e = (GEN)ex[i]; n = dethnf_i(I);
660: ne = powgi(n,e);
661: de = egalii(n, gcoeff(I,1,1))? ne: powgi(gcoeff(I,1,1), e);
662: N = mulii(N, ne);
663: d = mulii(d, de);
664: }
665: *pd = d; return N;
666: }
667:
668: /* try to split ideal / (some integer) on FB */
669: static long
670: factorgensimple(GEN nf,GEN ideal)
671: {
672: long N,i,v,av1 = avma,lo, L = lg(vectbase);
673: GEN x;
674: if (typ(ideal) != t_MAT) ideal = (GEN)ideal[1]; /* idele */
675: x = dethnf_i(ideal);
676: N = lg(ideal)-1;
677: if (gcmp1(x)) { avma=av1; primfact[0]=0; return 1; }
678: for (lo=0, i=1; i<L; i++)
679: {
680: GEN q,y, P=(GEN)vectbase[i], p=(GEN)P[1];
681: long vx0, vx, sum_ef, e,f,lo0,i0;
682:
683: vx = pvaluation(x,p,&y);
684: if (!vx) continue;
685: /* p | x = Nideal */
686: vx0 = vx; sum_ef = 0; lo0 =lo; i0 = i;
687: for(;;)
688: {
689: e = itos((GEN)P[3]);
690: f = itos((GEN)P[4]); sum_ef += e * f;
691: v = idealval(nf,ideal,P);
692: if (v)
693: {
694: lo++; primfact[lo]=i; expoprimfact[lo]=v;
695: vx -= v * f; if (!vx) break;
696: }
697: i++; if (i == L) break;
698: P = (GEN)vectbase[i]; q = (GEN)P[1];
699: if (!egalii(p,q)) break;
700: }
701: if (vx)
702: { /* divisible by P | p not in FB */
703: long k,l,n;
704: k = N - sum_ef;
705: if (vx % k) break;
706: k = vx / k; /* ideal / p^k may factor on FB */
707: for (l = lo0+1; l <= lo; l++)
708: {
709: P = (GEN)vectbase[primfact[l]];
710: expoprimfact[l] -= k * itos((GEN)P[3]); /* may become 0 */
711: }
712: n = lo0+1;
713: for (l = i0; l < i; l++) /* update exponents for other P | p */
714: {
715: if (n <= lo && primfact[n] == l) { n++; continue; }
716: lo++; primfact[lo] = l; P = (GEN)vectbase[l];
717: expoprimfact[lo] = - k * itos((GEN)P[3]);
718: }
719: /* check that ideal / p^k / (tentative factorisation) is integral */
720: {
721: GEN z = ideal;
722: for (l = lo0+1; l <= lo; l++)
723: {
724: GEN n = stoi(- expoprimfact[l]);
725: z = idealmulpowprime(nf, z, (GEN)vectbase[primfact[l]], n);
726: }
727: z = gdiv(z, gpowgs(p, k));
728: if (!gcmp1(denom(z))) { avma=av1; return 0; }
729: ideal = z;
730: }
731: }
732: if (gcmp1(y)) { avma=av1; primfact[0]=lo; return 1; }
733: x = y;
734: }
735: avma=av1; return 0;
736: }
737:
738: static void
739: add_to_fact(long l, long v, long e)
740: {
741: long i,lo;
742: if (!e) return;
743: for (i=1; i<=l && primfact[i] < v; i++)/*empty*/;
744: if (i <= l && primfact[i] == v)
745: expoprimfact[i] += e;
746: else
747: {
748: lo = ++primfact[0];
749: primfact[lo] = v;
750: expoprimfact[lo] = e;
751: }
752: }
753:
754: static void
755: init_sub(long l, GEN perm, GEN *v, GEN *ex)
756: {
757: long i;
758: *v = cgetg(l,t_VECSMALL);
759: *ex= cgetg(l,t_VECSMALL);
760: for (i=1; i<l; i++) (*v)[i] = itos((GEN)perm[i]);
761: }
762:
763: /* factor x = x0[1] on vectbase (modulo principal ideals).
764: * We may have x0[2] = NULL (principal part not wanted), in which case,
765: * don't compute archimedean parts */
766: static GEN
767: split_ideal(GEN nf, GEN x0, long prec, GEN vperm)
768: {
769: GEN p1,vdir,id,z,v,ex,y, x = (GEN)x0[1];
770: long nbtest_lim,nbtest,bou,i,ru, lgsub;
771: int flag = (gexpo(gcoeff(x,1,1)) < 100);
772:
773: y = x0;
774: if (flag && factorgensimple(nf,y)) return y;
775:
776: y = ideallllred(nf,x0,NULL,prec);
777: if (flag && ((!x0[2] && gegal((GEN)y[1], (GEN)x[1]))
778: || (x0[2] && gcmp0((GEN)y[2])))) flag = 0; /* y == x0 */
779: if (flag && factorgensimple(nf,y)) return y;
780:
781: z = init_idele(nf); ru = lg(z[2]);
782: if (!x0[2]) { z[2] = 0; x0 = x; } /* stop cheating */
783: vdir = cgetg(ru,t_VEC);
784: for (i=2; i<ru; i++) vdir[i]=zero;
785: for (i=1; i<ru; i++)
786: {
787: vdir[i]=lstoi(10);
788: y = ideallllred(nf,x0,vdir,prec);
789: if (factorgensimple(nf,y)) return y;
790: vdir[i]=zero;
791: }
792: nbtest = 0; nbtest_lim = (ru-1) << 2; lgsub = 3;
793: init_sub(lgsub, vperm, &v, &ex);
794: for(;;)
795: {
796: int non0 = 0;
797: id = x0;
798: for (i=1; i<lgsub; i++)
799: {
800: ex[i] = mymyrand() >> randshift;
801: if (ex[i])
802: { /* don't let id become too large as lgsub gets bigger: avoid
803: prec problems */
804: if (non0) id = ideallllred(nf,id,NULL,prec);
805: non0++;
806: z[1]=vectbase[v[i]]; p1=idealpowred(nf,z,stoi(ex[i]),prec);
807: id = idealmulh(nf,id,p1);
808: }
809: }
810: if (id == x0) continue;
811:
812: for (i=1; i<ru; i++) vdir[i] = lstoi(mymyrand() >> randshift);
813: for (bou=1; bou<ru; bou++)
814: {
815: if (bou>1)
816: {
817: for (i=1; i<ru; i++) vdir[i]=zero;
818: vdir[bou]=lstoi(10);
819: }
820: nbtest++;
821: y = ideallllred(nf,id,vdir,prec);
822: if (DEBUGLEVEL>2)
823: fprintferr("nbtest = %ld, ideal = %Z\n",nbtest,y[1]);
824: if (factorgensimple(nf,y))
825: {
826: long l = primfact[0];
827: for (i=1; i<lgsub; i++) add_to_fact(l,v[i],-ex[i]);
828: return y;
829: }
830: }
831: if (nbtest > nbtest_lim)
832: {
833: nbtest = 0;
834: if (lgsub < 7)
835: {
836: lgsub++; nbtest_lim <<= 2;
837: init_sub(lgsub, vperm, &v, &ex);
838: }
839: else nbtest_lim = VERYBIGINT; /* don't increase further */
840: if (DEBUGLEVEL)
841: fprintferr("split_ideal: increasing factor base [%ld]\n",lgsub);
842: }
843: }
844: }
845:
846: static void
847: get_split_expo(GEN xalpha, GEN yalpha, GEN vperm)
848: {
849: long i, colW = lg(xalpha)-1;
850: GEN vinvperm = new_chunk(lg(vectbase));
851: for (i=1; i<lg(vectbase); i++) vinvperm[itos((GEN)vperm[i])]=i;
852: for (i=1; i<=primfact[0]; i++)
853: {
854: long k = vinvperm[primfact[i]];
855: long l = k - colW;
856: if (l <= 0) xalpha[k]=lstoi(expoprimfact[i]);
857: else yalpha[l]=lstoi(expoprimfact[i]);
858: }
859: }
860:
861: GEN
862: init_red_mod_units(GEN bnf, long prec)
863: {
864: GEN z, s = gzero, p1,s1,mat, matunit = (GEN)bnf[3];
865: long i,j, RU = lg(matunit);
866:
867: if (RU == 1) return NULL;
868: mat = cgetg(RU,t_MAT);
869: for (j=1; j<RU; j++)
870: {
871: p1=cgetg(RU+1,t_COL); mat[j]=(long)p1;
872: s1=gzero;
873: for (i=1; i<RU; i++)
874: {
875: p1[i] = lreal(gcoeff(matunit,i,j));
876: s1 = gadd(s1, gsqr((GEN)p1[i]));
877: }
878: p1[RU]=zero; if (gcmp(s1,s) > 0) s = s1;
879: }
880: s = gsqrt(gmul2n(s,RU),prec);
881: if (gcmpgs(s,100000000) < 0) s = stoi(100000000);
882: z = cgetg(3,t_VEC);
883: z[1] = (long)mat;
884: z[2] = (long)s; return z;
885: }
886:
887: /* z computed above. Return unit exponents that would reduce col (arch) */
888: GEN
889: red_mod_units(GEN col, GEN z, long prec)
890: {
891: long i,RU;
892: GEN x,mat,N2;
893:
894: if (!z) return NULL;
895: mat= (GEN)z[1];
896: N2 = (GEN)z[2];
897: RU = lg(mat); x = cgetg(RU+1,t_COL);
898: for (i=1; i<RU; i++) x[i]=lreal((GEN)col[i]);
899: x[RU] = (long)N2;
900: x = lllintern(concatsp(mat,x), 1,prec);
901: if (!x) return NULL;
902: x = (GEN)x[RU];
903: if (signe(x[RU]) < 0) x = gneg_i(x);
904: if (!gcmp1((GEN)x[RU])) err(bugparier,"red_mod_units");
905: setlg(x,RU); return x;
906: }
907:
908: /* clg2 format changed for version 2.0.21 (contained ideals, now archs)
909: * Compatibility mode: old clg2 format */
910: static GEN
911: get_Garch(GEN nf, GEN gen, GEN clg2, long prec)
912: {
913: long i,c;
914: GEN g,z,J,Garch, clorig = (GEN)clg2[3];
915:
916: c = lg(gen); Garch = cgetg(c,t_MAT);
917: for (i=1; i<c; i++)
918: {
919: g = (GEN)gen[i];
920: z = (GEN)clorig[i]; J = (GEN)z[1];
921: if (!gegal(g,J))
922: {
923: z = idealinv(nf,z); J = (GEN)z[1];
924: J = gmul(J,denom(J));
925: if (!gegal(g,J))
926: {
927: z = ideallllred(nf,z,NULL,prec); J = (GEN)z[1];
928: if (!gegal(g,J))
929: err(bugparier,"isprincipal (incompatible bnf generators)");
930: }
931: }
932: Garch[i] = z[2];
933: }
934: return Garch;
935: }
936:
937: /* [x] archimedian components, A column vector. return [x] A
938: * x may be a translated GEN (y + k) */
939: static GEN
940: act_arch(GEN A, GEN x)
941: {
942: GEN z, a;
943: long i,l = lg(A);
944: if (typ(A) == t_MAT)
945: {
946: /* assume lg(x) >= l */
947: z = cgetg(l, t_VEC);
948: for (i=1; i<l; i++) z[i] = (long)act_arch((GEN)A[i], x);
949: return z;
950: }
951: /* A a t_COL of t_INT. Assume lg(A)==lg(x) */
952: l = lg(A); z = cgetg(l, t_VEC);
953: if (l==1) return z;
954: a = gmul((GEN)A[1], (GEN)x[1]);
955: for (i=2; i<l; i++)
956: if (signe(A[i])) a = gadd(a, gmul((GEN)A[i], (GEN)x[i]));
957: settyp(a, t_VEC); return a;
958: }
959:
960: static long
961: prec_arch(GEN bnf)
962: {
963: GEN a = (GEN)bnf[4];
964: long i, l = lg(a), prec;
965:
966: for (i=1; i<l; i++)
967: if ( (prec = gprecision((GEN)a[i])) ) return prec;
968: return DEFAULTPREC;
969: }
970:
971: /* col = archimedian components of x, Nx = kNx^e its norm, dx a bound for its
972: * denominator. Return x or NULL (fail) */
973: GEN
974: isprincipalarch(GEN bnf, GEN col, GEN kNx, GEN e, GEN dx, long *pe)
975: {
976: GEN nf, x, matunit, s;
977: long N, R1, RU, i, prec = gprecision(col);
978: bnf = checkbnf(bnf); nf = checknf(bnf);
979: if (!prec) prec = prec_arch(bnf);
980: matunit = (GEN)bnf[3];
981: N = degpol(nf[1]);
982: R1 = itos(gmael(nf,2,1));
983: RU = (N + R1)>>1;
984: col = cleancol(col,N,prec); settyp(col, t_COL);
985: if (RU > 1)
986: { /* reduce mod units */
987: GEN u, z = init_red_mod_units(bnf,prec);
988: u = red_mod_units(col,z,prec);
989: if (!u && z) return NULL;
990: if (u) col = gadd(col, gmul(matunit, u));
991: }
992: s = gdivgs(gmul(e, glog(kNx,prec)), N);
993: for (i=1; i<=R1; i++) col[i] = lexp(gadd(s, (GEN)col[i]),prec);
994: for ( ; i<=RU; i++) col[i] = lexp(gadd(s, gmul2n((GEN)col[i],-1)),prec);
995: /* d.alpha such that x = alpha \prod gj^ej */
996: x = grndtoi(gmul(dx, gauss_realimag(nf,col)), pe);
997: return (*pe > -5)? NULL: gdiv(x, dx);
998: }
999:
1000: /* y = C \prod g[i]^e[i] ? */
1001: static int
1002: fact_ok(GEN nf, GEN y, GEN C, GEN g, GEN e)
1003: {
1004: long i, c = lg(e);
1005: GEN z = C? C: gun;
1006: for (i=1; i<c; i++)
1007: if (signe(e[i])) z = idealmul(nf, z, idealpow(nf, (GEN)g[i], (GEN)e[i]));
1008: if (typ(z) != t_MAT) z = idealhermite(nf,z);
1009: if (typ(y) != t_MAT) y = idealhermite(nf,y);
1010: return gegal(y, z);
1011: }
1012:
1013: /* assume x in HNF. cf class_group_gen for notations */
1014: static GEN
1015: isprincipalall0(GEN bnf, GEN x, long *ptprec, long flag)
1016: {
1017: long i,lW,lB,e,c, prec = *ptprec;
1018: GEN Q,xar,Wex,Bex,U,y,p1,gen,cyc,xc,ex,d,col,A;
1019: GEN W = (GEN)bnf[1];
1020: GEN B = (GEN)bnf[2];
1021: GEN WB_C = (GEN)bnf[4];
1022: GEN vperm = (GEN)bnf[6];
1023: GEN nf = (GEN)bnf[7];
1024: GEN clg2 = (GEN)bnf[9];
1025: int old_format = (typ(clg2[2]) == t_MAT);
1026:
1027: U = (GEN)clg2[1];
1028: cyc = gmael3(bnf,8,1,2); c = lg(cyc)-1;
1029: gen = gmael3(bnf,8,1,3);
1030:
1031: vectbase = (GEN)bnf[5]; /* needed by factorgensimple */
1032:
1033: /* factor x */
1034: xc = content(x); if (!gcmp1(xc)) x = gdiv(x,xc);
1035:
1036: p1 = init_idele(nf); p1[1] = (long)x;
1037: if (!(flag & nf_GEN)) p1[2] = 0; /* don't compute archimedean part */
1038: xar = split_ideal(nf,p1,prec,vperm);
1039:
1040: lW = lg(W)-1; Wex = zerocol(lW);
1041: lB = lg(B)-1; Bex = zerocol(lB); get_split_expo(Wex,Bex,vperm);
1042:
1043: /* x = g_W Wex + g_B Bex - [xar]
1044: * = g_W A - [xar] + [C_B]Bex since g_W B + g_B = [C_B] */
1045: A = gsub(Wex, gmul(B,Bex));
1046: if (old_format) U = ginv(U);
1047: Q = gmul(U, A);
1048: ex = cgetg(c+1,t_COL);
1049: for (i=1; i<=c; i++)
1050: Q[i] = (long)truedvmdii((GEN)Q[i],(GEN)cyc[i],(GEN*)(ex+i));
1051: if (!(flag & nf_GEN)) return gcopy(ex);
1052:
1053: /* compute arch component of the missing principal ideal */
1054: if (old_format)
1055: {
1056: GEN Garch, V = (GEN)clg2[2];
1057: p1 = c? concatsp(gmul(V,Q), Bex): Bex;
1058: col = act_arch(p1, WB_C);
1059: if (c)
1060: {
1061: Garch = get_Garch(nf,gen,clg2,prec);
1062: col = gadd(col, act_arch(ex, Garch));
1063: }
1064: }
1065: else
1066: { /* g A = G Ur A + [ga]A, Ur A = D Q + R as above (R = ex)
1067: = G R + [GD]Q + [ga]A */
1068: GEN ga = (GEN)clg2[2], GD = (GEN)clg2[3];
1069: col = act_arch(Bex, WB_C + lW);
1070: if (lW) col = gadd(col, act_arch(A, ga));
1071: if (c) col = gadd(col, act_arch(Q, GD));
1072: }
1073: col = gsub(col, (GEN)xar[2]);
1074:
1075: /* find coords on Zk; Q = N (x / \prod gj^ej) = N(alpha), denom(alpha) | d */
1076: Q = gdiv(dethnf_i(x), get_norm_fact(gen, ex, &d));
1077: col = isprincipalarch(bnf, col, Q, gun, d, &e);
1078: if (col && !fact_ok(nf,x, col,gen,ex)) col = NULL;
1079: if (!col)
1080: {
1081: *ptprec = prec + (e >> TWOPOTBITS_IN_LONG) + (MEDDEFAULTPREC-2);
1082: if (flag & nf_FORCE)
1083: {
1084: if (DEBUGLEVEL) err(warner,"precision too low for generators, e = %ld",e);
1085: return NULL;
1086: }
1087: err(warner,"precision too low for generators, not given");
1088: e = 0;
1089: }
1090: y = cgetg(4,t_VEC);
1091: y[1] = lcopy(ex);
1092: y[2] = e? lmul(xc,col): lgetg(1,t_COL);
1093: y[3] = lstoi(-e); return y;
1094: }
1095:
1096: static GEN
1097: triv_gen(GEN nf, GEN x, long c, long flag)
1098: {
1099: GEN y;
1100: if (!(flag & nf_GEN)) return zerocol(c);
1101: y = cgetg(4,t_VEC);
1102: y[1] = (long)zerocol(c);
1103: x = (typ(x) == t_COL)? gcopy(x): algtobasis(nf,x);
1104: y[2] = (long)x;
1105: y[3] = lstoi(BIGINT); return y;
1106: }
1107:
1108: GEN
1109: isprincipalall(GEN bnf,GEN x,long flag)
1110: {
1111: long av = avma,c,pr, tx = typ(x);
1112: GEN nf;
1113:
1114: bnf = checkbnf(bnf); nf = (GEN)bnf[7];
1115: if (tx == t_POLMOD)
1116: {
1117: if (!gegal((GEN)x[1],(GEN)nf[1]))
1118: err(talker,"not the same number field in isprincipal");
1119: x = (GEN)x[2]; tx = t_POL;
1120: }
1121: if (tx == t_POL || tx == t_COL)
1122: {
1123: if (gcmp0(x)) err(talker,"zero ideal in isprincipal");
1124: return triv_gen(nf, x, lg(mael3(bnf,8,1,2))-1, flag);
1125: }
1126: x = idealhermite(nf,x);
1127: if (lg(x)==1) err(talker,"zero ideal in isprincipal");
1128: if (lgef(nf[1])==4)
1129: return gerepileupto(av, triv_gen(nf, gcoeff(x,1,1), 0, flag));
1130:
1131: pr = prec_arch(bnf); /* precision of unit matrix */
1132: c = getrand();
1133: for (;;)
1134: {
1135: long av1 = avma;
1136: GEN y = isprincipalall0(bnf,x,&pr,flag);
1137: if (y) return gerepileupto(av,y);
1138:
1139: if (DEBUGLEVEL) err(warnprec,"isprincipalall0",pr);
1140: avma = av1; bnf = bnfnewprec(bnf,pr); setrand(c);
1141: }
1142: }
1143:
1144: /* isprincipal for C * \prod P[i]^e[i] (C omitted if NULL) */
1145: GEN
1146: isprincipalfact(GEN bnf,GEN P, GEN e, GEN C, long flag)
1147: {
1148: long av = avma, l = lg(e), i,prec,c;
1149: GEN id,id2, nf = checknf(bnf), z = NULL; /* gcc -Wall */
1150: int gen = flag & (nf_GEN | nf_GENMAT);
1151:
1152: prec = prec_arch(bnf);
1153: if (gen)
1154: {
1155: z = cgetg(3,t_VEC);
1156: z[2] = (flag & nf_GENMAT)? lgetg(1, t_MAT): lmodulcp(gun,(GEN)nf[1]);
1157: }
1158: id = C;
1159: for (i=1; i<l; i++) /* compute prod P[i]^e[i] */
1160: if (signe(e[i]))
1161: {
1162: if (gen) z[1] = P[i]; else z = (GEN)P[i];
1163: id2 = idealpowred(bnf,z, (GEN)e[i],prec);
1164: id = id? idealmulred(nf,id,id2,prec): id2;
1165: }
1166: if (id == C)
1167: {
1168: if (!C) id = gun;
1169: return isprincipalall(bnf, id, flag);
1170: }
1171: c = getrand();
1172: for (;;)
1173: {
1174: long av1 = avma;
1175: GEN y = isprincipalall0(bnf, gen? (GEN)id[1]: id,&prec,flag);
1176: if (y)
1177: {
1178: if (gen && typ(y) == t_VEC)
1179: {
1180: GEN u = lift_intern(basistoalg(nf, (GEN)y[2]));
1181: if (flag & nf_GENMAT)
1182: y[2] = (gcmp1(u)&&lg(id[2])>1)? id[2]: (long)arch_mul((GEN)id[2], u);
1183: else
1184: y[2] = (long)algtobasis(nf, gmul((GEN)id[2], u));
1185: y = gcopy(y);
1186: }
1187: return gerepileupto(av,y);
1188: }
1189:
1190: if (flag & nf_GIVEPREC)
1191: {
1192: if (DEBUGLEVEL)
1193: err(warner,"insufficient precision for generators, not given");
1194: avma = av; return stoi(prec);
1195: }
1196: if (DEBUGLEVEL) err(warnprec,"isprincipalall0",prec);
1197: avma = av1; bnf = bnfnewprec(bnf,prec); setrand(c);
1198: }
1199: }
1200:
1201: GEN
1202: isprincipal(GEN bnf,GEN x)
1203: {
1204: return isprincipalall(bnf,x,nf_REGULAR);
1205: }
1206:
1207: GEN
1208: isprincipalgen(GEN bnf,GEN x)
1209: {
1210: return isprincipalall(bnf,x,nf_GEN);
1211: }
1212:
1213: GEN
1214: isprincipalforce(GEN bnf,GEN x)
1215: {
1216: return isprincipalall(bnf,x,nf_FORCE);
1217: }
1218:
1219: GEN
1220: isprincipalgenforce(GEN bnf,GEN x)
1221: {
1222: return isprincipalall(bnf,x,nf_GEN | nf_FORCE);
1223: }
1224:
1225: GEN
1226: isunit(GEN bnf,GEN x)
1227: {
1228: long av=avma,tetpil,tx = typ(x),i,R1,RU,n;
1229: GEN p1,logunit,y,ex,nf,z,pi2_sur_w,gn,emb;
1230:
1231: bnf = checkbnf(bnf); nf=(GEN)bnf[7];
1232: logunit=(GEN)bnf[3]; RU=lg(logunit);
1233: p1 = gmael(bnf,8,4); /* roots of 1 */
1234: gn = (GEN)p1[1]; n = itos(gn);
1235: z = (GEN)p1[2];
1236: switch(tx)
1237: {
1238: case t_INT: case t_FRAC: case t_FRACN:
1239: if (!gcmp1(x) && !gcmp_1(x)) return cgetg(1,t_COL);
1240: y = zerocol(RU); i = (gsigne(x) > 0)? 0: n>>1;
1241: y[RU] = (long)gmodulss(i, n); return y;
1242:
1243: case t_POLMOD:
1244: if (!gegal((GEN)nf[1],(GEN)x[1]))
1245: err(talker,"not the same number field in isunit");
1246: x = (GEN)x[2]; /* fall through */
1247: case t_POL:
1248: p1 = x; x = algtobasis(bnf,x); break;
1249:
1250: case t_COL:
1251: if (lgef(nf[1])-2 == lg(x)) { p1 = basistoalg(nf,x); break; }
1252:
1253: default:
1254: err(talker,"not an algebraic number in isunit");
1255: }
1256: if (!gcmp1(denom(x))) { avma = av; return cgetg(1,t_COL); }
1257: if (typ(p1) != t_POLMOD) p1 = gmodulcp(p1,(GEN)nf[1]);
1258: p1 = gnorm(p1);
1259: if (!is_pm1(p1)) { avma = av; return cgetg(1,t_COL); }
1260:
1261: R1 = itos(gmael(nf,2,1)); p1 = cgetg(RU+1,t_COL);
1262: for (i=1; i<=R1; i++) p1[i] = un;
1263: for ( ; i<=RU; i++) p1[i] = deux;
1264: logunit = concatsp(logunit,p1);
1265: /* ex = fundamental units exponents */
1266: {
1267: GEN rx, rlog = greal(logunit);
1268: long e, prec = nfgetprec(nf);
1269: for (i=1;;)
1270: {
1271: rx = get_arch_real(nf,x,&emb, MEDDEFAULTPREC);
1272: if (rx)
1273: {
1274: ex = grndtoi(gauss(rlog, rx), &e);
1275: if (gcmp0((GEN)ex[RU]) && e < -4) break;
1276: }
1277:
1278: if (++i > 4) err(precer,"isunit");
1279: prec = (prec-1)<<1;
1280: if (DEBUGLEVEL) err(warnprec,"isunit",prec);
1281: nf = nfnewprec(nf, prec);
1282: }
1283: }
1284:
1285: setlg(ex, RU);
1286: setlg(p1, RU); settyp(p1, t_VEC);
1287: for (i=1; i<RU; i++) p1[i] = coeff(logunit, 1, i);
1288: p1 = gneg(gimag(gmul(p1,ex))); if (!R1) p1 = gmul2n(p1, -1);
1289: p1 = gadd(garg((GEN)emb[1],DEFAULTPREC), p1);
1290: /* p1 = arg(the missing root of 1) */
1291:
1292: pi2_sur_w = divrs(mppi(DEFAULTPREC), n>>1);
1293: p1 = ground(gdiv(p1, pi2_sur_w));
1294: if (n > 2)
1295: {
1296: GEN ro = gmael(nf,6,1);
1297: GEN p2 = ground(gdiv(garg(poleval(z,ro), DEFAULTPREC), pi2_sur_w));
1298: p1 = mulii(p1, mpinvmod(p2, gn));
1299: }
1300:
1301: tetpil = avma; y = cgetg(RU+1,t_COL);
1302: for (i=1; i<RU; i++) y[i] = lcopy((GEN)ex[i]);
1303: y[RU] = lmodulcp(p1, gn); return gerepile(av,tetpil,y);
1304: }
1305:
1306: GEN
1307: signunits(GEN bnf)
1308: {
1309: long av,i,j,R1,RU,mun;
1310: GEN matunit,y,p1,p2,nf,pi;
1311:
1312: bnf = checkbnf(bnf); nf = (GEN)bnf[7];
1313: matunit = (GEN)bnf[3]; RU = lg(matunit);
1314: R1=itos(gmael(nf,2,1)); pi=mppi(MEDDEFAULTPREC);
1315: y=cgetg(RU,t_MAT); mun = lnegi(gun);
1316: for (j=1; j<RU; j++)
1317: {
1318: p1=cgetg(R1+1,t_COL); y[j]=(long)p1; av=avma;
1319: for (i=1; i<=R1; i++)
1320: {
1321: p2 = ground(gdiv(gimag(gcoeff(matunit,i,j)),pi));
1322: p1[i] = mpodd(p2)? mun: un;
1323: }
1324: avma=av;
1325: }
1326: return y;
1327: }
1328:
1329: /* LLL-reduce ideal and return T2 | ideal */
1330: static GEN
1331: red_ideal(GEN *ideal,GEN T2vec,GEN prvec)
1332: {
1333: long i;
1334: for (i=1; i<lg(T2vec); i++)
1335: {
1336: GEN p1,base, T2 = (GEN)T2vec[i];
1337: long prec = prvec[i];
1338:
1339: p1 = qf_base_change(T2, *ideal, 1);
1340: base = lllgramintern(p1,100,1,prec);
1341: if (base)
1342: {
1343: p1 = sqred1intern(qf_base_change(p1,base,1),1);
1344: if (p1) { *ideal = gmul(*ideal,base); return p1; }
1345: }
1346: if (DEBUGLEVEL) err(warner, "prec too low in red_ideal[%ld]: %ld",i,prec);
1347: }
1348: return NULL;
1349: }
1350:
1351: static double
1352: get_minkovski(long N, long R1, GEN D, GEN gborne)
1353: {
1354: const long prec = DEFAULTPREC;
1355: GEN p1,p2, pi = mppi(prec);
1356: double bound;
1357:
1358: p1 = gsqrt(gmulsg(N,gmul2n(pi,1)),prec);
1359: p1 = gdiv(p1, gexp(stoi(N),prec));
1360: p1 = gmulsg(N, gpui(p1, dbltor(2./(double)N),prec));
1361: p2 = gpui(gdivsg(4,pi), gdivgs(stoi(N-R1),N),prec);
1362: p1 = gmul(p1,p2);
1363: bound = gtodouble(gmul(p1, gpui(absi(D), dbltor(1./(double)N),prec)));
1364: bound *= gtodouble(gborne);
1365: if (DEBUGLEVEL) { fprintferr("Bound for norms = %.0f\n",bound); flusherr(); }
1366: return bound;
1367: }
1368:
1369: static void
1370: wr_rel(long *col)
1371: {
1372: long i;
1373: fprintferr("\nrel = ");
1374: for (i=1; i<=KC; i++)
1375: if (col[i]) fprintferr("%ld^%ld ",i,col[i]);
1376: fprintferr("\n");
1377: }
1378:
1379: static void
1380: set_log_embed(GEN col, GEN x, long R1, long prec)
1381: {
1382: long i, l = lg(x);
1383: for (i=1; i<=R1; i++) gaffect( glog((GEN)x[i],prec) , (GEN)col[i]);
1384: for ( ; i < l; i++) gaffect(gmul2n(glog((GEN)x[i],prec), 1), (GEN)col[i]);
1385: }
1386:
1387: static void
1388: set_fact(GEN c)
1389: {
1390: long i;
1391: c[0] = primfact[1]; /* for already_found_relation */
1392: for (i=1; i<=primfact[0]; i++) c[primfact[i]] = expoprimfact[i];
1393: }
1394:
1395: static void
1396: unset_fact(GEN c)
1397: {
1398: long i;
1399: for (i=1; i<=primfact[0]; i++) c[primfact[i]] = 0;
1400: }
1401:
1402: #define MAXTRY 20
1403:
1404: /* return -1 in case of precision problems. t = current number of relations */
1405: static long
1406: small_norm_for_buchall(long cglob,GEN *mat,GEN matarch,long LIMC, long PRECREG,
1407: GEN nf,GEN gborne,long nbrelpid,GEN invp,GEN L,GEN LLLnf)
1408: {
1409: const double eps = 0.000001;
1410: double *y,*z,**q,*v, MINKOVSKI_BOUND,BOUND;
1411: ulong av = avma, av1,av2,limpile;
1412: long i,j,k,noideal, nbrel = lg(mat)-1;
1413: long alldep = 0, nbsmallnorm,nbfact,R1, N = degpol(nf[1]);
1414: GEN x,xembed,M,T2,r,cbase,invcbase,T2vec,prvec;
1415:
1416: if (gsigne(gborne)<=0) return cglob;
1417: if (DEBUGLEVEL)
1418: fprintferr("\n#### Looking for %ld relations (small norms)\n",nbrel);
1419: xembed = NULL; /* gcc -Wall */
1420: nbsmallnorm = nbfact = 0;
1421: R1 = itos(gmael(nf,2,1));
1422: T2 = (GEN)LLLnf[1];
1423: M = (GEN)LLLnf[2];
1424: cbase =(GEN)LLLnf[3];
1425: invcbase=(GEN)LLLnf[4];
1426:
1427: prvec = cgetg(3,t_VECSMALL);
1428: prvec[1] = MEDDEFAULTPREC;
1429: prvec[2] = PRECREG;
1430: T2vec = cgetg(3,t_VEC);
1431: T2vec[1] = (long)gprec_w(T2,prvec[1]);
1432: T2vec[2] = (long)T2;
1433: minim_alloc(N+1, &q, &x, &y, &z, &v);
1434: av1 = avma;
1435: MINKOVSKI_BOUND = get_minkovski(N,R1,(GEN)nf[3],gborne);
1436: for (noideal=KC; noideal; noideal--)
1437: {
1438: long nbrelideal=0, dependent = 0;
1439: GEN IDEAL, ideal = (GEN)vectbase[noideal];
1440: GEN normideal = idealnorm(nf,ideal);
1441:
1442: if (alldep > 2*N) break;
1443:
1444: if (DEBUGLEVEL>1) fprintferr("\n*** Ideal no %ld: %Z\n", noideal, ideal);
1445: ideal = prime_to_ideal(nf,ideal);
1446: IDEAL = invcbase? gmul(invcbase, ideal): ideal;
1447: IDEAL = gmul(IDEAL, lllint(IDEAL)); /* should be almost T2-reduced */
1448: r = red_ideal(&IDEAL,T2vec,prvec);
1449: if (!r) return -1; /* precision problem */
1450:
1451: for (k=1; k<=N; k++)
1452: {
1453: v[k]=gtodouble(gcoeff(r,k,k));
1454: for (j=1; j<k; j++) q[j][k]=gtodouble(gcoeff(r,j,k));
1455: if (DEBUGLEVEL>3) fprintferr("v[%ld]=%.0f ",k,v[k]);
1456: }
1457:
1458: BOUND = MINKOVSKI_BOUND * pow(gtodouble(normideal),2./(double)N);
1459: if (DEBUGLEVEL>1)
1460: {
1461: if (DEBUGLEVEL>3) fprintferr("\n");
1462: fprintferr("BOUND = %.0f\n",BOUND); flusherr();
1463: }
1464: BOUND += eps; av2=avma; limpile = stack_lim(av2,1);
1465: k = N; y[N]=z[N]=0; x[N]= (long) sqrt(BOUND/v[N]);
1466: for(;; x[1]--)
1467: {
1468: ulong av3 = avma;
1469: double p;
1470: GEN col;
1471:
1472: for(;;) /* looking for primitive element of small norm */
1473: { /* cf minim00 */
1474: if (k>1)
1475: {
1476: long l = k-1;
1477: z[l] = 0;
1478: for (j=k; j<=N; j++) z[l] += q[l][j]*x[j];
1479: p = x[k]+z[k];
1480: y[l] = y[k]+p*p*v[k];
1481: x[l] = (long) floor(sqrt((BOUND-y[l])/v[l])-z[l]);
1482: k = l;
1483: }
1484: for(;;)
1485: {
1486: p = x[k]+z[k];
1487: if (y[k] + p*p*v[k] <= BOUND) break;
1488: k++; x[k]--;
1489: }
1490: if (k==1) /* element complete */
1491: {
1492: if (!x[1] && y[1]<=eps) goto ENDIDEAL;
1493: if (ccontent(x)==1) /* primitive */
1494: {
1495: GEN Nx, gx = gmul_mati_smallvec(IDEAL,x);
1496: long av4;
1497: if (!isnfscalar(gx))
1498: {
1499: xembed = gmul(M,gx); av4 = avma; nbsmallnorm++;
1500: Nx = ground(norm_by_embed(R1,xembed)); setsigne(Nx, 1);
1501: if (factorelt(nf,cbase,gx,Nx,KCZ,LIMC)) { avma = av4; break; }
1502: if (DEBUGLEVEL > 1) { fprintferr("."); flusherr(); }
1503: }
1504: avma = av3;
1505: }
1506: x[1]--;
1507: }
1508: }
1509: cglob++; col = mat[cglob];
1510: set_fact(col);
1511: if (cglob > 1 && cglob <= KC && ! addcolumntomatrix(col,invp,L))
1512: { /* Q-dependent from previous ones: forget it */
1513: cglob--; unset_fact(col);
1514: if (DEBUGLEVEL>1) { fprintferr("*"); flusherr(); }
1515: if (++dependent > MAXTRY) { alldep++; break; }
1516: avma = av3; continue;
1517: }
1518: /* record archimedean part */
1519: set_log_embed((GEN)matarch[cglob], xembed, R1, PRECREG);
1520: alldep = dependent = 0;
1521:
1522: if (DEBUGLEVEL)
1523: {
1524: if (DEBUGLEVEL==1) fprintferr("%4ld",cglob);
1525: else { fprintferr("cglob = %ld. ",cglob); wr_rel(mat[cglob]); }
1526: flusherr(); nbfact++;
1527: }
1528: if (cglob >= nbrel) goto END; /* we have enough */
1529: if (++nbrelideal == nbrelpid) break;
1530:
1531: if (low_stack(limpile, stack_lim(av2,1)))
1532: {
1533: if(DEBUGMEM>1) err(warnmem,"small_norm");
1534: invp = gerepilecopy(av2, invp);
1535: }
1536: }
1537: ENDIDEAL:
1538: invp = gerepilecopy(av1, invp);
1539: if (DEBUGLEVEL>1) msgtimer("for this ideal");
1540: }
1541: END:
1542: if (DEBUGLEVEL)
1543: {
1544: fprintferr("\n");
1545: msgtimer("small norm relations");
1546: if (DEBUGLEVEL>1)
1547: {
1548: GEN p1,tmp=cgetg(cglob+1,t_MAT);
1549:
1550: fprintferr("Elements of small norm gave %ld relations.\n",cglob);
1551: fprintferr("Computing rank: "); flusherr();
1552: for(j=1;j<=cglob;j++)
1553: {
1554: p1=cgetg(KC+1,t_COL); tmp[j]=(long)p1;
1555: for(i=1;i<=KC;i++) p1[i]=lstoi(mat[j][i]);
1556: }
1557: tmp = (GEN)sindexrank(tmp)[2]; k=lg(tmp)-1;
1558: fprintferr("%ld; independent columns: %Z\n",k,tmp);
1559: }
1560: if(nbsmallnorm)
1561: fprintferr("\nnb. fact./nb. small norm = %ld/%ld = %.3f\n",
1562: nbfact,nbsmallnorm,((double)nbfact)/nbsmallnorm);
1563: else
1564: fprintferr("\nnb. small norm = 0\n");
1565: }
1566: avma = av; return cglob;
1567: }
1568: #undef MAXTRY
1569:
1570: /* I assumed to be integral HNF, T2 a weighted T2 matrix */
1571: static GEN
1572: ideallllredpart1(GEN I, GEN T2, long prec)
1573: {
1574: GEN y,m,idealpro;
1575:
1576: y = lllgramintern(qf_base_change(T2,I,1),100,1,prec+1);
1577: if (!y) return NULL;
1578:
1579: /* I, m, y integral */
1580: m = gmul(I,(GEN)y[1]);
1581: if (isnfscalar(m)) m = gmul(I,(GEN)y[2]);
1582:
1583: idealpro = cgetg(3,t_VEC);
1584: idealpro[1] = (long)I;
1585: idealpro[2] = (long)m; /* irrational element of small T2 norm in I */
1586: if (DEBUGLEVEL>5) fprintferr("\nidealpro = %Z\n",idealpro);
1587: return idealpro;
1588: }
1589:
1590: static void
1591: dbg_newrel(long jideal, long jdir, long phase, long cglob, long *col,
1592: GEN colarch, long lim)
1593: {
1594: fprintferr("\n++++ cglob = %ld: new relation (need %ld)", cglob, lim);
1595: wr_rel(col);
1596: if (DEBUGLEVEL>3)
1597: {
1598: fprintferr("(jideal=%ld,jdir=%ld,phase=%ld)", jideal,jdir,phase);
1599: msgtimer("for this relation");
1600: }
1601: if (DEBUGLEVEL>6) fprintferr("archimedian part = %Z\n",colarch);
1602: flusherr() ;
1603: }
1604:
1605: static void
1606: dbg_cancelrel(long jideal,long jdir,long phase, long *col)
1607: {
1608: fprintferr("rel. cancelled. phase %ld: %ld ",phase);
1609: if (DEBUGLEVEL>3) fprintferr("(jideal=%ld,jdir=%ld)",jideal,jdir);
1610: wr_rel(col); flusherr();
1611: }
1612:
1613: static void
1614: dbg_outrel(long phase,long cglob, GEN vperm,GEN *mat,GEN maarch)
1615: {
1616: long i,j;
1617: GEN p1,p2;
1618:
1619: if (phase == 0)
1620: {
1621: ulong av = avma; p2=cgetg(cglob+1,t_MAT);
1622: for (j=1; j<=cglob; j++)
1623: {
1624: p1=cgetg(KC+1,t_COL); p2[j]=(long)p1;
1625: for (i=1; i<=KC; i++) p1[i]=lstoi(mat[j][i]);
1626: }
1627: fprintferr("\nRank = %ld\n", rank(p2));
1628: if (DEBUGLEVEL>3)
1629: {
1630: fprintferr("relations = \n");
1631: for (j=1; j <= cglob; j++) wr_rel(mat[j]);
1632: fprintferr("\nmatarch = %Z\n",maarch);
1633: }
1634: avma = av;
1635: }
1636: else if (DEBUGLEVEL>6)
1637: {
1638: fprintferr("before hnfadd:\nvectbase[vperm[]] = \n");
1639: fprintferr("[");
1640: for (i=1; i<=KC; i++)
1641: {
1642: bruterr((GEN)vectbase[vperm[i]],'g',-1);
1643: if (i<KC) fprintferr(",");
1644: }
1645: fprintferr("]~\n");
1646: }
1647: flusherr();
1648: }
1649:
1650: /* Check if we already have a column mat[l] equal to mat[s]
1651: * General check for colinearity useless since exceedingly rare */
1652: static long
1653: already_found_relation(GEN *mat,long s)
1654: {
1655: GEN coll, cols = mat[s];
1656: long l,bs;
1657:
1658: bs = 1; while (bs<=KC && !cols[bs]) bs++;
1659: if (bs > KC) return s; /* zero relation */
1660:
1661: for (l=s-1; l; l--)
1662: {
1663: coll = mat[l];
1664: if (bs == coll[0]) /* = index of first non zero elt in coll */
1665: {
1666: long b = bs;
1667: do b++; while (b<=KC && cols[b] == coll[b]);
1668: if (b > KC) return l;
1669: }
1670: }
1671: cols[0] = bs; return 0;
1672: }
1673:
1674: /* I integral ideal in HNF form */
1675: static GEN
1676: remove_content(GEN I)
1677: {
1678: long N = lg(I)-1;
1679: if (!gcmp1(gcoeff(I,N,N))) { GEN y=content(I); if (!gcmp1(y)) I=gdiv(I,y); }
1680: return I;
1681: }
1682:
1683: /* if phase != 1 re-initialize static variables. If <0 return immediately */
1684: static long
1685: random_relation(long phase,long cglob,long LIMC,long PRECLLL,long PRECREG,
1686: GEN nf,GEN subFB,GEN vecT2,GEN *mat,GEN matarch,GEN list_jideal)
1687: {
1688: static long jideal, jdir;
1689: long lim,i,av,av1,cptzer,nbT2,lgsub,r1, jlist = 1;
1690: GEN arch,col,colarch,ideal,idealpro,P,ex;
1691:
1692: if (phase != 1) { jideal=jdir=1; if (phase<0) return 0; }
1693:
1694: r1 = nf_get_r1(nf);
1695: lim = lg(mat)-1; /* requested number of relations */
1696: nbT2 = lg(vecT2)-1;
1697: lgsub = lg(subFB); ex = cgetg(lgsub, t_VECSMALL);
1698: cptzer = 0;
1699: if (DEBUGLEVEL && list_jideal)
1700: fprintferr("looking hard for %Z\n",list_jideal);
1701: P = NULL; /* gcc -Wall */
1702: for (av = avma;;)
1703: {
1704: if (list_jideal && jlist < lg(list_jideal) && jdir <= nbT2)
1705: jideal = list_jideal[jlist++];
1706: if (!list_jideal || jdir <= nbT2)
1707: {
1708: avma = av;
1709: P = prime_to_ideal(nf, (GEN)vectbase[jideal]);
1710: }
1711: ideal = P;
1712: do {
1713: for (i=1; i<lgsub; i++)
1714: {
1715: ex[i] = mymyrand()>>randshift;
1716: if (ex[i]) ideal = idealmulh(nf,ideal, gmael3(powsubFB,i,ex[i],1));
1717: }
1718: }
1719: while (ideal == P); /* If ex = 0, try another */
1720: ideal = remove_content(ideal);
1721:
1722: if (phase != 1) jdir = 1; else phase = 2;
1723: for (av1 = avma; jdir <= nbT2; jdir++, avma = av1)
1724: { /* reduce along various directions */
1725: if (DEBUGLEVEL>2)
1726: fprintferr("phase=%ld,jideal=%ld,jdir=%ld,rand=%ld\n",
1727: phase,jideal,jdir,getrand());
1728: idealpro = ideallllredpart1(ideal,(GEN)vecT2[jdir],PRECLLL);
1729: if (!idealpro) return -2;
1730: if (!factorgen(nf,idealpro,KCZ,LIMC))
1731: {
1732: if (DEBUGLEVEL>1) { fprintferr("."); flusherr(); }
1733: continue;
1734: }
1735: /* can factor ideal, record relation */
1736: cglob++; col = mat[cglob];
1737: set_fact(col); col[jideal]--;
1738: for (i=1; i<lgsub; i++) col[subFB[i]] -= ex[i];
1739: if (already_found_relation(mat, cglob))
1740: { /* Already known: forget it */
1741: if (DEBUGLEVEL>1) dbg_cancelrel(jideal,jdir,phase,col);
1742: cglob--; unset_fact(col); col[jideal] = 0;
1743: for (i=1; i<lgsub; i++) col[subFB[i]] = 0;
1744:
1745: if (++cptzer > MAXRELSUP)
1746: {
1747: if (list_jideal) { cptzer -= 10; break; }
1748: return -1;
1749: }
1750: continue;
1751: }
1752:
1753: /* Compute and record archimedian part */
1754: cptzer = 0; arch = NULL;
1755: for (i=1; i<lgsub; i++)
1756: if (ex[i])
1757: {
1758: GEN p1 = gmael3(powsubFB,i,ex[i],2);
1759: arch = arch? vecmul(arch, p1): p1;
1760: }
1761: colarch = (GEN)matarch[cglob];
1762: /* arch = archimedean component (MULTIPLICATIVE form) of ideal */
1763: arch = vecdiv(arch, gmul(gmael(nf,5,1), (GEN)idealpro[2]));
1764: set_log_embed(colarch, arch, r1, PRECREG);
1765: if (DEBUGLEVEL) dbg_newrel(jideal,jdir,phase,cglob,col,colarch,lim);
1766:
1767: /* Need more, try next P */
1768: if (cglob < lim) break;
1769:
1770: /* We have found enough. Return */
1771: if (phase)
1772: {
1773: jdir = 1;
1774: if (jideal == KC) jideal=1; else jideal++;
1775: }
1776: if (DEBUGLEVEL>2)
1777: fprintferr("Upon exit: jideal=%ld,jdir=%ld\n",jideal,jdir);
1778: avma = av; return cglob;
1779: }
1780: if (!list_jideal)
1781: {
1782: if (jideal == KC) jideal=1; else jideal++;
1783: }
1784: }
1785: }
1786:
1787: static long
1788: be_honest(GEN nf,GEN subFB,long PRECLLL)
1789: {
1790: ulong av;
1791: GEN MC = gmael(nf,5,2), M = gmael(nf,5,1), D = (GEN)nf[3];
1792: long ex,i,j,J,k,iz,nbtest, lgsub = lg(subFB), ru = lg(MC);
1793: GEN P,ideal,idealpro, exu = cgetg(ru, t_VECSMALL), MCtw= cgetg(ru, t_MAT);
1794:
1795: if (DEBUGLEVEL)
1796: {
1797: fprintferr("Be honest for primes from %ld to %ld\n", FB[KCZ+1],FB[KCZ2]);
1798: flusherr();
1799: }
1800: av = avma;
1801: for (iz=KCZ+1; iz<=KCZ2; iz++, avma = av)
1802: {
1803: if (DEBUGLEVEL>1) fprintferr("%ld ", FB[iz]);
1804: P = idealbase[numFB[FB[iz]]]; J = lg(P);
1805: /* if unramified, check all but 1 */
1806: if (J > 1 && !divise(D, gmael(P,1,1))) J--;
1807: for (j=1; j<J; j++)
1808: {
1809: GEN ideal0 = prime_to_ideal(nf,(GEN)P[j]);
1810: ulong av2 = avma;
1811: for(nbtest=0;;)
1812: {
1813: ideal = ideal0;
1814: for (i=1; i<lgsub; i++)
1815: {
1816: ex = mymyrand()>>randshift;
1817: if (ex) ideal = idealmulh(nf,ideal,gmael3(powsubFB,i,ex,1));
1818: }
1819: ideal = remove_content(ideal);
1820: for (k=1; k<ru; k++)
1821: {
1822: if (k==1)
1823: for (i=1; i<ru; i++) exu[i] = mymyrand()>>randshift;
1824: else
1825: {
1826: for (i=1; i<ru; i++) exu[i] = 0;
1827: exu[k] = 10;
1828: }
1829: for (i=1; i<ru; i++)
1830: MCtw[i] = exu[i]? lmul2n((GEN)MC[i],exu[i]<<1): MC[i];
1831: idealpro = ideallllredpart1(ideal, mulmat_real(MCtw,M), PRECLLL);
1832: if (idealpro && factorgen(nf,idealpro,iz-1,FB[iz-1])) break;
1833: nbtest++; if (nbtest==200) return 0;
1834: }
1835: avma = av2; if (k < ru) break;
1836: }
1837: }
1838: }
1839: if (DEBUGLEVEL)
1840: {
1841: if (DEBUGLEVEL>1) fprintferr("\n");
1842: msgtimer("be honest");
1843: }
1844: avma = av; return 1;
1845: }
1846:
1847: int
1848: trunc_error(GEN x)
1849: {
1850: return typ(x)==t_REAL && signe(x)
1851: && (expo(x)>>TWOPOTBITS_IN_LONG) + 3 > lg(x);
1852: }
1853:
1854: /* xarch = complex logarithmic embeddings of units (u_j) found so far */
1855: static GEN
1856: compute_multiple_of_R(GEN xarch,long RU,long N,GEN *ptsublambda)
1857: {
1858: GEN v,mdet,mdet_t,Im_mdet,kR,sublambda,lambda,xreal;
1859: GEN *gptr[2];
1860: long av = avma, i,j, sreg = lg(xarch)-1, R1 = 2*RU - N;
1861:
1862: if (DEBUGLEVEL) { fprintferr("\n#### Computing regulator\n"); flusherr(); }
1863: xreal=greal(xarch); v=cgetg(RU+1,t_COL); /* xreal = (log |sigma_i(u_j)|) */
1864: for (i=1; i<=R1; i++) v[i]=un;
1865: for ( ; i<=RU; i++) v[i]=deux;
1866: mdet=cgetg(sreg+2,t_MAT); mdet[1]=(long)v;
1867: for (j=2; j<=sreg+1; j++) mdet[j]=xreal[j-1]; /* det(Span(mdet)) = N * R */
1868:
1869: i = gprecision(mdet); /* truncate to avoid "near dependant" vectors */
1870: mdet_t = (i <= 4)? mdet: gprec_w(mdet,i-1);
1871: v = (GEN)sindexrank(mdet_t)[2]; /* list of independant column indices */
1872: /* check we have full rank for units */
1873: if (lg(v) != RU+1) { avma=av; return NULL; }
1874:
1875: Im_mdet = vecextract_p(mdet,v);
1876: /* integral multiple of R: the cols we picked form a Q-basis, they have an
1877: * index in the full lattice */
1878: kR = gdivgs(det2(Im_mdet), N);
1879: /* R > 0.2 uniformly */
1880: if (gexpo(kR) < -3) { avma=av; return NULL; }
1881:
1882: kR = mpabs(kR);
1883: sublambda = cgetg(sreg+1,t_MAT);
1884: lambda = gauss(Im_mdet,xreal); /* rational entries */
1885: for (i=1; i<=sreg; i++)
1886: {
1887: GEN p1 = cgetg(RU,t_COL), p2 = (GEN)lambda[i];
1888: sublambda[i] = (long)p1;
1889: for (j=1; j<RU; j++)
1890: {
1891: p1[j] = p2[j+1];
1892: if (trunc_error((GEN)p1[j])) { *ptsublambda = NULL; return gzero; }
1893: }
1894: }
1895: *ptsublambda = sublambda;
1896: gptr[0]=ptsublambda; gptr[1]=&kR;
1897: gerepilemany(av,gptr,2); return kR;
1898: }
1899:
1900: /* Assuming enough relations, c = Rz is close to an even integer, according
1901: * to Dirichlet's formula. Otherwise, close to a multiple.
1902: * Compute a tentative regulator (not a multiple this time) */
1903: static GEN
1904: compute_check(GEN sublambda, GEN z, GEN *parch, GEN *reg)
1905: {
1906: long av = avma, av2, tetpil;
1907: GEN p1,c,den, R = *reg; /* multiple of regulator */
1908:
1909: if (DEBUGLEVEL) { fprintferr("\n#### Computing check\n"); flusherr(); }
1910: c = gmul(R,z);
1911: sublambda = bestappr(sublambda,c); den = denom(sublambda);
1912: if (gcmp(den,c) > 0)
1913: {
1914: if (DEBUGLEVEL) fprintferr("c = %Z\nden = %Z\n",c,den);
1915: avma=av; return NULL;
1916: }
1917:
1918: p1 = gmul(sublambda,den); tetpil=avma;
1919: *parch = lllint(p1);
1920:
1921: av2=avma; p1 = det2(gmul(sublambda,*parch));
1922: affrr(mpabs(gmul(R,p1)), R); avma=av2;
1923:
1924: if (DEBUGLEVEL) msgtimer("bestappr/regulator");
1925: *parch = gerepile(av,tetpil,*parch); return gmul(R,z);
1926: }
1927:
1928: /* find the smallest (wrt norm) among I, I^-1 and red(I^-1) */
1929: static GEN
1930: inverse_if_smaller(GEN nf, GEN I, long prec)
1931: {
1932: GEN J, d, dmin, I1;
1933: int inv;
1934: J = (GEN)I[1];
1935: dmin = dethnf_i(J); I1 = idealinv(nf,I);
1936: J = (GEN)I1[1]; J = gmul(J,denom(J)); I1[1] = (long)J;
1937: d = dethnf_i(J);
1938: if (cmpii(d,dmin) < 0) {inv=1; I=I1; dmin=d;}
1939: else inv=0;
1940: /* try reducing (often _increases_ the norm) */
1941: I1 = ideallllred(nf,I1,NULL,prec);
1942: J = (GEN)I1[1];
1943: d = dethnf_i(J);
1944: if (cmpii(d,dmin) < 0) {inv=1; I=I1;}
1945: return I;
1946: }
1947:
1948: /* in place */
1949: static void
1950: neg_row(GEN U, long i)
1951: {
1952: GEN c = U + lg(U)-1;
1953: for (; c>U; c--) coeff(c,i,0) = lneg(gcoeff(c,i,0));
1954: }
1955:
1956: static void
1957: setlg_col(GEN U, long l)
1958: {
1959: GEN c = U + lg(U)-1;
1960: for (; c>U; c--) setlg(*c, l);
1961: }
1962:
1963: /* compute class group (clg1) + data for isprincipal (clg2) */
1964: static void
1965: class_group_gen(GEN nf,GEN W,GEN C,GEN vperm,GEN *ptclg1,GEN *ptclg2,long prec)
1966: {
1967: GEN z,G,Ga,ga,GD,cyc,X,Y,D,U,V,Ur,Ui,Uir,I,J,Vbase;
1968: long i,j,s,lo,lo0;
1969:
1970: if (DEBUGLEVEL)
1971: { fprintferr("\n#### Computing class group generators\n"); timer2(); }
1972: z = smith2(W); /* U W V = D, D diagonal, G = g Ui (G=new gens, g=old) */
1973: U = (GEN)z[1]; Ui = ginv(U);
1974: V = (GEN)z[2];
1975: D = (GEN)z[3];
1976: lo0 = lo = lg(D);
1977: /* we could set lo = lg(cyc) and truncate all matrices below
1978: * setlg_col(D && U && Y, lo) + setlg(D && V && X && Ui, lo)
1979: * but it's not worth the complication:
1980: * 1) gain is negligible (avoid computing z^0 if lo < lo0)
1981: * 2) when computing ga, the products XU and VY use the original matrices
1982: */
1983: Ur = reducemodHNF(U, D, &Y);
1984: Uir = reducemodHNF(Ui,W, &X);
1985: /* [x] = logarithmic embedding of x (arch. component)
1986: * NB: z = idealred(I) --> I = y z[1], with [y] = - z[2]
1987: * P invertible diagonal matrix (\pm 1) which is only implicitly defined
1988: * G = g Uir P + [Ga], Uir = Ui + WX
1989: * g = G P Ur + [ga], Ur = U + DY */
1990: Vbase = cgetg(lo0,t_VEC);
1991: if (typ(vperm) == t_VECSMALL)
1992: for (i=1; i<lo0; i++) Vbase[i] = vectbase[vperm[i]];
1993: else
1994: for (i=1; i<lo0; i++) Vbase[i] = vectbase[itos((GEN)vperm[i])];
1995:
1996: G = cgetg(lo,t_VEC);
1997: Ga= cgetg(lo,t_VEC);
1998: z = init_idele(nf);
1999: for (j=1; j<lo; j++)
2000: {
2001: GEN p1 = gcoeff(Uir,1,j);
2002: z[1]=Vbase[1]; I = idealpowred(nf,z,p1,prec);
2003: if (signe(p1)<0) I[1] = lmul((GEN)I[1],denom((GEN)I[1]));
2004: for (i=2; i<lo0; i++)
2005: {
2006: p1 = gcoeff(Uir,i,j); s = signe(p1);
2007: if (s)
2008: {
2009: z[1]=Vbase[i]; J = idealpowred(nf,z,p1,prec);
2010: if (s<0) J[1] = lmul((GEN)J[1],denom((GEN)J[1]));
2011: I = idealmulh(nf,I,J);
2012: I = ideallllred(nf,I,NULL,prec);
2013: }
2014: }
2015: J = inverse_if_smaller(nf, I, prec);
2016: if (J != I)
2017: { /* update wrt P */
2018: neg_row(Y ,j); V[j] = lneg((GEN)V[j]);
2019: neg_row(Ur,j); X[j] = lneg((GEN)X[j]);
2020: }
2021: G[j] = (long)J[1]; /* generator, order cyc[j] */
2022: Ga[j]= (long)J[2];
2023: }
2024: /* at this point Y = PY, Ur = PUr, V = VP, X = XP */
2025:
2026: /* G D =: [GD] = g (UiP + W XP) D + [Ga]D = g W (VP + XP D) + [Ga]D
2027: * NB: DP = PD and Ui D = W V. gW is given by (first lo0-1 cols of) C
2028: */
2029: GD = gadd(act_arch(gadd(V, gmul(X,D)), C),
2030: act_arch(D, Ga));
2031: /* -[ga] = [GD]PY + G PU - g = [GD]PY + [Ga] PU + gW XP PU
2032: = gW (XP PUr + VP PY) + [Ga]PUr */
2033: ga = gadd(act_arch(gadd(gmul(X,Ur), gmul(V,Y)), C),
2034: act_arch(Ur, Ga));
2035: ga = gneg(ga);
2036: /* TODO: could (LLL)reduce ga and GD mod units ? */
2037:
2038: cyc = cgetg(lo,t_VEC); /* elementary divisors */
2039: for (j=1; j<lo; j++)
2040: {
2041: cyc[j] = coeff(D,j,j);
2042: if (gcmp1((GEN)cyc[j]))
2043: { /* strip useless components */
2044: lo = j; setlg(cyc,lo); setlg_col(Ur,lo);
2045: setlg(G,lo); setlg(Ga,lo); setlg(GD,lo); break;
2046: }
2047: }
2048: z = cgetg(4,t_VEC); *ptclg1 = z;
2049: z[1]=(long)dethnf_i(W);
2050: z[2]=(long)cyc;
2051: z[3]=(long)G;
2052: z = cgetg(4,t_VEC); *ptclg2 = z;
2053: z[1]=(long)Ur;
2054: z[2]=(long)ga;
2055: z[3]=(long)GD;
2056: if (DEBUGLEVEL) msgtimer("classgroup generators");
2057: }
2058:
2059: /* real(MC * M), where columns a and b of MC have been multiplied by 2^20+1 */
2060: static GEN
2061: shift_t2(GEN T2, GEN M, GEN MC, long a, long b)
2062: {
2063: long i,j,N = lg(T2);
2064: GEN z, t2 = cgetg(N,t_MAT);
2065: for (j=1; j<N; j++)
2066: {
2067: t2[j] = lgetg(N,t_COL);
2068: for (i=1; i<=j; i++)
2069: {
2070: z = mul_real(gcoeff(MC,i,a), gcoeff(M,a,j));
2071: if (a!=b) z = gadd(z, mul_real(gcoeff(MC,i,b), gcoeff(M,b,j)));
2072: coeff(t2,j,i) = coeff(t2,i,j) = ladd(gcoeff(T2,i,j), gmul2n(z,20));
2073: }
2074: }
2075: return t2;
2076: }
2077:
2078: static GEN
2079: compute_vecT2(long RU,GEN nf)
2080: {
2081: GEN vecT2, M = gmael(nf,5,1), MC = gmael(nf,5,2), T2 = gmael(nf,5,3);
2082: long i,j,n = min(RU,9), ind = 1;
2083:
2084: vecT2=cgetg(1 + n*(n+1)/2,t_VEC);
2085: for (j=1; j<=n; j++)
2086: for (i=1; i<=j; i++) vecT2[ind++] = (long)shift_t2(T2,M,MC,i,j);
2087: if (DEBUGLEVEL) msgtimer("weighted T2 matrices");
2088: return vecT2;
2089: }
2090:
2091: /* cf. relationrank()
2092: *
2093: * If V depends linearly from the columns of the matrix, return 0.
2094: * Otherwise, update INVP and L and return 1. No GC.
2095: */
2096: int
2097: addcolumntomatrix(GEN V, GEN invp, GEN L)
2098: {
2099: GEN a = gmul_mat_smallvec(invp,V);
2100: long i,j,k, n = lg(invp);
2101:
2102: if (DEBUGLEVEL>6)
2103: {
2104: fprintferr("adding vector = %Z\n",V);
2105: fprintferr("vector in new basis = %Z\n",a);
2106: fprintferr("list = %Z\n",L);
2107: fprintferr("base change matrix =\n"); outerr(invp);
2108: }
2109: k = 1; while (k<n && (L[k] || gcmp0((GEN)a[k]))) k++;
2110: if (k == n) return 0;
2111: L[k] = 1;
2112: for (i=k+1; i<n; i++) a[i] = ldiv(gneg_i((GEN)a[i]),(GEN)a[k]);
2113: for (j=1; j<=k; j++)
2114: {
2115: GEN c = (GEN)invp[j], ck = (GEN)c[k];
2116: if (gcmp0(ck)) continue;
2117: c[k] = ldiv(ck, (GEN)a[k]);
2118: if (j==k)
2119: for (i=k+1; i<n; i++)
2120: c[i] = lmul((GEN)a[i], ck);
2121: else
2122: for (i=k+1; i<n; i++)
2123: c[i] = ladd((GEN)c[i], gmul((GEN)a[i], ck));
2124: }
2125: return 1;
2126: }
2127:
2128: /* compute the rank of A in M_n,r(Z) (C long), where rank(A) = r <= n.
2129: * Starting from the identity (canonical basis of Q^n), we obtain a base
2130: * change matrix P by taking the independant columns of A and vectors from
2131: * the canonical basis. Update invp = 1/P, and L in {0,1}^n, with L[i] = 0
2132: * of P[i] = e_i, and 1 if P[i] = A_i (i-th column of A)
2133: */
2134: static GEN
2135: relationrank(GEN *A, long r, GEN L)
2136: {
2137: long i, n = lg(L)-1, av = avma, lim = stack_lim(av,1);
2138: GEN invp = idmat(n);
2139:
2140: if (!r) return invp;
2141: if (r>n) err(talker,"incorrect matrix in relationrank");
2142: for (i=1; i<=r; i++)
2143: {
2144: if (! addcolumntomatrix(A[i],invp,L) && i == r)
2145: err(talker,"not a maximum rank matrix in relationrank");
2146: if (low_stack(lim, stack_lim(av,1)))
2147: {
2148: if(DEBUGMEM>1) err(warnmem,"relationrank");
2149: invp = gerepilecopy(av, invp);
2150: }
2151: }
2152: return gerepilecopy(av, invp);
2153: }
2154:
2155: static GEN
2156: codeprime(GEN bnf, GEN pr)
2157: {
2158: long j,av=avma,tetpil;
2159: GEN p,al,fa,p1;
2160:
2161: p=(GEN)pr[1]; al=(GEN)pr[2]; fa=primedec(bnf,p);
2162: for (j=1; j<lg(fa); j++)
2163: if (gegal(al,gmael(fa,j,2)))
2164: {
2165: p1=mulsi(lg(al)-1,p); tetpil=avma;
2166: return gerepile(av,tetpil,addsi(j-1,p1));
2167: }
2168: err(talker,"bug in codeprime/smallbuchinit");
2169: return NULL; /* not reached */
2170: }
2171:
2172: static GEN
2173: decodeprime(GEN nf, GEN co)
2174: {
2175: long n,indi,av=avma;
2176: GEN p,rem,p1;
2177:
2178: n=lg(nf[7])-1; p=dvmdis(co,n,&rem); indi=itos(rem)+1;
2179: p1=primedec(nf,p);
2180: return gerepilecopy(av,(GEN)p1[indi]);
2181: }
2182:
2183: /* v = bnf[10] */
2184: GEN
2185: get_matal(GEN v)
2186: {
2187: if (typ(v) == t_VEC)
2188: {
2189: GEN ma = (GEN)v[1];
2190: if (typ(ma) != t_INT) return ma;
2191: }
2192: return NULL;
2193: }
2194:
2195: GEN
2196: get_cycgen(GEN v)
2197: {
2198: if (typ(v) == t_VEC)
2199: {
2200: GEN h = (GEN)v[2];
2201: if (typ(h) == t_VEC) return h;
2202: }
2203: return NULL;
2204: }
2205:
2206: /* compute principal ideals corresponding to (gen[i]^cyc[i]) */
2207: static GEN
2208: makecycgen(GEN bnf)
2209: {
2210: GEN cyc,gen,h,nf,y,GD,D;
2211: long e,i,l;
2212:
2213: h = get_cycgen((GEN)bnf[10]);
2214: if (h) return h;
2215:
2216: nf = checknf(bnf);
2217: cyc = gmael3(bnf,8,1,2); D = diagonal(cyc);
2218: gen = gmael3(bnf,8,1,3); GD = gmael(bnf,9,3);
2219: l = lg(gen); h = cgetg(l, t_VEC);
2220: for (i=1; i<l; i++)
2221: {
2222: if (cmpis((GEN)cyc[i], 16) < 0)
2223: {
2224: GEN N = dethnf_i((GEN)gen[i]);
2225: y = isprincipalarch(bnf,(GEN)GD[i], N, (GEN)cyc[i], gun, &e);
2226: if (y && !fact_ok(nf,y,NULL,gen,(GEN)D[i])) y = NULL;
2227: if (y) { h[i] = (long)to_famat_all(y,gun); continue; }
2228: }
2229: y = isprincipalfact(bnf, gen, (GEN)D[i], NULL,
2230: nf_GEN|nf_GENMAT|nf_FORCE|nf_GIVEPREC);
2231: if (typ(y) != t_INT)
2232: h[i] = y[2];
2233: else
2234: {
2235: y = idealpow(nf, (GEN)gen[i], (GEN)cyc[i]);
2236: h[i] = isprincipalgenforce(bnf,y)[2];
2237: }
2238: }
2239: return h;
2240: }
2241:
2242: /* compute principal ideals corresponding to bnf relations */
2243: static GEN
2244: makematal(GEN bnf)
2245: {
2246: GEN W,B,pFB,vp,nf,ma, WB_C;
2247: long lW,lma,j,prec;
2248:
2249: ma = get_matal((GEN)bnf[10]);
2250: if (ma) return ma;
2251:
2252: W = (GEN)bnf[1];
2253: B = (GEN)bnf[2];
2254: WB_C= (GEN)bnf[4];
2255: vp = (GEN)bnf[6];
2256: nf = (GEN)bnf[7];
2257: lW=lg(W)-1; lma=lW+lg(B);
2258: pFB=cgetg(lma,t_VEC); ma=(GEN)bnf[5]; /* reorder factor base */
2259: for (j=1; j<lma; j++) pFB[j] = ma[itos((GEN)vp[j])];
2260: ma = cgetg(lma,t_MAT);
2261:
2262: prec = prec_arch(bnf);
2263: for (j=1; j<lma; j++)
2264: {
2265: long c = getrand(), e;
2266: GEN ex = (j<=lW)? (GEN)W[j]: (GEN)B[j-lW];
2267: GEN C = (j<=lW)? NULL: (GEN)pFB[j];
2268: GEN dx, Nx = get_norm_fact_primes(pFB, ex, C, &dx);
2269: GEN y = isprincipalarch(bnf,(GEN)WB_C[j], Nx,gun, dx, &e);
2270: if (y && !fact_ok(nf,y,C,pFB,ex)) y = NULL;
2271: if (y)
2272: {
2273: if (DEBUGLEVEL>1) fprintferr("*%ld ",j);
2274: ma[j] = (long)y; continue;
2275: }
2276:
2277: if (!y) y = isprincipalfact(bnf,pFB,ex,C, nf_GEN|nf_FORCE|nf_GIVEPREC);
2278: if (typ(y) != t_INT)
2279: {
2280: if (DEBUGLEVEL>1) fprintferr("%ld ",j);
2281: ma[j] = y[2]; continue;
2282: }
2283:
2284: prec = itos(y); j--;
2285: if (DEBUGLEVEL) err(warnprec,"makematal",prec);
2286: nf = nfnewprec(nf,prec);
2287: bnf = bnfinit0(nf,1,NULL,prec); setrand(c);
2288: }
2289: if (DEBUGLEVEL>1) fprintferr("\n");
2290: return ma;
2291: }
2292:
2293: /* insert O in bnf at index K
2294: * K = 1: matal
2295: * K = 2: cycgen */
2296: static void
2297: bnfinsert(GEN bnf, GEN O, long K)
2298: {
2299: GEN v = (GEN)bnf[10];
2300: if (typ(v) != t_VEC)
2301: {
2302: GEN w = cgetg(3, t_VEC);
2303: long i;
2304: for (i = 1; i < 3; i++)
2305: w[i] = (i==K)? (long)O: zero;
2306: w = gclone(w);
2307: bnf[10] = (long)w;
2308: }
2309: else
2310: v[K] = lclone(O);
2311: }
2312:
2313: GEN
2314: check_and_build_cycgen(GEN bnf)
2315: {
2316: GEN cycgen = get_cycgen((GEN)bnf[10]);
2317: if (!cycgen)
2318: {
2319: long av = avma;
2320: if (DEBUGLEVEL) err(warner,"completing bnf (building cycgen)");
2321: bnfinsert(bnf, makecycgen(bnf), 2); avma = av;
2322: cycgen = get_cycgen((GEN)bnf[10]);
2323: }
2324: return cycgen;
2325: }
2326:
2327: GEN
2328: check_and_build_matal(GEN bnf)
2329: {
2330: GEN matal = get_matal((GEN)bnf[10]);
2331: if (!matal)
2332: {
2333: long av = avma;
2334: if (DEBUGLEVEL) err(warner,"completing bnf (building matal)");
2335: bnfinsert(bnf, makematal(bnf), 1); avma = av;
2336: matal = get_matal((GEN)bnf[10]);
2337: }
2338: return matal;
2339: }
2340:
2341: GEN
2342: smallbuchinit(GEN pol,GEN gcbach,GEN gcbach2,GEN gRELSUP,GEN gborne,long nbrelpid,long minsFB,long prec)
2343: {
2344: long av=avma,av1,k;
2345: GEN y,bnf,pFB,vp,nf,mas,res,uni,v1,v2,v3;
2346:
2347: if (typ(pol)==t_VEC) bnf = checkbnf(pol);
2348: else
2349: {
2350: bnf=buchall(pol,gcbach,gcbach2,gRELSUP,gborne,nbrelpid,minsFB,-3,prec);
2351: bnf = checkbnf_discard(bnf);
2352: }
2353: pFB=(GEN)bnf[5]; vp=(GEN)bnf[6]; nf=(GEN)bnf[7];
2354: mas=(GEN)nf[5]; res=(GEN)bnf[8]; uni=(GEN)res[5];
2355:
2356: y=cgetg(13,t_VEC); y[1]=lcopy((GEN)nf[1]); y[2]=lcopy(gmael(nf,2,1));
2357: y[3]=lcopy((GEN)nf[3]); y[4]=lcopy((GEN)nf[7]);
2358: y[5]=lcopy((GEN)nf[6]); y[6]=lcopy((GEN)mas[5]);
2359: y[7]=lcopy((GEN)bnf[1]); y[8]=lcopy((GEN)bnf[2]);
2360: v1=cgetg(lg(pFB),t_VEC); y[9]=(long)v1;
2361: for (k=1; k<lg(pFB); k++)
2362: v1[k]=(long)codeprime(bnf,(GEN)pFB[itos((GEN)vp[k])]);
2363: v2=cgetg(3,t_VEC); y[10]=(long)v2;
2364: v2[1]=lcopy(gmael(res,4,1));
2365: v2[2]=(long)algtobasis(bnf,gmael(res,4,2));
2366: v3=cgetg(lg(uni),t_VEC); y[11]=(long)v3;
2367: for (k=1; k<lg(uni); k++)
2368: v3[k] = (long)algtobasis(bnf,(GEN)uni[k]);
2369: av1 = avma;
2370: y[12] = gcmp0((GEN)bnf[10])? lpilecopy(av1, makematal(bnf))
2371: : lcopy((GEN)bnf[10]);
2372: return gerepileupto(av,y);
2373: }
2374:
2375: static GEN
2376: get_regulator(GEN mun,long prec)
2377: {
2378: long av,tetpil;
2379: GEN p1;
2380:
2381: if (lg(mun)==1) return gun;
2382: av=avma; p1 = gtrans(greal(mun));
2383: setlg(p1,lg(p1)-1); p1 = det(p1);
2384: tetpil=avma; return gerepile(av,tetpil,gabs(p1,prec));
2385: }
2386:
2387: static GEN
2388: log_poleval(GEN P, GEN *ro, long i, GEN nf, long prec0)
2389: {
2390: GEN x = poleval(P, (GEN)(*ro)[i]);
2391: long prec = prec0, k = 0;
2392: while (gcmp0(x) || ((k = gprecision(x)) && k < DEFAULTPREC))
2393: {
2394: prec = (prec-2)<<1;
2395: if (DEBUGLEVEL) err(warnprec,"log_poleval",prec);
2396: *ro = get_roots((GEN)nf[1], itos(gmael(nf,2,1)),lg(*ro)-1,prec);
2397: x = poleval(P, (GEN)(*ro)[i]);
2398: }
2399: if (k > prec0) x = gprec_w(x,prec0);
2400: return glog(x, prec0);
2401: }
2402:
2403: /* return corrected archimedian components
2404: * (= log(sigma_i(x)) - log(|Nx|)/[K:Q]) for a (matrix of elements) */
2405: static GEN
2406: get_arch2_i(GEN nf, GEN a, long prec, int units)
2407: {
2408: GEN M,N, ro = dummycopy((GEN)nf[6]);
2409: long j,k, la = lg(a), ru = lg(ro), r1 = nf_get_r1(nf);
2410:
2411: M = cgetg(la,t_MAT); if (la == 1) return M;
2412: if (typ(a[1]) == t_COL) a = gmul((GEN)nf[7], a);
2413: if (units) N = NULL; /* no correction necessary */
2414: else
2415: {
2416: GEN pol = (GEN)nf[1];
2417: N = cgetg(la,t_VEC);
2418: for (k=1; k<la; k++) N[k] = (long)gabs(subres(pol,(GEN)a[k]),0);
2419: N = gdivgs(glog(N,prec), - (degpol(pol))); /* - log(|norm|) / [K:Q] */
2420: }
2421: for (k=1; k<la; k++)
2422: {
2423: GEN p, c = cgetg(ru,t_COL); M[k] = (long)c;
2424: for (j=1; j<ru; j++)
2425: {
2426: p = log_poleval((GEN)a[k],&ro,j,nf,prec);
2427: if (N) p = gadd(p, (GEN)N[k]);
2428: c[j]=(j<=r1)? (long) p: lmul2n(p,1);
2429: }
2430: }
2431: return M;
2432: }
2433:
2434: static GEN
2435: get_arch2(GEN nf, GEN a, long prec, int units)
2436: {
2437: long av = avma;
2438: return gerepilecopy(av, get_arch2_i(nf,a,prec,units));
2439: }
2440:
2441: static void
2442: my_class_group_gen(GEN bnf, GEN *ptcl, GEN *ptcl2, long prec)
2443: {
2444: GEN W=(GEN)bnf[1], C=(GEN)bnf[4], vperm=(GEN)bnf[6], nf=(GEN)bnf[7], *gptr[2];
2445: long av = avma;
2446:
2447: vectbase = (GEN)bnf[5]; /* needed by class_group_gen */
2448: class_group_gen(nf,W,C,vperm,ptcl,ptcl2, prec);
2449: gptr[0]=ptcl; gptr[1]=ptcl2; gerepilemany(av,gptr,2);
2450: }
2451:
2452: GEN
2453: bnfnewprec(GEN bnf, long prec)
2454: {
2455: GEN nf,ro,res,p1,funits,mun,matal,clgp,clgp2,y;
2456: long r1,r2,ru,pl1,pl2,prec1;
2457:
2458: bnf = checkbnf(bnf);
2459: if (prec <= 0) return nfnewprec(checknf(bnf),prec);
2460: y = cgetg(11,t_VEC);
2461: funits = check_units(bnf,"bnfnewprec");
2462: nf = (GEN)bnf[7];
2463: ro = (GEN)nf[6]; ru = lg(ro)-1;
2464: r1 = nf_get_r1(nf);
2465: r2 = (r1 + ru) >> 1;
2466: pl1 = (ru == 1 && r1 == 0)? 0: gexpo(funits);
2467: pl2 = gexpo(ro);
2468: prec1 = prec;
2469: prec += ((ru + r2 - 1) * (pl1 + (ru + r2) * pl2)) >> TWOPOTBITS_IN_LONG;
2470: nf = nfnewprec((GEN)bnf[7],prec);
2471: res = cgetg(7,t_VEC);
2472: ro = (GEN)nf[6];
2473: mun = get_arch2(nf,funits,prec,1);
2474: if (prec != prec1) { mun = gprec_w(mun,prec1); prec = prec1; }
2475: res[2]=(long)get_regulator(mun,prec);
2476: p1 = (GEN)bnf[8];
2477: res[3]=lcopy((GEN)p1[3]);
2478: res[4]=lcopy((GEN)p1[4]);
2479: res[5]=lcopy((GEN)p1[5]);
2480: res[6]=lcopy((GEN)p1[6]);
2481:
2482: y[1]=lcopy((GEN)bnf[1]);
2483: y[2]=lcopy((GEN)bnf[2]);
2484: y[3]=(long)mun;
2485: matal = check_and_build_matal(bnf);
2486: y[4]=(long)get_arch2(nf,matal,prec,0);
2487: y[5]=lcopy((GEN)bnf[5]);
2488: y[6]=lcopy((GEN)bnf[6]);
2489: y[7]=(long)nf;
2490: y[8]=(long)res;
2491: my_class_group_gen(y,&clgp,&clgp2,prec);
2492: res[1]=(long)clgp;
2493: y[9]=(long)clgp2;
2494: y[10]=lcopy((GEN)bnf[10]); return y;
2495: }
2496:
2497: GEN
2498: bnrnewprec(GEN bnr, long prec)
2499: {
2500: GEN y = cgetg(7,t_VEC);
2501: long i;
2502: checkbnr(bnr);
2503: y[1] = (long)bnfnewprec((GEN)bnr[1],prec);
2504: for (i=2; i<7; i++) y[i]=lcopy((GEN)bnr[i]);
2505: return y;
2506: }
2507:
2508: GEN
2509: bnfmake(GEN sbnf, long prec)
2510: {
2511: long av = avma, j,k,n,r1,r2,ru,lpf;
2512: GEN p1,x,bas,ro,nf,mun,funits,index;
2513: GEN pfc,vp,mc,clgp,clgp2,res,y,W,racu,reg,matal;
2514:
2515: if (typ(sbnf)!=t_VEC || lg(sbnf)!=13)
2516: err(talker,"incorrect sbnf in bnfmake");
2517: x=(GEN)sbnf[1]; bas=(GEN)sbnf[4]; n=lg(bas)-1;
2518: r1=itos((GEN)sbnf[2]); r2=(n-r1)>>1; ru=r1+r2;
2519: ro=(GEN)sbnf[5];
2520: if (prec > gprecision(ro)) ro=get_roots(x,r1,ru,prec);
2521: index = gun;
2522: for (j=2; j<=n; j++) index = mulii(index, denom(leading_term(bas[j])));
2523:
2524: nf=cgetg(10,t_VEC);
2525: nf[1]=sbnf[1]; p1=cgetg(3,t_VEC); p1[1]=lstoi(r1); p1[2]=lstoi(r2);
2526: nf[2]=(long)p1;
2527: nf[3]=sbnf[3];
2528: nf[4]=(long)index;
2529: nf[6]=(long)ro;
2530: nf[7]=(long)bas;
2531: get_nf_matrices(nf, 0);
2532:
2533: funits=cgetg(ru,t_VEC); p1 = (GEN)sbnf[11];
2534: for (k=1; k < lg(p1); k++)
2535: funits[k] = lmul(bas,(GEN)p1[k]);
2536: mun = get_arch2_i(nf,funits,prec,1);
2537:
2538: prec=gprecision(ro); if (prec<DEFAULTPREC) prec=DEFAULTPREC;
2539: matal = get_matal((GEN)sbnf[12]);
2540: if (!matal) matal = (GEN)sbnf[12];
2541: mc = get_arch2_i(nf,matal,prec,0);
2542:
2543: pfc=(GEN)sbnf[9]; lpf=lg(pfc);
2544: vectbase=cgetg(lpf,t_COL); vp=cgetg(lpf,t_COL);
2545: for (j=1; j<lpf; j++)
2546: {
2547: vp[j]=lstoi(j);
2548: vectbase[j]=(long)decodeprime(nf,(GEN)pfc[j]);
2549: }
2550: W = (GEN)sbnf[7];
2551: class_group_gen(nf,W,mc,vp,&clgp,&clgp2, prec); /* uses vectbase */
2552:
2553: reg = get_regulator(mun,prec);
2554: p1=cgetg(3,t_VEC); racu=(GEN)sbnf[10];
2555: p1[1]=racu[1]; p1[2]=lmul(bas,(GEN)racu[2]);
2556: racu=p1;
2557:
2558: res=cgetg(7,t_VEC);
2559: res[1]=(long)clgp; res[2]=(long)reg; res[3]=(long)dbltor(1.0);
2560: res[4]=(long)racu; res[5]=(long)funits; res[6]=lstoi(1000);
2561:
2562: y=cgetg(11,t_VEC);
2563: y[1]=(long)W; y[2]=sbnf[8]; y[3]=(long)mun;
2564: y[4]=(long)mc; y[5]=(long)vectbase; y[6]=(long)vp;
2565: y[7]=(long)nf; y[8]=(long)res; y[9]=(long)clgp2; y[10]=sbnf[12];
2566: return gerepilecopy(av,y);
2567: }
2568:
2569: static GEN
2570: classgroupall(GEN P, GEN data, long flag, long prec)
2571: {
2572: long court[3],doubl[4];
2573: long av=avma,flun,lx, minsFB=3,nbrelpid=4;
2574: GEN bach=doubl,bach2=doubl,RELSUP=court,borne=gun;
2575:
2576: if (!data) lx=1;
2577: else
2578: {
2579: lx = lg(data);
2580: if (typ(data)!=t_VEC || lx > 7)
2581: err(talker,"incorrect parameters in classgroup");
2582: }
2583: court[0]=evaltyp(t_INT) | evallg(3); affsi(5,court);
2584: doubl[0]=evaltyp(t_REAL)| evallg(4); affrr(dbltor(0.3),doubl);
2585: avma=av;
2586: switch(lx)
2587: {
2588: case 7: minsFB = itos((GEN)data[6]);
2589: case 6: nbrelpid= itos((GEN)data[5]);
2590: case 5: borne = (GEN)data[4];
2591: case 4: RELSUP = (GEN)data[3];
2592: case 3: bach2 = (GEN)data[2];
2593: case 2: bach = (GEN)data[1];
2594: }
2595: switch(flag)
2596: {
2597: case 0: flun=-2; break;
2598: case 1: flun=-3; break;
2599: case 2: flun=-1; break;
2600: case 3: return smallbuchinit(P,bach,bach2,RELSUP,borne,nbrelpid,minsFB,prec);
2601: case 4: flun=2; break;
2602: case 5: flun=3; break;
2603: case 6: flun=0; break;
2604: default: err(flagerr,"classgroupall");
2605: return NULL; /* not reached */
2606: }
2607: return buchall(P,bach,bach2,RELSUP,borne,nbrelpid,minsFB,flun,prec);
2608: }
2609:
2610: GEN
2611: bnfclassunit0(GEN P, long flag, GEN data, long prec)
2612: {
2613: if (typ(P)==t_INT) return quadclassunit0(P,0,data,prec);
2614: if (flag < 0 || flag > 2) err(flagerr,"bnfclassunit");
2615: return classgroupall(P,data,flag+4,prec);
2616: }
2617:
2618: GEN
2619: bnfinit0(GEN P, long flag, GEN data, long prec)
2620: {
2621: #if 0
2622: /* TODO: synchronize quadclassunit output with bnfinit */
2623: if (typ(P)==t_INT)
2624: {
2625: if (flag<4) err(impl,"specific bnfinit for quadratic fields");
2626: return quadclassunit0(P,0,data,prec);
2627: }
2628: #endif
2629: if (flag < 0 || flag > 3) err(flagerr,"bnfinit");
2630: return classgroupall(P,data,flag,prec);
2631: }
2632:
2633: GEN
2634: classgrouponly(GEN P, GEN data, long prec)
2635: {
2636: ulong av = avma;
2637: GEN z;
2638:
2639: if (typ(P)==t_INT)
2640: {
2641: z=quadclassunit0(P,0,data,prec); setlg(z,4);
2642: return gerepilecopy(av,z);
2643: }
2644: z=(GEN)classgroupall(P,data,6,prec)[1];
2645: return gerepilecopy(av,(GEN)z[5]);
2646: }
2647:
2648: GEN
2649: regulator(GEN P, GEN data, long prec)
2650: {
2651: ulong av = avma;
2652: GEN z;
2653:
2654: if (typ(P)==t_INT)
2655: {
2656: if (signe(P)>0)
2657: {
2658: z=quadclassunit0(P,0,data,prec);
2659: return gerepilecopy(av,(GEN)z[4]);
2660: }
2661: return gun;
2662: }
2663: z=(GEN)classgroupall(P,data,6,prec)[1];
2664: return gerepilecopy(av,(GEN)z[6]);
2665: }
2666:
2667: #ifdef INLINE
2668: INLINE
2669: #endif
2670: GEN
2671: col_0(long n)
2672: {
2673: GEN c = (GEN) gpmalloc(sizeof(long)*(n+1));
2674: long i;
2675: for (i=1; i<=n; i++) c[i]=0;
2676: c[0] = evaltyp(t_VECSMALL) | evallg(n);
2677: return c;
2678: }
2679:
2680: GEN
2681: cgetc_col(long n, long prec)
2682: {
2683: GEN c = cgetg(n+1,t_COL);
2684: long i;
2685: for (i=1; i<=n; i++) c[i] = (long)cgetc(prec);
2686: return c;
2687: }
2688:
2689: static GEN
2690: buchall_end(GEN nf,GEN CHANGE,long fl,long k, GEN fu, GEN clg1, GEN clg2,
2691: GEN reg, GEN c_1, GEN zu, GEN W, GEN B,
2692: GEN xarch, GEN matarch, GEN vectbase, GEN vperm)
2693: {
2694: long i, l = labs(fl)>1? 11: fl? 9: 8;
2695: GEN p1,z, RES = cgetg(11,t_COL);
2696:
2697: setlg(RES,l);
2698: RES[5]=(long)clg1;
2699: RES[6]=(long)reg;
2700: RES[7]=(long)c_1;
2701: RES[8]=(long)zu;
2702: RES[9]=(long)fu;
2703: RES[10]=lstoi(k);
2704: if (fl>=0)
2705: {
2706: RES[1]=nf[1];
2707: RES[2]=nf[2]; p1=cgetg(3,t_VEC); p1[1]=nf[3]; p1[2]=nf[4];
2708: RES[3]=(long)p1;
2709: RES[4]=nf[7];
2710: z=cgetg(2,t_MAT); z[1]=lcopy(RES); return z;
2711: }
2712: z=cgetg(11,t_VEC);
2713: z[1]=(long)W;
2714: z[2]=(long)B;
2715: z[3]=(long)xarch;
2716: z[4]=(long)matarch;
2717: z[5]=(long)vectbase;
2718: for (i=lg(vperm)-1; i>0; i--) vperm[i]=lstoi(vperm[i]);
2719: settyp(vperm, t_VEC);
2720: z[6]=(long)vperm;
2721: z[7]=(long)nf; RES+=4; RES[0]=evaltyp(t_VEC) | evallg(l-4);
2722: z[8]=(long)RES;
2723: z[9]=(long)clg2;
2724: z[10]=zero; /* dummy: we MUST have lg(bnf) != lg(nf) */
2725: if (CHANGE) { p1=cgetg(3,t_VEC); p1[1]=(long)z; p1[2]=(long)CHANGE; z=p1; }
2726: return gcopy(z);
2727: }
2728:
2729: static GEN
2730: buchall_for_degree_one_pol(GEN nf, GEN CHANGE, long flun)
2731: {
2732: long av = avma, k = EXP220;
2733: GEN W,B,xarch,matarch,vectbase,vperm;
2734: GEN fu=cgetg(1,t_VEC), reg=gun, c_1=gun, zu=cgetg(3,t_VEC);
2735: GEN clg1=cgetg(4,t_VEC), clg2=cgetg(4,t_VEC);
2736:
2737: clg1[1]=un; clg1[2]=clg1[3]=clg2[3]=lgetg(1,t_VEC);
2738: clg2[1]=clg2[2]=lgetg(1,t_MAT);
2739: zu[1]=deux; zu[2]=lnegi(gun);
2740: W=B=xarch=matarch=cgetg(1,t_MAT);
2741: vectbase=cgetg(1,t_COL); vperm=cgetg(1,t_VEC);
2742:
2743: return gerepileupto(av, buchall_end(nf,CHANGE,flun,k,fu,clg1,clg2,reg,c_1,zu,W,B,xarch,matarch,vectbase,vperm));
2744: }
2745:
2746: static GEN
2747: get_LLLnf(GEN nf, long prec)
2748: {
2749: GEN M = gmael(nf,5,1);
2750: GEN T2 = gmael(nf,5,3);
2751: GEN cbase = lllgramintern(T2,100,1,prec);
2752: GEN v = cgetg(5,t_VEC);
2753: if (!cbase) return NULL;
2754: if (gegal(cbase, idmat(lg(T2)-1))) cbase = NULL;
2755: v[1] = (long) (cbase? qf_base_change(T2, cbase, 1): T2);
2756: v[2] = (long) (cbase? gmul(M, cbase): M);
2757: v[3] = (long) cbase;
2758: v[4] = (long) (cbase? ZM_inv(cbase,gun): NULL); return v;
2759: }
2760:
2761: GEN
2762: buchall(GEN P,GEN gcbach,GEN gcbach2,GEN gRELSUP,GEN gborne,long nbrelpid,
2763: long minsFB,long flun,long prec)
2764: {
2765: ulong av = avma,av0,av1,limpile;
2766: long N,R1,R2,RU,PRECREG,PRECLLL,KCCO,RELSUP,LIMC,LIMC2,lim;
2767: long nlze,sreg,nrelsup,nreldep,phase,matmax,i,j,k,ss,cglob;
2768: long first=1, sfb_increase=0, sfb_trials=0, precdouble=0, precadd=0;
2769: double cbach,cbach2,drc,LOGD2;
2770: GEN p1,vecT2,fu,zu,nf,LLLnf,D,xarch,W,reg,lfun,z,clh,vperm,subFB;
2771: GEN B,C,c1,sublambda,pdep,parch,liste,invp,clg1,clg2, *mat;
2772: GEN CHANGE=NULL, extramat=NULL, extraC=NULL, list_jideal=NULL;
2773: char *precpb = NULL;
2774:
2775: if (DEBUGLEVEL) timer2();
2776:
2777: if (typ(P)==t_POL) nf = NULL; else { nf = checknf(P); P = (GEN)nf[1]; }
2778: if (typ(gRELSUP) != t_INT) gRELSUP = gtrunc(gRELSUP);
2779: RELSUP = itos(gRELSUP);
2780: if (RELSUP<=0) err(talker,"not enough relations in bnfxxx");
2781:
2782: /* Initializations */
2783: fu = NULL; /* gcc -Wall */
2784: N = degpol(P); PRECREG = max(BIGDEFAULTPREC,prec);
2785: if (!nf)
2786: {
2787: nf = initalgall0(P, nf_REGULAR, PRECREG);
2788: /* P was non-monic and nfinit changed it ? */
2789: if (lg(nf)==3) { CHANGE = (GEN)nf[2]; nf = (GEN)nf[1]; }
2790: }
2791: if (N<=1) return buchall_for_degree_one_pol(nf,CHANGE,flun);
2792: zu = rootsof1(nf);
2793: zu[2] = lmul((GEN)nf[7],(GEN)zu[2]);
2794: if (DEBUGLEVEL) msgtimer("initalg & rootsof1");
2795:
2796: R1 = itos(gmael(nf,2,1)); R2 = (N-R1)>>1; RU = R1+R2;
2797: D = (GEN)nf[3]; drc = fabs(gtodouble(D));
2798: LOGD2 = log(drc); LOGD2 = LOGD2*LOGD2;
2799: lim = (long) (exp(-(double)N) * sqrt(2*PI*N*drc) * pow(4/PI,(double)R2));
2800: if (lim < 3) lim = 3;
2801: cbach = min(12., gtodouble(gcbach)); cbach /= 2;
2802: cbach2 = gtodouble(gcbach2);
2803: if (DEBUGLEVEL)
2804: {
2805: fprintferr("N = %ld, R1 = %ld, R2 = %ld, RU = %ld\n",N,R1,R2,RU);
2806: fprintferr("D = %Z\n",D);
2807: }
2808: av0 = avma; mat = NULL; FB = NULL;
2809:
2810: START:
2811: avma = av0;
2812: if (first) first = 0; else desallocate(&mat);
2813: if (precpb)
2814: {
2815: precdouble++;
2816: if (precadd)
2817: PRECREG += precadd;
2818: else
2819: PRECREG = (PRECREG<<1)-2;
2820: if (DEBUGLEVEL)
2821: {
2822: char str[64]; sprintf(str,"buchall (%s)",precpb);
2823: err(warnprec,str,PRECREG);
2824: }
2825: precpb = NULL;
2826: nf = nfnewprec(nf,PRECREG); av0 = avma;
2827: }
2828: else
2829: cbach = check_bach(cbach,12.);
2830: precadd = 0;
2831:
2832: /* T2-LLL-reduce nf.zk */
2833: LLLnf = get_LLLnf(nf, PRECREG);
2834: if (!LLLnf) { precpb = "LLLnf"; goto START; }
2835:
2836: LIMC = (long)(cbach*LOGD2); if (LIMC < 20) LIMC = 20;
2837: LIMC2= max(3 * N, (long)(cbach2*LOGD2));
2838: if (LIMC2 < LIMC) LIMC2 = LIMC;
2839: if (DEBUGLEVEL) { fprintferr("LIMC = %ld, LIMC2 = %ld\n",LIMC,LIMC2); }
2840:
2841: /* initialize FB, [sub]vperm */
2842: lfun = FBgen(nf,LIMC2,LIMC);
2843: if (!lfun) goto START;
2844: vperm = cgetg(lg(vectbase), t_VECSMALL);
2845: sfb_trials = sfb_increase = 0;
2846: subFB = subFBgen(N,min(lim,LIMC2), minsFB, vperm, &ss);
2847: if (!subFB) goto START;
2848: nreldep = nrelsup = 0;
2849:
2850: /* PRECLLL = prec for LLL-reductions (idealred)
2851: * PRECREG = prec for archimedean components */
2852: PRECLLL = DEFAULTPREC
2853: + ((expi(D)*(lg(subFB)-2)+((N*N)>>2))>>TWOPOTBITS_IN_LONG);
2854: if (!precdouble) PRECREG = prec+1;
2855: if (PRECREG < PRECLLL) PRECREG = PRECLLL;
2856: KCCO = KC+RU-1 + max(ss,RELSUP); /* expected number of needed relations */
2857: if (DEBUGLEVEL)
2858: fprintferr("relsup = %ld, ss = %ld, KCZ = %ld, KC = %ld, KCCO = %ld\n",
2859: RELSUP,ss,KCZ,KC,KCCO);
2860: matmax = 10*KCCO + MAXRELSUP; /* make room for lots of relations */
2861: mat = (GEN*)gpmalloc(sizeof(GEN)*(matmax + 1));
2862: setlg(mat, KCCO+1);
2863: C = cgetg(KCCO+1,t_MAT);
2864: /* trivial relations (p) = prod P^e */
2865: cglob = 0; z = zerocol(RU);
2866: for (i=1; i<=KCZ; i++)
2867: {
2868: GEN P = idealbase[i];
2869: if (isclone(P))
2870: { /* all prime divisors in FB */
2871: unsetisclone(P); cglob++;
2872: C[cglob] = (long)z; /* 0 */
2873: mat[cglob] = p1 = col_0(KC);
2874: k = numideal[FB[i]];
2875: p1[0] = k+1; /* for already_found_relation */
2876: p1 += k;
2877: for (j=lg(P)-1; j; j--) p1[j] = itos(gmael(P,j,3));
2878: }
2879: }
2880: if (DEBUGLEVEL) fprintferr("After trivial relations, cglob = %ld\n",cglob);
2881: /* initialize for other relations */
2882: for (i=cglob+1; i<=KCCO; i++)
2883: {
2884: mat[i] = col_0(KC);
2885: C[i] = (long)cgetc_col(RU,PRECREG);
2886: }
2887: av1 = avma; liste = cgetg(KC+1,t_VECSMALL);
2888: for (i=1; i<=KC; i++) liste[i] = 0;
2889: invp = relationrank(mat,cglob,liste);
2890:
2891: /* relations through elements of small norm */
2892: cglob = small_norm_for_buchall(cglob,mat,C,(long)LIMC,PRECREG,
2893: nf,gborne,nbrelpid,invp,liste,LLLnf);
2894: if (cglob < 0) { precpb = "small_norm"; goto START; }
2895: avma = av1; limpile = stack_lim(av1,1);
2896:
2897: phase = 0;
2898: nlze = matmax = 0; /* for lint */
2899: vecT2 = NULL;
2900:
2901: /* random relations */
2902: if (cglob == KCCO) /* enough rels, but init random_relations just in case */
2903: ((void(*)(long))random_relation)(-1);
2904: else
2905: {
2906: GEN matarch;
2907:
2908: if (DEBUGLEVEL) fprintferr("\n#### Looking for random relations\n");
2909: MORE:
2910: if (sfb_increase)
2911: {
2912: if (DEBUGLEVEL) fprintferr("*** Increasing sub factor base\n");
2913: sfb_increase = 0;
2914: if (++sfb_trials >= SFB_MAX) goto START;
2915: subFB = subFBgen(N,min(lim,LIMC2), lg(subFB)-1+SFB_STEP,NULL,&ss);
2916: if (!subFB) goto START;
2917: nreldep = nrelsup = 0;
2918: }
2919:
2920: if (phase == 0) matarch = C;
2921: else
2922: { /* reduced the relation matrix at least once */
2923: long lgex = max(nlze, MIN_EXTRA); /* number of new relations sought */
2924: long slim; /* total number of relations (including lgex new ones) */
2925: setlg(extraC, lgex+1);
2926: setlg(extramat,lgex+1); /* were allocated after hnfspec */
2927: slim = cglob+lgex;
2928: if (slim > matmax)
2929: {
2930: mat = (long**)gprealloc(mat, (2*slim+1) * sizeof(long*),
2931: (matmax+1) * sizeof(long*));
2932: matmax = 2 * slim;
2933: }
2934: setlg(mat, slim+1);
2935: if (DEBUGLEVEL)
2936: fprintferr("\n(need %ld more relation%s)\n", lgex, (lgex==1)?"":"s");
2937: for (j=cglob+1; j<=slim; j++) mat[j] = col_0(KC);
2938: matarch = extraC - cglob; /* start at 0, the others at cglob */
2939: }
2940: if (!vecT2)
2941: {
2942: vecT2 = compute_vecT2(RU,nf);
2943: av1 = avma;
2944: }
2945: if (!powsubFB)
2946: {
2947: powsubFBgen(nf,subFB,CBUCHG+1,PRECREG);
2948: av1 = avma;
2949: }
2950: ss = random_relation(phase,cglob,(long)LIMC,PRECLLL,PRECREG,
2951: nf,subFB,vecT2,mat,matarch,list_jideal);
2952: if (ss < 0)
2953: { /* could not find relations */
2954: if (ss != -1) precpb = "random_relation"; /* precision pb */
2955: goto START;
2956: }
2957: if (DEBUGLEVEL > 2) dbg_outrel(phase,cglob,vperm,mat,matarch);
2958: if (phase)
2959: for (j=lg(extramat)-1; j>0; j--)
2960: {
2961: GEN c = mat[cglob+j], *g = (GEN*) extramat[j];
2962: for (k=1; k<=KC; k++) g[k] = stoi(c[vperm[k]]);
2963: }
2964: cglob = ss;
2965: }
2966:
2967: /* reduce relation matrices */
2968: if (phase == 0)
2969: { /* never reduced before */
2970: long lgex;
2971: W = hnfspec(mat,vperm,&pdep,&B,&C, lg(subFB)-1);
2972: phase = 1;
2973: nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1;
2974: if (nlze)
2975: {
2976: list_jideal = cgetg(nlze+1, t_VECSMALL);
2977: for (i=1; i<=nlze; i++) list_jideal[i] = vperm[i];
2978: }
2979: lgex = max(nlze, MIN_EXTRA); /* set lgex > 0, in case regulator is 0 */
2980: /* allocate now, reduce dimension later (setlg) when lgex goes down */
2981: extramat= cgetg(lgex+1,t_MAT);
2982: extraC = cgetg(lgex+1,t_MAT);
2983: for (j=1; j<=lgex; j++)
2984: {
2985: extramat[j]=lgetg(KC+1,t_COL);
2986: extraC[j] = (long)cgetc_col(RU,PRECREG);
2987: }
2988: }
2989: else
2990: {
2991: if (low_stack(limpile, stack_lim(av1,1)))
2992: {
2993: GEN *gptr[6];
2994: if(DEBUGMEM>1) err(warnmem,"buchall");
2995: gptr[0]=&W; gptr[1]=&C; gptr[2]=&B; gptr[3]=&pdep;
2996: gptr[4]=&extramat; gptr[5]=&extraC; gerepilemany(av1,gptr,6);
2997: }
2998: list_jideal = NULL;
2999: W = hnfadd(W,vperm,&pdep,&B,&C, extramat,extraC);
3000: nlze = lg(pdep)>1? lg(pdep[1])-1: lg(B[1])-1;
3001: if (nlze && ++nreldep > MAXRELSUP) { sfb_increase=1; goto MORE; }
3002: }
3003: if (nlze) goto MORE; /* dependent rows */
3004:
3005: /* first attempt at regulator for the check */
3006: sreg = (lg(mat)-1) - (lg(B)-1) - (lg(W)-1); /* = zc (cf hnffinal) */
3007: xarch = cgetg(sreg+1,t_MAT); /* cols corresponding to units */
3008: for (j=1; j<=sreg; j++) xarch[j] = C[j];
3009: reg = compute_multiple_of_R(xarch,RU,N,&sublambda);
3010: if (!reg)
3011: { /* not full rank for units */
3012: if (DEBUGLEVEL) fprintferr("regulator is zero.\n");
3013: if (++nrelsup > MAXRELSUP) goto START;
3014: nlze = MIN_EXTRA; goto MORE;
3015: }
3016: /* anticipate precision problems */
3017: if (!sublambda) { precpb = "bestappr"; goto START; }
3018:
3019: /* class number */
3020: clh = dethnf_i(W);
3021: if (DEBUGLEVEL) fprintferr("\n#### Tentative class number: %Z\n", clh);
3022:
3023: /* check */
3024: z = mulrr(lfun,gmul(gmul2n(gpuigs(shiftr(mppi(DEFAULTPREC),1),-R2),-R1),
3025: gsqrt(absi(D),DEFAULTPREC)));
3026: z = mulri(z,(GEN)zu[1]);
3027: /* z = Res (zeta_K, s = 1) * w D^(1/2) / [ 2^r1 (2pi)^r2 ] = h R */
3028: p1 = gmul2n(divir(clh,z), 1);
3029: /* c1 should be close to 2, and not much smaller */
3030: c1 = compute_check(sublambda,p1,&parch,®);
3031: /* precision problems? */
3032: if (!c1 || gcmpgs(gmul2n(c1,1),3) < 0)
3033: { /* has to be a precision problem unless we cheat on Bach constant */
3034: if (!precdouble) precpb = "compute_check";
3035: goto START;
3036: }
3037: if (gcmpgs(c1,3) > 0)
3038: {
3039: if (++nrelsup <= MAXRELSUP)
3040: {
3041: if (DEBUGLEVEL) fprintferr("\n ***** check = %f\n",gtodouble(c1)/2);
3042: nlze = MIN_EXTRA; goto MORE;
3043: }
3044: if (cbach < 11.99) { sfb_increase = 1; goto MORE; }
3045: err(warner,"suspicious check. Try to increase extra relations");
3046: }
3047:
3048: if (KCZ2 > KCZ)
3049: { /* "be honest" */
3050: if (!powsubFB) powsubFBgen(nf,subFB,CBUCHG+1,PRECREG);
3051: if (!be_honest(nf,subFB,PRECLLL)) goto START;
3052: }
3053:
3054: /* fundamental units */
3055: if (flun < 0 || flun > 1)
3056: {
3057: xarch = cleancol(gmul(xarch,parch),N,PRECREG);
3058: if (DEBUGLEVEL) msgtimer("cleancol");
3059: }
3060: if (labs(flun) > 1)
3061: {
3062: fu = getfu(nf,&xarch,reg,flun,&k,PRECREG);
3063: if (k <= 0 && labs(flun) > 2)
3064: {
3065: if (k < 0)
3066: {
3067: precadd = (DEFAULTPREC-2) + ((-k) >> TWOPOTBITS_IN_LONG);
3068: if (precadd <= 0) precadd = (DEFAULTPREC-2);
3069: }
3070: precpb = "getfu"; goto START;
3071: }
3072: }
3073:
3074: /* class group generators */
3075: i = lg(C)-sreg; C += sreg; C[0] = evaltyp(t_MAT)|evallg(i);
3076: C = cleancol(C,N,PRECREG);
3077: class_group_gen(nf,W,C,vperm, &clg1, &clg2, PRECREG);
3078:
3079: c1 = gdiv(gmul(reg,clh), z);
3080: desallocate(&mat);
3081:
3082: return gerepileupto(av, buchall_end(nf,CHANGE,flun,k,fu,clg1,clg2,reg,c1,zu,W,B,xarch,C,vectbase,vperm));
3083: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to buch2.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