1 #include "Global.h"
2 #include "Subpoly.h"
3 #include "Rat.h"
4
5 #define subl_int LLong
6 #define UnAided_IP_CHECK (POLY_Dmax>4)
7 #define SIMPLE_CTH (0)
8 #define INCOMPLETE_SL_REDUCTION (0)
9 #define TEST_Aided_IP_CHECK (0)
10 #ifndef CEQ_Nmax
11 #define CEQ_Nmax EQUA_Nmax
12 #endif
13
14 #define IMPROVE_SL_COORD (0) /* 0=no 1=SL(old) 2=GL(new) */
15 #define IMPROVE_SL_REGCD (0)
16
17 typedef struct {int nk, k[VERT_Nmax];} KeepList;
18 typedef struct {int ne; Equation e[CEQ_Nmax];} CEqList;
19
20 /* ========== Auxiliary routines from other modules ========== */
21
22 Long CompareEq(Equation *X, Equation *Y, int n);
23 int IsGoodCEq(Equation *_E, PolyPointList *_P, VertexNumList *_V);
24 int Finish_IP_Check(PolyPointList *_P, VertexNumList *_V, EqList *_F,
25 CEqList *_CEq, INCI *F_I, INCI *CEq_I); /* from Vertex.c */
26 int Improve_Coords(PolyPointList *_P,VertexNumList *_V); /* from Polynf.c */
27
New_Improve_Coords(PolyPointList * _P,VertexNumList * _V)28 int New_Improve_Coords(PolyPointList *_P,VertexNumList *_V)
29 { switch(IMPROVE_SL_COORD)
30 { case 0: return 1; case 1: return Improve_Coords(_P,_V);
31 case 2: Make_Poly_UTriang(_P); return 1; default: assert(0);}
32 }
33
34 /* ====================================================================== */
35 /* ========== ========== */
36 /* ========== S U B P O L Y S ========== */
37 /* ========== ========== */
38 /* ====================================================================== */
39
IsNewEq(Equation * E,CEqList * _C,EqList * _F,int * _n)40 int IsNewEq(Equation *E, CEqList *_C, EqList *_F, int *_n){
41 int i=_C->ne;
42 while(i--) if(!CompareEq(E,&_C->e[i],*_n)) return 0;
43 i=_F->ne;
44 while(i--) if(!CompareEq(E,&_F->e[i],*_n)) return 0;
45 return 1;
46 }
47
Irrel(INCI * _INCI,INCI * INCI_List,int * n_irrel)48 int Irrel(INCI *_INCI, INCI *INCI_List, int *n_irrel){
49 int i;
50 for (i=0;i<*n_irrel;i++) if (INCI_EQ(*_INCI, INCI_List[i])) return 1;
51 return 0;
52 }
53
INCI_To_VertexNumList(INCI X,VertexNumList * _V,int n)54 void INCI_To_VertexNumList(INCI X, VertexNumList *_V, int n){
55 _V->nv=0;
56 while(!INCI_EQ_0(X)) {n--; if(INCI_M2(X)) _V->v[_V->nv++]=n; X=INCI_D2(X);}
57 }
58
Remove_INCI_From_List(INCI New_INCI,INCI * INCI_List,int * n_INCI)59 void Remove_INCI_From_List(INCI New_INCI, INCI *INCI_List, int *n_INCI){
60 int i;
61 for (i=*n_INCI-1;i>=0;i--) {
62 if(INCI_LE(INCI_List[i],New_INCI)) INCI_List[i]=INCI_List[--(*n_INCI)];}
63 }
64
Add_INCI_To_List(INCI New_INCI,INCI * INCI_List,int * n_INCI)65 void Add_INCI_To_List(INCI New_INCI, INCI *INCI_List, int *n_INCI){
66 int i;
67 for (i=*n_INCI-1;i>=0;i--) {
68 if(INCI_LE(New_INCI,INCI_List[i])) return;
69 if(INCI_LE(INCI_List[i],New_INCI)) INCI_List[i]=INCI_List[--(*n_INCI)];}
70 assert(*n_INCI < CD2F_Nmax);
71 INCI_List[(*n_INCI)++]=New_INCI;
72 }
73 /*
74 An INCI_List corresponds to a list of facets of some polytope, which might
75 be the original polytope or one of its faces. This implies that no two
76 INCIs x,y of the INCI_List fulfill INCI_LE(x,y).
77 Consequently, Add_INCI_To_List enhances INCI_List by New_INCI only if there is
78 no y in INCI_List such that INCI_LE(New_INCI,y) is fulfilled;
79 at the same time it removes any y with INCI_LE(y,New_INCI) from INCI_List.
80 Remove_INCI_From_List removes any y with INCI_LE(y,New_INCI) from INCI_List.
81 */
82
83 INCI List_Complete(int n, int n_polys, int n_facets,
84 INCI *polys, INCI *facets);
85
Poly_Complete(int n,int * n_polyfacets,INCI * poly,INCI * polyfacets)86 INCI Poly_Complete(int n, int *n_polyfacets, INCI *poly, INCI *polyfacets){
87 int j, k, n_cd2_faces=0;
88 INCI cd2_faces[CD2F_Nmax];
89 INCI X=INCI_0();
90 for (j=0;j<*n_polyfacets;j++) X=INCI_OR(X,polyfacets[j]);
91 if (n==1) return INCI_XOR(X,*poly);
92 for (j=0;j<*n_polyfacets;j++) for (k=j+1;k<*n_polyfacets;k++)
93 Add_INCI_To_List(INCI_AND(polyfacets[j],polyfacets[k]),
94 cd2_faces,&n_cd2_faces);
95 return INCI_OR(INCI_XOR(X,*poly),
96 List_Complete(n-1,*n_polyfacets,n_cd2_faces,polyfacets,cd2_faces));
97 }
98
List_Complete(int n,int n_polys,int n_facets,INCI * polys,INCI * facets)99 INCI List_Complete(int n, int n_polys, int n_facets,
100 INCI *polys, INCI *facets){
101 int i, j, n_polyfacets;
102 INCI polyfacets[CEQ_Nmax];
103 INCI X=INCI_0();
104 for (i=0;i<n_polys;i++){
105 n_polyfacets=0;
106 for (j=0;j<n_facets;j++) if (INCI_LE(facets[j],polys[i]))
107 polyfacets[n_polyfacets++]=facets[j];
108 X=INCI_OR(X,Poly_Complete(n, &n_polyfacets, &(polys[i]), polyfacets));}
109 return X;
110 }
111
FE_Close_the_Hole(PolyPointList * _P,VertexNumList * _V,EqList * _E,CEqList * _CEq,int n_old_v,INCI * CEq_INCI,INCI Hole_Vert_INCI)112 void FE_Close_the_Hole(PolyPointList *_P, VertexNumList *_V, EqList *_E,
113 CEqList *_CEq, int n_old_v, INCI *CEq_INCI, INCI Hole_Vert_INCI){
114 static PolyPointList P;
115 VertexNumList Hole_Verts;
116 EqList BVE;
117 int i,j,n_new_Eq=0;
118 INCI_To_VertexNumList(Hole_Vert_INCI, &Hole_Verts, n_old_v);
119 if(Hole_Verts.nv<_P->n){
120 printf("Hole_Verts.nv=%d\n",Hole_Verts.nv);
121 puts("Hole_Verts:");
122 Print_INCI(Hole_Vert_INCI);
123 puts("\nCEq_INCI:");
124 for(i=0;i<_CEq->ne;i++) {Print_INCI(CEq_INCI[i]); puts("");}
125 Print_VL(_P,_V,"");
126 Print_EL((EqList *) _E,&_P->n,0,"");
127 Print_EL((EqList *) _CEq,&_P->n, 0,"");
128 exit(0);}
129 P.n=_P->n;
130 P.np=Hole_Verts.nv;
131 for (i=0;i<Hole_Verts.nv;i++) for (j=0;j<_P->n;j++)
132 P.x[i][j]=_P->x[Hole_Verts.v[i]][j];
133 Find_Equations(&P,&Hole_Verts,&BVE);
134 assert(BVE.ne);
135 for (i=0;i<BVE.ne;i++) if(IsGoodCEq(&BVE.e[i],_P,_V))
136 if(IsNewEq(&BVE.e[i],_CEq,_E,&_P->n)){
137 n_new_Eq++;
138 assert(_CEq->ne<CEQ_Nmax);
139 CEq_INCI[_CEq->ne]=Eq_To_INCI(&BVE.e[i],_P,_V);
140 _CEq->e[_CEq->ne++]=BVE.e[i];}
141 assert(n_new_Eq);
142 }
143
Close_the_Hole(PolyPointList * _P,VertexNumList * _V,EqList * _E,CEqList * _CEq,int old_ne,int n_old_v,int n_Hole_Faces,INCI * E_INCI,INCI * CEq_INCI,INCI Hole_Verts,INCI * Hole_Faces)144 void Close_the_Hole(PolyPointList *_P, VertexNumList *_V, EqList *_E,
145 CEqList *_CEq, int old_ne, int n_old_v, int n_Hole_Faces,
146 INCI *E_INCI, INCI *CEq_INCI, INCI Hole_Verts, INCI *Hole_Faces){
147
148 INCI Bad_Vert_INCI, cd2_Faces[CD2F_Nmax];
149 int n_cd2_Faces=n_Hole_Faces, np=_P->np, i, j;
150
151 /* puts("CTH");
152 Print_PPL(_P);
153 fprintf(outFILE,"_E:\n");
154 for (i=0;i<_E->ne;i++){
155 for (j=0;j<_P->n;j++) fprintf(outFILE,"%d ", _E->e[i].a[j]);
156 fprintf(outFILE," %d\n", _E->e[i].c);}
157 fprintf(outFILE,"n_irrel=%d, n_old_v=%d\n", n_irrel, n_old_v);
158 fprintf(outFILE,"E_INCI: ");
159 for (i=0;i<_E->ne;i++) Print_INCI(E_INCI[i]);
160 fprintf(outFILE,"\nHole_Faces: ");
161 for (i=0;i<n_Hole_Faces;i++) Print_INCI(Hole_Faces[i]);
162 fprintf(outFILE,"\n");
163 fflush(0); */
164
165 for (i=0;i<n_cd2_Faces;i++) cd2_Faces[i]=Hole_Faces[i];
166 _P->np=n_old_v;
167 while(n_Hole_Faces--){
168 int not_found=1;
169 VertexNumList Other_Verts;
170 Equation Eq, E1, E2;
171 /* Connect the codim 2 faces along the hole with other vertices
172 bounding the hole; calculate corresponding equations and
173 check if they are good and new */
174
175 for (j=0;(j<_E->ne)&¬_found;j++)
176 if (INCI_LE(Hole_Faces[n_Hole_Faces],E_INCI[j])){
177 E1=_E->e[j]; not_found=0;}
178 if (not_found)
179 fprintf(outFILE,"Equation1 not found in Close_The_Hole!\n");
180 not_found=1;
181 for (j=_E->ne;(j<old_ne)&¬_found;j++)
182 if (INCI_LE(Hole_Faces[n_Hole_Faces],E_INCI[j])){
183 E2=_E->e[j]; not_found=0;}
184 if (not_found)
185 fprintf(outFILE,"Equation2 not found in Close_The_Hole!\n");
186 INCI_To_VertexNumList(INCI_XOR(Hole_Faces[n_Hole_Faces],Hole_Verts),
187 &Other_Verts, n_old_v);
188 for (i=0;i<Other_Verts.nv;i++){
189 Eq=EEV_To_Equation(&E1,&E2,_P->x[Other_Verts.v[i]],_P->n);
190 if (IsGoodCEq(&Eq,_P,_V)) {
191 CEq_INCI[_CEq->ne]=Eq_To_INCI(&Eq,_P,_V);
192 if (!Irrel(&(CEq_INCI[_CEq->ne]),E_INCI,&_E->ne)){
193 assert(_CEq->ne<CEQ_Nmax);
194 _CEq->e[_CEq->ne]=Eq;
195 Remove_INCI_From_List(CEq_INCI[_CEq->ne],
196 Hole_Faces,&n_Hole_Faces);
197 _CEq->ne++;
198 break;}}}}
199
200 /* for (i=0;i<_CEq->ne;i++){
201 fprintf(outFILE,"_CEq->e[%d]: ",i);
202 for (j=0;j<_P->n;j++) fprintf(outFILE," %d",_CEq->e[i].a[j]);
203 fprintf(outFILE," %d\n",_CEq->e[i].c);}
204 for (i=0;i<_E->ne;i++){
205 fprintf(outFILE,"_E->e[%d]: ",i);
206 for (j=0;j<_P->n;j++) fprintf(outFILE," %d",_E->e[i].a[j]);
207 fprintf(outFILE," %d\n",_E->e[i].c);}
208 fflush(0); */
209
210 for (i=0;i<_CEq->ne;i++) for (j=0;j<i;j++) {
211 INCI New_Face=INCI_AND(CEq_INCI[i],CEq_INCI[j]);
212 int k;
213 if (INCI_abs(New_Face)<_P->n-1) continue;
214 for (k=0;k<_CEq->ne;k++) if (INCI_LE(New_Face,CEq_INCI[k]))
215 if((k!=i)&&(k!=j)) break;
216 if (k<_CEq->ne) continue;
217 for (k=0;k<_E->ne;k++) if (INCI_LE(New_Face,E_INCI[k])) break;
218 if (k<_E->ne) continue;
219 assert(n_cd2_Faces<CD2F_Nmax);
220 cd2_Faces[n_cd2_Faces++]=New_Face; }
221 Bad_Vert_INCI=List_Complete(_P->n-1, _CEq->ne, n_cd2_Faces,
222 CEq_INCI, cd2_Faces);
223 if (!INCI_EQ_0(Bad_Vert_INCI))
224 FE_Close_the_Hole(_P, _V, _E, _CEq, n_old_v, CEq_INCI, Bad_Vert_INCI);
225 _P->np=np;
226 }
227
Aided_IP_Check(PolyPointList * _P,VertexNumList * _V,EqList * _E,int n_irrel,int n_old_v)228 int Aided_IP_Check(PolyPointList *_P, VertexNumList *_V, EqList *_E,
229 int n_irrel, int n_old_v){
230
231 /* The first n_old_v entries of _P are the old vertices;
232 the first n_irrel entries of _E are the irrelevant old facets (which
233 remain facets), the remaining facets are the relevant facets */
234
235 int i, j, n_Hole_Faces=0, old_ne=_E->ne;
236 INCI CEq_INCI[CEQ_Nmax], Hole_Verts;
237 INCI E_INCI[EQUA_Nmax];
238 INCI Hole_Faces[CD2F_Nmax];
239 CEqList CEq;
240
241 if (n_irrel==0) {
242 fprintf(outFILE, "n_irrel=0 in Aided_IP_Check!"); exit(0); }
243
244 /* Create E_INCI: Incidences between _E and old vertices;
245 Intersect relevant and irrelevant facets to get the codim 2 faces
246 that bound the hole;
247 Create Hole_Verts as the union of Hole_Faces */
248
249 for (i=0;i<_E->ne;i++){
250 E_INCI[i]=INCI_0();
251 for (j=0;j<n_old_v;j++) E_INCI[i]=
252 INCI_PN(E_INCI[i],Eval_Eq_on_V(&(_E->e[i]),_P->x[j],_P->n));}
253
254 for (i=0;i<n_irrel;i++) for (j=n_irrel;j<_E->ne;j++) {
255 INCI New_Face=INCI_AND(E_INCI[i],E_INCI[j]);
256 int k;
257 if (INCI_abs(New_Face)<_P->n-1) continue;
258 for (k=0;k<_E->ne;k++) if (INCI_LE(New_Face,E_INCI[k]))
259 if((k!=i)&&(k!=j)) break;
260 if (k<_E->ne) continue;
261 assert(n_Hole_Faces<CD2F_Nmax);
262 Hole_Faces[n_Hole_Faces++]=New_Face;}
263
264 /* puts("AIP:");
265 Print_PPL(_P);
266 Print_EL(_E,&_P->n,0);
267 puts("E_INCI:");
268 for(i=0;i<_E->ne;i++) {Print_INCI(E_INCI[i]); puts("");}
269 printf("n_irrel=%d\n", n_irrel); */
270
271 Hole_Verts=INCI_0();
272 for (i=n_irrel;i<_E->ne;i++) Hole_Verts=INCI_OR(Hole_Verts, E_INCI[i]);
273 for (i=0;i<n_old_v;i++) _V->v[i]=i;
274 _V->nv=n_old_v;
275 _E->ne=n_irrel;
276
277 if (n_irrel==1) {
278 /* Find vertex of maximal distance from the irrelevant facet and make
279 Candidate equations by connecting it with the Hole_Faces */
280 Long maxdist;
281 int newvert=n_old_v;
282 if (n_old_v==_P->np) return 0;
283 maxdist=Eval_Eq_on_V(&(_E->e[0]),_P->x[n_old_v],_P->n);
284 for (i=n_old_v+1; i<_P->np; i++){
285 Long dist=Eval_Eq_on_V(&(_E->e[0]),_P->x[i],_P->n);
286 if (dist>maxdist) {
287 newvert=i;
288 maxdist=dist;
289 continue;}
290 if (dist==maxdist) if (Vec_Greater_Than(_P->x[i],_P->x[newvert],_P->n)){
291 newvert=i;
292 maxdist=dist;}}
293 _V->v[_V->nv++]=newvert;
294
295 for (i=0;i<n_Hole_Faces;i++){
296 Long dist;
297 int not_found=1;
298 Equation E1, E2;
299 for (j=0;(j<n_irrel)&¬_found;j++)
300 if (INCI_LE(Hole_Faces[i],E_INCI[j])){ E1=_E->e[j]; not_found=0;}
301 if (not_found)
302 fprintf(outFILE,"Equation1 not found in Aided_IP_Check!\n");
303 not_found=1;
304 for (j=n_irrel;(j<old_ne)&¬_found;j++)
305 if (INCI_LE(Hole_Faces[i],E_INCI[j])){ E2=_E->e[j]; not_found=0;}
306 if (not_found)
307 fprintf(outFILE,"Equation2 not found in Aided_IP_Check!\n");
308 CEq.e[i]=EEV_To_Equation(&E1,&E2,_P->x[newvert],_P->n);
309 for (j=0;j<n_old_v;j++)
310 if((dist=Eval_Eq_on_V(&(CEq.e[i]),_P->x[j],_P->n))){
311 if (dist<0){
312 int k;
313 CEq.e[i].c=-CEq.e[i].c;
314 for(k=0;k<_P->n;k++) CEq.e[i].a[k]=-CEq.e[i].a[k]; }
315 break;}}
316
317 E_INCI[0]=INCI_PN(E_INCI[0],1);
318 for (i=0;i<n_Hole_Faces;i++) CEq_INCI[i]=INCI_PN(Hole_Faces[i],0);
319
320 CEq.ne=n_Hole_Faces; }
321
322 else {
323 CEq.ne=0;
324 #if SIMPLE_CTH
325 FE_Close_the_Hole(_P,_V,_E,&CEq,n_old_v,CEq_INCI,Hole_Verts);}
326 #else
327 Close_the_Hole(_P,_V,_E,&CEq,old_ne,n_old_v,n_Hole_Faces,
328 E_INCI, CEq_INCI, Hole_Verts, Hole_Faces);}
329 #endif
330
331 /*Print_PPL(_P);
332 Print_VL(_P,_V,"Verts:");
333 for (i=0;i<CEq.ne;i++){
334 fprintf(outFILE,"CEq.e[%d]: ",i);
335 for (j=0;j<_P->n;j++) fprintf(outFILE," %d",CEq.e[i].a[j]);
336 fprintf(outFILE," %d\n",CEq.e[i].c);}
337 fprintf(outFILE,"CEq_INCI: ");
338 for (i=0;i<CEq.ne;i++) Print_INCI(CEq_INCI[i]);
339 fprintf(outFILE,"\n");
340 for (i=0;i<_E->ne;i++){
341 fprintf(outFILE,"_E->e[%d]: ",i);
342 for (j=0;j<_P->n;j++) fprintf(outFILE," %d",_E->e[i].a[j]);
343 fprintf(outFILE," %d\n",_E->e[i].c);}
344 fprintf(outFILE,"E_INCI: ");
345 for (i=0;i<_E->ne;i++) Print_INCI(E_INCI[i]);
346 fflush(0); */
347
348 return Finish_IP_Check(_P, _V, _E, &CEq, E_INCI, CEq_INCI);
349 }
350
351 void Make_All_Subpolys(PolyPointList *_P, EqList *_E, VertexNumList *_V,
352 KeepList *_KL, NF_List *_NFL);
353
Relevant(Long * X,int * n,EqList * _E,int * n_irrel)354 int Relevant(Long *X, int *n, EqList *_E, int *n_irrel){
355 int i;
356 for (i=0;i<*n_irrel;i++)
357 if (!Eval_Eq_on_V(&(_E->e[i]),X,*n)) return 0;
358 for (i=0;i<*n;i++) if (X[i]) return 1;
359 return 0;
360 }
361
kept(int i,KeepList * _KL)362 int kept(int i, KeepList *_KL){
363 int j;
364 for(j=0;j<_KL->nk;j++) if(_KL->k[j]==i) return j;
365 return -1;
366 }
367
Drop_and_Keep(PolyPointList * _P,VertexNumList * _V,EqList * _E,KeepList * _KL,int drop_num,NF_List * _NFL)368 void Drop_and_Keep(PolyPointList *_P, VertexNumList *_V, EqList *_E,
369 KeepList *_KL, int drop_num, NF_List *_NFL){
370
371 int i, j, n_irrel=0, IP;
372 int *new2old =(int *) malloc(sizeof(int)*POINT_Nmax);
373 Long drop_point[POLY_Dmax];
374 VertexNumList red_V;
375 EqList *new_E=(EqList *) malloc(sizeof(EqList));
376 PolyPointList *red_P = (PolyPointList *) malloc(sizeof (PolyPointList));
377
378 if(red_P==NULL) {puts("Unable to allocate space for red_P"); exit(0);}
379 if(new2old==NULL) {puts("Unable to allocate space for new2old"); exit(0);}
380
381 for (j=0;j<_P->n;j++) drop_point[j]=_P->x[_V->v[drop_num]][j];
382
383 /* Create new_E: Same as *_E, but irrelevant facets first */
384
385 new_E->ne=_E->ne;
386 for (i=0;i<_E->ne;i++){
387 if (Eval_Eq_on_V(&(_E->e[i]),drop_point,_P->n)){
388 for (j=0;j<_P->n;j++) new_E->e[n_irrel].a[j]=_E->e[i].a[j];
389 new_E->e[n_irrel++].c=_E->e[i].c;}
390 else {
391 for (j=0;j<_P->n;j++) new_E->e[_E->ne+n_irrel-i-1].a[j]=_E->e[i].a[j];
392 new_E->e[_E->ne+n_irrel-i-1].c=_E->e[i].c;}}
393
394 /* Create red_P: old vertices, points that are not on irrelevant facets */
395 red_P->np=0;
396 red_P->n=_P->n;
397 for (i=0;i<_V->nv;i++) if (i!=drop_num){
398 for (j=0;j<_P->n;j++) red_P->x[red_P->np][j]=_P->x[_V->v[i]][j];
399 new2old[red_P->np++]=_V->v[i];}
400 for (i=0;i<_P->np;i++)
401 if ((i!=_V->v[drop_num])&&Relevant(_P->x[i],&_P->n,new_E,&n_irrel)){
402 for (j=0;j<_P->n;j++) red_P->x[red_P->np][j]=_P->x[i][j];
403 new2old[red_P->np++]=i;}
404
405 /* Print_PPL(_P);
406 fprintf(outFILE,"new_E:\n");
407 for (i=0;i<new_E->ne;i++){
408 for (j=0;j<_P->n;j++) fprintf(outFILE,"%d ", new_E->e[i].a[j]);
409 fprintf(outFILE," %d\n", new_E->e[i].c);}
410 fprintf(outFILE,"\n");
411 Print_PPL(red_P);
412 fprintf(outFILE,"\n");
413 printf("drop_point: ");
414 for (i=0;i<_P->n;i++) printf("%d ", (int) drop_point[i]);
415 puts("");*/
416
417 _NFL->nIP++;
418
419 #if UnAided_IP_CHECK
420 IP=IP_Check(red_P,&red_V,new_E);
421 #else
422 IP=Aided_IP_Check(red_P,&red_V,new_E,n_irrel,_V->nv-1);
423 #endif
424
425 /* puts("After AIP (in D&K):");
426 Print_VL(red_P,&red_V);
427 Print_EL(new_E,&red_P->n,0); */
428
429 free(red_P);
430
431 if(IP){
432
433 VertexNumList new_V;
434 #if TEST_Aided_IP_CHECK
435 VertexNumList test_V;
436 EqList test_E;
437 #endif
438
439 /* Create new_V from red_V */
440 new_V.nv=red_V.nv;
441 for (i=0;i<red_V.nv;i++) if((new_V.v[i]=new2old[red_V.v[i]])==_P->np-1)
442 new_V.v[i]=_V->v[drop_num];
443
444 /* Drop the vertex _V->v[drop_num]: */
445 _P->np--;
446 for (j=0;j<_P->n;j++) _P->x[_V->v[drop_num]][j]=_P->x[_P->np][j];
447 j=kept(_P->np,_KL);
448 if (j>=0) _KL->k[j]=_V->v[drop_num];
449
450 #if TEST_Aided_IP_CHECK
451 if(!IP_Check(_P,&test_V,&test_E)||(test_V.nv!=red_V.nv)||
452 (test_E.ne!=new_E->ne)){
453 int k;
454 fprintf(outFILE,"_V: ");
455 for (i=0;i<_V->nv;i++) fprintf(outFILE,"%d ",_V->v[i]);
456 fprintf(outFILE,"\n");
457 fprintf(outFILE,"_E:\n");
458 for (i=0;i<_E->ne;i++){
459 for (k=0;k<_P->n;k++) fprintf(outFILE,"%d ", (int) _E->e[i].a[k]);
460 fprintf(outFILE," %d\n", (int) _E->e[i].c);}
461 fprintf(outFILE,"\n");
462 fprintf(outFILE,"drop_num: %d\n",drop_num);
463 fprintf(outFILE,"drop_point: ");
464 for (i=0;i<_P->n;i++) fprintf(outFILE,"%d ", (int) drop_point[i]);
465 fprintf(outFILE,"\n");
466 fprintf(outFILE,"_KL: ");
467 for (i=0;i<_KL->nk;i++) fprintf(outFILE,"%d ", _KL->k[i]);
468 fprintf(outFILE,"\n");
469
470 fprintf(outFILE,"red_V: ");
471 for (i=0;i<red_V.nv;i++) fprintf(outFILE,"%d ",red_V.v[i]);
472 fprintf(outFILE,"\n");
473 fprintf(outFILE,"new_E:\n");
474 for (i=0;i<new_E->ne;i++){
475 for (k=0;k<_P->n;k++) fprintf(outFILE,"%d ", (int) new_E->e[i].a[k]);
476 fprintf(outFILE," %d\n", new_E->e[i].c);}
477 fprintf(outFILE,"\n");
478 fprintf(outFILE,"n_irrel: %d\n",n_irrel);
479 fprintf(outFILE,"new2old:");
480 for (i=0;i<red_P->np;i++) fprintf(outFILE," %d", new2old[i]);
481 fprintf(outFILE,"\n");
482
483 fprintf(outFILE,"test_V: ");
484 for (i=0;i<test_V.nv;i++) fprintf(outFILE,"%d ",test_V.v[i]);
485 fprintf(outFILE,"\n");
486 fprintf(outFILE,"test_E:\n");
487 for (i=0;i<test_E.ne;i++){
488 for (k=0;k<_P->n;k++) fprintf(outFILE,"%d ", (int) test_E.e[i].a[k]);
489 fprintf(outFILE," %d\n", (int) test_E.e[i].c);}
490 fprintf(outFILE,"\n");
491
492 exit(0);}
493 #endif
494
495 Make_All_Subpolys(_P, new_E, &new_V, _KL, _NFL);
496
497 /* Reconstruct _P */
498 if (j>=0) _KL->k[j]=_P->np;
499 for (j=0;j<_P->n;j++){
500 _P->x[_P->np][j]=_P->x[_V->v[drop_num]][j];
501 _P->x[_V->v[drop_num]][j]=drop_point[j];}
502 _P->np++;}
503
504 /* Attach a keep label to the vertex _V->v[drop_num]: */
505 _KL->k[_KL->nk]=_V->v[drop_num];
506 _KL->nk++;
507 free(new2old); free(new_E);
508 }
509
510
511
Start_Make_All_Subpolys(PolyPointList * _P,NF_List * _NFL)512 void Start_Make_All_Subpolys(PolyPointList *_P, NF_List *_NFL){
513 int i,j;
514 VertexNumList V, new_V;
515 EqList E, new_E;
516 KeepList KL;
517 Long drop_point[POLY_Dmax];
518 _NFL->nIP++;
519
520 if (_NFL->rf) _NFL->rf=_NFL->rd;
521 _NFL->rd=0;
522
523 if(!IP_Check(_P,&V,&E)) {
524 fprintf(outFILE,"IP_check negative!\n"); return;}
525
526 if(_NFL->Nmin==0) {PairMat PM; Make_VEPM(_P,&V,&E,PM); /* complete poly if */
527 Complete_Poly(PM, &E, V.nv, _P); _NFL->Nmin=_P->np;} /* non-CWS input */
528
529 /* fprintf(outFILE,"StartMASP:\n");
530 for (i=0;i<V.nv;i++){
531 for (j=0;j<_P->n;j++) fprintf(outFILE,"%d ",(int) _P->x[V.v[i]][j]);
532 fprintf(outFILE,"\n");}
533 for (i=0;i<E.ne;i++){
534 for (j=0;j<_P->n;j++) fprintf(outFILE,"%d ", (int) E.e[i].a[j]);
535 fprintf(outFILE," %d\n", (int) E.e[i].c);}
536 fflush(0);*/
537
538 KL.nk=0;
539
540 if(_NFL->kf){
541 fprintf(outFILE,"The vertices are:\n");
542 for (i=0;i<V.nv;i++){
543 fprintf(outFILE,"%d. ",i);
544 for (j=0;j<_P->n;j++) fprintf(outFILE,"%d ",(int) _P->x[V.v[i]][j]);
545 fprintf(outFILE,"\n"); }
546 fprintf(outFILE,"How many of them do you want to keep?\n");
547 fscanf(inFILE,"%d",&(KL.nk));
548 fprintf(outFILE,"Which %d of them do you want to keep?\n",KL.nk);
549 for (i=0; i<KL.nk;i++){
550 fscanf(inFILE,"%d",&j);
551 KL.k[i]=V.v[j]; }
552 fprintf(outFILE,"Keeping\n");
553 for (i=0; i<KL.nk;i++) {
554 for(j=0;j<_P->n;j++) fprintf(outFILE,"%d ",(int) _P->x[KL.k[i]][j]);
555 fprintf(outFILE,"\n"); }}
556
557 else for(i=0;i<V.nv;i++) {
558 /* Try dropping the vertex V.v[i]: */
559 _P->np--;
560 for (j=0;j<_P->n;j++){
561 drop_point[j]=_P->x[V.v[i]][j];
562 _P->x[V.v[i]][j]=_P->x[_P->np][j]; }
563 if (!IP_Check(_P,&new_V,&new_E)) KL.k[KL.nk++]=V.v[i];
564 /* Reconstruct _P */
565 for (j=0;j<_P->n;j++){
566 _P->x[_P->np][j]=_P->x[V.v[i]][j];
567 _P->x[V.v[i]][j]=drop_point[j];}
568 _P->np++; }
569
570 Make_All_Subpolys(_P, &E, &V, &KL, _NFL);
571 }
572
573
IREgcd(Long * vec_in,int * d,Long * vec_out)574 void IREgcd(Long *vec_in, int *d, Long *vec_out){
575 int i, j, n_rem_perm, n_perm=1, a, b, perm_j, inv_perm[POLY_Dmax];
576 Long perm_vec_in[POLY_Dmax], perm_vec_out[POLY_Dmax], abs_sum, max_abs_sum=0;
577 for (i=1;i<=*d;i++) n_perm*=i;
578 for (i=0;i<n_perm;i++){
579 b=i;
580 n_rem_perm=n_perm;
581 for (j=0;j<*d;j++) inv_perm[j]=-1;
582 for (j=*d;j>0;j--){
583 n_rem_perm/=j;
584 a=b/n_rem_perm;
585 b=b%n_rem_perm;
586 perm_j=*d;
587 while (a>=0){
588 perm_j--;
589 if (inv_perm[perm_j]<0) a--;}
590 perm_vec_in[j-1]=vec_in[perm_j];
591 inv_perm[perm_j]=j-1;}
592 REgcd(perm_vec_in, d, perm_vec_out);
593 abs_sum=0;
594 for (j=0;j<*d;j++) abs_sum+=abs(perm_vec_out[j]);
595 if ((i==0)||(abs_sum<max_abs_sum)){
596 max_abs_sum=abs_sum;
597 for (j=0;j<*d;j++) vec_out[j]=perm_vec_out[inv_perm[j]]; } }
598 }
599
600
Make_RedVec(Long n,Long * badFE,Long * RedVec)601 int Make_RedVec(Long n, Long *badFE, Long *RedVec){
602 int i=0, j;
603 Long RV[POLY_Dmax], FE[POLY_Dmax], dist=0;
604 for (j=0;j<n;j++) if (badFE[j]) FE[i++]=badFE[j];
605 if(i==1) RV[0]=1;
606 else IREgcd(FE,&i,RV);
607 for (j=0;j<i;j++) dist+=FE[j]*RV[j];
608 if ((dist!=1)&&(dist!=-1)) {
609 fprintf(outFILE,"FE: ");
610 for (j=0;j<i;j++) fprintf(outFILE,"%d ",(int) FE[j]);
611 fprintf(outFILE,"\n"); fprintf(outFILE,"RV: ");
612 for (j=0;j<i;j++) fprintf(outFILE,"%d ",(int) RV[j]);
613 fprintf(outFILE,"\n");
614 fprintf(outFILE,"dist= %d\n", (int) dist);
615 return 0;}
616 i=0;
617 for (j=0;j<n;j++) if (badFE[j]) RedVec[j]=-dist*RV[i++];
618 else RedVec[j]=0;
619 return 1;
620 }
621
Reduce_Poly(PolyPointList * _P,EqList * _E,KeepList * _KL,NF_List * _NFL,int badfacet)622 void Reduce_Poly(PolyPointList *_P, EqList *_E, KeepList *_KL,
623 NF_List *_NFL, int badfacet){
624
625 /* reduce _P to the sublattice where badfacet has distance 1;
626 call Make_All_Subpolys for the reduced polyhedron new_P; */
627
628 VertexNumList V;
629 PolyPointList *new_P= (PolyPointList *) malloc(sizeof (PolyPointList));
630 KeepList new_KL;
631 int i, j;
632 Long RedVec[POLY_Dmax];
633 Equation BE=_E->e[badfacet];
634 Rat VPM_checksum_old, VPM_checksum_new;
635
636 /*Generate new_P (in old coordinates!): */
637
638 new_P->np=0;
639 new_KL.nk=0;
640 new_P->n=_P->n;
641 for (i=0;i<_P->np;i++){
642 Long dist=0;
643 for (j=0;j<_P->n;j++) dist-=BE.a[j]*_P->x[i][j];
644 if (!(dist%BE.c)){
645 for (j=0;j<_P->n;j++)
646 new_P->x[new_P->np][j]=_P->x[i][j];
647 if (kept(i,_KL)>=0) new_KL.k[new_KL.nk++]=new_P->np;
648 new_P->np++; } }
649
650 /* Print_PPL(new_P);
651 fprintf(outFILE,"kept: ");
652 for (j=0;j<new_KL.nk;j++) fprintf(outFILE,"%d ",new_KL.k[j]);
653 fprintf(outFILE,"\n\n"); */
654
655 if(new_KL.nk!=_KL->nk) {free(new_P); return;}
656 _NFL->nIP++;
657
658 if(!IP_Check(new_P,&V,_E)) {free(new_P); return;}
659
660 VPM_checksum_old=rI(0);
661 for (i=0;i<_E->ne;i++) {
662 Long s=0;
663 for (j=0;j<V.nv;j++)
664 s+=Eval_Eq_on_V(&(_E->e[i]),new_P->x[V.v[j]], new_P->n);
665 VPM_checksum_old=rS(VPM_checksum_old, rQ(rI(s),rI(_E->e[i].c)));}
666
667 /* Generate RedVec: */
668 if (!Make_RedVec(_P->n, BE.a, RedVec)){
669 Print_PPL(_P,"");
670 fprintf(outFILE,"Bad Facet: ");
671 for (j=0;j<_P->n;j++) fprintf(outFILE,"%d ",(int) BE.a[j]);
672 fprintf(outFILE," %d\n",(int) BE.c);
673 fprintf(outFILE,"kept: ");
674 for (j=0;j<_KL->nk;j++) fprintf(outFILE,"%d ",_KL->k[j]);
675 fprintf(outFILE,"\n");
676 fprintf(outFILE,"RedVec: ");
677 for (j=0;j<_P->n;j++) fprintf(outFILE,"%d ",(int) RedVec[j]);
678 fprintf(outFILE,"\n");
679 exit(0);}
680
681 /*Coordinates of new_P: */
682
683 for (i=0;i<new_P->np;i++){
684 Long dist=0;
685 for (j=0;j<new_P->n;j++) dist-=BE.a[j]*new_P->x[i][j];
686 if (!(dist%BE.c)){
687 for (j=0;j<new_P->n;j++)
688 new_P->x[i][j]-=(dist-dist/BE.c)*RedVec[j]; } }
689
690 /* Print_PPL(new_P);
691 fprintf(outFILE,"\n\n"); */
692
693 if(!New_Improve_Coords(new_P,&V)) {
694 Print_PPL(_P,"");
695 fprintf(outFILE,"Bad Facet: ");
696 for (j=0;j<_P->n;j++) fprintf(outFILE,"%d ",(int) BE.a[j]);
697 fprintf(outFILE," %d\n",(int) BE.c);
698 fprintf(outFILE,"kept: ");
699 for (j=0;j<_KL->nk;j++) fprintf(outFILE,"%d ",_KL->k[j]);
700 fprintf(outFILE,"\n");
701 fprintf(outFILE,"RedVec: ");
702 for (j=0;j<_P->n;j++) fprintf(outFILE,"%d ",(int) RedVec[j]);
703 fprintf(outFILE,"\n");
704 puts("WARNING: INCOMPLETE SL REDUCTION"); fflush(stdout);
705 #if INCOMPLETE_SL_REDUCTION
706 free(new_P); return;
707 #else
708 exit(0);
709 #endif
710 }
711
712 if(!IP_Check(new_P,&V,_E)) {
713 fprintf(outFILE,"Trouble in Reduce_Poly!"); exit(0);}
714
715 VPM_checksum_new=rI(0);
716 for (i=0;i<_E->ne;i++) {
717 Long s=0;
718 for (j=0;j<V.nv;j++)
719 s+=Eval_Eq_on_V(&(_E->e[i]),new_P->x[V.v[j]],new_P->n);
720 VPM_checksum_new=rS(VPM_checksum_new, rQ(rI(s),rI(_E->e[i].c)));}
721 if(VPM_checksum_new.D*VPM_checksum_old.N-
722 VPM_checksum_new.N*VPM_checksum_old.D) {
723 fprintf(outFILE,"Checksums don't match in Reduce_Poly!"); exit(0); }
724
725 Make_All_Subpolys(new_P,_E,&V,&new_KL, _NFL);
726 free(new_P);
727 }
728
729
Make_All_Subpolys(PolyPointList * _P,EqList * _E,VertexNumList * _V,KeepList * _KL,NF_List * _NFL)730 void Make_All_Subpolys(PolyPointList *_P, EqList *_E, VertexNumList *_V,
731 KeepList *_KL, NF_List *_NFL){
732
733 /* Creates all subpolyhedra of an IP-polyhedron given by _P, _E, _V.
734 The basic structure is:
735 Search for bad facet;
736 If (reflexive) if (!Add_NF_to_List) return;
737 Drop_and_Keep vertices (on bad facet / everywhere);
738 If (bad facet exists) Reduce_Poly; */
739
740 int i, j, badfacet=-1, nbadvert=0, maxnkept=-1;
741 KeepList new_KL=*_KL;
742
743 /* fprintf(outFILE,"MASP:\n");
744 Print_PPL(_P);
745 Print_VL(_P,_V);
746 Print_EL(_E,&_P->n,0);
747 fflush(0);*/
748
749 if (_V->nv>_NFL->VN) _NFL->VN=_V->nv;
750 if (_E->ne>_NFL->FN) _NFL->FN=_E->ne;
751
752 /* Choose bad facet with largest number of kept vertices: */
753
754 for (i=0;i<_E->ne;i++) if (_E->e[i].c>1){
755 int nkept=0;
756 for (j=0;j<_KL->nk;j++)
757 if (!Eval_Eq_on_V(&(_E->e[i]),_P->x[_KL->k[j]],_P->n)) nkept++;
758 if (nkept>maxnkept){
759 maxnkept=nkept;
760 badfacet=i; }}
761
762 if (!_NFL->rf) if (badfacet==-1){ /* not recover mode && reflexive */
763 if (_NFL->of<=-2) { /* flags -o0 (break on reflexive) */
764 if (_NFL->rd) { /* or -oc (recover missing) */
765 if(_NFL->of==-2){Add_NF_to_List(_P,_V,_E,_NFL); return;} /* flag -o0 */
766 else {if(!Add_NF_to_List(_P,_V,_E,_NFL)) return;} /* flag -oc */
767 }}
768 else if(!Add_NF_to_List(_P,_V,_E,_NFL)) return;} /* standard */
769
770 if(0 < _NFL->of) if(_NFL->of <= _NFL->rd) return; /* limit rd */
771
772 if (badfacet!=-1) {
773 /* Find vertices on bad facet: */
774 nbadvert=0;
775 for (i=0;i<_V->nv;i++)
776 if (!Eval_Eq_on_V(&(_E->e[badfacet]),_P->x[_V->v[i]],_P->n)){
777 j=_V->v[nbadvert];
778 _V->v[nbadvert]=_V->v[i];
779 _V->v[i]=j;
780 nbadvert++; } }
781
782 if (_NFL->rf&&(_NFL->rf<_NFL->rd)) fprintf(outFILE,"_NFL->rf<_NFL->rd\n");
783 if (_NFL->rf==_NFL->rd) _NFL->rf=0;
784 if (_NFL->rf) for (i=0;i<_NFL->b[_NFL->rd];i++)
785 if((kept(_V->v[i],&new_KL)<0)) new_KL.k[new_KL.nk++]=_V->v[i];
786 _NFL->rd++;
787
788 for(i=0;i<((badfacet==-1)?_V->nv:nbadvert);i++)
789 if((kept(_V->v[i],&new_KL)<0)){
790 _NFL->b[_NFL->rd-1]=i;
791 Drop_and_Keep(_P, _V, _E, &new_KL, i, _NFL); }
792
793 /* If a bad facet is kept, reduce the poly: */
794 if((_NFL->of==0)||(_NFL->of==-3)) if (badfacet!=-1) {/* -oc recover[of=-3] */
795 int currentSL=_NFL->SL;
796 _NFL->SL=1;
797 _NFL->b[_NFL->rd-1]=nbadvert;
798 Reduce_Poly(_P, _E, &new_KL, _NFL, badfacet);
799 _NFL->SL=currentSL; }
800
801 _NFL->rd--;
802 }
803
Ascii_to_Binary(CWS * W,PolyPointList * P,char * dbin,char * polyi,char * polyo)804 void Ascii_to_Binary(CWS *W, PolyPointList *P,
805 char *dbin, char *polyi, char *polyo){
806 NF_List *_NFL=(NF_List *) malloc(sizeof(NF_List));
807 VertexNumList V;
808 EqList F;
809 assert(_NFL!=NULL);
810 if(!(*polyo)) {
811 puts("You have to specify an output file via -po in -a-mode!\n");
812 printf("For more help type use option `-h'\n");
813 exit(0);}
814 _NFL->of=0; _NFL->rf=0;
815 _NFL->iname=polyi; _NFL->oname=polyo; _NFL->dbname=dbin;
816 Init_NF_List(_NFL);
817 _NFL->SL=0;
818
819 while(Read_CWS_PP(W,P)) {
820 _NFL->hc = 0;
821 _NFL->V= _NFL->F= _NFL->VN= _NFL->FN = _NFL->Xnuc = _NFL->Xdif = 0;
822 _NFL->Nmin=P->np; _NFL->Nmax=0;
823 if(_NFL->d==0) _NFL->d=P->n;
824 else if(_NFL->d-P->n) {puts("different dim!"); exit(0);}
825 if (!IP_Check(P,&V,&F)){
826 puts("IP_Check failed in Ascii_to_Binary!\n"); exit(0);}
827 if (Add_NF_to_List(P,&V,&F,_NFL)) if (outFILE!=stdout){
828 int i, j;
829 for(i=0;i<W->nw;i++) {
830 fprintf(outFILE,"%d ",(int) W->d[i]);
831 for(j=0;j<W->N;j++) fprintf(outFILE,"%d ",(int) W->W[i][j]);
832 if(i+1<W->nw) fprintf(outFILE," "); else fprintf(outFILE,"\n"); }
833 fflush(0);} }
834 Write_List_2_File(polyo,_NFL);
835 free(_NFL);
836 }
837
Do_the_Classification(CWS * W,PolyPointList * P,int oFlag,int rFlag,int kFlag,char * polyi,char * polyo,char * dbin)838 void Do_the_Classification(CWS *W, PolyPointList *P, /* char *fn, */
839 int oFlag, int rFlag, int kFlag, char *polyi, char *polyo, char *dbin) {
840
841 /* static int nw; */
842 NF_List *_NFL=(NF_List *) malloc(sizeof(NF_List));
843 time_t W_SAVE_TIME=time(NULL);
844 if(!(*polyo)) {
845 puts("You have to specify an output file via -po!\n");
846 printf("For more help use option '-h'\n");
847 exit(0);}
848 assert(_NFL!=NULL); _NFL->of=oFlag; _NFL->rf=rFlag; _NFL->kf=kFlag;
849 _NFL->iname=polyi; _NFL->oname=polyo; _NFL->dbname=dbin;
850 Init_NF_List(_NFL);
851 rFlag=0; /* now used as "read flag" */
852
853 while(Read_CWS_PP(W,P)) { /* make subpolys */
854 if(W->nw>0) _NFL->Nmin=P->np; else _NFL->Nmin=0 /* :: Complete_Poly() */ ;
855 _NFL->Nmax=0; _NFL->SL = _NFL->hc = 0;
856 _NFL->V= _NFL->F= _NFL->VN= _NFL->FN = _NFL->Xnuc = _NFL->Xdif = 0;
857 if(_NFL->d==0) _NFL->d=P->n;
858 else if(_NFL->d-P->n) {puts("different dim!"); exit(0);}
859 if(rFlag) {Read_File_2_List(polyo,_NFL); rFlag=0;}
860 Start_Make_All_Subpolys(P, _NFL);
861 Print_Weight_Info(W,_NFL);
862 if((WRITE_DIM <= P->n) && (MIN_NEW <= _NFL->NP))
863 if((int)difftime(time(NULL),W_SAVE_TIME) > MIN_W_SAVE_TIME) {
864 Write_List_2_File(polyo,_NFL);
865 rFlag=1;
866 _NFL->SAVE=W_SAVE_TIME=time(NULL); } }
867
868 if(rFlag==0) Write_List_2_File(polyo,_NFL);
869 _NFL->TIME=time(NULL);
870 fputs(ctime(&_NFL->TIME),stdout);
871 free(_NFL);
872 }
873
874
875 /* ====================================================================== */
876 /* ========== ========== */
877 /* ========== S U B L A T T I C E S ========== */
878 /* ========== ========== */
879 /* ====================================================================== */
880
Write_VPM(FILE * F,int * v,int * f,subl_int x[][VERT_Nmax])881 void Write_VPM(FILE *F,int *v, int *f, subl_int x[][VERT_Nmax]){
882 int i,j;
883 fprintf(F,"%d %d\n",*v,*f);
884 for(i=0;i<*f;i++) {
885 for(j=0;j<*v-1;j++) fprintf(F,"%lld ",(long long) x[i][j]);
886 fprintf(F,"%lld\n",(long long) x[i][j]); }
887 }
888
Write_M(FILE * F,int * v,int * f,subl_int x[][POLY_Dmax])889 void Write_M(FILE *F,int *v, int *f, subl_int x[][POLY_Dmax])
890 { int i,j; fprintf(F,"%d %d\n",*v,*f);
891 for(i=0;i<*f;i++)
892 { for(j=0;j<*v-1;j++) fprintf(F,"%lld ",(long long) x[i][j]);
893 fprintf(F,"%lld\n",(long long) x[i][j]);
894 }
895 }
896
897
Make_All_Sublat(NF_List * _L,int n,int v,subl_int diag[POLY_Dmax],subl_int u[][VERT_Nmax],char * mFlag,PolyPointList * _P)898 void Make_All_Sublat(NF_List *_L, int n, int v, subl_int diag[POLY_Dmax],
899 subl_int u[][VERT_Nmax], char *mFlag, PolyPointList *_P){
900 /* create all inequivalent decompositions diag=s*t into upper
901 triangular matrices s,t;
902 t*(first lines of u) becomes the poly P;
903 choose the elements of t (columns rising, line # falling),
904 calculate corresponding element of s at the same time;
905 finally calculate P and add to list */
906
907
908 VertexNumList V;
909 EqList F;
910 subl_int s[POLY_Dmax][POLY_Dmax], t[POLY_Dmax][POLY_Dmax], ts_lincol;
911 int i,j,k, lin=0, col=0, err=0;
912
913 _L->d=_P->n=n;
914 _P->np=V.nv=v;
915 for(i=0;i<n;i++) for (j=0;j<n;j++) t[i][j]=(i>=j)-1;
916 while (col>=0){
917 /* fprintf(outFILE,"t[%d][%d]=%lld\n",lin, col, t[lin][col]); */
918 if (lin==col){
919 do
920 t[lin][col]++;
921 while ((diag[lin]%t[lin][col])&&(t[lin][col]<=diag[lin]));
922 if (t[lin][col]>diag[lin]){
923 if (col==0) return;
924 t[lin][col]=0;
925 lin=0;
926 col--;}
927 else{
928 s[lin][col]=diag[lin]/t[lin][col];
929 if (col==0) lin=col=1;
930 else lin--; } }
931 else{
932 /* (t*s)[lin][col]=sum_j t[lin][j]*s[j][col]=0 */
933 ts_lincol=0;
934 for (j=col;j>lin;j--) ts_lincol+=t[lin][j]*s[j][col];
935 do{
936 t[lin][col]++;
937 ts_lincol+=s[col][col];}
938 while ((ts_lincol%t[lin][lin])&&(t[lin][col]<t[lin][lin]));
939 if (t[lin][col]>=t[lin][lin]){
940 t[lin][col]=-1;
941 lin++;}
942 else{
943 s[lin][col]=-ts_lincol/t[lin][lin];
944 if (lin) lin--;
945 else if (col<n-1) {col++; lin=col;}
946 else{
947 for (i=0;i<n;i++) for (j=0;j<v;j++){
948 /* _P->x[j][i]=s[i][i]*u[i][j];
949 for (k=i+1;k<n;k++) _P->x[j][i]+=s[i][k]*u[k][j]; } */
950 subl_int Pji=s[i][i]*u[i][j];
951 for (k=i+1;k<n;k++) Pji+=s[i][k]*u[k][j];
952 _P->x[j][i]=Pji;
953 if (labs(Pji)>20000000) err=1;}
954 _P->np=v;
955 if (err){
956 Print_PPL(_P,"");
957 fprintf(outFILE,"u: "); Write_VPM(outFILE,&v,&v,u);
958 fprintf(outFILE,"s: "); Write_M(outFILE,&n,&n,s);
959 fprintf(outFILE,"t: "); Write_M(outFILE,&n,&n,t);}
960 for(i=0;i<v;i++) V.v[i]=i;
961 New_Improve_Coords(_P,&V);
962
963 if (!IP_Check(_P,&V,&F)){
964 puts("IP_Check failed in Make_All_Sublat!\n"); exit(0);}
965 for (i=0;i<F.ne;i++) if(F.e[i].c!=1) {
966 fprintf(outFILE,"Not reflexive in Make_All_Sublat!\n");
967 exit(0);}
968 if (*mFlag != 'r') Add_NF_to_List(_P,&V,&F,_L);
969 else if (Poly_Max_check(_P,&V,&F)) if(Add_NF_to_List(_P,&V,&F,_L)){
970 Print_PPL(_P,""); }}}}}
971 }
972
Line_Recomb(int i,int line,int v,int f,subl_int x[][VERT_Nmax])973 int Line_Recomb(int i, int line, int v, int f, subl_int x[][VERT_Nmax]){
974 int j, k, changes=0;
975 subl_int y11, y12, y21, y22, g, new_;
976 if (abs(x[i][i])==1) return 0;
977 for (j=line;j<f;j++){
978 /* take a combination of i'th and j'th line such that x[i][i] becomes
979 gcd(x[i][i],x[j][i]) */
980 if (!x[j][i]) continue;
981 if (!x[i][i]||(x[j][i]%x[i][i])){
982 changes=1;
983 if ((g=LEgcd(x[i][i],x[j][i],&y11,&y12))<0) {y11*=-1; y12*=-1; g*=-1;}
984 y21=-x[j][i]/g;
985 y22=x[i][i]/g;
986 for (k=i;k<v;k++){
987 new_=y11*x[i][k]+y12*x[j][k];
988 x[j][k]=y21*x[i][k]+y22*x[j][k];
989 x[i][k]=new_; } } }
990 return changes;
991 }
992
Col_Recomb(int i,int col,int v,int f,subl_int x[][VERT_Nmax],subl_int u[][VERT_Nmax])993 void Col_Recomb(int i, int col, int v, int f, subl_int x[][VERT_Nmax],
994 subl_int u[][VERT_Nmax] ){
995 /* take a combination of i'th and col'th column such that x[i][i] becomes
996 gcd(x[i][i],x[i][col]); change u accordingly */
997 int k;
998 subl_int y11, y12, y21, y22, g, new_; /* "new" conflicts with g++ => new_ */
999 if ((g=LEgcd(x[i][i],x[i][col],&y11,&y12))<0) {y11*=-1; y12*=-1; g*=-1;}
1000 y21=-x[i][col]/g;
1001 y22=x[i][i]/g;
1002 for (k=i;k<f;k++){
1003 new_=y11*x[k][i]+y12*x[k][col];
1004 x[k][col]=y21*x[k][i]+y22*x[k][col];
1005 x[k][i]=new_; }
1006 for (k=0;k<v;k++){
1007 new_=y22*u[i][k]-y21*u[col][k];
1008 u[col][k]=-y12*u[i][k]+y11*u[col][k];
1009 u[i][k]=new_; }
1010 }
1011
lin_col_swap(int i,int lin,int col,int v,int f,subl_int x[][VERT_Nmax],subl_int u[][VERT_Nmax])1012 void lin_col_swap(int i, int lin, int col, int v, int f,
1013 subl_int x[][VERT_Nmax], subl_int u[][VERT_Nmax] ){
1014 /* swap i'th and lin'th line, i'th and col'th column, change u accordingly */
1015 int k;
1016 subl_int swap_vec[VERT_Nmax];
1017 if (lin!=i) for (k=i;k<v;k++){
1018 swap_vec[k]=x[lin][k];
1019 x[lin][k]=x[i][k];
1020 x[i][k]=swap_vec[k];}
1021 if (col!=i){
1022 for (k=i;k<f;k++){
1023 swap_vec[k]=x[k][col];
1024 x[k][col]=x[k][i];
1025 x[k][i]=swap_vec[k];}
1026 for (k=0;k<v;k++){
1027 swap_vec[k]=u[col][k];
1028 u[col][k]=u[i][k];
1029 u[i][k]=swap_vec[k];}}
1030 }
1031
1032
MakePolyOnSublat(NF_List * _L,subl_int x[VERT_Nmax][VERT_Nmax],int v,int f,int * max_order,char * mFlag,PolyPointList * _P)1033 void MakePolyOnSublat(NF_List *_L, subl_int x[VERT_Nmax][VERT_Nmax],
1034 int v, int f, int *max_order, char *mFlag,
1035 PolyPointList *_P){
1036
1037 /* Decompose the VPM x as x=w*diag*u, where w and u are SL(Z);
1038 the first lines of u are the coordinates on the coarsest lattices */
1039
1040 int i, j, k, col, lin;
1041 subl_int diag[POLY_Dmax];
1042 subl_int fac, order=1;
1043 subl_int u[VERT_Nmax][VERT_Nmax];
1044
1045 for (i=0;i<v;i++) for (j=0;j<v;j++) u[i][j]=(i==j);
1046
1047 for (i=0;i<VERT_Nmax;i++) {
1048 /* put gcd of remaining matrix elements at x[i][i], make zeroes at
1049 x[i][j] and x[j][i] for all j unequal to i */
1050 int chosen_lin=0, chosen_col=-1;
1051 subl_int min_entry, min_col_max=0, min_lin_gcd=0, min_lin_max=0, rem_gcd=0;
1052
1053 for (lin=i;lin<f;lin++){
1054 subl_int lin_gcd=0, lin_max=0;
1055 for (col=i;col<v;col++) if(x[lin][col]){
1056 if (!lin_gcd) lin_gcd=abs(x[lin][col]);
1057 else lin_gcd=NNgcd(lin_gcd,x[lin][col]);
1058 if (abs(x[lin][col])>lin_max) lin_max=abs(x[lin][col]);}
1059 if (lin_max) if ((!min_lin_gcd)||(lin_gcd<min_lin_gcd)||
1060 ((lin_gcd==min_lin_gcd)&&(lin_max<min_lin_max))){
1061 chosen_lin=lin;
1062 min_lin_gcd=lin_gcd;
1063 min_lin_max=lin_max;}
1064 rem_gcd=NNgcd(rem_gcd, lin_gcd);}
1065
1066 if (!min_lin_max) break;
1067
1068 min_entry=min_lin_max+1;
1069 for (col=i;col<v;col++)
1070 if(x[chosen_lin][col]) if(abs(x[chosen_lin][col])<=min_entry){
1071 subl_int col_max=0;
1072 for (lin=i;lin<f;lin++) if (abs(x[lin][col])>col_max)
1073 col_max=abs(x[lin][col]);
1074 if (!col_max){printf("col_max==0!!!"); exit(0);}
1075 if ((chosen_col==-1)||(abs(x[chosen_lin][col])<min_entry)||
1076 ((abs(x[chosen_lin][col])==min_entry)&&(col_max<min_col_max))){
1077 chosen_col=col;
1078 min_col_max=col_max;}}
1079
1080 lin_col_swap(i, chosen_lin, chosen_col, v, f, x, u);
1081
1082 /* put gcd of remaining matrix elements at x[i][i] */
1083 while(abs(x[i][i])>rem_gcd){
1084 for (col=i+1;col<v;col++) if ((x[i][col]%x[i][i]))
1085 Col_Recomb(i,col,v,f,x,u);
1086 if (Line_Recomb(i,i+1,v,f,x)) continue;
1087 if (abs(x[i][i])<=rem_gcd) break;
1088 for (col=i+1;col<v;col++){
1089 for(j=i+1;(j<f)&&(!(x[j][col]%x[i][i]));j++);
1090 if (j<f){
1091 /* set i'th element of j'th line to zero, then add j'th line to
1092 i'th line, change w accordingly */
1093 fac=x[j][i]/x[i][i];
1094 for (k=i;k<v;k++){
1095 x[j][k]-=fac*x[i][k];
1096 x[i][k]+=x[j][k];}
1097 break;}}}
1098
1099 /* make zeroes */
1100 for (col=i+1;col<v;col++) {
1101 fac=x[i][col]/x[i][i];
1102 for (k=i;k<f;k++) x[k][col]-=fac*x[k][i];
1103 for (k=0;k<v;k++) u[i][k]+=fac*u[col][k]; }
1104 for (k=i+1;k<v;k++) if (x[i][k]){
1105 printf("error in MakePolyOnSublat!!!\n"); exit(0);}}
1106
1107 for (j=0;j<i;j++) order*=(diag[j]=abs(x[j][j]));
1108 if (order>*max_order) *max_order=order;
1109 if (i>POLY_Dmax) {printf("diag has %d entries!!!\n", i); exit(0);}
1110 if (order>1) Make_All_Sublat(_L, i, v, diag, u, mFlag, _P);
1111 }
1112
uc_nf_to_P(PolyPointList * _P,int * MS,int * d,int * v,int * nuc,unsigned char * uc)1113 void uc_nf_to_P(PolyPointList *_P, int *MS, int *d, int *v, int *nuc,
1114 unsigned char *uc){
1115 Long tNF[POLY_Dmax][VERT_Nmax]; int i, j;
1116 UCnf2vNF(d, v, nuc, uc, tNF, MS);
1117 _P->n=*d; _P->np=*v; (*MS)%=4;
1118 for(i=0;i<*v;i++) for(j=0;j<*d;j++) _P->x[i][j]=tNF[j][i];
1119 }
1120
Find_Sublat_Polys(char mFlag,char * dbin,char * polyi,char * polyo,PolyPointList * _P)1121 void Find_Sublat_Polys(char mFlag, char *dbin, char *polyi, char *polyo,
1122 PolyPointList *_P){
1123 NF_List *_NFL=(NF_List *) malloc(sizeof(NF_List));
1124 VertexNumList Vnl; EqList Fel;
1125 subl_int x[VERT_Nmax][VERT_Nmax], y[VERT_Nmax][VERT_Nmax];
1126 time_t Tstart;
1127 int v, nu, i, j, k, max_order=1;
1128
1129 assert(_NFL!=NULL);
1130 if(!(*polyo)) {
1131 puts("You have to specify an output file via -po in -sm-mode.");
1132 printf("For more help use option `-h'\n");
1133 exit(0);}
1134 _NFL->of=0; _NFL->rf=0;
1135 _NFL->iname=polyi; _NFL->oname=polyo; _NFL->dbname=dbin;
1136 Init_NF_List(_NFL);
1137 _NFL->SL=0;
1138
1139 if (*dbin){
1140 DataBase *DB=&(_NFL->DB);
1141 char *dbname = (char *) malloc(1+strlen(dbin)+File_Ext_NCmax), *fx;
1142 unsigned char uc_poly[NUC_Nmax];
1143 int MS;
1144
1145 strcpy(dbname,dbin);
1146 strcat(dbname,".");
1147 fx=&dbname[strlen(dbin)+1];
1148
1149 printf("Reading DB-files, calculating sublattices:\n"); fflush(0);
1150
1151 /* read the DB-files and calculate Sublattices */
1152 for (v=2;v<=DB->nVmax;v++) if(DB->nNUC[v]){
1153 FILE *dbfile;
1154 char ext[4]={'v',0,0,0};
1155 ext[1]='0' + v / 10; ext[2]='0' + v % 10;
1156 Tstart=time(NULL);
1157 strcpy(fx,ext);
1158 dbfile=fopen(dbname,"rb");
1159 assert(dbfile!=NULL);
1160 if (!mFlag) {printf("Reading %s\n", dbname); fflush(0);}
1161 for (nu=0;nu<=DB->NUCmax;nu++) for (i=0; i<DB->NFnum[v][nu]; i++){
1162 for (j=0; j<nu; j++) uc_poly[j]=fgetc(dbfile);
1163 uc_nf_to_P(_P, &MS, &(_NFL->d), &v, &nu, uc_poly);
1164 if (MS>3) {printf("MS=%d!!!\n",MS); exit(0);}
1165 assert(IP_Check(_P,&Vnl,&Fel));
1166 assert(v==Vnl.nv);
1167 /* compute VPM */
1168 for(j=0;j<Fel.ne;j++) for(k=0;k<Vnl.nv;k++)
1169 x[j][k]=Eval_Eq_on_V(&(Fel.e[j]),_P->x[Vnl.v[k]],_P->n)-1;
1170 for(j=0;j<Fel.ne;j++) for(k=0;k<Vnl.nv;k++) y[k][j]=x[j][k];
1171 if (MS!=2)
1172 MakePolyOnSublat(_NFL, x, Vnl.nv, Fel.ne, &max_order, &mFlag, _P);
1173 if (MS>1)
1174 MakePolyOnSublat(_NFL, y, Fel.ne, Vnl.nv, &max_order, &mFlag, _P); }
1175
1176 if(ferror(dbfile)) {printf("File error in %s\n",dbname); exit(0);}
1177 fclose(dbfile);
1178 printf(" %dp (%ds)\n", (int)_NFL->NP, (int) difftime(time(NULL),Tstart));
1179 fflush(0); }}
1180
1181 else {
1182 CWS W;
1183 while(Read_CWS_PP(&W,_P)){
1184 assert(IP_Check(_P,&Vnl,&Fel));
1185 /* compute VPM */
1186 for(j=0;j<Fel.ne;j++) for(i=0;i<Vnl.nv;i++)
1187 x[j][i]=Eval_Eq_on_V(&(Fel.e[j]),_P->x[Vnl.v[i]],_P->n)-1;
1188 MakePolyOnSublat(_NFL, x, Vnl.nv, Fel.ne, &max_order, &mFlag, _P);
1189 if (!mFlag) Print_Weight_Info(&W,_NFL);}}
1190
1191 printf("max_order=%d\n", max_order);
1192 Write_List_2_File(polyo,_NFL);
1193 _NFL->TIME=time(NULL); fputs(ctime(&_NFL->TIME),stdout);
1194 free(_NFL);
1195 }
1196
1197
1198
1199 /* ====================================================================== */
1200 /* ========== ========== */
1201 /* ========== I R , V I R & M A X ========== */
1202 /* ========== ========== */
1203 /* ====================================================================== */
1204
irred(PolyPointList * _P)1205 int irred(PolyPointList *_P)
1206 {
1207 int i,j,drop_point[POLY_Dmax];
1208 VertexNumList V,NV; EqList E;
1209
1210 for (j=0;j<_P->n;j++) drop_point[j]=0; /* just to silence the compiler */
1211 IP_Check(_P,&V,&E);
1212 for(i=0;i<V.nv;i++){
1213
1214 /* Dropping the vertex V.v[i]: */
1215 _P->np--;
1216 for (j=0;j<_P->n;j++){
1217 drop_point[j]=_P->x[V.v[i]][j];
1218 _P->x[V.v[i]][j]=_P->x[_P->np][j]; }
1219 if (IP_Check(_P,&NV,&E)) return 0;
1220
1221 /* Reconstruct _P */
1222 for (j=0;j<_P->n;j++){
1223 _P->x[_P->np][j]=_P->x[V.v[i]][j];
1224 _P->x[V.v[i]][j]=drop_point[j];}
1225 _P->np++; }
1226 return 1;
1227 }
1228
DPircheck(CWS * _W,PolyPointList * _P)1229 void DPircheck(CWS *_W, PolyPointList *_P){
1230 int i, j, k;
1231 EqList E,DE;
1232 VertexNumList V,DV;
1233 EqList *B=&_W->B;
1234 PolyPointList *_RDP= (PolyPointList *) malloc(sizeof (PolyPointList));
1235 PolyPointList *_PD = (PolyPointList *) malloc(sizeof (PolyPointList));
1236 if(_RDP==NULL) {puts("Unable to allocate space for _P"); exit(0);}
1237 if(_PD==NULL) {printf("Unable to allocate _PD\n"); exit(0);}
1238
1239 IP_Check(_P,&V,&E);
1240 Make_Dual_Poly(_P,&V,&E,_PD);
1241 _RDP->n=_P->n;
1242 _RDP->np=B->ne;
1243 for (i=0;i<B->ne;i++) for (j=0;j<_P->n;j++) _RDP->x[i][j]=B->e[i].a[j];
1244 IP_Check(_RDP,&DV,&DE);
1245 _RDP->np=0;
1246 for (i=0;i<_PD->np;i++){
1247 int Good=1;
1248 for (k=0;k<DE.ne;k++){
1249 int Pairing=DE.e[k].c;
1250 for (j=0;j<_P->n;j++) Pairing+=DE.e[k].a[j]*_PD->x[i][j];
1251 if (Pairing<0) Good=0;}
1252 if(Good){
1253 for (j=0;j<_P->n;j++) _RDP->x[_RDP->np][j]=_PD->x[i][j];
1254 _RDP->np++;}}
1255 /* Print_PPL(_RDP); */
1256 if (irred(_RDP)){
1257 /* Print_PPL(_RDP); */
1258 for(i=0;i<_W->nw;i++){
1259 fprintf(outFILE,"%d ",(int) _W->d[i]);
1260 for(j=0;j<_W->N;j++) fprintf(outFILE,"%d ",(int) _W->W[i][j]);
1261 if(i+1<_W->nw) fprintf(outFILE," "); }
1262 fprintf(outFILE,"\n"); }
1263 free(_PD); free(_RDP);
1264 }
1265
1266
virred(PolyPointList * _P,EqList * B)1267 int virred(PolyPointList *_P, EqList *B)
1268 {
1269 int i,j,k,drop_point[POLY_Dmax],equal;
1270 VertexNumList V; EqList E;
1271 assert(Ref_Check(_P,&V,&E));
1272 for (j=0;j<_P->n;j++) drop_point[j]=0; /* just to silence the compiler */
1273 /* for(i=0;(i<V.nv);i++) {
1274 for (j=0;(j<_P->n);j++) fprintf(outFILE,"%d ",(int) _P->x[V.v[i]][j]);
1275 fprintf(outFILE,"\n");}
1276 fprintf(outFILE,"_CL->B:\n");
1277 for(k=0;k<B->ne;k++){
1278 for (j=0;(j<_P->n);j++) fprintf(outFILE,"%d ",(int) _CL->B[k][j]);
1279 fprintf(outFILE,"\n");} */
1280 for(k=0;k<B->ne;k++){
1281 equal=0;
1282 for(i=0;(i<V.nv)&&(!equal);i++){
1283 equal=1;
1284 for (j=0;(j<_P->n)&&equal;j++)
1285 if (_P->x[V.v[i]][j]-B->e[k].a[j]) equal=0; }
1286 /* if (equal) break;}*/
1287 i--;
1288 if (!equal){
1289 fprintf(outFILE,"Vertex not found!\n");
1290 for(i=0;(i<V.nv);i++) {
1291 for (j=0;(j<_P->n);j++) fprintf(outFILE,"%d ",(int) _P->x[V.v[i]][j]);
1292 fprintf(outFILE,"\n");}
1293 for (j=0;j<_P->n;j++) fprintf(outFILE,"%d ",(int) _P->x[V.v[i]][j]);
1294 fprintf(outFILE,"\n");
1295 for (j=0;j<_P->n;j++) fprintf(outFILE,"%d ",(int) B->e[k].a[j]);
1296 fprintf(outFILE,"\n");
1297 return 1;}
1298 /* Dropping the vertex V.v[i]: */
1299 _P->np--;
1300 for (j=0;j<_P->n;j++){
1301 drop_point[j]=_P->x[V.v[i]][j];
1302 _P->x[V.v[i]][j]=_P->x[_P->np][j]; }
1303 if (IP_Check(_P,&V,&E)) return 0;
1304
1305 /* Reconstruct _P */
1306 for (j=0;j<_P->n;j++){
1307 _P->x[_P->np][j]=_P->x[V.v[i]][j];
1308 _P->x[V.v[i]][j]=drop_point[j];}
1309 _P->np++;
1310 assert(Ref_Check(_P,&V,&E));}
1311 return 1;
1312 }
1313
DPvircheck(CWS * _W,PolyPointList * _P)1314 void DPvircheck(CWS *_W, PolyPointList *_P) {
1315 int i, j;
1316 EqList E; EqList *B=&_W->B;
1317 VertexNumList V;
1318 PolyPointList *_PD = (PolyPointList *) malloc(sizeof (PolyPointList));
1319 if(_PD==NULL) {printf("Unable to allocate _PD\n"); exit(0);}
1320
1321 Ref_Check(_P,&V,&E);
1322 Make_Dual_Poly(_P,&V,&E,_PD);
1323 if (virred(_PD,B)){
1324 /* Print_PPL(_PD); */
1325 for(i=0;i<_W->nw;i++){
1326 fprintf(outFILE,"%d ",(int) _W->d[i]);
1327 for(j=0;j<_W->N;j++) fprintf(outFILE,"%d ",(int) _W->W[i][j]);
1328 if(i+1<_W->nw) fprintf(outFILE," "); }
1329 fprintf(outFILE,"\n"); }
1330 free(_PD);
1331 }
1332
1333 int Find_Ref_Subpoly(PolyPointList *_P, EqList *_E, VertexNumList *_V,
1334 KeepList *_KL, int *rd);
1335
Find_RSP_Drop_and_Keep(PolyPointList * _P,VertexNumList * _V,EqList * _E,KeepList * _KL,int drop_num,int * rd)1336 int Find_RSP_Drop_and_Keep(PolyPointList *_P, VertexNumList *_V, EqList *_E,
1337 KeepList *_KL, int drop_num, int *rd){
1338
1339 int i, j, n_irrel=0, IP/*, test_IP*/;
1340 int *new2old =(int *) malloc(sizeof(int)*POINT_Nmax);
1341 Long drop_point[POLY_Dmax];
1342 VertexNumList red_V/*, test_V*/;
1343 EqList new_E/*, test_E*/;
1344 PolyPointList *red_P = (PolyPointList *) malloc(sizeof (PolyPointList));
1345
1346 if(red_P==NULL) {puts("Unable to allocate space for red_P"); exit(0);}
1347 if(new2old==NULL) {puts("Unable to allocate space for new2old"); exit(0);}
1348
1349 for (j=0;j<_P->n;j++) drop_point[j]=_P->x[_V->v[drop_num]][j];
1350
1351 /* Create new_E: Same as *_E, but irrelevant facets first */
1352
1353 new_E.ne=_E->ne;
1354 for (i=0;i<_E->ne;i++){
1355 if (Eval_Eq_on_V(&(_E->e[i]),drop_point,_P->n)){
1356 for (j=0;j<_P->n;j++) new_E.e[n_irrel].a[j]=_E->e[i].a[j];
1357 new_E.e[n_irrel++].c=_E->e[i].c;}
1358 else {
1359 for (j=0;j<_P->n;j++) new_E.e[_E->ne+n_irrel-i-1].a[j]=_E->e[i].a[j];
1360 new_E.e[_E->ne+n_irrel-i-1].c=_E->e[i].c;}}
1361
1362 /* Create red_P: old vertices, points that are not on irrelevant facets */
1363 red_P->np=0;
1364 red_P->n=_P->n;
1365 for (i=0;i<_V->nv;i++) if (i!=drop_num){
1366 for (j=0;j<_P->n;j++) red_P->x[red_P->np][j]=_P->x[_V->v[i]][j];
1367 new2old[red_P->np++]=_V->v[i];}
1368 for (i=0;i<_P->np;i++)
1369 if ((i!=_V->v[drop_num])&&Relevant(_P->x[i],&_P->n,&new_E,&n_irrel)){
1370 for (j=0;j<_P->n;j++) red_P->x[red_P->np][j]=_P->x[i][j];
1371 new2old[red_P->np++]=i;}
1372
1373 IP=Aided_IP_Check(red_P,&red_V,&new_E,n_irrel,_V->nv-1);
1374
1375 /* test_IP=IP_Check(red_P,&test_V,&test_E);
1376 printf("%d %d %d\n", _P->np, red_P->np, *rd); fflush(0);
1377 printf("%d %d %d\n", test_V.nv, test_E.ne, test_IP);
1378 printf("%d %d \n", red_V.nv, new_E.ne); fflush(0);
1379 if (IP) if ((!test_IP)||(test_V.nv!=red_V.nv)||
1380 (test_E.ne!=new_E.ne)){
1381 int k;
1382 fprintf(outFILE,"red_P:\n");fflush(0);
1383 Print_PPL(red_P);
1384 fprintf(outFILE,"\n");
1385 fprintf(outFILE,"red_V: ");
1386 for (i=0;i<red_V.nv;i++) fprintf(outFILE,"%d ",red_V.v[i]);
1387 fprintf(outFILE,"\n");
1388 fprintf(outFILE,"new_E:\n");
1389 for (i=0;i<new_E.ne;i++){
1390 for (k=0;k<_P->n;k++) fprintf(outFILE,"%d ", (int) new_E.e[i].a[k]);
1391 fprintf(outFILE," %d\n", (int) new_E.e[i].c);}
1392 fprintf(outFILE,"\n");
1393 fprintf(outFILE,"n_irrel: %d\n",n_irrel);
1394 fprintf(outFILE,"new2old:");
1395 for (i=0;i<red_P->np;i++) fprintf(outFILE," %d", new2old[i]);
1396 fprintf(outFILE,"\n");
1397 fprintf(outFILE,"test_V: ");
1398 for (i=0;i<test_V.nv;i++) fprintf(outFILE,"%d ",test_V.v[i]);
1399 fprintf(outFILE,"\n");
1400 fprintf(outFILE,"test_E:\n");
1401 for (i=0;i<test_E.ne;i++){
1402 for (k=0;k<_P->n;k++) fprintf(outFILE,"%d ", (int) test_E.e[i].a[k]);
1403 fprintf(outFILE," %d\n", (int) test_E.e[i].c);}
1404 fprintf(outFILE,"\n");
1405 exit(0);} */
1406
1407
1408 free(red_P);
1409
1410 if(IP){
1411 VertexNumList new_V;
1412
1413 /* Create new_V from red_V */
1414
1415 new_V.nv=red_V.nv;
1416 for (i=0;i<red_V.nv;i++) if((new_V.v[i]=new2old[red_V.v[i]])==_P->np-1)
1417 new_V.v[i]=_V->v[drop_num];
1418
1419 /* Drop the vertex _V->v[drop_num]: */
1420 _P->np--;
1421 for (j=0;j<_P->n;j++) _P->x[_V->v[drop_num]][j]=_P->x[_P->np][j];
1422 j=kept(_P->np,_KL);
1423 if (j>=0) _KL->k[j]=_V->v[drop_num];
1424
1425 if (Find_Ref_Subpoly(_P, &new_E, &new_V, _KL, rd)) {
1426 free(new2old);
1427 return 1;}
1428
1429 /* Reconstruct _P */
1430 if (j>=0) _KL->k[j]=_P->np;
1431 for (j=0;j<_P->n;j++){
1432 _P->x[_P->np][j]=_P->x[_V->v[drop_num]][j];
1433 _P->x[_V->v[drop_num]][j]=drop_point[j];}
1434 _P->np++;}
1435
1436 /* Attach a keep label to the vertex _V->v[drop_num]: */
1437 _KL->k[_KL->nk]=_V->v[drop_num];
1438 _KL->nk++;
1439 free(new2old);
1440 return 0;
1441 }
1442
1443
Start_Find_Ref_Subpoly(PolyPointList * _P)1444 int Start_Find_Ref_Subpoly(PolyPointList *_P){
1445 int i,j;
1446 VertexNumList V, new_V;
1447 EqList E, new_E;
1448 KeepList KL;
1449 Long drop_point[POLY_Dmax];
1450 int rd=0;
1451
1452 if(!IP_Check(_P,&V,&E)) {
1453 fprintf(outFILE,"IP_check negative!\n"); exit(0);}
1454
1455 KL.nk=0;
1456 for (j=0;j<_P->n;j++) drop_point[j]=0; /* just to silence the compiler */
1457
1458 for(i=0;i<V.nv;i++) {
1459
1460 /* Try dropping the vertex V.v[i]: */
1461 _P->np--;
1462 for (j=0;j<_P->n;j++){
1463 drop_point[j]=_P->x[V.v[i]][j];
1464 _P->x[V.v[i]][j]=_P->x[_P->np][j]; }
1465
1466 if (!IP_Check(_P,&new_V,&new_E)) KL.k[KL.nk++]=V.v[i];
1467
1468 /* Reconstruct _P */
1469
1470 for (j=0;j<_P->n;j++){
1471 _P->x[_P->np][j]=_P->x[V.v[i]][j];
1472 _P->x[V.v[i]][j]=drop_point[j];}
1473 _P->np++; }
1474
1475 return Find_Ref_Subpoly(_P, &E, &V, &KL, &rd);
1476 }
1477
1478
Find_Ref_Subpoly(PolyPointList * _P,EqList * _E,VertexNumList * _V,KeepList * _KL,int * rd)1479 int Find_Ref_Subpoly(PolyPointList *_P, EqList *_E, VertexNumList *_V,
1480 KeepList *_KL, int *rd){
1481
1482 /* Creates all subpolyhedra of an IP-polyhedron given by _P, _E, _V.
1483 The basic structure is:
1484 Search for bad facet;
1485 If (reflexive) if (!Add_NF_to_List) return;
1486 Drop_and_Keep vertices (on bad facet / everywhere);
1487 If (bad facet exists) Reduce_Poly; */
1488
1489 int i, j, badfacet=-1, nbadvert=0, maxnkept=-1;
1490 KeepList new_KL=*_KL;
1491
1492 /* Choose bad facet with largest number of kept vertices: */
1493
1494 for (i=0;i<_E->ne;i++) if (_E->e[i].c>1){
1495 int nkept=0;
1496 for (j=0;j<_KL->nk;j++)
1497 if (!Eval_Eq_on_V(&(_E->e[i]),_P->x[_KL->k[j]],_P->n)) nkept++;
1498 if (nkept>maxnkept){
1499 maxnkept=nkept;
1500 badfacet=i; }}
1501
1502 if (*rd>0) if (badfacet==-1) return 1;
1503
1504 if (badfacet!=-1) {
1505 /* Find vertices on bad facet: */
1506 nbadvert=0;
1507 for (i=0;i<_V->nv;i++)
1508 if (!Eval_Eq_on_V(&(_E->e[badfacet]),_P->x[_V->v[i]],_P->n)){
1509 j=_V->v[nbadvert];
1510 _V->v[nbadvert]=_V->v[i];
1511 _V->v[i]=j;
1512 nbadvert++; } }
1513
1514 (*rd)++;
1515
1516 for(i=0;i<((badfacet==-1)?_V->nv:nbadvert);i++)
1517 if((kept(_V->v[i],&new_KL)<0))
1518 if (Find_RSP_Drop_and_Keep(_P, _V, _E, &new_KL, i, rd)) return 1;
1519
1520 (*rd)--;
1521 return 0;
1522 }
1523
Max_check(CWS * _W,PolyPointList * _P)1524 void Max_check(CWS *_W, PolyPointList *_P) {
1525 int i, j;
1526 EqList E;
1527 VertexNumList V;
1528 IP_Check(_P,&V,&E);
1529 if (Poly_Max_check(_P,&V,&E)){
1530 for(i=0;i<_W->nw;i++){
1531 fprintf(outFILE,"%d ",(int) _W->d[i]);
1532 for(j=0;j<_W->N;j++) fprintf(outFILE,"%d ",(int) _W->W[i][j]);
1533 if(i+1<_W->nw) fprintf(outFILE," "); }
1534 fprintf(outFILE,"\n"); }
1535 }
1536
Poly_Max_check(PolyPointList * _P,VertexNumList * _V,EqList * _E)1537 int Poly_Max_check(PolyPointList *_P, VertexNumList *_V, EqList *_E){
1538 int rm;
1539 PolyPointList *_PD = (PolyPointList *) malloc(sizeof(PolyPointList));
1540 assert(_PD!=NULL);
1541 Make_Dual_Poly(_P,_V,_E,_PD);
1542 rm=!Start_Find_Ref_Subpoly(_PD);
1543 free(_PD);
1544 return rm;
1545 }
1546
Poly_Min_check(PolyPointList * _P,VertexNumList * _V,EqList * _E)1547 int Poly_Min_check(PolyPointList *_P, VertexNumList *_V, EqList *_E){
1548 PairMat PM;
1549 Make_VEPM(_P,_V,_E,PM);
1550 Complete_Poly(PM, _E, _V->nv, _P);
1551 return !Start_Find_Ref_Subpoly(_P);
1552 }
1553
Overall_check(CWS * _W,PolyPointList * _P)1554 void Overall_check(CWS *_W, PolyPointList *_P) {
1555 int i, j, k, span, lpm=0, vm=0, r=0;
1556 EqList E, DE;
1557 VertexNumList V, DV;
1558 EqList *B=&_W->B;
1559 PolyPointList *_RDP= (PolyPointList *) malloc(sizeof (PolyPointList));
1560 PolyPointList *_PD = (PolyPointList *) malloc(sizeof (PolyPointList));
1561 PolyPointList *_PD2 = (PolyPointList *) malloc(sizeof (PolyPointList));
1562 if(_RDP==NULL) {puts("Unable to allocate space for _P"); exit(0);}
1563 if(_PD==NULL) {printf("Unable to allocate _PD\n"); exit(0);}
1564 if(_PD2==NULL) {printf("Unable to allocate _PD\n"); exit(0);}
1565
1566 if ((_P->n<5) ? !IP_Check(_P,&V,&E) : !Ref_Check(_P,&V,&E)) {
1567 free(_PD); free(_PD2); free(_RDP); return;}
1568
1569 /* assert(Ref_Equations(&E)); */
1570
1571 for(i=0;i<_W->nw;i++){
1572 fprintf(outFILE,"%d ",(int) _W->d[i]);
1573 for(j=0;j<_W->N;j++) fprintf(outFILE,"%d ",(int) _W->W[i][j]);
1574 if(i+1<_W->nw) fprintf(outFILE," "); }
1575 fflush(0);
1576
1577 span=Span_Check(&E, B, &_P->n);
1578
1579 Make_Dual_Poly(_P,&V,&E,_PD);
1580 *_PD2=*_PD;
1581
1582 _RDP->n=_P->n;
1583 _RDP->np=B->ne;
1584 for (i=0;i<B->ne;i++) for (j=0;j<_P->n;j++) _RDP->x[i][j]=B->e[i].a[j];
1585 IP_Check(_RDP,&DV,&DE);
1586 _RDP->np=0;
1587 for (i=0;i<_PD->np;i++){
1588 int Good=1;
1589 for (k=0;k<DE.ne;k++){
1590 int Pairing=DE.e[k].c;
1591 for (j=0;j<_P->n;j++) Pairing+=DE.e[k].a[j]*_PD->x[i][j];
1592 if (Pairing<0) Good=0;}
1593 if(Good){
1594 for (j=0;j<_P->n;j++) _RDP->x[_RDP->np][j]=_PD->x[i][j];
1595 _RDP->np++;}}
1596 if (irred(_RDP)) lpm=1;
1597
1598 if (!Start_Find_Ref_Subpoly(_PD)) r=1;
1599
1600 if (span) if (virred(_PD2,B)) vm=1;
1601
1602 if ((!span&&vm)||(!lpm&&vm)||(r!=vm))
1603 fprintf(outFILE,"span:%d lpm:%d vm:%d r:%d\n", span, lpm, vm, r);
1604 else{
1605 if(r) fprintf(outFILE,"r");
1606 if(lpm) fprintf(outFILE,"l");
1607 if(span) fprintf(outFILE,"s");
1608 if(!r&&!lpm&&!span) fprintf(outFILE,"-");
1609 fprintf(outFILE,"\n");}
1610 fflush(0);
1611 free(_PD); free(_PD2); free(_RDP);
1612 }
1613