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