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);
222 L
= Rays2Polyhedron(M
, NbMaxCons
);
224 M2
= Matrix_Alloc(dim
+1, dim
+2);
225 Vector_Set(M2
->p
[0]+1, 0, dim
);
226 value_set_si(M2
->p
[0][0], 1);
227 value_set_si(M2
->p
[0][dim
+1], 1);
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);
238 L
= Rays2Polyhedron(M
, NbMaxCons
);
241 assert(t
<= MAX_TRY
);
246 for (i
= 0; i
< L
->NbConstraints
; ++i
) {
247 if (value_negz_p(L
->Constraint
[i
][dim
+1]))
249 if (value_notzero_p(L
->Constraint
[i
][dim
+2]))
251 for (j
= 1, r
= 1; j
< M
->NbRows
; ++j
) {
252 Inner_Product(M
->p
[j
]+1, L
->Constraint
[i
]+1, dim
+1, &tmp
);
253 if (value_notzero_p(tmp
))
257 Vector_Copy(M
->p
[j
]+1, M2
->p
[r
]+1, dim
);
258 value_set_si(M2
->p
[r
][0], 1);
259 value_set_si(M2
->p
[r
][dim
+1], 0);
263 T
= Rays2Polyhedron(M2
, P
->NbConstraints
+1);
277 void check_triangulization(Polyhedron
*P
, Polyhedron
*T
)
279 Polyhedron
*C
, *D
, *E
, *F
, *G
, *U
;
280 for (C
= T
; C
; C
= C
->next
) {
284 U
= DomainConvex(DomainUnion(U
, C
, 100), 100);
285 for (D
= C
->next
; D
; D
= D
->next
) {
290 E
= DomainIntersection(C
, D
, 600);
291 assert(E
->NbRays
== 0 || E
->NbEq
>= 1);
297 assert(PolyhedronIncludes(U
, P
));
298 assert(PolyhedronIncludes(P
, U
));
301 void Euclid(Value a
, Value b
, Value
*x
, Value
*y
, Value
*g
)
303 Value c
, d
, e
, f
, tmp
;
310 value_absolute(c
, a
);
311 value_absolute(d
, b
);
314 while(value_pos_p(d
)) {
315 value_division(tmp
, c
, d
);
316 value_multiply(tmp
, tmp
, f
);
317 value_substract(e
, e
, tmp
);
318 value_division(tmp
, c
, d
);
319 value_multiply(tmp
, tmp
, d
);
320 value_substract(c
, c
, tmp
);
327 else if (value_pos_p(a
))
329 else value_oppose(*x
, e
);
333 value_multiply(tmp
, a
, *x
);
334 value_substract(tmp
, c
, tmp
);
335 value_division(*y
, tmp
, b
);
344 Matrix
* unimodular_complete(Vector
*row
)
346 Value g
, b
, c
, old
, tmp
;
355 m
= Matrix_Alloc(row
->Size
, row
->Size
);
356 for (j
= 0; j
< row
->Size
; ++j
) {
357 value_assign(m
->p
[0][j
], row
->p
[j
]);
359 value_assign(g
, row
->p
[0]);
360 for (i
= 1; value_zero_p(g
) && i
< row
->Size
; ++i
) {
361 for (j
= 0; j
< row
->Size
; ++j
) {
363 value_set_si(m
->p
[i
][j
], 1);
365 value_set_si(m
->p
[i
][j
], 0);
367 value_assign(g
, row
->p
[i
]);
369 for (; i
< row
->Size
; ++i
) {
370 value_assign(old
, g
);
371 Euclid(old
, row
->p
[i
], &c
, &b
, &g
);
373 for (j
= 0; j
< row
->Size
; ++j
) {
375 value_multiply(tmp
, row
->p
[j
], b
);
376 value_division(m
->p
[i
][j
], tmp
, old
);
378 value_assign(m
->p
[i
][j
], c
);
380 value_set_si(m
->p
[i
][j
], 0);
392 * Returns a full-dimensional polyhedron with the same number
393 * of integer points as P
395 Polyhedron
*remove_equalities(Polyhedron
*P
)
399 Polyhedron
*p
= Polyhedron_Copy(P
), *q
;
400 unsigned dim
= p
->Dimension
;
405 while (p
->NbEq
> 0) {
407 Vector_Gcd(p
->Constraint
[0]+1, dim
+1, &g
);
408 Vector_AntiScale(p
->Constraint
[0]+1, p
->Constraint
[0]+1, g
, dim
+1);
409 Vector_Gcd(p
->Constraint
[0]+1, dim
, &g
);
410 if (value_notone_p(g
) && value_notmone_p(g
)) {
412 p
= Empty_Polyhedron(0);
415 v
= Vector_Alloc(dim
);
416 Vector_Copy(p
->Constraint
[0]+1, v
->p
, dim
);
417 m1
= unimodular_complete(v
);
418 m2
= Matrix_Alloc(dim
, dim
+1);
419 for (i
= 0; i
< dim
-1 ; ++i
) {
420 Vector_Copy(m1
->p
[i
+1], m2
->p
[i
], dim
);
421 value_set_si(m2
->p
[i
][dim
], 0);
423 Vector_Set(m2
->p
[dim
-1], 0, dim
);
424 value_set_si(m2
->p
[dim
-1][dim
], 1);
425 q
= Polyhedron_Image(p
, m2
, p
->NbConstraints
+1+p
->NbRays
);
438 * Returns a full-dimensional polyhedron with the same number
439 * of integer points as P
440 * nvar specifies the number of variables
441 * The remaining dimensions are assumed to be parameters
443 * factor is NbEq x (nparam+2) matrix, containing stride constraints
444 * on the parameters; column nparam is the constant;
445 * column nparam+1 is the stride
447 Polyhedron
*remove_equalities_p(Polyhedron
*P
, unsigned nvar
, Matrix
**factor
)
451 Polyhedron
*p
= P
, *q
;
452 unsigned dim
= p
->Dimension
;
457 f
= Matrix_Alloc(p
->NbEq
, dim
-nvar
+2);
461 while (nvar
> 0 && p
->NbEq
- skip
> 0) {
464 while (value_zero_p(p
->Constraint
[skip
][0]) &&
465 First_Non_Zero(p
->Constraint
[skip
]+1, nvar
) == -1)
470 Vector_Gcd(p
->Constraint
[skip
]+1, dim
+1, &g
);
471 Vector_AntiScale(p
->Constraint
[skip
]+1, p
->Constraint
[skip
]+1, g
, dim
+1);
472 Vector_Gcd(p
->Constraint
[skip
]+1, nvar
, &g
);
473 Vector_Copy(p
->Constraint
[skip
]+1+nvar
, f
->p
[j
], dim
-nvar
+1);
474 value_assign(f
->p
[j
][dim
-nvar
+1], g
);
475 v
= Vector_Alloc(dim
);
476 Vector_AntiScale(p
->Constraint
[skip
]+1, v
->p
, g
, nvar
);
477 Vector_Set(v
->p
+nvar
, 0, dim
-nvar
);
478 m1
= unimodular_complete(v
);
479 m2
= Matrix_Alloc(dim
, dim
+1);
480 for (i
= 0; i
< dim
-1 ; ++i
) {
481 Vector_Copy(m1
->p
[i
+1], m2
->p
[i
], dim
);
482 value_set_si(m2
->p
[i
][dim
], 0);
484 Vector_Set(m2
->p
[dim
-1], 0, dim
);
485 value_set_si(m2
->p
[dim
-1][dim
], 1);
486 q
= Polyhedron_Image(p
, m2
, p
->NbConstraints
+1+p
->NbRays
);
505 static void free_singles(int **singles
, int dim
)
508 for (i
= 0; i
< dim
; ++i
)
513 static int **find_singles(Polyhedron
*P
, int dim
, int max
, int *nsingle
)
516 int **singles
= (int **) malloc(dim
* sizeof(int *));
519 for (i
= 0; i
< dim
; ++i
) {
520 singles
[i
] = (int *) malloc((max
+ 1) *sizeof(int));
524 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
525 for (j
= 0, prev
= -1; j
< dim
; ++j
) {
526 if (value_notzero_p(P
->Constraint
[i
][j
+1])) {
531 singles
[prev
][0] = -1;
537 if (prev
>= 0 && singles
[prev
][0] >= 0)
538 singles
[prev
][++singles
[prev
][0]] = i
;
541 for (j
= 0; j
< dim
; ++j
)
542 if (singles
[j
][0] > 0)
545 free_singles(singles
, dim
);
552 * The number of points in P is equal to factor time
553 * the number of points in the polyhedron returned.
554 * The return value is zero if no reduction can be found.
556 Polyhedron
* Polyhedron_Reduce(Polyhedron
*P
, Value
* factor
)
558 int i
, j
, nsingle
, k
, p
;
559 unsigned dim
= P
->Dimension
;
562 value_set_si(*factor
, 1);
564 assert (P
->NbEq
== 0);
566 singles
= find_singles(P
, dim
, 2, &nsingle
);
573 Matrix
*m
= Matrix_Alloc((dim
-nsingle
)+1, dim
+1);
579 for (i
= 0, j
= 0; i
< dim
; ++i
) {
580 if (singles
[i
][0] != 2)
581 value_set_si(m
->p
[j
++][i
], 1);
583 for (k
= 0; k
<= 1; ++k
) {
585 value_oppose(tmp
, P
->Constraint
[p
][dim
+1]);
586 if (value_pos_p(P
->Constraint
[p
][i
+1]))
587 mpz_cdiv_q(pos
, tmp
, P
->Constraint
[p
][i
+1]);
589 mpz_fdiv_q(neg
, tmp
, P
->Constraint
[p
][i
+1]);
591 value_substract(tmp
, neg
, pos
);
592 value_increment(tmp
, tmp
);
593 value_multiply(*factor
, *factor
, tmp
);
596 value_set_si(m
->p
[dim
-nsingle
][dim
], 1);
597 P
= Polyhedron_Image(P
, m
, P
->NbConstraints
);
599 free_singles(singles
, dim
);
609 struct section
{ Polyhedron
* D
; evalue E
; };
611 static Polyhedron
* ParamPolyhedron_Reduce_mod(Polyhedron
*P
, unsigned nvar
,
616 unsigned dim
= P
->Dimension
;
618 singles
= find_singles(P
, nvar
, P
->NbConstraints
, &nsingle
);
625 Matrix
*m
= Matrix_Alloc((dim
-nsingle
)+1, dim
+1);
629 evalue_set_si(&mone
, -1, 1);
633 for (i
= 0, j
= 0; i
< dim
; ++i
) {
634 if (i
>= nvar
|| singles
[i
][0] < 2)
635 value_set_si(m
->p
[j
++][i
], 1);
643 /* put those with positive coefficients first; number: p */
644 for (p
= 0, n
= singles
[i
][0]-1; p
<= n
; ) {
645 while (value_pos_p(P
->Constraint
[singles
[i
][1+p
]][i
+1]))
647 while (value_neg_p(P
->Constraint
[singles
[i
][1+n
]][i
+1]))
650 int t
= singles
[i
][1+p
];
651 singles
[i
][1+p
] = singles
[i
][1+n
];
658 assert (p
>= 1 && n
>= 1);
659 s
= (struct section
*) malloc(p
* n
* sizeof(struct section
));
660 M
= Matrix_Alloc((p
-1) + (n
-1), dim
-nvar
+2);
661 for (k
= 0; k
< p
; ++k
) {
662 for (k2
= 0; k2
< p
; ++k2
) {
666 value_oppose(tmp
, P
->Constraint
[singles
[i
][1+k2
]][i
+1]);
667 value_set_si(M
->p
[q
][0], 1);
668 Vector_Combine(P
->Constraint
[singles
[i
][1+k
]]+1+nvar
,
669 P
->Constraint
[singles
[i
][1+k2
]]+1+nvar
,
672 P
->Constraint
[singles
[i
][1+k
]][i
+1],
675 for (l
= p
; l
< p
+n
; ++l
) {
676 value_oppose(tmp
, P
->Constraint
[singles
[i
][1+l
]][i
+1]);
677 for (l2
= p
; l2
< p
+n
; ++l2
) {
681 value_set_si(M
->p
[q
][0], 1);
682 Vector_Combine(P
->Constraint
[singles
[i
][1+l2
]]+1+nvar
,
683 P
->Constraint
[singles
[i
][1+l
]]+1+nvar
,
686 P
->Constraint
[singles
[i
][1+l2
]][i
+1],
689 s
[nd
].D
= Constraints2Polyhedron(M
, P
->NbRays
);
690 if (emptyQ(s
[nd
].D
)) {
691 Polyhedron_Free(s
[nd
].D
);
694 L
= ceil3(P
->Constraint
[singles
[i
][1+k
]]+1+nvar
,
696 P
->Constraint
[singles
[i
][1+k
]][i
+1]);
697 U
= ceil3(P
->Constraint
[singles
[i
][1+l
]]+1+nvar
,
699 P
->Constraint
[singles
[i
][1+l
]][i
+1]);
714 value_set_si(F
.d
, 0);
715 F
.x
.p
= new_enode(partition
, 2*nd
, -1);
716 for (k
= 0; k
< nd
; ++k
) {
717 EVALUE_SET_DOMAIN(F
.x
.p
->arr
[2*k
], s
[k
].D
);
718 value_clear(F
.x
.p
->arr
[2*k
+1].d
);
719 F
.x
.p
->arr
[2*k
+1] = s
[k
].E
;
724 free_evalue_refs(&F
);
727 value_set_si(m
->p
[dim
-nsingle
][dim
], 1);
728 P
= Polyhedron_Image(P
, m
, P
->NbConstraints
);
730 free_singles(singles
, nvar
);
734 free_evalue_refs(&mone
);
737 reduce_evalue(factor
);
743 Polyhedron
* ParamPolyhedron_Reduce(Polyhedron
*P
, unsigned nvar
,
746 return ParamPolyhedron_Reduce_mod(P
, nvar
, factor
);
749 Polyhedron
* ParamPolyhedron_Reduce(Polyhedron
*P
, unsigned nvar
,
752 Polyhedron
*R
= ParamPolyhedron_Reduce_mod(P
, nvar
, factor
);
753 evalue_mod2table(factor
, P
->Dimension
- nvar
);
754 reduce_evalue(factor
);
759 Bool
isIdentity(Matrix
*M
)
762 if (M
->NbRows
!= M
->NbColumns
)
765 for (i
= 0;i
< M
->NbRows
; i
++)
766 for (j
= 0; j
< M
->NbColumns
; j
++)
768 if(value_notone_p(M
->p
[i
][j
]))
771 if(value_notzero_p(M
->p
[i
][j
]))
777 void Param_Polyhedron_Print(FILE* DST
, Param_Polyhedron
*PP
, char **param_names
)
782 for(P
=PP
->D
;P
;P
=P
->next
) {
784 /* prints current val. dom. */
785 printf( "---------------------------------------\n" );
786 printf( "Domain :\n");
787 Print_Domain( stdout
, P
->Domain
, param_names
);
789 /* scan the vertices */
790 printf( "Vertices :\n");
791 FORALL_PVertex_in_ParamPolyhedron(V
,P
,PP
) {
793 /* prints each vertex */
794 Print_Vertex( stdout
, V
->Vertex
, param_names
);
797 END_FORALL_PVertex_in_ParamPolyhedron
;
801 void Enumeration_Print(FILE *Dst
, Enumeration
*en
, char **params
)
803 for (; en
; en
= en
->next
) {
804 Print_Domain(Dst
, en
->ValidityDomain
, params
);
805 print_evalue(Dst
, &en
->EP
, params
);
809 void Enumeration_mod2table(Enumeration
*en
, unsigned nparam
)
811 for (; en
; en
= en
->next
)
812 evalue_mod2table(&en
->EP
, nparam
);
815 void Free_ParamNames(char **params
, int m
)
822 int DomainIncludes(Polyhedron
*Pol1
, Polyhedron
*Pol2
)
825 for ( ; Pol1
; Pol1
= Pol1
->next
) {
826 for (P2
= Pol2
; P2
; P2
= P2
->next
)
827 if (!PolyhedronIncludes(Pol1
, P2
))