File: [local] / OpenXM_contrib / pari / src / basemath / Attic / subgroup.c (download)
Revision 1.1.1.1 (vendor branch), Sun Jan 9 17:35:31 2000 UTC (24 years, 8 months ago) by maekawa
Branch: PARI_GP
CVS Tags: maekawa-ipv6, VERSION_2_0_17_BETA, RELEASE_20000124, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, RELEASE_1_1_3, RELEASE_1_1_2 Changes since 1.1: +0 -0
lines
Import PARI/GP 2.0.17 beta.
|
/* $Id: subgroup.c,v 1.1.1.1 1999/09/16 13:47:37 karim Exp $ */
#include "pari.h"
GEN hnf0(GEN x, long r);
void push_val(entree *ep, GEN a);
void pop_val(entree *ep);
/* SUBGROUPS
* Assume: G = Gp x I, with Gp a p-group and (|I|,p)=1, and I small.
* Compute subgroups of I by recursive calls
* Loop through subgroups Hp of Gp using Birkhoff's algorithm.
* If (I is non trivial)
* lift Hp to G (mul by exponent of I)
* for each subgp of I, lift it to G (mult by exponent of Gp)
* consider the group generated by the two subgroups (concat)
*/
static long *powerlist, *mmu, *lam, *c, *maxc, *a, *maxa, **g, **maxg;
static GEN **H, subq, subqpart, hnfgroup;
static GEN BINMAT;
static long countsub, expoI;
static long *available, indexbound, lsubq, lsubqpart;
static char *gpch;
static entree *ep;
static void(*treatsub_fun)(GEN);
typedef struct slist {
struct slist *next;
long *data;
} slist;
static slist *sublist;
void
printtyp(long *typ)
{
long i;
for (i=1; i<=typ[0]; i++) fprintferr(" %ld ",typ[i]);
fprintferr("\n");
}
/* compute conjugate partition of typ */
static long*
conjugate(long *typ)
{
long *t, i, k = typ[0], l, last;
if (!k) { t = new_chunk(1); t[0]=0; return t; }
l = typ[1]; t = new_chunk(l+2);
t[1] = k; last = k;
for (i=2; i<=l; i++)
{
while (typ[last] < i) last--;
t[i] = last;
}
t[i] = 0; t[0] = l;
return t;
}
static void
std_fun(GEN x)
{
ep->value = (void*)x;
lisseq(gpch); countsub++;
}
void
addcell(GEN H)
{
long *pt,i,j, k = 0, n = lg(H)-1;
slist *cell = (slist*) gpmalloc(sizeof(slist) + n*(n+1)/2 * sizeof(long));
sublist->next = cell; cell->data = pt = (long*) (cell + 1);
for (j=1; j<=n; j++)
for(i=1; i<=j; i++) pt[k++] = itos(gcoeff(H,i,j));
sublist = cell;
}
static void
list_fun(GEN x)
{
addcell(hnf(concatsp(hnfgroup,x))); countsub++;
}
/* treat subgroup Hp (not in HNF, treatsub_fun should do it if desired) */
static void
treatsub(GEN Hp)
{
long i;
if (!subq) treatsub_fun(Hp);
else
{ /* not a p group, add the trivial part */
Hp = gmulsg(expoI,Hp); /* lift Hp to G */
for (i=1; i<lsubqpart; i++)
treatsub_fun(concatsp(Hp, (GEN)subqpart[i]));
}
}
/* assume t>0 and l>1 */
static void
dogroup()
{
long av = avma,av1, e,i,j,k,r,n,t2,ind, t = mmu[0], l = lam[0];
t2 = (l==t)? t-1: t;
n = t2 * l - (t2*(t2+1))/2; /* number of gamma_ij */
for (i=1, r=t+1; ; i++)
{
if (available[i]) c[r++] = i;
if (r > l) break;
}
if (DEBUGLEVEL>2) { fprintferr(" column selection:"); printtyp(c); }
/* a/g and maxa/maxg access the same data indexed differently */
for (ind=0,i=1; i<=t; ind+=(l-i),i++)
{
maxg[i] = maxa + (ind - (i+1)); /* only access maxg[i][i+1..l] */
g[i] = a + (ind - (i+1));
for (r=i+1; r<=l; r++)
if (c[r] < c[i]) maxg[i][r] = powerlist[mmu[i]-mmu[r]-1];
else if (lam[c[r]] < mmu[i]) maxg[i][r] = powerlist[lam[c[r]]-mmu[r]];
else maxg[i][r] = powerlist[mmu[i]-mmu[r]];
}
av1=avma; a[n-1]=0; for (i=0; i<n-1; i++) a[i]=1;
for(;;)
{
a[n-1]++;
if (a[n-1] > maxa[n-1])
{
j=n-2; while (j>=0 && a[j]==maxa[j]) j--;
if (j<0) { avma=av; return; }
a[j]++; for (k=j+1; k<n; k++) a[k]=1;
}
for (i=1; i<=t; i++)
{
for (r=1; r<i; r++) affsi(0, H[i][c[r]]);
affsi(powerlist[lam[c[r]]-mmu[r]], H[r][c[r]]);
for (r=i+1; r<=l; r++)
{
if (c[r] < c[i])
e = g[i][r]*powerlist[lam[c[r]]-mmu[i]+1];
else
if (lam[c[r]] < mmu[i]) e = g[i][r];
else e = g[i][r]*powerlist[lam[c[r]]-mmu[i]];
affsi(e, H[i][c[r]]);
}
}
treatsub((GEN)H); avma=av1;
}
}
/* c[1],...,c[r-1] filled */
static void
loop(long r)
{
long j;
if (r > mmu[0]) { dogroup(); return; }
if (r!=1 && (mmu[r-1] == mmu[r])) j = c[r-1]+1; else j = 1;
for ( ; j<=maxc[r]; j++)
if (available[j])
{
c[r] = j; available[j] = 0;
loop(r+1); available[j] = 1;
}
}
static void
dopsubtyp()
{
long av = avma, i,r, l = lam[0], t = mmu[0];
if (!t)
{
GEN p1 = cgetg(2,t_MAT);
p1[1] = (long)zerocol(l);
treatsub(p1); avma=av; return;
}
if (l==1) /* imply t = 1 */
{
GEN p1 = gtomat(stoi(powerlist[lam[1]-mmu[1]]));
treatsub(p1); avma=av; return;
}
c = new_chunk(l+1); c[0] = l;
maxc = new_chunk(l+1);
available = new_chunk(l+1);
a = new_chunk(l*(t+1));
maxa=new_chunk(l*(t+1));
g = (long**)new_chunk(t+1);
maxg = (long**)new_chunk(t+1);
if (DEBUGLEVEL) { fprintferr(" subgroup:"); printtyp(mmu); }
for (i=1; i<=t; i++)
{
for (r=1; r<=l; r++)
if (mmu[i] > lam[r]) break;
maxc[i] = r-1;
}
H = (GEN**)cgetg(t+1, t_MAT);
for (i=1; i<=t; i++)
{
H[i] = (GEN*)cgetg(l+1, t_COL);
for (r=1; r<=l; r++) H[i][r] = cgeti(3);
}
for (i=1; i<=l; i++) available[i]=1;
for (i=1; i<=t; i++) c[i]=0;
/* go through all column selections */
loop(1); avma=av; return;
}
static long
weight(long *typ)
{
long i,w = 0;
for (i=1; i<=typ[0]; i++) w += typ[i];
return w;
}
static long
dopsub(long p, long *gtyp, long *indexsubq)
{
long w,i,j,k,n, wg, wmin = 0, count = 0;
if (DEBUGLEVEL) { fprintferr("\ngroup:"); printtyp(gtyp); }
if (indexbound)
{
wg = weight(gtyp);
wmin = (long) (wg - (log((double)indexbound) / log((double)p)));
if (cmpis(gpuigs(stoi(p), wg - wmin), indexbound) > 0) wmin++;
}
lam = gtyp; n = lam[0]; mmu = new_chunk(n+1);
mmu[1] = -1; for (i=2; i<=n; i++) mmu[i]=0;
for(;;) /* go through all vectors mu_{i+1} <= mu_i <= lam_i */
{
mmu[1]++;
if (mmu[1] > lam[1])
{
for (j=2; j<=n; j++)
if (mmu[j] < lam[j] && mmu[j] < mmu[j-1]) break;
if (j>n) return count;
mmu[j]++; for (k=1; k<j; k++) mmu[k]=mmu[j];
}
for (j=1; j<=n; j++)
if (!mmu[j]) break;
mmu[0] = j-1; w = weight(mmu);
if (w >= wmin)
{
GEN p1 = gun;
if (subq) /* G not a p-group */
{
if (indexbound)
{
long indexH = itos(gpuigs(stoi(p), wg - w));
long bound = indexbound / indexH;
subqpart = cgetg(lsubq, t_VEC);
lsubqpart = 1;
for (i=1; i<lsubq; i++)
if (indexsubq[i] <= bound) subqpart[lsubqpart++] = subq[i];
}
else { subqpart = subq; lsubqpart = lsubq; }
}
if (DEBUGLEVEL)
{
long *lp = conjugate(lam);
long *mp = conjugate(mmu);
if (DEBUGLEVEL > 3)
{
fprintferr(" lambda = "); printtyp(lam);
fprintferr(" lambda' = "); printtyp(lp);
fprintferr(" mu = "); printtyp(mmu);
fprintferr(" mu' = "); printtyp(mp);
}
for (j=1; j<=mp[0]; j++)
{
p1 = mulii(p1, gpuigs(stoi(p), mp[j+1]*(lp[j]-mp[j])));
p1 = mulii(p1, gcoeff(BINMAT, lp[j]-mp[j+1]+1, mp[j]-mp[j+1]+1));
}
fprintferr(" alpha_lambda(mu,p) = %Z\n",p1);
}
countsub = 0;
dopsubtyp();
count += countsub;
if (DEBUGLEVEL)
{
fprintferr(" countsub = %ld\n", countsub);
msgtimer("for this type");
if (subq) p1 = mulis(p1,lsubqpart-1);
if (cmpis(p1,countsub))
{
fprintferr(" alpha = %Z\n",p1);
err(bugparier,"forsubgroup (alpha != countsub)");
}
}
}
}
}
static GEN
expand_sub(GEN x, long n)
{
long i,j, m = lg(x);
GEN p = idmat(n-1), q,c;
for (i=1; i<m; i++)
{
q = (GEN)p[i]; c = (GEN)x[i];
for (j=1; j<m; j++) q[j] = c[j];
for ( ; j<n; j++) q[j] = zero;
}
return p;
}
GEN matqpascal(long n, GEN q);
static long
subgroup_engine(GEN cyc, long bound)
{
long av=avma,i,j,k,imax,nbprim,count, n = lg(cyc);
GEN gtyp,fa,junk,primlist,p,listgtyp,indexsubq = NULL;
long oindexbound = indexbound;
long oexpoI = expoI;
long *opowerlist = powerlist;
GEN osubq = subq;
GEN oBINMAT = BINMAT;
long olsubq = lsubq;
long *ommu=mmu, *olam=lam, *oc=c, *omaxc=maxc, *oa=a, *omaxa=maxa, **og=g, **omaxg=maxg;
GEN **oH=H, osubqpart=subqpart;
long ocountsub=countsub;
long *oavailable=available, olsubqpart=lsubqpart;
if (typ(cyc) != t_VEC)
{
if (typ(cyc) != t_MAT) err(typeer,"forsubgroup");
cyc = mattodiagonal(cyc);
}
for (i=1; i<n-1; i++)
if (!divise((GEN)cyc[i], (GEN)cyc[i+1]))
err(talker,"not a group in forsubgroup");
if (n == 1 || gcmp1((GEN)cyc[1])) { treatsub(cyc); return 1; }
if (!signe(cyc[1]))
err(talker,"infinite group in forsubgroup");
if (DEBUGLEVEL) timer2();
indexbound = bound;
fa = factor((GEN)cyc[1]); primlist = (GEN)fa[1];
nbprim = lg(primlist);
listgtyp = new_chunk(n); k = 0;
for (i=1; i<nbprim; i++)
{
gtyp = new_chunk(n); p = (GEN)primlist[i];
for (j=1; j<n; j++)
{
gtyp[j] = pvaluation((GEN)cyc[j], p, &junk);
if (!gtyp[j]) break;
}
j--; gtyp[0] = j;
if (j > k) { k = j; imax = i; }
listgtyp[i] = (long)gtyp;
}
gtyp = (GEN)listgtyp[imax]; p = (GEN)primlist[imax];
k = gtyp[1];
powerlist = new_chunk(k+1); powerlist[0] = 1;
powerlist[1] = itos(p);
for (j=1; j<=k; j++) powerlist[j] = powerlist[1] * powerlist[j-1];
if (DEBUGLEVEL) BINMAT = matqpascal(gtyp[0]+1, p);
if (nbprim == 2) subq = NULL;
else
{ /* not a p-group */
GEN cyc2 = dummycopy(cyc);
GEN ohnfgroup = hnfgroup;
for (i=1; i<n; i++)
{
cyc2[i] = ldivis((GEN)cyc2[i], powerlist[gtyp[i]]);
if (gcmp1((GEN)cyc2[i])) break;
}
setlg(cyc2, i);
if (is_bigint(cyc[1]))
err(impl,"subgrouplist for large cyclic factors");
expoI = itos((GEN)cyc2[1]);
hnfgroup = diagonal(cyc2);
subq = subgrouplist(cyc2, bound);
hnfgroup = ohnfgroup;
lsubq = lg(subq);
for (i=1; i<lsubq; i++)
subq[i] = (long)expand_sub((GEN)subq[i], n);
if (indexbound)
{
indexsubq = new_chunk(lsubq);
for (i=1; i<lsubq; i++)
indexsubq[i] = itos(dethnf((GEN)subq[i]));
}
/* lift subgroups of I to G */
for (i=1; i<lsubq; i++)
subq[i] = lmulsg(powerlist[k],(GEN)subq[i]);
if (DEBUGLEVEL>2)
{
fprintferr("(lifted) subgp of prime to %Z part:\n",p);
outbeaut(subq);
}
}
count = dopsub(powerlist[1],gtyp,indexsubq);
if (DEBUGLEVEL) fprintferr("nb subgroup = %ld\n",count);
mmu=ommu; lam=olam; c=oc; maxc=omaxc; a=oa; maxa=omaxa; g=og; maxg=omaxg;
H=oH; subqpart=osubqpart,
countsub=ocountsub;
available=oavailable; lsubqpart=olsubqpart;
indexbound = oindexbound;
expoI = oexpoI;
powerlist = opowerlist;
subq = osubq;
BINMAT = oBINMAT;
lsubq = olsubq;
avma=av; return count;
}
void
forsubgroup(entree *oep, GEN cyc, long bound, char *och)
{
entree *saveep = ep;
char *savech = gpch;
void(*savefun)(GEN) = treatsub_fun;
long n, N = lg(cyc)-1;
treatsub_fun = std_fun;
cyc = dummycopy(cyc);
for (n = N; n > 1; n--) /* take care of trailing 1s */
if (!gcmp1((GEN)cyc[n])) break;
setlg(cyc, n+1);
gpch = och;
ep = oep;
push_val(ep, gzero);
(void)subgroup_engine(cyc,bound);
pop_val(ep);
treatsub_fun = savefun;
gpch = savech;
ep = saveep;
}
GEN
subgrouplist(GEN cyc, long bound)
{
void(*savefun)(GEN) = treatsub_fun;
slist *olist = sublist, *list;
long ii,i,j,k,nbsub,n, N = lg(cyc)-1, av = avma;
GEN z,H;
GEN ohnfgroup = hnfgroup;
sublist = list = (slist*) gpmalloc(sizeof(slist));
treatsub_fun = list_fun;
cyc = dummycopy(cyc);
for (n = N; n > 1; n--) /* take care of trailing 1s */
if (!gcmp1((GEN)cyc[n])) break;
setlg(cyc, n+1);
hnfgroup = diagonal(cyc);
nbsub = subgroup_engine(cyc,bound);
hnfgroup = ohnfgroup; avma = av;
z = cgetg(nbsub+1,t_VEC); sublist = list;
for (ii=1; ii<=nbsub; ii++)
{
list = sublist; sublist = list->next; free(list);
H = cgetg(N+1,t_MAT); z[ii]=(long)H; k=0;
for (j=1; j<=n; j++)
{
H[j] = lgetg(N+1, t_COL);
for (i=1; i<=j; i++) coeff(H,i,j) = lstoi(sublist->data[k++]);
for ( ; i<=N; i++) coeff(H,i,j) = zero;
}
for ( ; j<=N; j++)
{
H[j] = lgetg(N+1, t_COL);
for (i=1; i<=N; i++) coeff(H,i,j) = (i==j)? un: zero;
}
}
free(sublist); sublist = olist;
treatsub_fun = savefun; return z;
}