1 #include "defs.h"
2 #include "permfns.h"
3
4 extern char inf[],outf[];
5 extern short mp,mxexp,mb,mnpt,
6 perm[],svp2[],svg2[],cp[],orb[],
7 expperm[],adpt[],expcp[],lorbg[],lorbp[],lorbdef[],
8 base[],ntfpt[],ntbpt[],start[],invbase[],pno[],
9 facord[],orno[],lporb[],deftime[],orbperm[],
10 genorb[],*pptr[],*svgptr[],*svpptr[],*expptr[];
11 extern int psp,expsp,svsp;
12
13 short npt,p,im,cb,adno,*tp;
14 char ok;
15 FILE *ip,*op;
16 /* This program and normrun involve backtrack searches. These can be speeded up
17 by storing the perms for each coset rep in the stabilizer chain, which can
18 be computed using Schreier vectors. As many such perms as space allows are
19 stored in this way, in array expperm, using perm pointer expptr, and expcp
20 corresponds to cp.
21 */
22
23 int
sylprog(void)24 sylprog (void)
25 { char nontriv,bt,seek,b,incadno;
26 short i,j,k,l,m,n,lnt,fnt,mxp,mnb,mexp,nperms,nb,ct,np2,bno,stp,
27 lexp,*z,*itp,*ap,*sva;
28 int quot;
29 if ((ip=fopen(inf,"r"))==0)
30 { fprintf(stderr,"Cannot open %s.\n",inf); return(-1);}
31 printf("Input prime! "); scanf("%hd",&p);
32 fscanf(ip,"%hd%hd%hd%hd",&npt,&nperms,&nb,&l);
33 if (npt>mnpt) {fprintf(stderr,"npt too big. Increase NPT.\n"); return(-1);}
34 if (nb>mb) {fprintf(stderr,"nb too big. Increase MB.\n"); return(-1);}
35 if (nb*npt>svsp)
36 { fprintf(stderr,"Not enough sv space.Increase SVSP.\n");return(-1);}
37 if (l<=0) { fprintf(stderr,"Wrong input format.\n"); return(-1); }
38 readbaselo(nb,base,lorbg);
39 for (i=1;i<=npt;i++) invbase[i]=0;
40 for (i=1;i<=nb;i++) invbase[base[i]]=i;
41 lnt=0; fnt=0;
42 /* Determine order of sylp subgroup . Links ntfpt and ntbpt are used to bypass
43 trivial indices in the stab chain.
44 */
45 for (i=nb;i>=1;i--)
46 { j=lorbg[i]; lorbp[i]=1;
47 if (j>1)
48 { if (lnt==0) { lnt=i; k=i; ntfpt[lnt]=nb+1; }
49 else { ntbpt[k]=i; ntfpt[i]=k; k=i; }
50 while (j%p==0) {lorbp[i]*=p; j/=p; }
51 if (lorbp[i]>1) fnt=i;
52 }
53 else invbase[base[i]]=0;
54 }
55 ntbpt[fnt]=0;
56 if (fnt==0) {fprintf(stderr,"Sylp-group is trivial!\n"); return(-1); }
57 printf("Orbit lengths of Sylp-group:\n");
58 for (i=1;i<=nb;i++) printf("%4d",lorbp[i]); printf(".\n");
59
60 quot=psp/(npt+1); if (quot>mp) quot=mp; mxp=quot;
61 quot=expsp/npt; if (quot>mxexp) quot=mxexp; mexp=quot;
62 np2=2*nperms-1;
63 if (np2>=mxp)
64 {fprintf(stderr,"Out of space. Increase PSP (or MP).\n"); return(-1); }
65 for (i=0;i<mxp;i++) pptr[i]=perm+(npt+1)*i-1;
66 for (i=1;i<=nb;i++) svgptr[i]=svg2+npt*(i-1)-1;
67 for (i=1;i<=nb;i++) svpptr[i]=svp2+npt*(i-1)-1;
68 readpsv(0,nb,nperms,svgptr);
69 fclose(ip);
70
71 /* We now determine how much space we have to store cosetrep perms, and
72 calculate these with the function exprep
73 */
74 lexp=0; ct=0;
75 for (i=nb;i>fnt;i--)
76 { ct+=(lorbg[i]-1); if (ct+i>=mexp) {lexp=i; break; } }
77 if (lexp==nb)
78 {fprintf(stderr,"Not enough expanding space. Increase MXEXP.\n");return(-1);}
79 if (lexp==0) lexp=fnt;
80 z=expperm-1; i=0;
81 while (i!=lexp)
82 { i= (i==0) ? fnt : ntfpt[i]; expptr[i]=z;
83 z+=npt; start[i]=i-1;
84 }
85 printf("lexp=%d.\n",lexp);
86 start[lexp+1]=lexp; k=lexp;
87 for (i=lexp+1;i<=lnt;i++)
88 { l=lorbg[i]; start[i+1]=start[i]+l-1;
89 if (l>1) for (j=1;j<=npt;j++) if (svgptr[i][j]>0)
90 { k++; expptr[k]=z; z+=npt; exprep(j,k,svgptr[i]); }
91 }
92 nontriv=0; *pno=0; deftime[0]=0; stp=np2+1;
93 for (i=1;i<=nb;i++)
94 { for (j=1;j<=npt;j++) {svpptr[i][j]=0; svpptr[i][base[i]]= -1; }
95 lorbdef[i]=1;
96 }
97
98 /* Starting search. lorbdef records the length of the orbit ofthe Sylp group
99 found so far. tp is a perm that we are currently considering as a
100 possible new element of Sylp. nontriv is set true as soon as the first
101 element of Sylp has been found. adno is then number in the stabilizer chain
102 for which the coset reps are currently being advanced in the search.
103 */
104 for (bno=lnt;bno>=fnt;bno=ntbpt[bno])
105 { printf("bno=%d.\n",bno);
106 orb[1]=base[bno]; adno=bno;
107 while (lorbdef[bno]<lorbp[bno])
108 { *expcp=0; tp=pptr[np2+1]; itp=pptr[np2+2];
109 if (np2+2>=mxp)
110 { fprintf(stderr,"Out of space. Increase PSP.\n"); return(-1); }
111 for (i=1;i<=npt;i++) tp[i]=i; tp[npt+1]=bno;
112
113 /* If nontriv, then we compute all orbits of Sylp. The new element must
114 permute these orbits.
115 */
116 if (nontriv)
117 { allorbs(lporb,orno); bt=0;
118 for (i=1;i<= *lporb;i++) {orbperm[i]=0; deftime[i]=0; }
119 }
120 seek=1;
121 while (seek)
122 { for (i=1;i<=npt;i++) itp[i]=0;
123 b=1; ok=1;
124 while(1)
125 { k=expcp[*expcp];
126 /* Advance to next element in search. ok is set false if this element is
127 eliminated from membership of Sylp, and we advance again.
128 */
129 if (b && (*expcp==0 || k<=start[adno]))
130 { (*expcp)++;
131 if (adno<=lexp) {adpt[adno]=0; expcp[*expcp]=adno;}
132 else expcp[*expcp]=start[adno];
133 b=0;
134 }
135 else
136 { if (adno<=lexp)
137 { sva=svgptr[adno]; ap=adpt+adno; (*ap)++;
138 while (sva[*ap]<=0 && *ap<=npt) (*ap)++;
139 if (*ap<=npt) { exprep(*ap,adno,sva); break; }
140 }
141 else if (k<start[adno+1]) {expcp[*expcp]++; break;}
142 if (*expcp==1)
143 { fprintf(stderr,
144 "Premature end of search. Probably %d is not prime.\n",p);
145 return(-1);
146 }
147 (*expcp)--; bt=1; adno=ntbpt[adno]; b=1;
148 if (adno==bno && nontriv)
149 { im=expimage(base[adno]); l=orno[im];
150 if (lporb[l]>1) for (i=1;i<=npt;i++)
151 if (orno[i]==l) svpptr[adno][i]= -2;
152 /* Setting sv= -2 when bno=adno avoids testing more elements than necessary
153 which map the main orbit of base[bno] to another.
154 */
155 }
156 }
157 }
158 incadno=1;
159 while (incadno)
160 { cb=base[adno]; im=expimage(cb);
161 if (adno==bno && svpptr[adno][im]!=0) {ok=0;break;}
162 else
163 { j=im; k=invbase[j]; ct=1;
164 while (k>0 && k<adno)
165 { j=tp[j]; k=invbase[j]; ct++; }
166 if (k==adno && ppower(ct)==0) {ok=0;break;}
167 /* Because this element cannot be a p-element */
168 }
169 if (nontriv)
170 { if (bt) for (i=1;i<=*lporb;i++)
171 { j=orbperm[i]; if (deftime[j]>=adno) {orbperm[i]=0;deftime[j]=0;}}
172 bt=0; deforbperm();
173 if (ok==0) break;
174 /* This means that the element does not permute the orbits */
175 }
176 else tp[cb]=im;
177 if (adno==lnt)
178 { incadno=0;
179 for (cb=1;cb<=npt;cb++) if (invbase[cb]==0)
180 { im=expimage(cb);
181 if (nontriv)
182 { deforbperm(); if (ok==0) { bt=1; break; } }
183 else tp[cb]=im;
184 }
185 }
186 else adno=ntfpt[adno];
187 }
188 if (ok) for (i=1;i<=npt;i++) if (itp[i]==0)
189 { j=i; l=1; k=tp[j]; itp[k]=j;
190 while (k!=i) { l++; j=k;k=tp[j]; itp[k]=j; }
191 if (ppower(l)==0) { ok=0;bt=1;break;}
192 }
193 /* if ok then tp is a p-element permuting the orbits. The final test is to
194 check that it normalizes Sylp.
195 */
196 if (ok && nontriv) for (i=stp;i<np2 && ok;i+=2)
197 { *cp=0;
198 for (j=fnt;j<=lnt;j=ntfpt[j])
199 { k=image(tp[pptr[i][itp[base[j]]]]);
200 if (svpptr[j][k]!=0) addsv(k,svpptr[j]);
201 else {ok=0; bt=1; break; }
202 }
203 }
204 /* if ok then tp is the new element of Sylp */
205 if (ok)
206 { (*pno)++; pno[*pno]=np2+1; k=lorbdef[bno];
207 l=orbitsv(base[bno],svpptr[bno],k);
208 facord[np2+1]=l/k; lorbdef[bno]=l;
209 np2+=2; adno=bno; nontriv=1; seek=0;
210 printf("Found element in Sylp-group. Lorbdef=%d.\n",l);
211 for (i=1;i<=npt;i++) if (svpptr[bno][i]== -2) svpptr[bno][i]=0;
212 }
213 }
214 }
215 for (i=1;i<=npt;i++) if (svpptr[bno][i]== -2) svpptr[bno][i]=0;
216 }
217 op=fopen(outf,"w");
218 fprintf(op,"%4d%4d%4d%4d\n",npt,*pno,nb,4);
219 printbaselo(nb,base,lorbp); printpsv(nb,pno,svpptr);
220 fprintf(op,"%3d",p);
221 for (i=stp;i<np2;i+=2) fprintf(op,"%4d",facord[i]); fprintf(op,"\n");
222 return(0);
223 }
224
225 int
ppower(int x)226 ppower (int x)
227 /* This tests whether x is a power of p */
228 { while (x!=1)
229 { if (x%p==0) x /=p; else return(0); }
230 return(1);
231 }
232
233 int
deforbperm(void)234 deforbperm (void)
235 /* This tests whether the fact that tp sends cb (=base[adno]) to im
236 contradicts permutation of the orbits. If so then ok is set false.
237 */
238 { short i,j,k; i=orno[cb]; j=orno[im]; k=orbperm[i];
239 if (k==0)
240 { if (deftime[j]==0) if (lporb[i]==lporb[j])
241 { orbperm[i]=j; deftime[j]=adno; }
242 else ok=0; else ok=0;
243 }
244 else if (k!=j) ok=0;
245 if (ok) tp[cb]=im;
246 return(0);
247 }
248