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