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