1 #include "defs.h"
2 #include "permfns.h"
3 
4 extern char wrd,nt,isbase,inf[],outf1[],outf2[],fixed[];
5 extern int  perm[],sv[],cp[],actgen[],orb[],
6        base[],lorb[],order[],pno[], *pptr[],*svptr[],
7        mp,mb,mnpt;
8 extern int psp,svsp;
9 int npt;
10 FILE *ip,*op;
11 
12 /* The data structures in this program are similar to most permutation group
13    programs. Permutations are numbered 0,1,2,3,... (where 2x+1 is usually the
14    inverse of 2x) and stored in the array perm. Perm no x is pointed to by
15    pptr[x]. npt=no. of points. Usually pptr[i][npt+1] gives the number of the
16    group in the stabilizer chain in which the perm lies. So this no is i if the
17    perm fixes the first i-1 base points. nb=no of base points.
18    The base and the lengths of the orbits in the stab chain are stored in base
19    and lorb.
20    Schreier vectors are stored in sv, and pointed to by svptr[i], i=1,..,nb.
21    pno is a list of *pno (=pno[0]) perm nos.
22    cp is a similar list of length *cp, but it is used to represent the product
23    of the perms  cp[1]cp[2]..cp[*cp]. This product can be evaluated by the
24    procedure image in permfns.c
25 */
26 int
gpprog(void)27 gpprog (void)
28 { int i,j,k,l,m,n,nperms,nb,jobt,np2,seek,cord,ocord,given,ordknown,trivrel,
29         lsv,u,v,w,x,y,z,mxp,mnb,bpt,bno,id;
30   int quot;
31   float grporder;
32 
33   if ((ip=fopen(inf,"r"))== 0)
34   { fprintf(stderr,"Cannot open file %s\n",inf); return(-1); }
35   fscanf(ip,"%d%d%d%d",&npt,&nperms,&nb,&jobt);
36   if (npt>mnpt) { fprintf(stderr,"npt too big. Increase NPT.\n"); return(-1); }
37   if (jobt>0) { fprintf(stderr,"Wrong input format.\n"); return(-1); }
38   quot=psp/(npt+1); if (quot>mp) quot=mp; mxp=quot;
39   quot=svsp/npt; if (quot>mb) quot=mb; mnb=quot; if (mnb>=npt) mnb=npt-1;
40   if (nb>mnb)
41   { fprintf(stderr,"nb to big. Increase SVSP (or MB).\n"); return(-1); }
42 /* pptr[i] is the i th permutation, svptr[i] is the i the Schreier vector. */
43   for (i=0;i<mxp;i++) pptr[i]=perm+i*(npt+1)-1;
44   for (i=1;i<=mnb;i++) svptr[i]=sv+(i-1)*npt-1;
45 
46 /* Next we read any base points and orbit lengths */
47   if (nb!=0)
48   { for (i=1;i<=npt;i++) orb[i]=0;
49     for (i=1;i<=nb;i++)
50     { fscanf(ip,"%d",base+i); if (orb[base[i]]!=0)
51       { fprintf(stderr,"Repeated base element.\n"); return(-1); }
52       orb[base[i]]=1;
53     }
54   }
55   if (jobt<0)
56   { jobt= -jobt; if (jobt>nb) {fprintf(stderr,"jobt too big.\n"); return(-1); }
57     seek=jobt; given=jobt; ordknown=1;
58     for (i=1;i<=jobt;i++)  fscanf(ip,"%d",order+i);
59   }
60   else ordknown=0;
61 
62 /* Now we read the generating permutations */
63   np2=2*nperms-2;
64   for (i=0;i<=np2;i+=2)
65   { k=i/2 +1; j=readperm(pptr[i]);
66     if (j==2)
67     { fprintf(stderr,"Generator no %d is not a permutation.\n",k); return(-1); }
68     if (j==1)
69     { fprintf(stderr,"Generator no %d is the identity.\n",k);
70       if (nt) return(-1);
71       nperms-=1; np2-=2; i-=2;
72     }
73     else
74     { invert(pptr[i],pptr[i+1]); actgen[i]=1; x=1;
75       for (z=base[x];x<=nb && pptr[i][z]==z;z=base[x]) x++;
76       pptr[i][npt+1]=x;
77       if (x>nb)
78       { if (isbase)
79         { fprintf(stderr,"Given base is not a base!\n"); return(-1);}
80         nb++;  for (z=1;pptr[i][z]==z;z++); base[nb]=z;
81         printf("New base element no %d is %d.\n",nb,z);
82       }
83       if (x==1) printf("Generator no %d moves first base point.\n",k);
84       else printf("Generator no %d fixes first %d base point(s).\n",k,x-1);
85     }
86     if (nperms==0) { fprintf(stderr,"Trivial group!\n"); return(-1); }
87   }
88   fclose(ip);
89 
90   if (wrd) {op=fopen(outf2,"w"); fprintf(op,"%4d\n",nperms); }
91   bno=nb; for (i=0;i<=mnb;i++) lorb[i]=0;
92 
93 /* Now the main algorithm begins */
94 loop:
95   *pno=0; if (ordknown) ocord= (bno==nb) ? 1 : order[bno+1];
96 /* We make a list of the perm nos that fix the first bno-1 base pts */
97   for (i=0;i<=np2;i+=2)
98   { if (pptr[i][npt+1]>=bno && actgen[i]<=bno)
99     { (*pno)++; pno[*pno]=i; }
100   }
101   lorb[bno]=orbitsv(base[bno],svptr[bno],0);
102   if (ordknown)
103   { cord= (bno>=given) ? ocord*lorb[bno] : lorb[bno];
104     if (bno==seek && cord==order[bno])
105     { seek--; printf("Order is now as given for bno = %d.\n",bno);
106       bno--; if (bno==0) goto foundorder; goto loop;
107     }
108   }
109 /* Now we start generating the Schreier generators that fix the firat bno
110    base points, test them for membership, and add them as new generators
111    if necessary.
112 */
113   if (*pno!=0)
114   { nperms++; np2+=2; y=np2+1;
115     if (y>=mxp)
116     {fprintf(stderr,"Out of space. Increase PSP (or MP).\n"); return(-1); }
117     for (i=1;i<=npt;i++) fixed[i]=0;
118     for (i=1;i<bno;i++)
119     { u=base[i]; fixed[u]=1; pptr[np2][u]=u; pptr[y][u]=u; }
120     for (i=1;i<=lorb[bno];i++)
121     { *cp=0; addsv(orb[i],svptr[bno]);
122       for (w=1,x= *cp;w<=x;w++,x--)
123       { if (w==x) cp[w]-=1; else {u=cp[w]; cp[w]=cp[x]-1; cp[x]=u-1;}}
124       lsv= *cp;
125       for (j=1;j<=*pno;j++)
126       { *cp=lsv; u=pno[j]+1;
127         trivrel = (*cp>0) ? cp[*cp]==(u-1) : 0;
128         if (trivrel==0)
129         { (*cp)++; cp[*cp]=u; id=0;
130           for (l=bno;l<=nb;l++) fixed[base[l]]=0;
131           for (l=bno;l<=nb;l++)
132           { v=base[l]; u=image(v);
133             if (svptr[l][u]==0) goto newgen;
134             addsv(u,svptr[l]); pptr[np2][v]=v; pptr[y][v]=v; fixed[v]=1;
135           }
136           l=nb+1; id=1;
137 newgen:   if (isbase==0)
138           for (k=1;k<=npt;k++)
139           { if (fixed[k]==0)
140             { u=image(k); pptr[np2][k]=u; pptr[y][u]=k;
141               if (id && k!=u)
142               { id=0; nb++;
143   if (nb>mnb) { fprintf(stderr,"nb to big. Increase SVSP (or MB).\n"); return(-1); }
144                 base[nb]=k;
145                 printf("New base point no %d is %d.\n",nb,k);
146               }
147             }
148           }
149           if (id==0)
150           { pptr[np2][npt+1]=l; actgen[np2]=bno+1;
151             if (isbase) for (k=1;k<=npt;k++)
152             { u=image(k); pptr[np2][k]=u; pptr[y][u]=k;}
153             printf("New generator no %d, fixing first %d base pts is:\n",
154                     nperms,l-1);
155             for (m=1;m<=*cp;m++)
156             { z=cp[m]; y=z/2; x= (z==2*y) ? y+1 : -(y+1);
157               printf("%3d",x); if (m== *cp) printf("\n"); else printf("*");
158             }
159             if (wrd)
160             {for (m=0;m<=*cp;m++) fprintf(op,"%4d",cp[m]);fprintf(op,"\n");}
161             bno=l; goto loop;
162           }
163         }
164       }
165     }
166     nperms--; np2-=2;
167   }
168   if (ordknown && bno>given) order[bno]=cord;
169   bno--;
170 
171 foundorder:
172   if (bno==0)
173   { printf("The order of the group is:\n");
174     for (i=1;i<=nb;i++)
175     { printf("%3d",lorb[i]);
176       if (i==nb) printf(" =\n"); else printf("*");
177     }
178     grporder=1; for (i=1;i<=nb;i++) grporder *= lorb[i];
179     printf("%g\n",grporder);
180   }
181   else goto loop;
182   if (wrd) { fprintf(op,"%d\n",-1); fclose(op); }
183 
184 /* Now we output the generating perms and Schreier vectors */
185   op=fopen(outf1,"w");
186   fprintf(op,"%4d%4d%4d%4d\n",npt,nperms,nb,1);
187   printbaselo(nb,base,lorb);  *pno=0;
188   for (i=1;i<=nperms;i++) {(*pno)++; pno[*pno]=2*(i-1);}
189   printpsv(nb,pno,svptr);
190   fclose(op); return(0);
191 }
192