1 #include "defs.h"
2
3 extern char slg,con,check,inf1[],inf2[],inf3[],inf4[],inf5[],inf6[],outf1[];
4 extern int msp,psp,svsp;
5 extern short mm,mv,mp,mpt,mpr,mdim,mb,
6 mspace[],*vec[],**mat[],pinv[],perm[],sv[],cp[],
7 *pptr[],*svptr[],base[],lorb[];
8 short prime,dim,*spv,**spm,npt,nb,maxv,maxm,maxp,nmat,nm2;
9 FILE *ip,*op;
10
seeknln(void)11 void seeknln (void) { while (getc(ip)!='\n'); }
12
mcprog(void)13 int mcprog (void)
14 { short i,j,k,l,m,n,nperms,np2,*p,**q;
15 int quot;
16 if ((ip=fopen(inf1,"r"))==0)
17 { fprintf(stderr,"Cannot open %s.\n",inf1); return(-1);}
18 fscanf(ip,"%hd%hd%hd",&prime,&dim,&nmat);
19 if (prime>mpr)
20 { fprintf(stderr,"prime too big. Increase MPR.\n"); return(-1); }
21 if (dim>mdim)
22 { fprintf(stderr,"dim too big. Increase MDIM.\n"); return(-1); }
23 setpinv();
24 /* Set pointers to define vectors and matrices. Matrix elements are stored
25 in mspace, as vectors of vectors.
26 */
27 quot=msp/dim; if (quot>mv) quot=mv; maxv=quot;
28 for (i=0;i<maxv;i++) vec[i]=mspace-1+i*dim;
29 maxm=maxv/dim-1; if (maxm>=mm) maxm=mm-1;
30 for (i=0;i<=maxm;i++) mat[i]=vec-1+i*dim;
31 /* spm is a spare matrix used for inverting, etc. similarly spv */
32 spm=mat[maxm]; spv=spm[1]; nm2=nmat*2-1;
33 if (nm2>=maxm)
34 { fprintf(stderr,"Too many mats. Increase MM.\n"); return(-1); }
35 for (i=0;i<nm2;i+=2)
36 { readmat(mat[i]);
37 if (inv(mat[i],mat[i+1]) == -1) return(-1);
38 }
39 fclose(ip);
40 if ((ip=fopen(inf2,"r"))==0)
41 { fprintf(stderr,"Cannot open %s.\n",inf2); return(-1); }
42 fscanf(ip,"%hd",&i);
43 if (i!=nmat)
44 { fprintf(stderr,"nmat does not agree in %s.\n",inf2); return(-1); }
45 /* Compute matrices for strong generators of G. Then save only those in
46 gpname.sg (unless -o set)
47 */
48 while(1)
49 { fscanf(ip,"%hd",cp); if (*cp== -1) break;
50 for (i=1;i<=*cp;i++) fscanf(ip,"%hd",cp+i);
51 nmat++; nm2+=2;
52 if (nm2>=maxm)
53 { fprintf(stderr,"Too many mats. Increase MM.\n"); return(-1); }
54 prod(cp,mat[nm2-1]); inv(mat[nm2-1],mat[nm2]);
55 }
56 fclose(ip);
57 if ((ip=fopen(inf3,"r"))==0)
58 { fprintf(stderr,"Cannot open %s.\n",inf3); return(-1); }
59 fscanf(ip,"%hd%hd%hd%hd",&npt,&nperms,&nb,&l);
60 if (npt>mpt) {fprintf(stderr,"npt too big. Increase MPT.\n"); return(-1); }
61 if (nb>=mb) {fprintf(stderr,"nb too big. Increase MB.\n"); return(-1); }
62 if (nb*npt>svsp)
63 { fprintf(stderr,"svsp too small. Increase SVSP.\n"); return(-1);}
64 np2=2*nperms-1;
65 if (l<=0 || (l!=3 && slg) || np2>nm2 || (slg==0 && np2!=nm2))
66 { fprintf(stderr,"%s has invalid data nm2,np2,l=%d,%d,%d.\n",inf3,nm2,np2,l);
67 return(-1);
68 }
69 quot=psp/npt; if (quot>mp) quot=mp; maxp=quot;
70 if (np2>=maxp)
71 { fprintf(stderr,"Need more perm space. Increase PSP (or MP).\n");
72 return(-1);
73 }
74 for (i=0;i<maxp;i++) pptr[i]=perm+i*(npt+1)-1;
75 for (i=1;i<=nb;i++) svptr[i]=sv+(i-1)*npt-1;
76 readbaselo(nb,base,lorb); readpsv(0,nb,nperms,svptr);
77 if (slg)
78 { p=mspace+msp-1-nmat; for (i=1;i<=nmat;i++) p[i]=0;
79 for (i=1;i<=nperms;i++) { fscanf(ip,"%hd",cp+i); p[cp[i]]=1; }
80 fclose(ip);
81 k=0;
82 for (i=1;i<=nperms;i++)
83 { q=vec-1+(2*cp[i]-2)*dim; mat[k]=q; mat[k+1]=q+dim; k+=2; }
84 for (i=1;i<=nmat;i++) if (p[i]==0)
85 { q=vec-1+(2*i-2)*dim; mat[k]=q; mat[k+1]=q+dim; k+=2; }
86 nmat=nperms; nm2=np2;
87 /* If -t set, we check that the matrices satisfy the relations of G */
88 if (check)
89 { if ((ip=fopen(inf6,"r"))==0)
90 { fprintf(stderr,"Cannot open %s.\n",inf6); return(-1);}
91 fscanf(ip,"%hd%hd",&i,&n); seeknln();
92 for (i=1;i<=n;i++)
93 { fscanf(ip,"%hd",cp);
94 for (j=1;j<= *cp;j++) fscanf(ip,"%hd",cp+j);
95 q=mat[nm2+1]; prod(cp,q);
96 for (k=1;k<=dim;k++) for (l=1;l<=dim;l++)
97 if ((l==k && q[k][l]!=1) || (l!=k && q[k][l]!=0))
98 { fprintf(stderr,"Matrices do not satisfy group relations.\n");
99 fprintf(stderr," (Relation %d in %s.)\n",i,inf6);
100 fclose(ip); return(-1);
101 }
102 }
103 fclose(ip);
104 }
105
106 }
107 return(0);
108 }
109
110 int
conprog(int con)111 conprog (int con)
112 { short i,j,k,n,y,fac,lm3,bn,lmat,nperms,class,*d1,*d2,*wt,
113 *p,*q,*r,*cv,*cve,*cc,*sw,*bcf,*ibcf,**cbm,**icbm,**newmat;
114 int sum;
115 char id;
116 if (con==1)
117 { for (i=0;i<nm2;i+=2) {trans(mat[i],mat[i+1]); inv(mat[i+1],mat[i]);}}
118 if ((ip=fopen(inf4,"r"))==0)
119 { fprintf(stderr,"Cannot open %s.\n",inf4); return(-1); }
120 fscanf(ip,"%hd%hd%hd%hd",&i,&nperms,&j,&k); seeknln();
121 if (i!=npt) { fprintf(stderr,"npt error in %s.\n",inf4); return(-1); }
122 y= con==1 ? nperms : 1;
123 if (nm2+y>=maxp)
124 { fprintf(stderr,"Need more perm space. Increase PSP (or MP).\n");
125 return(-1);
126 }
127 y= con==1 ? nperms+4 : 1;
128 if (nm2+y>=maxm)
129 { fprintf(stderr,"Need more mat space. Increase MSP (or MM or MV).\n");
130 return(-1);
131 }
132 if (j>0) seeknln(); if (k!=0) seeknln();
133 /* In case con=1, we reverse the order of the generators, since
134 the generators of P as permutation group are in reverse order from
135 PCP generators.
136 We deal with the case con=1 separately, since we need to read in all
137 of the generators from inf4, and compute their matrices, in order to
138 compute the necessary basis changes.
139 */
140 if (con==1)
141 { for (i=nperms;i>=1;i--) {readvec(pptr[nm2+i],0); seeknln(); }
142 fclose(ip);
143
144 /* Now we express the perms in terms of the generators of G, and compute
145 their matrices.
146 */
147 for (i=1;i<=nperms;i++)
148 { p=pptr[nm2+i]; *cp=0;
149 for (j=1;j<=nb;j++) addsv(image(p[base[j]]),svptr[j]);
150 q=cp+ *cp; p=q+1; *p= *cp;
151 while (q>cp) {*(++p)= (*q%2==0) ? *q+1 : *q-1; q--;}
152 p=cp+ *cp+1; prod(p,mat[nm2+i]);
153 }
154 lmat=nm2+nperms;
155 /* Now do the basis changes.
156 The base change matrices will be output to inf6, for use later in nqrun.
157 */
158 op=fopen(inf6,"w");
159 fprintf(op,"%4d%4d%4d\n",prime,dim,2);
160 cbm=mat[lmat+1];
161 if (dim==1) { cbm[1][1]=id=1; }
162 else
163 { icbm=mat[lmat+2]; bcf=mat[lmat+4][1]; lm3=lmat+3;
164 for (i=1;i<=dim;i++) bcf[i]=0; ibcf=mat[lmat+4][2]; id=1;
165 for (n=dim;n>=1;n--)
166 { cv=mat[lm3][1]; cve=cv+dim; cc=mat[lm3][2];
167 for (i=dim;i>=1;i--) if (bcf[i]==0)
168 if (n==1) {p=cbm[1]; q=p+dim; while (++p<=q) *p=0; cbm[1][i]=1; break;}
169 else
170 { bn=i; p=cv;
171 while (++p<=cve) *p=0;
172 cv[i]=1;
173 while(1)
174 { for (j=1;j<=nperms;j++)
175 { comm(cv,cc,mat[nm2+j]);
176 for (k=dim;k>n;k--) if ((fac=cc[ibcf[k]])!=0)
177 { p=cc; r=p+dim; q=cbm[k]+1;
178 while (++p<=r)
179 { sum= *p-(fac* *q); *p=sum%prime; if (*p<0) *p+=prime; q++; }
180 }
181 for (k=dim;k>0;k--) if (cc[k]!=0)
182 { bn=k; id=0; sw=cv; cv=cc; cc=sw; cve=cv+dim; break;}
183 if (k>0) break;
184 }
185 if (j>nperms) break;
186 }
187 bcf[bn]=n; ibcf[n]=bn; p=cbm[n]; q=p+dim; r=cv+1; fac=pinv[cv[bn]];
188 while (++p<=q) {sum= *r*fac; *p=sum%prime; r++; }
189 break;
190 }
191 }
192 }
193 printmat(cbm);
194 if (id) printf("No first basis change.\n"); else
195 { /* printf("First dual basis change matrix (new in terms of old):\n");
196 for (i=1;i<=dim;i++)
197 {for (j=1;j<=dim;j++) printf("%3d",cbm[i][j]); printf("\n"); } */
198 inv(cbm,icbm); newmat=mat[lm3]; *cp=3; cp[1]=lmat+1; cp[3]=lmat+2;
199 for (i=0;i<=lmat;i++)
200 { cp[2]=i; prod(cp,newmat);
201 mat[lm3]=mat[i]; mat[i]=newmat; newmat=mat[lm3];
202 }
203 }
204 if (dim==1)
205 { d1=mat[lmat+2][1]; d2=mat[lmat+3][1];
206 wt=mat[lmat+4][1]; wt[1]=class=1; d1[1]=d2[1]=0;
207 }
208 else
209 { d1=bcf; d2=ibcf; wt= (dim==2) ? mat[lmat+5][1] : mat[lmat+4][3];
210 id=cbdef(nm2+1,nm2+nperms,lmat+1,d1,d2,wt,&class);
211 }
212 printmat(cbm);
213 fclose(op);
214 if (id) printf("No second basis change.\n"); else
215 { /* printf("Second dual basis change matrix (new in terms of old):\n");
216 for (i=1;i<=dim;i++)
217 {for (j=1;j<=dim;j++) printf("%3d",cbm[i][j]); printf("\n"); } */
218 inv(cbm,icbm); newmat=mat[lm3]; *cp=3; cp[1]=lmat+1; cp[3]=lmat+2;
219 for (i=0;i<=lmat;i++)
220 { cp[2]=i; prod(cp,newmat);
221 mat[lm3]=mat[i]; mat[i]=newmat; newmat=mat[lm3];
222 }
223 }
224 op=fopen(outf1,"w");
225 fprintf(op,"%4d%4d%4d%3d\n",prime,dim,nperms,class);
226 for (i=1;i<=dim;i++) fprintf(op,"%3d ",wt[i]); fprintf(op,"\n");
227 for (i=1;i<=dim;i++) fprintf(op,"%3d ",d1[i]); fprintf(op,"\n");
228 for (i=1;i<=dim;i++) fprintf(op,"%3d ",d2[i]); fprintf(op,"\n");
229 for (i=1;i<=nperms;i++) printmat(mat[nm2+i]);
230 fclose(op);
231 } /* con=1 */
232 else
233 { op=fopen(outf1,"w");
234 fprintf(op,"%4d%4d%4d\n",prime,dim,nperms);
235 for (i=1;i<=nperms;i++)
236 { p=pptr[nm2+1]; readvec(p,0); *cp=0; seeknln();
237 for (j=1;j<=nb;j++) addsv(image(p[base[j]]),svptr[j]);
238 if (con==0)
239 { q=cp+ *cp; p=q+1; *p= *cp;
240 while (q>cp) {*(++p)= (*q%2==0) ? *q+1 : *q-1; q--;}
241 p=cp+ *cp+1;
242 }
243 else p=cp;
244 prod(p,mat[nm2+1]); printmat(mat[nm2+1]);
245 }
246 fclose(ip); fclose(op);
247 }
248 return(0);
249 }
250
251 /* We repeat a few procedures from permfns.c, since that file is not loaded
252 in matcalc.
253 */
254 int
addsv(int pt,short * sv)255 addsv (int pt, short *sv)
256 { short pn;
257 pn=sv[pt];
258 while (pn!= -1)
259 { (*cp)++; cp[*cp]=pn; pt=pptr[pn][pt]; pn=sv[pt]; }
260 }
261
262 int
invert(short * ptr1,short * ptr2)263 invert (short *ptr1, short *ptr2)
264 { short i;
265 for (i=1;i<=npt;i++) ptr2[ptr1[i]]=i;
266 }
267
268 int
image(int pt)269 image (int pt)
270 { short i;
271 for (i=1;i<= *cp;i++) pt=pptr[cp[i]][pt];
272 return(pt);
273 }
274
275 int
readvec(short * ptr,int e)276 readvec (short *ptr, int e)
277 { short i; for (i=1;i<=npt+e;i++) fscanf(ip,"%hd",ptr+i); }
278
279 int
readbaselo(int nb,short * base,short * lorb)280 readbaselo (int nb, short *base, short *lorb)
281 { short i;
282 for (i=1;i<=nb;i++) fscanf(ip,"%hd",base+i);
283 for (i=1;i<=nb;i++) fscanf(ip,"%hd",lorb+i);
284 }
285
286 int
readpsv(int e,int nb,int nperms,short ** svptr)287 readpsv (int e, int nb, int nperms, short **svptr)
288 { short i,j,*k;
289 for (i=1;i<=nperms;i++)
290 { j=e+2*i-2; readvec(pptr[j],1); invert(pptr[j],pptr[j+1]); }
291 for (i=1;i<=nb;i++)
292 { readvec(svptr[i],0);
293 for (j=1;j<=npt;j++) { k=svptr[i]+j; if (*k>0) *k +=e;}
294 }
295 }
296
setpinv(void)297 int setpinv (void)
298 { short i,j;
299 for (i=0;i<prime;i++) pinv[i]=0;
300 for (i=1;i<prime;i++) if (pinv[i]==0) for (j=1;j<prime;j++)
301 if (i*j % prime ==1) { pinv[i]=j; pinv[j]=i; break; }
302 }
303