3 #include <barvinok/util.h>
4 #include <barvinok/options.h>
5 #include <polylib/ranking.h>
8 #define ALLOC(type) (type*)malloc(sizeof(type))
9 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
12 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
14 #define NALLOC(p,n) p = (void *)malloc((n) * sizeof(*p))
17 void manual_count(Polyhedron
*P
, Value
* result
)
19 Polyhedron
*U
= Universe_Polyhedron(0);
20 Enumeration
*en
= Polyhedron_Enumerate(P
,U
,1024,NULL
);
21 Value
*v
= compute_poly(en
,NULL
);
22 value_assign(*result
, *v
);
29 #include <barvinok/evalue.h>
30 #include <barvinok/util.h>
31 #include <barvinok/barvinok.h>
33 /* Return random value between 0 and max-1 inclusive
35 int random_int(int max
) {
36 return (int) (((double)(max
))*rand()/(RAND_MAX
+1.0));
39 Polyhedron
*Polyhedron_Read(unsigned MaxRays
)
42 unsigned NbRows
, NbColumns
;
47 while (fgets(s
, sizeof(s
), stdin
)) {
50 if (strncasecmp(s
, "vertices", sizeof("vertices")-1) == 0)
52 if (sscanf(s
, "%u %u", &NbRows
, &NbColumns
) == 2)
57 M
= Matrix_Alloc(NbRows
,NbColumns
);
60 P
= Rays2Polyhedron(M
, MaxRays
);
62 P
= Constraints2Polyhedron(M
, MaxRays
);
67 /* Inplace polarization
69 void Polyhedron_Polarize(Polyhedron
*P
)
71 unsigned NbRows
= P
->NbConstraints
+ P
->NbRays
;
75 q
= (Value
**)malloc(NbRows
* sizeof(Value
*));
77 for (i
= 0; i
< P
->NbRays
; ++i
)
79 for (; i
< NbRows
; ++i
)
80 q
[i
] = P
->Constraint
[i
-P
->NbRays
];
81 P
->NbConstraints
= NbRows
- P
->NbConstraints
;
82 P
->NbRays
= NbRows
- P
->NbRays
;
85 P
->Ray
= q
+ P
->NbConstraints
;
89 * Rather general polar
90 * We can optimize it significantly if we assume that
93 * Also, we calculate the polar as defined in Schrijver
94 * The opposite should probably work as well and would
95 * eliminate the need for multiplying by -1
97 Polyhedron
* Polyhedron_Polar(Polyhedron
*P
, unsigned NbMaxRays
)
101 unsigned dim
= P
->Dimension
+ 2;
102 Matrix
*M
= Matrix_Alloc(P
->NbRays
, dim
);
106 value_set_si(mone
, -1);
107 for (i
= 0; i
< P
->NbRays
; ++i
) {
108 Vector_Scale(P
->Ray
[i
], M
->p
[i
], mone
, dim
);
109 value_multiply(M
->p
[i
][0], M
->p
[i
][0], mone
);
110 value_multiply(M
->p
[i
][dim
-1], M
->p
[i
][dim
-1], mone
);
112 P
= Constraints2Polyhedron(M
, NbMaxRays
);
120 * Returns the supporting cone of P at the vertex with index v
122 Polyhedron
* supporting_cone(Polyhedron
*P
, int v
)
127 unsigned char *supporting
= (unsigned char *)malloc(P
->NbConstraints
);
128 unsigned dim
= P
->Dimension
+ 2;
130 assert(v
>=0 && v
< P
->NbRays
);
131 assert(value_pos_p(P
->Ray
[v
][dim
-1]));
135 for (i
= 0, n
= 0; i
< P
->NbConstraints
; ++i
) {
136 Inner_Product(P
->Constraint
[i
] + 1, P
->Ray
[v
] + 1, dim
- 1, &tmp
);
137 if ((supporting
[i
] = value_zero_p(tmp
)))
140 assert(n
>= dim
- 2);
142 M
= Matrix_Alloc(n
, dim
);
144 for (i
= 0, j
= 0; i
< P
->NbConstraints
; ++i
)
146 value_set_si(M
->p
[j
][dim
-1], 0);
147 Vector_Copy(P
->Constraint
[i
], M
->p
[j
++], dim
-1);
150 P
= Constraints2Polyhedron(M
, P
->NbRays
+1);
156 void value_lcm(const Value i
, const Value j
, Value
* lcm
)
160 value_multiply(aux
,i
,j
);
162 value_division(*lcm
,aux
,*lcm
);
166 unsigned char *supporting_constraints(Polyhedron
*P
, Param_Vertices
*v
, int *n
)
168 Value lcm
, tmp
, tmp2
;
169 unsigned dim
= P
->Dimension
+ 2;
170 unsigned nparam
= v
->Vertex
->NbColumns
- 2;
171 unsigned nvar
= dim
- nparam
- 2;
172 unsigned char *supporting
= (unsigned char *)malloc(P
->NbConstraints
);
177 row
= Vector_Alloc(nparam
+1);
182 value_set_si(lcm
, 1);
183 for (i
= 0, *n
= 0; i
< P
->NbConstraints
; ++i
) {
184 Vector_Set(row
->p
, 0, nparam
+1);
185 for (j
= 0 ; j
< nvar
; ++j
) {
186 value_set_si(tmp
, 1);
187 value_assign(tmp2
, P
->Constraint
[i
][j
+1]);
188 if (value_ne(lcm
, v
->Vertex
->p
[j
][nparam
+1])) {
189 value_assign(tmp
, lcm
);
190 value_lcm(lcm
, v
->Vertex
->p
[j
][nparam
+1], &lcm
);
191 value_division(tmp
, lcm
, tmp
);
192 value_multiply(tmp2
, tmp2
, lcm
);
193 value_division(tmp2
, tmp2
, v
->Vertex
->p
[j
][nparam
+1]);
195 Vector_Combine(row
->p
, v
->Vertex
->p
[j
], row
->p
,
196 tmp
, tmp2
, nparam
+1);
198 value_set_si(tmp
, 1);
199 Vector_Combine(row
->p
, P
->Constraint
[i
]+1+nvar
, row
->p
, tmp
, lcm
, nparam
+1);
200 for (j
= 0; j
< nparam
+1; ++j
)
201 if (value_notzero_p(row
->p
[j
]))
203 if ((supporting
[i
] = (j
== nparam
+ 1)))
215 Polyhedron
* supporting_cone_p(Polyhedron
*P
, Param_Vertices
*v
)
218 unsigned dim
= P
->Dimension
+ 2;
219 unsigned nparam
= v
->Vertex
->NbColumns
- 2;
220 unsigned nvar
= dim
- nparam
- 2;
222 unsigned char *supporting
;
224 supporting
= supporting_constraints(P
, v
, &n
);
225 M
= Matrix_Alloc(n
, nvar
+2);
227 for (i
= 0, j
= 0; i
< P
->NbConstraints
; ++i
)
229 value_set_si(M
->p
[j
][nvar
+1], 0);
230 Vector_Copy(P
->Constraint
[i
], M
->p
[j
++], nvar
+1);
233 P
= Constraints2Polyhedron(M
, P
->NbRays
+1);
239 Polyhedron
* triangulate_cone(Polyhedron
*P
, unsigned NbMaxCons
)
241 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
242 options
->MaxRays
= NbMaxCons
;
243 P
= triangulate_cone_with_options(P
, options
);
244 barvinok_options_free(options
);
248 Polyhedron
* triangulate_cone_with_options(Polyhedron
*P
,
249 struct barvinok_options
*options
)
251 const static int MAX_TRY
=10;
254 unsigned dim
= P
->Dimension
;
255 Matrix
*M
= Matrix_Alloc(P
->NbRays
+1, dim
+3);
257 Polyhedron
*L
, *R
, *T
;
258 assert(P
->NbEq
== 0);
264 Vector_Set(M
->p
[0]+1, 0, dim
+1);
265 value_set_si(M
->p
[0][0], 1);
266 value_set_si(M
->p
[0][dim
+2], 1);
267 Vector_Set(M
->p
[P
->NbRays
]+1, 0, dim
+2);
268 value_set_si(M
->p
[P
->NbRays
][0], 1);
269 value_set_si(M
->p
[P
->NbRays
][dim
+1], 1);
271 for (i
= 0, r
= 1; i
< P
->NbRays
; ++i
) {
272 if (value_notzero_p(P
->Ray
[i
][dim
+1]))
274 Vector_Copy(P
->Ray
[i
], M
->p
[r
], dim
+1);
275 value_set_si(M
->p
[r
][dim
+2], 0);
279 M2
= Matrix_Alloc(dim
+1, dim
+2);
282 if (options
->try_Delaunay_triangulation
) {
283 /* Delaunay triangulation */
284 for (r
= 1; r
< P
->NbRays
; ++r
) {
285 Inner_Product(M
->p
[r
]+1, M
->p
[r
]+1, dim
, &tmp
);
286 value_assign(M
->p
[r
][dim
+1], tmp
);
289 L
= Rays2Polyhedron(M3
, options
->MaxRays
);
294 /* Usually R should still be 0 */
297 for (r
= 1; r
< P
->NbRays
; ++r
) {
298 value_set_si(M
->p
[r
][dim
+1], random_int((t
+1)*dim
*P
->NbRays
)+1);
301 L
= Rays2Polyhedron(M3
, options
->MaxRays
);
305 assert(t
<= MAX_TRY
);
310 POL_ENSURE_FACETS(L
);
311 for (i
= 0; i
< L
->NbConstraints
; ++i
) {
312 /* Ignore perpendicular facets, i.e., facets with 0 z-coordinate */
313 if (value_negz_p(L
->Constraint
[i
][dim
+1]))
315 if (value_notzero_p(L
->Constraint
[i
][dim
+2]))
317 for (j
= 1, r
= 1; j
< M
->NbRows
; ++j
) {
318 Inner_Product(M
->p
[j
]+1, L
->Constraint
[i
]+1, dim
+1, &tmp
);
319 if (value_notzero_p(tmp
))
323 Vector_Copy(M
->p
[j
]+1, M2
->p
[r
]+1, dim
);
324 value_set_si(M2
->p
[r
][0], 1);
325 value_set_si(M2
->p
[r
][dim
+1], 0);
329 Vector_Set(M2
->p
[0]+1, 0, dim
);
330 value_set_si(M2
->p
[0][0], 1);
331 value_set_si(M2
->p
[0][dim
+1], 1);
332 T
= Rays2Polyhedron(M2
, P
->NbConstraints
+1);
346 void check_triangulization(Polyhedron
*P
, Polyhedron
*T
)
348 Polyhedron
*C
, *D
, *E
, *F
, *G
, *U
;
349 for (C
= T
; C
; C
= C
->next
) {
353 U
= DomainConvex(DomainUnion(U
, C
, 100), 100);
354 for (D
= C
->next
; D
; D
= D
->next
) {
359 E
= DomainIntersection(C
, D
, 600);
360 assert(E
->NbRays
== 0 || E
->NbEq
>= 1);
366 assert(PolyhedronIncludes(U
, P
));
367 assert(PolyhedronIncludes(P
, U
));
370 /* Computes x, y and g such that g = gcd(a,b) and a*x+b*y = g */
371 void Extended_Euclid(Value a
, Value b
, Value
*x
, Value
*y
, Value
*g
)
373 Value c
, d
, e
, f
, tmp
;
380 value_absolute(c
, a
);
381 value_absolute(d
, b
);
384 while(value_pos_p(d
)) {
385 value_division(tmp
, c
, d
);
386 value_multiply(tmp
, tmp
, f
);
387 value_subtract(e
, e
, tmp
);
388 value_division(tmp
, c
, d
);
389 value_multiply(tmp
, tmp
, d
);
390 value_subtract(c
, c
, tmp
);
397 else if (value_pos_p(a
))
399 else value_oppose(*x
, e
);
403 value_multiply(tmp
, a
, *x
);
404 value_subtract(tmp
, c
, tmp
);
405 value_division(*y
, tmp
, b
);
414 static int unimodular_complete_1(Matrix
*m
)
416 Value g
, b
, c
, old
, tmp
;
425 value_assign(g
, m
->p
[0][0]);
426 for (i
= 1; value_zero_p(g
) && i
< m
->NbColumns
; ++i
) {
427 for (j
= 0; j
< m
->NbColumns
; ++j
) {
429 value_set_si(m
->p
[i
][j
], 1);
431 value_set_si(m
->p
[i
][j
], 0);
433 value_assign(g
, m
->p
[0][i
]);
435 for (; i
< m
->NbColumns
; ++i
) {
436 value_assign(old
, g
);
437 Extended_Euclid(old
, m
->p
[0][i
], &c
, &b
, &g
);
439 for (j
= 0; j
< m
->NbColumns
; ++j
) {
441 value_multiply(tmp
, m
->p
[0][j
], b
);
442 value_division(m
->p
[i
][j
], tmp
, old
);
444 value_assign(m
->p
[i
][j
], c
);
446 value_set_si(m
->p
[i
][j
], 0);
458 int unimodular_complete(Matrix
*M
, int row
)
465 return unimodular_complete_1(M
);
467 left_hermite(M
, &H
, &Q
, &U
);
469 for (r
= 0; ok
&& r
< row
; ++r
)
470 if (value_notone_p(H
->p
[r
][r
]))
473 for (r
= row
; r
< M
->NbRows
; ++r
)
474 Vector_Copy(Q
->p
[r
], M
->p
[r
], M
->NbColumns
);
480 * Returns a full-dimensional polyhedron with the same number
481 * of integer points as P
483 Polyhedron
*remove_equalities(Polyhedron
*P
, unsigned MaxRays
)
485 Polyhedron
*Q
= Polyhedron_Copy(P
);
486 unsigned dim
= P
->Dimension
;
493 Q
= DomainConstraintSimplify(Q
, MaxRays
);
497 m1
= Matrix_Alloc(dim
, dim
);
498 for (i
= 0; i
< Q
->NbEq
; ++i
)
499 Vector_Copy(Q
->Constraint
[i
]+1, m1
->p
[i
], dim
);
501 /* m1 may not be unimodular, but we won't be throwing anything away */
502 unimodular_complete(m1
, Q
->NbEq
);
504 m2
= Matrix_Alloc(dim
+1-Q
->NbEq
, dim
+1);
505 for (i
= Q
->NbEq
; i
< dim
; ++i
)
506 Vector_Copy(m1
->p
[i
], m2
->p
[i
-Q
->NbEq
], dim
);
507 value_set_si(m2
->p
[dim
-Q
->NbEq
][dim
], 1);
510 P
= Polyhedron_Image(Q
, m2
, MaxRays
);
518 * Returns a full-dimensional polyhedron with the same number
519 * of integer points as P
520 * nvar specifies the number of variables
521 * The remaining dimensions are assumed to be parameters
523 * factor is NbEq x (nparam+2) matrix, containing stride constraints
524 * on the parameters; column nparam is the constant;
525 * column nparam+1 is the stride
527 * if factor is NULL, only remove equalities that don't affect
528 * the number of points
530 Polyhedron
*remove_equalities_p(Polyhedron
*P
, unsigned nvar
, Matrix
**factor
,
535 unsigned dim
= P
->Dimension
;
542 m1
= Matrix_Alloc(nvar
, nvar
);
543 P
= DomainConstraintSimplify(P
, MaxRays
);
545 f
= Matrix_Alloc(P
->NbEq
, dim
-nvar
+2);
549 for (i
= 0, j
= 0; i
< P
->NbEq
; ++i
) {
550 if (First_Non_Zero(P
->Constraint
[i
]+1, nvar
) == -1)
553 Vector_Gcd(P
->Constraint
[i
]+1, nvar
, &g
);
554 if (!factor
&& value_notone_p(g
))
558 Vector_Copy(P
->Constraint
[i
]+1+nvar
, f
->p
[j
], dim
-nvar
+1);
559 value_assign(f
->p
[j
][dim
-nvar
+1], g
);
562 Vector_Copy(P
->Constraint
[i
]+1, m1
->p
[j
], nvar
);
568 unimodular_complete(m1
, j
);
570 m2
= Matrix_Alloc(dim
+1-j
, dim
+1);
571 for (i
= 0; i
< nvar
-j
; ++i
)
572 Vector_Copy(m1
->p
[i
+j
], m2
->p
[i
], nvar
);
574 for (i
= nvar
-j
; i
<= dim
-j
; ++i
)
575 value_set_si(m2
->p
[i
][i
+j
], 1);
577 Q
= Polyhedron_Image(P
, m2
, MaxRays
);
584 void Line_Length(Polyhedron
*P
, Value
*len
)
590 assert(P
->Dimension
== 1);
596 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
597 value_oppose(tmp
, P
->Constraint
[i
][2]);
598 if (value_pos_p(P
->Constraint
[i
][1])) {
599 mpz_cdiv_q(tmp
, tmp
, P
->Constraint
[i
][1]);
600 if (!p
|| value_gt(tmp
, pos
))
601 value_assign(pos
, tmp
);
603 } else if (value_neg_p(P
->Constraint
[i
][1])) {
604 mpz_fdiv_q(tmp
, tmp
, P
->Constraint
[i
][1]);
605 if (!n
|| value_lt(tmp
, neg
))
606 value_assign(neg
, tmp
);
610 value_subtract(tmp
, neg
, pos
);
611 value_increment(*len
, tmp
);
613 value_set_si(*len
, -1);
622 * Factors the polyhedron P into polyhedra Q_i such that
623 * the number of integer points in P is equal to the product
624 * of the number of integer points in the individual Q_i
626 * If no factors can be found, NULL is returned.
627 * Otherwise, a linked list of the factors is returned.
629 * If there are factors and if T is not NULL, then a matrix will be
630 * returned through T expressing the old variables in terms of the
631 * new variables as they appear in the sequence of factors.
633 * The algorithm works by first computing the Hermite normal form
634 * and then grouping columns linked by one or more constraints together,
635 * where a constraints "links" two or more columns if the constraint
636 * has nonzero coefficients in the columns.
638 Polyhedron
* Polyhedron_Factor(Polyhedron
*P
, unsigned nparam
, Matrix
**T
,
642 Matrix
*M
, *H
, *Q
, *U
;
643 int *pos
; /* for each column: row position of pivot */
644 int *group
; /* group to which a column belongs */
645 int *cnt
; /* number of columns in the group */
646 int *rowgroup
; /* group to which a constraint belongs */
647 int nvar
= P
->Dimension
- nparam
;
648 Polyhedron
*F
= NULL
;
656 NALLOC(rowgroup
, P
->NbConstraints
);
658 M
= Matrix_Alloc(P
->NbConstraints
, nvar
);
659 for (i
= 0; i
< P
->NbConstraints
; ++i
)
660 Vector_Copy(P
->Constraint
[i
]+1, M
->p
[i
], nvar
);
661 left_hermite(M
, &H
, &Q
, &U
);
665 for (i
= 0; i
< P
->NbConstraints
; ++i
)
667 for (i
= 0, j
= 0; i
< H
->NbColumns
; ++i
) {
668 for ( ; j
< H
->NbRows
; ++j
)
669 if (value_notzero_p(H
->p
[j
][i
]))
671 assert (j
< H
->NbRows
);
674 for (i
= 0; i
< nvar
; ++i
) {
678 for (i
= 0; i
< H
->NbColumns
&& cnt
[0] < nvar
; ++i
) {
679 if (rowgroup
[pos
[i
]] == -1)
680 rowgroup
[pos
[i
]] = i
;
681 for (j
= pos
[i
]+1; j
< H
->NbRows
; ++j
) {
682 if (value_zero_p(H
->p
[j
][i
]))
684 if (rowgroup
[j
] != -1)
686 rowgroup
[j
] = group
[i
];
687 for (k
= i
+1; k
< H
->NbColumns
&& j
>= pos
[k
]; ++k
) {
692 if (group
[k
] != group
[i
] && value_notzero_p(H
->p
[j
][k
])) {
693 assert(cnt
[group
[k
]] != 0);
694 assert(cnt
[group
[i
]] != 0);
695 if (group
[i
] < group
[k
]) {
696 cnt
[group
[i
]] += cnt
[group
[k
]];
700 cnt
[group
[k
]] += cnt
[group
[i
]];
709 if (cnt
[0] != nvar
) {
710 /* Extract out pure context constraints separately */
711 Polyhedron
**next
= &F
;
714 *T
= Matrix_Alloc(nvar
, nvar
);
715 for (i
= nparam
? -1 : 0; i
< nvar
; ++i
) {
719 for (j
= 0, k
= 0; j
< P
->NbConstraints
; ++j
)
720 if (rowgroup
[j
] == -1) {
721 if (First_Non_Zero(P
->Constraint
[j
]+1+nvar
,
734 for (j
= 0, k
= 0; j
< P
->NbConstraints
; ++j
)
735 if (rowgroup
[j
] >= 0 && group
[rowgroup
[j
]] == i
) {
742 for (j
= 0; j
< nvar
; ++j
) {
744 for (l
= 0, m
= 0; m
< d
; ++l
) {
747 value_assign((*T
)->p
[j
][tot_d
+m
++], U
->p
[j
][l
]);
751 M
= Matrix_Alloc(k
, d
+nparam
+2);
752 for (j
= 0, k
= 0; j
< P
->NbConstraints
; ++j
) {
754 if (rowgroup
[j
] != i
)
756 value_assign(M
->p
[k
][0], P
->Constraint
[j
][0]);
757 for (l
= 0, m
= 0; m
< d
; ++l
) {
760 value_assign(M
->p
[k
][1+m
++], H
->p
[j
][l
]);
762 Vector_Copy(P
->Constraint
[j
]+1+nvar
, M
->p
[k
]+1+m
, nparam
+1);
765 *next
= Constraints2Polyhedron(M
, NbMaxRays
);
766 next
= &(*next
)->next
;
781 * Project on final dim dimensions
783 Polyhedron
* Polyhedron_Project(Polyhedron
*P
, int dim
)
786 int remove
= P
->Dimension
- dim
;
790 if (P
->Dimension
== dim
)
791 return Polyhedron_Copy(P
);
793 T
= Matrix_Alloc(dim
+1, P
->Dimension
+1);
794 for (i
= 0; i
< dim
+1; ++i
)
795 value_set_si(T
->p
[i
][i
+remove
], 1);
796 I
= Polyhedron_Image(P
, T
, P
->NbConstraints
);
801 /* Constructs a new constraint that ensures that
802 * the first constraint is (strictly) smaller than
805 static void smaller_constraint(Value
*a
, Value
*b
, Value
*c
, int pos
, int shift
,
806 int len
, int strict
, Value
*tmp
)
808 value_oppose(*tmp
, b
[pos
+1]);
809 value_set_si(c
[0], 1);
810 Vector_Combine(a
+1+shift
, b
+1+shift
, c
+1, *tmp
, a
[pos
+1], len
-shift
-1);
812 value_decrement(c
[len
-shift
-1], c
[len
-shift
-1]);
813 ConstraintSimplify(c
, c
, len
-shift
, tmp
);
817 /* For each pair of lower and upper bounds on the first variable,
818 * calls fn with the set of constraints on the remaining variables
819 * where these bounds are active, i.e., (stricly) larger/smaller than
820 * the other lower/upper bounds, the lower and upper bound and the
823 * If the first variable is equal to an affine combination of the
824 * other variables then fn is called with both lower and upper
825 * pointing to the corresponding equality.
827 void for_each_lower_upper_bound(Polyhedron
*P
, for_each_lower_upper_bound_fn fn
,
830 unsigned dim
= P
->Dimension
;
837 if (value_zero_p(P
->Constraint
[0][0]) &&
838 value_notzero_p(P
->Constraint
[0][1])) {
839 M
= Matrix_Alloc(P
->NbConstraints
-1, dim
-1+2);
840 for (i
= 1; i
< P
->NbConstraints
; ++i
) {
841 value_assign(M
->p
[i
-1][0], P
->Constraint
[i
][0]);
842 Vector_Copy(P
->Constraint
[i
]+2, M
->p
[i
-1]+1, dim
);
844 fn(M
, P
->Constraint
[0], P
->Constraint
[0], cb_data
);
850 pos
= ALLOCN(int, P
->NbConstraints
);
852 for (i
= 0, z
= 0; i
< P
->NbConstraints
; ++i
)
853 if (value_zero_p(P
->Constraint
[i
][1]))
854 pos
[P
->NbConstraints
-1 - z
++] = i
;
855 /* put those with positive coefficients first; number: p */
856 for (i
= 0, p
= 0, n
= P
->NbConstraints
-z
-1; i
< P
->NbConstraints
; ++i
)
857 if (value_pos_p(P
->Constraint
[i
][1]))
859 else if (value_neg_p(P
->Constraint
[i
][1]))
861 n
= P
->NbConstraints
-z
-p
;
862 assert (p
>= 1 && n
>= 1);
864 M
= Matrix_Alloc((p
-1) + (n
-1) + z
+ 1, dim
-1+2);
865 for (i
= 0; i
< z
; ++i
) {
866 value_assign(M
->p
[i
][0], P
->Constraint
[pos
[P
->NbConstraints
-1 - i
]][0]);
867 Vector_Copy(P
->Constraint
[pos
[P
->NbConstraints
-1 - i
]]+2,
870 for (k
= 0; k
< p
; ++k
) {
871 for (k2
= 0; k2
< p
; ++k2
) {
874 q
= 1 + z
+ k2
- (k2
> k
);
876 P
->Constraint
[pos
[k
]],
877 P
->Constraint
[pos
[k2
]],
878 M
->p
[q
], 0, 1, dim
+2, k2
> k
, &g
);
880 for (l
= p
; l
< p
+n
; ++l
) {
881 for (l2
= p
; l2
< p
+n
; ++l2
) {
884 q
= 1 + z
+ l2
-1 - (l2
> l
);
886 P
->Constraint
[pos
[l2
]],
887 P
->Constraint
[pos
[l
]],
888 M
->p
[q
], 0, 1, dim
+2, l2
> l
, &g
);
890 smaller_constraint(P
->Constraint
[pos
[k
]],
891 P
->Constraint
[pos
[l
]],
892 M
->p
[z
], 0, 1, dim
+2, 0, &g
);
893 fn(M
, P
->Constraint
[pos
[k
]], P
->Constraint
[pos
[l
]], cb_data
);
902 struct section
{ Polyhedron
* D
; evalue E
; };
912 static void PLL_cb(Matrix
*M
, Value
*lower
, Value
*upper
, void *cb_data
)
914 struct PLL_data
*data
= (struct PLL_data
*)cb_data
;
915 unsigned dim
= M
->NbColumns
-1;
921 T
= Constraints2Polyhedron(M2
, data
->MaxRays
);
923 data
->s
[data
->nd
].D
= DomainIntersection(T
, data
->C
, data
->MaxRays
);
926 POL_ENSURE_VERTICES(data
->s
[data
->nd
].D
);
927 if (emptyQ(data
->s
[data
->nd
].D
)) {
928 Polyhedron_Free(data
->s
[data
->nd
].D
);
931 L
= bv_ceil3(lower
+1+1, dim
-1+1, lower
[0+1], data
->s
[data
->nd
].D
);
932 U
= bv_ceil3(upper
+1+1, dim
-1+1, upper
[0+1], data
->s
[data
->nd
].D
);
934 eadd(&data
->mone
, U
);
935 emul(&data
->mone
, U
);
936 data
->s
[data
->nd
].E
= *U
;
943 static evalue
*ParamLine_Length_mod(Polyhedron
*P
, Polyhedron
*C
, unsigned MaxRays
)
945 unsigned dim
= P
->Dimension
;
946 unsigned nvar
= dim
- C
->Dimension
;
947 int ssize
= (P
->NbConstraints
+1) * (P
->NbConstraints
+1) / 4;
948 struct PLL_data data
;
954 value_init(data
.mone
.d
);
955 evalue_set_si(&data
.mone
, -1, 1);
957 data
.s
= ALLOCN(struct section
, ssize
);
959 data
.MaxRays
= MaxRays
;
961 for_each_lower_upper_bound(P
, PLL_cb
, &data
);
965 value_set_si(F
->d
, 0);
966 F
->x
.p
= new_enode(partition
, 2*data
.nd
, dim
-nvar
);
967 for (k
= 0; k
< data
.nd
; ++k
) {
968 EVALUE_SET_DOMAIN(F
->x
.p
->arr
[2*k
], data
.s
[k
].D
);
969 value_clear(F
->x
.p
->arr
[2*k
+1].d
);
970 F
->x
.p
->arr
[2*k
+1] = data
.s
[k
].E
;
974 free_evalue_refs(&data
.mone
);
979 evalue
* ParamLine_Length(Polyhedron
*P
, Polyhedron
*C
,
980 struct barvinok_options
*options
)
983 tmp
= ParamLine_Length_mod(P
, C
, options
->MaxRays
);
984 if (options
->lookup_table
) {
985 evalue_mod2table(tmp
, C
->Dimension
);
991 Bool
isIdentity(Matrix
*M
)
994 if (M
->NbRows
!= M
->NbColumns
)
997 for (i
= 0;i
< M
->NbRows
; i
++)
998 for (j
= 0; j
< M
->NbColumns
; j
++)
1000 if(value_notone_p(M
->p
[i
][j
]))
1003 if(value_notzero_p(M
->p
[i
][j
]))
1009 void Param_Polyhedron_Print(FILE* DST
, Param_Polyhedron
*PP
, char **param_names
)
1014 for(P
=PP
->D
;P
;P
=P
->next
) {
1016 /* prints current val. dom. */
1017 fprintf(DST
, "---------------------------------------\n");
1018 fprintf(DST
, "Domain :\n");
1019 Print_Domain(DST
, P
->Domain
, param_names
);
1021 /* scan the vertices */
1022 fprintf(DST
, "Vertices :\n");
1023 FORALL_PVertex_in_ParamPolyhedron(V
,P
,PP
) {
1025 /* prints each vertex */
1026 Print_Vertex(DST
, V
->Vertex
, param_names
);
1029 END_FORALL_PVertex_in_ParamPolyhedron
;
1033 void Enumeration_Print(FILE *Dst
, Enumeration
*en
, const char * const *params
)
1035 for (; en
; en
= en
->next
) {
1036 Print_Domain(Dst
, en
->ValidityDomain
, params
);
1037 print_evalue(Dst
, &en
->EP
, params
);
1041 void Enumeration_Free(Enumeration
*en
)
1047 free_evalue_refs( &(en
->EP
) );
1048 Domain_Free( en
->ValidityDomain
);
1055 void Enumeration_mod2table(Enumeration
*en
, unsigned nparam
)
1057 for (; en
; en
= en
->next
) {
1058 evalue_mod2table(&en
->EP
, nparam
);
1059 reduce_evalue(&en
->EP
);
1063 size_t Enumeration_size(Enumeration
*en
)
1067 for (; en
; en
= en
->next
) {
1068 s
+= domain_size(en
->ValidityDomain
);
1069 s
+= evalue_size(&en
->EP
);
1074 void Free_ParamNames(char **params
, int m
)
1081 /* Check whether every set in D2 is included in some set of D1 */
1082 int DomainIncludes(Polyhedron
*D1
, Polyhedron
*D2
)
1084 for ( ; D2
; D2
= D2
->next
) {
1086 for (P1
= D1
; P1
; P1
= P1
->next
)
1087 if (PolyhedronIncludes(P1
, D2
))
1095 int line_minmax(Polyhedron
*I
, Value
*min
, Value
*max
)
1100 value_oppose(I
->Constraint
[0][2], I
->Constraint
[0][2]);
1101 /* There should never be a remainder here */
1102 if (value_pos_p(I
->Constraint
[0][1]))
1103 mpz_fdiv_q(*min
, I
->Constraint
[0][2], I
->Constraint
[0][1]);
1105 mpz_fdiv_q(*min
, I
->Constraint
[0][2], I
->Constraint
[0][1]);
1106 value_assign(*max
, *min
);
1107 } else for (i
= 0; i
< I
->NbConstraints
; ++i
) {
1108 if (value_zero_p(I
->Constraint
[i
][1])) {
1113 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
1114 if (value_pos_p(I
->Constraint
[i
][1]))
1115 mpz_cdiv_q(*min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
1117 mpz_fdiv_q(*max
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
1125 PROCEDURES TO COMPUTE ENUMERATION. recursive procedure, recurse for
1128 @param pos index position of current loop index (1..hdim-1)
1129 @param P loop domain
1130 @param context context values for fixed indices
1131 @param exist number of existential variables
1132 @return the number of integer points in this
1136 void count_points_e (int pos
, Polyhedron
*P
, int exist
, int nparam
,
1137 Value
*context
, Value
*res
)
1142 value_set_si(*res
, 0);
1146 value_init(LB
); value_init(UB
); value_init(k
);
1150 if (lower_upper_bounds(pos
,P
,context
,&LB
,&UB
) !=0) {
1151 /* Problem if UB or LB is INFINITY */
1152 value_clear(LB
); value_clear(UB
); value_clear(k
);
1153 if (pos
> P
->Dimension
- nparam
- exist
)
1154 value_set_si(*res
, 1);
1156 value_set_si(*res
, -1);
1163 for (value_assign(k
,LB
); value_le(k
,UB
); value_increment(k
,k
)) {
1164 fprintf(stderr
, "(");
1165 for (i
=1; i
<pos
; i
++) {
1166 value_print(stderr
,P_VALUE_FMT
,context
[i
]);
1167 fprintf(stderr
,",");
1169 value_print(stderr
,P_VALUE_FMT
,k
);
1170 fprintf(stderr
,")\n");
1175 value_set_si(context
[pos
],0);
1176 if (value_lt(UB
,LB
)) {
1177 value_clear(LB
); value_clear(UB
); value_clear(k
);
1178 value_set_si(*res
, 0);
1183 value_set_si(*res
, 1);
1185 value_subtract(k
,UB
,LB
);
1186 value_add_int(k
,k
,1);
1187 value_assign(*res
, k
);
1189 value_clear(LB
); value_clear(UB
); value_clear(k
);
1193 /*-----------------------------------------------------------------*/
1194 /* Optimization idea */
1195 /* If inner loops are not a function of k (the current index) */
1196 /* i.e. if P->Constraint[i][pos]==0 for all P following this and */
1198 /* Then CNT = (UB-LB+1)*count_points(pos+1, P->next, context) */
1199 /* (skip the for loop) */
1200 /*-----------------------------------------------------------------*/
1203 value_set_si(*res
, 0);
1204 for (value_assign(k
,LB
);value_le(k
,UB
);value_increment(k
,k
)) {
1205 /* Insert k in context */
1206 value_assign(context
[pos
],k
);
1207 count_points_e(pos
+1, P
->next
, exist
, nparam
, context
, &c
);
1208 if(value_notmone_p(c
))
1209 value_addto(*res
, *res
, c
);
1211 value_set_si(*res
, -1);
1214 if (pos
> P
->Dimension
- nparam
- exist
&&
1221 fprintf(stderr
,"%d\n",CNT
);
1225 value_set_si(context
[pos
],0);
1226 value_clear(LB
); value_clear(UB
); value_clear(k
);
1228 } /* count_points_e */
1230 int DomainContains(Polyhedron
*P
, Value
*list_args
, int len
,
1231 unsigned MaxRays
, int set
)
1236 if (P
->Dimension
== len
)
1237 return in_domain(P
, list_args
);
1239 assert(set
); // assume list_args is large enough
1240 assert((P
->Dimension
- len
) % 2 == 0);
1242 for (i
= 0; i
< P
->Dimension
- len
; i
+= 2) {
1244 for (j
= 0 ; j
< P
->NbEq
; ++j
)
1245 if (value_notzero_p(P
->Constraint
[j
][1+len
+i
]))
1247 assert(j
< P
->NbEq
);
1248 value_absolute(m
, P
->Constraint
[j
][1+len
+i
]);
1249 k
= First_Non_Zero(P
->Constraint
[j
]+1, len
);
1251 assert(First_Non_Zero(P
->Constraint
[j
]+1+k
+1, len
- k
- 1) == -1);
1252 mpz_fdiv_q(list_args
[len
+i
], list_args
[k
], m
);
1253 mpz_fdiv_r(list_args
[len
+i
+1], list_args
[k
], m
);
1257 return in_domain(P
, list_args
);
1260 Polyhedron
*DomainConcat(Polyhedron
*head
, Polyhedron
*tail
)
1265 for (S
= head
; S
->next
; S
= S
->next
)
1271 evalue
*barvinok_lexsmaller_ev(Polyhedron
*P
, Polyhedron
*D
, unsigned dim
,
1272 Polyhedron
*C
, unsigned MaxRays
)
1275 Polyhedron
*RC
, *RD
, *Q
;
1276 unsigned nparam
= dim
+ C
->Dimension
;
1280 RC
= LexSmaller(P
, D
, dim
, C
, MaxRays
);
1284 exist
= RD
->Dimension
- nparam
- dim
;
1285 CA
= align_context(RC
, RD
->Dimension
, MaxRays
);
1286 Q
= DomainIntersection(RD
, CA
, MaxRays
);
1287 Polyhedron_Free(CA
);
1289 Polyhedron_Free(RC
);
1292 for (Q
= RD
; Q
; Q
= Q
->next
) {
1294 Polyhedron
*next
= Q
->next
;
1297 t
= barvinok_enumerate_e(Q
, exist
, nparam
, MaxRays
);
1303 free_evalue_refs(t
);
1315 Enumeration
*barvinok_lexsmaller(Polyhedron
*P
, Polyhedron
*D
, unsigned dim
,
1316 Polyhedron
*C
, unsigned MaxRays
)
1318 evalue
*EP
= barvinok_lexsmaller_ev(P
, D
, dim
, C
, MaxRays
);
1320 return partition2enumeration(EP
);
1323 /* "align" matrix to have nrows by inserting
1324 * the necessary number of rows and an equal number of columns in front
1326 Matrix
*align_matrix(Matrix
*M
, int nrows
)
1329 int newrows
= nrows
- M
->NbRows
;
1330 Matrix
*M2
= Matrix_Alloc(nrows
, newrows
+ M
->NbColumns
);
1331 for (i
= 0; i
< newrows
; ++i
)
1332 value_set_si(M2
->p
[i
][i
], 1);
1333 for (i
= 0; i
< M
->NbRows
; ++i
)
1334 Vector_Copy(M
->p
[i
], M2
->p
[newrows
+i
]+newrows
, M
->NbColumns
);
1338 static void print_varlist(FILE *out
, int n
, char **names
)
1342 for (i
= 0; i
< n
; ++i
) {
1345 fprintf(out
, "%s", names
[i
]);
1350 static void print_term(FILE *out
, Value v
, int pos
, int dim
, int nparam
,
1351 char **iter_names
, char **param_names
, int *first
)
1353 if (value_zero_p(v
)) {
1354 if (first
&& *first
&& pos
>= dim
+ nparam
)
1360 if (!*first
&& value_pos_p(v
))
1364 if (pos
< dim
+ nparam
) {
1365 if (value_mone_p(v
))
1367 else if (!value_one_p(v
))
1368 value_print(out
, VALUE_FMT
, v
);
1370 fprintf(out
, "%s", iter_names
[pos
]);
1372 fprintf(out
, "%s", param_names
[pos
-dim
]);
1374 value_print(out
, VALUE_FMT
, v
);
1377 char **util_generate_names(int n
, const char *prefix
)
1380 int len
= (prefix
? strlen(prefix
) : 0) + 10;
1381 char **names
= ALLOCN(char*, n
);
1383 fprintf(stderr
, "ERROR: memory overflow.\n");
1386 for (i
= 0; i
< n
; ++i
) {
1387 names
[i
] = ALLOCN(char, len
);
1389 fprintf(stderr
, "ERROR: memory overflow.\n");
1393 snprintf(names
[i
], len
, "%d", i
);
1395 snprintf(names
[i
], len
, "%s%d", prefix
, i
);
1401 void util_free_names(int n
, char **names
)
1404 for (i
= 0; i
< n
; ++i
)
1409 void Polyhedron_pprint(FILE *out
, Polyhedron
*P
, int dim
, int nparam
,
1410 char **iter_names
, char **param_names
)
1415 assert(dim
+ nparam
== P
->Dimension
);
1421 print_varlist(out
, nparam
, param_names
);
1422 fprintf(out
, " -> ");
1424 print_varlist(out
, dim
, iter_names
);
1425 fprintf(out
, " : ");
1428 fprintf(out
, "FALSE");
1429 else for (i
= 0; i
< P
->NbConstraints
; ++i
) {
1431 int v
= First_Non_Zero(P
->Constraint
[i
]+1, P
->Dimension
);
1432 if (v
== -1 && value_pos_p(P
->Constraint
[i
][0]))
1435 fprintf(out
, " && ");
1436 if (v
== -1 && value_notzero_p(P
->Constraint
[i
][1+P
->Dimension
]))
1437 fprintf(out
, "FALSE");
1438 else if (value_pos_p(P
->Constraint
[i
][v
+1])) {
1439 print_term(out
, P
->Constraint
[i
][v
+1], v
, dim
, nparam
,
1440 iter_names
, param_names
, NULL
);
1441 if (value_zero_p(P
->Constraint
[i
][0]))
1442 fprintf(out
, " = ");
1444 fprintf(out
, " >= ");
1445 for (j
= v
+1; j
<= dim
+nparam
; ++j
) {
1446 value_oppose(tmp
, P
->Constraint
[i
][1+j
]);
1447 print_term(out
, tmp
, j
, dim
, nparam
,
1448 iter_names
, param_names
, &first
);
1451 value_oppose(tmp
, P
->Constraint
[i
][1+v
]);
1452 print_term(out
, tmp
, v
, dim
, nparam
,
1453 iter_names
, param_names
, NULL
);
1454 fprintf(out
, " <= ");
1455 for (j
= v
+1; j
<= dim
+nparam
; ++j
)
1456 print_term(out
, P
->Constraint
[i
][1+j
], j
, dim
, nparam
,
1457 iter_names
, param_names
, &first
);
1461 fprintf(out
, " }\n");
1466 /* Construct a cone over P with P placed at x_d = 1, with
1467 * x_d the coordinate of an extra dimension
1469 * It's probably a mistake to depend so much on the internal
1470 * representation. We should probably simply compute the
1471 * vertices/facets first.
1473 Polyhedron
*Cone_over_Polyhedron(Polyhedron
*P
)
1475 unsigned NbConstraints
= 0;
1476 unsigned NbRays
= 0;
1480 if (POL_HAS(P
, POL_INEQUALITIES
))
1481 NbConstraints
= P
->NbConstraints
+ 1;
1482 if (POL_HAS(P
, POL_POINTS
))
1483 NbRays
= P
->NbRays
+ 1;
1485 C
= Polyhedron_Alloc(P
->Dimension
+1, NbConstraints
, NbRays
);
1486 if (POL_HAS(P
, POL_INEQUALITIES
)) {
1488 for (i
= 0; i
< P
->NbConstraints
; ++i
)
1489 Vector_Copy(P
->Constraint
[i
], C
->Constraint
[i
], P
->Dimension
+2);
1491 value_set_si(C
->Constraint
[P
->NbConstraints
][0], 1);
1492 value_set_si(C
->Constraint
[P
->NbConstraints
][1+P
->Dimension
], 1);
1494 if (POL_HAS(P
, POL_POINTS
)) {
1495 C
->NbBid
= P
->NbBid
;
1496 for (i
= 0; i
< P
->NbRays
; ++i
)
1497 Vector_Copy(P
->Ray
[i
], C
->Ray
[i
], P
->Dimension
+2);
1499 value_set_si(C
->Ray
[P
->NbRays
][0], 1);
1500 value_set_si(C
->Ray
[P
->NbRays
][1+C
->Dimension
], 1);
1502 POL_SET(C
, POL_VALID
);
1503 if (POL_HAS(P
, POL_INEQUALITIES
))
1504 POL_SET(C
, POL_INEQUALITIES
);
1505 if (POL_HAS(P
, POL_POINTS
))
1506 POL_SET(C
, POL_POINTS
);
1507 if (POL_HAS(P
, POL_VERTICES
))
1508 POL_SET(C
, POL_VERTICES
);
1512 /* Returns a (dim+nparam+1)x((dim-n)+nparam+1) matrix
1513 * mapping the transformed subspace back to the original space.
1514 * n is the number of equalities involving the variables
1515 * (i.e., not purely the parameters).
1516 * The remaining n coordinates in the transformed space would
1517 * have constant (parametric) values and are therefore not
1518 * included in the variables of the new space.
1520 Matrix
*compress_variables(Matrix
*Equalities
, unsigned nparam
)
1522 unsigned dim
= (Equalities
->NbColumns
-2) - nparam
;
1523 Matrix
*M
, *H
, *Q
, *U
, *C
, *ratH
, *invH
, *Ul
, *T1
, *T2
, *T
;
1528 for (n
= 0; n
< Equalities
->NbRows
; ++n
)
1529 if (First_Non_Zero(Equalities
->p
[n
]+1, dim
) == -1)
1532 return Identity(dim
+nparam
+1);
1534 value_set_si(mone
, -1);
1535 M
= Matrix_Alloc(n
, dim
);
1536 C
= Matrix_Alloc(n
+1, nparam
+1);
1537 for (i
= 0; i
< n
; ++i
) {
1538 Vector_Copy(Equalities
->p
[i
]+1, M
->p
[i
], dim
);
1539 Vector_Scale(Equalities
->p
[i
]+1+dim
, C
->p
[i
], mone
, nparam
+1);
1541 value_set_si(C
->p
[n
][nparam
], 1);
1542 left_hermite(M
, &H
, &Q
, &U
);
1547 ratH
= Matrix_Alloc(n
+1, n
+1);
1548 invH
= Matrix_Alloc(n
+1, n
+1);
1549 for (i
= 0; i
< n
; ++i
)
1550 Vector_Copy(H
->p
[i
], ratH
->p
[i
], n
);
1551 value_set_si(ratH
->p
[n
][n
], 1);
1552 ok
= Matrix_Inverse(ratH
, invH
);
1556 T1
= Matrix_Alloc(n
+1, nparam
+1);
1557 Matrix_Product(invH
, C
, T1
);
1560 if (value_notone_p(T1
->p
[n
][nparam
])) {
1561 for (i
= 0; i
< n
; ++i
) {
1562 if (!mpz_divisible_p(T1
->p
[i
][nparam
], T1
->p
[n
][nparam
])) {
1567 /* compress_params should have taken care of this */
1568 for (j
= 0; j
< nparam
; ++j
)
1569 assert(mpz_divisible_p(T1
->p
[i
][j
], T1
->p
[n
][nparam
]));
1570 Vector_AntiScale(T1
->p
[i
], T1
->p
[i
], T1
->p
[n
][nparam
], nparam
+1);
1572 value_set_si(T1
->p
[n
][nparam
], 1);
1574 Ul
= Matrix_Alloc(dim
+1, n
+1);
1575 for (i
= 0; i
< dim
; ++i
)
1576 Vector_Copy(U
->p
[i
], Ul
->p
[i
], n
);
1577 value_set_si(Ul
->p
[dim
][n
], 1);
1578 T2
= Matrix_Alloc(dim
+1, nparam
+1);
1579 Matrix_Product(Ul
, T1
, T2
);
1583 T
= Matrix_Alloc(dim
+nparam
+1, (dim
-n
)+nparam
+1);
1584 for (i
= 0; i
< dim
; ++i
) {
1585 Vector_Copy(U
->p
[i
]+n
, T
->p
[i
], dim
-n
);
1586 Vector_Copy(T2
->p
[i
], T
->p
[i
]+dim
-n
, nparam
+1);
1588 for (i
= 0; i
< nparam
+1; ++i
)
1589 value_set_si(T
->p
[dim
+i
][(dim
-n
)+i
], 1);
1590 assert(value_one_p(T2
->p
[dim
][nparam
]));
1597 /* Computes the left inverse of an affine embedding M and, if Eq is not NULL,
1598 * the equalities that define the affine subspace onto which M maps
1601 Matrix
*left_inverse(Matrix
*M
, Matrix
**Eq
)
1604 Matrix
*L
, *H
, *Q
, *U
, *ratH
, *invH
, *Ut
, *inv
;
1607 if (M
->NbColumns
== 1) {
1608 inv
= Matrix_Alloc(1, M
->NbRows
);
1609 value_set_si(inv
->p
[0][M
->NbRows
-1], 1);
1611 *Eq
= Matrix_Alloc(M
->NbRows
-1, 1+(M
->NbRows
-1)+1);
1612 for (i
= 0; i
< M
->NbRows
-1; ++i
) {
1613 value_oppose((*Eq
)->p
[i
][1+i
], M
->p
[M
->NbRows
-1][0]);
1614 value_assign((*Eq
)->p
[i
][1+(M
->NbRows
-1)], M
->p
[i
][0]);
1621 L
= Matrix_Alloc(M
->NbRows
-1, M
->NbColumns
-1);
1622 for (i
= 0; i
< L
->NbRows
; ++i
)
1623 Vector_Copy(M
->p
[i
], L
->p
[i
], L
->NbColumns
);
1624 right_hermite(L
, &H
, &U
, &Q
);
1627 t
= Vector_Alloc(U
->NbColumns
);
1628 for (i
= 0; i
< U
->NbColumns
; ++i
)
1629 value_oppose(t
->p
[i
], M
->p
[i
][M
->NbColumns
-1]);
1631 *Eq
= Matrix_Alloc(H
->NbRows
- H
->NbColumns
, 2 + U
->NbColumns
);
1632 for (i
= 0; i
< H
->NbRows
- H
->NbColumns
; ++i
) {
1633 Vector_Copy(U
->p
[H
->NbColumns
+i
], (*Eq
)->p
[i
]+1, U
->NbColumns
);
1634 Inner_Product(U
->p
[H
->NbColumns
+i
], t
->p
, U
->NbColumns
,
1635 (*Eq
)->p
[i
]+1+U
->NbColumns
);
1638 ratH
= Matrix_Alloc(H
->NbColumns
+1, H
->NbColumns
+1);
1639 invH
= Matrix_Alloc(H
->NbColumns
+1, H
->NbColumns
+1);
1640 for (i
= 0; i
< H
->NbColumns
; ++i
)
1641 Vector_Copy(H
->p
[i
], ratH
->p
[i
], H
->NbColumns
);
1642 value_set_si(ratH
->p
[ratH
->NbRows
-1][ratH
->NbColumns
-1], 1);
1644 ok
= Matrix_Inverse(ratH
, invH
);
1647 Ut
= Matrix_Alloc(invH
->NbRows
, U
->NbColumns
+1);
1648 for (i
= 0; i
< Ut
->NbRows
-1; ++i
) {
1649 Vector_Copy(U
->p
[i
], Ut
->p
[i
], U
->NbColumns
);
1650 Inner_Product(U
->p
[i
], t
->p
, U
->NbColumns
, &Ut
->p
[i
][Ut
->NbColumns
-1]);
1654 value_set_si(Ut
->p
[Ut
->NbRows
-1][Ut
->NbColumns
-1], 1);
1655 inv
= Matrix_Alloc(invH
->NbRows
, Ut
->NbColumns
);
1656 Matrix_Product(invH
, Ut
, inv
);
1662 /* Check whether all rays are revlex positive in the parameters
1664 int Polyhedron_has_revlex_positive_rays(Polyhedron
*P
, unsigned nparam
)
1667 for (r
= 0; r
< P
->NbRays
; ++r
) {
1669 if (value_notzero_p(P
->Ray
[r
][P
->Dimension
+1]))
1671 for (i
= P
->Dimension
-1; i
>= P
->Dimension
-nparam
; --i
) {
1672 if (value_neg_p(P
->Ray
[r
][i
+1]))
1674 if (value_pos_p(P
->Ray
[r
][i
+1]))
1677 /* A ray independent of the parameters */
1678 if (i
< P
->Dimension
-nparam
)
1684 static Polyhedron
*Recession_Cone(Polyhedron
*P
, unsigned nparam
, unsigned MaxRays
)
1687 unsigned nvar
= P
->Dimension
- nparam
;
1688 Matrix
*M
= Matrix_Alloc(P
->NbConstraints
, 1 + nvar
+ 1);
1690 for (i
= 0; i
< P
->NbConstraints
; ++i
)
1691 Vector_Copy(P
->Constraint
[i
], M
->p
[i
], 1+nvar
);
1692 R
= Constraints2Polyhedron(M
, MaxRays
);
1697 int Polyhedron_is_unbounded(Polyhedron
*P
, unsigned nparam
, unsigned MaxRays
)
1701 Polyhedron
*R
= Recession_Cone(P
, nparam
, MaxRays
);
1702 POL_ENSURE_VERTICES(R
);
1704 for (i
= 0; i
< R
->NbRays
; ++i
)
1705 if (value_zero_p(R
->Ray
[i
][1+R
->Dimension
]))
1707 is_unbounded
= R
->NbBid
> 0 || i
< R
->NbRays
;
1709 return is_unbounded
;
1712 void Vector_Oppose(Value
*p1
, Value
*p2
, unsigned len
)
1716 for (i
= 0; i
< len
; ++i
)
1717 value_oppose(p2
[i
], p1
[i
]);
1720 /* perform transposition inline; assumes M is a square matrix */
1721 void Matrix_Transposition(Matrix
*M
)
1725 assert(M
->NbRows
== M
->NbColumns
);
1726 for (i
= 0; i
< M
->NbRows
; ++i
)
1727 for (j
= i
+1; j
< M
->NbColumns
; ++j
)
1728 value_swap(M
->p
[i
][j
], M
->p
[j
][i
]);