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 /* Delaunay triangulation */
213 for (i
= 0, r
= 1; i
< P
->NbRays
; ++i
) {
214 if (value_notzero_p(P
->Ray
[i
][dim
+1]))
216 Vector_Copy(P
->Ray
[i
], M
->p
[r
], dim
+1);
217 Inner_Product(M
->p
[r
]+1, M
->p
[r
]+1, dim
, &tmp
);
218 value_assign(M
->p
[r
][dim
+1], tmp
);
219 value_set_si(M
->p
[r
][dim
+2], 0);
224 L
= Rays2Polyhedron(M3
, NbMaxCons
);
227 M2
= Matrix_Alloc(dim
+1, dim
+2);
232 /* Usually R should still be 0 */
235 for (r
= 1; r
< P
->NbRays
; ++r
) {
236 value_set_si(M
->p
[r
][dim
+1], random_int((t
+1)*dim
)+1);
239 L
= Rays2Polyhedron(M3
, NbMaxCons
);
243 assert(t
<= MAX_TRY
);
248 for (i
= 0; i
< L
->NbConstraints
; ++i
) {
249 /* Ignore perpendicular facets, i.e., facets with 0 z-coordinate */
250 if (value_negz_p(L
->Constraint
[i
][dim
+1]))
252 if (value_notzero_p(L
->Constraint
[i
][dim
+2]))
254 for (j
= 1, r
= 1; j
< M
->NbRows
; ++j
) {
255 Inner_Product(M
->p
[j
]+1, L
->Constraint
[i
]+1, dim
+1, &tmp
);
256 if (value_notzero_p(tmp
))
260 Vector_Copy(M
->p
[j
]+1, M2
->p
[r
]+1, dim
);
261 value_set_si(M2
->p
[r
][0], 1);
262 value_set_si(M2
->p
[r
][dim
+1], 0);
266 Vector_Set(M2
->p
[0]+1, 0, dim
);
267 value_set_si(M2
->p
[0][0], 1);
268 value_set_si(M2
->p
[0][dim
+1], 1);
269 T
= Rays2Polyhedron(M2
, P
->NbConstraints
+1);
283 void check_triangulization(Polyhedron
*P
, Polyhedron
*T
)
285 Polyhedron
*C
, *D
, *E
, *F
, *G
, *U
;
286 for (C
= T
; C
; C
= C
->next
) {
290 U
= DomainConvex(DomainUnion(U
, C
, 100), 100);
291 for (D
= C
->next
; D
; D
= D
->next
) {
296 E
= DomainIntersection(C
, D
, 600);
297 assert(E
->NbRays
== 0 || E
->NbEq
>= 1);
303 assert(PolyhedronIncludes(U
, P
));
304 assert(PolyhedronIncludes(P
, U
));
307 void Euclid(Value a
, Value b
, Value
*x
, Value
*y
, Value
*g
)
309 Value c
, d
, e
, f
, tmp
;
316 value_absolute(c
, a
);
317 value_absolute(d
, b
);
320 while(value_pos_p(d
)) {
321 value_division(tmp
, c
, d
);
322 value_multiply(tmp
, tmp
, f
);
323 value_substract(e
, e
, tmp
);
324 value_division(tmp
, c
, d
);
325 value_multiply(tmp
, tmp
, d
);
326 value_substract(c
, c
, tmp
);
333 else if (value_pos_p(a
))
335 else value_oppose(*x
, e
);
339 value_multiply(tmp
, a
, *x
);
340 value_substract(tmp
, c
, tmp
);
341 value_division(*y
, tmp
, b
);
350 Matrix
* unimodular_complete(Vector
*row
)
352 Value g
, b
, c
, old
, tmp
;
361 m
= Matrix_Alloc(row
->Size
, row
->Size
);
362 for (j
= 0; j
< row
->Size
; ++j
) {
363 value_assign(m
->p
[0][j
], row
->p
[j
]);
365 value_assign(g
, row
->p
[0]);
366 for (i
= 1; value_zero_p(g
) && i
< row
->Size
; ++i
) {
367 for (j
= 0; j
< row
->Size
; ++j
) {
369 value_set_si(m
->p
[i
][j
], 1);
371 value_set_si(m
->p
[i
][j
], 0);
373 value_assign(g
, row
->p
[i
]);
375 for (; i
< row
->Size
; ++i
) {
376 value_assign(old
, g
);
377 Euclid(old
, row
->p
[i
], &c
, &b
, &g
);
379 for (j
= 0; j
< row
->Size
; ++j
) {
381 value_multiply(tmp
, row
->p
[j
], b
);
382 value_division(m
->p
[i
][j
], tmp
, old
);
384 value_assign(m
->p
[i
][j
], c
);
386 value_set_si(m
->p
[i
][j
], 0);
398 * Returns a full-dimensional polyhedron with the same number
399 * of integer points as P
401 Polyhedron
*remove_equalities(Polyhedron
*P
)
405 Polyhedron
*p
= Polyhedron_Copy(P
), *q
;
406 unsigned dim
= p
->Dimension
;
411 while (p
->NbEq
> 0) {
413 Vector_Gcd(p
->Constraint
[0]+1, dim
+1, &g
);
414 Vector_AntiScale(p
->Constraint
[0]+1, p
->Constraint
[0]+1, g
, dim
+1);
415 Vector_Gcd(p
->Constraint
[0]+1, dim
, &g
);
416 if (value_notone_p(g
) && value_notmone_p(g
)) {
418 p
= Empty_Polyhedron(0);
421 v
= Vector_Alloc(dim
);
422 Vector_Copy(p
->Constraint
[0]+1, v
->p
, dim
);
423 m1
= unimodular_complete(v
);
424 m2
= Matrix_Alloc(dim
, dim
+1);
425 for (i
= 0; i
< dim
-1 ; ++i
) {
426 Vector_Copy(m1
->p
[i
+1], m2
->p
[i
], dim
);
427 value_set_si(m2
->p
[i
][dim
], 0);
429 Vector_Set(m2
->p
[dim
-1], 0, dim
);
430 value_set_si(m2
->p
[dim
-1][dim
], 1);
431 q
= Polyhedron_Image(p
, m2
, p
->NbConstraints
+1+p
->NbRays
);
444 * Returns a full-dimensional polyhedron with the same number
445 * of integer points as P
446 * nvar specifies the number of variables
447 * The remaining dimensions are assumed to be parameters
449 * factor is NbEq x (nparam+2) matrix, containing stride constraints
450 * on the parameters; column nparam is the constant;
451 * column nparam+1 is the stride
453 * if factor is NULL, only remove equalities that don't affect
454 * the number of points
456 Polyhedron
*remove_equalities_p(Polyhedron
*P
, unsigned nvar
, Matrix
**factor
)
460 Polyhedron
*p
= P
, *q
;
461 unsigned dim
= p
->Dimension
;
467 f
= Matrix_Alloc(p
->NbEq
, dim
-nvar
+2);
472 while (nvar
> 0 && p
->NbEq
- skip
> 0) {
475 while (value_zero_p(p
->Constraint
[skip
][0]) &&
476 First_Non_Zero(p
->Constraint
[skip
]+1, nvar
) == -1)
481 Vector_Gcd(p
->Constraint
[skip
]+1, dim
+1, &g
);
482 Vector_AntiScale(p
->Constraint
[skip
]+1, p
->Constraint
[skip
]+1, g
, dim
+1);
483 Vector_Gcd(p
->Constraint
[skip
]+1, nvar
, &g
);
484 if (!factor
&& value_notone_p(g
) && value_notmone_p(g
)) {
489 Vector_Copy(p
->Constraint
[skip
]+1+nvar
, f
->p
[j
], dim
-nvar
+1);
490 value_assign(f
->p
[j
][dim
-nvar
+1], g
);
492 v
= Vector_Alloc(dim
);
493 Vector_AntiScale(p
->Constraint
[skip
]+1, v
->p
, g
, nvar
);
494 Vector_Set(v
->p
+nvar
, 0, dim
-nvar
);
495 m1
= unimodular_complete(v
);
496 m2
= Matrix_Alloc(dim
, dim
+1);
497 for (i
= 0; i
< dim
-1 ; ++i
) {
498 Vector_Copy(m1
->p
[i
+1], m2
->p
[i
], dim
);
499 value_set_si(m2
->p
[i
][dim
], 0);
501 Vector_Set(m2
->p
[dim
-1], 0, dim
);
502 value_set_si(m2
->p
[dim
-1][dim
], 1);
503 q
= Polyhedron_Image(p
, m2
, p
->NbConstraints
+1+p
->NbRays
);
522 static void free_singles(int **singles
, int dim
)
525 for (i
= 0; i
< dim
; ++i
)
530 static int **find_singles(Polyhedron
*P
, int dim
, int max
, int *nsingle
)
533 int **singles
= (int **) malloc(dim
* sizeof(int *));
536 for (i
= 0; i
< dim
; ++i
) {
537 singles
[i
] = (int *) malloc((max
+ 1) *sizeof(int));
541 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
542 for (j
= 0, prev
= -1; j
< dim
; ++j
) {
543 if (value_notzero_p(P
->Constraint
[i
][j
+1])) {
548 singles
[prev
][0] = -1;
554 if (prev
>= 0 && singles
[prev
][0] >= 0)
555 singles
[prev
][++singles
[prev
][0]] = i
;
558 for (j
= 0; j
< dim
; ++j
)
559 if (singles
[j
][0] > 0)
562 free_singles(singles
, dim
);
569 * The number of points in P is equal to factor time
570 * the number of points in the polyhedron returned.
571 * The return value is zero if no reduction can be found.
573 Polyhedron
* Polyhedron_Reduce(Polyhedron
*P
, Value
* factor
)
575 int i
, j
, nsingle
, k
, p
;
576 unsigned dim
= P
->Dimension
;
579 value_set_si(*factor
, 1);
581 assert (P
->NbEq
== 0);
583 singles
= find_singles(P
, dim
, 2, &nsingle
);
590 Matrix
*m
= Matrix_Alloc((dim
-nsingle
)+1, dim
+1);
596 for (i
= 0, j
= 0; i
< dim
; ++i
) {
597 if (singles
[i
][0] != 2)
598 value_set_si(m
->p
[j
++][i
], 1);
600 for (k
= 0; k
<= 1; ++k
) {
602 value_oppose(tmp
, P
->Constraint
[p
][dim
+1]);
603 if (value_pos_p(P
->Constraint
[p
][i
+1]))
604 mpz_cdiv_q(pos
, tmp
, P
->Constraint
[p
][i
+1]);
606 mpz_fdiv_q(neg
, tmp
, P
->Constraint
[p
][i
+1]);
608 value_substract(tmp
, neg
, pos
);
609 value_increment(tmp
, tmp
);
610 value_multiply(*factor
, *factor
, tmp
);
613 value_set_si(m
->p
[dim
-nsingle
][dim
], 1);
614 P
= Polyhedron_Image(P
, m
, P
->NbConstraints
);
616 free_singles(singles
, dim
);
627 * Replaces constraint a x >= c by x >= ceil(c/a)
628 * where "a" is a common factor in the coefficients
629 * old is the constraint; v points to an initialized
630 * value that this procedure can use.
631 * Return non-zero if something changed.
632 * Result is placed in new.
634 int ConstraintSimplify(Value
*old
, Value
*new, int len
, Value
* v
)
636 Vector_Gcd(old
+1, len
- 2, v
);
641 Vector_AntiScale(old
+1, new+1, *v
, len
-2);
642 mpz_fdiv_q(new[len
-1], old
[len
-1], *v
);
647 * Project on final dim dimensions
649 Polyhedron
* Polyhedron_Project(Polyhedron
*P
, int dim
)
652 int remove
= P
->Dimension
- dim
;
656 if (P
->Dimension
== dim
)
657 return Polyhedron_Copy(P
);
659 T
= Matrix_Alloc(dim
+1, P
->Dimension
+1);
660 for (i
= 0; i
< dim
+1; ++i
)
661 value_set_si(T
->p
[i
][i
+remove
], 1);
662 I
= Polyhedron_Image(P
, T
, P
->NbConstraints
);
667 /* Constructs a new constraint that ensures that
668 * the first constraint is (strictly) smaller than
671 void smaller_constraint(Value
*a
, Value
*b
, Value
*c
, int pos
, int shift
,
672 int len
, int strict
, Value
*tmp
)
674 value_oppose(*tmp
, b
[pos
+1]);
675 value_set_si(c
[0], 1);
676 Vector_Combine(a
+1+shift
, b
+1+shift
, c
+1, *tmp
, a
[pos
+1], len
-shift
-1);
678 value_decrement(c
[len
-shift
-1], c
[len
-shift
-1]);
679 ConstraintSimplify(c
, c
, len
-shift
, tmp
);
682 struct section
{ Polyhedron
* D
; evalue E
; };
684 static Polyhedron
* ParamPolyhedron_Reduce_mod(Polyhedron
*P
, unsigned nvar
,
689 unsigned dim
= P
->Dimension
;
691 singles
= find_singles(P
, nvar
, P
->NbConstraints
, &nsingle
);
699 Matrix
*m
= Matrix_Alloc((dim
-nsingle
)+1, dim
+1);
703 evalue_set_si(&mone
, -1, 1);
704 C
= Polyhedron_Project(P
, dim
-nvar
);
708 for (i
= 0, j
= 0; i
< dim
; ++i
) {
709 if (i
>= nvar
|| singles
[i
][0] < 2)
710 value_set_si(m
->p
[j
++][i
], 1);
718 /* put those with positive coefficients first; number: p */
719 for (p
= 0, n
= singles
[i
][0]-1; p
<= n
; ) {
720 while (value_pos_p(P
->Constraint
[singles
[i
][1+p
]][i
+1]))
722 while (value_neg_p(P
->Constraint
[singles
[i
][1+n
]][i
+1]))
725 int t
= singles
[i
][1+p
];
726 singles
[i
][1+p
] = singles
[i
][1+n
];
733 assert (p
>= 1 && n
>= 1);
734 s
= (struct section
*) malloc(p
* n
* sizeof(struct section
));
735 M
= Matrix_Alloc((p
-1) + (n
-1), dim
-nvar
+2);
736 for (k
= 0; k
< p
; ++k
) {
737 for (k2
= 0; k2
< p
; ++k2
) {
742 P
->Constraint
[singles
[i
][1+k
]],
743 P
->Constraint
[singles
[i
][1+k2
]],
744 M
->p
[q
], i
, nvar
, dim
+2, k2
> k
, &g
);
746 for (l
= p
; l
< p
+n
; ++l
) {
747 for (l2
= p
; l2
< p
+n
; ++l2
) {
752 P
->Constraint
[singles
[i
][1+l2
]],
753 P
->Constraint
[singles
[i
][1+l
]],
754 M
->p
[q
], i
, nvar
, dim
+2, l2
> l
, &g
);
757 s
[nd
].D
= Constraints2Polyhedron(M2
, P
->NbRays
);
759 if (emptyQ(s
[nd
].D
)) {
760 Polyhedron_Free(s
[nd
].D
);
763 T
= DomainIntersection(s
[nd
].D
, C
, 0);
764 L
= bv_ceil3(P
->Constraint
[singles
[i
][1+k
]]+1+nvar
,
766 P
->Constraint
[singles
[i
][1+k
]][i
+1], T
);
767 U
= bv_ceil3(P
->Constraint
[singles
[i
][1+l
]]+1+nvar
,
769 P
->Constraint
[singles
[i
][1+l
]][i
+1], T
);
785 value_set_si(F
.d
, 0);
786 F
.x
.p
= new_enode(partition
, 2*nd
, dim
-nvar
);
787 for (k
= 0; k
< nd
; ++k
) {
788 EVALUE_SET_DOMAIN(F
.x
.p
->arr
[2*k
], s
[k
].D
);
789 value_clear(F
.x
.p
->arr
[2*k
+1].d
);
790 F
.x
.p
->arr
[2*k
+1] = s
[k
].E
;
795 free_evalue_refs(&F
);
798 value_set_si(m
->p
[dim
-nsingle
][dim
], 1);
799 P
= Polyhedron_Image(P
, m
, P
->NbConstraints
);
801 free_singles(singles
, nvar
);
805 free_evalue_refs(&mone
);
809 reduce_evalue(factor
);
815 Polyhedron
* ParamPolyhedron_Reduce(Polyhedron
*P
, unsigned nvar
,
818 return ParamPolyhedron_Reduce_mod(P
, nvar
, factor
);
821 Polyhedron
* ParamPolyhedron_Reduce(Polyhedron
*P
, unsigned nvar
,
827 evalue_set_si(&tmp
, 1, 1);
828 R
= ParamPolyhedron_Reduce_mod(P
, nvar
, &tmp
);
829 evalue_mod2table(&tmp
, P
->Dimension
- nvar
);
832 free_evalue_refs(&tmp
);
837 Bool
isIdentity(Matrix
*M
)
840 if (M
->NbRows
!= M
->NbColumns
)
843 for (i
= 0;i
< M
->NbRows
; i
++)
844 for (j
= 0; j
< M
->NbColumns
; j
++)
846 if(value_notone_p(M
->p
[i
][j
]))
849 if(value_notzero_p(M
->p
[i
][j
]))
855 void Param_Polyhedron_Print(FILE* DST
, Param_Polyhedron
*PP
, char **param_names
)
860 for(P
=PP
->D
;P
;P
=P
->next
) {
862 /* prints current val. dom. */
863 printf( "---------------------------------------\n" );
864 printf( "Domain :\n");
865 Print_Domain( stdout
, P
->Domain
, param_names
);
867 /* scan the vertices */
868 printf( "Vertices :\n");
869 FORALL_PVertex_in_ParamPolyhedron(V
,P
,PP
) {
871 /* prints each vertex */
872 Print_Vertex( stdout
, V
->Vertex
, param_names
);
875 END_FORALL_PVertex_in_ParamPolyhedron
;
879 void Enumeration_Print(FILE *Dst
, Enumeration
*en
, char **params
)
881 for (; en
; en
= en
->next
) {
882 Print_Domain(Dst
, en
->ValidityDomain
, params
);
883 print_evalue(Dst
, &en
->EP
, params
);
887 void Enumeration_mod2table(Enumeration
*en
, unsigned nparam
)
889 for (; en
; en
= en
->next
) {
890 evalue_mod2table(&en
->EP
, nparam
);
891 reduce_evalue(&en
->EP
);
895 size_t Enumeration_size(Enumeration
*en
)
899 for (; en
; en
= en
->next
) {
900 s
+= domain_size(en
->ValidityDomain
);
901 s
+= evalue_size(&en
->EP
);
906 void Free_ParamNames(char **params
, int m
)
913 int DomainIncludes(Polyhedron
*Pol1
, Polyhedron
*Pol2
)
916 for ( ; Pol1
; Pol1
= Pol1
->next
) {
917 for (P2
= Pol2
; P2
; P2
= P2
->next
)
918 if (!PolyhedronIncludes(Pol1
, P2
))
926 static Polyhedron
*p_simplify_constraints(Polyhedron
*P
, Vector
*row
,
927 Value
*g
, unsigned MaxRays
)
929 Polyhedron
*T
, *R
= P
;
930 int len
= P
->Dimension
+2;
933 /* Also look at equalities.
934 * If an equality can be "simplified" then there
935 * are no integer solutions anyway and the following loop
936 * will add a conflicting constraint
938 for (r
= 0; r
< R
->NbConstraints
; ++r
) {
939 if (ConstraintSimplify(R
->Constraint
[r
], row
->p
, len
, g
)) {
941 R
= AddConstraints(row
->p
, 1, R
, MaxRays
);
953 * Replaces constraint a x >= c by x >= ceil(c/a)
954 * where "a" is a common factor in the coefficients
955 * Destroys P and returns a newly allocated Polyhedron
956 * or just returns P in case no changes were made
958 Polyhedron
*DomainConstraintSimplify(Polyhedron
*P
, unsigned MaxRays
)
961 int len
= P
->Dimension
+2;
962 Vector
*row
= Vector_Alloc(len
);
964 Polyhedron
*R
= P
, *N
;
965 value_set_si(row
->p
[0], 1);
968 for (prev
= &R
; P
; P
= N
) {
971 T
= p_simplify_constraints(P
, row
, &g
, MaxRays
);
973 if (emptyQ(T
) && prev
!= &R
) {
985 if (R
->next
&& emptyQ(R
)) {
996 int line_minmax(Polyhedron
*I
, Value
*min
, Value
*max
)
1001 value_oppose(I
->Constraint
[0][2], I
->Constraint
[0][2]);
1002 /* There should never be a remainder here */
1003 if (value_pos_p(I
->Constraint
[0][1]))
1004 mpz_fdiv_q(*min
, I
->Constraint
[0][2], I
->Constraint
[0][1]);
1006 mpz_fdiv_q(*min
, I
->Constraint
[0][2], I
->Constraint
[0][1]);
1007 value_assign(*max
, *min
);
1008 } else for (i
= 0; i
< I
->NbConstraints
; ++i
) {
1009 if (value_zero_p(I
->Constraint
[i
][1])) {
1014 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
1015 if (value_pos_p(I
->Constraint
[i
][1]))
1016 mpz_cdiv_q(*min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
1018 mpz_fdiv_q(*max
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
1026 PROCEDURES TO COMPUTE ENUMERATION. recursive procedure, recurse for
1029 @param pos index position of current loop index (1..hdim-1)
1030 @param P loop domain
1031 @param context context values for fixed indices
1032 @param exist number of existential variables
1033 @return the number of integer points in this
1037 void count_points_e (int pos
, Polyhedron
*P
, int exist
, int nparam
,
1038 Value
*context
, Value
*res
)
1043 value_set_si(*res
, 0);
1047 value_init(LB
); value_init(UB
); value_init(k
);
1051 if (lower_upper_bounds(pos
,P
,context
,&LB
,&UB
) !=0) {
1052 /* Problem if UB or LB is INFINITY */
1053 value_clear(LB
); value_clear(UB
); value_clear(k
);
1054 if (pos
> P
->Dimension
- nparam
- exist
)
1055 value_set_si(*res
, 1);
1057 value_set_si(*res
, -1);
1064 for (value_assign(k
,LB
); value_le(k
,UB
); value_increment(k
,k
)) {
1065 fprintf(stderr
, "(");
1066 for (i
=1; i
<pos
; i
++) {
1067 value_print(stderr
,P_VALUE_FMT
,context
[i
]);
1068 fprintf(stderr
,",");
1070 value_print(stderr
,P_VALUE_FMT
,k
);
1071 fprintf(stderr
,")\n");
1076 value_set_si(context
[pos
],0);
1077 if (value_lt(UB
,LB
)) {
1078 value_clear(LB
); value_clear(UB
); value_clear(k
);
1079 value_set_si(*res
, 0);
1084 value_set_si(*res
, 1);
1086 value_substract(k
,UB
,LB
);
1087 value_add_int(k
,k
,1);
1088 value_assign(*res
, k
);
1090 value_clear(LB
); value_clear(UB
); value_clear(k
);
1094 /*-----------------------------------------------------------------*/
1095 /* Optimization idea */
1096 /* If inner loops are not a function of k (the current index) */
1097 /* i.e. if P->Constraint[i][pos]==0 for all P following this and */
1099 /* Then CNT = (UB-LB+1)*count_points(pos+1, P->next, context) */
1100 /* (skip the for loop) */
1101 /*-----------------------------------------------------------------*/
1104 value_set_si(*res
, 0);
1105 for (value_assign(k
,LB
);value_le(k
,UB
);value_increment(k
,k
)) {
1106 /* Insert k in context */
1107 value_assign(context
[pos
],k
);
1108 count_points_e(pos
+1, P
->next
, exist
, nparam
, context
, &c
);
1109 if(value_notmone_p(c
))
1110 value_addto(*res
, *res
, c
);
1112 value_set_si(*res
, -1);
1115 if (pos
> P
->Dimension
- nparam
- exist
&&
1122 fprintf(stderr
,"%d\n",CNT
);
1126 value_set_si(context
[pos
],0);
1127 value_clear(LB
); value_clear(UB
); value_clear(k
);
1129 } /* count_points_e */
1131 int DomainContains(Polyhedron
*P
, Value
*list_args
, int len
,
1132 unsigned MaxRays
, int set
)
1137 if (P
->Dimension
== len
)
1138 return in_domain(P
, list_args
);
1140 assert(set
); // assume list_args is large enough
1141 assert((P
->Dimension
- len
) % 2 == 0);
1143 for (i
= 0; i
< P
->Dimension
- len
; i
+= 2) {
1145 for (j
= 0 ; j
< P
->NbEq
; ++j
)
1146 if (value_notzero_p(P
->Constraint
[j
][1+len
+i
]))
1148 assert(j
< P
->NbEq
);
1149 value_absolute(m
, P
->Constraint
[j
][1+len
+i
]);
1150 k
= First_Non_Zero(P
->Constraint
[j
]+1, len
);
1152 assert(First_Non_Zero(P
->Constraint
[j
]+1+k
+1, len
- k
- 1) == -1);
1153 mpz_fdiv_q(list_args
[len
+i
], list_args
[k
], m
);
1154 mpz_fdiv_r(list_args
[len
+i
+1], list_args
[k
], m
);
1158 return in_domain(P
, list_args
);