3 #include <polylib/polylib.h>
7 Both this software and its documentation are
9 Copyright 1993 by IRISA /Universite de Rennes I - France,
10 Copyright 1995,1996 by BYU, Provo, Utah
19 /*---------------------------------------------------------------------*/
20 /* int exist_points(pos,P,context) */
21 /* pos : index position of current loop index (0..hdim-1) */
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
) {
32 value_init(LB
); value_init(UB
);
33 value_init(k
); value_init(tmp
);
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");
46 value_set_si(context
[pos
],0);
55 value_subtract(tmp
,UB
,LB
);
56 value_increment(tmp
,tmp
);
60 return (value_pos_p(tmp
));
63 for (value_assign(k
,LB
);value_le(k
,UB
);value_increment(k
,k
)) {
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
);
74 value_set_si(context
[pos
],0);
75 value_clear(UB
); value_clear(LB
);
76 value_clear(k
); value_clear(tmp
);
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
) {
89 /* Create a context vector size dim+2 and set it to all zeros */
90 context
= (Value
*) malloc((P
->Dimension
+2)*sizeof(Value
));
92 /* Initialize array 'context' */
93 for (i
=0;i
<(P
->Dimension
+2);i
++)
94 value_init(context
[i
]);
96 Vector_Set(context
,0,(P
->Dimension
+2));
98 /* Set context[P->Dimension+1] = 1 (the constant) */
99 value_set_si(context
[P
->Dimension
+1],1);
101 L
= Polyhedron_Scan(P
,C
,MAXRAYS
);
102 res
= exist_points(1,L
,context
);
105 /* Clear array 'context' */
106 for (i
=0;i
<(P
->Dimension
+2);i
++)
107 value_clear(context
[i
]);
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
) {
119 int res
, dim
, i
, j
, k
;
120 Polyhedron
*Q1
, *Q2
, *Q3
, *Q
;
125 if (Pol1
->next
|| Pol2
->next
) {
126 errormsg1("PolyhedronLTQ", "compoly", "Can only compare polyhedra");
129 if (Pol1
->Dimension
!= Pol2
->Dimension
) {
130 errormsg1("PolyhedronLTQ","diffdim","Polyhedra are not same dimension");
133 dim
= Pol1
->Dimension
+2;
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
);
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 */
147 Q1
= AddRays(Mat
->p
[0],1,Pol1
,NbMaxConstrs
);
148 Q2
= AddRays(Mat
->p
[0],1,Pol2
,NbMaxConstrs
);
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
);
158 Q
= DomainIntersection(Q1
,Q2
,NbMaxConstrs
);
161 fprintf(stdout
, "Q\n");
162 Polyhedron_Print(stdout
,P_VALUE_FMT
,Q
);
168 if (emptyQ(Q
)) res
= 0; /* not comparable */
170 Q1
= DomainIntersection(Pol1
,Q
,NbMaxConstrs
);
171 Q2
= DomainIntersection(Pol2
,Q
,NbMaxConstrs
);
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
);
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
);
185 /* Choose a contraint from Q1 */
186 for (i
=0; i
<Q1
->NbConstraints
; i
++) {
187 if (value_zero_p(Q1
->Constraint
[i
][0])) {
190 if (value_zero_p(Q1
->Constraint
[i
][INDEX
])) {
192 /* Ignore side constraint (they are in Q) */
195 else if (value_neg_p(Q1
->Constraint
[i
][INDEX
])) {
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
]);
204 /* Copy constraint to Mat */
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
]);
211 else if(value_neg_p(Q1
->Constraint
[i
][INDEX
])) {
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
]);
220 /* Lower or side bound -- ignore it */
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
])) {
229 /* Ignore side constraint (they are in Q) */
232 else if (value_pos_p(Q2
->Constraint
[j
][INDEX
])) {
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
]);
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
]);
247 else if (value_pos_p(Q2
->Constraint
[j
][INDEX
])) {
249 /* Lower bound -- make an upper bound from it */
250 value_set_si(Mat
->p
[1][0],1);
252 value_oppose(Mat
->p
[1][k
],Q2
->Constraint
[j
][k
]);
256 /* Upper or side bound -- ignore it */
261 fprintf(stdout
, "i=%d j=%d M=\n", i
+1, j
+1);
262 Matrix_Print(stdout
,P_VALUE_FMT
,Mat
);
265 /* Add Mat to Q and see if anything is made */
266 Q3
= AddConstraints(Mat
->p
[0],2,Q
,NbMaxConstrs
);
269 fprintf(stdout
, "Q3\n");
270 Polyhedron_Print(stdout
,P_VALUE_FMT
,Q3
);
277 fprintf(stdout
, "not empty\n");
283 fprintf(stdout
,"empty\n");
297 fprintf(stdout
, "res = %d\n", res
);
301 } /* PolyhedronLTQ */
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
308 int GaussSimplify(Matrix
*Mat1
,Matrix
*Mat2
) {
310 int NbRows
= Mat1
->NbRows
;
311 int NbCols
= Mat1
->NbColumns
;
313 int i
, j
, k
, n
, pivot
, Rank
;
316 column_index
=(int *)malloc(NbCols
* sizeof(int));
318 errormsg1("GaussSimplify", "outofmem", "out of memory space\n");
323 /* Initialize all the 'Value' variables */
324 value_init(gcd
); value_init(tmp
);
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 */
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
);
335 /* Normalize the pivot row */
336 Vector_Gcd(Mat1
->p
[Rank
],NbCols
,&gcd
);
340 if (value_ge(gcd
,tmp
)) {
342 for (k
=0; k
<NbCols
; k
++,cp
++)
343 value_division(*cp
,*cp
,gcd
);
345 if (value_neg_p(Mat1
->p
[Rank
][j
])) {
347 for (k
=0; k
<NbCols
; k
++,cp
++)
348 value_oppose(*cp
,*cp
);
350 /* End of normalize */
353 for (i
=0;i
<NbRows
;i
++) /* Zero out the rest of the column */
355 if (value_notzero_p(Mat1
->p
[i
][j
])) {
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
);
364 Vector_Combine(Mat1
->p
[i
],Mat1
->p
[Rank
],Mat1
->p
[i
],a2
,
366 Vector_Normalize(Mat1
->p
[i
],NbCols
);
367 value_clear(a
); value_clear(a1
); value_clear(a2
);
370 column_index
[Rank
]=j
;
373 } /* end of Gauss elimination */
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
++) {
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
])) {
390 /* Remove dependency of i' on j */
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
);
399 if (value_one_p(a2
)) {
400 Vector_Combine(Mat2
->p
[i
],Mat1
->p
[k
],Mat2
->p
[i
],a2
,
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
);
407 else if ((i
==j
) && value_zero_p(Mat2
->p
[i
][j
])) {
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 */
413 Vector_Combine(Mat2
->p
[i
],Mat1
->p
[k
],Mat2
->p
[i
],tmp
,
416 } /* if 'i' depends on just a constant, then leave it alone.*/
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");
429 if (value_notone_p(Mat2
->p
[Mat2
->NbRows
-1][NbCols
-1])) {
430 errormsg1("GaussSimplify", "corrtrans", "Corrupted transformation\n");
433 value_clear(gcd
); value_clear(tmp
);
436 } /* GaussSimplify */
442 Matrix
*a
=NULL
, *b
=NULL
, *c
, *d
, *e
, *f
;
443 Polyhedron
*A
, *B
, *C
, *D
, *last
, *tmp
;
444 int i
, nbPol
, nbMat
, func
;
446 fgets(s
, 128, stdin
);
449 ((sscanf(s
, "D %d", &nbPol
)<1) && (sscanf(s
, "M %d", &nbMat
)<1)) )
450 fgets(s
, 128, stdin
);
452 for (i
=0, A
=last
=(Polyhedron
*)0; i
<nbPol
; i
++) {
454 tmp
= Constraints2Polyhedron(a
,600);
456 if (!last
) A
= last
= tmp
;
469 ((sscanf(s
, "D %d", &nbPol
)<1) && (sscanf(s
, "M %d", &nbMat
)<1)) )
470 fgets(s
, 128, stdin
);
472 for (i
=0, B
=last
=(Polyhedron
*)0; i
<nbPol
; i
++) {
474 tmp
= Constraints2Polyhedron(b
,200);
476 if (!last
) B
= last
= tmp
;
486 fgets(s
, 128, stdin
);
487 while ((*s
=='#') || (sscanf(s
, "F %d", &func
)<1))
488 fgets(s
, 128, stdin
);
492 C
= DomainUnion(A
, B
, 200);
493 D
= DomainConvex(C
, 200);
494 d
= Polyhedron2Constraints(D
);
495 Matrix_Print(stdout
,P_VALUE_FMT
,d
);
501 D
= DomainSimplify(A
, B
, 200);
502 d
= Polyhedron2Constraints(D
);
503 Matrix_Print(stdout
,P_VALUE_FMT
,d
);
508 a
= Polyhedron2Constraints(A
);
509 Matrix_Print(stdout
,P_VALUE_FMT
,a
);
510 b
= Polyhedron2Constraints(B
);
511 Matrix_Print(stdout
,P_VALUE_FMT
,b
);
514 a
= Polyhedron2Rays(A
);
515 Matrix_Print(stdout
,P_VALUE_FMT
,a
);
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
);
528 f
= Matrix_Alloc(d
->NbRows
,a
->NbColumns
);
529 Matrix_Product(d
,a
,f
);
530 Matrix_Print(stdout
,P_VALUE_FMT
,f
);
532 f
= Matrix_Alloc(e
->NbRows
, d
->NbColumns
);
533 Matrix_Product(e
,d
,f
);
534 Matrix_Print(stdout
,P_VALUE_FMT
,f
);
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
);
547 f
= Matrix_Alloc(a
->NbRows
, d
->NbColumns
);
548 Matrix_Product(a
,d
,f
);
549 Matrix_Print(stdout
,P_VALUE_FMT
,f
);
551 f
= Matrix_Alloc(d
->NbRows
, e
->NbColumns
);
552 Matrix_Product(d
,e
,f
);
553 Matrix_Print(stdout
,P_VALUE_FMT
,f
);
557 /* Polyhedron_Print(stdout,"%5d", A); */
558 /* Matrix_Print(stdout,"%4d", b); */
560 C
= Polyhedron_Image(A
, b
, 400);
561 Polyhedron_Print(stdout
,P_VALUE_FMT
,C
);
566 Polyhedron_Not_Empty(A
,B
,600) ? "Not Empty" : "Empty");
570 i
= PolyhedronLTQ(A
,B
,600);
572 i
==-1 ? "A<B" : i
==1 ? "A>B" : i
==0 ? "A><B" : "error");
573 i
= PolyhedronLTQ(B
,A
,600);
575 i
==-1 ? "A<B" : i
==1 ? "A>B" : i
==0 ? "A><B" : "error");
578 printf("? unknown function\n");