Annotation of OpenXM_contrib/pari/src/basemath/subgroup.c, Revision 1.1.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