MLIR  19.0.0git
Barvinok.cpp
Go to the documentation of this file.
1 //===- Barvinok.cpp - Barvinok's Algorithm ----------------------*- C++ -*-===//
2 //
3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6 //
7 //===----------------------------------------------------------------------===//
8
12 #include <algorithm>
13 #include <bitset>
14
15 using namespace mlir;
16 using namespace presburger;
17 using namespace mlir::presburger::detail;
18
19 /// Assuming that the input cone is pointed at the origin,
20 /// converts it to its dual in V-representation.
21 /// Essentially we just remove the all-zeroes constant column.
23  unsigned numIneq = cone.getNumInequalities();
24  unsigned numVar = cone.getNumCols() - 1;
25  ConeV dual(numIneq, numVar, 0, 0);
26  // Assuming that an inequality of the form
27  // a1*x1 + ... + an*xn + b ≥ 0
28  // is represented as a row [a1, ..., an, b]
29  // and that b = 0.
30
31  for (auto i : llvm::seq<int>(0, numIneq)) {
32  assert(cone.atIneq(i, numVar) == 0 &&
33  "H-representation of cone is not centred at the origin!");
34  for (unsigned j = 0; j < numVar; ++j) {
35  dual.at(i, j) = cone.atIneq(i, j);
36  }
37  }
38
39  // Now dual is of the form [ [a1, ..., an] , ... ]
40  // which is the V-representation of the dual.
41  return dual;
42 }
43
44 /// Converts a cone in V-representation to the H-representation
45 /// of its dual, pointed at the origin (not at the original vertex).
46 /// Essentially adds a column consisting only of zeroes to the end.
48  unsigned rows = cone.getNumRows();
49  unsigned columns = cone.getNumColumns();
50  ConeH dual = defineHRep(columns);
51  // Add a new column (for constants) at the end.
52  // This will be initialized to zero.
53  cone.insertColumn(columns);
54
55  for (unsigned i = 0; i < rows; ++i)
57
58  // Now dual is of the form [ [a1, ..., an, 0] , ... ]
59  // which is the H-representation of the dual.
60  return dual;
61 }
62
63 /// Find the index of a cone in V-representation.
65  if (cone.getNumRows() > cone.getNumColumns())
66  return MPInt(0);
67
68  return cone.determinant();
69 }
70
71 /// Compute the generating function for a unimodular cone.
72 /// This consists of a single term of the form
73 /// sign * x^num / prod_j (1 - x^den_j)
74 ///
75 /// sign is either +1 or -1.
76 /// den_j is defined as the set of generators of the cone.
77 /// num is computed by expressing the vertex as a weighted
78 /// sum of the generators, and then taking the floor of the
79 /// coefficients.
82  ParamPoint vertex, int sign, const ConeH &cone) {
83  // Consider a cone with H-representation [0 -1].
84  // [-1 -2]
85  // Let the vertex be given by the matrix [ 2 2 0], with 2 params.
86  // [-1 -1/2 1]
87
88  // cone must be unimodular.
89  assert(abs(getIndex(getDual(cone))) == 1 && "input cone is not unimodular!");
90
91  unsigned numVar = cone.getNumVars();
92  unsigned numIneq = cone.getNumInequalities();
93
94  // Thus its ray matrix, U, is the inverse of the
95  // transpose of its inequality matrix, cone.
96  // The last column of the inequality matrix is null,
97  // so we remove it to obtain a square matrix.
99  transp.removeRow(numVar);
100
101  FracMatrix generators(numVar, numIneq);
102  transp.determinant(/*inverse=*/&generators); // This is the U-matrix.
103  // Thus the generators are given by U = [2 -1].
104  // [-1 0]
105
106  // The powers in the denominator of the generating
107  // function are given by the generators of the cone,
108  // i.e., the rows of the matrix U.
109  std::vector<Point> denominator(numIneq);
110  ArrayRef<Fraction> row;
111  for (auto i : llvm::seq<int>(0, numVar)) {
112  row = generators.getRow(i);
113  denominator[i] = Point(row);
114  }
115
116  // The vertex is v \in Z^{d x (n+1)}
117  // We need to find affine functions of parameters λ_i(p)
118  // such that v = Σ λ_i(p)*u_i,
119  // where u_i are the rows of U (generators)
120  // The λ_i are given by the columns of Λ = v^T U^{-1}, and
121  // we have transp = U^{-1}.
122  // Then the exponent in the numerator will be
123  // Σ -floor(-λ_i(p))*u_i.
124  // Thus we store the (exponent of the) numerator as the affine function -Λ,
125  // since the generators u_i are already stored as the exponent of the
126  // denominator. Note that the outer -1 will have to be accounted for, as it is
127  // not stored. See end for an example.
128
129  unsigned numColumns = vertex.getNumColumns();
130  unsigned numRows = vertex.getNumRows();
131  ParamPoint numerator(numColumns, numRows);
132  SmallVector<Fraction> ithCol(numRows);
133  for (auto i : llvm::seq<int>(0, numColumns)) {
134  for (auto j : llvm::seq<int>(0, numRows))
135  ithCol[j] = vertex(j, i);
136  numerator.setRow(i, transp.preMultiplyWithRow(ithCol));
137  numerator.negateRow(i);
138  }
139  // Therefore Λ will be given by [ 1 0 ] and the negation of this will be
140  // [ 1/2 -1 ]
141  // [ -1 -2 ]
142  // stored as the numerator.
143  // Algebraically, the numerator exponent is
144  // [ -2 ⌊ - N - M/2 + 1 ⌋ + 1 ⌊ 0 + M + 2 ⌋ ] -> first COLUMN of U is [2, -1]
145  // [ 1 ⌊ - N - M/2 + 1 ⌋ + 0 ⌊ 0 + M + 2 ⌋ ] -> second COLUMN of U is [-1, 0]
146
147  return GeneratingFunction(numColumns - 1, SmallVector<int>(1, sign),
148  std::vector({numerator}),
149  std::vector({denominator}));
150 }
151
152 /// We use Gaussian elimination to find the solution to a set of d equations
153 /// of the form
154 /// a_1 x_1 + ... + a_d x_d + b_1 m_1 + ... + b_p m_p + c = 0
155 /// where x_i are variables,
156 /// m_i are parameters and
157 /// a_i, b_i, c are rational coefficients.
158 ///
159 /// The solution expresses each x_i as an affine function of the m_i, and is
160 /// therefore represented as a matrix of size d x (p+1).
161 /// If there is no solution, we return null.
162 std::optional<ParamPoint>
164  // equations is a d x (d + p + 1) matrix.
165  // Each row represents an equation.
166  unsigned d = equations.getNumRows();
167  unsigned numCols = equations.getNumColumns();
168
169  // If the determinant is zero, there is no unique solution.
170  // Thus we return null.
171  if (FracMatrix(equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/d - 1,
172  /*fromColumn=*/0,
173  /*toColumn=*/d - 1))
174  .determinant() == 0)
175  return std::nullopt;
176
177  // Perform row operations to make each column all zeros except for the
178  // diagonal element, which is made to be one.
179  for (unsigned i = 0; i < d; ++i) {
180  // First ensure that the diagonal element is nonzero, by swapping
181  // it with a row that is non-zero at column i.
182  if (equations(i, i) != 0)
183  continue;
184  for (unsigned j = i + 1; j < d; ++j) {
185  if (equations(j, i) == 0)
186  continue;
187  equations.swapRows(j, i);
188  break;
189  }
190
191  Fraction diagElement = equations(i, i);
192
193  // Apply row operations to make all elements except the diagonal to zero.
194  for (unsigned j = 0; j < d; ++j) {
195  if (i == j)
196  continue;
197  if (equations(j, i) == 0)
198  continue;
199  // Apply row operations to make element (j, i) zero by subtracting the
200  // ith row, appropriately scaled.
201  Fraction currentElement = equations(j, i);
203  /*scale=*/-currentElement / diagElement);
204  }
205  }
206
207  // Rescale diagonal elements to 1.
208  for (unsigned i = 0; i < d; ++i)
209  equations.scaleRow(i, 1 / equations(i, i));
210
211  // Now we have reduced the equations to the form
212  // x_i + b_1' m_1 + ... + b_p' m_p + c' = 0
213  // i.e. each variable appears exactly once in the system, and has coefficient
214  // one.
215  //
216  // Thus we have
217  // x_i = - b_1' m_1 - ... - b_p' m_p - c
218  // and so we return the negation of the last p + 1 columns of the matrix.
219  //
220  // We copy these columns and return them.
221  ParamPoint vertex =
222  equations.getSubMatrix(/*fromRow=*/0, /*toRow=*/d - 1,
223  /*fromColumn=*/d, /*toColumn=*/numCols - 1);
224  vertex.negateMatrix();
225  return vertex;
226 }
227
228 /// This is an implementation of the Clauss-Loechner algorithm for chamber
229 /// decomposition.
230 ///
231 /// We maintain a list of pairwise disjoint chambers and the generating
232 /// functions corresponding to each one. We iterate over the list of regions,
233 /// each time adding the current region's generating function to the chambers
234 /// where it is active and separating the chambers where it is not.
235 ///
236 /// Given the region each generating function is active in, for each subset of
237 /// generating functions the region that (the sum of) precisely this subset is
238 /// in, is the intersection of the regions that these are active in,
239 /// intersected with the complements of the remaining regions.
240 std::vector<std::pair<PresburgerSet, GeneratingFunction>>
242  unsigned numSymbols, ArrayRef<std::pair<PresburgerSet, GeneratingFunction>>
243  regionsAndGeneratingFunctions) {
244  assert(!regionsAndGeneratingFunctions.empty() &&
245  "there must be at least one chamber!");
246  // We maintain a list of regions and their associated generating function
247  // initialized with the universe and the empty generating function.
248  std::vector<std::pair<PresburgerSet, GeneratingFunction>> chambers = {
250  GeneratingFunction(numSymbols, {}, {}, {})}};
251
252  // We iterate over the region list.
253  //
254  // For each activity region R_j (corresponding to the generating function
255  // gf_j), we examine all the current chambers R_i.
256  //
257  // If R_j has a full-dimensional intersection with an existing chamber R_i,
258  // then that chamber is replaced by two new ones:
259  // 1. the intersection R_i \cap R_j, where the generating function is
260  // gf_i + gf_j.
261  // 2. the difference R_i - R_j, where the generating function is gf_i.
262  //
263  // At each step, we define a new chamber list after considering gf_j,
264  // replacing and appending chambers as discussed above.
265  //
266  // The loop has the invariant that the union over all the chambers gives the
267  // universe at every step.
268  for (const auto &[region, generatingFunction] :
269  regionsAndGeneratingFunctions) {
270  std::vector<std::pair<PresburgerSet, GeneratingFunction>> newChambers;
271
272  for (const auto &[currentRegion, currentGeneratingFunction] : chambers) {
273  PresburgerSet intersection = currentRegion.intersect(region);
274
275  // If the intersection is not full-dimensional, we do not modify
276  // the chamber list.
277  if (!intersection.isFullDim()) {
278  newChambers.emplace_back(currentRegion, currentGeneratingFunction);
279  continue;
280  }
281
282  // If it is, we add the intersection and the difference as chambers.
283  newChambers.emplace_back(intersection,
284  currentGeneratingFunction + generatingFunction);
285  newChambers.emplace_back(currentRegion.subtract(region),
286  currentGeneratingFunction);
287  }
288  chambers = std::move(newChambers);
289  }
290
291  return chambers;
292 }
293
294 /// For a polytope expressed as a set of n inequalities, compute the generating
295 /// function corresponding to the lattice points included in the polytope. This
296 /// algorithm has three main steps:
297 /// 1. Enumerate the vertices, by iterating over subsets of inequalities and
298 /// checking for satisfiability. For each d-subset of inequalities (where d
299 /// is the number of variables), we solve to obtain the vertex in terms of
300 /// the parameters, and then check for the region in parameter space where
301 /// this vertex satisfies the remaining (n - d) inequalities.
302 /// 2. For each vertex, identify the tangent cone and compute the generating
303 /// function corresponding to it. The generating function depends on the
304 /// parametric expression of the vertex and the (non-parametric) generators
305 /// of the tangent cone.
306 /// 3. [Clauss-Loechner decomposition] Identify the regions in parameter space
307 /// (chambers) where each vertex is active, and accordingly compute the
308 /// GF of the polytope in each chamber.
309 ///
310 /// Verdoolaege, Sven, et al. "Counting integer points in parametric
311 /// polytopes using Barvinok's rational functions." Algorithmica 48 (2007):
312 /// 37-66.
313 std::vector<std::pair<PresburgerSet, GeneratingFunction>>
315  const PolyhedronH &poly) {
316  unsigned numVars = poly.getNumRangeVars();
317  unsigned numSymbols = poly.getNumSymbolVars();
318  unsigned numIneqs = poly.getNumInequalities();
319
320  // We store a list of the computed vertices.
321  std::vector<ParamPoint> vertices;
322  // For each vertex, we store the corresponding active region and the
323  // generating functions of the tangent cone, in order.
324  std::vector<std::pair<PresburgerSet, GeneratingFunction>>
325  regionsAndGeneratingFunctions;
326
327  // We iterate over all subsets of inequalities with cardinality numVars,
328  // using permutations of numVars 1's and (numIneqs - numVars) 0's.
329  //
330  // For a given permutation, we consider a subset which contains
331  // the i'th inequality if the i'th bit in the bitset is 1.
332  //
333  // We start with the permutation that takes the last numVars inequalities.
334  SmallVector<int> indicator(numIneqs);
335  for (unsigned i = numIneqs - numVars; i < numIneqs; ++i)
336  indicator[i] = 1;
337
338  do {
339  // Collect the inequalities corresponding to the bits which are set
340  // and the remaining ones.
341  auto [subset, remainder] = poly.getInequalities().splitByBitset(indicator);
342  // All other inequalities are stored in a2 and b2c2.
343  //
344  // These are column-wise splits of the inequalities;
345  // a2 stores the coefficients of the variables, and
346  // b2c2 stores the coefficients of the parameters and the constant term.
347  FracMatrix a2(numIneqs - numVars, numVars);
348  FracMatrix b2c2(numIneqs - numVars, numSymbols + 1);
349  a2 = FracMatrix(
350  remainder.getSubMatrix(0, numIneqs - numVars - 1, 0, numVars - 1));
351  b2c2 = FracMatrix(remainder.getSubMatrix(0, numIneqs - numVars - 1, numVars,
352  numVars + numSymbols));
353
354  // Find the vertex, if any, corresponding to the current subset of
355  // inequalities.
356  std::optional<ParamPoint> vertex =
357  solveParametricEquations(FracMatrix(subset)); // d x (p+1)
358
359  if (!vertex)
360  continue;
361  if (std::find(vertices.begin(), vertices.end(), vertex) != vertices.end())
362  continue;
363  // If this subset corresponds to a vertex that has not been considered,
364  // store it.
365  vertices.push_back(*vertex);
366
367  // If a vertex is formed by the intersection of more than d facets, we
368  // assume that any d-subset of these facets can be solved to obtain its
369  // expression. This assumption is valid because, if the vertex has two
370  // distinct parametric expressions, then a nontrivial equality among the
371  // parameters holds, which is a contradiction as we know the parameter
372  // space to be full-dimensional.
373
374  // Let the current vertex be [X | y], where
375  // X represents the coefficients of the parameters and
376  // y represents the constant term.
377  //
378  // The region (in parameter space) where this vertex is active is given
379  // by substituting the vertex into the *remaining* inequalities of the
380  // polytope (those which were not collected into subset), i.e., into the
381  // inequalities [A2 | B2 | c2].
382  //
383  // Thus, the coefficients of the parameters after substitution become
384  // (A2 • X + B2)
385  // and the constant terms become
386  // (A2 • y + c2).
387  //
388  // The region is therefore given by
389  // (A2 • X + B2) p + (A2 • y + c2) ≥ 0
390  //
391  // This is equivalent to A2 • [X | y] + [B2 | c2].
392  //
393  // Thus we premultiply [X | y] with each row of A2
394  // and add each row of [B2 | c2].
395  FracMatrix activeRegion(numIneqs - numVars, numSymbols + 1);
396  for (unsigned i = 0; i < numIneqs - numVars; i++) {
397  activeRegion.setRow(i, vertex->preMultiplyWithRow(a2.getRow(i)));
399  }
400
401  // We convert the representation of the active region to an integers-only
402  // form so as to store it as a PresburgerSet.
403  IntegerPolyhedron activeRegionRel(
404  PresburgerSpace::getRelationSpace(0, numSymbols, 0, 0), activeRegion);
405
406  // Now, we compute the generating function at this vertex.
407  // We collect the inequalities corresponding to each vertex to compute
408  // the tangent cone at that vertex.
409  //
410  // We only need the coefficients of the variables (NOT the parameters)
411  // as the generating function only depends on these.
412  // We translate the cones to be pointed at the origin by making the
413  // constant terms zero.
414  ConeH tangentCone = defineHRep(numVars);
415  for (unsigned j = 0, e = subset.getNumRows(); j < e; ++j) {
416  SmallVector<MPInt> ineq(numVars + 1);
417  for (unsigned k = 0; k < numVars; ++k)
418  ineq[k] = subset(j, k);
420  }
421  // We assume that the tangent cone is unimodular, so there is no need
422  // to decompose it.
423  //
424  // In the general case, the unimodular decomposition may have several
425  // cones.
426  GeneratingFunction vertexGf(numSymbols, {}, {}, {});
427  SmallVector<std::pair<int, ConeH>, 4> unimodCones = {{1, tangentCone}};
428  for (const std::pair<int, ConeH> &signedCone : unimodCones) {
429  auto [sign, cone] = signedCone;
430  vertexGf = vertexGf +
431  computeUnimodularConeGeneratingFunction(*vertex, sign, cone);
432  }
433  // We store the vertex we computed with the generating function of its
434  // tangent cone.
435  regionsAndGeneratingFunctions.emplace_back(PresburgerSet(activeRegionRel),
436  vertexGf);
437  } while (std::next_permutation(indicator.begin(), indicator.end()));
438
439  // Now, we use Clauss-Loechner decomposition to identify regions in parameter
440  // space where each vertex is active. These regions (chambers) have the
441  // property that no two of them have a full-dimensional intersection, i.e.,
442  // they may share "facets" or "edges", but their intersection can only have
443  // up to numVars - 1 dimensions.
444  //
445  // In each chamber, we sum up the generating functions of the active vertices
446  // to find the generating function of the polytope.
447  return computeChamberDecomposition(numSymbols, regionsAndGeneratingFunctions);
448 }
449
450 /// We use an iterative procedure to find a vector not orthogonal
451 /// to a given set, ignoring the null vectors.
452 /// Let the inputs be {x_1, ..., x_k}, all vectors of length n.
453 ///
454 /// In the following,
455 /// vs[:i] means the elements of vs up to and including the i'th one,
456 /// <vs, us> means the dot product of vs and us,
457 /// vs ++ [v] means the vector vs with the new element v appended to it.
458 ///
459 /// We proceed iteratively; for steps d = 0, ... n-1, we construct a vector
460 /// which is not orthogonal to any of {x_1[:d], ..., x_n[:d]}, ignoring
461 /// the null vectors.
462 /// At step d = 0, we let vs = [1]. Clearly this is not orthogonal to
463 /// any vector in the set {x_1[0], ..., x_n[0]}, except the null ones,
464 /// which we ignore.
465 /// At step d > 0 , we need a number v
466 /// s.t. <x_i[:d], vs++[v]> != 0 for all i.
467 /// => <x_i[:d-1], vs> + x_i[d]*v != 0
468 /// => v != - <x_i[:d-1], vs> / x_i[d]
469 /// We compute this value for all x_i, and then
470 /// set v to be the maximum element of this set plus one. Thus
471 /// v is outside the set as desired, and we append it to vs
472 /// to obtain the result of the d'th step.
474  ArrayRef<Point> vectors) {
475  unsigned dim = vectors[0].size();
476  assert(
477  llvm::all_of(vectors,
478  [&](const Point &vector) { return vector.size() == dim; }) &&
479  "all vectors need to be the same size!");
480
481  SmallVector<Fraction> newPoint = {Fraction(1, 1)};
482  Fraction maxDisallowedValue = -Fraction(1, 0),
483  disallowedValue = Fraction(0, 1);
484
485  for (unsigned d = 1; d < dim; ++d) {
486  // Compute the disallowed values - <x_i[:d-1], vs> / x_i[d] for each i.
487  maxDisallowedValue = -Fraction(1, 0);
488  for (const Point &vector : vectors) {
489  if (vector[d] == 0)
490  continue;
491  disallowedValue =
492  -dotProduct(ArrayRef(vector).slice(0, d), newPoint) / vector[d];
493
494  // Find the biggest such value
495  maxDisallowedValue = std::max(maxDisallowedValue, disallowedValue);
496  }
497  newPoint.push_back(maxDisallowedValue + 1);
498  }
499  return newPoint;
500 }
501
502 /// We use the following recursive formula to find the coefficient of
503 /// s^power in the rational function given by P(s)/Q(s).
504 ///
505 /// Let P[i] denote the coefficient of s^i in the polynomial P(s).
506 /// (P/Q)[r] =
507 /// if (r == 0) then
508 /// P[0]/Q[0]
509 /// else
510 /// (P[r] - {Σ_{i=1}^r (P/Q)[r-i] * Q[i])}/(Q[0])
511 /// We therefore recursively call getCoefficientInRationalFunction on
512 /// all i \in [0, power).
513 ///
514 /// https://math.ucdavis.edu/~deloera/researchsummary/
515 /// barvinokalgorithm-latte1.pdf, p. 1285
517  unsigned power, ArrayRef<QuasiPolynomial> num, ArrayRef<Fraction> den) {
518  assert(!den.empty() && "division by empty denominator in rational function!");
519
520  unsigned numParam = num[0].getNumInputs();
521  // We use the isEqual method of PresburgerSpace, which QuasiPolynomial
522  // inherits from.
523  assert(
524  llvm::all_of(
525  num, [&](const QuasiPolynomial &qp) { return num[0].isEqual(qp); }) &&
526  "the quasipolynomials should all belong to the same space!");
527
528  std::vector<QuasiPolynomial> coefficients;
529  coefficients.reserve(power + 1);
530
531  coefficients.push_back(num[0] / den[0]);
532  for (unsigned i = 1; i <= power; ++i) {
533  // If the power is not there in the numerator, the coefficient is zero.
534  coefficients.push_back(i < num.size() ? num[i]
535  : QuasiPolynomial(numParam, 0));
536
537  // After den.size(), the coefficients are zero, so we stop
538  // subtracting at that point (if it is less than i).
539  unsigned limit = std::min<unsigned long>(i, den.size() - 1);
540  for (unsigned j = 1; j <= limit; ++j)
541  coefficients[i] = coefficients[i] -
542  coefficients[i - j] * QuasiPolynomial(numParam, den[j]);
543
544  coefficients[i] = coefficients[i] / den[0];
545  }
546  return coefficients[power].simplify();
547 }
548
549 /// Substitute x_i = t^μ_i in one term of a generating function, returning
550 /// a quasipolynomial which represents the exponent of the numerator
551 /// of the result, and a vector which represents the exponents of the
552 /// denominator of the result.
553 /// If the returned value is {num, dens}, it represents the function
554 /// t^num / \prod_j (1 - t^dens[j]).
555 /// v represents the affine functions whose floors are multiplied by the
556 /// generators, and ds represents the list of generators.
557 std::pair<QuasiPolynomial, std::vector<Fraction>>
558 substituteMuInTerm(unsigned numParams, const ParamPoint &v,
559  const std::vector<Point> &ds, const Point &mu) {
560  unsigned numDims = mu.size();
561 #ifndef NDEBUG
562  for (const Point &d : ds)
563  assert(d.size() == numDims &&
564  "μ has to have the same number of dimensions as the generators!");
565 #endif
566
567  // First, the exponent in the numerator becomes
568  // - (μ • u_1) * (floor(first col of v))
569  // - (μ • u_2) * (floor(second col of v)) - ...
570  // - (μ • u_d) * (floor(d'th col of v))
571  // So we store the negation of the dot products.
572
573  // We have d terms, each of whose coefficient is the negative dot product.
574  SmallVector<Fraction> coefficients;
575  coefficients.reserve(numDims);
576  for (const Point &d : ds)
577  coefficients.push_back(-dotProduct(mu, d));
578
579  // Then, the affine function is a single floor expression, given by the
580  // corresponding column of v.
581  ParamPoint vTranspose = v.transpose();
582  std::vector<std::vector<SmallVector<Fraction>>> affine;
583  affine.reserve(numDims);
584  for (unsigned j = 0; j < numDims; ++j)
585  affine.push_back({SmallVector<Fraction>(vTranspose.getRow(j))});
586
587  QuasiPolynomial num(numParams, coefficients, affine);
588  num = num.simplify();
589
590  std::vector<Fraction> dens;
591  dens.reserve(ds.size());
592  // Similarly, each term in the denominator has exponent
593  // given by the dot product of μ with u_i.
594  for (const Point &d : ds) {
595  // This term in the denominator is
596  // (1 - t^dens.back())
597  dens.push_back(dotProduct(d, mu));
598  }
599
600  return {num, dens};
601 }
602
603 /// Normalize all denominator exponents dens to their absolute values
604 /// by multiplying and dividing by the inverses, in a function of the form
605 /// sign * t^num / prod_j (1 - t^dens[j]).
606 /// Here, sign = ± 1,
607 /// num is a QuasiPolynomial, and
608 /// each dens[j] is a Fraction.
610  std::vector<Fraction> &dens) {
611  // We track the number of exponents that are negative in the
612  // denominator, and convert them to their absolute values.
613  unsigned numNegExps = 0;
614  Fraction sumNegExps(0, 1);
615  for (const auto &den : dens) {
616  if (den < 0) {
617  numNegExps += 1;
618  sumNegExps += den;
619  }
620  }
621
622  // If we have (1 - t^-c) in the denominator, for positive c,
623  // multiply and divide by t^c.
624  // We convert all negative-exponent terms at once; therefore
625  // we multiply and divide by t^sumNegExps.
626  // Then we get
627  // -(1 - t^c) in the denominator,
628  // increase the numerator by c, and
629  // flip the sign of the function.
630  if (numNegExps % 2 == 1)
631  sign = -sign;
632  num = num - QuasiPolynomial(num.getNumInputs(), sumNegExps);
633 }
634
635 /// Compute the binomial coefficients nCi for 0 ≤ i ≤ r,
636 /// where n is a QuasiPolynomial.
637 std::vector<QuasiPolynomial> getBinomialCoefficients(const QuasiPolynomial &n,
638  unsigned r) {
639  unsigned numParams = n.getNumInputs();
640  std::vector<QuasiPolynomial> coefficients;
641  coefficients.reserve(r + 1);
642  coefficients.emplace_back(numParams, 1);
643  for (unsigned j = 1; j <= r; ++j)
644  // We use the recursive formula for binomial coefficients here and below.
645  coefficients.push_back(
646  (coefficients[j - 1] * (n - QuasiPolynomial(numParams, j - 1)) /
647  Fraction(j, 1))
648  .simplify());
649  return coefficients;
650 }
651
652 /// Compute the binomial coefficients nCi for 0 ≤ i ≤ r,
653 /// where n is a QuasiPolynomial.
654 std::vector<Fraction> getBinomialCoefficients(const Fraction &n,
655  const Fraction &r) {
656  std::vector<Fraction> coefficients;
657  coefficients.reserve((int64_t)floor(r));
658  coefficients.emplace_back(1);
659  for (unsigned j = 1; j <= r; ++j)
660  coefficients.push_back(coefficients[j - 1] * (n - (j - 1)) / (j));
661  return coefficients;
662 }
663
664 /// We have a generating function of the form
665 /// f_p(x) = \sum_i sign_i * (x^n_i(p)) / (\prod_j (1 - x^d_{ij})
666 ///
667 /// where sign_i is ±1,
668 /// n_i \in Q^p -> Q^d is the sum of the vectors d_{ij}, weighted by the
669 /// floors of d affine functions on p parameters.
670 /// d_{ij} \in Q^d are vectors.
671 ///
672 /// We need to find the number of terms of the form x^t in the expansion of
673 /// this function.
674 /// However, direct substitution (x = (1, ..., 1)) causes the denominator
675 /// to become zero.
676 ///
677 /// We therefore use the following procedure instead:
678 /// 1. Substitute x_i = (s+1)^μ_i for some vector μ. This makes the generating
679 /// function a function of a scalar s.
680 /// 2. Write each term in this function as P(s)/Q(s), where P and Q are
681 /// polynomials. P has coefficients as quasipolynomials in d parameters, while
682 /// Q has coefficients as scalars.
683 /// 3. Find the constant term in the expansion of each term P(s)/Q(s). This is
684 /// equivalent to substituting s = 0.
685 ///
686 /// Verdoolaege, Sven, et al. "Counting integer points in parametric
687 /// polytopes using Barvinok's rational functions." Algorithmica 48 (2007):
688 /// 37-66.
691  // Step (1) We need to find a μ such that we can substitute x_i =
692  // (s+1)^μ_i. After this substitution, the exponent of (s+1) in the
693  // denominator is (μ_i • d_{ij}) in each term. Clearly, this cannot become
694  // zero. Hence we find a vector μ that is not orthogonal to any of the
695  // d_{ij} and substitute x accordingly.
696  std::vector<Point> allDenominators;
697  for (ArrayRef<Point> den : gf.getDenominators())
698  allDenominators.insert(allDenominators.end(), den.begin(), den.end());
699  Point mu = getNonOrthogonalVector(allDenominators);
700
701  unsigned numParams = gf.getNumParams();
702  const std::vector<std::vector<Point>> &ds = gf.getDenominators();
703  QuasiPolynomial totalTerm(numParams, 0);
704  for (unsigned i = 0, e = ds.size(); i < e; ++i) {
705  int sign = gf.getSigns()[i];
706
707  // Compute the new exponents of (s+1) for the numerator and the
708  // denominator after substituting μ.
709  auto [numExp, dens] =
710  substituteMuInTerm(numParams, gf.getNumerators()[i], ds[i], mu);
711  // Now the numerator is (s+1)^numExp
712  // and the denominator is \prod_j (1 - (s+1)^dens[j]).
713
714  // Step (2) We need to express the terms in the function as quotients of
715  // polynomials. Each term is now of the form
716  // sign_i * (s+1)^numExp / (\prod_j (1 - (s+1)^dens[j]))
717  // For the i'th term, we first normalize the denominator to have only
718  // positive exponents. We convert all the dens[j] to their
719  // absolute values and change the sign and exponent in the numerator.
720  normalizeDenominatorExponents(sign, numExp, dens);
721
722  // Then, using the formula for geometric series, we replace each (1 -
723  // (s+1)^(dens[j])) with
724  // (-s)(\sum_{0 ≤ k < dens[j]} (s+1)^k).
725  for (auto &j : dens)
726  j = abs(j) - 1;
727  // Note that at this point, the semantics of dens[j] changes to mean
728  // a term (\sum_{0 ≤ k ≤ dens[j]} (s+1)^k). The denominator is, as before,
729  // a product of these terms.
730
731  // Since the -s are taken out, the sign changes if there is an odd number
732  // of such terms.
733  unsigned r = dens.size();
734  if (dens.size() % 2 == 1)
735  sign = -sign;
736
737  // Thus the term overall now has the form
738  // sign'_i * (s+1)^numExp /
739  // (s^r * \prod_j (\sum_{0 ≤ k < dens[j]} (s+1)^k)).
740  // This means that
741  // the numerator is a polynomial in s, with coefficients as
742  // quasipolynomials (given by binomial coefficients), and the denominator
743  // is a polynomial in s, with integral coefficients (given by taking the
744  // convolution over all j).
745
746  // Step (3) We need to find the constant term in the expansion of each
747  // term. Since each term has s^r as a factor in the denominator, we avoid
748  // substituting s = 0 directly; instead, we find the coefficient of s^r in
749  // sign'_i * (s+1)^numExp / (\prod_j (\sum_k (s+1)^k)),
750  // Letting P(s) = (s+1)^numExp and Q(s) = \prod_j (...),
751  // we need to find the coefficient of s^r in P(s)/Q(s),
752  // for which we use the getCoefficientInRationalFunction() function.
753
754  // First, we compute the coefficients of P(s), which are binomial
755  // coefficients.
756  // We only need the first r+1 of these, as higher-order terms do not
757  // contribute to the coefficient of s^r.
758  std::vector<QuasiPolynomial> numeratorCoefficients =
759  getBinomialCoefficients(numExp, r);
760
761  // Then we compute the coefficients of each individual term in Q(s),
762  // which are (dens[i]+1) C (k+1) for 0 ≤ k ≤ dens[i].
763  std::vector<std::vector<Fraction>> eachTermDenCoefficients;
764  std::vector<Fraction> singleTermDenCoefficients;
765  eachTermDenCoefficients.reserve(r);
766  for (const Fraction &den : dens) {
767  singleTermDenCoefficients = getBinomialCoefficients(den + 1, den + 1);
768  eachTermDenCoefficients.push_back(
769  ArrayRef<Fraction>(singleTermDenCoefficients).slice(1));
770  }
771
772  // Now we find the coefficients in Q(s) itself
773  // by taking the convolution of the coefficients
774  // of all the terms.
775  std::vector<Fraction> denominatorCoefficients;
776  denominatorCoefficients = eachTermDenCoefficients[0];
777  for (unsigned j = 1, e = eachTermDenCoefficients.size(); j < e; ++j)
778  denominatorCoefficients = multiplyPolynomials(denominatorCoefficients,
779  eachTermDenCoefficients[j]);
780
781  totalTerm =
782  totalTerm + getCoefficientInRationalFunction(r, numeratorCoefficients,
783  denominatorCoefficients) *
784  QuasiPolynomial(numParams, sign);
785  }
786
788 }
std::pair< QuasiPolynomial, std::vector< Fraction > > substituteMuInTerm(unsigned numParams, const ParamPoint &v, const std::vector< Point > &ds, const Point &mu)
Substitute x_i = t^μ_i in one term of a generating function, returning a quasipolynomial which repres...
Definition: Barvinok.cpp:558
std::vector< QuasiPolynomial > getBinomialCoefficients(const QuasiPolynomial &n, unsigned r)
Compute the binomial coefficients nCi for 0 ≤&#160;i ≤&#160;r, where n is a QuasiPolynomial.
Definition: Barvinok.cpp:637
void normalizeDenominatorExponents(int &sign, QuasiPolynomial &num, std::vector< Fraction > &dens)
Normalize all denominator exponents dens to their absolute values by multiplying and dividing by the ...
Definition: Barvinok.cpp:609
static Value max(ImplicitLocOpBuilder &builder, Value value, Value bound)
Fraction determinant(FracMatrix *inverse=nullptr) const
Definition: Matrix.cpp:587
MPInt determinant(IntMatrix *inverse=nullptr) const
Definition: Matrix.cpp:553
An IntegerPolyhedron represents the set of points from a PresburgerSpace that satisfy a list of affin...
An IntegerRelation represents the set of points from a PresburgerSpace that satisfy a list of affine ...
Adds an inequality (>= 0) from the coefficients specified in inEq.
MPInt atIneq(unsigned i, unsigned j) const
Returns the value at the specified inequality row and column.
unsigned getNumCols() const
Returns the number of columns in the constraint system.
This class provides support for multi-precision arithmetic.
Definition: MPInt.h:87
unsigned getNumRows() const
Definition: Matrix.h:85
Matrix< T > getSubMatrix(unsigned fromRow, unsigned toRow, unsigned fromColumn, unsigned toColumn) const
Definition: Matrix.cpp:386
void scaleRow(unsigned row, const T &scale)
Multiply the specified row by a factor of scale.
Definition: Matrix.cpp:315
void insertColumn(unsigned pos)
Definition: Matrix.cpp:150
MutableArrayRef< T > getRow(unsigned row)
Get a [Mutable]ArrayRef corresponding to the specified row.
Definition: Matrix.cpp:132
void setRow(unsigned row, ArrayRef< T > elems)
Set the specified row to elems.
Definition: Matrix.cpp:142
std::pair< Matrix< T >, Matrix< T > > splitByBitset(ArrayRef< int > indicator)
Split the rows of a matrix into two matrices according to which bits are 1 and which are 0 in a given...
Definition: Matrix.cpp:414
void removeRow(unsigned pos)
Definition: Matrix.cpp:232
unsigned getNumColumns() const
Definition: Matrix.h:87
Matrix< T > transpose() const
Definition: Matrix.cpp:82
SmallVector< T, 8 > preMultiplyWithRow(ArrayRef< T > rowVec) const
The given vector is interpreted as a row vector v.
Definition: Matrix.cpp:348
void negateMatrix()
Negate the entire matrix.
Definition: Matrix.cpp:342
void swapRows(unsigned row, unsigned otherRow)
Swap the given rows.
Definition: Matrix.cpp:112
void addToRow(unsigned sourceRow, unsigned targetRow, const T &scale)
Add scale multiples of the source row to the target row.
Definition: Matrix.cpp:301
void negateRow(unsigned row)
Negate the specified row.
Definition: Matrix.cpp:336
T & at(unsigned row, unsigned column)
Access the element at the specified row and column.
Definition: Matrix.h:61
bool isFullDim() const
Return whether the given PresburgerRelation is full-dimensional.
PresburgerSet intersect(const PresburgerRelation &set) const
static PresburgerSet getUniverse(const PresburgerSpace &space)
Return a universe set of the specified type that contains all points.
static PresburgerSpace getSetSpace(unsigned numDims=0, unsigned numSymbols=0, unsigned numLocals=0)
static PresburgerSpace getRelationSpace(unsigned numDomain=0, unsigned numRange=0, unsigned numSymbols=0, unsigned numLocals=0)
std::vector< ParamPoint > getNumerators() const
std::vector< std::vector< Point > > getDenominators() const
std::vector< std::pair< PresburgerSet, GeneratingFunction > > computeChamberDecomposition(unsigned numSymbols, ArrayRef< std::pair< PresburgerSet, GeneratingFunction >> regionsAndGeneratingFunctions)
Given a list of possibly intersecting regions (PresburgerSet) and the generating functions active in ...
Definition: Barvinok.cpp:241
std::optional< ParamPoint > solveParametricEquations(FracMatrix equations)
Find the solution of a set of equations that express affine constraints between a set of variables an...
Definition: Barvinok.cpp:163
SlowMPInt abs(const SlowMPInt &x)
Redeclarations of friend declarations above to make it discoverable by lookups.
SmallVector< Fraction > Point
ConeV getDual(ConeH cone)
Given a cone in H-representation, return its dual.
Definition: Barvinok.cpp:22
QuasiPolynomial computeNumTerms(const GeneratingFunction &gf)
Find the number of terms in a generating function, as a quasipolynomial in the parameter space of the...
Definition: Barvinok.cpp:690
Point getNonOrthogonalVector(ArrayRef< Point > vectors)
Find a vector that is not orthogonal to any of the given vectors, i.e., has nonzero dot product with ...
Definition: Barvinok.cpp:473
GeneratingFunction computeUnimodularConeGeneratingFunction(ParamPoint vertex, int sign, const ConeH &cone)
Compute the generating function for a unimodular cone.
Definition: Barvinok.cpp:81
PolyhedronH defineHRep(int numVars, int numSymbols=0)
Definition: Barvinok.h:52
std::vector< std::pair< PresburgerSet, GeneratingFunction > > computePolytopeGeneratingFunction(const PolyhedronH &poly)
Compute the generating function corresponding to a polytope.
Definition: Barvinok.cpp:314
MPInt getIndex(const ConeV &cone)
Get the index of a cone, i.e., the volume of the parallelepiped spanned by its generators,...
Definition: Barvinok.cpp:64
QuasiPolynomial getCoefficientInRationalFunction(unsigned power, ArrayRef< QuasiPolynomial > num, ArrayRef< Fraction > den)
Find the coefficient of a given power of s in a rational function given by P(s)/Q(s),...
Definition: Barvinok.cpp:516
Fraction dotProduct(ArrayRef< Fraction > a, ArrayRef< Fraction > b)
Compute the dot product of two vectors.
Definition: Utils.cpp:533
std::vector< Fraction > multiplyPolynomials(ArrayRef< Fraction > a, ArrayRef< Fraction > b)
Find the product of two polynomials, each given by an array of coefficients.
Definition: Utils.cpp:544
MPInt floor(const Fraction &f)
Definition: Fraction.h:74
Include the generated interface declarations.
A class to represent fractions.
Definition: Fraction.h:28
Eliminates variable at the specified position using Fourier-Motzkin variable elimination.