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 Polyhedron
*remove_equalities_p(Polyhedron
*P
, unsigned nvar
, Matrix
**factor
)
455 Polyhedron
*p
= P
, *q
;
456 unsigned dim
= p
->Dimension
;
461 f
= Matrix_Alloc(p
->NbEq
, dim
-nvar
+2);
465 while (nvar
> 0 && p
->NbEq
- skip
> 0) {
468 while (value_zero_p(p
->Constraint
[skip
][0]) &&
469 First_Non_Zero(p
->Constraint
[skip
]+1, nvar
) == -1)
474 Vector_Gcd(p
->Constraint
[skip
]+1, dim
+1, &g
);
475 Vector_AntiScale(p
->Constraint
[skip
]+1, p
->Constraint
[skip
]+1, g
, dim
+1);
476 Vector_Gcd(p
->Constraint
[skip
]+1, nvar
, &g
);
477 Vector_Copy(p
->Constraint
[skip
]+1+nvar
, f
->p
[j
], dim
-nvar
+1);
478 value_assign(f
->p
[j
][dim
-nvar
+1], g
);
479 v
= Vector_Alloc(dim
);
480 Vector_AntiScale(p
->Constraint
[skip
]+1, v
->p
, g
, nvar
);
481 Vector_Set(v
->p
+nvar
, 0, dim
-nvar
);
482 m1
= unimodular_complete(v
);
483 m2
= Matrix_Alloc(dim
, dim
+1);
484 for (i
= 0; i
< dim
-1 ; ++i
) {
485 Vector_Copy(m1
->p
[i
+1], m2
->p
[i
], dim
);
486 value_set_si(m2
->p
[i
][dim
], 0);
488 Vector_Set(m2
->p
[dim
-1], 0, dim
);
489 value_set_si(m2
->p
[dim
-1][dim
], 1);
490 q
= Polyhedron_Image(p
, m2
, p
->NbConstraints
+1+p
->NbRays
);
509 static void free_singles(int **singles
, int dim
)
512 for (i
= 0; i
< dim
; ++i
)
517 static int **find_singles(Polyhedron
*P
, int dim
, int max
, int *nsingle
)
520 int **singles
= (int **) malloc(dim
* sizeof(int *));
523 for (i
= 0; i
< dim
; ++i
) {
524 singles
[i
] = (int *) malloc((max
+ 1) *sizeof(int));
528 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
529 for (j
= 0, prev
= -1; j
< dim
; ++j
) {
530 if (value_notzero_p(P
->Constraint
[i
][j
+1])) {
535 singles
[prev
][0] = -1;
541 if (prev
>= 0 && singles
[prev
][0] >= 0)
542 singles
[prev
][++singles
[prev
][0]] = i
;
545 for (j
= 0; j
< dim
; ++j
)
546 if (singles
[j
][0] > 0)
549 free_singles(singles
, dim
);
556 * The number of points in P is equal to factor time
557 * the number of points in the polyhedron returned.
558 * The return value is zero if no reduction can be found.
560 Polyhedron
* Polyhedron_Reduce(Polyhedron
*P
, Value
* factor
)
562 int i
, j
, nsingle
, k
, p
;
563 unsigned dim
= P
->Dimension
;
566 value_set_si(*factor
, 1);
568 assert (P
->NbEq
== 0);
570 singles
= find_singles(P
, dim
, 2, &nsingle
);
577 Matrix
*m
= Matrix_Alloc((dim
-nsingle
)+1, dim
+1);
583 for (i
= 0, j
= 0; i
< dim
; ++i
) {
584 if (singles
[i
][0] != 2)
585 value_set_si(m
->p
[j
++][i
], 1);
587 for (k
= 0; k
<= 1; ++k
) {
589 value_oppose(tmp
, P
->Constraint
[p
][dim
+1]);
590 if (value_pos_p(P
->Constraint
[p
][i
+1]))
591 mpz_cdiv_q(pos
, tmp
, P
->Constraint
[p
][i
+1]);
593 mpz_fdiv_q(neg
, tmp
, P
->Constraint
[p
][i
+1]);
595 value_substract(tmp
, neg
, pos
);
596 value_increment(tmp
, tmp
);
597 value_multiply(*factor
, *factor
, tmp
);
600 value_set_si(m
->p
[dim
-nsingle
][dim
], 1);
601 P
= Polyhedron_Image(P
, m
, P
->NbConstraints
);
603 free_singles(singles
, dim
);
614 * Replaces constraint a x >= c by x >= ceil(c/a)
615 * where "a" is a common factor in the coefficients
616 * old is the constraint; v points to an initialized
617 * value that this procedure can use.
618 * Return non-zero if something changed.
619 * Result is placed in new.
621 int ConstraintSimplify(Value
*old
, Value
*new, int len
, Value
* v
)
623 Vector_Gcd(old
+1, len
- 2, v
);
628 Vector_AntiScale(old
+1, new+1, *v
, len
-2);
629 mpz_fdiv_q(new[len
-1], old
[len
-1], *v
);
634 * Project on final dim dimensions
636 static Polyhedron
* Polyhedron_Project(Polyhedron
*P
, int dim
)
639 int remove
= P
->Dimension
- dim
;
643 if (P
->Dimension
== dim
)
644 return Polyhedron_Copy(P
);
646 T
= Matrix_Alloc(dim
+1, P
->Dimension
+1);
647 for (i
= 0; i
< dim
+1; ++i
)
648 value_set_si(T
->p
[i
][i
+remove
], 1);
649 I
= Polyhedron_Image(P
, T
, P
->NbConstraints
);
654 /* Constructs a new constraint that ensures that
655 * the first constraint is (strictly) smaller than
658 void smaller_constraint(Value
*a
, Value
*b
, Value
*c
, int pos
, int shift
,
659 int len
, int strict
, Value
*tmp
)
661 value_oppose(*tmp
, b
[pos
+1]);
662 value_set_si(c
[0], 1);
663 Vector_Combine(a
+1+shift
, b
+1+shift
, c
+1, *tmp
, a
[pos
+1], len
-shift
-1);
665 value_decrement(c
[len
-shift
-1], c
[len
-shift
-1]);
666 ConstraintSimplify(c
, c
, len
-shift
, tmp
);
669 struct section
{ Polyhedron
* D
; evalue E
; };
671 static Polyhedron
* ParamPolyhedron_Reduce_mod(Polyhedron
*P
, unsigned nvar
,
676 unsigned dim
= P
->Dimension
;
678 singles
= find_singles(P
, nvar
, P
->NbConstraints
, &nsingle
);
686 Matrix
*m
= Matrix_Alloc((dim
-nsingle
)+1, dim
+1);
690 evalue_set_si(&mone
, -1, 1);
691 C
= Polyhedron_Project(P
, dim
-nvar
);
695 for (i
= 0, j
= 0; i
< dim
; ++i
) {
696 if (i
>= nvar
|| singles
[i
][0] < 2)
697 value_set_si(m
->p
[j
++][i
], 1);
705 /* put those with positive coefficients first; number: p */
706 for (p
= 0, n
= singles
[i
][0]-1; p
<= n
; ) {
707 while (value_pos_p(P
->Constraint
[singles
[i
][1+p
]][i
+1]))
709 while (value_neg_p(P
->Constraint
[singles
[i
][1+n
]][i
+1]))
712 int t
= singles
[i
][1+p
];
713 singles
[i
][1+p
] = singles
[i
][1+n
];
720 assert (p
>= 1 && n
>= 1);
721 s
= (struct section
*) malloc(p
* n
* sizeof(struct section
));
722 M
= Matrix_Alloc((p
-1) + (n
-1), dim
-nvar
+2);
723 for (k
= 0; k
< p
; ++k
) {
724 for (k2
= 0; k2
< p
; ++k2
) {
729 P
->Constraint
[singles
[i
][1+k
]],
730 P
->Constraint
[singles
[i
][1+k2
]],
731 M
->p
[q
], i
, nvar
, dim
+2, k2
> k
, &g
);
733 for (l
= p
; l
< p
+n
; ++l
) {
734 for (l2
= p
; l2
< p
+n
; ++l2
) {
739 P
->Constraint
[singles
[i
][1+l2
]],
740 P
->Constraint
[singles
[i
][1+l
]],
741 M
->p
[q
], i
, nvar
, dim
+2, l2
> l
, &g
);
744 s
[nd
].D
= Constraints2Polyhedron(M2
, P
->NbRays
);
746 if (emptyQ(s
[nd
].D
)) {
747 Polyhedron_Free(s
[nd
].D
);
750 T
= DomainIntersection(s
[nd
].D
, C
, 0);
751 L
= bv_ceil3(P
->Constraint
[singles
[i
][1+k
]]+1+nvar
,
753 P
->Constraint
[singles
[i
][1+k
]][i
+1], T
);
754 U
= bv_ceil3(P
->Constraint
[singles
[i
][1+l
]]+1+nvar
,
756 P
->Constraint
[singles
[i
][1+l
]][i
+1], T
);
772 value_set_si(F
.d
, 0);
773 F
.x
.p
= new_enode(partition
, 2*nd
, dim
-nvar
);
774 for (k
= 0; k
< nd
; ++k
) {
775 EVALUE_SET_DOMAIN(F
.x
.p
->arr
[2*k
], s
[k
].D
);
776 value_clear(F
.x
.p
->arr
[2*k
+1].d
);
777 F
.x
.p
->arr
[2*k
+1] = s
[k
].E
;
782 free_evalue_refs(&F
);
785 value_set_si(m
->p
[dim
-nsingle
][dim
], 1);
786 P
= Polyhedron_Image(P
, m
, P
->NbConstraints
);
788 free_singles(singles
, nvar
);
792 free_evalue_refs(&mone
);
796 reduce_evalue(factor
);
802 Polyhedron
* ParamPolyhedron_Reduce(Polyhedron
*P
, unsigned nvar
,
805 return ParamPolyhedron_Reduce_mod(P
, nvar
, factor
);
808 Polyhedron
* ParamPolyhedron_Reduce(Polyhedron
*P
, unsigned nvar
,
814 evalue_set_si(&tmp
, 1, 1);
815 R
= ParamPolyhedron_Reduce_mod(P
, nvar
, &tmp
);
816 evalue_mod2table(&tmp
, P
->Dimension
- nvar
);
819 free_evalue_refs(&tmp
);
824 Bool
isIdentity(Matrix
*M
)
827 if (M
->NbRows
!= M
->NbColumns
)
830 for (i
= 0;i
< M
->NbRows
; i
++)
831 for (j
= 0; j
< M
->NbColumns
; j
++)
833 if(value_notone_p(M
->p
[i
][j
]))
836 if(value_notzero_p(M
->p
[i
][j
]))
842 void Param_Polyhedron_Print(FILE* DST
, Param_Polyhedron
*PP
, char **param_names
)
847 for(P
=PP
->D
;P
;P
=P
->next
) {
849 /* prints current val. dom. */
850 printf( "---------------------------------------\n" );
851 printf( "Domain :\n");
852 Print_Domain( stdout
, P
->Domain
, param_names
);
854 /* scan the vertices */
855 printf( "Vertices :\n");
856 FORALL_PVertex_in_ParamPolyhedron(V
,P
,PP
) {
858 /* prints each vertex */
859 Print_Vertex( stdout
, V
->Vertex
, param_names
);
862 END_FORALL_PVertex_in_ParamPolyhedron
;
866 void Enumeration_Print(FILE *Dst
, Enumeration
*en
, char **params
)
868 for (; en
; en
= en
->next
) {
869 Print_Domain(Dst
, en
->ValidityDomain
, params
);
870 print_evalue(Dst
, &en
->EP
, params
);
874 void Enumeration_mod2table(Enumeration
*en
, unsigned nparam
)
876 for (; en
; en
= en
->next
) {
877 evalue_mod2table(&en
->EP
, nparam
);
878 reduce_evalue(&en
->EP
);
882 size_t Enumeration_size(Enumeration
*en
)
886 for (; en
; en
= en
->next
) {
887 s
+= domain_size(en
->ValidityDomain
);
888 s
+= evalue_size(&en
->EP
);
893 void Free_ParamNames(char **params
, int m
)
900 int DomainIncludes(Polyhedron
*Pol1
, Polyhedron
*Pol2
)
903 for ( ; Pol1
; Pol1
= Pol1
->next
) {
904 for (P2
= Pol2
; P2
; P2
= P2
->next
)
905 if (!PolyhedronIncludes(Pol1
, P2
))
913 static Polyhedron
*p_simplify_constraints(Polyhedron
*P
, Vector
*row
,
914 Value
*g
, unsigned MaxRays
)
916 Polyhedron
*T
, *R
= P
;
917 int len
= P
->Dimension
+2;
920 /* Also look at equalities.
921 * If an equality can be "simplified" then there
922 * are no integer solutions anyway and the following loop
923 * will add a conflicting constraint
925 for (r
= 0; r
< R
->NbConstraints
; ++r
) {
926 if (ConstraintSimplify(R
->Constraint
[r
], row
->p
, len
, g
)) {
928 R
= AddConstraints(row
->p
, 1, R
, MaxRays
);
940 * Replaces constraint a x >= c by x >= ceil(c/a)
941 * where "a" is a common factor in the coefficients
942 * Destroys P and returns a newly allocated Polyhedron
943 * or just returns P in case no changes were made
945 Polyhedron
*DomainConstraintSimplify(Polyhedron
*P
, unsigned MaxRays
)
948 int len
= P
->Dimension
+2;
949 Vector
*row
= Vector_Alloc(len
);
951 Polyhedron
*R
= P
, *N
;
952 value_set_si(row
->p
[0], 1);
955 for (prev
= &R
; P
; P
= N
) {
958 T
= p_simplify_constraints(P
, row
, &g
, MaxRays
);
960 if (emptyQ(T
) && prev
!= &R
) {
972 if (R
->next
&& emptyQ(R
)) {
983 int line_minmax(Polyhedron
*I
, Value
*min
, Value
*max
)
988 value_oppose(I
->Constraint
[0][2], I
->Constraint
[0][2]);
989 /* There should never be a remainder here */
990 if (value_pos_p(I
->Constraint
[0][1]))
991 mpz_fdiv_q(*min
, I
->Constraint
[0][2], I
->Constraint
[0][1]);
993 mpz_fdiv_q(*min
, I
->Constraint
[0][2], I
->Constraint
[0][1]);
994 value_assign(*max
, *min
);
995 } else for (i
= 0; i
< I
->NbConstraints
; ++i
) {
996 if (value_zero_p(I
->Constraint
[i
][1])) {
1001 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
1002 if (value_pos_p(I
->Constraint
[i
][1]))
1003 mpz_cdiv_q(*min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
1005 mpz_fdiv_q(*max
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
1013 PROCEDURES TO COMPUTE ENUMERATION. recursive procedure, recurse for
1016 @param pos index position of current loop index (1..hdim-1)
1017 @param P loop domain
1018 @param context context values for fixed indices
1019 @param exist number of existential variables
1020 @return the number of integer points in this
1024 void count_points_e (int pos
, Polyhedron
*P
, int exist
, int nparam
,
1025 Value
*context
, Value
*res
)
1030 value_set_si(*res
, 0);
1034 value_init(LB
); value_init(UB
); value_init(k
);
1038 if (lower_upper_bounds(pos
,P
,context
,&LB
,&UB
) !=0) {
1039 /* Problem if UB or LB is INFINITY */
1040 value_clear(LB
); value_clear(UB
); value_clear(k
);
1041 if (pos
> P
->Dimension
- nparam
- exist
)
1042 value_set_si(*res
, 1);
1044 value_set_si(*res
, -1);
1051 for (value_assign(k
,LB
); value_le(k
,UB
); value_increment(k
,k
)) {
1052 fprintf(stderr
, "(");
1053 for (i
=1; i
<pos
; i
++) {
1054 value_print(stderr
,P_VALUE_FMT
,context
[i
]);
1055 fprintf(stderr
,",");
1057 value_print(stderr
,P_VALUE_FMT
,k
);
1058 fprintf(stderr
,")\n");
1063 value_set_si(context
[pos
],0);
1064 if (value_lt(UB
,LB
)) {
1065 value_clear(LB
); value_clear(UB
); value_clear(k
);
1066 value_set_si(*res
, 0);
1071 value_set_si(*res
, 1);
1073 value_substract(k
,UB
,LB
);
1074 value_add_int(k
,k
,1);
1075 value_assign(*res
, k
);
1077 value_clear(LB
); value_clear(UB
); value_clear(k
);
1081 /*-----------------------------------------------------------------*/
1082 /* Optimization idea */
1083 /* If inner loops are not a function of k (the current index) */
1084 /* i.e. if P->Constraint[i][pos]==0 for all P following this and */
1086 /* Then CNT = (UB-LB+1)*count_points(pos+1, P->next, context) */
1087 /* (skip the for loop) */
1088 /*-----------------------------------------------------------------*/
1091 value_set_si(*res
, 0);
1092 for (value_assign(k
,LB
);value_le(k
,UB
);value_increment(k
,k
)) {
1093 /* Insert k in context */
1094 value_assign(context
[pos
],k
);
1095 count_points_e(pos
+1, P
->next
, exist
, nparam
, context
, &c
);
1096 if(value_notmone_p(c
))
1097 value_addto(*res
, *res
, c
);
1099 value_set_si(*res
, -1);
1102 if (pos
> P
->Dimension
- nparam
- exist
&&
1109 fprintf(stderr
,"%d\n",CNT
);
1113 value_set_si(context
[pos
],0);
1114 value_clear(LB
); value_clear(UB
); value_clear(k
);
1116 } /* count_points_e */
1118 int DomainContains(Polyhedron
*P
, Value
*list_args
, int len
,
1119 unsigned MaxRays
, int set
)
1124 if (P
->Dimension
== len
)
1125 return in_domain(P
, list_args
);
1127 assert(set
); // assume list_args is large enough
1128 assert((P
->Dimension
- len
) % 2 == 0);
1130 for (i
= 0; i
< P
->Dimension
- len
; i
+= 2) {
1132 for (j
= 0 ; j
< P
->NbEq
; ++j
)
1133 if (value_notzero_p(P
->Constraint
[j
][1+len
+i
]))
1135 assert(j
< P
->NbEq
);
1136 value_absolute(m
, P
->Constraint
[j
][1+len
+i
]);
1137 k
= First_Non_Zero(P
->Constraint
[j
]+1, len
);
1139 assert(First_Non_Zero(P
->Constraint
[j
]+1+k
+1, len
- k
- 1) == -1);
1140 mpz_fdiv_q(list_args
[len
+i
], list_args
[k
], m
);
1141 mpz_fdiv_r(list_args
[len
+i
+1], list_args
[k
], m
);
1145 return in_domain(P
, list_args
);