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