1 #include "defs.h"
2
3 # define MSP 100000
4 # define MM 101
5 # define MV 1000
6 # define MPT 2000
7 # define MPR 2000
8 # define MDIM 300
9 # define MARG 3000
10
11 int msp=MSP;
12 short prime,dim,*spv,**spm,nmat,ordim,npt,*fptr,*nfptr,**sp,**imsp,*endsp,
13 mspace[MSP],*vec[MV],**mat[MM],pinv[MPR],*pptr[MM],*spptr[MPT+1],
14 mm=MM,mv=MV,mpt=MPT,mpr=MPR,mdim=MDIM,marg=MARG;
15 char orvec;
16 FILE *ip,*op;
17
18 int
main(int argc,char * argv[])19 main (int argc, char *argv[])
20 { short nm,nv,i,j,pt,perm,*p1,*p2,*q1,arg;
21 char inf[80],outf[80],err; int x;
22 /* Defaults: inf=gpname.inmat outf=gpname.inperm */
23 arg=1; err=0;
24 if (argc<=arg) { err=1; goto error;}
25 strcpy(inf,argv[1]); strcat(inf,"."); strcpy(outf,inf);
26 arg++; if (argc<=arg) strcat(inf,"inmat"); else strcat(inf,argv[arg]);
27 arg++; if (argc<=arg) strcat(outf,"inperm"); else strcat(outf,argv[arg]);
28 if ((ip=fopen(inf,"r"))==0)
29 { fprintf(stderr,"Cannot open %s.\n",inf); exit(1); }
30 fscanf(ip,"%hd%hd%hd",&prime,&dim,&nmat);
31 if (prime>mpr) {fprintf(stderr,"prime too big. Increase MPR.\n"); exit(1);}
32 if (dim>mdim) {fprintf(stderr,"dim too big. Increase MDIM.\n"); exit(1);}
33 if (nmat>mm) {fprintf(stderr,"Too many mats. Increase MM.\n"); exit(1); }
34 setpinv();
35 nm=nmat+1; nv=(nm+2)*dim;
36 if (nv>=mv) {fprintf(stderr,"Too many vectors. Increase MV.\n"); exit(1);}
37 x=(nv+1)*dim;
38 if (x>msp-marg) {fprintf(stderr,"Out of space. Increase MSP.\n"); exit(1);}
39 /* Set matrix pointers as in mcp.c */
40 for (i=0;i<=nv;i++) vec[i]=mspace-1+i*dim;
41 for (i=0;i<=nmat;i++) { mat[i]=vec-1+i*dim; if (i>0) readmat(mat[i]); }
42 fclose(ip);
43 endsp=mspace+msp-marg;
44 spm=mat[0]; spv=spm[1]; sp=vec-1+nm*dim; imsp=sp+dim; fptr=vec[nv];
45 for (i=1;i<=nmat;i++) {pptr[i]=fptr; fptr+=MPT; }
46 if (fptr>=endsp) {fprintf(stderr,"Out of space. Increase MSP.\n"); exit(1);}
47 reenter:
48 printf("Input dim of orbit space (-1 for orbit of vector): ");
49 scanf("%hd",&ordim);
50 if (ordim< -1 || ordim==0 || ordim>dim)
51 { printf("Inappropriate dimension.\n"); while (getchar()!='\n');goto reenter;}
52 else if (ordim== -1) { ordim=1; orvec=1; } else orvec=0;
53 printf("Input generating vectors of orbit space.\n");
54 for (i=1;i<=ordim;i++)
55 { p1=imsp[i]; for (j=1;j<=dim;j++) scanf("%hd",p1+j); }
56 npt=1;
57 if (orvec==0) normalize();
58 encode(); spptr[1]=fptr; fptr=nfptr;
59 if (fptr>=endsp) {fprintf(stderr,"Out of space. Increase MSP.\n"); exit(1);}
60
61 for (pt=1;pt<=npt;pt++)
62 { decode(pt);
63 for (perm=1;perm<=nmat;perm++)
64 { for (j=1;j<=ordim;j++) im(sp[j],imsp[j],mat[perm]);
65 if (orvec==0) normalize(); encode();
66 for (i=1;i<=npt;i++)
67 { p1=fptr; p2=nfptr; q1=spptr[i]+1;
68 while (++p1<=p2) { if (*p1!= *q1) break; q1++; }
69 if (p1>p2) break;
70 }
71 pptr[perm][pt]=i;
72 if (i>npt)
73 { npt++; if (npt>mpt) {fprintf(stderr,"Too many points.\n"); exit(1); }
74 spptr[npt]=fptr; fptr=nfptr;
75 if (fptr>=endsp)
76 {fprintf(stderr,"Out of space. Increase MSP.\n"); exit(1);}
77 if (npt%10==0) printf("npt=%d.\n",npt);
78 }
79 }
80 }
81 op=fopen(outf,"w"); fprintf(op,"%4d%4d%4d%4d\n",npt,nmat,0,0);
82 for (i=1;i<=nmat;i++)
83 if (npt>=1000)
84 { for (j=1;j<=npt;j++) fprintf(op,"%5d",pptr[i][j]); fprintf(op,"\n"); }
85 else
86 { for (j=1;j<=npt;j++) fprintf(op,"%4d",pptr[i][j]); fprintf(op,"\n"); }
87 error:
88 if (err)
89 { fprintf(stderr,"Usage: matperm gpname [inf] [outf]\n");
90 exit(1);
91 }
92 exit(0);
93 }
94
95 int
encode(void)96 encode (void)
97 { short i,*p1,*p2,*q; char c;
98 nfptr=fptr;
99 for (i=1;i<=ordim;i++)
100 { c=orvec; p1=imsp[i]; p2=p1+dim; q=p1;
101 while (++q<=p2) if (*q!=0)
102 if (c==0) { c=1; *(++nfptr)=q-p1;}
103 else {*(++nfptr)=q-p1; *(++nfptr)= *q; }
104 *(++nfptr)=0;
105 }
106 }
107
108 int
decode(int n)109 decode (int n)
110 { short i,*p1,*p2,*ptr,*q; char c;
111 ptr=spptr[n];
112 for (i=1;i<=ordim;i++)
113 { c=orvec; p1=sp[i]; p2=p1+dim; q=p1; while (++q<=p2) *q=0;
114 while (*(++ptr)!=0)
115 if (c==0) {p1[*ptr]=1; c=1;} else {p1[*ptr]= *(ptr+1); ptr++; }
116 }
117 }
118
119 int
normalize(void)120 normalize (void)
121 { short rst,i,j,k,fac,*p1,*p2,*q1;
122 int n;
123 rst=1;
124 for (i=1;i<=dim;i++) for (j=rst;j<=ordim;j++) if (imsp[j][i] !=0)
125 { p1=imsp[j]; if (j>rst) {imsp[j]=imsp[rst]; imsp[rst]=p1;}
126 fac=pinv[p1[i]]; p2=p1+dim+1;
127 while (--p2>p1) {n= *p2*fac; n%=prime; *p2=n; }
128 for (k=1;k<=ordim;k++)
129 { if (imsp[k][i]==0 || k==rst) continue;
130 fac=prime-imsp[k][i]; p2=p1+dim+1; q1=imsp[k]+dim;
131 while (--p2>p1) { n= *q1+fac* *p2; n%=prime; *q1=n; q1--; }
132 }
133 if (rst==ordim) return(0);
134 rst++; break;
135 }
136 }
137
138 int
setpinv(void)139 setpinv (void)
140 { int i,j;
141 for (i=0;i<prime;i++) pinv[i]=0;
142 for (i=1;i<prime;i++) if (pinv[i]==0) for (j=1;j<prime;j++)
143 if (i*j % prime ==1) { pinv[i]=j; pinv[j]=i; break; }
144 }
145