semigroup_holes.cc

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