6 #include <isl/val_gmp.h>
9 #include <isl_set_polylib.h>
10 #include <barvinok/polylib.h>
11 #include <barvinok/util.h>
12 #include <barvinok/options.h>
13 #include <polylib/ranking.h>
15 #include "lattice_point.h"
17 #define ALLOC(type) (type*)malloc(sizeof(type))
18 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
21 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
23 #define NALLOC(p,n) p = (void *)malloc((n) * sizeof(*p))
26 void manual_count(Polyhedron
*P
, Value
* result
)
28 isl_ctx
*ctx
= isl_ctx_alloc();
32 int nvar
= P
->Dimension
;
34 dim
= isl_space_set_alloc(ctx
, 0, nvar
);
35 set
= isl_set_new_from_polylib(P
, dim
);
37 v
= isl_set_count_val(set
);
38 isl_val_get_num_gmp(v
, *result
);
47 #include <barvinok/evalue.h>
48 #include <barvinok/util.h>
49 #include <barvinok/barvinok.h>
51 /* Return random value between 0 and max-1 inclusive
53 int random_int(int max
) {
54 return (int) (((double)(max
))*rand()/(RAND_MAX
+1.0));
57 Polyhedron
*Polyhedron_Read(unsigned MaxRays
)
60 unsigned NbRows
, NbColumns
;
65 while (fgets(s
, sizeof(s
), stdin
)) {
68 if (strncasecmp(s
, "vertices", sizeof("vertices")-1) == 0)
70 if (sscanf(s
, "%u %u", &NbRows
, &NbColumns
) == 2)
75 M
= Matrix_Alloc(NbRows
,NbColumns
);
78 P
= Rays2Polyhedron(M
, MaxRays
);
80 P
= Constraints2Polyhedron(M
, MaxRays
);
85 /* Inplace polarization
87 void Polyhedron_Polarize(Polyhedron
*P
)
94 POL_ENSURE_VERTICES(P
);
95 NbRows
= P
->NbConstraints
+ P
->NbRays
;
96 q
= (Value
**)malloc(NbRows
* sizeof(Value
*));
98 for (i
= 0; i
< P
->NbRays
; ++i
)
100 for (; i
< NbRows
; ++i
)
101 q
[i
] = P
->Constraint
[i
-P
->NbRays
];
102 P
->NbConstraints
= NbRows
- P
->NbConstraints
;
103 P
->NbRays
= NbRows
- P
->NbRays
;
106 P
->Ray
= q
+ P
->NbConstraints
;
110 * Rather general polar
111 * We can optimize it significantly if we assume that
114 * Also, we calculate the polar as defined in Schrijver
115 * The opposite should probably work as well and would
116 * eliminate the need for multiplying by -1
118 Polyhedron
* Polyhedron_Polar(Polyhedron
*P
, unsigned NbMaxRays
)
122 unsigned dim
= P
->Dimension
+ 2;
123 Matrix
*M
= Matrix_Alloc(P
->NbRays
, dim
);
127 value_set_si(mone
, -1);
128 for (i
= 0; i
< P
->NbRays
; ++i
) {
129 Vector_Scale(P
->Ray
[i
], M
->p
[i
], mone
, dim
);
130 value_multiply(M
->p
[i
][0], M
->p
[i
][0], mone
);
131 value_multiply(M
->p
[i
][dim
-1], M
->p
[i
][dim
-1], mone
);
133 P
= Constraints2Polyhedron(M
, NbMaxRays
);
141 * Returns the supporting cone of P at the vertex with index v
143 Polyhedron
* supporting_cone(Polyhedron
*P
, int v
)
148 unsigned char *supporting
= (unsigned char *)malloc(P
->NbConstraints
);
149 unsigned dim
= P
->Dimension
+ 2;
151 assert(v
>=0 && v
< P
->NbRays
);
152 assert(value_pos_p(P
->Ray
[v
][dim
-1]));
156 for (i
= 0, n
= 0; i
< P
->NbConstraints
; ++i
) {
157 Inner_Product(P
->Constraint
[i
] + 1, P
->Ray
[v
] + 1, dim
- 1, &tmp
);
158 if ((supporting
[i
] = value_zero_p(tmp
)))
161 assert(n
>= dim
- 2);
163 M
= Matrix_Alloc(n
, dim
);
165 for (i
= 0, j
= 0; i
< P
->NbConstraints
; ++i
)
167 value_set_si(M
->p
[j
][dim
-1], 0);
168 Vector_Copy(P
->Constraint
[i
], M
->p
[j
++], dim
-1);
171 P
= Constraints2Polyhedron(M
, P
->NbRays
+1);
177 #define INT_BITS (sizeof(unsigned) * 8)
179 unsigned *supporting_constraints(Matrix
*Constraints
, Param_Vertices
*v
, int *n
)
181 Value lcm
, tmp
, tmp2
;
182 unsigned dim
= Constraints
->NbColumns
;
183 unsigned nparam
= v
->Vertex
->NbColumns
- 2;
184 unsigned nvar
= dim
- nparam
- 2;
185 int len
= (Constraints
->NbRows
+INT_BITS
-1)/INT_BITS
;
186 unsigned *supporting
= (unsigned *)calloc(len
, sizeof(unsigned));
193 row
= Vector_Alloc(nparam
+1);
198 value_set_si(lcm
, 1);
199 for (i
= 0, *n
= 0, ix
= 0, bx
= MSB
; i
< Constraints
->NbRows
; ++i
) {
200 Vector_Set(row
->p
, 0, nparam
+1);
201 for (j
= 0 ; j
< nvar
; ++j
) {
202 value_set_si(tmp
, 1);
203 value_assign(tmp2
, Constraints
->p
[i
][j
+1]);
204 if (value_ne(lcm
, v
->Vertex
->p
[j
][nparam
+1])) {
205 value_assign(tmp
, lcm
);
206 value_lcm(lcm
, lcm
, v
->Vertex
->p
[j
][nparam
+1]);
207 value_division(tmp
, lcm
, tmp
);
208 value_multiply(tmp2
, tmp2
, lcm
);
209 value_division(tmp2
, tmp2
, v
->Vertex
->p
[j
][nparam
+1]);
211 Vector_Combine(row
->p
, v
->Vertex
->p
[j
], row
->p
,
212 tmp
, tmp2
, nparam
+1);
214 value_set_si(tmp
, 1);
215 Vector_Combine(row
->p
, Constraints
->p
[i
]+1+nvar
, row
->p
, tmp
, lcm
, nparam
+1);
216 for (j
= 0; j
< nparam
+1; ++j
)
217 if (value_notzero_p(row
->p
[j
]))
219 if (j
== nparam
+ 1) {
220 supporting
[ix
] |= bx
;
234 Polyhedron
* supporting_cone_p(Polyhedron
*P
, Param_Vertices
*v
)
237 unsigned dim
= P
->Dimension
+ 2;
238 unsigned nparam
= v
->Vertex
->NbColumns
- 2;
239 unsigned nvar
= dim
- nparam
- 2;
243 unsigned *supporting
;
246 Polyhedron_Matrix_View(P
, &View
, P
->NbConstraints
);
247 supporting
= supporting_constraints(&View
, v
, &n
);
248 M
= Matrix_Alloc(n
, nvar
+2);
250 for (i
= 0, j
= 0, ix
= 0, bx
= MSB
; i
< P
->NbConstraints
; ++i
) {
251 if (supporting
[ix
] & bx
) {
252 value_set_si(M
->p
[j
][nvar
+1], 0);
253 Vector_Copy(P
->Constraint
[i
], M
->p
[j
++], nvar
+1);
258 P
= Constraints2Polyhedron(M
, P
->NbRays
+1);
264 Polyhedron
* triangulate_cone(Polyhedron
*P
, unsigned NbMaxCons
)
266 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
267 options
->MaxRays
= NbMaxCons
;
268 P
= triangulate_cone_with_options(P
, options
);
269 barvinok_options_free(options
);
273 Polyhedron
* triangulate_cone_with_options(Polyhedron
*P
,
274 struct barvinok_options
*options
)
276 const static int MAX_TRY
=10;
279 unsigned dim
= P
->Dimension
;
280 Matrix
*M
= Matrix_Alloc(P
->NbRays
+1, dim
+3);
282 Polyhedron
*L
, *R
, *T
;
283 assert(P
->NbEq
== 0);
289 Vector_Set(M
->p
[0]+1, 0, dim
+1);
290 value_set_si(M
->p
[0][0], 1);
291 value_set_si(M
->p
[0][dim
+2], 1);
292 Vector_Set(M
->p
[P
->NbRays
]+1, 0, dim
+2);
293 value_set_si(M
->p
[P
->NbRays
][0], 1);
294 value_set_si(M
->p
[P
->NbRays
][dim
+1], 1);
296 for (i
= 0, r
= 1; i
< P
->NbRays
; ++i
) {
297 if (value_notzero_p(P
->Ray
[i
][dim
+1]))
299 Vector_Copy(P
->Ray
[i
], M
->p
[r
], dim
+1);
300 value_set_si(M
->p
[r
][dim
+2], 0);
304 M2
= Matrix_Alloc(dim
+1, dim
+2);
307 if (options
->try_Delaunay_triangulation
) {
308 /* Delaunay triangulation */
309 for (r
= 1; r
< P
->NbRays
; ++r
) {
310 Inner_Product(M
->p
[r
]+1, M
->p
[r
]+1, dim
, &tmp
);
311 value_assign(M
->p
[r
][dim
+1], tmp
);
314 L
= Rays2Polyhedron(M3
, options
->MaxRays
);
319 /* Usually R should still be 0 */
322 for (r
= 1; r
< P
->NbRays
; ++r
) {
323 value_set_si(M
->p
[r
][dim
+1], random_int((t
+1)*dim
*P
->NbRays
)+1);
326 L
= Rays2Polyhedron(M3
, options
->MaxRays
);
330 assert(t
<= MAX_TRY
);
335 POL_ENSURE_FACETS(L
);
336 for (i
= 0; i
< L
->NbConstraints
; ++i
) {
337 /* Ignore perpendicular facets, i.e., facets with 0 z-coordinate */
338 if (value_negz_p(L
->Constraint
[i
][dim
+1]))
340 if (value_notzero_p(L
->Constraint
[i
][dim
+2]))
342 for (j
= 1, r
= 1; j
< M
->NbRows
; ++j
) {
343 Inner_Product(M
->p
[j
]+1, L
->Constraint
[i
]+1, dim
+1, &tmp
);
344 if (value_notzero_p(tmp
))
348 Vector_Copy(M
->p
[j
]+1, M2
->p
[r
]+1, dim
);
349 value_set_si(M2
->p
[r
][0], 1);
350 value_set_si(M2
->p
[r
][dim
+1], 0);
354 Vector_Set(M2
->p
[0]+1, 0, dim
);
355 value_set_si(M2
->p
[0][0], 1);
356 value_set_si(M2
->p
[0][dim
+1], 1);
357 T
= Rays2Polyhedron(M2
, P
->NbConstraints
+1);
371 void check_triangulization(Polyhedron
*P
, Polyhedron
*T
)
373 Polyhedron
*C
, *D
, *E
, *F
, *G
, *U
;
374 for (C
= T
; C
; C
= C
->next
) {
378 U
= DomainConvex(DomainUnion(U
, C
, 100), 100);
379 for (D
= C
->next
; D
; D
= D
->next
) {
384 E
= DomainIntersection(C
, D
, 600);
385 assert(E
->NbRays
== 0 || E
->NbEq
>= 1);
391 assert(PolyhedronIncludes(U
, P
));
392 assert(PolyhedronIncludes(P
, U
));
395 /* Computes x, y and g such that g = gcd(a,b) and a*x+b*y = g */
396 void Extended_Euclid(Value a
, Value b
, Value
*x
, Value
*y
, Value
*g
)
398 Value c
, d
, e
, f
, tmp
;
405 value_absolute(c
, a
);
406 value_absolute(d
, b
);
409 while(value_pos_p(d
)) {
410 value_division(tmp
, c
, d
);
411 value_multiply(tmp
, tmp
, f
);
412 value_subtract(e
, e
, tmp
);
413 value_division(tmp
, c
, d
);
414 value_multiply(tmp
, tmp
, d
);
415 value_subtract(c
, c
, tmp
);
422 else if (value_pos_p(a
))
424 else value_oppose(*x
, e
);
428 value_multiply(tmp
, a
, *x
);
429 value_subtract(tmp
, c
, tmp
);
430 value_division(*y
, tmp
, b
);
439 static int unimodular_complete_1(Matrix
*m
)
441 Value g
, b
, c
, old
, tmp
;
450 value_assign(g
, m
->p
[0][0]);
451 for (i
= 1; value_zero_p(g
) && i
< m
->NbColumns
; ++i
) {
452 for (j
= 0; j
< m
->NbColumns
; ++j
) {
454 value_set_si(m
->p
[i
][j
], 1);
456 value_set_si(m
->p
[i
][j
], 0);
458 value_assign(g
, m
->p
[0][i
]);
460 for (; i
< m
->NbColumns
; ++i
) {
461 value_assign(old
, g
);
462 Extended_Euclid(old
, m
->p
[0][i
], &c
, &b
, &g
);
464 for (j
= 0; j
< m
->NbColumns
; ++j
) {
466 value_multiply(tmp
, m
->p
[0][j
], b
);
467 value_division(m
->p
[i
][j
], tmp
, old
);
469 value_assign(m
->p
[i
][j
], c
);
471 value_set_si(m
->p
[i
][j
], 0);
483 int unimodular_complete(Matrix
*M
, int row
)
490 return unimodular_complete_1(M
);
492 left_hermite(M
, &H
, &Q
, &U
);
494 for (r
= 0; ok
&& r
< row
; ++r
)
495 if (value_notone_p(H
->p
[r
][r
]))
498 for (r
= row
; r
< M
->NbRows
; ++r
)
499 Vector_Copy(Q
->p
[r
], M
->p
[r
], M
->NbColumns
);
505 * left_hermite may leave positive entries below the main diagonal in H.
506 * This function postprocesses the output of left_hermite to make
507 * the non-zero entries below the main diagonal negative.
509 void neg_left_hermite(Matrix
*A
, Matrix
**H_p
, Matrix
**Q_p
, Matrix
**U_p
)
514 left_hermite(A
, &H
, &Q
, &U
);
519 for (row
= 0, col
= 0; col
< H
->NbColumns
; ++col
, ++row
) {
520 while (value_zero_p(H
->p
[row
][col
]))
522 for (i
= 0; i
< col
; ++i
) {
523 if (value_negz_p(H
->p
[row
][i
]))
526 /* subtract column col from column i in H and U */
527 for (j
= 0; j
< H
->NbRows
; ++j
)
528 value_subtract(H
->p
[j
][i
], H
->p
[j
][i
], H
->p
[j
][col
]);
529 for (j
= 0; j
< U
->NbRows
; ++j
)
530 value_subtract(U
->p
[j
][i
], U
->p
[j
][i
], U
->p
[j
][col
]);
532 /* add row i to row col in Q */
533 for (j
= 0; j
< Q
->NbColumns
; ++j
)
534 value_addto(Q
->p
[col
][j
], Q
->p
[col
][j
], Q
->p
[i
][j
]);
540 * Returns a full-dimensional polyhedron with the same number
541 * of integer points as P
543 Polyhedron
*remove_equalities(Polyhedron
*P
, unsigned MaxRays
)
547 Polyhedron
*Q
= Polyhedron_Copy(P
);
552 Q
= DomainConstraintSimplify(Q
, MaxRays
);
556 Polyhedron_Matrix_View(Q
, &M
, Q
->NbEq
);
557 T
= compress_variables(&M
, 0);
562 P
= Polyhedron_Preimage(Q
, T
, MaxRays
);
572 * Returns a full-dimensional polyhedron with the same number
573 * of integer points as P
574 * nvar specifies the number of variables
575 * The remaining dimensions are assumed to be parameters
577 * factor is NbEq x (nparam+2) matrix, containing stride constraints
578 * on the parameters; column nparam is the constant;
579 * column nparam+1 is the stride
581 * if factor is NULL, only remove equalities that don't affect
582 * the number of points
584 Polyhedron
*remove_equalities_p(Polyhedron
*P
, unsigned nvar
, Matrix
**factor
,
589 unsigned dim
= P
->Dimension
;
596 m1
= Matrix_Alloc(nvar
, nvar
);
597 P
= DomainConstraintSimplify(P
, MaxRays
);
599 f
= Matrix_Alloc(P
->NbEq
, dim
-nvar
+2);
603 for (i
= 0, j
= 0; i
< P
->NbEq
; ++i
) {
604 if (First_Non_Zero(P
->Constraint
[i
]+1, nvar
) == -1)
607 Vector_Gcd(P
->Constraint
[i
]+1, nvar
, &g
);
608 if (!factor
&& value_notone_p(g
))
612 Vector_Copy(P
->Constraint
[i
]+1+nvar
, f
->p
[j
], dim
-nvar
+1);
613 value_assign(f
->p
[j
][dim
-nvar
+1], g
);
616 Vector_Copy(P
->Constraint
[i
]+1, m1
->p
[j
], nvar
);
622 unimodular_complete(m1
, j
);
624 m2
= Matrix_Alloc(dim
+1-j
, dim
+1);
625 for (i
= 0; i
< nvar
-j
; ++i
)
626 Vector_Copy(m1
->p
[i
+j
], m2
->p
[i
], nvar
);
628 for (i
= nvar
-j
; i
<= dim
-j
; ++i
)
629 value_set_si(m2
->p
[i
][i
+j
], 1);
631 Q
= Polyhedron_Image(P
, m2
, MaxRays
);
638 void Line_Length(Polyhedron
*P
, Value
*len
)
644 assert(P
->Dimension
== 1);
647 if (mpz_divisible_p(P
->Constraint
[0][2], P
->Constraint
[0][1]))
648 value_set_si(*len
, 1);
650 value_set_si(*len
, 0);
658 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
659 value_oppose(tmp
, P
->Constraint
[i
][2]);
660 if (value_pos_p(P
->Constraint
[i
][1])) {
661 mpz_cdiv_q(tmp
, tmp
, P
->Constraint
[i
][1]);
662 if (!p
|| value_gt(tmp
, pos
))
663 value_assign(pos
, tmp
);
665 } else if (value_neg_p(P
->Constraint
[i
][1])) {
666 mpz_fdiv_q(tmp
, tmp
, P
->Constraint
[i
][1]);
667 if (!n
|| value_lt(tmp
, neg
))
668 value_assign(neg
, tmp
);
672 value_subtract(tmp
, neg
, pos
);
673 value_increment(*len
, tmp
);
675 value_set_si(*len
, -1);
683 /* Update group[k] to the group column k belongs to.
684 * When merging two groups, only the group of the current
685 * group leader is changed. Here we change the group of
686 * the other members to also point to the group that the
687 * old group leader now points to.
689 static void update_group(int *group
, int *cnt
, int k
)
698 * Factors the polyhedron P into polyhedra Q_i such that
699 * the number of integer points in P is equal to the product
700 * of the number of integer points in the individual Q_i
702 * If no factors can be found, NULL is returned.
703 * Otherwise, a linked list of the factors is returned.
705 * If there are factors and if T is not NULL, then a matrix will be
706 * returned through T expressing the old variables in terms of the
707 * new variables as they appear in the sequence of factors.
709 * The algorithm works by first computing the Hermite normal form
710 * and then grouping columns linked by one or more constraints together,
711 * where a constraints "links" two or more columns if the constraint
712 * has nonzero coefficients in the columns.
714 Polyhedron
* Polyhedron_Factor(Polyhedron
*P
, unsigned nparam
, Matrix
**T
,
718 Matrix
*M
, *H
, *Q
, *U
;
719 int *pos
; /* for each column: row position of pivot */
720 int *group
; /* group to which a column belongs */
721 int *cnt
; /* number of columns in the group */
722 int *rowgroup
; /* group to which a constraint belongs */
723 int nvar
= P
->Dimension
- nparam
;
724 Polyhedron
*F
= NULL
;
732 NALLOC(rowgroup
, P
->NbConstraints
);
734 M
= Matrix_Alloc(P
->NbConstraints
, nvar
);
735 for (i
= 0; i
< P
->NbConstraints
; ++i
)
736 Vector_Copy(P
->Constraint
[i
]+1, M
->p
[i
], nvar
);
737 left_hermite(M
, &H
, &Q
, &U
);
741 for (i
= 0; i
< P
->NbConstraints
; ++i
)
743 for (i
= 0, j
= 0; i
< H
->NbColumns
; ++i
) {
744 for ( ; j
< H
->NbRows
; ++j
)
745 if (value_notzero_p(H
->p
[j
][i
]))
749 for (i
= 0; i
< nvar
; ++i
) {
753 for (i
= 0; i
< H
->NbColumns
&& cnt
[0] < nvar
; ++i
) {
754 if (pos
[i
] == H
->NbRows
)
755 continue; /* A line direction */
756 if (rowgroup
[pos
[i
]] == -1)
757 rowgroup
[pos
[i
]] = i
;
758 for (j
= pos
[i
]+1; j
< H
->NbRows
; ++j
) {
759 if (value_zero_p(H
->p
[j
][i
]))
761 if (rowgroup
[j
] != -1)
763 rowgroup
[j
] = group
[i
];
764 for (k
= i
+1; k
< H
->NbColumns
&& j
>= pos
[k
]; ++k
) {
765 update_group(group
, cnt
, k
);
766 update_group(group
, cnt
, i
);
767 if (group
[k
] != group
[i
] && value_notzero_p(H
->p
[j
][k
])) {
768 assert(cnt
[group
[k
]] != 0);
769 assert(cnt
[group
[i
]] != 0);
770 if (group
[i
] < group
[k
]) {
771 cnt
[group
[i
]] += cnt
[group
[k
]];
773 group
[group
[k
]] = group
[i
];
775 cnt
[group
[k
]] += cnt
[group
[i
]];
777 group
[group
[i
]] = group
[k
];
783 for (i
= 1; i
< nvar
; ++i
)
784 update_group(group
, cnt
, i
);
786 if (cnt
[0] != nvar
) {
787 /* Extract out pure context constraints separately */
788 Polyhedron
**next
= &F
;
791 *T
= Matrix_Alloc(nvar
, nvar
);
792 for (i
= nparam
? -1 : 0; i
< nvar
; ++i
) {
796 for (j
= 0, k
= 0; j
< P
->NbConstraints
; ++j
)
797 if (rowgroup
[j
] == -1) {
798 if (First_Non_Zero(P
->Constraint
[j
]+1+nvar
,
811 for (j
= 0, k
= 0; j
< P
->NbConstraints
; ++j
)
812 if (rowgroup
[j
] >= 0 && group
[rowgroup
[j
]] == i
) {
819 for (j
= 0; j
< nvar
; ++j
) {
821 for (l
= 0, m
= 0; m
< d
; ++l
) {
824 value_assign((*T
)->p
[j
][tot_d
+m
++], U
->p
[j
][l
]);
828 M
= Matrix_Alloc(k
, d
+nparam
+2);
829 for (j
= 0, k
= 0; j
< P
->NbConstraints
; ++j
) {
831 if (rowgroup
[j
] != i
)
833 value_assign(M
->p
[k
][0], P
->Constraint
[j
][0]);
834 for (l
= 0, m
= 0; m
< d
; ++l
) {
837 value_assign(M
->p
[k
][1+m
++], H
->p
[j
][l
]);
839 Vector_Copy(P
->Constraint
[j
]+1+nvar
, M
->p
[k
]+1+m
, nparam
+1);
842 *next
= Constraints2Polyhedron(M
, NbMaxRays
);
843 next
= &(*next
)->next
;
857 /* Computes the intersection of the contexts of a list of factors */
858 Polyhedron
*Factor_Context(Polyhedron
*F
, unsigned nparam
, unsigned MaxRays
)
861 Polyhedron
*C
= NULL
;
863 for (Q
= F
; Q
; Q
= Q
->next
) {
865 Polyhedron
*next
= Q
->next
;
868 if (Q
->Dimension
!= nparam
)
869 QC
= Polyhedron_Project(Q
, nparam
);
872 C
= Q
== QC
? Polyhedron_Copy(QC
) : QC
;
875 C
= DomainIntersection(C
, QC
, MaxRays
);
886 * Project on final dim dimensions
888 Polyhedron
* Polyhedron_Project(Polyhedron
*P
, int dim
)
891 int remove
= P
->Dimension
- dim
;
895 if (P
->Dimension
== dim
)
896 return Polyhedron_Copy(P
);
898 T
= Matrix_Alloc(dim
+1, P
->Dimension
+1);
899 for (i
= 0; i
< dim
+1; ++i
)
900 value_set_si(T
->p
[i
][i
+remove
], 1);
901 I
= Polyhedron_Image(P
, T
, P
->NbConstraints
);
906 /* Constructs a new constraint that ensures that
907 * the first constraint is (strictly) smaller than
910 static void smaller_constraint(Value
*a
, Value
*b
, Value
*c
, int pos
, int shift
,
911 int len
, int strict
, Value
*tmp
)
913 value_oppose(*tmp
, b
[pos
+1]);
914 value_set_si(c
[0], 1);
915 Vector_Combine(a
+1+shift
, b
+1+shift
, c
+1, *tmp
, a
[pos
+1], len
-shift
-1);
917 value_decrement(c
[len
-shift
-1], c
[len
-shift
-1]);
918 ConstraintSimplify(c
, c
, len
-shift
, tmp
);
922 /* For each pair of lower and upper bounds on the first variable,
923 * calls fn with the set of constraints on the remaining variables
924 * where these bounds are active, i.e., (stricly) larger/smaller than
925 * the other lower/upper bounds, the lower and upper bound and the
928 * If the first variable is equal to an affine combination of the
929 * other variables then fn is called with both lower and upper
930 * pointing to the corresponding equality.
932 * If there is no lower (or upper) bound, then NULL is passed
933 * as the corresponding bound.
935 void for_each_lower_upper_bound(Polyhedron
*P
,
936 for_each_lower_upper_bound_init init
,
937 for_each_lower_upper_bound_fn fn
,
940 unsigned dim
= P
->Dimension
;
947 if (value_zero_p(P
->Constraint
[0][0]) &&
948 value_notzero_p(P
->Constraint
[0][1])) {
949 M
= Matrix_Alloc(P
->NbConstraints
-1, dim
-1+2);
950 for (i
= 1; i
< P
->NbConstraints
; ++i
) {
951 value_assign(M
->p
[i
-1][0], P
->Constraint
[i
][0]);
952 Vector_Copy(P
->Constraint
[i
]+2, M
->p
[i
-1]+1, dim
);
956 fn(M
, P
->Constraint
[0], P
->Constraint
[0], cb_data
);
962 pos
= ALLOCN(int, P
->NbConstraints
);
964 for (i
= 0, z
= 0; i
< P
->NbConstraints
; ++i
)
965 if (value_zero_p(P
->Constraint
[i
][1]))
966 pos
[P
->NbConstraints
-1 - z
++] = i
;
967 /* put those with positive coefficients first; number: p */
968 for (i
= 0, p
= 0, n
= P
->NbConstraints
-z
-1; i
< P
->NbConstraints
; ++i
)
969 if (value_pos_p(P
->Constraint
[i
][1]))
971 else if (value_neg_p(P
->Constraint
[i
][1]))
973 n
= P
->NbConstraints
-z
-p
;
978 M
= Matrix_Alloc((p
? p
-1 : 0) + (n
? n
-1 : 0) + z
+ 1, dim
-1+2);
979 for (i
= 0; i
< z
; ++i
) {
980 value_assign(M
->p
[i
][0], P
->Constraint
[pos
[P
->NbConstraints
-1 - i
]][0]);
981 Vector_Copy(P
->Constraint
[pos
[P
->NbConstraints
-1 - i
]]+2,
984 for (k
= p
? 0 : -1; k
< p
; ++k
) {
985 for (k2
= 0; k2
< p
; ++k2
) {
988 q
= 1 + z
+ k2
- (k2
> k
);
990 P
->Constraint
[pos
[k
]],
991 P
->Constraint
[pos
[k2
]],
992 M
->p
[q
], 0, 1, dim
+2, k2
> k
, &g
);
994 for (l
= n
? p
: p
-1; l
< p
+n
; ++l
) {
997 for (l2
= p
; l2
< p
+n
; ++l2
) {
1000 q
= 1 + z
+ l2
-1 - (l2
> l
);
1002 P
->Constraint
[pos
[l2
]],
1003 P
->Constraint
[pos
[l
]],
1004 M
->p
[q
], 0, 1, dim
+2, l2
> l
, &g
);
1007 smaller_constraint(P
->Constraint
[pos
[k
]],
1008 P
->Constraint
[pos
[l
]],
1009 M
->p
[z
], 0, 1, dim
+2, 0, &g
);
1010 lower
= p
? P
->Constraint
[pos
[k
]] : NULL
;
1011 upper
= n
? P
->Constraint
[pos
[l
]] : NULL
;
1012 fn(M
, lower
, upper
, cb_data
);
1021 struct section
{ Polyhedron
* D
; evalue E
; };
1031 static void PLL_init(unsigned n
, void *cb_data
)
1033 struct PLL_data
*data
= (struct PLL_data
*)cb_data
;
1035 data
->s
= ALLOCN(struct section
, n
);
1038 /* Computes ceil(-coef/abs(d)) */
1039 static evalue
* bv_ceil3(Value
*coef
, int len
, Value d
, Polyhedron
*P
)
1043 Vector
*val
= Vector_Alloc(len
);
1046 Vector_Oppose(coef
, val
->p
, len
);
1047 value_absolute(t
, d
);
1049 EP
= ceiling(val
->p
, t
, len
-1, P
);
1057 static void PLL_cb(Matrix
*M
, Value
*lower
, Value
*upper
, void *cb_data
)
1059 struct PLL_data
*data
= (struct PLL_data
*)cb_data
;
1060 unsigned dim
= M
->NbColumns
-1;
1068 M2
= Matrix_Copy(M
);
1069 T
= Constraints2Polyhedron(M2
, data
->MaxRays
);
1071 data
->s
[data
->nd
].D
= DomainIntersection(T
, data
->C
, data
->MaxRays
);
1074 POL_ENSURE_VERTICES(data
->s
[data
->nd
].D
);
1075 if (emptyQ(data
->s
[data
->nd
].D
)) {
1076 Polyhedron_Free(data
->s
[data
->nd
].D
);
1079 L
= bv_ceil3(lower
+1+1, dim
-1+1, lower
[0+1], data
->s
[data
->nd
].D
);
1080 U
= bv_ceil3(upper
+1+1, dim
-1+1, upper
[0+1], data
->s
[data
->nd
].D
);
1082 eadd(&data
->mone
, U
);
1083 emul(&data
->mone
, U
);
1084 data
->s
[data
->nd
].E
= *U
;
1090 static evalue
*ParamLine_Length_mod(Polyhedron
*P
, Polyhedron
*C
, unsigned MaxRays
)
1092 unsigned dim
= P
->Dimension
;
1093 unsigned nvar
= dim
- C
->Dimension
;
1094 struct PLL_data data
;
1100 value_init(data
.mone
.d
);
1101 evalue_set_si(&data
.mone
, -1, 1);
1104 data
.MaxRays
= MaxRays
;
1106 for_each_lower_upper_bound(P
, PLL_init
, PLL_cb
, &data
);
1108 free_evalue_refs(&data
.mone
);
1112 return evalue_zero();
1117 value_set_si(F
->d
, 0);
1118 F
->x
.p
= new_enode(partition
, 2*data
.nd
, dim
-nvar
);
1119 for (k
= 0; k
< data
.nd
; ++k
) {
1120 EVALUE_SET_DOMAIN(F
->x
.p
->arr
[2*k
], data
.s
[k
].D
);
1121 value_clear(F
->x
.p
->arr
[2*k
+1].d
);
1122 F
->x
.p
->arr
[2*k
+1] = data
.s
[k
].E
;
1129 evalue
* ParamLine_Length(Polyhedron
*P
, Polyhedron
*C
,
1130 struct barvinok_options
*options
)
1133 tmp
= ParamLine_Length_mod(P
, C
, options
->MaxRays
);
1134 if (options
->lookup_table
) {
1135 evalue_mod2table(tmp
, C
->Dimension
);
1141 Bool
isIdentity(Matrix
*M
)
1144 if (M
->NbRows
!= M
->NbColumns
)
1147 for (i
= 0;i
< M
->NbRows
; i
++)
1148 for (j
= 0; j
< M
->NbColumns
; j
++)
1150 if(value_notone_p(M
->p
[i
][j
]))
1153 if(value_notzero_p(M
->p
[i
][j
]))
1159 void Param_Polyhedron_Print(FILE* DST
, Param_Polyhedron
*PP
,
1160 const char **param_names
)
1165 for(P
=PP
->D
;P
;P
=P
->next
) {
1167 /* prints current val. dom. */
1168 fprintf(DST
, "---------------------------------------\n");
1169 fprintf(DST
, "Domain :\n");
1170 Print_Domain(DST
, P
->Domain
, param_names
);
1172 /* scan the vertices */
1173 fprintf(DST
, "Vertices :\n");
1174 FORALL_PVertex_in_ParamPolyhedron(V
,P
,PP
) {
1176 /* prints each vertex */
1177 Print_Vertex(DST
, V
->Vertex
, param_names
);
1180 END_FORALL_PVertex_in_ParamPolyhedron
;
1184 void Enumeration_Print(FILE *Dst
, Enumeration
*en
, const char **params
)
1186 for (; en
; en
= en
->next
) {
1187 Print_Domain(Dst
, en
->ValidityDomain
, params
);
1188 print_evalue(Dst
, &en
->EP
, params
);
1192 void Enumeration_Free(Enumeration
*en
)
1198 free_evalue_refs( &(en
->EP
) );
1199 Domain_Free( en
->ValidityDomain
);
1206 void Enumeration_mod2table(Enumeration
*en
, unsigned nparam
)
1208 for (; en
; en
= en
->next
) {
1209 evalue_mod2table(&en
->EP
, nparam
);
1210 reduce_evalue(&en
->EP
);
1214 size_t Enumeration_size(Enumeration
*en
)
1218 for (; en
; en
= en
->next
) {
1219 s
+= domain_size(en
->ValidityDomain
);
1220 s
+= evalue_size(&en
->EP
);
1225 /* Check whether every set in D2 is included in some set of D1 */
1226 int DomainIncludes(Polyhedron
*D1
, Polyhedron
*D2
)
1228 for ( ; D2
; D2
= D2
->next
) {
1230 for (P1
= D1
; P1
; P1
= P1
->next
)
1231 if (PolyhedronIncludes(P1
, D2
))
1239 int line_minmax(Polyhedron
*I
, Value
*min
, Value
*max
)
1244 value_oppose(I
->Constraint
[0][2], I
->Constraint
[0][2]);
1245 /* There should never be a remainder here */
1246 if (value_pos_p(I
->Constraint
[0][1]))
1247 mpz_fdiv_q(*min
, I
->Constraint
[0][2], I
->Constraint
[0][1]);
1249 mpz_fdiv_q(*min
, I
->Constraint
[0][2], I
->Constraint
[0][1]);
1250 value_assign(*max
, *min
);
1251 } else for (i
= 0; i
< I
->NbConstraints
; ++i
) {
1252 if (value_zero_p(I
->Constraint
[i
][1])) {
1257 value_oppose(I
->Constraint
[i
][2], I
->Constraint
[i
][2]);
1258 if (value_pos_p(I
->Constraint
[i
][1]))
1259 mpz_cdiv_q(*min
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
1261 mpz_fdiv_q(*max
, I
->Constraint
[i
][2], I
->Constraint
[i
][1]);
1267 int DomainContains(Polyhedron
*P
, Value
*list_args
, int len
,
1268 unsigned MaxRays
, int set
)
1273 if (P
->Dimension
== len
)
1274 return in_domain(P
, list_args
);
1276 assert(set
); // assume list_args is large enough
1277 assert((P
->Dimension
- len
) % 2 == 0);
1279 for (i
= 0; i
< P
->Dimension
- len
; i
+= 2) {
1281 for (j
= 0 ; j
< P
->NbEq
; ++j
)
1282 if (value_notzero_p(P
->Constraint
[j
][1+len
+i
]))
1284 assert(j
< P
->NbEq
);
1285 value_absolute(m
, P
->Constraint
[j
][1+len
+i
]);
1286 k
= First_Non_Zero(P
->Constraint
[j
]+1, len
);
1288 assert(First_Non_Zero(P
->Constraint
[j
]+1+k
+1, len
- k
- 1) == -1);
1289 mpz_fdiv_q(list_args
[len
+i
], list_args
[k
], m
);
1290 mpz_fdiv_r(list_args
[len
+i
+1], list_args
[k
], m
);
1294 return in_domain(P
, list_args
);
1297 Polyhedron
*DomainConcat(Polyhedron
*head
, Polyhedron
*tail
)
1302 for (S
= head
; S
->next
; S
= S
->next
)
1308 evalue
*barvinok_lexsmaller_ev(Polyhedron
*P
, Polyhedron
*D
, unsigned dim
,
1309 Polyhedron
*C
, unsigned MaxRays
)
1312 Polyhedron
*RC
, *RD
, *Q
;
1313 unsigned nparam
= dim
+ C
->Dimension
;
1317 RC
= LexSmaller(P
, D
, dim
, C
, MaxRays
);
1321 exist
= RD
->Dimension
- nparam
- dim
;
1322 CA
= align_context(RC
, RD
->Dimension
, MaxRays
);
1323 Q
= DomainIntersection(RD
, CA
, MaxRays
);
1324 Polyhedron_Free(CA
);
1326 Polyhedron_Free(RC
);
1329 for (Q
= RD
; Q
; Q
= Q
->next
) {
1331 Polyhedron
*next
= Q
->next
;
1334 t
= barvinok_enumerate_e(Q
, exist
, nparam
, MaxRays
);
1351 Enumeration
*barvinok_lexsmaller(Polyhedron
*P
, Polyhedron
*D
, unsigned dim
,
1352 Polyhedron
*C
, unsigned MaxRays
)
1354 evalue
*EP
= barvinok_lexsmaller_ev(P
, D
, dim
, C
, MaxRays
);
1356 return partition2enumeration(EP
);
1359 /* "align" matrix to have nrows by inserting
1360 * the necessary number of rows and an equal number of columns in front
1362 Matrix
*align_matrix(Matrix
*M
, int nrows
)
1365 int newrows
= nrows
- M
->NbRows
;
1366 Matrix
*M2
= Matrix_Alloc(nrows
, newrows
+ M
->NbColumns
);
1367 for (i
= 0; i
< newrows
; ++i
)
1368 value_set_si(M2
->p
[i
][i
], 1);
1369 for (i
= 0; i
< M
->NbRows
; ++i
)
1370 Vector_Copy(M
->p
[i
], M2
->p
[newrows
+i
]+newrows
, M
->NbColumns
);
1374 static void print_varlist(FILE *out
, int n
, char **names
)
1378 for (i
= 0; i
< n
; ++i
) {
1381 fprintf(out
, "%s", names
[i
]);
1386 static void print_term(FILE *out
, Value v
, int pos
, int dim
, int nparam
,
1387 char **iter_names
, char **param_names
, int *first
)
1389 if (value_zero_p(v
)) {
1390 if (first
&& *first
&& pos
>= dim
+ nparam
)
1396 if (!*first
&& value_pos_p(v
))
1400 if (pos
< dim
+ nparam
) {
1401 if (value_mone_p(v
))
1403 else if (!value_one_p(v
))
1404 value_print(out
, VALUE_FMT
, v
);
1406 fprintf(out
, "%s", iter_names
[pos
]);
1408 fprintf(out
, "%s", param_names
[pos
-dim
]);
1410 value_print(out
, VALUE_FMT
, v
);
1413 char **util_generate_names(int n
, const char *prefix
)
1416 int len
= (prefix
? strlen(prefix
) : 0) + 10;
1417 char **names
= ALLOCN(char*, n
);
1419 fprintf(stderr
, "ERROR: memory overflow.\n");
1422 for (i
= 0; i
< n
; ++i
) {
1423 names
[i
] = ALLOCN(char, len
);
1425 fprintf(stderr
, "ERROR: memory overflow.\n");
1429 snprintf(names
[i
], len
, "%d", i
);
1431 snprintf(names
[i
], len
, "%s%d", prefix
, i
);
1437 void util_free_names(int n
, char **names
)
1440 for (i
= 0; i
< n
; ++i
)
1445 void Polyhedron_pprint(FILE *out
, Polyhedron
*P
, int dim
, int nparam
,
1446 char **iter_names
, char **param_names
)
1451 assert(dim
+ nparam
== P
->Dimension
);
1457 print_varlist(out
, nparam
, param_names
);
1458 fprintf(out
, " -> ");
1460 print_varlist(out
, dim
, iter_names
);
1461 fprintf(out
, " : ");
1464 fprintf(out
, "FALSE");
1465 else for (i
= 0; i
< P
->NbConstraints
; ++i
) {
1467 int v
= First_Non_Zero(P
->Constraint
[i
]+1, P
->Dimension
);
1468 if (v
== -1 && value_pos_p(P
->Constraint
[i
][0]))
1471 fprintf(out
, " && ");
1472 if (v
== -1 && value_notzero_p(P
->Constraint
[i
][1+P
->Dimension
]))
1473 fprintf(out
, "FALSE");
1474 else if (value_pos_p(P
->Constraint
[i
][v
+1])) {
1475 print_term(out
, P
->Constraint
[i
][v
+1], v
, dim
, nparam
,
1476 iter_names
, param_names
, NULL
);
1477 if (value_zero_p(P
->Constraint
[i
][0]))
1478 fprintf(out
, " = ");
1480 fprintf(out
, " >= ");
1481 for (j
= v
+1; j
<= dim
+nparam
; ++j
) {
1482 value_oppose(tmp
, P
->Constraint
[i
][1+j
]);
1483 print_term(out
, tmp
, j
, dim
, nparam
,
1484 iter_names
, param_names
, &first
);
1487 value_oppose(tmp
, P
->Constraint
[i
][1+v
]);
1488 print_term(out
, tmp
, v
, dim
, nparam
,
1489 iter_names
, param_names
, NULL
);
1490 fprintf(out
, " <= ");
1491 for (j
= v
+1; j
<= dim
+nparam
; ++j
)
1492 print_term(out
, P
->Constraint
[i
][1+j
], j
, dim
, nparam
,
1493 iter_names
, param_names
, &first
);
1497 fprintf(out
, " }\n");
1502 /* Construct a cone over P with P placed at x_d = 1, with
1503 * x_d the coordinate of an extra dimension
1505 * It's probably a mistake to depend so much on the internal
1506 * representation. We should probably simply compute the
1507 * vertices/facets first.
1509 Polyhedron
*Cone_over_Polyhedron(Polyhedron
*P
)
1511 unsigned NbConstraints
= 0;
1512 unsigned NbRays
= 0;
1516 if (POL_HAS(P
, POL_INEQUALITIES
))
1517 NbConstraints
= P
->NbConstraints
+ 1;
1518 if (POL_HAS(P
, POL_POINTS
))
1519 NbRays
= P
->NbRays
+ 1;
1521 C
= Polyhedron_Alloc(P
->Dimension
+1, NbConstraints
, NbRays
);
1522 if (POL_HAS(P
, POL_INEQUALITIES
)) {
1524 for (i
= 0; i
< P
->NbConstraints
; ++i
)
1525 Vector_Copy(P
->Constraint
[i
], C
->Constraint
[i
], P
->Dimension
+2);
1527 value_set_si(C
->Constraint
[P
->NbConstraints
][0], 1);
1528 value_set_si(C
->Constraint
[P
->NbConstraints
][1+P
->Dimension
], 1);
1530 if (POL_HAS(P
, POL_POINTS
)) {
1531 C
->NbBid
= P
->NbBid
;
1532 for (i
= 0; i
< P
->NbRays
; ++i
)
1533 Vector_Copy(P
->Ray
[i
], C
->Ray
[i
], P
->Dimension
+2);
1535 value_set_si(C
->Ray
[P
->NbRays
][0], 1);
1536 value_set_si(C
->Ray
[P
->NbRays
][1+C
->Dimension
], 1);
1538 POL_SET(C
, POL_VALID
);
1539 if (POL_HAS(P
, POL_INEQUALITIES
))
1540 POL_SET(C
, POL_INEQUALITIES
);
1541 if (POL_HAS(P
, POL_POINTS
))
1542 POL_SET(C
, POL_POINTS
);
1543 if (POL_HAS(P
, POL_VERTICES
))
1544 POL_SET(C
, POL_VERTICES
);
1548 /* Returns a (dim+nparam+1)x((dim-n)+nparam+1) matrix
1549 * mapping the transformed subspace back to the original space.
1550 * n is the number of equalities involving the variables
1551 * (i.e., not purely the parameters).
1552 * The remaining n coordinates in the transformed space would
1553 * have constant (parametric) values and are therefore not
1554 * included in the variables of the new space.
1556 Matrix
*compress_variables(Matrix
*Equalities
, unsigned nparam
)
1558 unsigned dim
= (Equalities
->NbColumns
-2) - nparam
;
1559 Matrix
*M
, *H
, *Q
, *U
, *C
, *ratH
, *invH
, *Ul
, *T1
, *T2
, *T
;
1564 for (n
= 0; n
< Equalities
->NbRows
; ++n
)
1565 if (First_Non_Zero(Equalities
->p
[n
]+1, dim
) == -1)
1568 return Identity(dim
+nparam
+1);
1570 value_set_si(mone
, -1);
1571 M
= Matrix_Alloc(n
, dim
);
1572 C
= Matrix_Alloc(n
+1, nparam
+1);
1573 for (i
= 0; i
< n
; ++i
) {
1574 Vector_Copy(Equalities
->p
[i
]+1, M
->p
[i
], dim
);
1575 Vector_Scale(Equalities
->p
[i
]+1+dim
, C
->p
[i
], mone
, nparam
+1);
1577 value_set_si(C
->p
[n
][nparam
], 1);
1578 left_hermite(M
, &H
, &Q
, &U
);
1583 ratH
= Matrix_Alloc(n
+1, n
+1);
1584 invH
= Matrix_Alloc(n
+1, n
+1);
1585 for (i
= 0; i
< n
; ++i
)
1586 Vector_Copy(H
->p
[i
], ratH
->p
[i
], n
);
1587 value_set_si(ratH
->p
[n
][n
], 1);
1588 ok
= Matrix_Inverse(ratH
, invH
);
1592 T1
= Matrix_Alloc(n
+1, nparam
+1);
1593 Matrix_Product(invH
, C
, T1
);
1596 if (value_notone_p(T1
->p
[n
][nparam
])) {
1597 for (i
= 0; i
< n
; ++i
) {
1598 if (!mpz_divisible_p(T1
->p
[i
][nparam
], T1
->p
[n
][nparam
])) {
1603 /* compress_params should have taken care of this */
1604 for (j
= 0; j
< nparam
; ++j
)
1605 assert(mpz_divisible_p(T1
->p
[i
][j
], T1
->p
[n
][nparam
]));
1606 Vector_AntiScale(T1
->p
[i
], T1
->p
[i
], T1
->p
[n
][nparam
], nparam
+1);
1608 value_set_si(T1
->p
[n
][nparam
], 1);
1610 Ul
= Matrix_Alloc(dim
+1, n
+1);
1611 for (i
= 0; i
< dim
; ++i
)
1612 Vector_Copy(U
->p
[i
], Ul
->p
[i
], n
);
1613 value_set_si(Ul
->p
[dim
][n
], 1);
1614 T2
= Matrix_Alloc(dim
+1, nparam
+1);
1615 Matrix_Product(Ul
, T1
, T2
);
1619 T
= Matrix_Alloc(dim
+nparam
+1, (dim
-n
)+nparam
+1);
1620 for (i
= 0; i
< dim
; ++i
) {
1621 Vector_Copy(U
->p
[i
]+n
, T
->p
[i
], dim
-n
);
1622 Vector_Copy(T2
->p
[i
], T
->p
[i
]+dim
-n
, nparam
+1);
1624 for (i
= 0; i
< nparam
+1; ++i
)
1625 value_set_si(T
->p
[dim
+i
][(dim
-n
)+i
], 1);
1626 assert(value_one_p(T2
->p
[dim
][nparam
]));
1633 /* Computes the left inverse of an affine embedding M and, if Eq is not NULL,
1634 * the equalities that define the affine subspace onto which M maps
1637 Matrix
*left_inverse(Matrix
*M
, Matrix
**Eq
)
1640 Matrix
*L
, *H
, *Q
, *U
, *ratH
, *invH
, *Ut
, *inv
;
1643 if (M
->NbColumns
== 1) {
1644 inv
= Matrix_Alloc(1, M
->NbRows
);
1645 value_set_si(inv
->p
[0][M
->NbRows
-1], 1);
1647 *Eq
= Matrix_Alloc(M
->NbRows
-1, 1+(M
->NbRows
-1)+1);
1648 for (i
= 0; i
< M
->NbRows
-1; ++i
) {
1649 value_oppose((*Eq
)->p
[i
][1+i
], M
->p
[M
->NbRows
-1][0]);
1650 value_assign((*Eq
)->p
[i
][1+(M
->NbRows
-1)], M
->p
[i
][0]);
1657 L
= Matrix_Alloc(M
->NbRows
-1, M
->NbColumns
-1);
1658 for (i
= 0; i
< L
->NbRows
; ++i
)
1659 Vector_Copy(M
->p
[i
], L
->p
[i
], L
->NbColumns
);
1660 right_hermite(L
, &H
, &U
, &Q
);
1663 t
= Vector_Alloc(U
->NbColumns
);
1664 for (i
= 0; i
< U
->NbColumns
; ++i
)
1665 value_oppose(t
->p
[i
], M
->p
[i
][M
->NbColumns
-1]);
1667 *Eq
= Matrix_Alloc(H
->NbRows
- H
->NbColumns
, 2 + U
->NbColumns
);
1668 for (i
= 0; i
< H
->NbRows
- H
->NbColumns
; ++i
) {
1669 Vector_Copy(U
->p
[H
->NbColumns
+i
], (*Eq
)->p
[i
]+1, U
->NbColumns
);
1670 Inner_Product(U
->p
[H
->NbColumns
+i
], t
->p
, U
->NbColumns
,
1671 (*Eq
)->p
[i
]+1+U
->NbColumns
);
1674 ratH
= Matrix_Alloc(H
->NbColumns
+1, H
->NbColumns
+1);
1675 invH
= Matrix_Alloc(H
->NbColumns
+1, H
->NbColumns
+1);
1676 for (i
= 0; i
< H
->NbColumns
; ++i
)
1677 Vector_Copy(H
->p
[i
], ratH
->p
[i
], H
->NbColumns
);
1678 value_set_si(ratH
->p
[ratH
->NbRows
-1][ratH
->NbColumns
-1], 1);
1680 ok
= Matrix_Inverse(ratH
, invH
);
1683 Ut
= Matrix_Alloc(invH
->NbRows
, U
->NbColumns
+1);
1684 for (i
= 0; i
< Ut
->NbRows
-1; ++i
) {
1685 Vector_Copy(U
->p
[i
], Ut
->p
[i
], U
->NbColumns
);
1686 Inner_Product(U
->p
[i
], t
->p
, U
->NbColumns
, &Ut
->p
[i
][Ut
->NbColumns
-1]);
1690 value_set_si(Ut
->p
[Ut
->NbRows
-1][Ut
->NbColumns
-1], 1);
1691 inv
= Matrix_Alloc(invH
->NbRows
, Ut
->NbColumns
);
1692 Matrix_Product(invH
, Ut
, inv
);
1698 /* Check whether all rays are revlex positive in the parameters
1700 int Polyhedron_has_revlex_positive_rays(Polyhedron
*P
, unsigned nparam
)
1703 for (r
= 0; r
< P
->NbRays
; ++r
) {
1705 if (value_notzero_p(P
->Ray
[r
][P
->Dimension
+1]))
1707 for (i
= P
->Dimension
-1; i
>= P
->Dimension
-nparam
; --i
) {
1708 if (value_neg_p(P
->Ray
[r
][i
+1]))
1710 if (value_pos_p(P
->Ray
[r
][i
+1]))
1713 /* A ray independent of the parameters */
1714 if (i
< P
->Dimension
-nparam
)
1720 static Polyhedron
*Recession_Cone(Polyhedron
*P
, unsigned nparam
, unsigned MaxRays
)
1723 unsigned nvar
= P
->Dimension
- nparam
;
1724 Matrix
*M
= Matrix_Alloc(P
->NbConstraints
, 1 + nvar
+ 1);
1726 for (i
= 0; i
< P
->NbConstraints
; ++i
)
1727 Vector_Copy(P
->Constraint
[i
], M
->p
[i
], 1+nvar
);
1728 R
= Constraints2Polyhedron(M
, MaxRays
);
1733 int Polyhedron_is_unbounded(Polyhedron
*P
, unsigned nparam
, unsigned MaxRays
)
1737 Polyhedron
*R
= Recession_Cone(P
, nparam
, MaxRays
);
1738 POL_ENSURE_VERTICES(R
);
1740 for (i
= 0; i
< R
->NbRays
; ++i
)
1741 if (value_zero_p(R
->Ray
[i
][1+R
->Dimension
]))
1743 is_unbounded
= R
->NbBid
> 0 || i
< R
->NbRays
;
1745 return is_unbounded
;
1748 static void SwapColumns(Value
**V
, int n
, int i
, int j
)
1752 for (r
= 0; r
< n
; ++r
)
1753 value_swap(V
[r
][i
], V
[r
][j
]);
1756 void Polyhedron_ExchangeColumns(Polyhedron
*P
, int Column1
, int Column2
)
1758 SwapColumns(P
->Constraint
, P
->NbConstraints
, Column1
, Column2
);
1759 SwapColumns(P
->Ray
, P
->NbRays
, Column1
, Column2
);
1762 Polyhedron_Matrix_View(P
, &M
, P
->NbConstraints
);
1763 Gauss(&M
, P
->NbEq
, P
->Dimension
+1);
1767 /* perform transposition inline; assumes M is a square matrix */
1768 void Matrix_Transposition(Matrix
*M
)
1772 assert(M
->NbRows
== M
->NbColumns
);
1773 for (i
= 0; i
< M
->NbRows
; ++i
)
1774 for (j
= i
+1; j
< M
->NbColumns
; ++j
)
1775 value_swap(M
->p
[i
][j
], M
->p
[j
][i
]);
1778 /* Matrix "view" of first rows rows */
1779 void Polyhedron_Matrix_View(Polyhedron
*P
, Matrix
*M
, unsigned rows
)
1782 M
->NbColumns
= P
->Dimension
+2;
1783 M
->p_Init
= P
->p_Init
;
1784 M
->p
= P
->Constraint
;
1787 int Last_Non_Zero(Value
*p
, unsigned len
)
1791 for (i
= len
- 1; i
>= 0; --i
)
1792 if (value_notzero_p(p
[i
]))