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