semigroup_holes.cc

   1 #include <assert.h>
   2 #include <iostream>
   3 #include <barvinok/barvinok.h>
   4 #include <barvinok/genfun.h>
   5
   6 /* This program computes the generating function of the holes in a semigroup.
   7  * The input is a matrix with the generators m_i of the semigroup as columns.
   8  * The output is the generating function for the set
   9  *
  10  *      (L \cap \cone { m_i }_i) \setminus S,
  11  *
  12  * with
  13  *      S = { \sum_i \lambda_i m_i | \lambda_i \in \Z_+ }
  14  * and
  15  *      L = { \sum_i \lambda_i m_i | \lambda_i \in \Z }
  16  */
  17
  18 using std::cout;
  19 using std::cerr;
  20 using std::endl;
  21
  22 /* Compute cone { m_i }_i, with m_i the columsn of generators.
  23  */
  24 static Polyhedron *compute_cone(Matrix *generators, barvinok_options *options)
  25 {
  26         Matrix *M;
  27         Polyhedron *cone;
  28
  29         M = Matrix_Alloc(generators->NbColumns + 1, 1 + generators->NbRows + 1);
  30         assert(M);
  31         value_set_si(M->p[0][0], 1);
  32         value_set_si(M->p[0][M->NbColumns - 1], 1);
  33         for (int i = 0; i < generators->NbColumns; ++i) {
  34                 value_set_si(M->p[1 + i][0], 1);
  35                 for (int j = 0; j < generators->NbRows; ++j)
  36                         value_assign(M->p[1 + i][1 + j], generators->p[j][i]);
  37         }
  38         cone = Rays2Polyhedron(M, options->MaxRays);
  39         Matrix_Free(M);
  40         return cone;
  41 }
  42
  43 /* Compute the generating function of L \cap \cone { m_i }_i, with
  44  * m_i the columns of generators and
  45  *
  46  *      L = { \sum_i \lambda_i m_i | \lambda_i \in \Z }
  47  *
  48  * We first compute the cone C = \cone { m_i }_i = { x | A x + b >= 0 }.
  49  * and we compute a minimal set of generators of the lattice by
  50  * computing the Hermite normal form H of generators.
  51  * Finally, we compute the generating function for the set
  52  *
  53  *      { x | \exists y : A x + b >= 0, x = H y }
  54  *
  55  * Since the number of columns in H is less that or equal to the number
  56  * of rows, no projection needs to be performed during the computation
  57  * of this generating function.
  58  */
  59 static gen_fun *compute_lattice_gf(Matrix *generators,
  60         barvinok_options *options)
  61 {
  62         Matrix *H;
  63         Matrix *M;
  64         Polyhedron *cone, *L;
  65         int col, row;
  66         gen_fun *gf;
  67
  68         cone = compute_cone(generators, options);
  69
  70         left_hermite(generators, &H, NULL, NULL);
  71
  72         for (row = col = 0; col < H->NbColumns; ++col) {
  73                 while (row < H->NbRows && value_zero_p(H->p[row][col]))
  74                         ++row;
  75                 if (row >= H->NbRows)
  76                         break;
  77         }
  78
  79         M = Matrix_Alloc(cone->NbConstraints + cone->Dimension,
  80                          1 + col + cone->Dimension + 1);
  81         assert(M);
  82         for (int i = 0; i < cone->NbConstraints; ++i) {
  83                 value_assign(M->p[i][0], cone->Constraint[i][0]);
  84                 for (int j = 0; j < cone->Dimension + 1; ++j)
  85                         value_assign(M->p[i][1 + col + j],
  86                                         cone->Constraint[i][1 + j]);
  87         }
  88         for (int i = 0; i < cone->Dimension; ++i) {
  89                 int row = cone->NbConstraints + i;
  90                 value_set_si(M->p[row][1 + col + i], -1);
  91                 for (int j = 0; j < col; ++j)
  92                         value_assign(M->p[row][1 + j], H->p[i][j]);
  93         }
  94         Matrix_Free(H);
  95         Polyhedron_Free(cone);
  96         L = Constraints2Polyhedron(M, options->MaxRays);
  97         Matrix_Free(M);
  98         gf = barvinok_enumerate_e_series(L, col, generators->NbRows, options);
  99         Polyhedron_Free(L);
 100         return gf;
 101 }
 102
 103 /* Compute the generating function of
 104  *
 105  *      S = { x | \exists \lambda_i: x = \sum_i \lambda_i m_i, \lambda_i >= 0 }
 106  *
 107  * with m_i the columns of generators.
 108  */
 109 static gen_fun *compute_semigroup_gf(Matrix *generators,
 110         barvinok_options *options)
 111 {
 112         Matrix *M;
 113         Polyhedron *S;
 114         gen_fun *gf;
 115
 116         M = Matrix_Alloc(generators->NbRows + generators->NbColumns,
 117                          1 + generators->NbColumns + generators->NbRows + 1);
 118         assert(M);
 119         for (int i = 0; i < generators->NbRows; ++i) {
 120                 value_set_si(M->p[i][1 + generators->NbColumns + i], -1);
 121                 for (int j = 0; j < generators->NbColumns; ++j)
 122                         value_assign(M->p[i][1 + j], generators->p[i][j]);
 123         }
 124         for (int i = 0; i < generators->NbColumns; ++i) {
 125                 int row = generators->NbRows + i;
 126                 value_set_si(M->p[row][0], 1);
 127                 value_set_si(M->p[row][1 + i], 1);
 128         }
 129         S = Constraints2Polyhedron(M, options->MaxRays);
 130         Matrix_Free(M);
 131         gf = barvinok_enumerate_e_series(S, generators->NbColumns,
 132                                                 generators->NbRows, options);
 133         Polyhedron_Free(S);
 134         return gf;
 135 }
 136
 137 /* Compute the generating function of
 138  *
 139  *      (L \cap \cone { m_i }_i) \setminus S,
 140  *
 141  * Since S \subset (L \cap \cone { m_i }_i), this is simply computed as
 142  *
 143  *      f_{(L \cap \cone { m_i }_i)} - f_S
 144  */
 145 static gen_fun *compute_holes_gf(Matrix *generators, barvinok_options *options)
 146 {
 147         gen_fun *lattice;
 148         gen_fun *semigroup;
 149         gen_fun *holes;
 150         QQ mone(-1, 1);
 151
 152         lattice = compute_lattice_gf(generators, options);
 153
 154         semigroup = compute_semigroup_gf(generators, options);
 155
 156         holes = lattice;
 157         holes->add(mone, semigroup, options);
 158         delete semigroup;
 159
 160         return holes;
 161 }
 162
 163 int main(int argc, char **argv)
 164 {
 165         Matrix *generators;
 166         gen_fun *holes;
 167         barvinok_options *options = barvinok_options_new_with_defaults();
 168
 169         argc = barvinok_options_parse(options, argc, argv, ISL_ARG_ALL);
 170
 171         generators = Matrix_Read();
 172         assert(generators);
 173         holes = compute_holes_gf(generators, options);
 174         Matrix_Free(generators);
 175
 176         holes->print(cout, 0, NULL);
 177         cout << endl;
 178
 179         barvinok_options_free(options);
 180         return 0;
 181 }