Annotation of OpenXM_contrib/pari/src/basemath/subgroup.c, Revision 1.1
1.1 ! maekawa 1: /* $Id: subgroup.c,v 1.1.1.1 1999/09/16 13:47:37 karim Exp $ */
! 2: #include "pari.h"
! 3: GEN hnf0(GEN x, long r);
! 4: void push_val(entree *ep, GEN a);
! 5: void pop_val(entree *ep);
! 6:
! 7: /* SUBGROUPS
! 8: * Assume: G = Gp x I, with Gp a p-group and (|I|,p)=1, and I small.
! 9: * Compute subgroups of I by recursive calls
! 10: * Loop through subgroups Hp of Gp using Birkhoff's algorithm.
! 11: * If (I is non trivial)
! 12: * lift Hp to G (mul by exponent of I)
! 13: * for each subgp of I, lift it to G (mult by exponent of Gp)
! 14: * consider the group generated by the two subgroups (concat)
! 15: */
! 16: static long *powerlist, *mmu, *lam, *c, *maxc, *a, *maxa, **g, **maxg;
! 17: static GEN **H, subq, subqpart, hnfgroup;
! 18: static GEN BINMAT;
! 19: static long countsub, expoI;
! 20: static long *available, indexbound, lsubq, lsubqpart;
! 21: static char *gpch;
! 22: static entree *ep;
! 23: static void(*treatsub_fun)(GEN);
! 24: typedef struct slist {
! 25: struct slist *next;
! 26: long *data;
! 27: } slist;
! 28:
! 29: static slist *sublist;
! 30:
! 31: void
! 32: printtyp(long *typ)
! 33: {
! 34: long i;
! 35: for (i=1; i<=typ[0]; i++) fprintferr(" %ld ",typ[i]);
! 36: fprintferr("\n");
! 37: }
! 38:
! 39: /* compute conjugate partition of typ */
! 40: static long*
! 41: conjugate(long *typ)
! 42: {
! 43: long *t, i, k = typ[0], l, last;
! 44:
! 45: if (!k) { t = new_chunk(1); t[0]=0; return t; }
! 46: l = typ[1]; t = new_chunk(l+2);
! 47: t[1] = k; last = k;
! 48: for (i=2; i<=l; i++)
! 49: {
! 50: while (typ[last] < i) last--;
! 51: t[i] = last;
! 52: }
! 53: t[i] = 0; t[0] = l;
! 54: return t;
! 55: }
! 56:
! 57: static void
! 58: std_fun(GEN x)
! 59: {
! 60: ep->value = (void*)x;
! 61: lisseq(gpch); countsub++;
! 62: }
! 63:
! 64: void
! 65: addcell(GEN H)
! 66: {
! 67: long *pt,i,j, k = 0, n = lg(H)-1;
! 68: slist *cell = (slist*) gpmalloc(sizeof(slist) + n*(n+1)/2 * sizeof(long));
! 69:
! 70: sublist->next = cell; cell->data = pt = (long*) (cell + 1);
! 71: for (j=1; j<=n; j++)
! 72: for(i=1; i<=j; i++) pt[k++] = itos(gcoeff(H,i,j));
! 73: sublist = cell;
! 74: }
! 75:
! 76: static void
! 77: list_fun(GEN x)
! 78: {
! 79: addcell(hnf(concatsp(hnfgroup,x))); countsub++;
! 80: }
! 81:
! 82: /* treat subgroup Hp (not in HNF, treatsub_fun should do it if desired) */
! 83: static void
! 84: treatsub(GEN Hp)
! 85: {
! 86: long i;
! 87: if (!subq) treatsub_fun(Hp);
! 88: else
! 89: { /* not a p group, add the trivial part */
! 90: Hp = gmulsg(expoI,Hp); /* lift Hp to G */
! 91: for (i=1; i<lsubqpart; i++)
! 92: treatsub_fun(concatsp(Hp, (GEN)subqpart[i]));
! 93: }
! 94: }
! 95:
! 96: /* assume t>0 and l>1 */
! 97: static void
! 98: dogroup()
! 99: {
! 100: long av = avma,av1, e,i,j,k,r,n,t2,ind, t = mmu[0], l = lam[0];
! 101:
! 102: t2 = (l==t)? t-1: t;
! 103: n = t2 * l - (t2*(t2+1))/2; /* number of gamma_ij */
! 104: for (i=1, r=t+1; ; i++)
! 105: {
! 106: if (available[i]) c[r++] = i;
! 107: if (r > l) break;
! 108: }
! 109: if (DEBUGLEVEL>2) { fprintferr(" column selection:"); printtyp(c); }
! 110: /* a/g and maxa/maxg access the same data indexed differently */
! 111: for (ind=0,i=1; i<=t; ind+=(l-i),i++)
! 112: {
! 113: maxg[i] = maxa + (ind - (i+1)); /* only access maxg[i][i+1..l] */
! 114: g[i] = a + (ind - (i+1));
! 115: for (r=i+1; r<=l; r++)
! 116: if (c[r] < c[i]) maxg[i][r] = powerlist[mmu[i]-mmu[r]-1];
! 117: else if (lam[c[r]] < mmu[i]) maxg[i][r] = powerlist[lam[c[r]]-mmu[r]];
! 118: else maxg[i][r] = powerlist[mmu[i]-mmu[r]];
! 119: }
! 120: av1=avma; a[n-1]=0; for (i=0; i<n-1; i++) a[i]=1;
! 121: for(;;)
! 122: {
! 123: a[n-1]++;
! 124: if (a[n-1] > maxa[n-1])
! 125: {
! 126: j=n-2; while (j>=0 && a[j]==maxa[j]) j--;
! 127: if (j<0) { avma=av; return; }
! 128: a[j]++; for (k=j+1; k<n; k++) a[k]=1;
! 129: }
! 130: for (i=1; i<=t; i++)
! 131: {
! 132: for (r=1; r<i; r++) affsi(0, H[i][c[r]]);
! 133: affsi(powerlist[lam[c[r]]-mmu[r]], H[r][c[r]]);
! 134: for (r=i+1; r<=l; r++)
! 135: {
! 136: if (c[r] < c[i])
! 137: e = g[i][r]*powerlist[lam[c[r]]-mmu[i]+1];
! 138: else
! 139: if (lam[c[r]] < mmu[i]) e = g[i][r];
! 140: else e = g[i][r]*powerlist[lam[c[r]]-mmu[i]];
! 141: affsi(e, H[i][c[r]]);
! 142: }
! 143: }
! 144: treatsub((GEN)H); avma=av1;
! 145: }
! 146: }
! 147:
! 148: /* c[1],...,c[r-1] filled */
! 149: static void
! 150: loop(long r)
! 151: {
! 152: long j;
! 153:
! 154: if (r > mmu[0]) { dogroup(); return; }
! 155:
! 156: if (r!=1 && (mmu[r-1] == mmu[r])) j = c[r-1]+1; else j = 1;
! 157: for ( ; j<=maxc[r]; j++)
! 158: if (available[j])
! 159: {
! 160: c[r] = j; available[j] = 0;
! 161: loop(r+1); available[j] = 1;
! 162: }
! 163: }
! 164:
! 165: static void
! 166: dopsubtyp()
! 167: {
! 168: long av = avma, i,r, l = lam[0], t = mmu[0];
! 169:
! 170: if (!t)
! 171: {
! 172: GEN p1 = cgetg(2,t_MAT);
! 173: p1[1] = (long)zerocol(l);
! 174: treatsub(p1); avma=av; return;
! 175: }
! 176: if (l==1) /* imply t = 1 */
! 177: {
! 178: GEN p1 = gtomat(stoi(powerlist[lam[1]-mmu[1]]));
! 179: treatsub(p1); avma=av; return;
! 180: }
! 181: c = new_chunk(l+1); c[0] = l;
! 182: maxc = new_chunk(l+1);
! 183: available = new_chunk(l+1);
! 184: a = new_chunk(l*(t+1));
! 185: maxa=new_chunk(l*(t+1));
! 186: g = (long**)new_chunk(t+1);
! 187: maxg = (long**)new_chunk(t+1);
! 188:
! 189: if (DEBUGLEVEL) { fprintferr(" subgroup:"); printtyp(mmu); }
! 190: for (i=1; i<=t; i++)
! 191: {
! 192: for (r=1; r<=l; r++)
! 193: if (mmu[i] > lam[r]) break;
! 194: maxc[i] = r-1;
! 195: }
! 196: H = (GEN**)cgetg(t+1, t_MAT);
! 197: for (i=1; i<=t; i++)
! 198: {
! 199: H[i] = (GEN*)cgetg(l+1, t_COL);
! 200: for (r=1; r<=l; r++) H[i][r] = cgeti(3);
! 201: }
! 202: for (i=1; i<=l; i++) available[i]=1;
! 203: for (i=1; i<=t; i++) c[i]=0;
! 204: /* go through all column selections */
! 205: loop(1); avma=av; return;
! 206: }
! 207:
! 208: static long
! 209: weight(long *typ)
! 210: {
! 211: long i,w = 0;
! 212: for (i=1; i<=typ[0]; i++) w += typ[i];
! 213: return w;
! 214: }
! 215:
! 216: static long
! 217: dopsub(long p, long *gtyp, long *indexsubq)
! 218: {
! 219: long w,i,j,k,n, wg, wmin = 0, count = 0;
! 220:
! 221: if (DEBUGLEVEL) { fprintferr("\ngroup:"); printtyp(gtyp); }
! 222: if (indexbound)
! 223: {
! 224: wg = weight(gtyp);
! 225: wmin = (long) (wg - (log((double)indexbound) / log((double)p)));
! 226: if (cmpis(gpuigs(stoi(p), wg - wmin), indexbound) > 0) wmin++;
! 227: }
! 228: lam = gtyp; n = lam[0]; mmu = new_chunk(n+1);
! 229: mmu[1] = -1; for (i=2; i<=n; i++) mmu[i]=0;
! 230: for(;;) /* go through all vectors mu_{i+1} <= mu_i <= lam_i */
! 231: {
! 232: mmu[1]++;
! 233: if (mmu[1] > lam[1])
! 234: {
! 235: for (j=2; j<=n; j++)
! 236: if (mmu[j] < lam[j] && mmu[j] < mmu[j-1]) break;
! 237: if (j>n) return count;
! 238: mmu[j]++; for (k=1; k<j; k++) mmu[k]=mmu[j];
! 239: }
! 240: for (j=1; j<=n; j++)
! 241: if (!mmu[j]) break;
! 242: mmu[0] = j-1; w = weight(mmu);
! 243: if (w >= wmin)
! 244: {
! 245: GEN p1 = gun;
! 246:
! 247: if (subq) /* G not a p-group */
! 248: {
! 249: if (indexbound)
! 250: {
! 251: long indexH = itos(gpuigs(stoi(p), wg - w));
! 252: long bound = indexbound / indexH;
! 253: subqpart = cgetg(lsubq, t_VEC);
! 254: lsubqpart = 1;
! 255: for (i=1; i<lsubq; i++)
! 256: if (indexsubq[i] <= bound) subqpart[lsubqpart++] = subq[i];
! 257: }
! 258: else { subqpart = subq; lsubqpart = lsubq; }
! 259: }
! 260: if (DEBUGLEVEL)
! 261: {
! 262: long *lp = conjugate(lam);
! 263: long *mp = conjugate(mmu);
! 264:
! 265: if (DEBUGLEVEL > 3)
! 266: {
! 267: fprintferr(" lambda = "); printtyp(lam);
! 268: fprintferr(" lambda' = "); printtyp(lp);
! 269: fprintferr(" mu = "); printtyp(mmu);
! 270: fprintferr(" mu' = "); printtyp(mp);
! 271: }
! 272: for (j=1; j<=mp[0]; j++)
! 273: {
! 274: p1 = mulii(p1, gpuigs(stoi(p), mp[j+1]*(lp[j]-mp[j])));
! 275: p1 = mulii(p1, gcoeff(BINMAT, lp[j]-mp[j+1]+1, mp[j]-mp[j+1]+1));
! 276: }
! 277: fprintferr(" alpha_lambda(mu,p) = %Z\n",p1);
! 278: }
! 279: countsub = 0;
! 280: dopsubtyp();
! 281: count += countsub;
! 282: if (DEBUGLEVEL)
! 283: {
! 284: fprintferr(" countsub = %ld\n", countsub);
! 285: msgtimer("for this type");
! 286: if (subq) p1 = mulis(p1,lsubqpart-1);
! 287: if (cmpis(p1,countsub))
! 288: {
! 289: fprintferr(" alpha = %Z\n",p1);
! 290: err(bugparier,"forsubgroup (alpha != countsub)");
! 291: }
! 292: }
! 293: }
! 294: }
! 295: }
! 296:
! 297: static GEN
! 298: expand_sub(GEN x, long n)
! 299: {
! 300: long i,j, m = lg(x);
! 301: GEN p = idmat(n-1), q,c;
! 302:
! 303: for (i=1; i<m; i++)
! 304: {
! 305: q = (GEN)p[i]; c = (GEN)x[i];
! 306: for (j=1; j<m; j++) q[j] = c[j];
! 307: for ( ; j<n; j++) q[j] = zero;
! 308: }
! 309: return p;
! 310: }
! 311:
! 312: GEN matqpascal(long n, GEN q);
! 313:
! 314: static long
! 315: subgroup_engine(GEN cyc, long bound)
! 316: {
! 317: long av=avma,i,j,k,imax,nbprim,count, n = lg(cyc);
! 318: GEN gtyp,fa,junk,primlist,p,listgtyp,indexsubq = NULL;
! 319: long oindexbound = indexbound;
! 320: long oexpoI = expoI;
! 321: long *opowerlist = powerlist;
! 322: GEN osubq = subq;
! 323: GEN oBINMAT = BINMAT;
! 324: long olsubq = lsubq;
! 325:
! 326: long *ommu=mmu, *olam=lam, *oc=c, *omaxc=maxc, *oa=a, *omaxa=maxa, **og=g, **omaxg=maxg;
! 327: GEN **oH=H, osubqpart=subqpart;
! 328: long ocountsub=countsub;
! 329: long *oavailable=available, olsubqpart=lsubqpart;
! 330:
! 331: if (typ(cyc) != t_VEC)
! 332: {
! 333: if (typ(cyc) != t_MAT) err(typeer,"forsubgroup");
! 334: cyc = mattodiagonal(cyc);
! 335: }
! 336: for (i=1; i<n-1; i++)
! 337: if (!divise((GEN)cyc[i], (GEN)cyc[i+1]))
! 338: err(talker,"not a group in forsubgroup");
! 339: if (n == 1 || gcmp1((GEN)cyc[1])) { treatsub(cyc); return 1; }
! 340: if (!signe(cyc[1]))
! 341: err(talker,"infinite group in forsubgroup");
! 342: if (DEBUGLEVEL) timer2();
! 343: indexbound = bound;
! 344: fa = factor((GEN)cyc[1]); primlist = (GEN)fa[1];
! 345: nbprim = lg(primlist);
! 346: listgtyp = new_chunk(n); k = 0;
! 347: for (i=1; i<nbprim; i++)
! 348: {
! 349: gtyp = new_chunk(n); p = (GEN)primlist[i];
! 350: for (j=1; j<n; j++)
! 351: {
! 352: gtyp[j] = pvaluation((GEN)cyc[j], p, &junk);
! 353: if (!gtyp[j]) break;
! 354: }
! 355: j--; gtyp[0] = j;
! 356: if (j > k) { k = j; imax = i; }
! 357: listgtyp[i] = (long)gtyp;
! 358: }
! 359: gtyp = (GEN)listgtyp[imax]; p = (GEN)primlist[imax];
! 360: k = gtyp[1];
! 361: powerlist = new_chunk(k+1); powerlist[0] = 1;
! 362: powerlist[1] = itos(p);
! 363: for (j=1; j<=k; j++) powerlist[j] = powerlist[1] * powerlist[j-1];
! 364:
! 365: if (DEBUGLEVEL) BINMAT = matqpascal(gtyp[0]+1, p);
! 366: if (nbprim == 2) subq = NULL;
! 367: else
! 368: { /* not a p-group */
! 369: GEN cyc2 = dummycopy(cyc);
! 370: GEN ohnfgroup = hnfgroup;
! 371: for (i=1; i<n; i++)
! 372: {
! 373: cyc2[i] = ldivis((GEN)cyc2[i], powerlist[gtyp[i]]);
! 374: if (gcmp1((GEN)cyc2[i])) break;
! 375: }
! 376: setlg(cyc2, i);
! 377: if (is_bigint(cyc[1]))
! 378: err(impl,"subgrouplist for large cyclic factors");
! 379: expoI = itos((GEN)cyc2[1]);
! 380: hnfgroup = diagonal(cyc2);
! 381: subq = subgrouplist(cyc2, bound);
! 382: hnfgroup = ohnfgroup;
! 383: lsubq = lg(subq);
! 384: for (i=1; i<lsubq; i++)
! 385: subq[i] = (long)expand_sub((GEN)subq[i], n);
! 386: if (indexbound)
! 387: {
! 388: indexsubq = new_chunk(lsubq);
! 389: for (i=1; i<lsubq; i++)
! 390: indexsubq[i] = itos(dethnf((GEN)subq[i]));
! 391: }
! 392: /* lift subgroups of I to G */
! 393: for (i=1; i<lsubq; i++)
! 394: subq[i] = lmulsg(powerlist[k],(GEN)subq[i]);
! 395: if (DEBUGLEVEL>2)
! 396: {
! 397: fprintferr("(lifted) subgp of prime to %Z part:\n",p);
! 398: outbeaut(subq);
! 399: }
! 400: }
! 401: count = dopsub(powerlist[1],gtyp,indexsubq);
! 402: if (DEBUGLEVEL) fprintferr("nb subgroup = %ld\n",count);
! 403:
! 404: mmu=ommu; lam=olam; c=oc; maxc=omaxc; a=oa; maxa=omaxa; g=og; maxg=omaxg;
! 405: H=oH; subqpart=osubqpart,
! 406: countsub=ocountsub;
! 407: available=oavailable; lsubqpart=olsubqpart;
! 408:
! 409: indexbound = oindexbound;
! 410: expoI = oexpoI;
! 411: powerlist = opowerlist;
! 412: subq = osubq;
! 413: BINMAT = oBINMAT;
! 414: lsubq = olsubq;
! 415: avma=av; return count;
! 416: }
! 417:
! 418: void
! 419: forsubgroup(entree *oep, GEN cyc, long bound, char *och)
! 420: {
! 421: entree *saveep = ep;
! 422: char *savech = gpch;
! 423: void(*savefun)(GEN) = treatsub_fun;
! 424: long n, N = lg(cyc)-1;
! 425:
! 426: treatsub_fun = std_fun;
! 427: cyc = dummycopy(cyc);
! 428: for (n = N; n > 1; n--) /* take care of trailing 1s */
! 429: if (!gcmp1((GEN)cyc[n])) break;
! 430: setlg(cyc, n+1);
! 431: gpch = och;
! 432: ep = oep;
! 433: push_val(ep, gzero);
! 434: (void)subgroup_engine(cyc,bound);
! 435: pop_val(ep);
! 436: treatsub_fun = savefun;
! 437: gpch = savech;
! 438: ep = saveep;
! 439: }
! 440:
! 441: GEN
! 442: subgrouplist(GEN cyc, long bound)
! 443: {
! 444: void(*savefun)(GEN) = treatsub_fun;
! 445: slist *olist = sublist, *list;
! 446: long ii,i,j,k,nbsub,n, N = lg(cyc)-1, av = avma;
! 447: GEN z,H;
! 448: GEN ohnfgroup = hnfgroup;
! 449:
! 450: sublist = list = (slist*) gpmalloc(sizeof(slist));
! 451: treatsub_fun = list_fun;
! 452: cyc = dummycopy(cyc);
! 453: for (n = N; n > 1; n--) /* take care of trailing 1s */
! 454: if (!gcmp1((GEN)cyc[n])) break;
! 455: setlg(cyc, n+1);
! 456: hnfgroup = diagonal(cyc);
! 457: nbsub = subgroup_engine(cyc,bound);
! 458: hnfgroup = ohnfgroup; avma = av;
! 459: z = cgetg(nbsub+1,t_VEC); sublist = list;
! 460: for (ii=1; ii<=nbsub; ii++)
! 461: {
! 462: list = sublist; sublist = list->next; free(list);
! 463: H = cgetg(N+1,t_MAT); z[ii]=(long)H; k=0;
! 464: for (j=1; j<=n; j++)
! 465: {
! 466: H[j] = lgetg(N+1, t_COL);
! 467: for (i=1; i<=j; i++) coeff(H,i,j) = lstoi(sublist->data[k++]);
! 468: for ( ; i<=N; i++) coeff(H,i,j) = zero;
! 469: }
! 470: for ( ; j<=N; j++)
! 471: {
! 472: H[j] = lgetg(N+1, t_COL);
! 473: for (i=1; i<=N; i++) coeff(H,i,j) = (i==j)? un: zero;
! 474: }
! 475: }
! 476: free(sublist); sublist = olist;
! 477: treatsub_fun = savefun; return z;
! 478: }
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 / src / basemath