2 This file is part of PolyLib.
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.
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.
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/>.
20 #include <polylib/polylib.h>
24 Both this software and its documentation are
26 Copyright 1993 by IRISA /Universite de Rennes I - France,
27 Copyright 1995,1996 by BYU, Provo, Utah
36 /*---------------------------------------------------------------------*/
37 /* int exist_points(pos,P,context) */
38 /* pos : index position of current loop index (0..hdim-1) */
40 /* context : context values for fixed indices */
41 /* recursive procedure, recurs for each imbriquation */
42 /* returns 1 if there exists any integer points in this polyhedron */
43 /* returns 0 if there are none */
44 /*---------------------------------------------------------------------*/
45 static int exist_points(int pos
,Polyhedron
*Pol
,Value
*context
) {
49 value_init(LB
); value_init(UB
);
50 value_init(k
); value_init(tmp
);
54 /* Problem if UB or LB is INFINITY */
55 if (lower_upper_bounds(pos
,Pol
,context
,&LB
,&UB
) !=0) {
56 errormsg1("exist_points", "infdom", "infinite domain");
63 value_set_si(context
[pos
],0);
72 value_subtract(tmp
,UB
,LB
);
73 value_increment(tmp
,tmp
);
77 return (value_pos_p(tmp
));
80 for (value_assign(k
,LB
);value_le(k
,UB
);value_increment(k
,k
)) {
82 /* insert k in context */
83 value_assign(context
[pos
],k
);
84 if (exist_points(pos
+1,Pol
->next
,context
) > 0 ) {
85 value_clear(LB
); value_clear(UB
);
86 value_clear(k
); value_clear(tmp
);
91 value_set_si(context
[pos
],0);
92 value_clear(UB
); value_clear(LB
);
93 value_clear(k
); value_clear(tmp
);
97 /*--------------------------------------------------------------*/
98 /* Test to see if there are any integral points in a polyhedron */
99 /*--------------------------------------------------------------*/
100 int Polyhedron_Not_Empty(Polyhedron
*P
,Polyhedron
*C
,int MAXRAYS
) {
106 /* Create a context vector size dim+2 and set it to all zeros */
107 context
= (Value
*) malloc((P
->Dimension
+2)*sizeof(Value
));
109 /* Initialize array 'context' */
110 for (i
=0;i
<(P
->Dimension
+2);i
++)
111 value_init(context
[i
]);
113 Vector_Set(context
,0,(P
->Dimension
+2));
115 /* Set context[P->Dimension+1] = 1 (the constant) */
116 value_set_si(context
[P
->Dimension
+1],1);
118 L
= Polyhedron_Scan(P
,C
,MAXRAYS
);
119 res
= exist_points(1,L
,context
);
122 /* Clear array 'context' */
123 for (i
=0;i
<(P
->Dimension
+2);i
++)
124 value_clear(context
[i
]);
129 /* PolyhedronLTQ ( P1, P2 ) */
130 /* P1 and P2 must be simple polyhedra */
131 /* result = 0 --> not comparable */
132 /* result = -1 --> P1 < P2 */
133 /* result = 1 --> P1 > P2 */
134 int PolyhedronLTQ (Polyhedron
*Pol1
,Polyhedron
*Pol2
,int NbMaxConstrs
) {
136 int res
, dim
, i
, j
, k
;
137 Polyhedron
*Q1
, *Q2
, *Q3
, *Q
;
142 if (Pol1
->next
|| Pol2
->next
) {
143 errormsg1("PolyhedronLTQ", "compoly", "Can only compare polyhedra");
146 if (Pol1
->Dimension
!= Pol2
->Dimension
) {
147 errormsg1("PolyhedronLTQ","diffdim","Polyhedra are not same dimension");
150 dim
= Pol1
->Dimension
+2;
153 fprintf(stdout
, "P1\n");
154 Polyhedron_Print(stdout
,P_VALUE_FMT
,Pol1
);
155 fprintf(stdout
, "P2\n");
156 Polyhedron_Print(stdout
,P_VALUE_FMT
,Pol2
);
159 /* Create the Line to add */
160 Mat
= Matrix_Alloc(1,dim
);
161 Vector_Set(Mat
->p_Init
,0,dim
);
162 value_set_si(Mat
->p
[0][INDEX
],1); /* line in INDEX dimension */
164 Q1
= AddRays(Mat
->p
[0],1,Pol1
,NbMaxConstrs
);
165 Q2
= AddRays(Mat
->p
[0],1,Pol2
,NbMaxConstrs
);
168 fprintf(stdout
, "Q1\n");
169 Polyhedron_Print(stdout
,P_VALUE_FMT
,Q1
);
170 fprintf(stdout
, "Q2\n");
171 Polyhedron_Print(stdout
,P_VALUE_FMT
,Q2
);
175 Q
= DomainIntersection(Q1
,Q2
,NbMaxConstrs
);
178 fprintf(stdout
, "Q\n");
179 Polyhedron_Print(stdout
,P_VALUE_FMT
,Q
);
185 if (emptyQ(Q
)) res
= 0; /* not comparable */
187 Q1
= DomainIntersection(Pol1
,Q
,NbMaxConstrs
);
188 Q2
= DomainIntersection(Pol2
,Q
,NbMaxConstrs
);
191 fprintf(stdout
, "Q1\n");
192 Polyhedron_Print(stdout
,P_VALUE_FMT
,Q1
);
193 fprintf(stdout
, "Q2\n");
194 Polyhedron_Print(stdout
,P_VALUE_FMT
,Q2
);
197 /* Test if Q1 < Q2 */
198 /* Build a constraint array for >= Q1 and <= Q2 */
199 Mat
= Matrix_Alloc(2,dim
);
200 Vector_Set(Mat
->p_Init
,0,2*dim
);
202 /* Choose a contraint from Q1 */
203 for (i
=0; i
<Q1
->NbConstraints
; i
++) {
204 if (value_zero_p(Q1
->Constraint
[i
][0])) {
207 if (value_zero_p(Q1
->Constraint
[i
][INDEX
])) {
209 /* Ignore side constraint (they are in Q) */
212 else if (value_neg_p(Q1
->Constraint
[i
][INDEX
])) {
214 /* copy -constraint to Mat */
215 value_set_si(Mat
->p
[0][0],1);
216 for (k
=1; k
<dim
; k
++)
217 value_oppose(Mat
->p
[0][k
],Q1
->Constraint
[i
][k
]);
221 /* Copy constraint to Mat */
223 value_set_si(Mat
->p
[0][0],1);
224 for (k
=1; k
<dim
; k
++)
225 value_assign(Mat
->p
[0][k
],Q1
->Constraint
[i
][k
]);
228 else if(value_neg_p(Q1
->Constraint
[i
][INDEX
])) {
230 /* Upper bound -- make a lower bound from it */
231 value_set_si(Mat
->p
[0][0],1);
232 for (k
=1; k
<dim
; k
++)
233 value_oppose(Mat
->p
[0][k
],Q1
->Constraint
[i
][k
]);
237 /* Lower or side bound -- ignore it */
241 /* Choose a constraint from Q2 */
242 for (j
=0; j
<Q2
->NbConstraints
; j
++) {
243 if (value_zero_p(Q2
->Constraint
[j
][0])) { /* equality */
244 if (value_zero_p(Q2
->Constraint
[j
][INDEX
])) {
246 /* Ignore side constraint (they are in Q) */
249 else if (value_pos_p(Q2
->Constraint
[j
][INDEX
])) {
251 /* Copy -constraint to Mat */
252 value_set_si(Mat
->p
[1][0],1);
253 for (k
=1; k
<dim
; k
++)
254 value_oppose(Mat
->p
[1][k
],Q2
->Constraint
[j
][k
]);
258 /* Copy constraint to Mat */
259 value_set_si(Mat
->p
[1][0],1);
260 for (k
=1; k
<dim
; k
++)
261 value_assign(Mat
->p
[1][k
],Q2
->Constraint
[j
][k
]);
264 else if (value_pos_p(Q2
->Constraint
[j
][INDEX
])) {
266 /* Lower bound -- make an upper bound from it */
267 value_set_si(Mat
->p
[1][0],1);
269 value_oppose(Mat
->p
[1][k
],Q2
->Constraint
[j
][k
]);
273 /* Upper or side bound -- ignore it */
278 fprintf(stdout
, "i=%d j=%d M=\n", i
+1, j
+1);
279 Matrix_Print(stdout
,P_VALUE_FMT
,Mat
);
282 /* Add Mat to Q and see if anything is made */
283 Q3
= AddConstraints(Mat
->p
[0],2,Q
,NbMaxConstrs
);
286 fprintf(stdout
, "Q3\n");
287 Polyhedron_Print(stdout
,P_VALUE_FMT
,Q3
);
294 fprintf(stdout
, "not empty\n");
300 fprintf(stdout
,"empty\n");
314 fprintf(stdout
, "res = %d\n", res
);
318 } /* PolyhedronLTQ */
321 Given Mat1, a matrix of equalities, performs Gaussian elimination.
322 Find a minimum basis, Returns the rank.
323 Mat1 is context, Mat2 is reduced in context of Mat1
325 int GaussSimplify(Matrix
*Mat1
,Matrix
*Mat2
) {
327 int NbRows
= Mat1
->NbRows
;
328 int NbCols
= Mat1
->NbColumns
;
330 int i
, j
, k
, n
, pivot
, Rank
;
333 column_index
=(int *)malloc(NbCols
* sizeof(int));
335 errormsg1("GaussSimplify", "outofmem", "out of memory space\n");
340 /* Initialize all the 'Value' variables */
341 value_init(gcd
); value_init(tmp
);
344 for (j
=0; j
<NbCols
; j
++) { /* for each column starting at */
345 for (i
=Rank
; i
<NbRows
; i
++) /* diagonal, look down to find */
346 if (value_notzero_p(Mat1
->p
[i
][j
])) /* the first non-zero entry */
348 if (i
!=NbRows
) { /* was one found ? */
349 if (i
!=Rank
) /* was it found below the diagonal?*/
350 Vector_Exchange(Mat1
->p
[Rank
],Mat1
->p
[i
],NbCols
);
352 /* Normalize the pivot row */
353 Vector_Gcd(Mat1
->p
[Rank
],NbCols
,&gcd
);
357 if (value_ge(gcd
,tmp
)) {
359 for (k
=0; k
<NbCols
; k
++,cp
++)
360 value_division(*cp
,*cp
,gcd
);
362 if (value_neg_p(Mat1
->p
[Rank
][j
])) {
364 for (k
=0; k
<NbCols
; k
++,cp
++)
365 value_oppose(*cp
,*cp
);
367 /* End of normalize */
370 for (i
=0;i
<NbRows
;i
++) /* Zero out the rest of the column */
372 if (value_notzero_p(Mat1
->p
[i
][j
])) {
374 value_init(a
); value_init(a1
); value_init(a2
);
375 value_absolute(a1
,Mat1
->p
[i
][j
]);
376 value_absolute(a2
,Mat1
->p
[Rank
][j
]);
377 value_gcd(a
, a1
, a2
));
378 value_divexact(a1
, a1
, a
);
379 value_divexact(a2
, a2
, a
);
381 Vector_Combine(Mat1
->p
[i
],Mat1
->p
[Rank
],Mat1
->p
[i
],a2
,
383 Vector_Normalize(Mat1
->p
[i
],NbCols
);
384 value_clear(a
); value_clear(a1
); value_clear(a2
);
387 column_index
[Rank
]=j
;
390 } /* end of Gauss elimination */
392 if (Mat2
) { /* Mat2 is a transformation matrix (i,j->f(i,j))....
393 can't scale it because can't scale both sides of -> */
394 /* normalizes an affine transformation */
395 /* priority of forms */
396 /* 1. i' -> i (identity) */
397 /* 2. i' -> i + constant (uniform) */
398 /* 3. i' -> constant (broadcast) */
399 /* 4. i' -> j (permutation) */
400 /* 5. i' -> j + constant ( ) */
401 /* 6. i' -> i + j + constant (non-uniform) */
402 for (k
=0; k
<Rank
; k
++) {
404 for (i
=0; i
<(Mat2
->NbRows
-1);i
++) { /* all but the last row 0...0 1 */
405 if ((i
!=j
) && value_notzero_p(Mat2
->p
[i
][j
])) {
407 /* Remove dependency of i' on j */
409 value_init(a
); value_init(a1
); value_init(a2
);
410 value_absolute(a1
,Mat2
->p
[i
][j
]);
411 value_absolute(a2
,Mat2
->p
[k
][j
]);
412 value_gcd(a
, a1
, a2
);
413 value_divexact(a1
, a1
, a
);
414 value_divexact(a2
, a2
, a
);
416 if (value_one_p(a2
)) {
417 Vector_Combine(Mat2
->p
[i
],Mat1
->p
[k
],Mat2
->p
[i
],a2
,
420 /* Vector_Normalize(Mat2->p[i],NbCols); -- can't do T */
421 } /* otherwise, can't do it without mult lhs prod (2i,3j->...) */
422 value_clear(a
); value_clear(a1
); value_clear(a2
);
424 else if ((i
==j
) && value_zero_p(Mat2
->p
[i
][j
])) {
426 /* 'i' does not depend on j */
427 for (n
=j
+1; n
< (NbCols
-1); n
++) {
428 if (value_notzero_p(Mat2
->p
[i
][n
])) { /* i' depends on some n */
430 Vector_Combine(Mat2
->p
[i
],Mat1
->p
[k
],Mat2
->p
[i
],tmp
,
433 } /* if 'i' depends on just a constant, then leave it alone.*/
439 /* Check last row of transformation Mat2 */
440 for (j
=0; j
<(NbCols
-1); j
++)
441 if (value_notzero_p(Mat2
->p
[Mat2
->NbRows
-1][j
])) {
442 errormsg1("GaussSimplify", "corrtrans", "Corrupted transformation\n");
446 if (value_notone_p(Mat2
->p
[Mat2
->NbRows
-1][NbCols
-1])) {
447 errormsg1("GaussSimplify", "corrtrans", "Corrupted transformation\n");
450 value_clear(gcd
); value_clear(tmp
);
453 } /* GaussSimplify */
459 Matrix
*a
=NULL
, *b
=NULL
, *c
, *d
, *e
, *f
;
460 Polyhedron
*A
, *B
, *C
, *D
, *last
, *tmp
;
461 int i
, nbPol
, nbMat
, func
;
463 fgets(s
, 128, stdin
);
466 ((sscanf(s
, "D %d", &nbPol
)<1) && (sscanf(s
, "M %d", &nbMat
)<1)) )
467 fgets(s
, 128, stdin
);
469 for (i
=0, A
=last
=(Polyhedron
*)0; i
<nbPol
; i
++) {
471 tmp
= Constraints2Polyhedron(a
,600);
473 if (!last
) A
= last
= tmp
;
486 ((sscanf(s
, "D %d", &nbPol
)<1) && (sscanf(s
, "M %d", &nbMat
)<1)) )
487 fgets(s
, 128, stdin
);
489 for (i
=0, B
=last
=(Polyhedron
*)0; i
<nbPol
; i
++) {
491 tmp
= Constraints2Polyhedron(b
,200);
493 if (!last
) B
= last
= tmp
;
503 fgets(s
, 128, stdin
);
504 while ((*s
=='#') || (sscanf(s
, "F %d", &func
)<1))
505 fgets(s
, 128, stdin
);
509 C
= DomainUnion(A
, B
, 200);
510 D
= DomainConvex(C
, 200);
511 d
= Polyhedron2Constraints(D
);
512 Matrix_Print(stdout
,P_VALUE_FMT
,d
);
518 D
= DomainSimplify(A
, B
, 200);
519 d
= Polyhedron2Constraints(D
);
520 Matrix_Print(stdout
,P_VALUE_FMT
,d
);
525 a
= Polyhedron2Constraints(A
);
526 Matrix_Print(stdout
,P_VALUE_FMT
,a
);
527 b
= Polyhedron2Constraints(B
);
528 Matrix_Print(stdout
,P_VALUE_FMT
,b
);
531 a
= Polyhedron2Rays(A
);
532 Matrix_Print(stdout
,P_VALUE_FMT
,a
);
536 /* a = ec , da = c , ed = 1 */
537 right_hermite(a
,&c
,&d
,&e
);
538 Matrix_Print(stdout
,P_VALUE_FMT
,c
);
539 Matrix_Print(stdout
,P_VALUE_FMT
,d
);
540 Matrix_Print(stdout
,P_VALUE_FMT
,e
);
541 f
= Matrix_Alloc(e
->NbRows
,c
->NbColumns
);
542 Matrix_Product(e
,c
,f
);
543 Matrix_Print(stdout
,P_VALUE_FMT
,f
);
545 f
= Matrix_Alloc(d
->NbRows
,a
->NbColumns
);
546 Matrix_Product(d
,a
,f
);
547 Matrix_Print(stdout
,P_VALUE_FMT
,f
);
549 f
= Matrix_Alloc(e
->NbRows
, d
->NbColumns
);
550 Matrix_Product(e
,d
,f
);
551 Matrix_Print(stdout
,P_VALUE_FMT
,f
);
555 /* a = ce , ad = c , de = 1 */
556 left_hermite(a
,&c
,&d
,&e
);
557 Matrix_Print(stdout
,P_VALUE_FMT
,c
);
558 Matrix_Print(stdout
,P_VALUE_FMT
,d
);
559 Matrix_Print(stdout
,P_VALUE_FMT
,e
);
560 f
= Matrix_Alloc(c
->NbRows
, e
->NbColumns
);
561 Matrix_Product(c
,e
,f
);
562 Matrix_Print(stdout
,P_VALUE_FMT
,f
);
564 f
= Matrix_Alloc(a
->NbRows
, d
->NbColumns
);
565 Matrix_Product(a
,d
,f
);
566 Matrix_Print(stdout
,P_VALUE_FMT
,f
);
568 f
= Matrix_Alloc(d
->NbRows
, e
->NbColumns
);
569 Matrix_Product(d
,e
,f
);
570 Matrix_Print(stdout
,P_VALUE_FMT
,f
);
574 /* Polyhedron_Print(stdout,"%5d", A); */
575 /* Matrix_Print(stdout,"%4d", b); */
577 C
= Polyhedron_Image(A
, b
, 400);
578 Polyhedron_Print(stdout
,P_VALUE_FMT
,C
);
583 Polyhedron_Not_Empty(A
,B
,600) ? "Not Empty" : "Empty");
587 i
= PolyhedronLTQ(A
,B
,600);
589 i
==-1 ? "A<B" : i
==1 ? "A>B" : i
==0 ? "A><B" : "error");
590 i
= PolyhedronLTQ(B
,A
,600);
592 i
==-1 ? "A<B" : i
==1 ? "A>B" : i
==0 ? "A><B" : "error");
595 printf("? unknown function\n");