hilbert.c

   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"
  13
  14 static ZSolveMatrix Matrix2zsolve(Matrix *M)
  15 {
  16     int i, j;
  17     ZSolveMatrix zmatrix;
  18
  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]);
  25         }
  26
  27     return zmatrix;
  28 }
  29
  30 static Matrix *VectorArray2Matrix(VectorArray array, unsigned cols)
  31 {
  32     int i, j;
  33     Matrix *M = Matrix_Alloc(array->Size, cols+1);
  34
  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);
  39     }
  40     return M;
  41 }
  42
  43 static void Polyhedron_Remove_Positivity_Constraint(Polyhedron *P)
  44 {
  45     int i;
  46
  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;
  56     }
  57 }
  58
  59 /* Return
  60  *             T 0
  61  *             0 1
  62  */
  63 static Matrix *homogenize(Matrix *T)
  64 {
  65     int i;
  66     Matrix *H = Matrix_Alloc(T->NbRows+1, T->NbColumns+1);
  67
  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;
  72 }
  73
  74 static Matrix *Polyhedron2standard_form(Polyhedron *P, Matrix **T)
  75 {
  76     int i, j;
  77     int rows;
  78     unsigned dim = P->Dimension;
  79     Matrix *M2;
  80     Matrix *H, *U;
  81     Matrix M;
  82
  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]));
  87
  88     Polyhedron_Matrix_View(P, &M, P->NbConstraints);
  89     H = standard_constraints(&M, 0, &rows, &U);
  90     *T = homogenize(U);
  91     Matrix_Free(U);
  92
  93     M2 = Matrix_Alloc(rows, 2+dim+rows);
  94
  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);
  98     }
  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;
 105     }
 106     Matrix_Free(H);
 107     return M2;
 108 }
 109
 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.
 113  */
 114 Matrix *Cone_Hilbert_Basis(Polyhedron *C, unsigned MaxRays)
 115 {
 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;
 124
 125     remove_all_equalities(&C, NULL, NULL, &CV, 0, MaxRays);
 126     dim = C->Dimension;
 127
 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);
 131
 132     M2 = Polyhedron2standard_form(C, &T);
 133     matrix = Matrix2zsolve(M2);
 134     Matrix_Free(M2);
 135
 136     rhs = createVector(matrix->Height);
 137     for (i = 0; i < matrix->Height; i++)
 138         rhs[i] = 0;
 139
 140     initialsystem = createLinearSystem();
 141     setLinearSystemMatrix(initialsystem, matrix);
 142     deleteMatrix(matrix);
 143
 144     setLinearSystemRHS(initialsystem, rhs);
 145     deleteVector(rhs);
 146
 147     setLinearSystemLimit(initialsystem, -1, 0, MAXINT, 0);
 148     setLinearSystemEquationType(initialsystem, -1, EQUATION_EQUAL, 0);
 149
 150     ctx = createZSolveContextFromSystem(initialsystem, NULL, 0, 0, NULL, NULL);
 151     zsolveSystem(ctx, 0);
 152
 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);
 160
 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;
 171     }
 172
 173     return M3;
 174 }
 175
 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.
 179  */
 180 static int lex_cmp(const void *va, const void *vb)
 181 {
 182     int i;
 183     const Value *a = *(const Value **)va;
 184     const Value *b = *(const Value **)vb;
 185
 186     for (i = 0; ; ++i) {
 187         int sign = mpz_cmp(a[i], b[i]);
 188         if (sign)
 189             return sign;
 190     }
 191 }
 192
 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.
 198  */
 199 Matrix *Cone_Hilbert_Integer_Hull(Polyhedron *C,
 200                                   struct barvinok_options *options)
 201 {
 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;
 208
 209     /* When checking for redundant points below, we want to
 210      * check if there are any _rational_ solutions.
 211      */
 212     POL_UNSET(options->MaxRays, POL_INTEGER);
 213
 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);
 220     }
 221
 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);
 225
 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];
 232     }
 233     hilbert->NbRows = k;
 234
 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;
 243
 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]);
 253             }
 254         }
 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);
 262         }
 263
 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);
 268
 269         result = constraints_opt(LP, obj->p, obj->p[0], lp_min, &tmp,
 270                                  options);
 271
 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;
 278         }
 279
 280         Matrix_Free(LP);
 281         Vector_Free(obj);
 282     }
 283     value_clear(tmp);
 284
 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);
 289     }
 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);
 293     }
 294     Matrix_Free(rays);
 295     Matrix_Free(hilbert);
 296
 297     options->MaxRays = MaxRays;
 298     return hull;
 299 }