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