1 #include "defs.h"
2 #include "permfns.h"
3
4 extern char mult,conv,gens,inf1[],inf2[],outf1[],outf2[],outf3[],outf4[];
5 extern short perm[],sv[],cp[],fpt[],orb[],igno[],
6 base[],lorb[],pno[],*pptr[],*svptr[],
7 gno[],ngno[],power[],wt[],d1[],d2[],
8 facord[],pinv[],rel[],
9 *relptr[],mp,mpt,mb,mexp;
10 extern int psp,svsp;
11 short npt,np,npt1,nb,exp,prime;
12 FILE *ip,*op;
13
14 int
resetsv(void)15 resetsv (void)
16 /* This recomputes the Schreier vectors, after a change of generators */
17 { short i,j,lo,bno,b;
18 for (bno=1;bno<=nb;bno++) if (lorb[bno]>1)
19 { *pno=0; lo=0; b=base[bno];
20 for (i=1;i<=np;i++)
21 { j=gno[i];
22 if (pptr[j][npt1]>=bno)
23 { (*pno)++; pno[*pno]=j; lo=orbitsv(b,svptr[bno],lo); }
24 }
25 if (lo!=lorb[bno]) {fprintf(stderr,"Orbit length wrong.\n"); return(-1); }
26 }
27 return(0);
28 }
29
30 int
setfixb(short * p)31 setfixb (short *p)
32 { short bno; bno=1;
33 while (bno<=nb && p[base[bno]]==base[bno]) bno++;
34 p[npt1]=bno;
35 }
36
37 int
pcprog(void)38 pcprog (void)
39 { short i,j,k,l,m,n,mxp,ngads,pairs,class,hgen,coeff,nf,nf1,nf2,pt,*pk,
40 *pnf,*pnf1,*pnf2,*ipnf,*ipnf1,*p1,*p2,*ip1,*ip2,nwt;
41 int quot;
42 char diff;
43 if ((ip=fopen(inf1,"r"))==0)
44 { fprintf(stderr,"Cannot open %s.\n",inf1); return(-1); }
45 fscanf(ip,"%hd%hd%hd%hd",&npt,&np,&nb,&l);
46 if (npt>mpt) {fprintf(stderr,"npt too big. Increase MPT.\n"); return(-1); }
47 if (nb>=mb) {fprintf(stderr,"nb too big. Increase MB.\n"); return(-1); }
48 if (nb*npt>svsp)
49 { fprintf(stderr,"Out of sv space. Increase SVSP.\n");return(-1);}
50 if (l!=4) { fprintf(stderr,"Wrong input format.\n"); return(-1); }
51 readbaselo(nb,base,lorb); npt1=npt+1;
52 quot=psp/npt1; if (quot>mp) quot=mp; mxp=quot; mxp=2*(mxp/2);
53 for (i=0;i<mxp;i++) pptr[i]=perm+i*npt1-1;
54 for (i=1;i<=nb;i++) svptr[i]=sv+(i-1)*npt-1;
55 if (2*(np+1)>mxp)
56 { fprintf(stderr,"Out of perm space. Increase PSP (or MP).\n"); return(-1); }
57 readpsv(0,nb,np,svptr);
58 fscanf(ip,"%hd",&prime);
59 for (i=1;i<=np;i++)
60 { fscanf(ip,"%hd",facord+i); gno[i]=2*(i-1); igno[2*i-1]=i; }
61 fclose(ip);
62
63 /* Since we will be changing generators, we link the perm nos with fpt.
64 gno[i] will be the current gen no i in the subnormal series for P.
65 */
66 nf=2*np; for (i=nf;i<mxp;i+=2) fpt[i]=i+2; fpt[mxp-2]= -1;
67
68 /* Now we start to convert the subnormal series to a central series */
69 for (i=1;i<=np;i++)
70 { p1=pptr[gno[i]]; m=facord[i];
71 /* If m>prime, we introduce p-th powers of this generator, until all indices
72 in the series have order p.
73 */
74 while (m>prime)
75 { pnf=pptr[nf]; ipnf=pnf+npt1;
76 for (n=1;n<=npt;n++)
77 { pt=n; for (l=1;l<=prime;l++) pt=p1[pt]; pnf[n]=pt; ipnf[pt]=n; }
78 for (n=np+1;n>=i+1;n--)
79 { l=gno[n-1]; gno[n]=l; igno[l+1]=n; facord[n+1]=facord[n]; }
80 setfixb(pnf); facord[i+1]=prime; m/=prime;
81 facord[i]=m; gno[i]=nf; igno[nf+1]=i; np++; p1=pnf; nf=fpt[nf];
82 if (nf== -1)
83 { fprintf(stderr,"Out of perm space. Increase PSP (or MP).\n");
84 return(-1);
85 }
86 resetsv();
87 }
88 ip1=p1+npt1;
89 /* In the main bit of the algorithm, we compute commutator [i,j], and change j
90 until this lies in G(j-1).
91 */
92 for (j=1;j<=(i-2);j++)
93 { p2=pptr[gno[j]];
94 while(1)
95 { ip2=p2+npt1; pnf=pptr[nf]; ipnf=pnf+npt1;
96 for (n=1;n<=npt;n++)
97 { l=p1[p2[ip1[ip2[n]]]]; pnf[n]=l; ipnf[l]=n; }
98 setfixb(pnf); firstgen(pnf,&hgen,&coeff); k=hgen;
99 if (k>=i) { fprintf(stderr,"Series is not subnormal.\n"); return(-1); }
100 if (k<j) break;
101 nf1=fpt[nf]; pnf1=pptr[nf1]; ipnf1=pnf1+npt1; nf2=fpt[nf1];
102 if (nf1== -1 || nf2== -1)
103 { fprintf(stderr,"Out of perm space. Increase PSP (or MP).\n");
104 return(-1);
105 }
106 pk=pptr[gno[k]];
107 for (n=1;n<=npt;n++)
108 { pt=n; for (l=1;l<=coeff;l++) pt=pk[pt]; pnf1[n]=pt; ipnf1[pt]=n; }
109 setfixb(pnf1); fpt[gno[k]]=nf2; nf2=gno[k]; fpt[nf1]=nf2;
110 pnf2=pptr[nf2];
111 while(1)
112 { for (n=1;n<=npt;n++) pnf2[n]=pnf[ipnf1[n]];
113 firstgen(pnf2,&hgen,&coeff);
114 if (hgen>=k) {fprintf(stderr,"Impossible error!\n"); return(-1); }
115 if (hgen<j)
116 { for (m=k;m>j;m--) {l=gno[m-1]; gno[m]=l; igno[l+1]=m; }
117 gno[j]=nf1; p2=pnf1; igno[nf1+1]=j;
118 fpt[nf]=nf2; resetsv(); break;
119 }
120 for (m=k;m>=(hgen+2);m--) {l=gno[m-1];gno[m]=l;igno[l+1]=m;}
121 if (pptr[gno[hgen]][npt1]>=pnf1[npt1])
122 { l=gno[hgen]; gno[hgen+1]=l; igno[l+1]=hgen+1; }
123 else
124 { gno[hgen+1]=nf1; igno[nf1+1]=hgen+1; nf1=gno[hgen]; fpt[nf1]=nf2;
125 fpt[nf]=nf1; pnf1=pptr[nf1];
126 }
127 pk=pptr[gno[hgen]+1];
128 for (n=1;n<=npt;n++)
129 { pt=n; for (l=1;l<=coeff;l++) pt=pk[pt];pt=ipnf1[pt]; pnf1[pt]=n;}
130 ipnf1=pnf1+npt1; invert(pnf1,ipnf1); setfixb(pnf1);
131 if (hgen==j)
132 { p2=pnf1; gno[j]=nf1; igno[nf1+1]=j; fpt[nf]=nf2; resetsv(); break;}
133 k=hgen;
134 } /* while(1) */
135 } /* while(1) */
136 } /* j loop */
137 } /* i loop */
138 printf("Central series found.\n");
139
140 /* We are now ready to compute the pcp, with definitions (d1,d2) and weigths
141 (wt). This involves two stages. In the first, we compute all commutators
142 and powers, to determine what the weights and definitios are, and change
143 some of the generators if necessary. We do not change the Schreier vector
144 generators - the new generators are used as PCP generators. Whenever we
145 change a generator, we have to completely restart this process. In the
146 second stage, the PCP is computed and output.
147 */
148 exp=np; setpinv();
149 if (exp>=mexp)
150 { fprintf(stderr,"exp too big. Increase MEXP.\n"); return(-1); }
151 for (i=1;i<=exp;i++)
152 { ngno[i]=gno[i]; wt[i]=1; power[i]=1; d1[i]=0; d2[i]=0; }
153 ngads=0; class=1;
154
155 restart:
156 for (i=exp;i>1;i--) for (j=exp;j>=i;j--)
157 { pnf=pptr[nf]; ipnf=pnf+npt1; p1=pptr[ngno[i]]; ip1=p1+npt1;
158 if (i==j)
159 { if (mult) continue;
160 nwt=wt[i]+1;
161 for (n=1;n<=npt;n++)
162 { pt=n; for (l=1;l<=prime;l++) pt=p1[pt]; pnf[n]=pt; ipnf[pt]=n; }
163 }
164 else
165 { nwt=wt[i]+wt[j]; p2=pptr[ngno[j]]; ip2=p2+npt1;
166 for (n=1;n<=npt;n++) {pt=p2[p1[ip2[ip1[n]]]]; pnf[n]=pt; ipnf[pt]=n; }
167 }
168 hgen=express(pnf,rel,nwt);
169 if (hgen>0)
170 { wt[hgen]=nwt; if (nwt>class) class=nwt; l=ngno[hgen];
171 diff=0; p1=pptr[l];
172 for (k=1;k<=npt;k++) if (pnf[k]!=p1[k]) {diff=1; break;}
173 if (diff)
174 { if (l!=gno[hgen]) { fpt[l]=fpt[nf]; fpt[nf]=l;} else ngads++;
175 ngno[hgen]=nf; pptr[nf][npt1]=hgen; nf=fpt[nf];
176 if (nf== -1)
177 { fprintf(stderr,"Out of perm space. Increase PSP (or MP).\n");
178 return(-1);
179 }
180 for (k=1;k<=exp;k++) { d1[k]=0; d2[k]=0;}
181 goto restart;
182 }
183 else { d1[hgen]= exp+1-i; d2[hgen]=exp+1-j;}
184 }
185 }
186
187 /* Now we are ready to compute the PCP. We output it as we go along. */
188 op=fopen(outf1,"w");
189 fprintf(op,"%4d%4d%4d%4d\n",npt,np,nb,2);
190 printbaselo(nb,base,lorb);
191 *pno=0; *gno=np; printpsv(nb,gno,svptr);
192 for (i=exp;i>0;i--) fprintf(op,"%4d",wt[i]); fprintf(op,"\n");
193 fprintf(op,"%4d%4d\n",prime,ngads);
194 for (i=1;i<=exp;i++) fprintf(op,"%4d",power[i]); fprintf(op,"\n");
195 for (i=1;i<=exp;i++)
196 { j=ngno[i]; if (j!=gno[i]) printvec(pptr[j],1);}
197 fclose(op);
198 /* Output pcp gens as permutations if required. */
199 if (gens)
200 { op=fopen(outf4,"w");
201 fprintf(op,"%4d%4d%4d%4d\n",npt,exp,0,0);
202 for (i=exp;i>=1;i--) { j=ngno[i]; printvec(pptr[j],0);}
203 fclose(op);
204 }
205 op=fopen(outf2,"w");
206 fprintf(op,"%4d%4d%4d%4d%4d%4d\n",prime,exp,exp,exp-1,class,mult);
207 for (i=exp;i>0;i--) fprintf(op,"%4d",wt[i]); fprintf(op,"\n");
208 for (i=exp;i>0;i--) fprintf(op,"%4d",d1[i]); fprintf(op,"\n");
209 for (i=exp;i>0;i--) fprintf(op,"%4d",d2[i]); fprintf(op,"\n");
210 for (i=2;i<exp;i++) for (j=1;j<i;j++)
211 { pnf=pptr[nf]; ipnf=pnf+npt1; p1=pptr[ngno[exp+1-i]]; ip1=p1+npt1;
212 p2=pptr[ngno[exp+1-j]]; ip2=p2+npt1;
213 for (n=1;n<=npt;n++) {pt=p2[p1[ip2[ip1[n]]]]; pnf[n]=pt; ipnf[pt]=n; }
214 express(pnf,rel+1,0);
215 l= *(rel+ *(rel+1)); m=1+exp-l;
216 if (d1[m]==i && d2[m]==j) *rel=l; else *rel=0;
217 for (k=0;k<= *(rel+1)+1;k++) fprintf(op,"%4d",rel[k]); fprintf(op,"\n");
218 }
219 for (i=1;i<exp;i++)
220 { pnf=pptr[nf]; ipnf=pnf+npt1; p1=pptr[ngno[exp+1-i]]; ip1=p1+npt1;
221 for (n=1;n<=npt;n++)
222 { pt=n; for (l=1;l<=prime;l++) pt=p1[pt]; pnf[n]=pt; ipnf[pt]=n; }
223 express(pnf,rel+1,0);
224 l= *(rel+ *(rel+1)); m=1+exp-l;
225 if (d1[m]==i && d2[m]==i) *rel=l; else *rel=0;
226 for (k=0;k<= *(rel+1)+1;k++) fprintf(op,"%4d",rel[k]); fprintf(op,"\n");
227 }
228 fclose(op);
229 if (conv)
230 { if ((ip=fopen(inf2,"r"))==0)
231 { fprintf(stderr,"Cannot open %s.\n",inf2); return(-1); }
232 op=fopen(outf3,"w");
233 fscanf(ip,"%hd%hd%hd%hd",&k,&l,&m,&n);
234 if (k!=npt || n<=0)
235 { fprintf(stderr,"inf2 has invalid npt or format.\n"); return(-1);}
236 for (i=1;i<=3;i++) while (getc(ip)!='\n');
237 fprintf(op,"%3d\n",l); pnf=pptr[nf]; ipnf=pnf+npt1;
238 for (i=1;i<=l;i++)
239 { readvec(pnf,1); invert(pnf,ipnf); express(pnf,rel,0);
240 fprintf(op,"%4d ",*rel);
241 for (j=1;j<= *rel;j++) fprintf(op,"%4d",rel[j]);
242 fprintf(op,"\n");
243 }
244 }
245 return(0);
246 }
247