1 #include "Global.h"
2 #include "Rat.h"
3
4 #undef TEST_Wbase
5 #undef USE_Old_Wbase
6 #define NO_COORD_IMPROVEMENT /* switch off weight permutation */
7
8 typedef struct {Long x[AMBI_Dmax][AMBI_Dmax]; int n, N;} CWLatticeBasis;
9
10 void Make_CWS_Points(CWS *_C, PolyPointList *_P);
11 void Make_RGC_Points(CWS *Cin, PolyPointList *_P);
12 void CWS_to_PermCWS(CWS *Cin, CWS *C, int *pi);
13
14 /* ========== I/O functions: ========== */
15
IsNextDigit()16 int IsNextDigit(){
17 char c; c=fgetc(inFILE); ungetc(c,inFILE);
18 if(c=='0') return -1;
19 if((c<'0') || ('9'<c)) return 0; else return 1;
20 }
21
Print_PPL(PolyPointList * _P,const char * comment)22 void Print_PPL(PolyPointList *_P, const char *comment){
23 int i,j;
24 if(_P->np>20){
25 fprintf(outFILE,"%d %d %s\n",_P->np,_P->n,comment);
26 for(i=0;i<_P->np;i++) {
27 for(j=0;j<_P->n;j++) fprintf(outFILE,"%d ",(int) _P->x[i][j]);
28 fprintf(outFILE,"\n");}}
29 else {
30 fprintf(outFILE,"%d %d %s\n",_P->n,_P->np,comment);
31 for(i=0;i<_P->n;i++) {
32 for(j=0;j<_P->np;j++) fprintf(outFILE," %4d",(int) _P->x[j][i]);
33 fprintf(outFILE,"\n");}}
34 }
35
Print_VL(PolyPointList * _P,VertexNumList * _V,const char * comment)36 void Print_VL(PolyPointList *_P, VertexNumList *_V, const char *comment){
37 int i,j;
38 if(_V->nv>20){
39 fprintf(outFILE,"%d %d %s\n",_V->nv,_P->n,comment);
40 for(i=0;i<_V->nv;i++) {
41 for(j=0;j<_P->n;j++) fprintf(outFILE,"%d ",(int) _P->x[_V->v[i]][j]);
42 fprintf(outFILE,"\n");}}
43 else {
44 fprintf(outFILE,"%d %d %s\n",_P->n,_V->nv,comment);
45 for(i=0;i<_P->n;i++) {
46 for(j=0;j<_V->nv;j++) fprintf(outFILE," %4d",(int) _P->x[_V->v[j]][i]);
47 fprintf(outFILE,"\n");}}
48 }
49
Print_EL(EqList * _E,int * n,int suppress_c,const char * comment)50 void Print_EL(EqList *_E, int *n, int suppress_c, const char *comment){
51 int i,j;
52 fprintf(outFILE,"%d %d %s\n",_E->ne,*n,comment);
53 for(i=0;i<_E->ne;i++) {
54 for(j=0;j<*n;j++) fprintf(outFILE," %3d",(int) _E->e[i].a[j]);
55 if (!suppress_c) fprintf(outFILE," %5d",(int) _E->e[i].c);
56 fprintf(outFILE,"\n");}
57 }
58
Print_Matrix(Long Matrix[][VERT_Nmax],int n_lines,int n_columns,const char * comment)59 void Print_Matrix(Long Matrix[][VERT_Nmax], int n_lines, int n_columns,
60 const char *comment){
61 int i,j;
62 fprintf(outFILE,"%d %d %s\n",n_lines, n_columns, comment);
63 for(i=0;i<n_lines;i++) {
64 for(j=0;j<n_columns;j++) fprintf(outFILE," %3d",(int) Matrix[i][j]);
65 fprintf(outFILE,"\n");}
66 }
auxString2Int(char * c,int * n)67 int auxString2Int(char *c,int *n)
68 { int j=0; *n=0; while(('0'<=c[j])&&(c[j]<='9')) *n = 10*(*n)+c[j++]-'0';
69 if(j) while(c[j]==' ') j++;
70 return j;
71 }
CWSZerror(char * c)72 void CWSZerror(char *c)
73 { printf("Format error %s in Read_CWS_Zinfo\n",c);exit(0);
74 }
Print_CWS_Zinfo(CWS * CW)75 void Print_CWS_Zinfo(CWS *CW)
76 { int i,j; if(CW->nw) for(i=0;i<CW->nz;i++)
77 { fprintf(outFILE,"/Z%d: ",CW->m[i]);
78 for(j=0;j<CW->N;j++) fprintf(outFILE,"%d ",CW->z[i][j]);
79 }
80 }
Read_CWS_Zinfo(FILE * inFILE,CWS * CW)81 int Read_CWS_Zinfo(FILE *inFILE,CWS *CW) /* return !EOF */
82 { int *nz=&CW->nz; int i=0,n; char c[999],b=' '; *nz=0;
83 for(n=0;n<999;n++)
84 { c[n]=fgetc(inFILE); if(feof(inFILE)) return 0; if(c[n]=='\n') break;
85 } if(n==999) {puts("Out of space in Read_CWS_Zinfo");exit(0);}
86 while(c[i]==b)i++;
87 if((c[i]=='=')&&(c[i+1]=='d'))i+=2;
88 while(c[i]==b)i++;
89 /* printf("Read_CWS_Zinfo n=%d -> ",n);
90 {int I; for(I=0;I<n;I++)printf("%c",c[I]);puts("");} */
91 while(i<n)
92 { int j, k, s; if((c[i]!='/')||(c[i+1]!='Z')) return 1;
93 i+=2; assert(*nz<POLY_Dmax);
94 if((j=auxString2Int(&c[i],&CW->m[*nz])))
95 { if(c[i+j]!=':') CWSZerror(":"); else {i+=j+1; while(c[i]==b)i++;}
96 for(k=0;k<CW->N;k++)
97 { if((j=auxString2Int(&c[i],&CW->z[*nz][k])))
98 {if((j)) i+=j; else CWSZerror("missing");} /* ????? */
99 } s=0;
100 for(k=0;k<CW->N;k++) s+=CW->z[*nz][k];
101 /* if(s % CW->m[*nz]) CWSZerror("det!=1"); */
102 (CW->nz)++;
103 } else CWSZerror("Order");
104 } return 1;
105 }
106
checkDimension(int polyDim,int codim,int index)107 void checkDimension(int polyDim, int codim, int index){
108 if (index == 1){
109 if (POLY_Dmax < (polyDim + codim - 1)){
110 if (codim > 1){
111 printf("Please increase POLY_Dmax to at least %d = %d + %d - 1\n",
112 (polyDim + codim - 1), polyDim, codim);
113 printf("(POLY_Dmax >= dim N + codim - 1 is required)\n");}
114 else printf("Please increase POLY_Dmax to at least %d\n", polyDim);
115 exit(0); }}
116 else if (POLY_Dmax < polyDim +1){
117 printf("Please increase POLY_Dmax to at least %d = %d + 1\n",
118 (polyDim + 1), polyDim);
119 printf("(option -G requires POLY_Dmax >= dim(cone) = dim(support) + 1)\n"
120 );
121 exit(0);}
122 }
123
ReadCwsPp(CWS * _CW,PolyPointList * _P,int codim,int index)124 int ReadCwsPp(CWS *_CW, PolyPointList *_P, int codim, int index)
125 /* _P is always an M-lattice polytope
126 codim = 1, index = 1: CY hypersurface
127 codim > 1, index = 1: _P reflexive, CICY with codimension codim
128 codim = 1, index > 1: _P a Gorenstein polytope with index index,
129 determines a reflexive Gorenstein cone */
130 { int i, j, FilterFlag=(inFILE==NULL);
131 int IN[AMBI_Dmax*(AMBI_Dmax+1)], S;
132 static int InputOK;
133 _CW->nw=_CW->N=_CW->nz=0;
134 _CW->index = index;
135 if(FilterFlag) inFILE=stdin;
136 else if(inFILE==stdin) {
137 puts("Degrees and weights `d1 w11 w12 ... d2 w21 w22 ...'");
138 puts(" or `#lines #columns' (= `PolyDim #Points' or `#Points PolyDim'):");}
139 for(i=0;i<AMBI_Dmax*(AMBI_Dmax+1);i++)
140 { char c;
141 while(' '==(c=fgetc(inFILE )));
142 ungetc(c,inFILE); /* read blanks */
143 if(IsNextDigit()) fscanf(inFILE,"%d",&IN[i]); else break;
144 }
145 if(i==0) { if(!InputOK){puts("-h gives you help\n"); exit(0);}
146 else return 0; } InputOK++;
147 if(i==1) { puts("Error in INPUT: need at least 2 numbers!"); exit(0);}
148 if(i==2) /* READ PolyPointList */
149 { int tr=0;
150 while('\n'!=fgetc(inFILE)); /* read to end of line */
151 if(IN[0]==IN[1]){
152 puts("The number of points must be larger than the dimension!");
153 exit(0);}
154 if(IN[0]>IN[1]) {tr=IN[0];IN[0]=IN[1];IN[1]=tr;} tr=!tr;
155 checkDimension(IN[0], codim, index);
156 if(IN[1]>POINT_Nmax) {puts("Please increase POINT_Nmax"); exit(0);}
157 if((inFILE==stdin)&&!FilterFlag)
158 printf("Type the %d coordinates as %s=%d lines with %s=%d columns:\n",
159 IN[0]*IN[1], tr ? "dim" : "#pts", tr ? IN[0] : IN[1],
160 tr ? "#pts" : "dim", tr ? IN[1] : IN[0]);
161 _P->n=IN[0]; _P->np=IN[1];
162 /* allow all numbers in one string (distributed over lines) or
163 * matrix blocks with trailing comments */
164 if(tr) for(i=0;i<IN[0];i++) for(j=0;j<IN[1];j++)
165 { int X; fscanf(inFILE,"%d",&X); _P->x[j][i]=X; }
166 else for(i=0;i<IN[1];i++) for(j=0;j<IN[0];j++)
167 { int X; fscanf(inFILE,"%d",&X); _P->x[i][j]=X; }
168 /* Finish_Poly_Points(_P); */
169 while(fgetc(inFILE )-'\n') if(feof(inFILE)) return 0;/* read to EOL */
170 if(FilterFlag) inFILE=NULL;
171 return 1;
172 } /* End of reading PolyPointList */
173 assert(i!=3);
174 S=IN[i-1];
175 for(j=0;j<i-1;j++) if((IN[j]==0)||(S<IN[j])) break;
176 if(j==i-1) /* Single Weights with: "w1 ... d" */
177 { _CW->nw=1; _CW->d[0]=S; _CW->N=i-1;
178 assert(_CW->N <= POLY_Dmax+1); /* Increase POLY_Dmax */
179 for(j=0;j<_CW->N;j++) {_CW->W[0][j]=IN[j]; assert(IN[j]>0);}
180 goto MAP;
181 }
182 _CW->d[0] = IN[0]; /* ASSIGN (COMBINED) WEIGHTS */
183 S = IN[0] * index;
184 for(j=1; j<i; j++)
185 { if(S<IN[j])
186 { if(S) { puts("Error in INPUT: degree vs. weights!"); exit(0);}
187 else
188 { _CW->d[++(_CW->nw)]=IN[j];
189 S = IN[j] * index;
190 if(1==_CW->nw)
191 { if(i%(1+_CW->N)) {puts("INPUT error: numbers?"); exit(0);}
192 } if(j%(1+_CW->N)) {puts("INPUT error: degrees?"); exit(0);}
193 }
194 }
195 else
196 { _CW->W[_CW->nw][j - 1 - (1 + _CW->N)*_CW->nw] = IN[j]; S-=IN[j];
197 if(0==_CW->nw) if(++(_CW->N)>AMBI_Dmax)
198 { puts("Increase AMBI_Dmax!"); exit(0);
199 }
200 }
201 }
202 ++(_CW->nw);
203 checkDimension(_CW->N - _CW->nw, codim, index);
204
205 MAP: for(i=0;i<_CW->nw;i++) /* check consistency of CWS */
206 { Long sum = _CW->d[i] * index, *w = _CW->W[i];
207 for(j=0;j<_CW->N;j++) { assert(w[j]>=0); sum-=w[j]; }
208 if((sum)&&(index==1)){
209 printf("Use option -l for (single) WeightSystems with ");
210 printf("d!=\\sum(w)\n(only Read_Weight makes the correct ");
211 puts("PolyPointList in that case)"); exit(0);}
212 }
213 _CW->nz=0;
214 if(!Read_CWS_Zinfo(inFILE,_CW)) /* read Z to EOL */
215 { if(!InputOK) puts("-h gives you help\n"); return 0;}
216 Make_CWS_Points(_CW, _P);
217 if(FilterFlag) inFILE=NULL;
218 return 1;
219 }
220
Read_CWS_PP(CWS * _CW,PolyPointList * _P)221 int Read_CWS_PP(CWS *_CW, PolyPointList *_P){
222 return ReadCwsPp(_CW, _P, 1, 1);
223 }
224
Read_PP(PolyPointList * _P)225 int Read_PP(PolyPointList *_P)
226 { int i, j, FilterFlag=(inFILE==NULL);
227 int IN[AMBI_Dmax*(AMBI_Dmax+1)];
228 static int InputOK;
229 /* _CW->nw=_CW->N=_CW->nz=0; */
230
231 if(FilterFlag) inFILE=stdin;
232 else if(inFILE==stdin) {
233 printf("`#lines #columns' (= `PolyDim #Points' or `#Points PolyDim'):\n");
234 };
235 for(i=0;i<AMBI_Dmax*(AMBI_Dmax+1);i++)
236 { char c;
237 while(' '==(c=fgetc(inFILE )));
238 ungetc(c,inFILE); /* read blanks */
239 if(IsNextDigit()) fscanf(inFILE,"%d",&IN[i]); else break;
240 }
241 if(i==0) { if(!InputOK){puts("-h gives you help\n"); exit(0);}
242 else return 0; } InputOK++;
243 if(i==1) { puts("Error in INPUT: need at least 2 numbers!"); exit(0);}
244 if(i==2) /* READ PolyPointList */
245 { int tr=0;
246 while('\n'!=fgetc(inFILE)); /* read to end of line */
247 if(IN[0]==IN[1]){
248 puts("The number of points must be larger than the dimension!");
249 exit(0);}
250 if(IN[0]>IN[1]) {tr=IN[0];IN[0]=IN[1];IN[1]=tr;} tr=!tr;
251 if(IN[0]>POLY_Dmax) {puts("increase POLY_Dmax!"); exit(0);}
252 if(IN[1]>POINT_Nmax) {puts("increase POINT_Nmax!"); exit(0);}
253 if((inFILE==stdin)&&!FilterFlag)
254 printf("Type the %d coordinates as %s=%d lines with %s=%d columns:\n",
255 IN[0]*IN[1], tr ? "dim" : "#pts", tr ? IN[0] : IN[1],
256 tr ? "#pts" : "dim", tr ? IN[1] : IN[0]);
257 _P->n=IN[0]; _P->np=IN[1];
258 /* allow all numbers in one string (distributed over lines) or
259 * matrix blocks with trailing comments */
260 if(tr) for(i=0;i<IN[0];i++) for(j=0;j<IN[1];j++)
261 { int X; fscanf(inFILE,"%d",&X); _P->x[j][i]=X; }
262 else for(i=0;i<IN[1];i++) for(j=0;j<IN[0];j++)
263 { int X; fscanf(inFILE,"%d",&X); _P->x[i][j]=X; }
264 /* Finish_Poly_Points(_P); */
265 while(fgetc(inFILE )-'\n') if(feof(inFILE)) return 0;/* read to EOL */
266 if(FilterFlag) inFILE=NULL;
267 return 1;
268 }
269 if(i>2){puts("Error: expected input format is matrix of polytope points!");exit(0);}
270 if(FilterFlag) inFILE=NULL;
271 return 1;
272 }
273
Read_CWS(CWS * _CW,PolyPointList * _P)274 int Read_CWS(CWS *_CW, PolyPointList *_P)
275 { int i, j, FilterFlag=(inFILE==NULL);
276 int IN[AMBI_Dmax*(AMBI_Dmax+1)], S;
277 static int InputOK;
278 _CW->nw=_CW->N=_CW->nz=0;
279 _CW->index = 1;
280
281 if(FilterFlag) inFILE=stdin;
282 else if(inFILE==stdin) {
283 printf("Degrees and weights `d1 w11 w12 ... d2 w21 w22 ...':\n");
284 };
285 for(i=0;i<AMBI_Dmax*(AMBI_Dmax+1);i++)
286 { char c;
287 while(' '==(c=fgetc(inFILE )));
288 ungetc(c,inFILE); /* read blanks */
289 if(IsNextDigit()) fscanf(inFILE,"%d",&IN[i]); else break;
290 }
291 if(i==0) { if(!InputOK){puts("-h gives you help\n"); exit(0);}
292 else return 0; } InputOK++;
293 if(i==1) { puts("Error in INPUT: need at least 2 numbers!"); exit(0);}
294 if(i==2) {puts("Error: expected input format is CWS!"); exit(0);}
295 assert(i!=3);
296 S=IN[i-1]; for(j=0;j<i-1;j++) if(S<IN[j]) break;
297 if(j==i-1) /* Single Weights with: "w1 ... d" */
298 { _CW->nw=1; _CW->d[0]=S; _CW->N=i-1;
299 assert(_CW->N <= POLY_Dmax+1); /* Increase POLY_Dmax */
300 for(j=0;j<_CW->N;j++) {_CW->W[0][j]=IN[j]; assert(IN[j]>0);}
301 goto MAP;
302 }
303 S=_CW->d[_CW->nw]=IN[0]; /* ASSIGN WEIGHTS */
304 for(j=1; j<i; j++)
305 { if(S<IN[j])
306 { if(S) { puts("Error in INPUT: degree vs. weights!"); exit(0);}
307 else
308 { S=_CW->d[++(_CW->nw)]=IN[j];
309 if(1==_CW->nw)
310 { if(i%(1+_CW->N)) {puts("INPUT error: numbers?"); exit(0);}
311 } if(j%(1+_CW->N)) {puts("INPUT error: degrees?"); exit(0);}
312 }
313 }
314 else
315 { _CW->W[_CW->nw][j - 1 - (1 + _CW->N)*_CW->nw] = IN[j]; S-=IN[j];
316 if(0==_CW->nw) if(++(_CW->N)>AMBI_Dmax)
317 { puts("Increase AMBI_Dmax!"); exit(0);
318 }
319 }
320 }
321 ++(_CW->nw);
322 if(_CW->N-_CW->nw > POLY_Dmax){
323 printf("Please increase POLY_Dmax to at least %d\n", _CW->N-_CW->nw);
324 exit(0);} /* increase POLY_Dmax */
325
326 MAP: for(i=0;i<_CW->nw;i++) /* check consistency of CWS */
327 { Long sum=_CW->d[i], *w=_CW->W[i];
328 for(j=0;j<_CW->N;j++) { assert(w[j]>=0); sum-=w[j]; }
329 if(sum){ /*printf("Use poly.x -w for (single) WeightSystems with ");
330 printf("d!=\\sum(w)\n(only Read_Weight makes the correct ");
331 puts("PolyPointList in that case)"); exit(0);*/
332 printf("cannot handle (single) WeightSystems with ");
333 printf("d!=\\sum(w)\n");
334 exit(0);}}
335 _CW->nz=0;
336 if(!Read_CWS_Zinfo(inFILE,_CW)) /* read Z to EOL */
337 { if(!InputOK) puts("-h gives you help\n"); return 0;}
338 Make_CWS_Points(_CW, _P); /* now make POLY */
339 if(FilterFlag) inFILE=NULL;
340 return 1;
341 }
342
Print_CWH(CWS * _W,BaHo * _BH)343 void Print_CWH(CWS *_W, BaHo *_BH){
344 int i, j;
345 for(i=0;i<_W->nw;i++) {
346 fprintf(outFILE,"%d ",(int) _W->d[i]);
347 for(j=0;j<_W->N;j++) fprintf(outFILE,"%d ",(int) _W->W[i][j]);
348 if(i+1<_W->nw) fprintf(outFILE," "); }
349 Print_CWS_Zinfo(_W);
350 if(_BH->np) {
351 fprintf(outFILE,"M:%d %d N:%d %d",_BH->mp,_BH->mv,_BH->np,_BH->nv);
352 if (_BH->n == 3) fprintf(outFILE," Pic:%d Cor:%d",_BH->pic, _BH->cor);
353 if (_BH->n > 3) {
354 fprintf(outFILE," H:%d",_BH->h1[1]);
355 for(i=2;i<_BH->n-1;i++) fprintf(outFILE,",%d",_BH->h1[i]);
356 if(_BH->n==4) fprintf(outFILE," [%d]",2*(_BH->h1[1]-_BH->h1[2]));
357 if(_BH->n==5)
358 fprintf(outFILE," [%d]",48+6*(_BH->h1[1]-_BH->h1[2]+_BH->h1[3]));}}
359 else if(_BH->mp)
360 fprintf(outFILE,"M:%d %d F:%d",_BH->mp,_BH->mv,_BH->nv);
361 else fprintf(outFILE,"V:%d F:%d",_BH->mv,_BH->nv);
362 fprintf(outFILE,"\n");
363 }
364
365 /* ========== END of I/O functions ========== */
366
367
368 /* ========== Solve W-Eq. -> triangular LatticeBasis
369 *
370 * B[0]:= (-n1, n0, 0, ...) / GCD(n0,n1);
371 * B[i]:= (0,...,0,g/G,0,...)- (ni/G) * ExtGCD.K(n0,...,n(i-1),0,...);
372 * with g=GCD(n0,...,n[i-1]); G=GCD(g,ni);
373 *
374 * CWS by iteration: NextWeight=W[]*B; NewB=Basis(NextW); B[i+1]=B[i]*NewB;
375 * */
376
PrintBasis(CWLatticeBasis * _B)377 void PrintBasis(CWLatticeBasis *_B)
378 { int i,j; puts("Basis:");
379 for(i=0;i<_B->n;i++)
380 { for(j=0;j<_B->N;j++) fprintf(outFILE,"%6d ",(int) _B->x[i][j]);
381 puts("");
382 } puts("End of Basis - -");
383 }
384
Orig_Solve_Next_WEq(Long * NW,CWLatticeBasis * _B)385 void Orig_Solve_Next_WEq(Long *NW, CWLatticeBasis *_B)
386 { int i, j, P=0, p[AMBI_Dmax]; Long W[AMBI_Dmax], G; _B->n=_B->N-1;
387 for(i=0;i<_B->N;i++)
388 { for(j=0;j<_B->n;j++) _B->x[j][i]=0; /* init B.x=0 */
389 if(NW[i]) {p[P]=i; W[P++]=NW[i];} /* non-zero weights */
390 }
391 if(P<2) puts("need two non-zero weights in >>Solve_Next_WEq<<");
392 for(i=0;i<p[0];i++) _B->x[i][i]=1;
393 while((++i)<p[1]) _B->x[i-1][i]=1;
394 G=Fgcd(W[0],W[1]); if(W[0]/G<0) G=-G;
395 _B->x[i-1][p[0]]=-W[1]/G; _B->x[i-1][p[1]]=W[0]/G;
396 j=2; while(++i<_B->N)
397 { if(NW[i])
398 { int k; Long *X=_B->x[i-1], K[AMBI_Dmax], g=REgcd(W, &j, K);
399 G=Fgcd(g,NW[i]); if(g/G<0) G=-G; X[i]= g/G; g=W[j]/G;
400 for(k=0;k<j;k++) X[p[k]]=-K[k]*g;
401 j++;
402 }
403 else _B->x[i-1][i]=1;
404 }
405 }
Solve_Next_WEq(Long * NW,CWLatticeBasis * _B)406 void Solve_Next_WEq(Long *NW, CWLatticeBasis *_B)
407 { Long W[AMBI_Dmax], *X[AMBI_Dmax], GLZ[AMBI_Dmax][AMBI_Dmax];
408 int i, j, P=0, p[AMBI_Dmax]; _B->n=_B->N-1;
409 #ifdef TEST_Wbase
410 Orig_Solve_Next_WEq(NW,_B); PrintBasis(_B);
411 #endif
412 for(i=0;i<_B->N;i++)
413 { for(j=0;j<_B->n;j++) _B->x[j][i]=0; /* init B.x=0 */
414 if(NW[i]) {p[P]=i; X[P]=GLZ[P]; W[P++]=NW[i];} /* non-zero weights */
415 } if(P>1) W_to_GLZ(W,&P,X); /* P>1, compute GLZ */
416 else {/* printf("P=%d W[0]=%d for W_to_GLZ\n",P,W[0]);exit(0);*/assert(P);
417 for(i=0;i<p[0];i++)_B->x[i][i]=1;
418 while((++i)<_B->N)_B->x[i-1][i]=1;
419 return; }
420 for(i=1;i<P;i++) if(X[i][i]<0) for(j=0;j<=i;j++) X[i][j] *= -1;
421 for(i=0;i<p[0];i++) _B->x[i][i]=1;
422 while((++i)<p[1]) _B->x[i-1][i]=1;
423 _B->x[i-1][p[0]]=X[1][0]; _B->x[i-1][p[1]]=X[1][1]; j=2;
424 while(++i<_B->N)
425 { if(NW[i])
426 { int k; Long *B=_B->x[i-1]; for(k=0;k<=j;k++) B[p[k]]=X[j][k]; j++;
427 }
428 else _B->x[i-1][i]=1;
429 }
430 #ifdef TEST_Wbase
431 /* for(i=0;i<_B->N;i++)printf("%d ",NW[i]);puts("=NW");for(i=0;i<P;i++)
432 {for(j=0;j<P;j++)printf("%4d ",GLZ[i][j]);puts("=GL");} */
433 printf("New version: "); PrintBasis(_B);
434 #endif
435 #ifdef USE_Old_Wbase
436 Orig_Solve_Next_WEq(NW,_B);
437 #endif
438 }
Make_CWS_Basis(CWS * _C,CWLatticeBasis * _B)439 void Make_CWS_Basis(CWS *_C, CWLatticeBasis *_B){
440 int i,j,k,l; Long W[AMBI_Dmax]; CWLatticeBasis NB, AuxB;
441 AuxB.N=_B->N=_C->N;
442 Solve_Next_WEq(_C->W[0],_B);
443 for(i=1;i<_C->nw;i++){
444 for(j=0;j<_B->n;j++){
445 W[j]=0;
446 for(k=0;k<_B->N;k++) W[j]+=_C->W[i][k]*_B->x[j][k]; }
447 NB.N=_B->n;
448 Solve_Next_WEq(W,&NB);
449 AuxB.n=NB.n;
450 for(j=0;j<AuxB.N;j++) for(k=0;k<AuxB.n;k++) {
451 AuxB.x[k][j]=0;
452 for(l=0;l<NB.N;l++) AuxB.x[k][j] += NB.x[k][l] * _B->x[l][j]; }
453 *_B=AuxB; }
454 /* init coordinate hyperplanes */
455 _C->B.ne=_B->N; /* == transposed of Basis */
456 for(i=0;i<_B->N;i++) {
457 for(j=0;j<_B->n;j++) _C->B.e[i].a[j]=_B->x[j][i];
458 _C->B.e[i].c=1;}
459 }
460
Poly_To_Ambi(CWLatticeBasis * _B,Long * x,Long * X)461 void Poly_To_Ambi(CWLatticeBasis *_B, Long *x, Long *X)
462 { int i, j;
463 for(i=0;i<_B->N;i++)
464 { X[i]=1; for(j=0;j<_B->n;j++) X[i] += _B->x[j][i] * x[j];
465 }
466 }
467 /* B.x[i][A] is lower triangular: zero for i<p and Amin[p] <= A <Amin[p+1]
468 * X^A=x^p B_p^A+1, hence for A \in [ Amin[p],Amin[p+1] ) and all k:
469 * x^iB_i^A \in [0,d_k/w_k^A] - 1 - \sum_{l<i} x^l B_l^A
470 * */
471
472 void QuotZ_2_SublatG(Long Z[][VERT_Nmax],int *zm,Long *M,int *d,
473 Long G[][POLY_Dmax]); /* normalize and diagonalize Z */
CWS_2_SublatZ(CWS * C,CWLatticeBasis * B,int * m,Long * M,Long G[POLY_Dmax][POLY_Dmax])474 void CWS_2_SublatZ(CWS *C,CWLatticeBasis *B, /* in */
475 int *m, Long *M, Long G[POLY_Dmax][POLY_Dmax]) /* out */
476 { int i,j,k,d=C->N-C->nw; Long Z[POLY_Dmax][VERT_Nmax]; *m=C->nz;
477 for(i=0;i<C->nz;i++){ M[i]=C->m[i]; for(j=0;j<d;j++)
478 { Z[i][j]=0; for(k=0;k<C->N;k++) Z[i][j]+=C->z[i][k]*B->x[j][k];}
479 } /* PrintBasis(B); */ QuotZ_2_SublatG(Z,m,M,&d,G);
480 }
Reduce_PPL_2_Sublat(PolyPointList * P,int * nm,Long * M,Long G[][POLY_Dmax])481 void Reduce_PPL_2_Sublat(PolyPointList *P,int *nm,Long *M,Long G[][POLY_Dmax])
482 { int i,j,n=0,N; for(N=0;N<P->np;N++)
483 { Long X[POLY_Dmax]; for(i=0;i<P->n;i++)
484 { X[i]=0; for(j=0;j<P->n;j++) X[i]+=G[i][j]*P->x[N][j];
485 } for(i=0;i<*nm;i++) if(X[i]%M[i]) break; if(i<*nm) continue;
486 for(i=0;i<*nm;i++) P->x[n][i]=X[i]/M[i];
487 for(i=*nm;i<P->n;i++) P->x[n][i]=X[i];
488 n++;
489 } P->np=n; /* Print_PPL(P,"");printf("nm=%d\n",*nm); */
490 }
PD_Floor(Long N,Long D)491 Long PD_Floor(Long N,Long D) /* assuming PosDenom D>0: F <= N/D < F+1 */
492 { Long F=N/D; return (F*D>N) ? F-1 : F;
493 }
494
Old_Make_CWS_Points(CWS * Cin,PolyPointList * _P)495 void Old_Make_CWS_Points(CWS *Cin, PolyPointList *_P)
496 { int i, j, Amin[POLY_Dmax+1]; Long *x=_P->x[_P->np=0], xmin[POLY_Dmax],
497 xmax[POLY_Dmax], Xmax[AMBI_Dmax], xaux[POLY_Dmax], L, R; CWS *_C=Cin;
498 CWLatticeBasis B; Long G[POLY_Dmax][POLY_Dmax],M[POLY_Dmax];int m=Cin->nz;
499 #ifndef NO_COORD_IMPROVEMENT /* ==== Perm Coord Improvement ==== */
500 int pi[AMBI_Dmax]; CWS Caux; _C=&Caux; CWS_to_PermCWS(Cin,_C, pi);
501 #endif /* = End of Perm Coord Improvement = */
502 Make_CWS_Basis(_C, &B); /* make `triangular' Basis */
503 #ifndef NO_COORD_IMPROVEMENT /* ==== Perm Coord Improvement ==== */
504 assert(_C->N==_C->B.ne);
505 for(i=0;i<Cin->N;i++) Cin->B.e[i]=_C->B.e[pi[i]]; Cin->B.ne=_C->N;
506 #endif /* = End of Perm Coord Improvement = */
507 i=_P->n=B.n; Amin[0]=0; Amin[B.n]=j=B.N; /* inversion structure */
508 while(--i) {while(!B.x[i-1][--j]); Amin[i]=++j;}
509 for(i=0;i<B.N;i++) /* compute Xmax */
510 { Xmax[i]=0; for(j=0;j<_C->nw;j++) if(_C->W[j][i])
511 { L=_C->d[j]/_C->W[j][i];
512 if(Xmax[i]) {if(L<Xmax[i]) Xmax[i]=L;}
513 else Xmax[i]=L;
514 }
515 }
516 j=B.n-1; i=Amin[j+1]-1; R=B.x[j][i]; /* compute xmin[j] and xmax[j] */
517 xmin[j]=-PD_Floor(1,R); xmax[j]=PD_Floor(Xmax[i]-1,R); /* since R > 0 */
518 /** /printf("R=%2d: %2d <= x[%d=&%d] <= %2d\n",R,xmin[j],j,i,xmax[j]);/ **/
519 while((i--)>Amin[j]) /* if(R=B.x[B.n-1][i]): R!=0 => new limits */
520 { Long Low=-1, Upp=Low+Xmax[i]; R=B.x[B.n-1][i];
521 if(R>0) {if(xmax[j]>(L=PD_Floor(Upp,R))) xmax[j]=L;
522 if(xmin[j]<(L=-PD_Floor(-Low,R))) xmin[j]=L;}
523 else {if(xmax[j]>(L=PD_Floor(-Low,-R))) xmax[j]=L;
524 if(xmin[j]<(L=-PD_Floor(Upp,-R))) xmin[j]=L;}
525 /** /printf("R=%2d: %2d <= x[%d=&%d] <= %2d\n",R,xmin[j],j,i,xmax[j]);/ **/
526 } /* this completes the limits for x[B.n-1] */
527 x[j]=xmin[j]; while(j<B.n)
528 { int k; if(x[j]>xmax[j]) {if(B.n==(++j)) break; else x[j]++;}
529 else /* compute limits[j-1] and initialize x[j-1] */
530 { Long Low=-1, Upp=Xmax[i=Amin[j--]-1]; int RangeFlag=0;
531 for(k=j+1;k<B.n;k++) Low-=x[k]*B.x[k][i]; /* compute offset */
532 Upp+=Low; R=B.x[j][i];
533 xmin[j]=-PD_Floor(-Low,R); xmax[j]=PD_Floor(Upp,R);
534 /** /printf("R=%2d: %2d <= x[%d=&%d] <= %2d\n",R,xmin[j],j,i,xmax[j]);/ **/
535 while((i--)>Amin[j]) if((R=B.x[j][i])) /* R!=0 => new limits */
536 { Low=-1; Upp=Xmax[i];
537 for(k=j+1;k<B.n;k++) Low-=x[k]*B.x[k][i];
538 Upp+=Low; R=B.x[j][i];
539 if(R>0) {if(xmax[j]>(L=PD_Floor(Upp,R))) xmax[j]=L;
540 if(xmin[j]<(L=-PD_Floor(-Low,R))) xmin[j]=L;}
541 else {if(xmax[j]>(L=PD_Floor(-Low,-R))) xmax[j]=L;
542 if(xmin[j]<(L=-PD_Floor(Upp,-R))) xmin[j]=L;}
543 /** /printf("R=%2d: %2d <= x[%d=&%d] <= %2d\n",R,xmin[j],j,i,xmax[j]);/ **/
544 } /* completes limits for x[] except */
545 else /* when R=0 and X out of Range: */
546 { Long X=1; for(k=j+1;k<B.n;k++) X+=x[k]*B.x[k][i];
547 if((X<0)||(X>Xmax[i])) RangeFlag=1;
548 }
549 if(RangeFlag) ++x[++j]; else x[j]=xmin[j];
550 if(j==0)
551 { while(x[0]<=xmax[0])
552 { Long *y;
553 if((++_P->np) < POINT_Nmax)
554 { y=(_P->x[_P->np]); for(k=0;k<B.n;k++) y[k]=x[k];
555 }
556 else if(_P->np == POINT_Nmax)
557 { y=xaux; for(k=0;k<B.n;k++) y[k]=x[k];
558 }
559 else {puts("Increase POINT_Nmax");exit(0);}
560 x=y; ++x[0];
561 }
562 x[j=1]++;
563 }
564 }
565 } if(m)CWS_2_SublatZ(_C,&B,&m,M,G);if(m)Reduce_PPL_2_Sublat(_P,&m,M,G);
566 }
567
Compute_X0(int N,CWS * _C,Long * X0)568 int Compute_X0(int N, CWS *_C, Long *X0){
569 /* finds a first X = X0 such that sum w_i X_i = d for every weight w in _C,
570 this becomes the origin after the change of coordinates from X to x */
571 int j;
572 Long Xmax=0;
573 if (!N){
574 for(j=0;j<_C->nw;j++)
575 if(_C->W[j][0]){
576 if (_C->d[j]%_C->W[j][N]) return 0;
577 if(Xmax) {if (Xmax != _C->d[j]/_C->W[j][N]) return 0;}
578 else Xmax = _C->d[j]/_C->W[j][N]; }
579 else if(_C->d[j]) return 0;
580 X0[0] = Xmax;
581 return 1;}
582 for(j=0;j<_C->nw;j++) if(_C->W[j][N]!=0){
583 Long L=_C->d[j]/_C->W[j][N];
584 if(Xmax) {if(L<Xmax) Xmax = L;}
585 else Xmax = L; }
586 for (X0[N]=0; X0[N] <= Xmax; X0[N]++){
587 if (Compute_X0(N-1, _C, X0)){
588 for(j=0; j<_C->nw; j++) _C->d[j] += X0[N] * _C->W[j][N]; /* reset Cin*/
589 return 1;}
590 for(j=0; j<_C->nw; j++) _C->d[j] -= _C->W[j][N];}
591 for(j=0; j<_C->nw; j++) _C->d[j] += (Xmax + 1) * _C->W[j][N]; /* reset Cin */
592 return 0;
593 }
594
Make_CWS_Points(CWS * Cin,PolyPointList * _P)595 void Make_CWS_Points(CWS *Cin, PolyPointList *_P)
596 { int i, j, Amin[POLY_Dmax+1], m=Cin->nz;
597 Long *x=_P->x[_P->np=0], xmin[POLY_Dmax],
598 xmax[POLY_Dmax], Xmax[AMBI_Dmax], X0[AMBI_Dmax], xaux[POLY_Dmax], L, R;
599 CWS *_C=Cin;
600 CWLatticeBasis B; Long G[POLY_Dmax][POLY_Dmax],M[POLY_Dmax];
601 #ifndef NO_COORD_IMPROVEMENT /* ==== Perm Coord Improvement ==== */
602 int pi[AMBI_Dmax]; CWS Caux; _C=&Caux; CWS_to_PermCWS(Cin,_C, pi);
603 #endif /* = End of Perm Coord Improvement = */
604 Make_CWS_Basis(_C, &B); /* make `triangular' Basis */
605 #ifndef NO_COORD_IMPROVEMENT /* ==== Perm Coord Improvement ==== */
606 assert(_C->N==_C->B.ne);
607 for(i=0;i<Cin->N;i++) Cin->B.e[i]=_C->B.e[pi[i]]; Cin->B.ne=_C->N;
608 #endif /* = End of Perm Coord Improvement = */
609 if (Cin->index == 1) for (i=0; i< Cin->N; i++) X0[i] = 1;
610 else if(!Compute_X0(Cin->N - 1, Cin, X0)){_P->n=0; puts("no X0!");return;}
611 /* X0 is the reference point in X-space that is transformed to the
612 origin of x-space for _P */
613 /*printf("\nX0: ");for (i=0;i<Cin->N;i++) printf("%ld ", X0[i]);puts("");
614 for(i=0;i<Cin->nw;i++) {
615 fprintf(outFILE,"%d ",(int) Cin->d[i]);
616 for(j=0;j<Cin->N;j++) fprintf(outFILE,"%d ",(int) Cin->W[i][j]);
617 if(i+1<Cin->nw) fprintf(outFILE," "); }*/
618 i=_P->n=B.n; Amin[0]=0; Amin[B.n]=j=B.N; /* inversion structure */
619 while(--i) {while(!B.x[i-1][--j]); Amin[i]=++j;}
620 for(i=0;i<B.N;i++) /* compute Xmax */
621 { Xmax[i]=0; for(j=0;j<_C->nw;j++) if(_C->W[j][i])
622 { L=_C->d[j]/_C->W[j][i];
623 if(Xmax[i]) {if(L<Xmax[i]) Xmax[i]=L;}
624 else Xmax[i]=L;
625 }
626 }
627 j=B.n-1; i=Amin[j+1]-1; R=B.x[j][i]; /* compute xmin[j] and xmax[j] */
628 xmin[j] = -PD_Floor(X0[i],R);
629 xmax[j] = PD_Floor(Xmax[i]-X0[i],R); /* since R > 0 */
630 /** /printf("R=%2d: %2d <= x[%d=&%d] <= %2d\n",R,xmin[j],j,i,xmax[j]);/ **/
631 while((i--)>Amin[j]) /* if(R=B.x[B.n-1][i]): R!=0 => new limits */
632 { Long Low=-X0[i], Upp=Low+Xmax[i]; R=B.x[B.n-1][i];
633 if(R>0) {if(xmax[j]>(L=PD_Floor(Upp,R))) xmax[j]=L;
634 if(xmin[j]<(L=-PD_Floor(-Low,R))) xmin[j]=L;}
635 else {if(xmax[j]>(L=PD_Floor(-Low,-R))) xmax[j]=L;
636 if(xmin[j]<(L=-PD_Floor(Upp,-R))) xmin[j]=L;}
637 /** /printf("R=%2d: %2d <= x[%d=&%d] <= %2d\n",R,xmin[j],j,i,xmax[j]);/ **/
638 } /* this completes the limits for x[B.n-1] */
639 x[j]=xmin[j];
640 while(j<B.n)
641 { int k; if(x[j]>xmax[j]) {if(B.n==(++j)) break; else x[j]++;}
642 else /* compute limits[j-1] and initialize x[j-1] */
643 { Long Upp=Xmax[i=Amin[j--]-1], Low=-X0[i]; int RangeFlag=0;
644 for(k=j+1;k<B.n;k++) Low-=x[k]*B.x[k][i]; /* compute offset */
645 Upp+=Low; R=B.x[j][i];
646 xmin[j]=-PD_Floor(-Low,R); xmax[j]=PD_Floor(Upp,R);
647 /** /printf("R=%2d: %2d <= x[%d=&%d] <= %2d\n",R,xmin[j],j,i,xmax[j]);/ **/
648 while((i--)>Amin[j]) if((R=B.x[j][i])) /* R!=0 => new limits */
649 { Low=-X0[i]; Upp=Xmax[i];
650 for(k=j+1;k<B.n;k++) Low-=x[k]*B.x[k][i];
651 Upp+=Low; R=B.x[j][i];
652 if(R>0) {if(xmax[j]>(L=PD_Floor(Upp,R))) xmax[j]=L;
653 if(xmin[j]<(L=-PD_Floor(-Low,R))) xmin[j]=L;}
654 else {if(xmax[j]>(L=PD_Floor(-Low,-R))) xmax[j]=L;
655 if(xmin[j]<(L=-PD_Floor(Upp,-R))) xmin[j]=L;}
656 /** /printf("R=%2d: %2d <= x[%d=&%d] <= %2d\n",R,xmin[j],j,i,xmax[j]);/ **/
657 } /* completes limits for x[] except */
658 else /* when R=0 and X out of Range: */
659 { Long X=1; for(k=j+1;k<B.n;k++) X+=x[k]*B.x[k][i];
660 if((X<0)||(X>Xmax[i])) RangeFlag=1;
661 }
662 if(RangeFlag) ++x[++j]; else x[j]=xmin[j];
663 if(j==0)
664 { while(x[0]<=xmax[0])
665 { Long *y;
666 if((++_P->np) < POINT_Nmax)
667 { y=(_P->x[_P->np]); for(k=0;k<B.n;k++) y[k]=x[k];
668 }
669 else if(_P->np == POINT_Nmax)
670 { y=xaux; for(k=0;k<B.n;k++) y[k]=x[k];
671 }
672 else {puts("Increase POINT_Nmax");exit(0);}
673 x=y; ++x[0];
674 }
675 x[j=1]++;
676 }
677 }
678 } if(m)CWS_2_SublatZ(_C,&B,&m,M,G);if(m)Reduce_PPL_2_Sublat(_P,&m,M,G);
679 }
680
681 /* ========== Coordinate improvement via CWS-Permutations ========== */
682
683 typedef struct {Long s, g, *N, **G; int *p;} Tri_GLZ_MPaux;
Tri_GLZ_Norm(int * d,Long ** S)684 Long Tri_GLZ_Norm(int *d, Long **S) /* sum-norm = default */
685 { Long norm=0; int i,j; for(i=0;i<*d;i++)for(j=0;j<*d;j++)
686 #ifdef MAX_NORM /* max(|*|) */
687 { Long x=S[i][j]; if(x<0) x=-x; if(x>norm) norm=x; } return norm;
688 #else /* sum(|*|) */
689 { Long x=S[i][j]; norm += (x>0) ? x : -x ; } return norm;
690 #endif
691 }
Tri_GLZ_Basis_Perm(int * d,int * pi,Tri_GLZ_MPaux * AP)692 void Tri_GLZ_Basis_Perm(int *d,int *pi, /* int *pinv, */ Tri_GLZ_MPaux * AP)
693 { Long g, norm, N[AMBI_Dmax], M[AMBI_Dmax][AMBI_Dmax], *S[AMBI_Dmax];
694 int i, j; /* if(pi[0]>pi[1]) return; *//* eliminate equiv. permut. */
695 for(i=0;i<*d;i++) { N[i]=AP->N[pi[i]]; S[i]=M[i];}
696 g=W_to_GLZ(N,d,S); norm=Tri_GLZ_Norm(d,S);
697 if(AP->s) assert(g==AP->g); else AP->g=g;
698 #ifdef TEST_MIN_GLZ
699 { Long err=0, max=0, tsum=0, pos; for(i=0;i<*d;i++)
700 { for(j=0;j<*d;j++) {pos=S[i][j]; if(pos<0)pos=-pos; tsum+=pos;
701 if(max<pos)max=pos;} if(err==0) {
702 for(j=0;j<*d;j++) err+=N[j]*S[i][j]; if(i==0)err-=g;}
703 } assert((norm==tsum)||(norm==max));
704 printf("max=%lld, sum=%lld\n",(long long) max, (long long) tsum);
705 if(err) { for(i=0;i<*d;i++) printf("%d",pi[i]);
706 printf(" g=%d norm=%d",(int)g,(int)norm); puts("");
707 for(i=0;i<*d;i++){printf("S[%d]=",i);for(j=0;j<*d;j++)printf(" %5d",
708 S[i][j]);if(!i)printf(" g=%d norm=%d",(int)g,(int)norm);puts("");}}
709 }
710 #endif
711 if((0 == AP->s) || (norm < AP->s)) /* init or improve AP->G */
712 { for(i=0;i<*d;i++)for(j=0;j<*d;j++)AP->G[i][pi[j]]=S[i][j]; AP->s=norm;
713 for(i=0;i<*d;i++) AP->p[i]=pi[i];
714 }
715 }
716
717 #ifndef NO_COORD_IMPROVEMENT
718 /* improved by permutation pi=pinv^-1 */
Wperm_to_GLZ(Long * W,int * d,Long ** G,int * P)719 Long Wperm_to_GLZ(Long *W, int *d, Long **G, int *P)
720 { Tri_GLZ_MPaux AS; int i, j, pi[AMBI_Dmax], pinv[AMBI_Dmax];
721 AS.s=0; AS.N=W; AS.G=G; AS.p=P;
722 Map_Permut(d,pi, /* pinv,*/ (Tri_GLZ_Basis_Perm),(void *) &AS);
723 for(i=0;i<*d;i++)for(j=0;j<*d;j++)G[i][j]=AS.G[i][j];
724 #ifdef TEST_MIN_GLZ
725 for(i=0;i<*d;i++){printf(" G[%d]=",i);for(j=0;j<*d;j++)printf(" %5d",
726 G[i][j]);if(!i)printf(" g=%d norm=%d",(int)AS.g,(int)AS.s);puts("");}
727 #endif
728 return AS.g;
729 }
CWS_to_PermCWS(CWS * Cin,CWS * C,int * pi)730 void CWS_to_PermCWS(CWS *Cin, CWS *C, int *pi)
731 { int i, j, N=C->N=Cin->N, n=C->nw=Cin->nw, l=0, A[AMBI_Dmax];
732 Long *X, *G[AMBI_Dmax], M[AMBI_Dmax][AMBI_Dmax], W[AMBI_Dmax];
733 for(j=1;j<n;j++) if(Cin->d[j]>Cin->d[l]) l=j; X=Cin->W[l];
734 j=0; for(i=0;i<N;i++) if(X[i]) {A[j]=i; W[j]=labs(X[i]); G[j]=M[j]; j++;}
735 Wperm_to_GLZ(W,&j,G,pi); for(i=0;i<j;i++) pi[i]=A[pi[i]];
736 for(i=0;i<N;i++) if(X[i]==0) pi[j++]=i;
737 C->d[0]=Cin->d[l]; for(i=0;i<N;i++) C->W[0][i]=Cin->W[l][pi[i]];
738 for(j=1;j<n;j++) {int L=j-(j<=l); C->d[j]=Cin->d[L];
739 for(i=0;i<N;i++) C->W[j][i]=Cin->W[L][pi[i]];}
740 C->nz=Cin->nz; for(j=0;j<C->nz;j++){C->m[j]=Cin->m[j];
741 for(i=0;i<N;i++) C->z[j][i]=Cin->z[j][pi[i]];}
742 }
743 #endif
744
745 /* ========== For WS for 5d-polytopes, Sep 2017 ========== */
int_ld(Long w)746 int int_ld(Long w){int i=-1; while (w) {w /= 2; i++;} return i;}
747
Initialize_C5S(C5stats * _C5S,int n)748 void Initialize_C5S(C5stats *_C5S, int n){
749 int k;
750 _C5S->n_nonIP = 0;
751 _C5S->n_IP_nonRef = 0;
752 _C5S->n_ref = 0;
753 _C5S->nr_max_mp = 0;
754 _C5S->nr_max_mv = 0;
755 _C5S->nr_max_nv = 0;
756 _C5S->nr_max_w = 0;
757 for (k=0; k<MAXLD; k++) {_C5S->n_w[k] = 0; _C5S->nr_n_w[k] = 0;}
758 _C5S->max_mp = 0;
759 _C5S->max_mv = 0;
760 _C5S->max_np = 0;
761 _C5S->max_nv = 0;
762 _C5S->max_h22 = 0;
763 _C5S->max_w = 0;
764 for (k=1; k<n-1; k++) _C5S->max_h1[k] = 0;
765 for (k=0; k<=n; k++) _C5S->max_nf[k] = 0;
766 _C5S->max_chi = -100000000;
767 _C5S->min_chi = 100000000;
768 }
769
Update_C5S(BaHo * _BH,int * nf,Long * W,C5stats * _C5S)770 void Update_C5S(BaHo *_BH, int *nf, Long *W, C5stats *_C5S){
771 if (_BH->np) { // reflexive case
772 int i, chi = 48+6*(_BH->h1[1]-_BH->h1[2]+_BH->h1[3]), ld = int_ld(W[5]);
773 assert(0<=ld); assert(ld<MAXLD);
774 if (_BH->mp > _C5S->max_mp) _C5S->max_mp = _BH->mp;
775 if (_BH->mv > _C5S->max_mv) _C5S->max_mv = _BH->mv;
776 if (_BH->np > _C5S->max_np) _C5S->max_np = _BH->np;
777 if (_BH->nv > _C5S->max_nv) _C5S->max_nv = _BH->nv;
778 if (_BH->h22 > _C5S->max_h22) _C5S->max_h22 = _BH->h22;
779 for (i=1; i<_BH->n-1; i++)
780 if (_BH->h1[i] > _C5S->max_h1[i]) _C5S->max_h1[i] = _BH->h1[i];
781 for (i=0; i<=_BH->n; i++)
782 if (nf[i] > _C5S->max_nf[i]) _C5S->max_nf[i] = nf[i];
783 if (chi > _C5S->max_chi) _C5S->max_chi = chi;
784 if (chi < _C5S->min_chi) _C5S->min_chi = chi;
785 _C5S->n_ref++;
786 _C5S->n_w[ld]++;
787 if(W[5] > _C5S->max_w) _C5S->max_w = W[5];}
788 else { // IP but non-reflexive case
789 if (_BH->mp > _C5S->nr_max_mp) _C5S->nr_max_mp = _BH->mp;
790 if (_BH->mv > _C5S->nr_max_mv) _C5S->nr_max_mv = _BH->mv;
791 if (_BH->nv > _C5S->nr_max_nv) _C5S->nr_max_nv = _BH->nv;
792 _C5S->n_IP_nonRef++;
793 _C5S->nr_n_w[int_ld(W[5])]++;
794 if(W[5] > _C5S->nr_max_w) _C5S->nr_max_w = W[5];}
795 }
796
Print_C5S(C5stats * _C5S)797 void Print_C5S(C5stats *_C5S){
798 int i;
799 printf("non-IP: #=%ld\n", _C5S->n_nonIP);
800 printf("IP, non-reflexive: #=%ld, max_mp=%d, max_mv=%d, max_nv=%d, max_w=%ld\n",
801 _C5S->n_IP_nonRef, _C5S->nr_max_mp, _C5S->nr_max_mv, _C5S->nr_max_nv,
802 _C5S->nr_max_w);
803 printf(" #(w5) of given ld: ");
804 for (i=0; i<MAXLD; i++) printf(" %d:%ld", i, _C5S->nr_n_w[i]);
805 puts("");
806 printf("reflexive: #=%ld, max_mp=%d, max_mv=%d, max_np=%d, max_nv=%d, max_w=%ld\n",
807 _C5S->n_ref, _C5S->max_mp, _C5S->max_mv, _C5S->max_np, _C5S->max_nv,
808 _C5S->max_w);
809 printf(" #(w5) of given ld: ");
810 for (i=0; i<MAXLD; i++) printf(" %d:%ld", i, _C5S->n_w[i]);
811 puts("");
812 printf(" max #(faces): %d %d %d %d %d\n", _C5S->max_nf[0], _C5S->max_nf[1],
813 _C5S->max_nf[2], _C5S->max_nf[3], _C5S->max_nf[4]);
814 printf(" h11<=%d, h12<=%d, h13<=%d, h22<=%d, %d<=chi<=%d\n",
815 _C5S->max_h1[1], _C5S->max_h1[2], _C5S->max_h1[3],
816 _C5S->max_h22, _C5S->min_chi, _C5S->max_chi);
817 }
818
819