//===- Barvinok.cpp - Barvinok's Algorithm ----------------------*- C++ -*-===// // // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. // See https://llvm.org/LICENSE.txt for license information. // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception // //===----------------------------------------------------------------------===// #include "mlir/Analysis/Presburger/Barvinok.h" #include "mlir/Analysis/Presburger/Utils.h" #include "llvm/ADT/Sequence.h" #include <algorithm> usingnamespacemlir; usingnamespacepresburger; usingnamespacemlir::presburger::detail; /// Assuming that the input cone is pointed at the origin, /// converts it to its dual in V-representation. /// Essentially we just remove the all-zeroes constant column. ConeV mlir::presburger::detail::getDual(ConeH cone) { … } /// Converts a cone in V-representation to the H-representation /// of its dual, pointed at the origin (not at the original vertex). /// Essentially adds a column consisting only of zeroes to the end. ConeH mlir::presburger::detail::getDual(ConeV cone) { … } /// Find the index of a cone in V-representation. DynamicAPInt mlir::presburger::detail::getIndex(const ConeV &cone) { … } /// Compute the generating function for a unimodular cone. /// This consists of a single term of the form /// sign * x^num / prod_j (1 - x^den_j) /// /// sign is either +1 or -1. /// den_j is defined as the set of generators of the cone. /// num is computed by expressing the vertex as a weighted /// sum of the generators, and then taking the floor of the /// coefficients. GeneratingFunction mlir::presburger::detail::computeUnimodularConeGeneratingFunction( ParamPoint vertex, int sign, const ConeH &cone) { … } /// We use Gaussian elimination to find the solution to a set of d equations /// of the form /// a_1 x_1 + ... + a_d x_d + b_1 m_1 + ... + b_p m_p + c = 0 /// where x_i are variables, /// m_i are parameters and /// a_i, b_i, c are rational coefficients. /// /// The solution expresses each x_i as an affine function of the m_i, and is /// therefore represented as a matrix of size d x (p+1). /// If there is no solution, we return null. std::optional<ParamPoint> mlir::presburger::detail::solveParametricEquations(FracMatrix equations) { … } /// This is an implementation of the Clauss-Loechner algorithm for chamber /// decomposition. /// /// We maintain a list of pairwise disjoint chambers and the generating /// functions corresponding to each one. We iterate over the list of regions, /// each time adding the current region's generating function to the chambers /// where it is active and separating the chambers where it is not. /// /// Given the region each generating function is active in, for each subset of /// generating functions the region that (the sum of) precisely this subset is /// in, is the intersection of the regions that these are active in, /// intersected with the complements of the remaining regions. std::vector<std::pair<PresburgerSet, GeneratingFunction>> mlir::presburger::detail::computeChamberDecomposition( unsigned numSymbols, ArrayRef<std::pair<PresburgerSet, GeneratingFunction>> regionsAndGeneratingFunctions) { … } /// For a polytope expressed as a set of n inequalities, compute the generating /// function corresponding to the lattice points included in the polytope. This /// algorithm has three main steps: /// 1. Enumerate the vertices, by iterating over subsets of inequalities and /// checking for satisfiability. For each d-subset of inequalities (where d /// is the number of variables), we solve to obtain the vertex in terms of /// the parameters, and then check for the region in parameter space where /// this vertex satisfies the remaining (n - d) inequalities. /// 2. For each vertex, identify the tangent cone and compute the generating /// function corresponding to it. The generating function depends on the /// parametric expression of the vertex and the (non-parametric) generators /// of the tangent cone. /// 3. [Clauss-Loechner decomposition] Identify the regions in parameter space /// (chambers) where each vertex is active, and accordingly compute the /// GF of the polytope in each chamber. /// /// Verdoolaege, Sven, et al. "Counting integer points in parametric /// polytopes using Barvinok's rational functions." Algorithmica 48 (2007): /// 37-66. std::vector<std::pair<PresburgerSet, GeneratingFunction>> mlir::presburger::detail::computePolytopeGeneratingFunction( const PolyhedronH &poly) { … } /// We use an iterative procedure to find a vector not orthogonal /// to a given set, ignoring the null vectors. /// Let the inputs be {x_1, ..., x_k}, all vectors of length n. /// /// In the following, /// vs[:i] means the elements of vs up to and including the i'th one, /// <vs, us> means the dot product of vs and us, /// vs ++ [v] means the vector vs with the new element v appended to it. /// /// We proceed iteratively; for steps d = 0, ... n-1, we construct a vector /// which is not orthogonal to any of {x_1[:d], ..., x_n[:d]}, ignoring /// the null vectors. /// At step d = 0, we let vs = [1]. Clearly this is not orthogonal to /// any vector in the set {x_1[0], ..., x_n[0]}, except the null ones, /// which we ignore. /// At step d > 0 , we need a number v /// s.t. <x_i[:d], vs++[v]> != 0 for all i. /// => <x_i[:d-1], vs> + x_i[d]*v != 0 /// => v != - <x_i[:d-1], vs> / x_i[d] /// We compute this value for all x_i, and then /// set v to be the maximum element of this set plus one. Thus /// v is outside the set as desired, and we append it to vs /// to obtain the result of the d'th step. Point mlir::presburger::detail::getNonOrthogonalVector( ArrayRef<Point> vectors) { … } /// We use the following recursive formula to find the coefficient of /// s^power in the rational function given by P(s)/Q(s). /// /// Let P[i] denote the coefficient of s^i in the polynomial P(s). /// (P/Q)[r] = /// if (r == 0) then /// P[0]/Q[0] /// else /// (P[r] - {Σ_{i=1}^r (P/Q)[r-i] * Q[i])}/(Q[0]) /// We therefore recursively call `getCoefficientInRationalFunction` on /// all i \in [0, power). /// /// https://math.ucdavis.edu/~deloera/researchsummary/ /// barvinokalgorithm-latte1.pdf, p. 1285 QuasiPolynomial mlir::presburger::detail::getCoefficientInRationalFunction( unsigned power, ArrayRef<QuasiPolynomial> num, ArrayRef<Fraction> den) { … } /// Substitute x_i = t^μ_i in one term of a generating function, returning /// a quasipolynomial which represents the exponent of the numerator /// of the result, and a vector which represents the exponents of the /// denominator of the result. /// If the returned value is {num, dens}, it represents the function /// t^num / \prod_j (1 - t^dens[j]). /// v represents the affine functions whose floors are multiplied by the /// generators, and ds represents the list of generators. std::pair<QuasiPolynomial, std::vector<Fraction>> substituteMuInTerm(unsigned numParams, const ParamPoint &v, const std::vector<Point> &ds, const Point &mu) { … } /// Normalize all denominator exponents `dens` to their absolute values /// by multiplying and dividing by the inverses, in a function of the form /// sign * t^num / prod_j (1 - t^dens[j]). /// Here, sign = ± 1, /// num is a QuasiPolynomial, and /// each dens[j] is a Fraction. void normalizeDenominatorExponents(int &sign, QuasiPolynomial &num, std::vector<Fraction> &dens) { … } /// Compute the binomial coefficients nCi for 0 ≤ i ≤ r, /// where n is a QuasiPolynomial. std::vector<QuasiPolynomial> getBinomialCoefficients(const QuasiPolynomial &n, unsigned r) { … } /// Compute the binomial coefficients nCi for 0 ≤ i ≤ r, /// where n is a QuasiPolynomial. std::vector<Fraction> getBinomialCoefficients(const Fraction &n, const Fraction &r) { … } /// We have a generating function of the form /// f_p(x) = \sum_i sign_i * (x^n_i(p)) / (\prod_j (1 - x^d_{ij}) /// /// where sign_i is ±1, /// n_i \in Q^p -> Q^d is the sum of the vectors d_{ij}, weighted by the /// floors of d affine functions on p parameters. /// d_{ij} \in Q^d are vectors. /// /// We need to find the number of terms of the form x^t in the expansion of /// this function. /// However, direct substitution (x = (1, ..., 1)) causes the denominator /// to become zero. /// /// We therefore use the following procedure instead: /// 1. Substitute x_i = (s+1)^μ_i for some vector μ. This makes the generating /// function a function of a scalar s. /// 2. Write each term in this function as P(s)/Q(s), where P and Q are /// polynomials. P has coefficients as quasipolynomials in d parameters, while /// Q has coefficients as scalars. /// 3. Find the constant term in the expansion of each term P(s)/Q(s). This is /// equivalent to substituting s = 0. /// /// Verdoolaege, Sven, et al. "Counting integer points in parametric /// polytopes using Barvinok's rational functions." Algorithmica 48 (2007): /// 37-66. QuasiPolynomial mlir::presburger::detail::computeNumTerms(const GeneratingFunction &gf) { … }
Codebase Browser
llvm