update pet for support for recent clangs
[barvinok.git] / semigroup_holes.cc
blob2662363bd569430d0e0764c7a2832f76ecd69ff3
1 #include <assert.h>
2 #include <iostream>
3 #include <barvinok/barvinok.h>
4 #include <barvinok/genfun.h>
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
10 * (L \cap \cone { m_i }_i) \setminus S,
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 }
18 using std::cout;
19 using std::cerr;
20 using std::endl;
22 /* Compute cone { m_i }_i, with m_i the columsn of generators.
24 static Polyhedron *compute_cone(Matrix *generators, barvinok_options *options)
26 Matrix *M;
27 Polyhedron *cone;
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]);
38 cone = Rays2Polyhedron(M, options->MaxRays);
39 Matrix_Free(M);
40 return cone;
43 /* Compute the generating function of L \cap \cone { m_i }_i, with
44 * m_i the columns of generators and
46 * L = { \sum_i \lambda_i m_i | \lambda_i \in \Z }
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
53 * { x | \exists y : A x + b >= 0, x = H y }
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.
59 static gen_fun *compute_lattice_gf(Matrix *generators,
60 barvinok_options *options)
62 Matrix *H;
63 Matrix *M;
64 Polyhedron *cone, *L;
65 int col, row;
66 gen_fun *gf;
68 cone = compute_cone(generators, options);
70 left_hermite(generators, &H, NULL, NULL);
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;
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]);
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]);
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;
103 /* Compute the generating function of
105 * S = { x | \exists \lambda_i: x = \sum_i \lambda_i m_i, \lambda_i >= 0 }
107 * with m_i the columns of generators.
109 static gen_fun *compute_semigroup_gf(Matrix *generators,
110 barvinok_options *options)
112 Matrix *M;
113 Polyhedron *S;
114 gen_fun *gf;
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]);
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);
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;
137 /* Compute the generating function of
139 * (L \cap \cone { m_i }_i) \setminus S,
141 * Since S \subset (L \cap \cone { m_i }_i), this is simply computed as
143 * f_{(L \cap \cone { m_i }_i)} - f_S
145 static gen_fun *compute_holes_gf(Matrix *generators, barvinok_options *options)
147 gen_fun *lattice;
148 gen_fun *semigroup;
149 gen_fun *holes;
150 QQ mone(-1, 1);
152 lattice = compute_lattice_gf(generators, options);
154 semigroup = compute_semigroup_gf(generators, options);
156 holes = lattice;
157 holes->add(mone, semigroup, options);
158 delete semigroup;
160 return holes;
163 int main(int argc, char **argv)
165 Matrix *generators;
166 gen_fun *holes;
167 barvinok_options *options = barvinok_options_new_with_defaults();
169 argc = barvinok_options_parse(options, argc, argv, ISL_ARG_ALL);
171 generators = Matrix_Read();
172 assert(generators);
173 holes = compute_holes_gf(generators, options);
174 Matrix_Free(generators);
176 holes->print(cout, 0, NULL);
177 cout << endl;
179 barvinok_options_free(options);
180 return 0;