?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
fixed a small bug in eval ehrhart
[polylib.git] / source / oldpolytest.c
blobb0c7b61ece208f7fe8bee6f468564677c158376b
1 /* polytest.c */
2 #include <stdio.h>
3 #include <polylib/polylib.h>
4
5 /* alpha.c
6 COPYRIGHT
7 Both this software and its documentation are
8
9 Copyright 1993 by IRISA /Universite de Rennes I - France,
10 Copyright 1995,1996 by BYU, Provo, Utah
11 all rights reserved.
12
13 Permission is granted to copy, use, and distribute
14 for any commercial or noncommercial purpose under the terms
15 of the GNU General Public license, version 2, June 1991
16 (see file : LICENSING).
17 */
18
19 #include <stdio.h>
20 #include <stdlib.h>
21 #include <string.h>
22
23 /*---------------------------------------------------------------------*/
24 /* int exist_points(pos,P,context) */
25 /* pos : index position of current loop index (0..hdim-1) */
26 /* P: loop domain */
27 /* context : context values for fixed indices */
28 /* recursive procedure, recurs for each imbriquation */
29 /* returns 1 if there exists any integer points in this polyhedron */
30 /* returns 0 if there are none */
31 /*---------------------------------------------------------------------*/
32 static int exist_points(int pos,Polyhedron *Pol,Value *context) {
33
34 Value LB, UB, k,tmp;
35
36 value_init(LB); value_init(UB);
37 value_init(k); value_init(tmp);
38 value_set_si(LB,0);
39 value_set_si(UB,0);
40
41 /* Problem if UB or LB is INFINITY */
42 if (lower_upper_bounds(pos,Pol,context,&LB,&UB) !=0) {
43 errormsg1("exist_points", "infdom", "infinite domain");
44 value_clear(LB);
45 value_clear(UB);
46 value_clear(k);
47 value_clear(tmp);
48 return -1;
49 }
50 value_set_si(context[pos],0);
51 if(value_lt(UB,LB)) {
52 value_clear(LB);
53 value_clear(UB);
54 value_clear(k);
55 value_clear(tmp);
56 return 0;
57 }
58 if (!Pol->next) {
59 value_substract(tmp,UB,LB);
60 value_increment(tmp,tmp);
61 value_clear(UB);
62 value_clear(LB);
63 value_clear(k);
64 return (value_pos_p(tmp));
65 }
66
67 for (value_assign(k,LB);value_le(k,UB);value_increment(k,k)) {
68
69 /* insert k in context */
70 value_assign(context[pos],k);
71 if (exist_points(pos+1,Pol->next,context) > 0 ) {
72 value_clear(LB); value_clear(UB);
73 value_clear(k); value_clear(tmp);
74 return 1;
75 }
76 }
77 /* Reset context */
78 value_set_si(context[pos],0);
79 value_clear(UB); value_clear(LB);
80 value_clear(k); value_clear(tmp);
81 return 0;
82 }
83
84 /*--------------------------------------------------------------*/
85 /* Test to see if there are any integral points in a polyhedron */
86 /*--------------------------------------------------------------*/
87 int Polyhedron_Not_Empty(Polyhedron *P,Polyhedron *C,int MAXRAYS) {
88
89 int res,i;
90 Value *context;
91 Polyhedron *L;
92
93 /* Create a context vector size dim+2 and set it to all zeros */
94 context = (Value *) malloc((P->Dimension+2)*sizeof(Value));
95
96 /* Initialize array 'context' */
97 for (i=0;i<(P->Dimension+2);i++)
98 value_init(context[i]);
99
100 Vector_Set(context,0,(P->Dimension+2));
101
102 /* Set context[P->Dimension+1] = 1 (the constant) */
103 value_set_si(context[P->Dimension+1],1);
104
105 L = Polyhedron_Scan(P,C,MAXRAYS);
106 res = exist_points(1,L,context);
107 Domain_Free(L);
108
109 /* Clear array 'context' */
110 for (i=0;i<(P->Dimension+2);i++)
111 value_clear(context[i]);
112 free(context);
113 return res;
114 }
115
116 /* PolyhedronLTQ ( P1, P2 ) */
117 /* P1 and P2 must be simple polyhedra */
118 /* result = 0 --> not comparable */
119 /* result = -1 --> P1 < P2 */
120 /* result = 1 --> P1 > P2 */
121 int PolyhedronLTQ (Polyhedron *Pol1,Polyhedron *Pol2,int NbMaxConstrs) {
122
123 int res, dim, i, j, k;
124 Polyhedron *Q1, *Q2, *Q3, *Q;
125 Matrix *Mat;
126
127 #define INDEX 1
128
129 if (Pol1->next || Pol2->next) {
130 errormsg1("PolyhedronLTQ", "compoly", "Can only compare polyhedra");
131 return 0;
132 }
133 if (Pol1->Dimension != Pol2->Dimension) {
134 errormsg1("PolyhedronLTQ","diffdim","Polyhedra are not same dimension");
135 return 0;
136 }
137 dim = Pol1->Dimension+2;
138
139 #ifdef DEBUG
140 fprintf(stdout, "P1\n");
141 Polyhedron_Print(stdout,P_VALUE_FMT,Pol1);
142 fprintf(stdout, "P2\n");
143 Polyhedron_Print(stdout,P_VALUE_FMT,Pol2);
144 #endif
145
146 /* Create the Line to add */
147 Mat = Matrix_Alloc(1,dim);
148 Vector_Set(Mat->p_Init,0,dim);
149 value_set_si(Mat->p[0][INDEX],1); /* line in INDEX dimension */
150
151 Q1 = AddRays(Mat->p[0],1,Pol1,NbMaxConstrs);
152 Q2 = AddRays(Mat->p[0],1,Pol2,NbMaxConstrs);
153
154 #ifdef DEBUG
155 fprintf(stdout, "Q1\n");
156 Polyhedron_Print(stdout,P_VALUE_FMT,Q1);
157 fprintf(stdout, "Q2\n");
158 Polyhedron_Print(stdout,P_VALUE_FMT,Q2);
159 #endif
160
161 Matrix_Free(Mat);
162 Q = DomainIntersection(Q1,Q2,NbMaxConstrs);
163
164 #ifdef DEBUG
165 fprintf(stdout, "Q\n");
166 Polyhedron_Print(stdout,P_VALUE_FMT,Q);
167 #endif
168
169 Domain_Free(Q1);
170 Domain_Free(Q2);
171
172 if (emptyQ(Q)) res = 0; /* not comparable */
173 else {
174 Q1 = DomainIntersection(Pol1,Q,NbMaxConstrs);
175 Q2 = DomainIntersection(Pol2,Q,NbMaxConstrs);
176
177 #ifdef DEBUG
178 fprintf(stdout, "Q1\n");
179 Polyhedron_Print(stdout,P_VALUE_FMT,Q1);
180 fprintf(stdout, "Q2\n");
181 Polyhedron_Print(stdout,P_VALUE_FMT,Q2);
182 #endif
183
184 /* Test if Q1 < Q2 */
185 /* Build a constraint array for >= Q1 and <= Q2 */
186 Mat = Matrix_Alloc(2,dim);
187 Vector_Set(Mat->p_Init,0,2*dim);
188
189 /* Choose a contraint from Q1 */
190 for (i=0; i<Q1->NbConstraints; i++) {
191 if (value_zero_p(Q1->Constraint[i][0])) {
192
193 /* Equality */
194 if (value_zero_p(Q1->Constraint[i][INDEX])) {
195
196 /* Ignore side constraint (they are in Q) */
197 continue;
198 }
199 else if (value_neg_p(Q1->Constraint[i][INDEX])) {
200
201 /* copy -constraint to Mat */
202 value_set_si(Mat->p[0][0],1);
203 for (k=1; k<dim; k++)
204 value_oppose(Mat->p[0][k],Q1->Constraint[i][k]);
205 }
206 else {
207
208 /* Copy constraint to Mat */
209
210 value_set_si(Mat->p[0][0],1);
211 for (k=1; k<dim; k++)
212 value_assign(Mat->p[0][k],Q1->Constraint[i][k]);
213 }
214 }
215 else if(value_neg_p(Q1->Constraint[i][INDEX])) {
216
217 /* Upper bound -- make a lower bound from it */
218 value_set_si(Mat->p[0][0],1);
219 for (k=1; k<dim; k++)
220 value_oppose(Mat->p[0][k],Q1->Constraint[i][k]);
221 }
222 else {
223
224 /* Lower or side bound -- ignore it */
225 continue;
226 }
227
228 /* Choose a constraint from Q2 */
229 for (j=0; j<Q2->NbConstraints; j++) {
230 if (value_zero_p(Q2->Constraint[j][0])) { /* equality */
231 if (value_zero_p(Q2->Constraint[j][INDEX])) {
232
233 /* Ignore side constraint (they are in Q) */
234 continue;
235 }
236 else if (value_pos_p(Q2->Constraint[j][INDEX])) {
237
238 /* Copy -constraint to Mat */
239 value_set_si(Mat->p[1][0],1);
240 for (k=1; k<dim; k++)
241 value_oppose(Mat->p[1][k],Q2->Constraint[j][k]);
242 }
243 else {
244
245 /* Copy constraint to Mat */
246 value_set_si(Mat->p[1][0],1);
247 for (k=1; k<dim; k++)
248 value_assign(Mat->p[1][k],Q2->Constraint[j][k]);
249 };
250 }
251 else if (value_pos_p(Q2->Constraint[j][INDEX])) {
252
253 /* Lower bound -- make an upper bound from it */
254 value_set_si(Mat->p[1][0],1);
255 for(k=1;k<dim;k++)
256 value_oppose(Mat->p[1][k],Q2->Constraint[j][k]);
257 }
258 else {
259
260 /* Upper or side bound -- ignore it */
261 continue;
262 };
263
264 #ifdef DEBUG
265 fprintf(stdout, "i=%d j=%d M=\n", i+1, j+1);
266 Matrix_Print(stdout,P_VALUE_FMT,Mat);
267 #endif
268
269 /* Add Mat to Q and see if anything is made */
270 Q3 = AddConstraints(Mat->p[0],2,Q,NbMaxConstrs);
271
272 #ifdef DEBUG
273 fprintf(stdout, "Q3\n");
274 Polyhedron_Print(stdout,P_VALUE_FMT,Q3);
275 #endif
276
277 if (!emptyQ(Q3)) {
278 Domain_Free(Q3);
279
280 #ifdef DEBUG
281 fprintf(stdout, "not empty\n");
282 #endif
283 res = -1;
284 goto LTQdone;
285 }
286 #ifdef DEBUG
287 fprintf(stdout,"empty\n");
288 #endif
289 Domain_Free(Q3);
290 } /* end for j */
291 } /* end for i */
292 res = 1;
293 LTQdone:
294 Matrix_Free(Mat);
295 Domain_Free(Q1);
296 Domain_Free(Q2);
297 }
298 Domain_Free(Q);
299
300 #ifdef DEBUG
301 fprintf(stdout, "res = %d\n", res);
302 #endif
303
304 return res;
305 } /* PolyhedronLTQ */
306
307 /* GaussSimplify --
308 Given Mat1, a matrix of equalities, performs Gaussian elimination.
309 Find a minimum basis, Returns the rank.
310 Mat1 is context, Mat2 is reduced in context of Mat1
311 */
312 int GaussSimplify(Matrix *Mat1,Matrix *Mat2) {
313
314 int NbRows = Mat1->NbRows;
315 int NbCols = Mat1->NbColumns;
316 int *column_index;
317 int i, j, k, n, pivot, Rank;
318 Value gcd, tmp, *cp;
319
320 column_index=(int *)malloc(NbCols * sizeof(int));
321 if (!column_index) {
322 errormsg1("GaussSimplify", "outofmem", "out of memory space\n");
323 Pol_status = 1;
324 return 0;
325 }
326
327 /* Initialize all the 'Value' variables */
328 value_init(gcd); value_init(tmp);
329
330 Rank=0;
331 for (j=0; j<NbCols; j++) { /* for each column starting at */
332 for (i=Rank; i<NbRows; i++) /* diagonal, look down to find */
333 if (value_notzero_p(Mat1->p[i][j])) /* the first non-zero entry */
334 break;
335 if (i!=NbRows) { /* was one found ? */
336 if (i!=Rank) /* was it found below the diagonal?*/
337 Vector_Exchange(Mat1->p[Rank],Mat1->p[i],NbCols);
338
339 /* Normalize the pivot row */
340 Vector_Gcd(Mat1->p[Rank],NbCols,&gcd);
341
342 /* If (gcd >= 2) */
343 value_set_si(tmp,2);
344 if (value_ge(gcd,tmp)) {
345 cp = Mat1->p[Rank];
346 for (k=0; k<NbCols; k++,cp++)
347 value_division(*cp,*cp,gcd);
348 }
349 if (value_neg_p(Mat1->p[Rank][j])) {
350 cp = Mat1->p[Rank];
351 for (k=0; k<NbCols; k++,cp++)
352 value_oppose(*cp,*cp);
353 }
354 /* End of normalize */
355
356 pivot=i;
357 for (i=0;i<NbRows;i++) /* Zero out the rest of the column */
358 if (i!=Rank) {
359 if (value_notzero_p(Mat1->p[i][j])) {
360 Value a, a1, a2;
361 value_init(a); value_init(a1); value_init(a2);
362 value_absolute(a1,Mat1->p[i][j]);
363 value_absolute(a2,Mat1->p[Rank][j]);
364 Gcd(a1,a2,&a));
365 value_division(a1,a1,a);
366 value_division(a2,a2,a);
367 value_oppose(a1,a1);
368 Vector_Combine(Mat1->p[i],Mat1->p[Rank],Mat1->p[i],a2,
369 a1,NbCols);
370 Vector_Normalize(Mat1->p[i],NbCols);
371 value_clear(a); value_clear(a1); value_clear(a2);
372 }
373 }
374 column_index[Rank]=j;
375 Rank++;
376 }
377 } /* end of Gauss elimination */
378
379 if (Mat2) { /* Mat2 is a transformation matrix (i,j->f(i,j))....
380 can't scale it because can't scale both sides of -> */
381 /* normalizes an affine transformation */
382 /* priority of forms */
383 /* 1. i' -> i (identity) */
384 /* 2. i' -> i + constant (uniform) */
385 /* 3. i' -> constant (broadcast) */
386 /* 4. i' -> j (permutation) */
387 /* 5. i' -> j + constant ( ) */
388 /* 6. i' -> i + j + constant (non-uniform) */
389 for (k=0; k<Rank; k++) {
390 j = column_index[k];
391 for (i=0; i<(Mat2->NbRows-1);i++) { /* all but the last row 0...0 1 */
392 if ((i!=j) && value_notzero_p(Mat2->p[i][j])) {
393
394 /* Remove dependency of i' on j */
395 Value a, a1, a2;
396 value_init(a); value_init(a1); value_init(a2);
397 value_absolute(a1,Mat2->p[i][j]);
398 value_absolute(a2,Mat2->p[k][j]);
399 Gcd(a1,a2,&a);
400 value_division(a1,a1,a);
401 value_division(a2,a2,a);
402 value_oppose(a1,a1);
403 if (value_one_p(a2)) {
404 Vector_Combine(Mat2->p[i],Mat1->p[k],Mat2->p[i],a2,
405 a1,NbCols);
406
407 /* Vector_Normalize(Mat2->p[i],NbCols); -- can't do T */
408 } /* otherwise, can't do it without mult lhs prod (2i,3j->...) */
409 value_clear(a); value_clear(a1); value_clear(a2);
410 }
411 else if ((i==j) && value_zero_p(Mat2->p[i][j])) {
412
413 /* 'i' does not depend on j */
414 for (n=j+1; n < (NbCols-1); n++) {
415 if (value_notzero_p(Mat2->p[i][n])) { /* i' depends on some n */
416 value_set_si(tmp,1);
417 Vector_Combine(Mat2->p[i],Mat1->p[k],Mat2->p[i],tmp,
418 tmp,NbCols);
419 break;
420 } /* if 'i' depends on just a constant, then leave it alone.*/
421 }
422 }
423 }
424 }
425
426 /* Check last row of transformation Mat2 */
427 for (j=0; j<(NbCols-1); j++)
428 if (value_notzero_p(Mat2->p[Mat2->NbRows-1][j])) {
429 errormsg1("GaussSimplify", "corrtrans", "Corrupted transformation\n");
430 break;
431 }
432
433 if (value_notone_p(Mat2->p[Mat2->NbRows-1][NbCols-1])) {
434 errormsg1("GaussSimplify", "corrtrans", "Corrupted transformation\n");
435 }
436 }
437 value_clear(gcd); value_clear(tmp);
438 free(column_index);
439 return Rank;
440 } /* GaussSimplify */
441
442 char s[128];
443
444 int main() {
445
446 Matrix *a=NULL, *b=NULL, *c, *d, *e, *f;
447 Polyhedron *A, *B, *C, *D, *last, *tmp;
448 int i, nbPol, nbMat, func;
449
450 fgets(s, 128, stdin);
451 nbPol = nbMat = 0;
452 while ((*s=='#') ||
453 ((sscanf(s, "D %d", &nbPol)<1) && (sscanf(s, "M %d", &nbMat)<1)) )
454 fgets(s, 128, stdin);
455
456 for (i=0, A=last=(Polyhedron *)0; i<nbPol; i++) {
457 a = Matrix_Read();
458 tmp = Constraints2Polyhedron(a,600);
459 Matrix_Free(a);
460 if (!last) A = last = tmp;
461 else {
462 last->next = tmp;
463 last = tmp;
464 }
465 }
466
467 if (nbMat) {
468 a = Matrix_Read();
469 }
470 fgets(s,128,stdin);
471 nbPol = nbMat = 0;
472 while ( (*s=='#') ||
473 ((sscanf(s, "D %d", &nbPol)<1) && (sscanf(s, "M %d", &nbMat)<1)) )
474 fgets(s, 128, stdin);
475
476 for (i=0, B=last=(Polyhedron *)0; i<nbPol; i++) {
477 b = Matrix_Read();
478 tmp = Constraints2Polyhedron(b,200);
479 Matrix_Free(b);
480 if (!last) B = last = tmp;
481 else {
482 last->next = tmp;
483 last = tmp;
484 }
485 }
486
487 if (nbMat) {
488 b = Matrix_Read();
489 }
490 fgets(s, 128, stdin);
491 while ((*s=='#') || (sscanf(s, "F %d", &func)<1))
492 fgets(s, 128, stdin);
493
494 switch (func) {
495 case 1:
496 C = DomainUnion(A, B, 200);
497 D = DomainConvex(C, 200);
498 d = Polyhedron2Constraints(D);
499 Matrix_Print(stdout,P_VALUE_FMT,d);
500 Matrix_Free(d);
501 Domain_Free(C);
502 Domain_Free(D);
503 break;
504 case 2:
505 D = DomainSimplify(A, B, 200);
506 d = Polyhedron2Constraints(D);
507 Matrix_Print(stdout,P_VALUE_FMT,d);
508 Matrix_Free(d);
509 Domain_Free(D);
510 break;
511 case 3:
512 a = Polyhedron2Constraints(A);
513 Matrix_Print(stdout,P_VALUE_FMT,a);
514 b = Polyhedron2Constraints(B);
515 Matrix_Print(stdout,P_VALUE_FMT,b);
516 break;
517 case 4:
518 a = Polyhedron2Rays(A);
519 Matrix_Print(stdout,P_VALUE_FMT,a);
520 break;
521 case 5:
522
523 /* a = ec , da = c , ed = 1 */
524 right_hermite(a,&c,&d,&e);
525 Matrix_Print(stdout,P_VALUE_FMT,c);
526 Matrix_Print(stdout,P_VALUE_FMT,d);
527 Matrix_Print(stdout,P_VALUE_FMT,e);
528 f = Matrix_Alloc(e->NbRows,c->NbColumns);
529 Matrix_Product(e,c,f);
530 Matrix_Print(stdout,P_VALUE_FMT,f);
531 Matrix_Free(f);
532 f = Matrix_Alloc(d->NbRows,a->NbColumns);
533 Matrix_Product(d,a,f);
534 Matrix_Print(stdout,P_VALUE_FMT,f);
535 Matrix_Free(f);
536 f = Matrix_Alloc(e->NbRows, d->NbColumns);
537 Matrix_Product(e,d,f);
538 Matrix_Print(stdout,P_VALUE_FMT,f);
539 break;
540 case 6:
541
542 /* a = ce , ad = c , de = 1 */
543 left_hermite(a,&c,&d,&e);
544 Matrix_Print(stdout,P_VALUE_FMT,c);
545 Matrix_Print(stdout,P_VALUE_FMT,d);
546 Matrix_Print(stdout,P_VALUE_FMT,e);
547 f = Matrix_Alloc(c->NbRows, e->NbColumns);
548 Matrix_Product(c,e,f);
549 Matrix_Print(stdout,P_VALUE_FMT,f);
550 Matrix_Free(f);
551 f = Matrix_Alloc(a->NbRows, d->NbColumns);
552 Matrix_Product(a,d,f);
553 Matrix_Print(stdout,P_VALUE_FMT,f);
554 Matrix_Free(f);
555 f = Matrix_Alloc(d->NbRows, e->NbColumns);
556 Matrix_Product(d,e,f);
557 Matrix_Print(stdout,P_VALUE_FMT,f);
558 break;
559 case 7:
560
561 /* Polyhedron_Print(stdout,"%5d", A); */
562 /* Matrix_Print(stdout,"%4d", b); */
563
564 C = Polyhedron_Image(A, b, 400);
565 Polyhedron_Print(stdout,P_VALUE_FMT,C);
566 break;
567 case 8:
568
569 printf("%s\n",
570 Polyhedron_Not_Empty(A,B,600) ? "Not Empty" : "Empty");
571 break;
572 case 9:
573
574 i = PolyhedronLTQ(A,B,600);
575 printf("%s\n",
576 i==-1 ? "A<B" : i==1 ? "A>B" : i==0 ? "A><B" : "error");
577 i = PolyhedronLTQ(B,A,600);
578 printf("%s\n",
579 i==-1 ? "A<B" : i==1 ? "A>B" : i==0 ? "A><B" : "error");
580 break;
581 default:
582 printf("? unknown function\n");
583 }
584
585 Domain_Free(A);
586 Domain_Free(B);
587
588 return 0;
589 }
590
591
Polyhedral library