1 #include "defs.h"
2
3 extern char inf[],outf[],temp[],temp2[],mspace[],*vec[],**mat[],
4 cvec[],pinv[],mstr[],full,intop,opt,aut;
5 extern short svec[],mdim,mv,mm,mpr,msl;
6 extern int space;
7 char prime;
8 short dim,maxnull;
9
10 FILE *ip,*op;
11
12 int
pkprog(void)13 pkprog (void)
14 { char ngens,neg,fac,there,saved,gottm,gotsgn,err,digexp,nc,provn,
15 trd,gotfac,heldc,len,**m0,**m1,*m01,irr;
16 short nvecs,i,j,rno1,rno2,nlty,fnz,adno,sdim,qdim,tdim;
17 int msp;
18 irr= opt ? 1 : 0;
19 if ((ip=fopen(inf,"r"))==0)
20 { fprintf(stderr,"Cannot open %s.\n",inf); return(-1);}
21 fscanf(ip,"%hd%hd%hd",&rno1,&dim,&rno2);
22 prime=rno1; ngens=rno2;
23 if (prime>mpr)
24 { fprintf(stderr,"prime too big. Nothing to be done.\n"); return(-1);}
25 if (dim>mdim)
26 { fprintf(stderr,"dim too big. Increase MDIM.\n"); return(-1);}
27 if (ngens>=mm)
28 { fprintf(stderr,"Too many generators. Nothing to be done.\n"); return(-1);}
29 nvecs=dim*(ngens+1);
30 if (nvecs>mv)
31 { fprintf(stderr,"Too many vectors. Increase MV.\n"); return(-1);}
32 msp=nvecs*dim;
33 if (msp>space)
34 { fprintf(stderr,"Not enough space. Increase SPACE.\n"); return(-1);}
35 seeknln(); setpinv();
36 /* Set vector and matrix pointers, as usual.
37 We take care only to use ngens+1 matrices, thus saving as much storage as
38 possible. This means using temporary files to store matrices occasionally.
39 */
40 for (i=0;i<nvecs;i++) vec[i]= mspace-1+i*dim;
41 for (i=0;i<=ngens;i++) mat[i]=vec-1+i*dim;
42 for (i=1;i<=ngens;i++) readmat(mat[i]);
43 fclose(ip);
44 if (aut)
45 { printf("Enter maximum value of nullity that should be handled: ");
46 scanf("%hd",&maxnull);
47 }
48 m0=mat[0]; m1=mat[1]; m01=m0[1];
49 /* The next section reads in the given expression in the generators, and
50 computes the corresponding matrix as m0. A temporary file is used to
51 remember the sum so far, whilst a product is being computed.
52 */
53 while (1)
54 { printf("Enter element of group algebra. End expression with $. Example:\n");
55 printf("Enter 3a[3]a[3]a[2] + a[2] - a[4] as 3.3*3*2 + 2 - 4$\n");
56 fflush(stdout);
57 there=0; saved=0; heldc=0; err=0;
58 while (heldc!='$')
59 { gotfac=0; gotsgn=0; gottm=0; provn=0; digexp=0; len=0;
60 while (gottm==0)
61 { if (heldc) { nc=heldc; heldc=0;} else nc=findc();
62 if (gotsgn==0)
63 { if (nc=='-') { neg=1; gotsgn=1; digexp=1; continue;}
64 else if (nc=='+') { neg=0; gotsgn=1; digexp=1; continue;}
65 else if (there==0) { neg=0; gotsgn=1; digexp=1;}
66 }
67 if ((nc=='+' || nc=='-' || nc=='$') && digexp==0 && (len!=0 || provn!=0))
68 { heldc=nc; gottm=1;
69 if (gotfac==0)
70 { if (provn)
71 { fac= neg ? prime-1 : 1; len=1; mstr[1]=provn-'0';}
72 }
73 }
74 else if (nc=='.' && digexp==0 && provn!=0 && gotfac==0)
75 { gotfac=1; fac=provn-'0'; fac%=prime; provn=0;
76 if (neg) fac=prime-fac; digexp=1;
77 }
78 else if (nc=='*' && digexp==0 && (provn!=0 || len!=0))
79 { if (gotfac==0)
80 { gotfac=1; fac= neg ? prime-1 : 1; len=1; mstr[1]=provn-'0';}
81 digexp=1;
82 }
83 else if (digit(nc) && digexp)
84 { if (gotfac==0) provn=nc;
85 else { len++; mstr[len]= nc-'0';}
86 digexp=0;
87 }
88 else
89 { fprintf(stderr,"Format error; Try again.\n"); err=1;
90 while (getchar()!='\n');
91 break;
92 }
93 }
94 if (err) break;
95 if (len==1)
96 { if (there) sum(fac,mat[mstr[1]],m0); else ncopy(fac,mat[mstr[1]],m0);}
97 else
98 { if (there)
99 /* We have a partial sum, and now another product to compute, so we must
100 store the partial sum in a temporary file.
101 */
102 { op=fopen(temp,"w"); printmat(m0); saved=1; fclose(op);}
103 ncopy(fac,mat[mstr[1]],m0);
104 for (i=2;i<=len;i++) prod(m0,mat[mstr[i]],m0);
105 if (saved)
106 { ip=fopen(temp,"r");
107 for (i=1;i<=dim;i++) rvecsum(1,m0[i]);
108 fclose(ip); saved=0;
109 }
110 }
111 there=1;
112 }
113 if (err) continue;
114 /* err means the expression input was garbled */
115 if (intop) for (i=1;i<=dim;i++)
116 { for (j=1;j<=dim;j++) printf("%3d",m0[i][j]); printf("\n"); }
117 /* Store the matrix m0 in a temporary file temp. Then compute its nullity.
118 (nlty writes generators of the null space to another temporary file, temp)
119 */
120 op=fopen(temp2,"w"); printmat(m0); fclose(op);
121 op=fopen(temp,"w"); nlty=null(m0,m1); fclose(op);
122 ip=fopen(inf,"r"); seeknln(); readmat(m1); fclose(ip);
123 if (nlty==0 || (aut && nlty>maxnull)) continue;
124 if (aut) break;
125 printf("Shall we use this matrix (y/n)? ");
126 nc=findc();
127 if (nc=='y') break;
128 }
129 /* Now we have chosen the matrix that we are going to use for the irreducibility
130 test. Next we arrange to consider the one-dimensional subspaces of the
131 null-space in order. cvec is used for this purpose.
132 */
133 for (i=1;i<=nlty;i++) cvec[i]=0;
134 fnz=nlty; adno=nlty; trd=0;
135 while (fnz!=0)
136 { cvec[adno]++;
137 for (i=1;i<=dim;i++) m01[i]=0;
138 /* Compute the current vector in the nullspace by reading from the file temp,
139 and store it in m01
140 */
141 ip=fopen(temp,"r");
142 for (i=1;i<=nlty;i++)
143 { if ((fac=cvec[i])==0) seeknln(); else {rvecsum(fac,m01); seeknln();}}
144 fclose(ip);
145 if (intop)
146 { for (j=1;j<=dim;j++) printf("%3d",m01[j]); printf("\n");
147 if (opt)
148 { printf("Shall we skip this one? "); nc=findc();
149 if (nc=='y') goto next;
150 }
151 }
152 /* Now compute the subspace it generates */
153 sdim=spgen(m0,ngens);
154 if (intop)
155 { for (i=1;i<=dim;i++) printf("%3d",svec[i]); printf("\n");
156 for (i=1;i<=dim;i++)
157 { for (j=1;j<=dim;j++) printf("%3d",m0[i][j]); printf("\n");}
158 }
159 if (sdim<dim)
160 /* sdim<dim means the module is reducible */
161 if (opt==0) break;
162 else
163 { printf("Shall we compute (y/n)? ");
164 nc=findc();
165 if (nc=='y') { opt=0; break;}
166 irr=0;
167 }
168 next:
169 /* The last vector generated the whole space, so go on to the next */
170 adno=nlty;
171 while (cvec[adno]==prime-1) { cvec[adno]=0; adno--;}
172 if (adno<=fnz) { if (adno==fnz) {cvec[adno]=0; adno--;} fnz--;}
173 }
174 if (opt || sdim==dim)
175 /* Each vector in the null-space generated the whole space. Finally we have to
176 transpose everything, and try just one vector in the null-space.
177 */
178 { printf("Transforming.\n"); fflush(stdout); trd=1;
179 ip=fopen(temp2,"r"); readmat(m0); fclose(ip);
180 trans(m0,m1); op=fopen(temp,"w"); null1(m1,m0); fclose(op);
181 /* That ruined m1 (the first generator), so we must re-read it */
182 ip=fopen(inf,"r"); seeknln(); readmat(m1); fclose(ip);
183 for (i=1;i<=ngens;i++) {trans(mat[i],m0); copy(m0,mat[i]);}
184 ip=fopen(temp,"r");
185 for (i=1;i<=dim;i++) {fscanf(ip,"%hd",&rno1); m01[i]=rno1;}
186 fclose(ip);
187 sdim=spgen(m0,ngens);
188 if (opt && sdim<dim)
189 { printf("Shall we compute (y/n)? ");
190 nc=findc();
191 if (nc=='y') opt=0; else irr=0;
192 }
193 if (intop)
194 { for (i=1;i<=dim;i++) printf("%3d",svec[i]); printf("\n");
195 for (i=1;i<=dim;i++)
196 { for (j=1;j<=dim;j++) printf("%3d",m0[i][j]); printf("\n");}
197 }
198 }
199
200 if (opt && irr==0) { unlink(temp); unlink(temp2); return(0);}
201 if (sdim==dim)
202 { printf("Module is irreducible.\n");
203 unlink(temp); unlink(temp2);return(0);
204 }
205 /* Space is reducible. Finally, we compute the generators of the subspace
206 and quotient space
207 */
208 qdim=dim-sdim;
209 strcpy(outf,inf);
210 if (full) { strcat(outf,"f"); tdim=dim;}
211 else { strcat(outf,"s"); tdim=sdim;}
212 /* if we have transposed, then we first output to the temporary file,
213 and then later read it back in and transpose it back again.
214 */
215 if (trd) op=fopen(temp,"w");
216 else { op=fopen(outf,"w"); fprintf(op,"%3d%5d%3d\n",prime,tdim,ngens);}
217 opnmat(m0,ngens,tdim,1);
218 if (trd==0|| full) fclose(op);
219 if (full==0)
220 { strcpy(outf,inf); strcat(outf,"q");
221 if (trd==0)
222 { op=fopen(outf,"w"); fprintf(op,"%3d%5d%3d\n",prime,qdim,ngens);}
223 opnmat(m0,ngens,dim,sdim+1);
224 fclose(op);
225 }
226 if (trd)
227 { if (full==0) dim=sdim;
228 /* Retranspose */
229 ip=fopen(temp,"r"); op=fopen(outf,"w");
230 fprintf(op,"%3d%5d%3d\n",prime,dim,ngens);
231 for (i=1;i<=ngens;i++)
232 { readmat(m0); trans(m0,m1); printmat(m1);}
233 fclose(op);
234 if (full==0)
235 { dim=qdim; strcpy(outf,inf); strcat(outf,"s");
236 op=fopen(outf,"w"); fprintf(op,"%3d%5d%3d\n",prime,dim,ngens);
237 for (i=1;i<=ngens;i++)
238 { readmat(m0); trans(m0,m1); printmat(m1);}
239 fclose(ip); fclose(op);
240 }
241 }
242 unlink(temp); unlink(temp2); return(0);
243 }
244
seeknln(void)245 void seeknln (void) { while (getc(ip)!='\n');}
246
setpinv(void)247 void setpinv (void)
248 { int i,j;
249 for (i=0;i<prime;i++) pinv[i]=0;
250 for (i=1;i<prime;i++) if (pinv[i]==0) for (j=1;j<prime;j++)
251 if (i*j % prime ==1) { pinv[i]=j; pinv[j]=i; break; }
252 }
253
findc(void)254 int findc (void)
255 { char c; while ((c=getchar())==' ' || c=='\n'); return(c); }
256
digit(int c)257 int digit (int c)
258 { if (c>='1' && c<='9') return(1); else return(0);}
259