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