Annotation of OpenXM_contrib/pari-2.2/src/basemath/buch2.c, Revision 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