1 #include <barvinok/util.h>
2 #include <barvinok/basis_reduction.h>
3 #include <barvinok/sample.h>
4 #include <barvinok/options.h>
7 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
9 /* If P has no rays, then we return NULL.
10 * Otherwise, look for the coordinate axis with the smallest maximal non-zero
11 * coefficient over all rays and a constraint that bounds the values on
12 * this axis to the maximal value over the vertices plus the above maximal
13 * non-zero coefficient times the number of rays minus 1.
14 * Any integer point x outside this region is the sum of a point inside
15 * the region and an integer multiple of the rays.
16 * Write x = \sum_i a_i v_i + \sum_j b_j r_j
17 * with \sum_i a_i = 1.
18 * Then x = \sum_i a_i v_i + \sum_j {b_j} r_j + \sum_j [b_j] r_j
19 * with y = \sum_i a_i v_i + \sum_j {b_j} r_j a point inside the region.
21 static Polyhedron
*remove_ray(Polyhedron
*P
, unsigned MaxRays
)
24 Vector
*min
, *max
, *c
;
31 POL_ENSURE_VERTICES(P
);
33 for (; r
< P
->NbRays
; ++r
)
34 if (value_zero_p(P
->Ray
[r
][P
->Dimension
+1]))
36 if (P
->NbBid
== 0 && r
== P
->NbRays
)
39 max
= Vector_Alloc(P
->Dimension
);
40 min
= Vector_Alloc(P
->Dimension
);
41 for (r
= 0; r
< P
->NbBid
; ++r
)
42 for (i
= 0 ; i
< P
->Dimension
; ++i
)
43 if (value_abs_gt(P
->Ray
[r
][1+i
], max
->p
[i
]))
44 value_absolute(max
->p
[i
], P
->Ray
[r
][1+i
]);
46 for (i
= 0 ; i
< P
->Dimension
; ++i
)
47 value_oppose(min
->p
[i
], max
->p
[i
]);
50 for (r
= P
->NbBid
; r
< P
->NbRays
; ++r
) {
51 if (value_notzero_p(P
->Ray
[r
][P
->Dimension
+1]))
53 for (i
= 0 ; i
< P
->Dimension
; ++i
) {
54 if (value_gt(P
->Ray
[r
][1+i
], max
->p
[i
]))
55 value_assign(max
->p
[i
], P
->Ray
[r
][1+i
]);
56 if (value_lt(P
->Ray
[r
][1+i
], min
->p
[i
]))
57 value_assign(min
->p
[i
], P
->Ray
[r
][1+i
]);
66 for (i
= 0 ; i
< P
->Dimension
; ++i
) {
67 if (value_notzero_p(min
->p
[i
]) &&
68 (value_zero_p(s
) || value_abs_lt(min
->p
[i
], s
))) {
69 value_assign(s
, min
->p
[i
]);
72 if (value_notzero_p(max
->p
[i
]) &&
73 (value_zero_p(s
) || value_abs_lt(max
->p
[i
], s
))) {
74 value_assign(s
, max
->p
[i
]);
79 for (r
= P
->NbBid
; r
< P
->NbRays
; ++r
)
80 if (value_notzero_p(P
->Ray
[r
][P
->Dimension
+1]))
84 mpz_cdiv_q(v
, P
->Ray
[r
][1+pos
], P
->Ray
[r
][P
->Dimension
+1]);
86 mpz_fdiv_q(v
, P
->Ray
[r
][1+pos
], P
->Ray
[r
][P
->Dimension
+1]);
88 for ( ; r
< P
->NbRays
; ++r
) {
89 if (value_zero_p(P
->Ray
[r
][P
->Dimension
+1]))
93 mpz_cdiv_q(tmp
, P
->Ray
[r
][1+pos
], P
->Ray
[r
][P
->Dimension
+1]);
97 mpz_fdiv_q(tmp
, P
->Ray
[r
][1+pos
], P
->Ray
[r
][P
->Dimension
+1]);
103 c
= Vector_Alloc(1+P
->Dimension
+1);
105 value_set_si(tmp
, rays
);
106 value_addmul(v
, tmp
, s
);
107 value_set_si(c
->p
[0], 1);
108 if (value_pos_p(s
)) {
109 value_set_si(c
->p
[1+pos
], -1);
110 value_assign(c
->p
[1+P
->Dimension
], v
);
112 value_set_si(c
->p
[1+pos
], 1);
113 value_oppose(c
->p
[1+P
->Dimension
], v
);
115 value_decrement(c
->p
[1+P
->Dimension
], c
->p
[1+P
->Dimension
]);
117 R
= AddConstraints(c
->p
, 1, P
, MaxRays
);
131 static void print_minmax(Polyhedron
*P
)
134 POL_ENSURE_VERTICES(P
);
135 Polyhedron_Print(stderr
, P_VALUE_FMT
, P
);
136 for (i
= 0; i
< P
->Dimension
; ++i
) {
142 mpz_cdiv_q(min
, P
->Ray
[0][1+i
], P
->Ray
[0][1+P
->Dimension
]);
143 mpz_fdiv_q(max
, P
->Ray
[0][1+i
], P
->Ray
[0][1+P
->Dimension
]);
145 for (j
= 1; j
< P
->NbRays
; ++j
) {
146 mpz_cdiv_q(tmp
, P
->Ray
[j
][1+i
], P
->Ray
[j
][1+P
->Dimension
]);
147 if (value_lt(tmp
, min
))
148 value_assign(min
, tmp
);
149 mpz_fdiv_q(tmp
, P
->Ray
[j
][1+i
], P
->Ray
[j
][1+P
->Dimension
]);
150 if (value_gt(tmp
, max
))
151 value_assign(max
, tmp
);
153 fprintf(stderr
, "i: %d, min: ", i
);
154 value_print(stderr
, VALUE_FMT
, min
);
155 fprintf(stderr
, ", max: ");
156 value_print(stderr
, VALUE_FMT
, max
);
157 fprintf(stderr
, "\n");
165 /* Remove coordinates that have a fixed value and return the matrix
166 * that adds these fixed coordinates again through T.
168 static Polyhedron
*Polyhedron_RemoveFixedColumns(Polyhedron
*P
, Matrix
**T
)
171 int dim
= P
->Dimension
;
172 int *remove
= ALLOCN(int, dim
);
176 assert(POL_HAS(P
, POL_INEQUALITIES
));
177 for (i
= 0; i
< dim
; ++i
)
180 for (i
= 0; i
< P
->NbEq
; ++i
) {
181 int pos
= First_Non_Zero(P
->Constraint
[i
]+1, dim
);
182 if (First_Non_Zero(P
->Constraint
[i
]+1+pos
+1, dim
-pos
-1) != -1)
188 Q
= Polyhedron_Alloc(P
->Dimension
-NbEq
, P
->NbConstraints
-NbEq
, P
->NbRays
);
189 Q
->NbEq
= P
->NbEq
- NbEq
;
190 for (i
= 0, k
= 0; i
< P
->NbConstraints
; ++i
) {
192 int pos
= First_Non_Zero(P
->Constraint
[i
]+1, dim
);
193 if (First_Non_Zero(P
->Constraint
[i
]+1+pos
+1, dim
-pos
-1) == -1)
196 value_assign(Q
->Constraint
[k
][0], P
->Constraint
[i
][0]);
197 for (j
= 0, n
= 0; j
< P
->Dimension
; ++j
) {
201 value_assign(Q
->Constraint
[k
][1+j
-n
], P
->Constraint
[i
][1+j
]);
203 value_assign(Q
->Constraint
[k
][1+j
-n
], P
->Constraint
[i
][1+j
]);
206 for (i
= 0; i
< Q
->NbRays
; ++i
) {
207 value_assign(Q
->Ray
[i
][0], P
->Ray
[i
][0]);
208 for (j
= 0, n
= 0; j
< P
->Dimension
; ++j
) {
212 value_assign(Q
->Ray
[i
][1+j
-n
], P
->Ray
[i
][1+j
]);
214 value_assign(Q
->Ray
[i
][1+j
-n
], P
->Ray
[i
][1+j
]);
216 *T
= Matrix_Alloc(P
->Dimension
+1, Q
->Dimension
+1);
217 for (i
= 0, n
= 0; i
< P
->Dimension
; ++i
) {
218 if (remove
[i
] != -1) {
219 value_oppose((*T
)->p
[i
][Q
->Dimension
],
220 P
->Constraint
[remove
[i
]][1+P
->Dimension
]);
223 value_set_si((*T
)->p
[i
][i
-n
], 1);
225 value_set_si((*T
)->p
[i
][i
-n
], 1);
226 POL_SET(Q
, POL_VALID
);
227 if (POL_HAS(P
, POL_INEQUALITIES
))
228 POL_SET(Q
, POL_INEQUALITIES
);
229 if (POL_HAS(P
, POL_FACETS
))
230 POL_SET(Q
, POL_FACETS
);
231 if (POL_HAS(P
, POL_POINTS
))
232 POL_SET(Q
, POL_POINTS
);
233 if (POL_HAS(P
, POL_VERTICES
))
234 POL_SET(Q
, POL_VERTICES
);
239 static Polyhedron
*remove_all_equalities(Polyhedron
*P
, Matrix
**T
,
242 /* Matrix "view" of equalities */
245 M
.NbColumns
= P
->Dimension
+2;
246 M
.p_Init
= P
->p_Init
;
249 *T
= compress_variables(&M
, 0);
253 P
= Polyhedron_Preimage(P
, *T
, MaxRays
);
258 static Vector
*product_sample(Polyhedron
*P
, Matrix
*T
,
259 struct barvinok_options
*options
)
262 Vector
*sample
= NULL
;
263 Vector
*tmp
= Vector_Alloc(T
->NbRows
);
265 for (; P
; P
= P
->next
) {
267 Polyhedron
*next
= P
->next
;
269 P_sample
= Polyhedron_Sample(P
, options
);
275 Vector_Copy(P_sample
->p
, tmp
->p
+i
, P
->Dimension
);
276 Vector_Free(P_sample
);
281 sample
= Vector_Alloc(T
->NbRows
+ 1);
282 Matrix_Vector_Product(T
, tmp
->p
, sample
->p
);
283 value_set_si(sample
->p
[T
->NbRows
], 1);
289 /* This function implements the algorithm described in
290 * "An Implementation of the Generalized Basis Reduction Algorithm
291 * for Integer Programming" of Cook el al. to find an integer point
293 * If the polyhedron is unbounded, we first remove its rays.
295 Vector
*Polyhedron_Sample(Polyhedron
*P
, struct barvinok_options
*options
)
298 Vector
*sample
= NULL
, *obj
= NULL
;
310 if (P
->Dimension
== 0) {
311 sample
= Vector_Alloc(1);
312 value_set_si(sample
->p
[0], 1);
316 if (P
->Dimension
== 1)
317 POL_ENSURE_VERTICES(P
);
319 for (i
= 0; i
< P
->NbRays
; ++i
)
320 if (value_one_p(P
->Ray
[i
][1+P
->Dimension
])) {
321 sample
= Vector_Alloc(P
->Dimension
+1);
322 Vector_Copy(P
->Ray
[i
]+1, sample
->p
, P
->Dimension
+1);
328 Polyhedron
*Q
= remove_all_equalities(P
, &T
, options
->MaxRays
);
331 Q_sample
= Polyhedron_Sample(Q
, options
);
334 sample
= Vector_Alloc(P
->Dimension
+ 1);
335 Matrix_Vector_Product(T
, Q_sample
->p
, sample
->p
);
336 Vector_Free(Q_sample
);
342 Q
= Polyhedron_Factor(P
, 0, &T
, options
->MaxRays
);
344 sample
= product_sample(Q
, T
, options
);
353 obj
= Vector_Alloc(P
->Dimension
+1);
354 value_set_si(obj
->p
[0], 1);
355 res
= polyhedron_range(P
, obj
->p
, obj
->p
[0], &min
, &max
, options
);
356 if (res
== lp_unbounded
) {
362 Q
= remove_ray(P
, options
->MaxRays
);
364 sample
= Polyhedron_Sample(Q
, options
);
368 if (res
== lp_empty
) {
374 assert(res
== lp_ok
);
376 if (value_eq(min
, max
)) {
380 Matrix
*basis
= Polyhedron_Reduced_Basis(P
, options
);
384 T
= Matrix_Alloc(P
->Dimension
+1, P
->Dimension
+1);
385 inv
= Matrix_Alloc(P
->Dimension
+1, P
->Dimension
+1);
386 for (i
= 0; i
< P
->Dimension
; ++i
)
387 for (j
= 0; j
< P
->Dimension
; ++j
)
388 value_assign(T
->p
[i
][j
], basis
->p
[i
][j
]);
389 value_set_si(T
->p
[P
->Dimension
][P
->Dimension
], 1);
393 ok
= Matrix_Inverse(M
, inv
);
397 Q
= Polyhedron_Image(P
, T
, options
->MaxRays
);
398 res
= polyhedron_range(Q
, obj
->p
, obj
->p
[0], &min
, &max
, options
);
399 assert(res
== lp_ok
);
406 v
= Vector_Alloc(1+Q
->Dimension
+1);
407 value_set_si(v
->p
[1], -1);
409 for (value_assign(tmp
, min
); value_le(tmp
, max
); value_increment(tmp
, tmp
)) {
413 value_assign(v
->p
[1+Q
->Dimension
], tmp
);
415 R
= AddConstraints(v
->p
, 1, Q
, options
->MaxRays
);
416 R
= DomainConstraintSimplify(R
, options
->MaxRays
);
422 S
= Polyhedron_RemoveFixedColumns(R
, &T
);
424 S_sample
= Polyhedron_Sample(S
, options
);
427 Vector
*Q_sample
= obj
;
429 Matrix_Vector_Product(T
, S_sample
->p
, Q_sample
->p
);
431 Vector_Free(S_sample
);
435 sample
= Vector_Alloc(P
->Dimension
+ 1);
436 Matrix_Vector_Product(inv
, Q_sample
->p
, sample
->p
);
437 Vector_Free(Q_sample
);