1 #include <polylib/polylibgmp.h>
6 #ifndef HAVE_ENUMERATE4
7 #define Polyhedron_Enumerate(a,b,c,d) Polyhedron_Enumerate(a,b,c)
10 void manual_count(Polyhedron
*P
, Value
* result
)
12 Polyhedron
*U
= Universe_Polyhedron(0);
13 Enumeration
*en
= Polyhedron_Enumerate(P
,U
,1024,NULL
);
14 Value
*v
= compute_poly(en
,NULL
);
15 value_assign(*result
, *v
);
22 #include "ev_operations.h"
26 /* Return random value between 0 and max-1 inclusive
28 int random_int(int max
) {
29 return (int) (((double)(max
))*rand()/(RAND_MAX
+1.0));
32 /* Inplace polarization
34 void Polyhedron_Polarize(Polyhedron
*P
)
36 unsigned NbRows
= P
->NbConstraints
+ P
->NbRays
;
40 q
= (Value
**)malloc(NbRows
* sizeof(Value
*));
42 for (i
= 0; i
< P
->NbRays
; ++i
)
44 for (; i
< NbRows
; ++i
)
45 q
[i
] = P
->Constraint
[i
-P
->NbRays
];
46 P
->NbConstraints
= NbRows
- P
->NbConstraints
;
47 P
->NbRays
= NbRows
- P
->NbRays
;
50 P
->Ray
= q
+ P
->NbConstraints
;
54 * Rather general polar
55 * We can optimize it significantly if we assume that
58 * Also, we calculate the polar as defined in Schrijver
59 * The opposite should probably work as well and would
60 * eliminate the need for multiplying by -1
62 Polyhedron
* Polyhedron_Polar(Polyhedron
*P
, unsigned NbMaxRays
)
66 unsigned dim
= P
->Dimension
+ 2;
67 Matrix
*M
= Matrix_Alloc(P
->NbRays
, dim
);
71 value_set_si(mone
, -1);
72 for (i
= 0; i
< P
->NbRays
; ++i
) {
73 Vector_Scale(P
->Ray
[i
], M
->p
[i
], mone
, dim
);
74 value_multiply(M
->p
[i
][0], M
->p
[i
][0], mone
);
75 value_multiply(M
->p
[i
][dim
-1], M
->p
[i
][dim
-1], mone
);
77 P
= Constraints2Polyhedron(M
, NbMaxRays
);
85 * Returns the supporting cone of P at the vertex with index v
87 Polyhedron
* supporting_cone(Polyhedron
*P
, int v
)
92 unsigned char *supporting
= (unsigned char *)malloc(P
->NbConstraints
);
93 unsigned dim
= P
->Dimension
+ 2;
95 assert(v
>=0 && v
< P
->NbRays
);
96 assert(value_pos_p(P
->Ray
[v
][dim
-1]));
100 for (i
= 0, n
= 0; i
< P
->NbConstraints
; ++i
) {
101 Inner_Product(P
->Constraint
[i
] + 1, P
->Ray
[v
] + 1, dim
- 1, &tmp
);
102 if ((supporting
[i
] = value_zero_p(tmp
)))
105 assert(n
>= dim
- 2);
107 M
= Matrix_Alloc(n
, dim
);
109 for (i
= 0, j
= 0; i
< P
->NbConstraints
; ++i
)
111 value_set_si(M
->p
[j
][dim
-1], 0);
112 Vector_Copy(P
->Constraint
[i
], M
->p
[j
++], dim
-1);
115 P
= Constraints2Polyhedron(M
, P
->NbRays
+1);
121 void value_lcm(Value i
, Value j
, Value
* lcm
)
125 value_multiply(aux
,i
,j
);
127 value_division(*lcm
,aux
,*lcm
);
131 Polyhedron
* supporting_cone_p(Polyhedron
*P
, Param_Vertices
*v
)
134 Value lcm
, tmp
, tmp2
;
135 unsigned char *supporting
= (unsigned char *)malloc(P
->NbConstraints
);
136 unsigned dim
= P
->Dimension
+ 2;
137 unsigned nparam
= v
->Vertex
->NbColumns
- 2;
138 unsigned nvar
= dim
- nparam
- 2;
143 row
= Vector_Alloc(nparam
+1);
148 value_set_si(lcm
, 1);
149 for (i
= 0, n
= 0; i
< P
->NbConstraints
; ++i
) {
150 Vector_Set(row
->p
, 0, nparam
+1);
151 for (j
= 0 ; j
< nvar
; ++j
) {
152 value_set_si(tmp
, 1);
153 value_assign(tmp2
, P
->Constraint
[i
][j
+1]);
154 if (value_ne(lcm
, v
->Vertex
->p
[j
][nparam
+1])) {
155 value_assign(tmp
, lcm
);
156 value_lcm(lcm
, v
->Vertex
->p
[j
][nparam
+1], &lcm
);
157 value_division(tmp
, lcm
, tmp
);
158 value_multiply(tmp2
, tmp2
, lcm
);
159 value_division(tmp2
, tmp2
, v
->Vertex
->p
[j
][nparam
+1]);
161 Vector_Combine(row
->p
, v
->Vertex
->p
[j
], row
->p
,
162 tmp
, tmp2
, nparam
+1);
164 value_set_si(tmp
, 1);
165 Vector_Combine(row
->p
, P
->Constraint
[i
]+1+nvar
, row
->p
, tmp
, lcm
, nparam
+1);
166 for (j
= 0; j
< nparam
+1; ++j
)
167 if (value_notzero_p(row
->p
[j
]))
169 if ((supporting
[i
] = (j
== nparam
+ 1)))
177 M
= Matrix_Alloc(n
, nvar
+2);
179 for (i
= 0, j
= 0; i
< P
->NbConstraints
; ++i
)
181 value_set_si(M
->p
[j
][nvar
+1], 0);
182 Vector_Copy(P
->Constraint
[i
], M
->p
[j
++], nvar
+1);
185 P
= Constraints2Polyhedron(M
, P
->NbRays
+1);
191 Polyhedron
* triangularize_cone(Polyhedron
*P
, unsigned NbMaxCons
)
193 const static int MAX_TRY
=10;
196 unsigned dim
= P
->Dimension
;
197 Matrix
*M
= Matrix_Alloc(P
->NbRays
+1, dim
+3);
199 Polyhedron
*L
, *R
, *T
;
200 assert(P
->NbEq
== 0);
205 Vector_Set(M
->p
[0]+1, 0, dim
+1);
206 value_set_si(M
->p
[0][0], 1);
207 value_set_si(M
->p
[0][dim
+2], 1);
208 Vector_Set(M
->p
[P
->NbRays
]+1, 0, dim
+2);
209 value_set_si(M
->p
[P
->NbRays
][0], 1);
210 value_set_si(M
->p
[P
->NbRays
][dim
+1], 1);
212 for (i
= 0, r
= 1; i
< P
->NbRays
; ++i
) {
213 if (value_notzero_p(P
->Ray
[i
][dim
+1]))
215 Vector_Copy(P
->Ray
[i
], M
->p
[r
], dim
+1);
216 Inner_Product(M
->p
[r
]+1, M
->p
[r
]+1, dim
, &tmp
);
217 value_assign(M
->p
[r
][dim
+1], tmp
);
218 value_set_si(M
->p
[r
][dim
+2], 0);
223 L
= Rays2Polyhedron(M3
, NbMaxCons
);
226 M2
= Matrix_Alloc(dim
+1, dim
+2);
231 /* Usually R should still be 0 */
234 for (r
= 1; r
< P
->NbRays
; ++r
) {
235 value_set_si(M
->p
[r
][dim
+1], random_int((t
+1)*dim
)+1);
238 L
= Rays2Polyhedron(M3
, NbMaxCons
);
242 assert(t
<= MAX_TRY
);
247 for (i
= 0; i
< L
->NbConstraints
; ++i
) {
248 if (value_negz_p(L
->Constraint
[i
][dim
+1]))
250 if (value_notzero_p(L
->Constraint
[i
][dim
+2]))
252 for (j
= 1, r
= 1; j
< M
->NbRows
; ++j
) {
253 Inner_Product(M
->p
[j
]+1, L
->Constraint
[i
]+1, dim
+1, &tmp
);
254 if (value_notzero_p(tmp
))
258 Vector_Copy(M
->p
[j
]+1, M2
->p
[r
]+1, dim
);
259 value_set_si(M2
->p
[r
][0], 1);
260 value_set_si(M2
->p
[r
][dim
+1], 0);
264 Vector_Set(M2
->p
[0]+1, 0, dim
);
265 value_set_si(M2
->p
[0][0], 1);
266 value_set_si(M2
->p
[0][dim
+1], 1);
267 T
= Rays2Polyhedron(M2
, P
->NbConstraints
+1);
281 void check_triangulization(Polyhedron
*P
, Polyhedron
*T
)
283 Polyhedron
*C
, *D
, *E
, *F
, *G
, *U
;
284 for (C
= T
; C
; C
= C
->next
) {
288 U
= DomainConvex(DomainUnion(U
, C
, 100), 100);
289 for (D
= C
->next
; D
; D
= D
->next
) {
294 E
= DomainIntersection(C
, D
, 600);
295 assert(E
->NbRays
== 0 || E
->NbEq
>= 1);
301 assert(PolyhedronIncludes(U
, P
));
302 assert(PolyhedronIncludes(P
, U
));
305 void Euclid(Value a
, Value b
, Value
*x
, Value
*y
, Value
*g
)
307 Value c
, d
, e
, f
, tmp
;
314 value_absolute(c
, a
);
315 value_absolute(d
, b
);
318 while(value_pos_p(d
)) {
319 value_division(tmp
, c
, d
);
320 value_multiply(tmp
, tmp
, f
);
321 value_substract(e
, e
, tmp
);
322 value_division(tmp
, c
, d
);
323 value_multiply(tmp
, tmp
, d
);
324 value_substract(c
, c
, tmp
);
331 else if (value_pos_p(a
))
333 else value_oppose(*x
, e
);
337 value_multiply(tmp
, a
, *x
);
338 value_substract(tmp
, c
, tmp
);
339 value_division(*y
, tmp
, b
);
348 Matrix
* unimodular_complete(Vector
*row
)
350 Value g
, b
, c
, old
, tmp
;
359 m
= Matrix_Alloc(row
->Size
, row
->Size
);
360 for (j
= 0; j
< row
->Size
; ++j
) {
361 value_assign(m
->p
[0][j
], row
->p
[j
]);
363 value_assign(g
, row
->p
[0]);
364 for (i
= 1; value_zero_p(g
) && i
< row
->Size
; ++i
) {
365 for (j
= 0; j
< row
->Size
; ++j
) {
367 value_set_si(m
->p
[i
][j
], 1);
369 value_set_si(m
->p
[i
][j
], 0);
371 value_assign(g
, row
->p
[i
]);
373 for (; i
< row
->Size
; ++i
) {
374 value_assign(old
, g
);
375 Euclid(old
, row
->p
[i
], &c
, &b
, &g
);
377 for (j
= 0; j
< row
->Size
; ++j
) {
379 value_multiply(tmp
, row
->p
[j
], b
);
380 value_division(m
->p
[i
][j
], tmp
, old
);
382 value_assign(m
->p
[i
][j
], c
);
384 value_set_si(m
->p
[i
][j
], 0);
396 * Returns a full-dimensional polyhedron with the same number
397 * of integer points as P
399 Polyhedron
*remove_equalities(Polyhedron
*P
)
403 Polyhedron
*p
= Polyhedron_Copy(P
), *q
;
404 unsigned dim
= p
->Dimension
;
409 while (p
->NbEq
> 0) {
411 Vector_Gcd(p
->Constraint
[0]+1, dim
+1, &g
);
412 Vector_AntiScale(p
->Constraint
[0]+1, p
->Constraint
[0]+1, g
, dim
+1);
413 Vector_Gcd(p
->Constraint
[0]+1, dim
, &g
);
414 if (value_notone_p(g
) && value_notmone_p(g
)) {
416 p
= Empty_Polyhedron(0);
419 v
= Vector_Alloc(dim
);
420 Vector_Copy(p
->Constraint
[0]+1, v
->p
, dim
);
421 m1
= unimodular_complete(v
);
422 m2
= Matrix_Alloc(dim
, dim
+1);
423 for (i
= 0; i
< dim
-1 ; ++i
) {
424 Vector_Copy(m1
->p
[i
+1], m2
->p
[i
], dim
);
425 value_set_si(m2
->p
[i
][dim
], 0);
427 Vector_Set(m2
->p
[dim
-1], 0, dim
);
428 value_set_si(m2
->p
[dim
-1][dim
], 1);
429 q
= Polyhedron_Image(p
, m2
, p
->NbConstraints
+1+p
->NbRays
);
442 * Returns a full-dimensional polyhedron with the same number
443 * of integer points as P
444 * nvar specifies the number of variables
445 * The remaining dimensions are assumed to be parameters
447 * factor is NbEq x (nparam+2) matrix, containing stride constraints
448 * on the parameters; column nparam is the constant;
449 * column nparam+1 is the stride
451 * if factor is NULL, only remove equalities that don't affect
452 * the number of points
454 Polyhedron
*remove_equalities_p(Polyhedron
*P
, unsigned nvar
, Matrix
**factor
)
458 Polyhedron
*p
= P
, *q
;
459 unsigned dim
= p
->Dimension
;
465 f
= Matrix_Alloc(p
->NbEq
, dim
-nvar
+2);
470 while (nvar
> 0 && p
->NbEq
- skip
> 0) {
473 while (value_zero_p(p
->Constraint
[skip
][0]) &&
474 First_Non_Zero(p
->Constraint
[skip
]+1, nvar
) == -1)
479 Vector_Gcd(p
->Constraint
[skip
]+1, dim
+1, &g
);
480 Vector_AntiScale(p
->Constraint
[skip
]+1, p
->Constraint
[skip
]+1, g
, dim
+1);
481 Vector_Gcd(p
->Constraint
[skip
]+1, nvar
, &g
);
482 if (!factor
&& value_notone_p(g
) && value_notmone_p(g
)) {
487 Vector_Copy(p
->Constraint
[skip
]+1+nvar
, f
->p
[j
], dim
-nvar
+1);
488 value_assign(f
->p
[j
][dim
-nvar
+1], g
);
490 v
= Vector_Alloc(dim
);
491 Vector_AntiScale(p
->Constraint
[skip
]+1, v
->p
, g
, nvar
);
492 Vector_Set(v
->p
+nvar
, 0, dim
-nvar
);
493 m1
= unimodular_complete(v
);
494 m2
= Matrix_Alloc(dim
, dim
+1);
495 for (i
= 0; i
< dim
-1 ; ++i
) {
496 Vector_Copy(m1
->p
[i
+1], m2
->p
[i
], dim
);
497 value_set_si(m2
->p
[i
][dim
], 0);
499 Vector_Set(m2
->p
[dim
-1], 0, dim
);
500 value_set_si(m2
->p
[dim
-1][dim
], 1);
501 q
= Polyhedron_Image(p
, m2
, p
->NbConstraints
+1+p
->NbRays
);
520 static void free_singles(int **singles
, int dim
)
523 for (i
= 0; i
< dim
; ++i
)
528 static int **find_singles(Polyhedron
*P
, int dim
, int max
, int *nsingle
)
531 int **singles
= (int **) malloc(dim
* sizeof(int *));
534 for (i
= 0; i
< dim
; ++i
) {
535 singles
[i
] = (int *) malloc((max
+ 1) *sizeof(int));
539 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
540 for (j
= 0, prev
= -1; j
< dim
; ++j
) {
541 if (value_notzero_p(P
->Constraint
[i
][j
+1])) {
546 singles
[prev
][0] = -1;
552 if (prev
>= 0 && singles
[prev
][0] >= 0)
553 singles
[prev
][++singles
[prev
][0]] = i
;
556 for (j
= 0; j
< dim
; ++j
)
557 if (singles
[j
][0] > 0)
560 free_singles(singles
, dim
);
567 * The number of points in P is equal to factor time
568 * the number of points in the polyhedron returned.
569 * The return value is zero if no reduction can be found.
571 Polyhedron
* Polyhedron_Reduce(Polyhedron
*P
, Value
* factor
)
573 int i
, j
, nsingle
, k
, p
;
574 unsigned dim
= P
->Dimension
;
577 value_set_si(*factor
, 1);
579 assert (P
->NbEq
== 0);
581 singles
= find_singles(P
, dim
, 2, &nsingle
);
588 Matrix
*m
= Matrix_Alloc((dim
-nsingle
)+1, dim
+1);
594 for (i
= 0, j
= 0; i
< dim
; ++i
) {
595 if (singles
[i
][0] != 2)
596 value_set_si(m
->p
[j
++][i
], 1);
598 for (k
= 0; k
<= 1; ++k
) {
600 value_oppose(tmp
, P
->Constraint
[p
][dim
+1]);
601 if (value_pos_p(P
->Constraint
[p
][i
+1]))
602 mpz_cdiv_q(pos
, tmp
, P
->Constraint
[p
][i
+1]);
604 mpz_fdiv_q(neg
, tmp
, P
->Constraint
[p
][i
+1]);
606 value_substract(tmp
, neg
, pos
);
607 value_increment(tmp
, tmp
);
608 value_multiply(*factor
, *factor
, tmp
);
611 value_set_si(m
->p
[dim
-nsingle
][dim
], 1);
612 P
= Polyhedron_Image(P
, m
, P
->NbConstraints
);
614 free_singles(singles
, dim
);
625 * Replaces constraint a x >= c by x >= ceil(c/a)
626 * where "a" is a common factor in the coefficients
627 * old is the constraint; v points to an initialized
628 * value that this procedure can use.
629 * Return non-zero if something changed.
630 * Result is placed in new.
632 int ConstraintSimplify(Value
*old
, Value
*new, int len
, Value
* v
)
634 Vector_Gcd(old
+1, len
- 2, v
);
639 Vector_AntiScale(old
+1, new+1, *v
, len
-2);
640 mpz_fdiv_q(new[len
-1], old
[len
-1], *v
);
645 * Project on final dim dimensions
647 static Polyhedron
* Polyhedron_Project(Polyhedron
*P
, int dim
)
650 int remove
= P
->Dimension
- dim
;
654 if (P
->Dimension
== dim
)
655 return Polyhedron_Copy(P
);
657 T
= Matrix_Alloc(dim
+1, P
->Dimension
+1);
658 for (i
= 0; i
< dim
+1; ++i
)
659 value_set_si(T
->p
[i
][i
+remove
], 1);
660 I
= Polyhedron_Image(P
, T
, P
->NbConstraints
);
665 /* Constructs a new constraint that ensures that
666 * the first constraint is (strictly) smaller than
669 void smaller_constraint(Value
*a
, Value
*b
, Value
*c
, int pos
, int shift
,
670 int len
, int strict
, Value
*tmp
)
672 value_oppose(*tmp
, b
[pos
+1]);
673 value_set_si(c
[0], 1);
674 Vector_Combine(a
+1+shift
, b
+1+shift
, c
+1, *tmp
, a
[pos
+1], len
-shift
-1);
676 value_decrement(c
[len
-shift
-1], c
[len
-shift
-1]);
677 ConstraintSimplify(c
, c
, len
-shift
, tmp
);
680 struct section
{ Polyhedron
* D
; evalue E
; };
682 static Polyhedron
* ParamPolyhedron_Reduce_mod(Polyhedron
*P
, unsigned nvar
,
687 unsigned dim
= P
->Dimension
;
689 singles
= find_singles(P
, nvar
, P
->NbConstraints
, &nsingle
);
697 Matrix
*m
= Matrix_Alloc((dim
-nsingle
)+1, dim
+1);
701 evalue_set_si(&mone
, -1, 1);
702 C
= Polyhedron_Project(P
, dim
-nvar
);
706 for (i
= 0, j
= 0; i
< dim
; ++i
) {
707 if (i
>= nvar
|| singles
[i
][0] < 2)
708 value_set_si(m
->p
[j
++][i
], 1);
716 /* put those with positive coefficients first; number: p */
717 for (p
= 0, n
= singles
[i
][0]-1; p
<= n
; ) {
718 while (value_pos_p(P
->Constraint
[singles
[i
][1+p
]][i
+1]))
720 while (value_neg_p(P
->Constraint
[singles
[i
][1+n
]][i
+1]))
723 int t
= singles
[i
][1+p
];
724 singles
[i
][1+p
] = singles
[i
][1+n
];
731 assert (p
>= 1 && n
>= 1);
732 s
= (struct section
*) malloc(p
* n
* sizeof(struct section
));
733 M
= Matrix_Alloc((p
-1) + (n
-1), dim
-nvar
+2);
734 for (k
= 0; k
< p
; ++k
) {
735 for (k2
= 0; k2
< p
; ++k2
) {
740 P
->Constraint
[singles
[i
][1+k
]],
741 P
->Constraint
[singles
[i
][1+k2
]],
742 M
->p
[q
], i
, nvar
, dim
+2, k2
> k
, &g
);
744 for (l
= p
; l
< p
+n
; ++l
) {
745 for (l2
= p
; l2
< p
+n
; ++l2
) {
750 P
->Constraint
[singles
[i
][1+l2
]],
751 P
->Constraint
[singles
[i
][1+l
]],
752 M
->p
[q
], i
, nvar
, dim
+2, l2
> l
, &g
);
755 s
[nd
].D
= Constraints2Polyhedron(M2
, P
->NbRays
);
757 if (emptyQ(s
[nd
].D
)) {
758 Polyhedron_Free(s
[nd
].D
);
761 T
= DomainIntersection(s
[nd
].D
, C
, 0);
762 L
= bv_ceil3(P
->Constraint
[singles
[i
][1+k
]]+1+nvar
,
764 P
->Constraint
[singles
[i
][1+k
]][i
+1], T
);
765 U
= bv_ceil3(P
->Constraint
[singles
[i
][1+l
]]+1+nvar
,
767 P
->Constraint
[singles
[i
][1+l
]][i
+1], T
);
783 value_set_si(F
.d
, 0);
784 F
.x
.p
= new_enode(partition
, 2*nd
, dim
-nvar
);
785 for (k
= 0; k
< nd
; ++k
) {
786 EVALUE_SET_DOMAIN(F
.x
.p
->arr
[2*k
], s
[k
].D
);
787 value_clear(F
.x
.p
->arr
[2*k
+1].d
);
788 F
.x
.p
->arr
[2*k
+1] = s
[k
].E
;
793 free_evalue_refs(&F
);
796 value_set_si(m
->p
[dim
-nsingle
][dim
], 1);
797 P
= Polyhedron_Image(P
, m
, P
->NbConstraints
);
799 free_singles(singles
, nvar
);
803 free_evalue_refs(&mone
);
807 reduce_evalue(factor
);
813 Polyhedron
* ParamPolyhedron_Reduce(Polyhedron
*P
, unsigned nvar
,
816 return ParamPolyhedron_Reduce_mod(P
, nvar
, factor
);
819 Polyhedron
* ParamPolyhedron_Reduce(Polyhedron
*P
, unsigned nvar
,
825 evalue_set_si(&tmp
, 1, 1);
826 R
= ParamPolyhedron_Reduce_mod(P
, nvar
, &tmp
);
827 evalue_mod2table(&tmp
, P
->Dimension
- nvar
);
830 free_evalue_refs(&tmp
);
835 Bool
isIdentity(Matrix
*M
)
838 if (M
->NbRows
!= M
->NbColumns
)
841 for (i
= 0;i
< M
->NbRows
; i
++)
842 for (j
= 0; j
< M
->NbColumns
; j
++)
844 if(value_notone_p(M
->p
[i
][j
]))
847 if(value_notzero_p(M
->p
[i
][j
]))
853 void Param_Polyhedron_Print(FILE* DST
, Param_Polyhedron
*PP
, char **param_names
)
858 for(P
=PP
->D
;P
;P
=P
->next
) {
860 /* prints current val. dom. */
861 printf( "---------------------------------------\n" );
862 printf( "Domain :\n");
863 Print_Domain( stdout
, P
->Domain
, param_names
);
865 /* scan the vertices */
866 printf( "Vertices :\n");
867 FORALL_PVertex_in_ParamPolyhedron(V
,P
,PP
) {
869 /* prints each vertex */
870 Print_Vertex( stdout
, V
->Vertex
, param_names
);
873 END_FORALL_PVertex_in_ParamPolyhedron
;
877 void Enumeration_Print(FILE *Dst
, Enumeration
*en
, char **params
)
879 for (; en
; en
= en
->next
) {
880 Print_Domain(Dst
, en
->ValidityDomain
, params
);
881 print_evalue(Dst
, &en
->EP
, params
);
885 void Enumeration_mod2table(Enumeration
*en
, unsigned nparam
)
887 for (; en
; en
= en
->next
) {
888 evalue_mod2table(&en
->EP
, nparam
);
889 reduce_evalue(&en
->EP
);
893 size_t Enumeration_size(Enumeration
*en
)
897 for (; en
; en
= en
->next
) {
898 s
+= domain_size(en
->ValidityDomain
);
899 s
+= evalue_size(&en
->EP
);
904 void Free_ParamNames(char **params
, int m
)
911 int DomainIncludes(Polyhedron
*Pol1
, Polyhedron
*Pol2
)
914 for ( ; Pol1
; Pol1
= Pol1
->next
) {
915 for (P2
= Pol2
; P2
; P2
= P2
->next
)
916 if (!PolyhedronIncludes(Pol1
, P2
))
924 static Polyhedron
*p_simplify_constraints(Polyhedron
*P
, Vector
*row
,
925 Value
*g
, unsigned MaxRays
)
927 Polyhedron
*T
, *R
= P
;
928 int len
= P
->Dimension
+2;
931 /* Also look at equalities.
932 * If an equality can be "simplified" then there
933 * are no integer solutions anyway and the following loop
934 * will add a conflicting constraint
936 for (r
= 0; r
< R
->NbConstraints
; ++r
) {
937 if (ConstraintSimplify(R
->Constraint
[r
], row
->p
, len
, g
)) {
939 R
= AddConstraints(row
->p
, 1, R
, MaxRays
);
951 * Replaces constraint a x >= c by x >= ceil(c/a)
952 * where "a" is a common factor in the coefficients
953 * Destroys P and returns a newly allocated Polyhedron
954 * or just returns P in case no changes were made
956 Polyhedron
*DomainConstraintSimplify(Polyhedron
*P
, unsigned MaxRays
)
959 int len
= P
->Dimension
+2;
960 Vector
*row
= Vector_Alloc(len
);
962 Polyhedron
*R
= P
, *N
;
963 value_set_si(row
->p
[0], 1);
966 for (prev
= &R
; P
; P
= N
) {
969 T
= p_simplify_constraints(P
, row
, &g
, MaxRays
);
971 if (emptyQ(T
) && prev
!= &R
) {
983 if (R
->next
&& emptyQ(R
)) {
994 int line_minmax(Polyhedron
*I
, Value
*min
, Value
*max
)
999 value_oppose(I
->Constraint
[0][2], I
->Constraint
[0][2]);
1000 /* There should never be a remainder here */
1001 if (value_pos_p(I
->Constraint
[0][1]))
1002 mpz_fdiv_q(*min
, I
->Constraint
[0][2], I
->Constraint
[0][1]);
1004 mpz_fdiv_q(*min
, I
->Constraint
[0][2], I
->Constraint
[0][1]);
1005 value_assign(*max
, *min
);
1006 } else for (i
= 0; i
< I
->NbConstraints
; ++i
) {
1007 if (value_zero_p(I
->Constraint
[i
][1])) {
1012 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
1013 if (value_pos_p(I
->Constraint
[i
][1]))
1014 mpz_cdiv_q(*min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
1016 mpz_fdiv_q(*max
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
1024 PROCEDURES TO COMPUTE ENUMERATION. recursive procedure, recurse for
1027 @param pos index position of current loop index (1..hdim-1)
1028 @param P loop domain
1029 @param context context values for fixed indices
1030 @param exist number of existential variables
1031 @return the number of integer points in this
1035 void count_points_e (int pos
, Polyhedron
*P
, int exist
, int nparam
,
1036 Value
*context
, Value
*res
)
1041 value_set_si(*res
, 0);
1045 value_init(LB
); value_init(UB
); value_init(k
);
1049 if (lower_upper_bounds(pos
,P
,context
,&LB
,&UB
) !=0) {
1050 /* Problem if UB or LB is INFINITY */
1051 value_clear(LB
); value_clear(UB
); value_clear(k
);
1052 if (pos
> P
->Dimension
- nparam
- exist
)
1053 value_set_si(*res
, 1);
1055 value_set_si(*res
, -1);
1062 for (value_assign(k
,LB
); value_le(k
,UB
); value_increment(k
,k
)) {
1063 fprintf(stderr
, "(");
1064 for (i
=1; i
<pos
; i
++) {
1065 value_print(stderr
,P_VALUE_FMT
,context
[i
]);
1066 fprintf(stderr
,",");
1068 value_print(stderr
,P_VALUE_FMT
,k
);
1069 fprintf(stderr
,")\n");
1074 value_set_si(context
[pos
],0);
1075 if (value_lt(UB
,LB
)) {
1076 value_clear(LB
); value_clear(UB
); value_clear(k
);
1077 value_set_si(*res
, 0);
1082 value_set_si(*res
, 1);
1084 value_substract(k
,UB
,LB
);
1085 value_add_int(k
,k
,1);
1086 value_assign(*res
, k
);
1088 value_clear(LB
); value_clear(UB
); value_clear(k
);
1092 /*-----------------------------------------------------------------*/
1093 /* Optimization idea */
1094 /* If inner loops are not a function of k (the current index) */
1095 /* i.e. if P->Constraint[i][pos]==0 for all P following this and */
1097 /* Then CNT = (UB-LB+1)*count_points(pos+1, P->next, context) */
1098 /* (skip the for loop) */
1099 /*-----------------------------------------------------------------*/
1102 value_set_si(*res
, 0);
1103 for (value_assign(k
,LB
);value_le(k
,UB
);value_increment(k
,k
)) {
1104 /* Insert k in context */
1105 value_assign(context
[pos
],k
);
1106 count_points_e(pos
+1, P
->next
, exist
, nparam
, context
, &c
);
1107 if(value_notmone_p(c
))
1108 value_addto(*res
, *res
, c
);
1110 value_set_si(*res
, -1);
1113 if (pos
> P
->Dimension
- nparam
- exist
&&
1120 fprintf(stderr
,"%d\n",CNT
);
1124 value_set_si(context
[pos
],0);
1125 value_clear(LB
); value_clear(UB
); value_clear(k
);
1127 } /* count_points_e */
1129 int DomainContains(Polyhedron
*P
, Value
*list_args
, int len
,
1130 unsigned MaxRays
, int set
)
1135 if (P
->Dimension
== len
)
1136 return in_domain(P
, list_args
);
1138 assert(set
); // assume list_args is large enough
1139 assert((P
->Dimension
- len
) % 2 == 0);
1141 for (i
= 0; i
< P
->Dimension
- len
; i
+= 2) {
1143 for (j
= 0 ; j
< P
->NbEq
; ++j
)
1144 if (value_notzero_p(P
->Constraint
[j
][1+len
+i
]))
1146 assert(j
< P
->NbEq
);
1147 value_absolute(m
, P
->Constraint
[j
][1+len
+i
]);
1148 k
= First_Non_Zero(P
->Constraint
[j
]+1, len
);
1150 assert(First_Non_Zero(P
->Constraint
[j
]+1+k
+1, len
- k
- 1) == -1);
1151 mpz_fdiv_q(list_args
[len
+i
], list_args
[k
], m
);
1152 mpz_fdiv_r(list_args
[len
+i
+1], list_args
[k
], m
);
1156 return in_domain(P
, list_args
);