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