doc: add another paper refering to the library
[barvinok.git] / hilbert.c
blobf1516d9ed3f77d0a9c63b46251baa1f4c1fd58ce
1 #include <assert.h>
2 #include <stdlib.h>
3 #define Vector ZSolveVector
4 #define Matrix ZSolveMatrix
5 #include "zsolve/libzsolve.h"
6 #undef Vector
7 #undef Matrix
8 #include <barvinok/options.h>
9 #include <barvinok/util.h>
10 #include "hilbert.h"
11 #include "normalization.h"
12 #include "polysign.h"
14 static ZSolveMatrix Matrix2zsolve(Matrix *M)
16 int i, j;
17 ZSolveMatrix zmatrix;
19 zmatrix = createMatrix(M->NbColumns-2, M->NbRows);
20 for (i = 0; i < M->NbRows; ++i)
21 for (j = 0; j < M->NbColumns-2; ++j) {
22 assert(mpz_cmp_si(M->p[i][1+j], -MAXINT) > 0);
23 assert(mpz_cmp_si(M->p[i][1+j], MAXINT) < 0);
24 zmatrix->Data[i*zmatrix->Width+j] = mpz_get_si(M->p[i][1+j]);
27 return zmatrix;
30 static Matrix *VectorArray2Matrix(VectorArray array, unsigned cols)
32 int i, j;
33 Matrix *M = Matrix_Alloc(array->Size, cols+1);
35 for (i = 0; i < array->Size; ++i) {
36 for (j = 0; j < cols; ++j)
37 value_set_si(M->p[i][j], array->Data[i][j]);
38 value_set_si(M->p[i][cols], 1);
40 return M;
43 static void Polyhedron_Remove_Positivity_Constraint(Polyhedron *P)
45 int i;
47 for (i = 0; i < P->NbConstraints; ++i) {
48 if (First_Non_Zero(P->Constraint[i]+1, P->Dimension) != -1)
49 continue;
50 if (i < P->NbConstraints-1)
51 Vector_Exchange(P->Constraint[i],
52 P->Constraint[P->NbConstraints-1],
53 P->Dimension+2);
54 P->NbConstraints--;
55 --i;
59 /* Return
60 * T 0
61 * 0 1
63 static Matrix *homogenize(Matrix *T)
65 int i;
66 Matrix *H = Matrix_Alloc(T->NbRows+1, T->NbColumns+1);
68 for (i = 0; i < T->NbRows; ++i)
69 Vector_Copy(T->p[i], H->p[i], T->NbColumns);
70 value_set_si(H->p[T->NbRows][T->NbColumns], 1);
71 return H;
74 static Matrix *Polyhedron2standard_form(Polyhedron *P, Matrix **T)
76 int i, j;
77 int rows;
78 unsigned dim = P->Dimension;
79 Matrix *M2;
80 Matrix *H, *U;
81 Matrix M;
83 assert(P->NbEq == 0);
84 Polyhedron_Remove_Positivity_Constraint(P);
85 for (i = 0; i < P->NbConstraints; ++i)
86 assert(value_zero_p(P->Constraint[i][1+dim]));
88 Polyhedron_Matrix_View(P, &M, P->NbConstraints);
89 H = standard_constraints(&M, 0, &rows, &U);
90 *T = homogenize(U);
91 Matrix_Free(U);
93 M2 = Matrix_Alloc(rows, 2+dim+rows);
95 for (i = dim; i < H->NbRows; ++i) {
96 Vector_Copy(H->p[i], M2->p[i-dim]+1, dim);
97 value_set_si(M2->p[i-dim][1+i], -1);
99 for (i = 0, j = H->NbRows-dim; i < dim; ++i) {
100 if (First_Non_Zero(H->p[i], i) == -1)
101 continue;
102 Vector_Oppose(H->p[i], M2->p[j]+1, dim);
103 value_set_si(M2->p[j][1+j+dim], 1);
104 ++j;
106 Matrix_Free(H);
107 return M2;
110 /* Assumes C is a linear cone (i.e. with apex zero).
111 * All equalities are removed first to speed up the computation
112 * in zsolve.
114 Matrix *Cone_Hilbert_Basis(Polyhedron *C, unsigned MaxRays)
116 unsigned dim;
117 int i;
118 Matrix *M2, *M3, *T;
119 Matrix *CV = NULL;
120 LinearSystem initialsystem;
121 ZSolveMatrix matrix;
122 ZSolveVector rhs;
123 ZSolveContext ctx;
125 remove_all_equalities(&C, NULL, NULL, &CV, 0, MaxRays);
126 dim = C->Dimension;
128 for (i = 0; i < C->NbConstraints; ++i)
129 assert(value_zero_p(C->Constraint[i][1+dim]) ||
130 First_Non_Zero(C->Constraint[i]+1, dim) == -1);
132 M2 = Polyhedron2standard_form(C, &T);
133 matrix = Matrix2zsolve(M2);
134 Matrix_Free(M2);
136 rhs = createVector(matrix->Height);
137 for (i = 0; i < matrix->Height; i++)
138 rhs[i] = 0;
140 initialsystem = createLinearSystem();
141 setLinearSystemMatrix(initialsystem, matrix);
142 deleteMatrix(matrix);
144 setLinearSystemRHS(initialsystem, rhs);
145 deleteVector(rhs);
147 setLinearSystemLimit(initialsystem, -1, 0, MAXINT, 0);
148 setLinearSystemEquationType(initialsystem, -1, EQUATION_EQUAL, 0);
150 ctx = createZSolveContextFromSystem(initialsystem, NULL, 0, 0, NULL, NULL);
151 zsolveSystem(ctx, 0);
153 M2 = VectorArray2Matrix(ctx->Homs, C->Dimension);
154 deleteZSolveContext(ctx, 1);
155 Matrix_Transposition(T);
156 M3 = Matrix_Alloc(M2->NbRows, M2->NbColumns);
157 Matrix_Product(M2, T, M3);
158 Matrix_Free(M2);
159 Matrix_Free(T);
161 if (CV) {
162 Matrix *T, *M;
163 T = Transpose(CV);
164 M = Matrix_Alloc(M3->NbRows, T->NbColumns);
165 Matrix_Product(M3, T, M);
166 Matrix_Free(M3);
167 Matrix_Free(CV);
168 Matrix_Free(T);
169 Polyhedron_Free(C);
170 M3 = M;
173 return M3;
176 /* Assumes no two arrays of values are the same, so we can just
177 * look for the first different elements, without knowing the
178 * length of the arrays.
180 static int lex_cmp(const void *va, const void *vb)
182 int i;
183 const Value *a = *(const Value **)va;
184 const Value *b = *(const Value **)vb;
186 for (i = 0; ; ++i) {
187 int sign = mpz_cmp(a[i], b[i]);
188 if (sign)
189 return sign;
193 /* Compute integer hull of truncated linear cone C, i.e., of C with
194 * the origin removed.
195 * Here, we do this by first computing the Hilbert basis of C
196 * and then discarding elements from this basis that are rational
197 * overconvex combinations of other elements in the basis.
199 Matrix *Cone_Hilbert_Integer_Hull(Polyhedron *C,
200 struct barvinok_options *options)
202 int i, j, k;
203 Matrix *hilbert = Cone_Hilbert_Basis(C, options->MaxRays);
204 Matrix *rays, *hull;
205 unsigned dim = C->Dimension;
206 Value tmp;
207 unsigned MaxRays = options->MaxRays;
209 /* When checking for redundant points below, we want to
210 * check if there are any _rational_ solutions.
212 POL_UNSET(options->MaxRays, POL_INTEGER);
214 POL_ENSURE_VERTICES(C);
215 rays = Matrix_Alloc(C->NbRays-1, C->Dimension);
216 for (i = 0, j = 0; i < C->NbRays; ++i) {
217 if (value_notzero_p(C->Ray[i][1+C->Dimension]))
218 continue;
219 Vector_Copy(C->Ray[i]+1, rays->p[j++], C->Dimension);
222 /* We only sort the pointers into the big Value array */
223 qsort(rays->p, rays->NbRows, sizeof(Value *), lex_cmp);
224 qsort(hilbert->p, hilbert->NbRows, sizeof(Value *), lex_cmp);
226 /* Remove rays from Hilbert basis */
227 for (i = 0, j = 0, k = 0; i < hilbert->NbRows && j < rays->NbRows; ++i) {
228 if (Vector_Equal(hilbert->p[i], rays->p[j], C->Dimension))
229 ++j;
230 else
231 hilbert->p[k++] = hilbert->p[i];
233 hilbert->NbRows = k;
235 /* Now remove points that are overconvex combinations of other points */
236 value_init(tmp);
237 for (i = 0; hilbert->NbRows > 1 && i < hilbert->NbRows; ++i) {
238 Matrix *LP;
239 Vector *obj;
240 int nray = rays->NbRows;
241 int npoint = hilbert->NbRows;
242 enum lp_result result;
244 LP = Matrix_Alloc(dim + 1 + nray + (npoint-1), 2 + nray + (npoint-1));
245 for (j = 0; j < dim; ++j) {
246 for (k = 0; k < nray; ++k)
247 value_assign(LP->p[j][k+1], rays->p[k][j]);
248 for (k = 0; k < npoint; ++k) {
249 if (k == i)
250 value_oppose(LP->p[j][1+nray+npoint-1], hilbert->p[k][j]);
251 else
252 value_assign(LP->p[j][1+nray+k-(k>i)], hilbert->p[k][j]);
255 value_set_si(LP->p[dim][0], 1);
256 for (k = 0; k < nray+npoint-1; ++k)
257 value_set_si(LP->p[dim][1+k], 1);
258 value_set_si(LP->p[dim][LP->NbColumns-1], -1);
259 for (k = 0; k < LP->NbColumns-2; ++k) {
260 value_set_si(LP->p[dim+1+k][0], 1);
261 value_set_si(LP->p[dim+1+k][1+k], 1);
264 /* Somewhat arbitrary objective function. */
265 obj = Vector_Alloc(LP->NbColumns-1);
266 value_set_si(obj->p[0], 1);
267 value_set_si(obj->p[obj->Size-1], 1);
269 result = constraints_opt(LP, obj->p, obj->p[0], lp_min, &tmp,
270 options);
272 /* If the LP is not empty, the point can be discarded */
273 if (result != lp_empty) {
274 hilbert->NbRows--;
275 if (i < hilbert->NbRows)
276 hilbert->p[i] = hilbert->p[hilbert->NbRows];
277 --i;
280 Matrix_Free(LP);
281 Vector_Free(obj);
283 value_clear(tmp);
285 hull = Matrix_Alloc(rays->NbRows + hilbert->NbRows, dim+1);
286 for (i = 0; i < rays->NbRows; ++i) {
287 Vector_Copy(rays->p[i], hull->p[i], dim);
288 value_set_si(hull->p[i][dim], 1);
290 for (i = 0; i < hilbert->NbRows; ++i) {
291 Vector_Copy(hilbert->p[i], hull->p[rays->NbRows+i], dim);
292 value_set_si(hull->p[rays->NbRows+i][dim], 1);
294 Matrix_Free(rays);
295 Matrix_Free(hilbert);
297 options->MaxRays = MaxRays;
298 return hull;