Annotation of OpenXM_contrib/pari-2.2/src/basemath/subgroup.c, Revision 1.2
1.2 ! noro 1: /* $Id: subgroup.c,v 1.25 2002/04/19 20:13:01 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: #include "pari.h"
17: extern GEN hnf0(GEN x, long r);
18: void push_val(entree *ep, GEN a);
19: void pop_val(entree *ep);
20:
1.2 ! noro 21: typedef struct {
! 22: entree *ep;
! 23: char *s;
! 24: } expr_t;
! 25:
! 26: typedef struct slist {
! 27: struct slist *next;
! 28: long *data;
! 29: } slist;
! 30:
! 31: typedef struct {
! 32: GEN hnfgroup;
! 33: GEN listKer;
! 34: ulong count;
! 35: slist *list;
! 36: } sublist_t;
! 37:
! 38: /* MAX: [G:H] <= bound, EXACT: [G:H] = bound, TYPE: type(H) = bound */
! 39: enum { b_NONE, b_MAX, b_EXACT, b_TYPE };
! 40:
1.1 noro 41: /* SUBGROUPS
1.2 ! noro 42: * G = Gp x I, with Gp a p-Sylow (I assumed small).
1.1 noro 43: * Compute subgroups of I by recursive calls
44: * Loop through subgroups Hp of Gp using Birkhoff's algorithm.
45: * If (I is non trivial)
46: * lift Hp to G (mul by exponent of I)
47: * for each subgp of I, lift it to G (mult by exponent of Gp)
48: * consider the group generated by the two subgroups (concat)
1.2 ! noro 49: *
! 50: * type(H) = mu --> H = Z/p^mu[1] x ... x Z/p^mu[len(mu)] */
! 51: typedef struct subgp_iter {
! 52: long *M, *L; /* mu = p-subgroup type, lambda = p-group type */
! 53: long *powlist; /* [i] = p^i, i = 0.. */
! 54: long *c, *maxc, *a, *maxa, **g, **maxg;
! 55: long *available;
! 56: GEN **H; /* p-subgroup of type mu, in matrix form */
! 57: GEN cyc; /* cyclic factors of G */
! 58: GEN subq;/* subgrouplist(I) */
! 59: GEN subqpart; /* J in subq s.t [I:J][Gp:Hp] <= indexbound */
! 60: GEN bound; /* if != NULL, impose a "bound" on [G:H] (see boundtype) */
! 61: int boundtype;
! 62: long countsub; /* number of subgroups of type M (so far) */
! 63: long count; /* number of p-subgroups so far [updated when M completed] */
! 64: GEN expoI; /* exponent of I */
! 65: void(*fun)(struct subgp_iter *, GEN); /* applied to each subgroup */
! 66: void *fundata; /* for fun */
! 67: } subgp_iter;
1.1 noro 68:
1.2 ! noro 69: #define len(x) (x)[0]
! 70: #define setlen(x,l) len(x)=(l)
1.1 noro 71:
1.2 ! noro 72: void
! 73: printtyp(const long *typ)
1.1 noro 74: {
1.2 ! noro 75: long i, l = len(typ);
! 76: for (i=1; i<=l; i++) fprintferr(" %ld ",typ[i]);
1.1 noro 77: fprintferr("\n");
78: }
79:
80: /* compute conjugate partition of typ */
81: static long*
82: conjugate(long *typ)
83: {
1.2 ! noro 84: long *t, i, k = len(typ), l, last;
1.1 noro 85:
1.2 ! noro 86: if (!k) { t = new_chunk(1); setlen(t,0); return t; }
1.1 noro 87: l = typ[1]; t = new_chunk(l+2);
88: t[1] = k; last = k;
89: for (i=2; i<=l; i++)
90: {
91: while (typ[last] < i) last--;
92: t[i] = last;
93: }
1.2 ! noro 94: t[i] = 0; setlen(t,l);
1.1 noro 95: return t;
96: }
1.2 ! noro 97: /* --- subgp_iter 'fun' associated to forsubgroup -------------- */
1.1 noro 98: static void
1.2 ! noro 99: std_fun(subgp_iter *T, GEN x)
1.1 noro 100: {
1.2 ! noro 101: expr_t *E = (expr_t *)T->fundata;
! 102: E->ep->value = (void*)x;
! 103: (void)lisseq(E->s); T->countsub++;
1.1 noro 104: }
1.2 ! noro 105: /* ----subgp_iter 'fun' associated to subgrouplist ------------- */
! 106: static void
! 107: addcell(sublist_t *S, GEN H)
1.1 noro 108: {
109: long *pt,i,j, k = 0, n = lg(H)-1;
110: slist *cell = (slist*) gpmalloc(sizeof(slist) + n*(n+1)/2 * sizeof(long));
111:
1.2 ! noro 112: S->list->next = cell; cell->data = pt = (long*) (cell + 1);
1.1 noro 113: for (j=1; j<=n; j++)
114: for(i=1; i<=j; i++) pt[k++] = itos(gcoeff(H,i,j));
1.2 ! noro 115: S->list = cell;
! 116: S->count++;
! 117: }
! 118:
! 119: extern int hnfdivide(GEN A, GEN B);
! 120:
! 121: /* 1 if h^(-1)*list[i] integral for some i, 0 otherwise */
! 122: static int
! 123: hnflistdivise(GEN list,GEN h)
! 124: {
! 125: long i, n = lg(list);
! 126: for (i=1; i<n; i++)
! 127: if (hnfdivide(h, (GEN)list[i])) break;
! 128: return i < n;
1.1 noro 129: }
130:
131: static void
1.2 ! noro 132: list_fun(subgp_iter *T, GEN x)
1.1 noro 133: {
1.2 ! noro 134: sublist_t *S = (sublist_t*)T->fundata;
! 135: GEN H = hnf(concatsp(S->hnfgroup,x));
! 136: if (S->listKer && hnflistdivise(S->listKer, H)) return; /* test conductor */
! 137: addcell(S, H); T->countsub++;
1.1 noro 138: }
1.2 ! noro 139: /* -------------------------------------------------------------- */
1.1 noro 140:
1.2 ! noro 141: /* treat subgroup Hp (not in HNF, T->fun should do it if desired) */
1.1 noro 142: static void
1.2 ! noro 143: treatsub(subgp_iter *T, GEN H)
1.1 noro 144: {
145: long i;
1.2 ! noro 146: if (!T->subq) T->fun(T, H);
1.1 noro 147: else
148: { /* not a p group, add the trivial part */
1.2 ! noro 149: GEN Hp = gmul(T->expoI, H); /* lift H to G */
! 150: long l = lg(T->subqpart);
! 151: for (i=1; i<l; i++)
! 152: T->fun(T, concatsp(Hp, (GEN)T->subqpart[i]));
1.1 noro 153: }
154: }
155:
156: /* assume t>0 and l>1 */
157: static void
1.2 ! noro 158: dogroup(subgp_iter *T)
1.1 noro 159: {
1.2 ! noro 160: const long *powlist = T->powlist;
! 161: long *M = T->M;
! 162: long *L = T->L;
! 163: long *c = T->c;
! 164: long *a = T->a, *maxa = T->maxa;
! 165: long **g = T->g, **maxg = T->maxg;
! 166: GEN **H = T->H;
! 167: gpmem_t av = avma;
! 168: long e,i,j,k,r,n,t2,ind, t = len(M), l = len(L);
1.1 noro 169:
170: t2 = (l==t)? t-1: t;
171: n = t2 * l - (t2*(t2+1))/2; /* number of gamma_ij */
172: for (i=1, r=t+1; ; i++)
173: {
1.2 ! noro 174: if (T->available[i]) c[r++] = i;
1.1 noro 175: if (r > l) break;
176: }
177: if (DEBUGLEVEL>2) { fprintferr(" column selection:"); printtyp(c); }
178: /* a/g and maxa/maxg access the same data indexed differently */
179: for (ind=0,i=1; i<=t; ind+=(l-i),i++)
180: {
181: maxg[i] = maxa + (ind - (i+1)); /* only access maxg[i][i+1..l] */
182: g[i] = a + (ind - (i+1));
183: for (r=i+1; r<=l; r++)
1.2 ! noro 184: if (c[r] < c[i]) maxg[i][r] = powlist[M[i]-M[r]-1];
! 185: else if (L[c[r]] < M[i]) maxg[i][r] = powlist[L[c[r]]-M[r]];
! 186: else maxg[i][r] = powlist[M[i]-M[r]];
1.1 noro 187: }
1.2 ! noro 188: av = avma; a[n-1]=0; for (i=0; i<n-1; i++) a[i]=1;
1.1 noro 189: for(;;)
190: {
191: a[n-1]++;
192: if (a[n-1] > maxa[n-1])
193: {
194: j=n-2; while (j>=0 && a[j]==maxa[j]) j--;
1.2 ! noro 195: if (j < 0) { avma = av; return; }
! 196:
1.1 noro 197: a[j]++; for (k=j+1; k<n; k++) a[k]=1;
198: }
199: for (i=1; i<=t; i++)
200: {
201: for (r=1; r<i; r++) affsi(0, H[i][c[r]]);
1.2 ! noro 202: affsi(powlist[L[c[r]] - M[r]], H[r][c[r]]);
1.1 noro 203: for (r=i+1; r<=l; r++)
204: {
1.2 ! noro 205: if (c[r] < c[i])
! 206: e = g[i][r] * powlist[L[c[r]] - M[i]+1];
1.1 noro 207: else
1.2 ! noro 208: if (L[c[r]] < M[i]) e = g[i][r];
! 209: else e = g[i][r] * powlist[L[c[r]] - M[i]];
1.1 noro 210: affsi(e, H[i][c[r]]);
211: }
212: }
1.2 ! noro 213: treatsub(T, (GEN)H); avma = av;
1.1 noro 214: }
215: }
216:
1.2 ! noro 217: /* T->c[1],...,T->c[r-1] filled */
1.1 noro 218: static void
1.2 ! noro 219: loop(subgp_iter *T, long r)
1.1 noro 220: {
221: long j;
222:
1.2 ! noro 223: if (r > len(T->M)) { dogroup(T); return; }
1.1 noro 224:
1.2 ! noro 225: if (r!=1 && (T->M[r-1] == T->M[r])) j = T->c[r-1]+1; else j = 1;
! 226: for ( ; j<=T->maxc[r]; j++)
! 227: if (T->available[j])
1.1 noro 228: {
1.2 ! noro 229: T->c[r] = j; T->available[j] = 0;
! 230: loop(T, r+1); T->available[j] = 1;
1.1 noro 231: }
232: }
233:
234: static void
1.2 ! noro 235: dopsubtyp(subgp_iter *T)
1.1 noro 236: {
1.2 ! noro 237: gpmem_t av = avma;
! 238: long i,r, l = len(T->L), t = len(T->M);
1.1 noro 239:
240: if (!t)
241: {
242: GEN p1 = cgetg(2,t_MAT);
1.2 ! noro 243: p1[1] = (long)zerocol(l);
! 244: treatsub(T, p1); avma = av; return;
1.1 noro 245: }
246: if (l==1) /* imply t = 1 */
247: {
1.2 ! noro 248: GEN p1 = gtomat(stoi(T->powlist[T->L[1]-T->M[1]]));
! 249: treatsub(T, p1); avma = av; return;
1.1 noro 250: }
1.2 ! noro 251: T->c = new_chunk(l+1); setlen(T->c, l);
! 252: T->maxc = new_chunk(l+1);
! 253: T->available = new_chunk(l+1);
! 254: T->a = new_chunk(l*(t+1));
! 255: T->maxa= new_chunk(l*(t+1));
! 256: T->g = (long**)new_chunk(t+1);
! 257: T->maxg = (long**)new_chunk(t+1);
1.1 noro 258:
1.2 ! noro 259: if (DEBUGLEVEL) { fprintferr(" subgroup:"); printtyp(T->M); }
1.1 noro 260: for (i=1; i<=t; i++)
261: {
262: for (r=1; r<=l; r++)
1.2 ! noro 263: if (T->M[i] > T->L[r]) break;
! 264: T->maxc[i] = r-1;
1.1 noro 265: }
1.2 ! noro 266: T->H = (GEN**)cgetg(t+1, t_MAT);
1.1 noro 267: for (i=1; i<=t; i++)
268: {
1.2 ! noro 269: T->H[i] = (GEN*)cgetg(l+1, t_COL);
! 270: for (r=1; r<=l; r++) T->H[i][r] = cgeti(3);
1.1 noro 271: }
1.2 ! noro 272: for (i=1; i<=l; i++) T->available[i]=1;
! 273: for (i=1; i<=t; i++) T->c[i]=0;
1.1 noro 274: /* go through all column selections */
1.2 ! noro 275: loop(T, 1); avma = av; return;
1.1 noro 276: }
277:
278: static long
279: weight(long *typ)
280: {
1.2 ! noro 281: long i, l = len(typ), w = 0;
! 282: for (i=1; i<=l; i++) w += typ[i];
1.1 noro 283: return w;
284: }
285:
1.2 ! noro 286: static void
! 287: dopsub(subgp_iter *T, GEN p, GEN indexsubq)
1.1 noro 288: {
1.2 ! noro 289: long *M, *L = T->L;
! 290: long w,i,j,k, wG = weight(L), wmin = 0, wmax = wG, n = len(L);
1.1 noro 291:
1.2 ! noro 292: if (DEBUGLEVEL) { fprintferr("\ngroup:"); printtyp(L); }
! 293: T->count = 0;
! 294: switch(T->boundtype)
! 295: {
! 296: case b_MAX: /* upper bound */
! 297: wmin = (long) (wG - (log(gtodouble(T->bound)) / log(gtodouble(p))));
! 298: if (cmpii(gpowgs(p, wG - wmin), T->bound) > 0) wmin++;
! 299: break;
! 300: case b_EXACT: /* exact value */
! 301: wmin = wmax = wG - ggval(T->bound, p);
! 302: break;
1.1 noro 303: }
1.2 ! noro 304: T->M = M = new_chunk(n+1);
! 305: M[1] = -1; for (i=2; i<=n; i++) M[i]=0;
1.1 noro 306: for(;;) /* go through all vectors mu_{i+1} <= mu_i <= lam_i */
307: {
1.2 ! noro 308: M[1]++;
! 309: if (M[1] > L[1])
1.1 noro 310: {
311: for (j=2; j<=n; j++)
1.2 ! noro 312: if (M[j] < L[j] && M[j] < M[j-1]) break;
! 313: if (j > n) return;
! 314:
! 315: M[j]++; for (k=1; k<j; k++) M[k]=M[j];
1.1 noro 316: }
317: for (j=1; j<=n; j++)
1.2 ! noro 318: if (!M[j]) break;
! 319: setlen(M, j-1); w = weight(M);
! 320: if (w >= wmin && w <= wmax)
! 321: {
1.1 noro 322: GEN p1 = gun;
323:
1.2 ! noro 324: if (T->subq) /* G not a p-group */
1.1 noro 325: {
1.2 ! noro 326: if (T->bound)
1.1 noro 327: {
1.2 ! noro 328: GEN indexH = gpowgs(p, wG - w);
! 329: GEN B = divii(T->bound, indexH);
! 330: long k, l = lg(T->subq);
! 331: T->subqpart = cgetg(l, t_VEC);
! 332: k = 1;
! 333: for (i=1; i<l; i++)
! 334: if (cmpii((GEN)indexsubq[i], B) <= 0)
! 335: T->subqpart[k++] = T->subq[i];
! 336: setlg(T->subqpart, k);
1.1 noro 337: }
1.2 ! noro 338: else T->subqpart = T->subq;
1.1 noro 339: }
340: if (DEBUGLEVEL)
341: {
1.2 ! noro 342: long *Lp = conjugate(L);
! 343: long *Mp = conjugate(M);
! 344: GEN BINMAT = matqpascal(len(L)+1, p);
1.1 noro 345:
346: if (DEBUGLEVEL > 3)
347: {
1.2 ! noro 348: fprintferr(" lambda = "); printtyp(L);
! 349: fprintferr(" lambda'= "); printtyp(Lp);
! 350: fprintferr(" mu = "); printtyp(M);
! 351: fprintferr(" mu'= "); printtyp(Mp);
1.1 noro 352: }
1.2 ! noro 353: for (j=1; j<=len(Mp); j++)
1.1 noro 354: {
1.2 ! noro 355: p1 = mulii(p1, gpuigs(p, Mp[j+1]*(Lp[j]-Mp[j])));
! 356: p1 = mulii(p1, gcoeff(BINMAT, Lp[j]-Mp[j+1]+1, Mp[j]-Mp[j+1]+1));
1.1 noro 357: }
358: fprintferr(" alpha_lambda(mu,p) = %Z\n",p1);
359: }
1.2 ! noro 360: T->countsub = 0;
! 361:
! 362: dopsubtyp(T);
! 363:
! 364: T->count += T->countsub;
1.1 noro 365: if (DEBUGLEVEL)
366: {
1.2 ! noro 367: fprintferr(" countsub = %ld\n", T->countsub);
1.1 noro 368: msgtimer("for this type");
1.2 ! noro 369: if (T->fun != list_fun || !((sublist_t*)(T->fundata))->listKer)
1.1 noro 370: {
1.2 ! noro 371: if (T->subq) p1 = mulis(p1,lg(T->subqpart)-1);
! 372: if (cmpis(p1,T->countsub))
! 373: {
! 374: fprintferr(" alpha = %Z\n",p1);
! 375: err(bugparier,"forsubgroup (alpha != countsub)");
! 376: }
1.1 noro 377: }
378: }
379: }
380: }
381: }
382:
1.2 ! noro 383: void
! 384: parse_bound(subgp_iter *T)
! 385: {
! 386: GEN b, B = T->bound;
! 387: if (!B) { T->boundtype = b_NONE; return; }
! 388: switch(typ(B))
! 389: {
! 390: case t_INT: /* upper bound */
! 391: T->boundtype = b_MAX;
! 392: break;
! 393: case t_VEC: /* exact value */
! 394: b = (GEN)B[1];
! 395: if (lg(B) != 2 || typ(b) != t_INT) err(typeer,"subgroup");
! 396: T->boundtype = b_EXACT;
! 397: T->bound = b;
! 398: break;
! 399: case t_COL: /* exact type */
! 400: if (lg(B) > len(T->L)+1) err(typeer,"subgroup");
! 401: err(impl,"exact type in subgrouplist");
! 402: T->boundtype = b_TYPE;
! 403: break;
! 404: default: err(typeer,"subgroup");
! 405: }
! 406: if (signe(T->bound) <= 0)
! 407: err(talker,"subgroup: index bound must be positive");
! 408: }
! 409:
1.1 noro 410: static GEN
411: expand_sub(GEN x, long n)
412: {
413: long i,j, m = lg(x);
414: GEN p = idmat(n-1), q,c;
415:
416: for (i=1; i<m; i++)
417: {
418: q = (GEN)p[i]; c = (GEN)x[i];
419: for (j=1; j<m; j++) q[j] = c[j];
420: for ( ; j<n; j++) q[j] = zero;
421: }
422: return p;
423: }
424:
425: extern GEN matqpascal(long n, GEN q);
1.2 ! noro 426: static GEN subgrouplist_i(GEN cyc, GEN bound, GEN expoI, GEN listKer);
! 427:
! 428: static GEN
! 429: init_powlist(long k, long p)
! 430: {
! 431: GEN z = new_chunk(k+1);
! 432: long i;
! 433: z[0] = 1; z[1] = (long)p;
! 434: for (i=1; i<=k; i++)
! 435: {
! 436: GEN t = muluu(p, z[i-1]);
! 437: z[i] = itos(t);
! 438: }
! 439: return z;
! 440: }
1.1 noro 441:
1.2 ! noro 442: static void
! 443: subgroup_engine(subgp_iter *T)
1.1 noro 444: {
1.2 ! noro 445: gpmem_t av = avma;
! 446: GEN B,L,fa,junk,primlist,p,listL,indexsubq = NULL;
! 447: GEN cyc = T->cyc;
! 448: long i,j,k,imax,nbprim, n = lg(cyc);
! 449:
! 450: if (typ(cyc) != t_VEC)
1.1 noro 451: {
452: if (typ(cyc) != t_MAT) err(typeer,"forsubgroup");
453: cyc = mattodiagonal(cyc);
454: }
455: for (i=1; i<n-1; i++)
456: if (!divise((GEN)cyc[i], (GEN)cyc[i+1]))
457: err(talker,"not a group in forsubgroup");
1.2 ! noro 458: if (n == 1) { T->fun(T, cyc); return; }
1.1 noro 459: if (!signe(cyc[1]))
460: err(talker,"infinite group in forsubgroup");
1.2 ! noro 461: if (DEBUGLEVEL) (void)timer2();
1.1 noro 462: fa = factor((GEN)cyc[1]); primlist = (GEN)fa[1];
463: nbprim = lg(primlist);
1.2 ! noro 464: listL = new_chunk(n); imax = k = 0;
1.1 noro 465: for (i=1; i<nbprim; i++)
466: {
1.2 ! noro 467: L = new_chunk(n); p = (GEN)primlist[i];
1.1 noro 468: for (j=1; j<n; j++)
469: {
1.2 ! noro 470: L[j] = pvaluation((GEN)cyc[j], p, &junk);
! 471: if (!L[j]) break;
1.1 noro 472: }
1.2 ! noro 473: j--; setlen(L, j);
1.1 noro 474: if (j > k) { k = j; imax = i; }
1.2 ! noro 475: listL[i] = (long)L;
1.1 noro 476: }
1.2 ! noro 477: L = (GEN)listL[imax]; p = (GEN)primlist[imax];
! 478: k = L[1];
! 479: T->L = L;
! 480: T->powlist = init_powlist(k, itos(p));
! 481: B = T->bound;
! 482: parse_bound(T);
1.1 noro 483:
1.2 ! noro 484: if (nbprim == 2)
! 485: {
! 486: T->subq = NULL;
! 487: if (T->boundtype == b_EXACT)
! 488: {
! 489: (void)pvaluation(T->bound,p,&B);
! 490: if (!gcmp1(B)) { avma = av; return; }
! 491: }
! 492: }
1.1 noro 493: else
494: { /* not a p-group */
1.2 ! noro 495: GEN cycI = dummycopy(cyc);
! 496: long lsubq;
1.1 noro 497: for (i=1; i<n; i++)
498: {
1.2 ! noro 499: cycI[i] = ldivis((GEN)cycI[i], T->powlist[L[i]]);
! 500: if (gcmp1((GEN)cycI[i])) break;
1.1 noro 501: }
1.2 ! noro 502: setlg(cycI, i); /* cyclic factors of I */
! 503: if (T->boundtype == b_EXACT)
1.1 noro 504: {
1.2 ! noro 505: (void)pvaluation(T->bound,p,&B);
! 506: B = _vec(B);
! 507: }
! 508: T->expoI = (GEN)cycI[1];
! 509: T->subq = subgrouplist_i(cycI, B, T->expoI, NULL);
! 510:
! 511: lsubq = lg(T->subq);
! 512: for (i=1; i<lsubq; i++)
! 513: T->subq[i] = (long)expand_sub((GEN)T->subq[i], n);
! 514: if (T->bound)
! 515: {
! 516: indexsubq = cgetg(lsubq,t_VEC);
1.1 noro 517: for (i=1; i<lsubq; i++)
1.2 ! noro 518: indexsubq[i] = (long)dethnf_i((GEN)T->subq[i]);
1.1 noro 519: }
520: /* lift subgroups of I to G */
1.2 ! noro 521: for (i=1; i<lsubq; i++)
! 522: T->subq[i] = lmulsg(T->powlist[k],(GEN)T->subq[i]);
1.1 noro 523: if (DEBUGLEVEL>2)
524: {
525: fprintferr("(lifted) subgp of prime to %Z part:\n",p);
1.2 ! noro 526: outbeaut(T->subq);
1.1 noro 527: }
528: }
1.2 ! noro 529: dopsub(T, p,indexsubq);
! 530: if (DEBUGLEVEL) fprintferr("nb subgroup = %ld\n",T->count);
! 531: avma = av;
1.1 noro 532: }
533:
534: static GEN
1.2 ! noro 535: get_snf(GEN x, long *N)
1.1 noro 536: {
537: GEN cyc;
538: long n;
539: switch(typ(x))
540: {
541: case t_MAT:
542: if (!isdiagonal(x)) return NULL;
543: cyc = mattodiagonal_i(x); break;
1.2 ! noro 544: case t_VEC: if (lg(x) == 4 && typ(x[2]) == t_VEC) x = (GEN)x[2];
1.1 noro 545: case t_COL: cyc = dummycopy(x); break;
546: default: return NULL;
547: }
1.2 ! noro 548: *N = lg(cyc)-1;
! 549: for (n = *N; n > 0; n--) /* take care of trailing 1s */
1.1 noro 550: {
551: GEN c = (GEN)cyc[n];
552: if (typ(c) != t_INT) return NULL;
553: if (!gcmp1(c)) break;
554: }
555: setlg(cyc, n+1);
556: for ( ; n > 0; n--)
557: {
558: GEN c = (GEN)cyc[n];
559: if (typ(c) != t_INT) return NULL;
560: }
561: return cyc;
562: }
563:
564: void
1.2 ! noro 565: forsubgroup(entree *ep, GEN cyc, GEN bound, char *ch)
1.1 noro 566: {
1.2 ! noro 567: subgp_iter T;
! 568: expr_t E;
! 569: long N;
1.1 noro 570:
1.2 ! noro 571: T.fun = &std_fun;
! 572: cyc = get_snf(cyc,&N);
1.1 noro 573: if (!cyc) err(typeer,"forsubgroup");
1.2 ! noro 574: T.bound = bound;
! 575: T.cyc = cyc;
! 576: E.s = ch;
! 577: E.ep= ep; T.fundata = (void*)&E;
1.1 noro 578: push_val(ep, gzero);
1.2 ! noro 579:
! 580: subgroup_engine(&T);
! 581:
1.1 noro 582: pop_val(ep);
583: }
584:
1.2 ! noro 585: static GEN
! 586: subgrouplist_i(GEN cyc, GEN bound, GEN expoI, GEN listKer)
1.1 noro 587: {
1.2 ! noro 588: gpmem_t av = avma;
! 589: subgp_iter T;
! 590: sublist_t S;
! 591: slist *list, *sublist;
! 592: long ii,i,j,k,nbsub,n,N;
1.1 noro 593: GEN z,H;
594:
1.2 ! noro 595: cyc = get_snf(cyc, &N);
1.1 noro 596: if (!cyc) err(typeer,"subgrouplist");
1.2 ! noro 597: n = lg(cyc)-1; /* not necessarily = N */
! 598:
! 599: S.list = sublist = (slist*) gpmalloc(sizeof(slist));
! 600: S.hnfgroup = diagonal(cyc);
! 601: S.listKer = listKer;
! 602: S.count = 0;
! 603: T.fun = &list_fun;
! 604: T.fundata = (void*)&S;
! 605:
! 606: T.cyc = cyc;
! 607: T.bound = bound;
! 608: T.expoI = expoI;
! 609:
! 610: subgroup_engine(&T);
! 611:
! 612: nbsub = S.count;
! 613: avma = av;
! 614: z = cgetg(nbsub+1,t_VEC);
1.1 noro 615: for (ii=1; ii<=nbsub; ii++)
616: {
617: list = sublist; sublist = list->next; free(list);
618: H = cgetg(N+1,t_MAT); z[ii]=(long)H; k=0;
619: for (j=1; j<=n; j++)
620: {
621: H[j] = lgetg(N+1, t_COL);
622: for (i=1; i<=j; i++) coeff(H,i,j) = lstoi(sublist->data[k++]);
623: for ( ; i<=N; i++) coeff(H,i,j) = zero;
624: }
625: for ( ; j<=N; j++)
626: {
627: H[j] = lgetg(N+1, t_COL);
628: for (i=1; i<=N; i++) coeff(H,i,j) = (i==j)? un: zero;
629: }
630: }
1.2 ! noro 631: free(sublist); return z;
! 632: }
! 633:
! 634: GEN
! 635: subgrouplist(GEN cyc, GEN bound)
! 636: {
! 637: return subgrouplist_i(cyc,bound,NULL,NULL);
! 638: }
! 639:
! 640: GEN
! 641: subgroupcondlist(GEN cyc, GEN bound, GEN listKer)
! 642: {
! 643: return subgrouplist_i(cyc,bound,NULL,listKer);
1.1 noro 644: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to subgroup.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