1 #include <stdio.h>
2 extern char inf1[],inf2[],inf3[],outf1[],outf2[],outf3[],temp1[],temp2[],
3 expg,exph,cr,dcr,triv;
4 extern int psp,trsp,svsp;
5 extern int mpt,mb,
6 tree[],perm[],gorb[],lorbg[],lorbh[],base[],
7 fpt[],bpt[],coh_index[],*cp[],*trad[],*sv[],
8 *tailad[],**cpad[],**svgptr[],**svhptr[];
9 int trptr;
10 int npt,**lcp,**ucp,*stop,dummy,*pst,*pend,nb,fnt,lnt,
11 ind,coind,oind,cind,bno,pt,**expp,npg;
12 FILE *fopen(),*ip,*op,*opx;
13
14 /* The data structures of this program differ from other permutation group
15 programs. Permutations are not numbered, but are located by their base
16 adresses. Consequently, arrays sv and cp are arrays of pointers. The
17 current permutation in cp is stored in the middle of cp (since it has to be
18 expanded in both directions), between adresses lcp and ucp of cp.
19 The next six procedures are the equivlents of the corresponding ones
20 in permfns.c, but addsv can go in both directions.
21 In this half of the program, coset reps of H in G are computed.
22 */
23
24 int
image(int pt)25 image (int pt)
26 { int **p;
27 p=lcp;
28 while (p<=ucp) {pt=(*p)[pt]; p++;}
29 return(pt);
30 }
31
32 int
addsvb(int pt,int ** sv)33 addsvb (int pt, int **sv)
34 { int *p;
35 while ((p=sv[pt])!=stop) {pt=p[pt]; lcp--; *lcp =p-npt; }
36 }
37
38 int
addsvf(int pt,int ** sv)39 addsvf (int pt, int **sv)
40 { int *p;
41 while ((p=sv[pt])!=stop) {pt=p[pt]; ucp++; *ucp=p; }
42 }
43
44 int
invert(int * a,int * b)45 invert (int *a, int *b)
46 {int i; for (i=1;i<=npt;i++) b[a[i]]=i; }
47
48 int
rdperm(int * a,int * b)49 rdperm (int *a, int *b)
50 { int i,*r;
51 for (i=1;i<=npt;i++) fscanf(ip,"%d",a+i);
52 r=pst+1; fscanf(ip,"%d",r); invert(a,b);
53 }
54
55 int
rdsv(int ** sv)56 rdsv (int **sv)
57 { int i,j;
58 for (i=1;i<=npt;i++)
59 { fscanf(ip,"%d",&j); sv[i]= (j==0) ? 0 : (j== -1) ? stop : cp[j];}
60 }
61
62 int
cnprg1(void)63 cnprg1 (void)
64 { int i,j,k,ad,sad,nph,*svpt,*p1,*p2,**csvg,**csvh,nexp;
65 char pthere,w,ok,allok; long fac;
66 stop = &dummy;
67 if ((ip=fopen(inf1,"r"))==0)
68 { fprintf(stderr,"Cannot open %s.\n",inf1); return(-1); }
69 fscanf(ip,"%d%d%d%d",&npt,&npg,&nb,&k);
70 if (npt>mpt) {fprintf(stderr,"npt too big. Increase mpt.\n"); return(-1); }
71 if (nb>=mb) {fprintf(stderr,"nb too big. Increase MB.\n"); return(-1); }
72 if (2*nb*npt>svsp)
73 { fprintf(stderr,"svsp too small. Increase SVSP.\n"); return(-1); }
74 if (k<=0) {fprintf(stderr,"Wrong input format for G.\n"); return(-1); }
75 for (i=1;i<=nb;i++) fscanf(ip,"%d",base+i);
76 for (i=1;i<=nb;i++) fscanf(ip,"%d",lorbg+i);
77 pst=perm-1; pend=perm+psp-1;
78 if (pst+2*npg*npt>pend)
79 { fprintf(stderr,"Out of space.Increase PSP.\n"); return(-1);}
80 for (i=1;i<=npg;i++)
81 { p1=pst; p2=pst+npt; pst+=(2*npt); rdperm(p1,p2); cp[2*i-1]=p2; }
82 for (i=1;i<=nb;i++) {svgptr[i]=sv+npt*(i-1)-1; rdsv(svgptr[i]); }
83 fclose(ip);
84 /* G is read. Now read the subgroup H */
85 if (triv)
86 { nph=0; for (i=1;i<=nb;i++)
87 { lorbh[i]=1; svhptr[i]=sv+npt*(nb+i-1)-1;
88 for (j=1;j<=npt;j++) svhptr[i][j]=0; svhptr[i][base[i]]=stop;
89 }
90 }
91 else
92 { if ((ip=fopen(inf2,"r"))==0)
93 { fprintf(stderr,"Cannot open %s.\n",inf2); return(-1); }
94 fscanf(ip,"%d%d%d%d",&i,&nph,&j,&k); w=0;
95 if (k<=0) { fprintf(stderr,"Wrong input format for H.\n"); return(-1); }
96 if (i!=npt || j!=nb) w=1;
97 else for (i=1;i<=nb;i++)
98 { fscanf(ip,"%d",lorbh+i); if (lorbh[i]!=base[i]) w=1; }
99 if (w)
100 { fprintf(stderr,"npt or base for G and H do not agree.\n");return(-1);}
101 for (i=1;i<=nb;i++) fscanf(ip,"%d",lorbh+i);
102 if (pst+2*npt*nph+2>pend)
103 { fprintf(stderr,"Out of space.Increase PSP.\n"); return(-1);}
104 for (i=1;i<=nph;i++)
105 { p2=pend-npt; p1=p2-npt; pend-=(2*npt); rdperm(p1,p2); cp[2*i-1]=p2;}
106 for (i=1;i<=nb;i++) {svhptr[i]=sv+npt*(nb+i-1)-1; rdsv(svhptr[i]); }
107 fclose(ip);
108 }
109 for (i=1;i<=nb;i++) fpt[i]=0;
110 p1=pst; p2=p1+npt; pst+=(2*npt); lcp=cp; expp= cp+2*npt-1;
111 cind=1; fac=1; coh_index[nb+1]=1; lnt=0;
112
113 /* Now we compute the indices of H in G in the stabilizer chain into the array
114 coh_index. The nontrivial indices are linked by fpt and bpt.
115 */
116 for (i=nb;i>=1;i--)
117 { j=lorbg[i]; fac=fac*j; j=lorbh[i]; fac=fac/j; ind=fac; coh_index[i]=ind;
118 if (ind>coh_index[i+1])
119 { if (lnt==0) lnt=i;
120 else { fpt[i]=fnt; bpt[fnt]=i; }
121 fnt=i;
122 }
123 }
124 bpt[fnt]= -1; fpt[lnt]= -1;
125 printf("Indices: ");
126 for (i=1;i<=nb;i++) printf("%5d",coh_index[i]); printf("\n");
127 if (ind==1) { fprintf(stderr,"G=H.\n"); return(-1); }
128 nexp=0; ind=1; oind=1; trptr=0;
129
130 /* Now we start the backtrack search for the left coset reps of H in G. The
131 necessary information about these is stored in the array tree. This is done
132 so as to make it as easy as possible to compute which coset an arbitrary
133 element of G lies in. If expg is true, then those perms in G which arise
134 are expanded and stored.
135 If dcr is true, then a small subset of the coset reps will need to be
136 output as perms, later. For this purpose, the base images of all coset reps
137 are output to a temporary file.
138 ind is the current coh_index coh_index[bno], and cind is the no of coset reps
139 found so far.
140 */
141 for (i=fnt;i!= -1;i=fpt[i])
142 { trad[i]=tree+trptr; tree[trptr]=base[i]; trptr+=2; }
143 if (dcr) opx=fopen(temp2,"w");
144 if (cr)
145 { op=fopen(outf3,"w");
146 fprintf(op,"%4d %d%4d%4d\n",npt,coh_index[1]-1,0,0);
147 }
148 sad=trptr; tree[trptr]=1; trptr++;
149 for (bno=lnt;bno!= -1;bno=bpt[bno])
150 { printf("bno=%d.\n",bno);
151 *gorb=1; gorb[1]=base[bno];
152 ind=coh_index[bno]; csvg=svgptr[bno]; csvh=svhptr[bno];
153 allok=(lorbh[bno]==1);
154 /* if allok, then all coset reps in this link of the stabilizer chain will be
155 coset reps of H in G.
156 */
157 if (expg) {for (i=1;i<=npt;i++) expp[i]=0; expp[base[bno]]=stop; }
158 for (pt=1;pt<=npt;pt++)
159 { svpt=csvg[pt];
160 if (svpt!=0 && svpt!=stop)
161 { ucp=lcp-1; addsvf(pt,csvg); pthere=0;
162 if (expg)
163 { for (i=1;i<=npt;i++) p1[i]=image(i);
164 invert(p1,p2); ucp=lcp; *ucp=p1;
165 }
166 j=sad;
167 for (i=fpt[bno];i!= -1;i=fpt[i])
168 { tailad[i]=tree+j; cpad[i]=ucp+1; j+=2; }
169 /* In the search, for each basic coset rep g in the stabilizer chain, we have
170 to try out gh as a possible coset rep of H in G, for each h in the list of
171 coset reps of H[bno-1] in G[bno-1]. coind and oind record the search
172 through this list.
173 */
174 coind=1;
175 while (coind<=oind)
176 { if (coind>1) advance(); ok=1;
177 if (allok==0)
178 for (i=1;i<=*gorb && ok;i++) if (csvh[image(gorb[i])]!=0) ok=0;
179 /* That was the membership test. ok true means we have found a new rep */
180 if (ok)
181 { if (pthere)
182 { ad=fpt[bno];
183 while (*tailad[ad]== *trad[ad]) ad=fpt[ad];
184 tree[trptr]= *tailad[ad];
185 }
186 else
187 { pthere=1; ad=bno; tree[trptr]=pt;
188 if (expg)
189 { expp[pt]=p1; p1=pst; p2=p1+npt; pst+=(2*npt); nexp++;
190 if (pst>pend)
191 { fprintf(stderr,"Out of perm space. Increase PSP.\n");
192 return(-1);
193 }
194 }
195 }
196 *(trad[ad]+1)=trptr; trad[ad]=tree+trptr; trptr+=2;
197 for (i=fpt[ad];i!= -1;i=fpt[i])
198 { *(trad[i]+1)= -1; tree[trptr]= *tailad[i];
199 trad[i]=tree+trptr; trptr+=2;
200 }
201 cind++; tree[trptr]=cind; trptr++;
202 if (trptr>trsp)
203 { fprintf(stderr,"Out of tree space. Increase TRSP.\n");
204 return(-1);
205 }
206 if (dcr) for (i=fnt;i!= -1;i=fpt[i]) fprintf(opx,"%5d",*trad[i]);
207 /* if cr, then we output each coset rep */
208 if (cr)
209 if (npt>=1000)
210 { for (i=1;i<=npt;i++) fprintf(op,"%5d",image(i));fprintf(op,"\n");}
211 else
212 { for (i=1;i<=npt;i++) fprintf(op,"%4d",image(i));fprintf(op,"\n");}
213 if (cind==ind) goto nextbno;
214 } /* if (ok) */
215 coind++;
216 } /* while (coind<=oind) */
217 (*gorb)++; gorb[*gorb]=pt;
218 } /* if (svpt!=0 ... */
219 } /* for (pt=1;pt<=npt;... */
220 fprintf(stderr,"Premature end of search.\n"); return(-1);
221 nextbno:
222 if (expg) for (i=1;i<=npt;i++) svgptr[bno][i]=expp[i];
223 sad-=2; oind=ind;
224 } /* main bno loop */
225 if (dcr) fclose(opx); if (cr) fclose(op);
226 for (i=fnt;i!= -1;i=fpt[i]) *(trad[i]+1)= -1;
227 printf("Tree space used = %d.\n",trptr);
228 if (expg) printf("Number of perms expanded = %d.\n",nexp);
229 return(0);
230 }
231
232 int
advance(void)233 advance (void)
234 /* This advances the element h in the search through elements gh */
235 { int ad,k,*p,*q;
236 ad=lnt;
237 while ((k= *(tailad[ad]+1)) == -1) ad=bpt[ad];
238 q=tree+k; ucp=cpad[ad]-1;
239 while (ad!= -1)
240 { tailad[ad]=q; cpad[ad]=ucp+1;
241 if (expg) { p=svgptr[ad][*q]; if (p!=stop) {ucp++; *ucp=p; }}
242 else addsvf(*q,svgptr[ad]);
243 q+=2; ad=fpt[ad];
244 }
245 }
246