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
9#include "mlir/Analysis/Presburger/Barvinok.h"
10#include "mlir/Analysis/Presburger/Utils.h"
11#include "llvm/ADT/Sequence.h"
12#include <algorithm>
13#include <bitset>
14
15using namespace mlir;
16using namespace presburger;
17using 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.
22ConeV mlir::presburger::detail::getDual(ConeH cone) {
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>(Begin: 0, End: 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(row: i, column: 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.
47ConeH mlir::presburger::detail::getDual(ConeV cone) {
48 unsigned rows = cone.getNumRows();
49 unsigned columns = cone.getNumColumns();
50 ConeH dual = defineHRep(numVars: columns);
51 // Add a new column (for constants) at the end.
52 // This will be initialized to zero.
53 cone.insertColumn(pos: columns);
54
55 for (unsigned i = 0; i < rows; ++i)
56 dual.addInequality(inEq: cone.getRow(row: 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.
64MPInt mlir::presburger::detail::getIndex(const ConeV &cone) {
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.
80GeneratingFunction
81mlir::presburger::detail::computeUnimodularConeGeneratingFunction(
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.
98 FracMatrix transp = FracMatrix(cone.getInequalities()).transpose();
99 transp.removeRow(pos: 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>(Begin: 0, End: numVar)) {
112 row = generators.getRow(row: 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>(Begin: 0, End: numColumns)) {
134 for (auto j : llvm::seq<int>(Begin: 0, End: numRows))
135 ithCol[j] = vertex(j, i);
136 numerator.setRow(row: i, elems: transp.preMultiplyWithRow(rowVec: ithCol));
137 numerator.negateRow(row: 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.
162std::optional<ParamPoint>
163mlir::presburger::detail::solveParametricEquations(FracMatrix equations) {
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(row: j, otherRow: 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);
202 equations.addToRow(/*sourceRow=*/i, /*targetRow=*/j,
203 /*scale=*/-currentElement / diagElement);
204 }
205 }
206
207 // Rescale diagonal elements to 1.
208 for (unsigned i = 0; i < d; ++i)
209 equations.scaleRow(row: i, scale: 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.
240std::vector<std::pair<PresburgerSet, GeneratingFunction>>
241mlir::presburger::detail::computeChamberDecomposition(
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 = {
249 {PresburgerSet::getUniverse(space: PresburgerSpace::getSetSpace(numDims: numSymbols)),
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(set: 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(args: currentRegion, args: currentGeneratingFunction);
279 continue;
280 }
281
282 // If it is, we add the intersection and the difference as chambers.
283 newChambers.emplace_back(args&: intersection,
284 args: currentGeneratingFunction + generatingFunction);
285 newChambers.emplace_back(args: currentRegion.subtract(set: region),
286 args: 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.
313std::vector<std::pair<PresburgerSet, GeneratingFunction>>
314mlir::presburger::detail::computePolytopeGeneratingFunction(
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(fromRow: 0, toRow: numIneqs - numVars - 1, fromColumn: 0, toColumn: numVars - 1));
351 b2c2 = FracMatrix(remainder.getSubMatrix(fromRow: 0, toRow: numIneqs - numVars - 1, fromColumn: numVars,
352 toColumn: numVars + numSymbols));
353
354 // Find the vertex, if any, corresponding to the current subset of
355 // inequalities.
356 std::optional<ParamPoint> vertex =
357 solveParametricEquations(equations: FracMatrix(subset)); // d x (p+1)
358
359 if (!vertex)
360 continue;
361 if (std::find(first: vertices.begin(), last: vertices.end(), val: 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(x: *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(row: i, elems: vertex->preMultiplyWithRow(rowVec: a2.getRow(row: i)));
398 activeRegion.addToRow(row: i, rowVec: b2c2.getRow(row: i), scale: 1);
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(numDomain: 0, numRange: numSymbols, numSymbols: 0, numLocals: 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);
419 tangentCone.addInequality(inEq: ineq);
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: *vertex, sign, cone);
432 }
433 // We store the vertex we computed with the generating function of its
434 // tangent cone.
435 regionsAndGeneratingFunctions.emplace_back(args: PresburgerSet(activeRegionRel),
436 args&: vertexGf);
437 } while (std::next_permutation(first: indicator.begin(), last: 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.
473Point mlir::presburger::detail::getNonOrthogonalVector(
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(a: ArrayRef(vector).slice(N: 0, M: d), b: newPoint) / vector[d];
493
494 // Find the biggest such value
495 maxDisallowedValue = std::max(a: maxDisallowedValue, b: disallowedValue);
496 }
497 newPoint.push_back(Elt: 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
516QuasiPolynomial mlir::presburger::detail::getCoefficientInRationalFunction(
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(n: power + 1);
530
531 coefficients.push_back(x: 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(x: 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>(a: i, b: 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.
557std::pair<QuasiPolynomial, std::vector<Fraction>>
558substituteMuInTerm(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(N: numDims);
576 for (const Point &d : ds)
577 coefficients.push_back(Elt: -dotProduct(a: mu, b: 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(n: numDims);
584 for (unsigned j = 0; j < numDims; ++j)
585 affine.push_back(x: {SmallVector<Fraction>(vTranspose.getRow(row: j))});
586
587 QuasiPolynomial num(numParams, coefficients, affine);
588 num = num.simplify();
589
590 std::vector<Fraction> dens;
591 dens.reserve(n: 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(x: dotProduct(a: d, b: 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.
609void normalizeDenominatorExponents(int &sign, QuasiPolynomial &num,
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.
637std::vector<QuasiPolynomial> getBinomialCoefficients(const QuasiPolynomial &n,
638 unsigned r) {
639 unsigned numParams = n.getNumInputs();
640 std::vector<QuasiPolynomial> coefficients;
641 coefficients.reserve(n: r + 1);
642 coefficients.emplace_back(args&: numParams, args: 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 x: (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.
654std::vector<Fraction> getBinomialCoefficients(const Fraction &n,
655 const Fraction &r) {
656 std::vector<Fraction> coefficients;
657 coefficients.reserve(n: (int64_t)floor(f: r));
658 coefficients.emplace_back(args: 1);
659 for (unsigned j = 1; j <= r; ++j)
660 coefficients.push_back(x: 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.
689QuasiPolynomial
690mlir::presburger::detail::computeNumTerms(const GeneratingFunction &gf) {
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(position: allDenominators.end(), first: den.begin(), last: den.end());
699 Point mu = getNonOrthogonalVector(vectors: 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, v: gf.getNumerators()[i], ds: 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, num&: 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(f: 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(n: 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(n: r);
766 for (const Fraction &den : dens) {
767 singleTermDenCoefficients = getBinomialCoefficients(n: den + 1, r: den + 1);
768 eachTermDenCoefficients.push_back(
769 x: ArrayRef<Fraction>(singleTermDenCoefficients).slice(N: 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(a: denominatorCoefficients,
779 b: eachTermDenCoefficients[j]);
780
781 totalTerm =
782 totalTerm + getCoefficientInRationalFunction(power: r, num: numeratorCoefficients,
783 den: denominatorCoefficients) *
784 QuasiPolynomial(numParams, sign);
785 }
786
787 return totalTerm.simplify();
788}
789

source code of mlir/lib/Analysis/Presburger/Barvinok.cpp