1 #include <polylib/polylibgmp.h>
2 #include <barvinok/util.h>
3 #include "basis_reduction.h"
6 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
8 /* If P has no rays, then we return NULL.
9 * Otherwise, transform the polyhedron such that one of the rays
10 * is the first unit vector and cut it off at a height that ensures
11 * that if the whole polyhedron has any points, then the remaining part
12 * has integer points. In particular we add the largest coefficient
13 * of a ray to the highest vertex (rounded up).
15 * The matrix that transforms the resulting polytope to part of the
16 * original polyhedron is returned through T.
18 static Polyhedron
*remove_ray(Polyhedron
*P
, Matrix
**T
, unsigned MaxRays
)
31 for (; r
< P
->NbRays
; ++r
)
32 if (value_zero_p(P
->Ray
[r
][P
->Dimension
+1]))
34 if (P
->NbBid
== 0 && r
== P
->NbRays
)
38 v
= Vector_Alloc(P
->Dimension
+1);
39 Vector_Gcd(P
->Ray
[r
]+1, P
->Dimension
, &g
);
40 Vector_AntiScale(P
->Ray
[r
]+1, v
->p
, g
, P
->Dimension
+1);
41 M
= unimodular_complete(v
);
42 value_set_si(M
->p
[P
->Dimension
][P
->Dimension
], 1);
45 P
= Polyhedron_Preimage(P
, M2
, 0);
54 value_set_si(size
, 0);
56 for (i
= 0; i
< P
->NbBid
; ++i
) {
57 value_absolute(tmp
, P
->Ray
[i
][1]);
58 if (value_gt(tmp
, size
))
59 value_assign(size
, tmp
);
61 for (i
= P
->NbBid
; i
< P
->NbRays
; ++i
) {
62 if (value_zero_p(P
->Ray
[i
][P
->Dimension
+1])) {
63 if (value_gt(P
->Ray
[i
][1], size
))
64 value_assign(size
, P
->Ray
[i
][1]);
67 mpz_cdiv_q(tmp
, P
->Ray
[i
][1], P
->Ray
[i
][P
->Dimension
+1]);
68 if (first
|| value_gt(tmp
, offset
)) {
69 value_assign(offset
, tmp
);
73 value_addto(offset
, offset
, size
);
77 v
= Vector_Alloc(P
->Dimension
+2);
78 value_set_si(v
->p
[0], 1);
79 value_set_si(v
->p
[1], -1);
80 value_assign(v
->p
[1+P
->Dimension
], offset
);
81 R
= AddConstraints(v
->p
, 1, P
, MaxRays
);
89 static void print_minmax(Polyhedron
*P
)
92 POL_ENSURE_VERTICES(P
);
93 Polyhedron_Print(stderr
, P_VALUE_FMT
, P
);
94 for (i
= 0; i
< P
->Dimension
; ++i
) {
100 mpz_cdiv_q(min
, P
->Ray
[0][1+i
], P
->Ray
[0][1+P
->Dimension
]);
101 mpz_fdiv_q(max
, P
->Ray
[0][1+i
], P
->Ray
[0][1+P
->Dimension
]);
103 for (j
= 1; j
< P
->NbRays
; ++j
) {
104 mpz_cdiv_q(tmp
, P
->Ray
[j
][1+i
], P
->Ray
[j
][1+P
->Dimension
]);
105 if (value_lt(tmp
, min
))
106 value_assign(min
, tmp
);
107 mpz_fdiv_q(tmp
, P
->Ray
[j
][1+i
], P
->Ray
[j
][1+P
->Dimension
]);
108 if (value_gt(tmp
, max
))
109 value_assign(max
, tmp
);
111 fprintf(stderr
, "i: %d, min: ", i
);
112 value_print(stderr
, VALUE_FMT
, min
);
113 fprintf(stderr
, ", max: ");
114 value_print(stderr
, VALUE_FMT
, max
);
115 fprintf(stderr
, "\n");
123 /* Remove coordinates that have a fixed value and return the matrix
124 * that adds these fixed coordinates again through T.
126 static Polyhedron
*Polyhedron_RemoveFixedColumns(Polyhedron
*P
, Matrix
**T
)
129 int dim
= P
->Dimension
;
130 int *remove
= ALLOCN(int, dim
);
133 assert(POL_HAS(P
, POL_INEQUALITIES
));
134 for (i
= 0; i
< dim
; ++i
)
136 for (i
= 0; i
< P
->NbEq
; ++i
) {
137 int pos
= First_Non_Zero(P
->Constraint
[i
]+1, dim
);
138 assert(First_Non_Zero(P
->Constraint
[i
]+1+pos
+1, dim
-pos
-1) == -1);
141 Q
= Polyhedron_Alloc(P
->Dimension
-P
->NbEq
, P
->NbConstraints
-P
->NbEq
, P
->NbRays
);
142 for (i
= 0; i
< Q
->NbConstraints
; ++i
) {
143 value_assign(Q
->Constraint
[i
][0], P
->Constraint
[P
->NbEq
+i
][0]);
144 for (j
= 0, n
= 0; j
< P
->Dimension
; ++j
) {
148 value_assign(Q
->Constraint
[i
][1+j
-n
], P
->Constraint
[P
->NbEq
+i
][1+j
]);
150 value_assign(Q
->Constraint
[i
][1+j
-n
], P
->Constraint
[P
->NbEq
+i
][1+j
]);
152 for (i
= 0; i
< Q
->NbRays
; ++i
) {
153 value_assign(Q
->Ray
[i
][0], P
->Ray
[i
][0]);
154 for (j
= 0, n
= 0; j
< P
->Dimension
; ++j
) {
158 value_assign(Q
->Ray
[i
][1+j
-n
], P
->Ray
[i
][1+j
]);
160 value_assign(Q
->Ray
[i
][1+j
-n
], P
->Ray
[i
][1+j
]);
162 *T
= Matrix_Alloc(P
->Dimension
+1, Q
->Dimension
+1);
163 for (i
= 0, n
= 0; i
< P
->Dimension
; ++i
) {
165 value_oppose((*T
)->p
[i
][Q
->Dimension
], P
->Constraint
[n
][1+P
->Dimension
]);
168 value_set_si((*T
)->p
[i
][i
-n
], 1);
170 value_set_si((*T
)->p
[i
][i
-n
], 1);
171 POL_SET(Q
, POL_VALID
);
172 if (POL_HAS(P
, POL_INEQUALITIES
))
173 POL_SET(Q
, POL_INEQUALITIES
);
174 if (POL_HAS(P
, POL_FACETS
))
175 POL_SET(Q
, POL_FACETS
);
176 if (POL_HAS(P
, POL_POINTS
))
177 POL_SET(Q
, POL_POINTS
);
178 if (POL_HAS(P
, POL_VERTICES
))
179 POL_SET(Q
, POL_VERTICES
);
184 /* This function implements the algorithm described in
185 * "An Implementation of the Generalized Basis Reduction Algorithm
186 * for Integer Programming" of Cook el al. to find an integer point
188 * If the polyhedron is unbounded, we first remove its rays.
190 Vector
*Polyhedron_Sample(Polyhedron
*P
, unsigned MaxRays
)
193 Vector
*sample
= NULL
;
200 POL_ENSURE_VERTICES(P
);
202 for (i
= 0; i
< P
->NbRays
; ++i
)
203 if (value_one_p(P
->Ray
[i
][1+P
->Dimension
])) {
204 sample
= Vector_Alloc(P
->Dimension
+1);
205 Vector_Copy(P
->Ray
[i
]+1, sample
->p
, P
->Dimension
+1);
210 assert(P
->NbEq
== 0);
212 Q
= remove_ray(P
, &T
, MaxRays
);
216 Q_sample
= Polyhedron_Sample(Q
, MaxRays
);
220 sample
= Vector_Alloc(P
->Dimension
+ 1);
221 Matrix_Vector_Product(T
, Q_sample
->p
, sample
->p
);
222 Vector_Free(Q_sample
);
229 Matrix
*basis
= reduced_basis(P
);
231 T
= Matrix_Alloc(P
->Dimension
+1, P
->Dimension
+1);
232 inv
= Matrix_Alloc(P
->Dimension
+1, P
->Dimension
+1);
233 for (i
= 0; i
< P
->Dimension
; ++i
)
234 for (j
= 0; j
< P
->Dimension
; ++j
)
235 value_assign(T
->p
[i
][j
], basis
->p
[i
][j
]);
236 value_set_si(T
->p
[P
->Dimension
][P
->Dimension
], 1);
240 ok
= Matrix_Inverse(M
, inv
);
244 Q
= Polyhedron_Image(P
, T
, MaxRays
);
250 mpz_cdiv_q(min
, Q
->Ray
[0][1], Q
->Ray
[0][1+Q
->Dimension
]);
251 mpz_fdiv_q(max
, Q
->Ray
[0][1], Q
->Ray
[0][1+Q
->Dimension
]);
253 for (j
= 1; j
< Q
->NbRays
; ++j
) {
254 mpz_cdiv_q(tmp
, Q
->Ray
[j
][1], Q
->Ray
[j
][1+Q
->Dimension
]);
255 if (value_lt(tmp
, min
))
256 value_assign(min
, tmp
);
257 mpz_fdiv_q(tmp
, Q
->Ray
[j
][1], Q
->Ray
[j
][1+Q
->Dimension
]);
258 if (value_gt(tmp
, max
))
259 value_assign(max
, tmp
);
262 v
= Vector_Alloc(1+Q
->Dimension
+1);
263 value_set_si(v
->p
[1], -1);
265 for (value_assign(tmp
, min
); value_le(tmp
, max
); value_increment(tmp
, tmp
)) {
269 value_assign(v
->p
[1+Q
->Dimension
], tmp
);
271 R
= AddConstraints(v
->p
, 1, Q
, MaxRays
);
272 S
= Polyhedron_RemoveFixedColumns(R
, &T
);
274 S_sample
= Polyhedron_Sample(S
, MaxRays
);
277 Vector
*Q_sample
= Vector_Alloc(Q
->Dimension
+ 1);
278 Matrix_Vector_Product(T
, S_sample
->p
, Q_sample
->p
);
280 Vector_Free(S_sample
);
281 sample
= Vector_Alloc(P
->Dimension
+ 1);
282 Matrix_Vector_Product(inv
, Q_sample
->p
, sample
->p
);
283 Vector_Free(Q_sample
);