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