Home

barvinok: User Guide

1. e e e e e e e S Figure 5 84 A polytope and its integer projections 2 L e e e S Figure 5 86 A transformed polytope and its integer projection 107 small In particular if the dimension m is fixed then the lattice width see sub section 5 23 of lattice point free polytopes is bounded by a constant w m Lagarias et al 1990 Banaszczyk et al 1999 This means that in the direction of the lattice width of a polytope the gaps will be not be larger than w m Barvinok and Woods 2003 Theorem 4 3 Otherwise we would be able to put a uniformly scaled down version of the polytope in the gap and it would contain no lattice points which would contradict the fact that its lattice width is bounded by w m contains such a scaled down copy of the original polytope However neither the hori zontal nor the vertical direction is a lattice width direction of this polytope The actual lattice width of this polytope was computed in Example as 3 with correspond ing direction 1 2 Figure 5 86 shows the result of applying the unimodular transformation 1 2 0 1 to the polytope Note that the horizontal direction now has gaps of width at most 1 After shifting subtracting and projecting in the vertical direction we therefore end up with a set S with gaps of width 1 and we then have to shift and
2. Figure 5 90 The parametric polytope from Example for different values of the parameter 113 x Figure 5 92 The transformed parametric polytope from Example forO lt p lt 5 114 The first step is to compute the lattice width of P p In Example we already obtained the decomposition of the parameter domain into Ci p 0 lt p lt 5 p 5 lt lt 1 with corresponding lattice widths and lattice width directions 0 1 We p 5 p 1 0 We p 5 p Note that in this example the gaps in both coordinate directions are 1 so in principle there is no need to perform any unimodular transformation here Nevertheless we will apply the transformation that would be applied by the algorithm On the first domain we extend the lattice width direction to the unimodular matrix After application to the existentially quantified variables x we obtain the parametric polytope 20 5 gt 0 2252 5 gt 0 4 2 5 gt 0 2x 5 gt 0 gt 0 shown in the top of Figure 5 92 We now temporarily remove the existential quan tification on 1 resulting in T p x1 Z 3o Z s x P Since there is only one existentially quantified variable left we know we only have to shift and subtract the set once to obt
3. _ Pommersheim J Pop S pos 23 pow power Preparata F P 93 724 primal space print print evalue projection theorem 34 Pugh W Pulleyblank W R 22 QQ quasi polynomial Ehrhart see Ehrhart quasi polynomial Queyranne M Rabl T Rambau J 87 24 ran range range map 16 rational generating function 60 Ray 22 B6 ray open see open ray Rays read 16 recession cone Redon X 13 reduce_evalue reduced basis generalized see generalized reduced basis reduction parameter regular triangulation relation 26 27 30 representation explicit see explicit representation respecting 18 revlex positive Rom W O 87 723 Rosser E 23 Rutherford T 27 Ryan C T 119 723 25 Sakellariou R 63 67 sample 2 Sannella D 119 scan Scarf H E 34 727 725 schedule schedule_forest 18 Schnorr C P 24 Seghir R 29 105 117 118 725 727 Shallcross D F 27 Shamos M I 93 724 Shmonin G 98 104 105 722 short rat 30 Shoup V Shpeisman 723 Silber G A 118 724 simple cone size 23 27 slack variable Smaoui H 118 724 Smith Normal Form SNF 45 solutions SolveDiophantine Sottile F 118 Source specialization standard form 85 Stanley P 44 53 723 Stefanov T stride Stroobandt D 117 sum 41
4. yp qeu ya y y l ew _ 1 x Grim xj n 0 xum d 5 b l 0 2490 lin m 2 n gt 0 nr 1 MI n gt 0 1 n no ni n2 1 5 37 52 0793 m 5 12 271 702 yi 1 1 1 zi is ba ia Oy ee s ON 3 q 497 4 4 99 199 The coefficient of y is then 1 371 51 1079 721 1 0 1 720 1212 22 720 240 which is the same value as the one we found in Example 5 50 To compare the difference between the old implementation described in subsec tion 5 15 and the new implementation described here we repeat the experiments of Ver doolaege and Bruynooghe 2008 The new experiments were performed on the same 84 1000 100 10 a o 1 E 0 1 P reduction to unweighted Barvinok 0 01 nested sums 4 Euler Maclaurin Laurent expansion old 0 001 Laurent expansion new 0 20 40 60 80 100 Figure 5 54 Execution times for summing a monomial over a difficult non parametric triangle hardware Intel Core2 as the old experiments but using barvinok version barvinok 9 It should be noted that the running times reported by Verdoolaege and Bruynooghe 2008 had been mistakenly scaled down by a factor of 0 6 The experiments were cut off at 10 minutes 600 seconds The new results are shown in Figures 5 54 and 5 55 The time measuring re
5. where L has the lines as columns and L an integral basis for the orthogonal comple ment Section 5 7 Note that the inverse transformation p L p 0 gt 14 has integral coefficients since L can be extended to a unimodular matrix If the parameter values p for which P contains integer points are restricted to a non standard lattice we first replace the parameters by a different set of parameters that lie on the standard lattice through parameter compression Meister 2004 The left inverse of can be computes as explained in Section giving We have to be careful to only apply this transformation when both the equalities puted in Section are satisfied and some additional divisibility constraints In par ticular if a d is a row of C with a Z and d Z the transformation can only be applied to parameter values p such that d divides a p The complete transformation is given by p CL p or L C p 57 5 9 Parametric Volume Computation The volume of a parametric polytope can serve as an approximation for the number of integer points in the polytope We basically follow the description of here except that we focus on volume computation for linearly parametrized polytopes which we exploit to determine the sign of the determinants we compute as explained below Note first that the vertices of a linearly parametrized polytope are affine expres si
6. Wet tees Pee Pee m vor eres e Es T 6 2 Publications Refering to the Library 120 128 List of Figures Ll Aloopnesf 1 2 The iteration domain of the loop nest in Figure 1 1 1 3 loop with non unit stridd 0 L4 A loop witha modulo condition 9 15 A program with three loop nests 22 quasi polynomial 1 2 p 3p 2 Ln 5 11 Possible locations of the vector w with respect to the rays of a 3 dimensional cone ee 5 17 Examples of decompositions in primal space 5 18 Possible locations of w with respect to u and projected onto a plane orthogonal to the other rays lt 0 5 19 Possible locations of w with respect to and projected onto a plane orthogonal to the other rays when ajaj gt 0 5 44 Sum of polynomial x over the integer points in a triangle 5 54 Execution times for summing a monomial over a difficult non triangle oo erpa aee a G a aae rss 5 55 Execution times for summing a monomial over a difficult parametric 5 68 The Hilbert basis and the integer hull of a truncated 5 70 integer hull of a truncated 5 73 A polytope and its candidate width direction 5 75 T
7. eos gt 52 gt gt 52 gt S1 lt lt S2 II EQUUS I 3 S SS 5 2 lt lt Mm 5 lt lt 55 lt lt 5 gt gt S2 gt gt the ith power of m if i is negative then the result is the i th power of the inverse of m compute an overapproximation of the transitive clo sure of m and return a list containing the overap proximation and a boolean that is true if the overap proximation is known to be exact the element at position i in the list is q obviously equal to q2 is f obviously equal to f2 is s equal to s2 is m equal to m2 15 s a subset of s2 15 m a subset of m5 is s a proper subset of s2 is m a proper subset of m5 is s a superset of s2 15 m a superset of m is s a proper superset of s2 is m a proper superset of m a map from s to s2 of those elements that live in the same space and such that the elements of s are lexicographically strictly smaller than those of s2 a map from the domain of m to the domain of mz of those elements such that their images live in the same space and such that the images of the el ements of m are lexicographically strictly smaller than those of m a map from s to s2 of those elements that live in the same space and such that the elements of s are lexicographically smaller than those of 52 a map f
8. isl give isl union pw qpolynomial isl union map apply union pw qpolynomial __151 take isl union map umap __isl_take isl union pw qpolynomial upwqp Compose the given map with the given piecewise quasipolynomial That is compute the sum over all elements in the intersection of the range of the map and the domain of the piecewise quasipolynomial as a function of an element in the domain of the map 1 2 Calculator The iscc calculator offers an interface to some of the functionality provided by the isl and barvinok libraries The language used by iscc is extremely simple The calcula tor supports operations on constants and dynamically typed variables and assignments 2 to those variables If the result of an expression is not used inside another ex pression and not assigned to a variable then this result is printed on the screen The operators are overloaded based on the types of the arguments which may be sets re lations piecewise quasipolynomials piecewise quasipolynomial folds lists strings or booleans The supported operations are shown in Table I Note that when an opera tion requires an argument of a certain type a binary list with the first element of the required type may also be used instead 1 2 4 Sets and Iteration Domains Within the polyhedral model for analysis and transformation of static affine programs the most basic kind of set is the iteration domain The iteration domain represents the iterations of
9. type polynomial size 3 pos 1 enode arr 0 d 2 type periodic x n 5 d 1 size 2 1 pos 1 x n 3 d 1 arr 2 8 arr 6 x n 1 1 1 x n 2 Figure 2 2 The quasi polynomial 1 2 p 3p Example 2 1 Figure 2 2 is skillful reconstruction of Figure 2 from Loechner 1999 It shows the contents of the enode structures representing the quasi polynomial 1 2 p 3p 3 2 2 2 Options The barvinok options structure contains various options that influence the behavior of the library struct barvinok options struct barvinok stats stats PolyLib options unsigned MaxRays NTL options LLL reduction parameter delta LLL a LLL b long LLL a long LLL b barvinok options define BV SPECIALIZATION BF 2 define BV SPECIALIZATION DF 1 define BV SPECIALIZATION RANDOM 6 define BV SPECIALIZATION TODD 3 int incremental specialization unsigned long max index int primal int lookup table 24 int count_sample_infinite int try_Delaunay_triangulation define BV_APPROX_SIGN_NONE 0 define BV_APPROX_SIGN_APPROX 1 define BV_APPROX_SIGN_LOWER 2 define BV_APPROX_SIGN_UPPER 3 int polynomial_approximation define BV_APPROX_NONE 0 define BV_APPROX_DROP 1 define BV_APPROX_SCALE 2 define BV_APPROX_VOLUME 3 define BV_APPROX_BERNOULLI 4 int approximation_method define BV_APPROX_SCALE_FAST 1 lt lt 8 define BV_APPROX_SCALE_NARROW 1 lt lt
10. barvinok bound devos pwqp U gt max 2 3 U lt 10 and U gt Q Options lower compute lower bound instead of upper bound 3 11 polytope minimize The program polytope minimize has been superseded by isl sisl polyhedron minimize 3 12 polyhedron integer hull The program polyhedron integer hull takes a polyhedron in PolyLib notation and prints its integer hull The integer hull is computed as explained in subsection 5 21 42 3 13 polytope_lattice_width The program polytope_lattice_width computes the lattice width of a parametric polytope The input is the same as that of barvinok_enumerate The lattice width is computed as explained in subsection 5 23 Options direction d printthe lattice width directions 43 4 polymake clients The barvinok distribution includes a couple of polymake Gawrilow and Joswig 2000 clients in the polymake subdir e lattice points file Computes the property LATTICE POINTS of a polytope the number of lattice points in the polytope e h star vector lt file gt Computes the property H STAR VECTOR of a lattice polytope the h vector of the polytope Stanley 1993 44 5 Implementation details 5 1 Aninterior point of a polyhedron We often need a point that lies in the interior of a polyhedron The function inner point implements the following algorithm Each polyhedron P can be written as the sum of a polytope P and a cone C the recession cone or characteri
11. qi G z gt Td i then xAvita FG 2 ST cy 1 The method gen_fun Hadamard_product computes the Hadamard product of the current rational generating function with the rational generating function gf as ex plained in Verdoolaege 2005 Section 4 5 2 32 2 6 Counting Functions Our library provides essentially three different counting functions one for non parametric polytopes one for parametric polytopes and one for parametric sets with existential variables The old versions of these functions have a MaxRays argument while the new versions have a more general barvinok options argument For more informa tion on barvinok options see Section 2 2 void barvinok count Polyhedron P Value result unsigned NbMaxCons void barvinok count with options Polyhedron P Value result struct barvinok options options The function barvinok count or barvinok count with options enumerates the non parametric polytope P and returns the result in the Value pointed to by result which needs to have been allocated and initialized If P is a union then only the first set in the union will be taken into account For the meaning of the argument NbMaxCons see the discussion on MaxRays in Section 2 2 The function barvinok enumerate for enumerating parametric polytopes was meant to be a drop in replacement of PolyLib s Polyhedron Enumerate function Unfortunately the latter has been changed to accept an extra argu
12. ROG A Finally we consider the contributions of the vertices For the vertex 3 4 we have r 1 0 1 1 Since v is integer we have ti t 0 Also Cy 1 1 2 y D and y D D Since the total degree of the polynomial xix is two we only need the coefficients of u K up ton n 2 ni n 010 ECT Di PED je Di 1 10 2 12 buc Dy ola aga yen ge EED Dp 2 0 ae ar Ca 1 1 G28 9 1 69 6 1 CD Di Da 0 2 SP 00 890 55 Cn Da We find 3 1 h x D D D D 2D D2 x s i xc js C 1D3 at js 3 1 B 7 5 ge E EE gee a QU Regg a The contribution of this vertex is therefore 1355 we have ti t 0 Also 1 1 2 y D and y2 D D We similarly For the vertex v 2 5 we have rj 0 1 and r2 1 1 Since v is integer find 7 5 5750 7 1152 1 1 h x sux xX a x2 X1 576 24 24 The contribution of this vertex is therefore 29 1067 76 For the vertex v 2 4 we have r 1 0 and rz 0 1 Since v is integer we have t t 0 The computations are easier in this case since Cz 0 y Di and D We find 2 1 1 1 1 h x 4 12 12 T 144 The contribution of this vertex is therefore 253 h 2 4 144 The
13. f xy f x y 1 The only remaining problem is that the operation in is fairly expensive In particular this operation is performed by first computing the Hadamard product of the two generating functions which corresponds to the intersection of the sets and then subtracting the resulting generating function from this first generating function The last operation is fairly cheap but the Hadamard product has a time complexity which while polynomial if the dimension in this case the maximum of k in 5 26 is fixed is exponential in this dimension Furthermore this dimension increases by max k d on each successive application of the Hadamard product while max k gt d as soon as some projection is performed which certainly happens in the recursive application of the algorithm Now the total number of Hadamard products is bounded by a constant w m and so the increase in dimension is also bounded by a constant so the whole operation is still polynomial for fixed dimension Nevertheless we do not want to perform any additional Hadamard products if we do not really have to That is we would like to be able to detect when we can stop performing these operations before reaching the upper bound Let f x z f z and let f x z be the result of applying the set differ ence in k times Denote the corresponding sets by and We want to find the smallest k such that f x y z f
14. gt i exists f 1 lt i lt n and 1 5 f lt 2 and 1 1 4 lt 5 lt 1 Since greatest integer parts occur relatively frequently there is a special notation for them in isl using The above set can therefore also be represented as n gt i 1 lt i lt i 1 5 Gi 1 5 22 Actually since modulos are pretty common too 151 also has a special notation for them and the set can therefore also be respresented as n gt i 1 lt lt and i 1 5 lt 2 It should be noted that always rounds down towards while integer division in C truncates towards 0 When translating conditions in C containing integer divisions and or modulos to isl constraints the user should therefore make sure the sign of the dividend is positive If not the integer division needs to be translated differently for positive and negative values for i 0 i lt n i T t i a i for i 0 i lt n i for j 0 j n i j F t j fCt jl tlj 1 for i 0 i lt n i B b i tlil Figure 1 5 A program with three loop nests Most programs involve more than one statement Although it is possible to work with different sets each representing the iteration domain of a statement it is usually more convenient to collect all iteration domains in a single set To be able to differen tiate between iterations of different statements with t
15. gt cat petr 46 1 1 1 1 10 11 1000 101 100 10010 1 38 0 3 n gt barvinok_enumerate series lt petr POLYHEDRON Dimension 4 Constraints 5 Equations 80 Rays 5 Lines 0 Constraints 5 6 Inequality 1 1 1 1 8 Inequality 1 1 0 0 98 Inequality 9 1 1 0 0 Inequality 0 0 1 1 Inequality 0 0 0 9 1 Rays 5 6 Ray 1 1 1 3 Ray 1 1 0 2 Ray 1 0 9 1 Ray 0 G9 9 1 Vertex 1 1 1 3 1 POLYHEDRON Dimension 1 Constraints 1 Equations 9 Rays 2 Lines 1 Constraints 1 3 Inequality 0 1 Rays 2 3 Line 1 Vertex 8 1 n 3 CC1 n 1 n 1 n72 1 n 3 Options floor f convert fractionals to floorings convert C convert fractionals to periodics series S compute rational generating function instead of piecewise step polynomial explicit e convert computed rational generating function to a piecewise step polynomial 3 3 barvinok enumerate e The program barvinok enumerate e enumerates a parametric projected set It takes a single polytope in PolyLib notation as input followed by two lines indicating the number or existential variables and the number of parameters and optionally followed by a list of parameter names The syntax for the line indicating the number of existen tial variables is the letter E followed by a space and the actual number For indicating the number of parameters the letter P is used The following example corresponds to Verdoola
16. of the other ujs since there exists a convex combination of and u on this hyperplane In particular since a and a have the same sign we have Qj uj t ucH for ojo gt 0 5 8 Qj t Qj t aj The corresponding cones K and K with Ko K therefore intersect in a common face F c H Let K pos tu m U w then any x K lies both in some cone with 0 lt i lt k and in some with k 1 lt i x d Just subtract an appropriate multiple of Equation 5 7 The cones ae and K E therefore both form a triangulation of and hence k d amp gt a Ej 5 9 i l je i k 1 jh or a K F 5 10 i l j with s 1forl lt i lt k s 1 fork 1 lt i lt d 6 1 1 and some lower dimensional faces Figure 5 11 shows the possible configurations in the case of a 3 dimensional cone As explained above there are several ways of avoiding the lower dimensional faces in 5 10 Here we will apply the following proposition 48 Proposition 5 12 and Verdoolaege 2008 Let 2 al Pil del Pi 0 5 13 ich be a finite linear identity of indicator functions of closed polyhedra P Q where the polyhedra with i I are full dimensional and those with i lower dimensional Let each closed polyhedron be given as P x b gt Bi for je Ji Let y Q be a vector such
17. that evaluates to the number of integer points in the dilation of P by a factor n The input is the same as that of barvinok count except that it may be followed by the variable name The functionality is the same as running barvinok enumerate on the cone over P placed at n 1 Options floor f convert fractionals to floorings convert c convert fractionals to periodics series s compute Ehrhart series instead of Ehrhart quasi polynomial 3 6 polyhedron sample The program polyhedron sample takes a polytope in PolyLib notation and prints an integer point in the polytope if there is one The point is computed using Polyhedron Sample 3 7 polytope scan The program polytope scan takes a polytope in PolyLib notation and prints a list of all integer points in the polytope Unless the direct options is given the order is based on the reduced basis computed with Polyhedron Reduced Basis Options direct d listthe points in the lexicographical order 3 8 lexmin The program lexmin implements an algorithm for performing PIP based on rational generating functions and provides an alternative for the technique Feautrier 1988 which is based on cutting planes The input is the same as that of the example program from piplib Feautrier 2006 except that the value for the big parameter needs to be 1 since there is no need for big parameters and it does not read any options from the input file 3 9 barvinok summate Gi
18. 1 1 lyjtys 1 3 0yiy3 0 Adding up all contributions in the final columns of these tables we obtain a grand total of 12 5 16 Summation through exponential substitution and Laurent ex pansions This section was inspired by Baldoni et al 2009 As in subsection 5 15 we want to compute n times the coefficient of the y term in the Taylor expansion of f eY However we do not write 5 46 i e e ediat sibi gle 7 1 ei as in 5 47 but as d 1 1 gSi DXb y Y _ iid Li les 82 instead The second term in each factor is of the form 1 e _ 1 1 xe x l e x xl e _1 b k AH ee T k 1 k ay REND g k 1 with b k u the Bernoulli polynomials The second line follows from 5 38 and the third line follows from b 0 1 We therefore have a 1 b k 1 s p e 11 255 Di tn j l k 0 b yn b 1 s p 1 mem J y ERIT VATIC IMP qi nj j l z i n p y 2 7 gt 0 n gt 0 where is the first non zero coefficient of b and the vector b contains the subse quent d f coefficients of b The expansion of the first term follows from 5 49 while the expansion of the second term follows from 5 48 Note however that unlike each factor is now a sum of two series instead of a product of two series In particular there are two kinds of non overlapping terms In the first kind of terms
19. 10 Off Qo The rows of span the orthogonal complement of the given subspace Since can be extended to a unimodular matrix these rows form an integral basis If the entries on the diagonal of H are all 1 then can be extended to a unimodular matrix by concatenating and The resulting matrix is unimodular since Mi H 0 Qi Q2 0 Ln mn m Q2 This method for extending a matrix of which only a few lines are known to uni modular matrix is more general than the method described by Bik 1996 which only considers extending a matrix given by a single row 0 5 8 Ensuring a polyhedron has only revlex positive rays The barvinok series with options function and all further gen fun manipula tions assume that the effective parameter domain has only revlex positive rays When used to compute rational generating functions the barvinok enumerate application will therefore transform the effective parameter domain of a problem if it has revlex negative rays It will then not compute the generating function f x Pon 7 pez but 80 WPry a n Zh e instead where p Tp t with T Z and t Z is an affine transformation that maps the transformed parameter space back to the original parameter space First assume that the parameter domain does not contain any lines and that there are no equalities in the description of Pp that force the values of p for which Pp contains integer
20. 7 dom 12 L5 domain domain 15 DomainConstraintSimplify DomainIncludes double description dual space eadd 28 29 Edmonds J 93 722 eequal 30 Ehrhart quasi polynomial Ehrhart series Ehrhart E Eisenbrand F Eisenschmidt E 87 722 emul 28 Enge A 27 enode 2224 27H29 Enumeration 23 eor esum 29 euler 42 Euler Maclaurin formula 129 local see local Euler Maclaurin for mula evalue 22 23 28130 33 evalue eval 29 evalue frac2floor 26 29 evalue sum 29 example exists 9 explicit representation exponential substitution 60 facet open see open facet Feautrier P Fern ndez F J 117 27 flooring 26 27 29 B0 BOH41 fractional 26130 BHI Fukuda K 92 118 727 fundamental parallelepiped Gale transform 87 Gallivan K 726 Garbervetsky D 017 727 Gawrilow E 44 722 Gelfand M 87 222 gen fun 31 32 56 gen fun add 32 gen fun coefficient 32 gen fun Hadamard product gen fun operator evalue gen fun substitute 32 generalized basis reduction generalized reduced basis 35 Girbal S 118 724 gist 8 GLPK 92 GMP 22 Gomory E 41 722 h vector H STAR VECTOR 44 h star vector 44 Haase C 018 720 Hadamard product 32 111 Hartmann M 27 Hartmann M E 93 722 Haws D 21 Hemmecke R 91 7271E723 H
21. In some situations we are given the generating function of some integer set and we would like to know if the set is infinite or not Typically we want to know if the set is empty or not but we cannot simply count the number of elements in the standard way since we may not have any guarantee that the set has only a finite number of elements We will consider the slightly more general case where we are given a rational generating function f x of the form such that 9x 5 80 seQnZ4 converges on some nonempty open subset of C Q is a pointed polyhedron and c s gt 0 and we want to compute S 5 81 seQnz4 where the sum may diverge i e S oo The following proposition shows that we can determine 5 in polynomial time For a sketch of an alternative technique see Woods 2005 Proof of Lemma 16 Proposition 5 82 Fix d and Given a rational generating function of the form with k k and a pointed polyhedron C 9 then there is a polynomial time algo rithm that determines for the corresponding function c s whether the sum 5 81 diverges and computes the value of S if it does not Proof Since Q is pointed the series converges a neighborhood of ef e eft for any such that rz lt 0 for any extremal ray of Q and such that bij for any bj in 5 26 Let and perform the substitution x The function g t 7 is then also a short rational gener
22. See Section 2 6 The following example was taken by Loechner 1999 from Loechner 1997 Chapter II 2 cat loechner Dimension of the matrix 7 Constraints ij cte 100000 0 lt i 100100 1 lt 010000 0 lt 1 10000 lt i 001000 0 lt 1 1 100 0 lt i j 1110 10 0 1 N V 2 parameters constraints 0 4 gt barvinok_enumerate lt loechner POLYHEDRON Dimension 5 Constraints 6 Equations 1 Rays 5 Lines 0 Constraints 6 7 Equality 1 1 1 0 1 0 Inequality 9 1 1 1 1 8 Inequality 9 1 0 0 6 98 Inequality 9 9 1 0 0 0 Inequality 89 2 2 9 1 0 Inequality 9 0 0 0 6 1 Rays 5 7 Ray 1 0 1 1 2 37 Ray 1 1 0 1 2 Vertex 0 0 0 0 0 1 1 Ray 1 0 1 0 0 1 1 POLYHEDRON Dimension 2 Constraints 1 Equations 9 Rays 3 Lines 2 Constraints 1 4 Inequality 0 0 1 Rays 3 4 Line 1 8 Line 0 1 Vertex 0 0 1 1 P Q gt 0 2P Q gt 0 1 gt 0 C 1 2 2 1 0 1 2 P C 3 8 072 C 1 2 1 2 0 0 1 4 Q C 5 4 1 72 0 0 1 gt 90 P Q 1 gt 0 1 gt 0 C 1 8 072 C 1 2 C 1 2 0 8 3 4 Q C 5 4 C 1 2 0 0 1 The output corresponds to 1P PO 3P 30 2 20 O 1 3 20 ifP lt Q lt 2P 10 3 5 20 2 20 if0 lt Q lt P 1 The following is an example of Petr Lison k
23. Analytical computation of Ehrhart polynomials and its applications for embedded systems Report CW 376 Depart ment of Computer Science K U Leuven Leuven Belgium URL http www cs kuleuven ac be publicaties rapporten cw CW376 abs html 1117 Verdoolaege S Seghir Beyls Loechner and Bruynooghe 20044 September Analytical computation of Ehrhart polynomials Enabling more compiler analyses and optimizations In Proceedings of International Confer ence on Compilers Architectures and Synthesis for Embedded Systems Wash ington D C pp 248 258 117 Verdoolaege S Beyls M Bruynooghe and Catthoor 20052 Experiences with enumeration of integer projections of parametric polytopes In R Bodik Ed Proceedings of 14th International Conference on Compiler Construc tion Edinburgh Scotland Volume 3443 of Lecture Notes in Computer Science Berlin pp 91 105 Springer Verlag Verdoolaege S K M Woods M Bruynooghe and R Cools 2005b Computation and manipulation of enumerators of integer projections of parametric polytopes Report CW 392 Dept of Computer Science K U Leuven Leuven Belgium 117 Verdoolaege S H Nikolov and T Stefanov 2006 March Improved derivation of process networks In 4th Workshop on Optimization for DSP and Embedded Systems ODES 4 118 126 Verdoolaege S H Nikolov and T Stefanov 2007a pn A tool for improved derivation of process networks EUR
24. bj 5 48 with the multinomial coefficients For the first factor we compute the Laurent expansion of 79 each of its factors 1 l Dive _ 1 b 1 Xi ra bin IYIN cl Sewer 4 4 n 140 b y 2 a JCD Em c 5 49 gt 0 if Of where bj is the first non zero coefficient of b and the vector b contains the subsequent d f coefficients of b Given a polynomial q y p Bm p y that we wish to sum over the integer points of a polytope P we perform the following operations for each unimodular cone in the decomposition of each vertex cone e For each m with Bm p 0 Compute all sums N 24 0 X nj 1 nj of exponents from such that N lt m and compute the corresponding coefficient yy in the prod uct of Laurent series by enumerating all combinations of n leading to the same N Note that there are only a finite number of N satisfying this con straint since gt d By reordering the variables such that the highest exponents occurs for the first variable the number of N can be reduced For each of these N Compute the coefficient g N p of y N in the product of Taylor ex pansions 5 48 e The contribution of this cone is the sum of m o By p YN m N D over all considered m and N Within each vertex cone computation the coefficients and g w p only need to be computed once Example 5 50 Consider
25. codegen m using s2 coefficients s 52 solutions 5 3 5 55 5 m3 cross s deltasm m deltas_map dom m dom q dom f domain m domain 4 domain f m domain map m M M M M M M s ranm s range ment simplify the representation of set s by trying to combine pairs of basic sets into a single basic set simplify the representation of map m by trying to combine pairs of basic maps into a single basic map simplify the representation of by trying to com bine pairs of basic sets in the domain of into a single basic set simplify the representation of f by trying to com bine pairs of basic sets in the domain of f into a single basic set generate code for the given domain generate code for the given schedule generate code for the domain s using the options rn generate code for the schedule m using the options nm The set of coefficients of valid constraints for 51 The set of elements satisfying the constraints with coefficients in Cartesian product of s and 52 Cartesian product of m and theset y x x yem the map x y y x xyem domain of map m domain of piecewise quasipolynomial 4 domain of piecewise quasipolynomial fold f domain of map m domain of piecewise quasipolynomial 4 domain of piecewise quasipolynomial fold f map from a wrapped copy of m to the domain of m range of map m range of map m continued on
26. jJ can be constructed as i j 0 lt i j lt 10 i i or as i j 0 lt i j lt 10 and not j Note that an occurrence of a relational operator in a set description may define several constraints one for each pair of arguments The elements in a list of arguments are separated by a comma If there are no constraints on the coordinates i e in case of a universe set the colon may be omitted as well For example 11 represents the entire unnamed zero dimensional space and should not be confused with which represents the empty set Returning to the iteration domain of the loop nest in Figure 1 1 we usually do not want to analyze such a program for a specific value of n but instead for all possible values of n at once A generic description of the iteration domain can be obtained through the introduction of a free parameter as in n gt i j 1 lt i lt and 1 lt j lt The optional parameters should be declared by placing them in a comma delimited list inside followed by an gt in front of the main set description The parameters are global and are identified by their names so the order inside the list is arbitrary This should be contrasted to the coordinates of a space the names of which are only relevant within the description of the set and which are instead identified by their positions That is n gt i j 1 lt i lt and 1 l
27. replace by 1 When used in the computation of the integer hull of some polytope it is useful to not only keep track of the best point so far but of all points found These points will all lie outside of the current approximation of the integer hull and adding them all instead of just one will typically get us to the complete integer hull quicker 94 12 3 4 5 6 7 8 9 10 11 12 13 14 15 16 tt e Figure 5 66 The integer points of a polytope projected on an objective function Example 5 67 Assume that the values of some objective function attained by the in teger points of some polytope are as shown in Figure 5 66 and assume we know that the optimal value lies between 1 and 16 In the first step we would look for a point attaining a value in the interval 1 16 Suppose this yields a point attaining the value 8 second line of the figure We record this point as the current best and update the search interval to 1 7 In the second step we look for a point attaining a value in the interval 1 4 but find nothing and set the search interval to 5 7 In the third step we consider the interval 5 7 and find a point attaining the value 6 We update the current best value and set the search interval to 5 5 In the fourth step we consider the interval 5 5 find no points and update the interval to 6 5 Since the lower bound is now larger than the upper bound the algorithm terminates returning the best
28. sum 16 supporting cone 69 Szarek S J 20 132 Tauzer J 21 Tawbi N 63 725 term 32 Todd polynomial 60 TOPCOM transverse cone type typeof 17 ub 17 under unimodular matrix unimodular transformation unimodular_complete union unrestricted sign unwrap upoly 16 upper bound 42 using 15 Value 22 29 B3 van Herick A W 118 720 726 Van Campenhout J 117 27 Van Engelen R A 63 65 726 Vasilache N 118 724 vec_QQ 30 vec ZZ 30 vector partition 87 Verdoolaege S AS SRE IRIS ON JAN gale um JBR 12 7 Vergne M 68 69 69 72 118 720 vertex 98 vertices 17 36 volume iov N Walsh 26 Walter M 723 width 92 Wilde D K 1241127 Wonnacott D 25 Woods M BA 54 104 105 112 17 723 25 27 wrap 17 write x n 23 x p 23 Yoshida R Zelevinsky A V 22 Zeron E S 118 724 zip 7 zsolve 77 30
29. value arr 0 0 value arr 1 9 If the size is 3 then the value is value arr 0 8 value arr 1 value arr 2 The type of arr 8 is typically fractional Finally the partition type is used to represent piecewise quasi polynomials We prefer to encode this information inside evalues themselves rather than using Enumerations since we want to perform the same kinds of operations on both quasi polynomials and piecewise quasi polynomials An enode of type partition may not be nested inside another enode The size of the array is twice the number of chambers Pointers to chambers are stored in the even slots whereas pointer to the associated quasi polynomials are stored in the odd slots To be able to store pointers to chambers the definition of evalue was changed as follows typedef struct _evalue Value d denominator union Value n numerator if denominator gt 0 struct enode p pointer if denominator 0 Polyhedron D domain if denominator 1 evalue Note that we allow a chamber to be a union of polyhedra as discussed in Verdoolaege 2005 Section 4 5 1 Chambers with extra variables i e those of Section 4 6 5 are only partially supported The field pos is set to the actual dimension i e the number of parameters 2 4 Operations on Quasi polynomials In this section we discuss some of the more important operations on evalues provided by the barvin
30. 1 define BV_APPROX_SCALE_NARROW2 1 lt lt 2 define BV_APPROX_SCALE_CHAMBER 1 lt lt 3 int scale_flags define BV_VOL_LIFT 0 define BV_VOL_VERTEX 1 define BV_VOL_BARYCENTER 2 int volume_triangulate basis reduction options define BV_GBR_NONE 0 define BV_GBR_GLPK 1 define BV_GBR_CDD 2 int gbr_lp_solver define BV_LP_POLYLIB 0 define BV_LP_GLPK 1 define BV_LP_CDD 2 define BV_LP_CDDF 3 int lp_solver define BV_HULL_GBR 0 define BV_HULL_HILBERT 1 int integer hull r struct barvinok_options barvinok_options_new_with_defaults The function barvinok_options_new_with_defaults can be used to create a barvinok_options structure with default values 25 e PolyLib options MaxRays The value of MaxRays is passed to various PolyLib functions and defines the maximum size of a table used in the double description computation in the PolyLib function Chernikova In earlier versions of PolyLib this parameter had to be conservatively set to a high number to ensure suc cessful operation resulting in significant memory overhead Our change to allow this table to grow dynamically is available in recent versions of PolyLib In these versions the value no longer indicates the maximal ta ble size but rather the size of the initial allocation This value may be set to 6 or left as set by barvinok options new with defaults e NTL options LLL a LLL b The values used for the reduction para
31. 1 lt j lt 52 n gt n 5 S1 for i 1 i lt n i 3 f S Figure 1 3 A loop with non unit stride for i 1 i lt n i if Ci 1 5 lt 2 JE S Figure 1 4 A loop with a modulo condition If a loop has a non unit stride as in then affine constraints on the coor dinates and the parameters are not sufficient to represent the iteration domain What is needed is a way to express that the value of i is equal to 1 plus 3 times some integer and this is where existentially quantified variables can be used Existentially quantified variables are introduced by the exists keyword and followed by a colon They can only be used within a single disjunct As an example the iteration domain of the loop in Figure 1 3 can be represented as n gt i exists a 1 lt i lt n andi 1 3a Existentially quantified variables are also useful to represent modulo constraints Consider for example the loop in Figure 1 4 The iterator values i for which the state ment S is executed lie between 1 and n and are such that the remainder of i 1 on division by 5 is less than or equal to 2 The constraint i 1 mod 5 lt 2 can be written as i 1 5 lt 2 where f is the greatest integer part of HL That is f is the unique integer value satisfying the constraints 5f lt i 1 and 5f gt i 1 4 The iteration domain of the statement in Figure 1 4 can therefore be represented as n
32. 2 5 1 4 5 T 3 1 2 3 1 3 4 2 3 5 3 4 5 1 2 43 2 4 5 88 We see that we have three chambers in the decomposition one with 8 vertices and two with 6 vertices Take the second vertex 1 2 3 of the first chamber This vertex corresponds to the saturation of the constraints j gt 0 k gt O and i p i e i j k p 0 0 Plugging in this vertex in the remaining constraints we see that it is only valid in case p gt 0 r gt 0 and 2p gt 0 For the remaining vertices of the first chamber we similarly find 0 1 2 0 0 0 p 0 q r gt O0andq gt 0 1 2 3 p 0 0 gt 0 gt 2 gt 0 0 1 4 0 p r qt r2 gt 0 p gt O0andq gt 0 1 3 4 p 0 gt 0 gt 2 gt 0 0 2 5 0 9 0 420 p20and q r20 2 3 5 p 2p q 0 gt 0 2 gt gt 0 9 4 5 0 4 9420 9 gt gt 0 3 4 5 p2p qr gt 0 2 gt gt 0 Combining these constraints with the initial constraints the problem on the parame ters p 0 q gt 0 we find the chamber p q r p OA ptr 0Aq 0 For the second chamber we have 1 2 3 p 0 0 pzOrzO0and2p q20 1 3 4 p 0 p 0 r gt 0and2p q20 2 3 5 p 2p 9 0 p202p qzO0andrz0 3 4 5 p 2p q r p20 2p q20andr20 1 2 5 4 0 0 q 2 0 2p q 2 Oand 2p q 2r gt 20 1 4 5 4 0 p 4 r q20 2p q 2r20and2p q420 The chamber is therefore
33. 3 Both arguments are assumed to correspond to indicator functions evalue esum evalue E int nvar evalue evalue sum evalue E int nvar unsigned MaxRays The function esum has been superseded by evalue sum The function evalue sum performs the summation operation from Verdoolaege 2005 Section 4 5 4 The piece wise step polynomial represented by E is summated over its first nvar variables Note that E must be zero or of type partition The function returns the result in a newly allocated evalue Note also that E needs to have been converted from fractionals to floorings using the function evalue frac2floor void evalue frac2floor evalue e This function also ensures that the arguments of the floorings are positive in the relevant chambers It currently assumes that the argument of each fractional in the original evalue has a minimum in the corresponding chamber double compute evalue const evalue e Value list args Value compute poly Enumeration en Value list args evalue evalue eval const evalue e Value values The functions compute evalue compute poly and evalue_eval evaluate a piece wise quasi polynomial at a certain point The argument list args points to an array of Values that is assumed to be long enough The double return value of compute evalue is inherited from PolyLib void print evalue FILE DST const evalue char pname The function print evalue dumps a human readable representat
34. 4 n 1 and 0 lt j n and 2 n 1 lt i j lt 4 n 2 and i lt 2 n 1 ub f u ub 0 u n gt 0 lt lt 18 Y m n gt i j gt i 1 j 1 1 lt i j lt n i j gt i 1 j 1 1 lt i lt n and 2 lt j lt m Gn 0 codegen gt A i gt 1 0 0 lt i lt N B i gt 1 1 1 lt i lt N Y k gt i 1 lt ic n gt 2 n n gt 1 m gt 1 lt c lt mj k gt max 3 k k 2 k gt 1 1 2 7 Comparison to other Interactive Polyhedral Tools Two related interactive polyhedral tools are the Omega calculator Kelly et al 1996 and SPPoC Boulet and Redon 2001 The syntax of iscc was very much inspired by that of the Omega calculator However the Omega calculator only knows sets and maps In particular it does not perform any form of counting An earlier version of barvinok came with a modified version of the Omega calculator that introduced an 13 operation for counting the number of elements in a set but it would simply print the result and not allow any further manipulations SPPoC does support counting but only the basic operation of counting the elements in a set In particular it does not support weighted counting nor the computation of upper bounds It also only supports single valued functions and not generic relations like the Omega calculator and iscc Inter nally SPPoC uses PolyL
35. 42 qi fi fixs S2 m si 42 qi s l f s compute schedule for the domains s that re spects all dependences in m and tries to minimize the dependences in compute a schedule for the domains in s that re spects all dependences in m and tries to minimize the dependences in and return the forest of bands in the schedule sum union union sum sum sum concatenation concatenation of stringification of o and S concatenation of and stringification of o difference set difference set difference difference product intersection intersection product product product product product intersect domain of m with s intersect domain of q with s intersect domain of f with s apply map m to set 51 apply q to each element in s and compute the sum of the results apply f to each element in s compute a bound of the results and return a list containing the bound and a boolean that is true if the bound is known to be tight continued on next page 18 Syntax Meaning m3 m no m m before ma mom m after BR qQ m 41 q2 m before qi q q m q after l m f m before f f m l f after m Im 5 gt 52 4 Qs q f s 53 1 S1 50 ma 96m m 5 42 41 h fi s join of m and m join of m and m join of m and m join of m and m join
36. Cones Poly Exp Cones Poly Exp hickerson 12 1 11625 9 24 8 90 7929 4 80 4 55 10 4251 432 4 19 803 0 66 0 62 100 980 1 42 1 35 84 0 13 0 12 200 550 1 00 0 92 76 0 12 0 12 300 474 0 93 0 86 58 0 12 0 10 500 410 0 90 0 83 42 0 10 0 10 1000 130 0 42 0 38 22 0 10 0 07 2000 10 0 10 0 10 22 0 10 0 09 5000 7 0142 0 11 7 0 12 0 10 hickerson 13 1 494836 489 463 483507 339 315 10 296151 325 309 55643 51 48 100 158929 203 192 9158 11 10 200 138296 184 173 6150 9 8 300 110438 168 157 4674 8 7 500 102403 163 151 3381 8 7 1000 83421 163 149 2490 8 7 2000 77055 170 153 1857 10 8 1 1 5000 57265 246 211 1488 13 1 10000 50963 319 269 1011 26 2 hickerson 14 1 1682743 2171 2064 552065 508 475 10 1027619 1453 1385 49632 62 59 100 455474 768 730 8470 14 13 200 406491 699 661 5554 11 10 300 328340 627 590 4332 11 9 500 303566 605 565 3464 11 9 1000 232626 581 532 2384 12 10 2000 195368 607 545 1792 14 12 5000 147496 785 682 1276 19 16 10000 128372 966 824 956 29 23 Table 2 Timing results of dual and primal decomposition with polynomial or expo nential substitution on the Hickerson examples 64 u amp gt Q X ag amp UR IR 1 5 32 1 I For the higher degree monomials we use the formula 1 1 1 b3 i Je m 5 33 k 0 id c with B the Bernoulli numbers which can be computed using the recurrence e m 1 2j o Bj0 1 5 34 po 7 Note that 5 33 is also
37. P3 2 8 ib x 2 x b ajb x ifa gt b i 0 x 0 x 0 effectively avoiding splintering if a given monomial contains a single integer part ex pression with argument of the form x b An argument of the form x c b can be handled through a variable substitution If the argument is of the form cx b with c 1 then Sakellariou 1996 3 27 proposes to rewrite the monomial as Ye mod Y mod b 0 0 mod b ye mod b x 0 y 0 y 0 67 and applying 5 35 However such an application results in an expression containing cx mod b y 0 which in turn leads to a polynomial of degree n 1 in mod D i e of degree one higher than the original expression Furthermore if the bound on x is rational then a itself contains a floor which on application of 5 35 results in a nested floor expression blocking the application of the same rule for the next variable Finally the case where a monomial contains multiple floor expressions either occurring in the input quasi polynomial or introduced by different variables having a rational bound with a non zero coefficient in the same variable is not handled Also note that if we disallow nested floor expressions then this rule will rarely be applicable since we try to eliminate variables with integer bounds first 5 14 Summation using local Euler Maclaurin formula In this section we provide some implementation details on using local Eu
38. T T v 9 k w v u U 1 i l j where u are the columns of A and Z ranges over 0 lt kj lt sj 45 Figure 5 4 The integer points in the fundamental parallelepiped of K Proof Since 0 lt x lt 1 it is clear that each such w lies inside the fundamental parallelepiped Furthermore v v547 A v vA aw kW vA A ez 71 4 ezixd 74 Finally if two such w equal 1 w Wo then w 5 K W p V SW with p ZZ k mods ie ky Since detS det Al we obtain all points in the fundamental parallelepiped by taking all Z4 satisfying 0 lt k lt o If the cone K is not closed then the coefficients of the open rays should be in 0 1 rather than in 0 1 In 5 3 we therefore need to replace the fractional part x x Lx by x x for the open rays Example 5 5 Let K be the cone 46 shown in Figure 5 4 Then _ 2 1 aim 1 2 a al le ol We have det A det S 2 and o o and 0 Therefore sfe off SLE EE w e up are E s 1 2 uallo qen o 5 3 Barvinok s decomposition of simple cones in primal space As described by De Loera et al 2004 the first implementation of Barvinok s count ing algorithm applied Barvinok s decomposition in the dual space Brion s polarizati
39. W 1991 Regular triangulations of convex polytopes Applied Geometry and Discrete Mathematics The Victor Klee Festschrift 4 443 456 58 Lepelley D A Louichi and H Smaoui 2008 April On Ehrhart polynomials and probability calculations in voting theory Social Choice and Welfare 30 3 363 383 118 Loechner V 1997 Contribution l tude des poly dres param tr s et applica tions en parall lisation automatique Ph D thesis University Louis Pasteur Strasbourg 137 Loechner V 1999 March Polylib A library for manipulating parameterized polyhedra Technical report ICPS Universit Louis Pasteur de Strasbourg France Loechner V and D K Wilde 1997 December Parameterized polyhedra and their vertices International Journal of Parallel Programming 25 6 525 549 22 Makhorin A 2006 July Gnu linear programming kit reference manual version 4 11 92 Meister B 2004 December Stating and Manipulating Periodicity in the Polytope Model Applications to Program Analysis and Optimization Ph D thesis ICPS Universit Louis Pasteur de Strasbourg France Meister B and S Verdoolaege 2008 April Polynomial approximations in the polytope model Bringing the power of quasi polynomials to the masses In J Sankaran and T Vander Aa Eds Digest of the 6th Workshop on Optimization for DSP and Embedded Systems ODES 6 117 Pfeifle J and J Rambau 2003 Computing triangulati
40. a statement in a loop nest Take for example the loop nest in and assume first that n has a fixed value say 5 The pairs of values of i and j for which statement S is executed are shown graphically in Figure 1 2 Mathematically this set of pairs can be represented as G eZ l1 lt i lt 5Al lt j lt ij and the isl notation is very similar i i j 153 55d 14313 In this notation the coordinates are separated by commas and enclosed in square brack ets This description of the space in which the set lives is followed by a colon and the constraints on the coordinates Assuming the iterators are incremented by one in every iterations of the loop a lower and upper bound on each loop iterator can be read off from the initialization and the test Note that in an iscc set the coordinates are as sumed to be integer by default For an iteration domain to be representable by such a set the iterators therefore need to be integers The constraints of a set need to be affine i e linear plus constant term These affine constraint may be combined through conjunctions and disjunctions or projections exists and negations not Note that the formula is immediately converted into Disjunctive Normal Form DNF so it may sometimes be more efficient to construct a set from smaller sets by applying basic operations such as intersection union and difference For example the following square with its diagonal removed G 5l0sijsl0 G
41. are two new facets but neither of them yields any further points We have therefore found the integer hull of the polytope 5 21 2 Optimization over the integer points of a polyhedron We assume that we want to find the minimum of some linear objective function When used in the computation of the integer hull of some polytope the objective function will therefore correspond to the inner normal of some facet During our search for an optimal integer point with respect to some objective func tion we will keep track of the best point so far as well as a lower bound and an upper bound such that the value at the optimal point if it is better than the current best lies between those two bounds Initially there is no best point yet and values for and u may be obtained from optimization over the linear relaxation When used in the com putation of the integer hull of some polytope the upper bound is one less than the value attained on the given facet of the current approximation As long as l lt u we perform the following steps e use the integer feasibility technique of subsection 5 20 to test whether there is any integer point with value in where u is uif the previous test for an integer point did not produce a point l if the previous test for an integer point did produce a point e if a point is found then remember it as the current best and replace u by the value at this point minus one e otherwise
42. automatique de poly dres pour la compilation Technique et Science Informatiques 20 8 1019 1048 13 Brion M 1988 Points entiers dans les poly dres convexes Annales Scientifiques de l cole Normale Sup rieure Quatri me S rie 21 4 653 663 47 Biieler A Enge Fukuda 2000 Exact volume computation for poly topes A practical study DMV Seminar Band 29 158 Clauss Loechner 1998 July Parametric analysis of polyhedral iteration spaces Journal of VLSI Signal Processing 19 2 179 194 22 Clauss P J Fernandez D Garbervetsky S Verdoolaege 2006 October Symbolic polynomial maximization over convex sets and its application to mem ory requirement estimation ICPS Research Report 06 04 Universit Louis Pas teur 117 Cohen J Hickey 1979 Two algorithms for determining volumes of convex polyhedra J ACM 26 3 401 414 58 Cook W M Hartmann R Kannan and C McDiarmid 1992 On integer points in polyhedra Combinatorica 12 1 27 37 93 Cook W Rutherford E Scarf D Shallcross 1993 An implemen tation of the generalized basis reduction algorithm for integer programming ORSA Journal on Computing 5 2 De Loera J A 1995 May Triangulations of Polytopes and Computational Alge bra Ph D thesis Cornell University 158 De Loera J A D Haws Hemmecke Huggins J Tauzer Yoshi
43. be cut by the affine hull of the facet The closedness of these facets is therefore predetermined by K For the other facets a choice will have to be made To be able to make the choice based on local information and without computing an explicit vector y we use the following convention We first assign an arbitrary total order to the rays If the affine hull of a facet separates the two rays not on the facet u and uj i e aja gt 0 5 8 then we choose y to lie on the side of the smallest ray according to the chosen order That is fij y gt 0 for the normal of the facet 49 pointing towards this smallest ray Otherwise i e if aja lt 0 the interior of will lie on one side of the facet and then we choose y to lie on the other side That is fiij y gt 0 for fi the normal of the facet pointing away from the cone Figure shows some example decompositions with an explicitly marked y To see that there is y satisfying the above constraints we need to show that RN S is non empty with 5 y fi j y gt 0 for all It will be easier to show this set is non empty when the u form an orthogonal basis Applying a non singular linear transformation T does not change the decomposition of w in terms of the u i e the o remain unchanged nor does this change any of the scalar products in the constraints that define RNS the normals are transformed by r3 Finding a vector y T RNS ensures that T
44. be found by combining say v p V2 p while the second direction can be found by combining say vy p and va p The parameter domain for the parametric polytope P p is 5 lt lt 5 The two closed chambers are therefore Ci peQ 5 ps lt 5 p C peQ 5 ps5 p To obtain a partition gives the internal point 0 0 which happens to meet the facets gt 0 and p gt 0 We therefore keep the facet with positive inner normal closed and open up the other facet The result is Ci p 0 lt p lt 5 5 lt lt 0 Since we are usually only interested in integer parameter values the latter chamber would become C 5 lt p lt 1 103 Our description differs slightly from that of of Eisenbrand and Shmonin 2007 First they consider all pairs of basic solutions instead of all pairs of vertices which means that they may consider basic solutions that are never feasible and that in case of a non simple polytope they may consider the same parametric vertex more than once The set of integer directions for a pair of vertices is the intersection of the sets of integer directions they obtain for each of the corresponding basic solutions Second they use a different method of creating a partition of partially open chambers which may lead to some lower dimensional chambers surviving and hence to a larger total number of chambers 5 24 Testing whether a set has an infinite number of points
45. c and the corre sponding widths we p we keep only a single direction of all those that yield the same width as an affine function of the parameters Then we construct the chambers where each of the widths is minimal 1 Ci ipe Q YJ Wa P lt wc p Note that many of the C may be empty or of lower dimension than Q and that the other C will intersect in common facets To obtain a partition of partially open full dimensional chambers we proceed as in subsection 5 4 99 Figure 5 75 The cone of directions C Example 5 74 Consider the non parametric polytope 3 50 gt 0 4 1 5x2 gt 0 xX 2x 3 gt 0 3x 40 3 gt 0 shown in Figure 5 75 The polytope has four vertices P 4x v 9 6 5 4 0 0 v4 5 3 The corresponding cones of directions for the given vertex to attain the minimum also shown in are pos 3 4 1 2 pos 4 5 1 2 pos 4 5 3 5 C pos 3 5 3 4 MM ee ee EN Let us now consider the directions in which v is minimal while v is maximal We find pos 4 5 3 4 V 0 100 Figure 5 76 The cone of directions Figure 5 77 The cone of directions C4 101 as shown in Figure 5 75 The vertices of the integer hull of C a
46. careful here since we allow the polyhedron to be unbounded and so we should apply the technique of with Q the projection of P on the remaining coordinates Note that testing whether we can stop is more expensive than applying the next iter ation since we have an extra 1 z factor in the denominator of one of the operands However we may save many iterations by stopping early and we will not needlessly replace a given representation of f x 2 by a more complex representation with more factors in the denominator An alternative way of checking whether we can stop is to test whether f7 x 2 f7 x y z as rational functions do so we would need to check that both REYD HYD fia e y 2 and 2 d z fa y 3 are zero and this Hadamard product is more expensive than the one above Example 5 91 Consider once more the parametric polytope 2xi t pt 520 2x p 520 20 5 gt 0 2x2 5 gt 0 from Woods 2004 Example 2 1 7 Example and assume we want to compute 3x eZ p x e PI 4 x That is we simply want to know for which values of p the polytope is non empty Now there are more efficient ways of answering this particular question but we will use it here as an example of the above algorithm The polytope P p is shown in Figure 5 90 for all integer value of the parameter for which the polytope is non empty 112
47. differences between image and domain deltas join in verse 1 and transitive closure The latter may result in an overapproximation The card operation computes the number of elements in a set or the number of elements in the image of a domain element of a map The operation is performed exactly and completely symbolically and the result is a piecewise quasipolynomial i e a subdivision of one or more spaces with a quasipolynomial associated to each cell in the subdivision As a trivial example the result of card A i gt B j 0 lt j lt is A i gt 1 1 i gt 0 Operations on piecewise quasipolynomials in clude sum and difference and the computation of an upper bound over the do main If the domain contains a pair of nested spaces then the upper bound is computed over the nested range As another trivial example the result of ube i gt j gt 72 1 lt lt 1 is i gt max i 2 i gt 8 True The first element in this list is the actual bound in the form of a piecewise quasipolynomial fold i e a maximum of quasipolynomials defined over cells The second indicates whether the bound is tight i e whether a maximum has been computed 1 2 5 Advanced Operations While the basic card operation simply counts the number of elements in an affine set it is also possible to assign a weight to each element of the set and to compute the sum of those weights over a
48. either uj or N uj and so the coefficient of one of these basis vectors will be strictly negative That is the sum of the normals will never be zero and the set R N S is non empty For each ray of cone Kj i e the cone with u replaced by w we now need to determine whether the facet not containing this ray is closed or not We denote the inward normal of this cone by n j Note that cone if it appears in 5 9 i e a 0 has the same facet opposite u and its normal nj will be equal to either n or nj depending on whether we are dealing with an external facet i e a facet of or an internal facet If on the other hand a 0 then nj noj If nij y gt 0 then the facet is closed Otherwise it is open It follows that the two or more occurrences of external facets are either all open or all closed while for internal facets exactly one is closed First consider the facet not containing uo w If o gt 0 then u and w are on the same side of the facet and so nj Otherwise Second if a 0 then replacing u by w does not change the affine hull of the facet and so nj Now consider the case that aja lt 0 i e u and u are on the same side of the hyperplane through the other rays If we project u u and w onto a plane orthogonal to the ridge through the other rays then the possible locations of w with respect to u and 50
49. evalue is a rational number with numerator x n and denominator d An enode is either a polynomial or a periodic depending on the value of type The length of the array arr is stored in size For a polynomial arr contains the coefficients For a periodic it contains the values for the different residue classes modulo the parameter indicated by pos For a polynomial pos refers to the variable of the polynomial The value of pos is 1 for the first parameter That is if the value of pos is 1 and the first parameter is p and if the length of the array is then in case it is a polynomial the enode represents 9 arr 1 p arr 2 p arr 1 1 p If it is a periodic then it represents 9 arr 1 arr 2 arr 1 1 Note that the elements of a periodic may themselves be other periodics or even polynomials In our library we only allow the elements of a periodic to be other periodics or rational numbers The chambers and their corresponding quasi polynomial are stored in Enumeration structures typedef struct enumeration Polyhedron ValidityDomain constraints on the parameters evalue EP dimension combined space struct enumeration next Ehrhart Polynomial corresponding to parameter values inside the domain ValidityDomain above f Enumeration For more information on these structures we refer to Loechner 1999 23 enode
50. integer hull Ph D thesis Ithaca NY USA 93 Hemmecke R 2002 On the computation of hilbert bases of cones World Scien tific 122 Hemmecke R Hemmecke M K ppe Malkin and M Walter 4ti2 a soft ware package for algebraic geometric and combinatorial problems on linear spaces Available at 91 Henrici P 1974 Applied and Computational Complex Analysis Pure and applied mathematics New York Wiley Interscience John Wiley amp Sons Volume 1 Power series integration conformal mapping location of zeros Pure and Applied Mathematics 59 Hubler S L and G Craciun 2012 Counting chemical compositions using ehrhart quasi polynomials Journal of Mathematical Chemistry 1 25 119 Huggins P 2006 iB4e A software framework for parametrizing specialized LP problems In A Iglesias and N Takayama Eds ICMS 2006 Proceedings of the Second International Congress on Mathematical Software Volume 4151 of Lecture Notes in Computer Science pp 245 247 Springer 93 Hung S and W Rom 1990 October An application of the hermite nor mal form in integer programming Linear Algebra and its Applications 140 15 163 179 87 Kannan R 1992 Lattice translates of a polytope and the Frobenius problem Com binatorica 12 2 161 177 Kelly W V Maslov W Pugh E Rosser T Shpeisman and D Wonnacott 1996 November The Omega calculator and library Technical
51. isl union map The union of the cells in the resulting isl pw qpolynomial is equal to the domain of the input isl map isl give isl pw qpolynomial isl pw qpolynomial sum isl take isl pw gpolynomial pwgp __isl_give isl union pw qpolynomial isl union pw gpolynomial sum isl take isl union pw qpolynomial upwqp Compute the sum of the given piecewise quasipolynomial over all integer points in the domain The result is a piecewise quasipolynomial that only involves the parameters If however the domain of the piecewise quasipolynomial wraps a relation then the sum is computed over all integer points in the range of that relation and the domain of the relation becomes the domain of the result __isl_give isl pw qpolynomial isl set apply pw qpolynomial isl take isl set set isl take isl pw qpolynomial pwqp __isl_give isl union pw qpolynomial isl union set apply union pw gpolynomial __151 take isl union set uset isl take isl union pw qpolynomial upwqp Compute the sum of the given piecewise quasipolynomial over all integer points in the intersection of the domain and the given set for i 1 i lt n i for j 1 j lt i j J S Figure 1 1 A loop nest Figure 1 2 The iteration domain of the loop nest in Figure 1 1 __isl_give isl pw gpolynomial isl map apply pw gpolynomial isl take isl map map isl take isl pw qpolynomial pwqp
52. looking for a point the moment curve oO 5 25 Enumerating integer projections of parametric polytopes In this section we are interested in computing c s t Z du eZ st u e P 5 83 with P c Q x Q4 x Q a rational pointed polyhedron such that P t u Q x Q I s t u PI is a polytope for any s This is equivalent to computing the number of points in the integer projection of a parametric polytope c s Z with Q7 x Q Q defined by t Exponential methods for computing are described by Verdoolaege et al 20053 and Seghir and Loechner 2006 Here we provide some implementation details for the polynomial method of Barvinok and Woods 2003 Theorem 1 7 for computing the generating function Y c s X which can then be converted into an explicit function c s Verdoolaege and Woods 2008 Corollary 1 11 Note that in contrast Bassin a WEN Themen L7 we may allow P to be an unbounded but still pointed polyhedron here as long as P is bounded since we replace their application of Kannan 1992 Lemma 3 1 by Eisenbrand and Shmonin 2007 Theorem 5 If there is only one existentially quantified variable m 1 then computing 5 83 is easy You simply shift P by 1 in the direction and subtract this shifted copy from the original D P P See e g Barvinok and Woods 2003 Figure 1 page 973 or Verdoolaege 2005 Figu
53. next page 15 Syntax Meaning m range map m identity s q lattice width s 5 lexmin s m lexmin 5 lexmax s m lexmax m 5 lift s 52 params s s params l parse file 5 S2 poly s m poly m q2 poly q q2 lpoly q q2 upoly q l pow m print o o read filename 52 sample s m sample m 5 Scan 5 m scan m source filename q2 sum qi a map from a wrapped copy of m to the range of m identity relation on s lattice width of the set s lexicographically minimal element of s lexicographically minimal image element lexicographically maximal element of s lexicographically maximal image element lift s to a space with extra dimensions correspond ing to the existentially quantified variables in s such that domain unwrap lift S is equal to S parameter domain of set s parameter domain of map m parse the file names and return a list consisting of the iteration domain the must write access rela tion the may write access relation the read access relation and the original schedule This operation is only available if pet support was compiled in polyhedral hull of s polyhedral hull of m polynomial approximation of 41 polynomial underapproximation of q polynomial overapproximation of q compute an overapproximation of the power of m and return a list containing the overapproximation and a boolean
54. of 0 Another alternative would be to reduce case B to x lt 1 and u x gt 1 while extending cases 4 and 5 to 1 lt 10 lt 1 and gt 1 u X 1 and lt l respectively with the remaining cases 1 lt l lt lt 1 having value 0 There does not appear to be a consistently better choice here as each of these three approaches seems to yield better results on some examples The last approach has the additional drawback that we would also have to deal with 5 cases even if the bounds are integers If at least one of the lower or upper bound is an integer affine expression then we can reduce the 3 or 5 cases to a single case case B by an affine substitution that ensure that the new lower or upper bound is zero In particular if X is an integer affine expression then we replace x by x and similarly for an upper bound 5 13 Exact Enumeration using Nested Sums The exact enumeration using nested sums proceeds in much the same way as the ap proximate enumeration from subsection 5 12 with the notable exception that we need to take the greatest or least integer part of any fractional bounds that may occur This has several consequences discussed below Since we will introduce floors during the recursive application of the procedure we may as well allow the weight p x in to be a piecewise quasipolynomial For the constant term becomes u amp 2 a ao X u K amp 1
55. once more the rational generating function 2 2 x 25 1 80 from Verdoolaege 2005 Example 39 and Example Assume we want to compute 2 XX yeroz2 We will need the following Bernoulli polynomials b 0 s 1 bd s 71 25 1 65 6s b 3 s l s 3 2 25 s 3 1 305 605 305 For the first term we have 2 0 2 1 0 0 1 1 and substitution yields 1 1 je 290 yi 2 e cxi hy l e 1 e t 17 1 seerd y i y T Ag 2 y T 2 14 Fe y Ge yi t yy Ce ee y ity We obtain the following results m ywyh m N m laBmynom_N 721 721 2 0 2 ly 4 0 Ly 2 0 2 0 4 0 740 120 179 179 2 2 ly 2 2 242 L 0 2 2 0 2 2 720212 360 3 1 3 1 211 211 GD p 60 ar 721 4 2 Iy 4 eo 81 For the second term we similarly obtain m m N 6mny N la 1 2 0 L D Dry GD E 1 1 2 0 1y 40 yt 2 0 1y CO aen 211 211 0 2 L D lyytyg 3 5p 179 179 2 ly 2 22 2 0 Np 2 2 7207 360 1 1 3 1 3 1 a 1 1 4 2 1ly y 4 gt Y 0 729 360 Finally for the third term we obtain m m N 6m ny m aBmynom_N 2 0 1 1 Lb 3 1 0 0 2
56. or all point s found 5 22 Computing the integer hull of a truncated cone In subsection 5 23 we will need to compute the integer hull of a cone with the origin removed C 0 5 22 1 Using the Hilbert basis of the cone As proposed by K ppe 2007 way of computing this integer hull is to first com pute the Hilbert basis of C see subsection 5 19 and to then remove from that Hilbert basis the points that are not vertices of the integer hull of C 0 The Hilbert basis of C is the minimal set of points b N 7 such that every integer point x C n Z4 can be written as a non negative integer combination of the b The vertices v of the integer hull of C V 0 are such that every integer point x C N Z4 0 can be written as s non negative rational combination of v Clearly any v is also a b since v can not be written as the sum of a rational convex combination of other integer points in C n Z V 0 and a non negative combination of the extremal rays r of C A fortiori it can therefore not be written as an integer combination of other integer points in C To obtain the v from the b we therefore simply need to remove first 0 0 and then those b that are not an extremal ray and that can be written as a combination bj f wihao gt Oand 1 k je 95 Figure 5 68 The Hilbert basis and the integer hull of a truncated cone Since the rg are al
57. orbarvinok series with options enumerates para metric polytopes in the form of a rational generating function The polyhedron P is assumed to have only revlex positive rays The function barvinok enumerate e series computes a generating function for the number of point in the parametric set defined by P with exist existentially quantified variables using the projection theorem as explained in subsection 5 25 The function barvinok enumerate scarf series computes a generating function for the number of point in the parametric set defined by P with exist existentially quantified variables which is assumed to be 2 This function implements the technique of Scarf and Woods stricted to problems with 3 or 4 constraints involving the existentially quantified vari ables 34 2 7 Auxiliary Functions In this section we briefly mention some auxiliary functions available in the barvinok library void Polyhedron_Polarize Polyhedron P The function Polyhedron Polarize polarizes its argument and is explained in Ver doolaege 2005 Section 4 4 2 int unimodular complete Matrix M int row The function unimodular complete extends the first row rows of M with an integral basis of the orthogonal complement as explained in Section 5 7 Returns non zero if the resulting matrix is unimodular int DomainIncludes Polyhedron D1 Polyhedron D2 The function DomainIncludes extends the function PolyhedronIncludes provided by
58. p q r 19 gt gt 0 Note that by intersecting with the initial constraints this chamber is no longer full dimensional and can therefore be discarded Finally for the third chamber we have 1 2 3 p 0 0 p2O0rz0and2p q20 1 3 4 p 0 pzOrzOand2p q20 2 3 5 p 2p 4 0 gt 0 2 gt gt 0 3 4 5 2 q r gt 0 2 gt gt 0 1 2 4 p r 0 0 p r2 0 r gt 0and2p q 2r 0 2 4 5 p n2p q 2n0 p rz0 2p q 2rzO0andrz 0 The chamber is therefore p q r p r Z0 qzZO0 rz0 j Now let us consider general parametric polytopes First note that we can follow the same procedure as above if we replace x by x c p in 5 60 i e if our problem has the form x gt Ax x b p Ac p as saturating a constraint x 0 is equivalent to saturating the constraint x c p and similarly aj x bj p is equivalent to aj x aj In the general case the problem has the form 5 62 Ax gt b p and then we apply the technique of subsection 5 17 Let A be a non singular square submatrix of A with the same number of columns and compute the left HNF H 89 A U with U unimodular and H lower triangular with non positive elements below the diagonal Replacing x by Ux we obtain Hx gt b p A U x lt b p with A the remaining rows of A and b p split in the same way If H happens to be the identit
59. quan tified variables such that this direction lies in the direction of one of the transformed variables Then the remaining existentially quantified variables are projected out by applying the technique recursively The resulting generating function will have gaps at most w m wide and so we have to subtract at most w m shifted copies of this generating function before we can plug in 1 to project out the selected and now only remaining existentially quantified variable We now look at each of these step in a bit more detail We are given a polyhedron P such that P is a polytope and we want to compute a 109 generating function f x y for the set T in 5 85 We first compute the lattice width directions of the m dimensional parametric polytope as in subsection 5 23 The result is a partition of the parameter domain i e the projection of P onto the first n d coordinates into partially open polyhedra Q together with the lattice width direction corresponding to each Q Since the generating functions only encode integer points we can replace each open facet a x b gt 0 by the closed facet a x b 1 gt 0 to obtain a collection of closed polyhedra Now let P PAO x Q and let f x be the generating function of the set T s t Ju Z s t u Pi Then clearly fa J fios y From now on we will consider a particular P with corresponding lattice width c and drop the i subscrip
60. report University of Maryland 13 Koeppe M M Queyranne and C T Ryan 2010 Parametric integer programming algorithm for bilevel mixed integer programs Journal of Optimization Theory and Applications 146 1 137 150 119 K ppe M 2006 LattE macchiato version 1 2 mk 0 7 1 an improved ver sion of De Loera et al s LattE program for counting integer points in polyhedra with variants of Barvinok s algorithm Available from URL http www math uni magdeburg de mkoeppe latte K ppe M 2007 A primal Barvinok algorithm based on irrational decompositions SIAM Journal on Discrete Mathematics 21 1 220 236 K ppe M 2007 June personal communication 95 K ppe M and S Verdoolaege 2008 Computing parametric rational generating functions with a primal Barvinok algorithm The Electronic Journal of Combi natorics 15 R16 K ppe M S Verdoolaege and K M Woods 2008 July An implementation of the barvinok woods integer projection algorithm In M Beck and T Stoll Eds The 2008 International Conference on Information Theory and Statisti cal Learning 117 123 Lagarias J C W Lenstra Jr and C P Schnorr 1990 Korkin zolotarev bases and successive minima of a lattice and its reciprocal lattice Combina torica 10 4 333 348 Lasserre J B and E S Zeron 2005 An alternative algorithm for counting lattice points in a convex polytope Math Oper Res 30 118 Lee C
61. that b y 0forallie Ij UD j Jj Foreachi we define the half open polyhedron x e Q bb gt Bi for j J with bj y gt 0 1 5 14 gt Bi for j Ji with bj y lt 0 Then a el 0 5 15 ich When applying this proposition to 5 10 we obtain x 2 5 Ki 5 16 where we start out from a given K which may be K itself i e a fully closed cone or the result of a previous application of the proposition either through a triangulation Section 5 4 or a previous decomposition In either case a suitable y is available either as an interior point of the cone or as the vector used in the previous application which may require a slight perturbation if it happens to lie on one of the new facets of the cones K We are however free to construct a new y on each application of the proposition In fact we will not even construct such a vector explicitly but rather apply a set of rules that is equivalent to a valid choice of y Below we will present an intuitive motivation for these rules For a more algebraic shorter and arguably simpler motivation we refer to K ppe and Verdoolaege 2008 The vector y has to satisfy b7 y gt 0 for normals b of closed facets and b7 y lt 0 for normals b of open facets of K These constraints delineate a non empty open cone R from which y should be selected For some of the new facets of the cones K j the cone R will not
62. that is true if the overapproximation is known to be exact print object read object from file a sample element of the set 5 a sample element of the map m the set s split into individual elements provided s contains a finite number of elements the map m split into individual elements provided m contains a finite number of elements read commands from file sum q over all integer points in the domain of or if the domain of q wraps a map over all integer points in the range of the wrapped map continued on next page 16 Syntax Meaning S typeof o l ubq vertices 5 s wrap m unwrap s write filename o mM zip m m4 m before m under m3 l last before m under m3 ms any last m before m3 under m4 a string representation of the type of compute an upper bound on the piecewise quasipolynomial q over all integer points in the do main of q and return a list containing the upper bound and a boolean that is true if the upper bound is known to be tight If the domain of q wraps a map then the upper bound is computed over all integer points in the range of the wrapped map instead a list of vertices of the rational polytope defined by the same constraints as s wrap the map in a set unwrap the map from the set write object to file the cross product of domain and range of i e w2yoa xo2zn woxo onem compute a map from any element x in the domain o
63. valid if m 0 i e 0 0 a fact that can be easily shown using Newton series Van Engelen et al 2006 Since we can only directly apply the summation formula when the lower bound is zero or one we need to consider several cases 1 21 u amp u amp I 1 T 9 xy 1 xy 1 n S UR 1 1 2 u 1 u amp K 1 mAN D S IR 1 3 1 lt O and u gt 0 u amp 1 2 exi zp AX CD x1 1 R x1 0 xy 1 as 8 5 UR 1 1 1 D If the coefficients in the lower and upper bound are all integer then the above 3 cases partition the integer points in the projection of P onto the remaining variables However if some of the coefficients are rational then the lower and upper bound can lie in the open interval 0 1 for some values of We may therefore also want to consider the following two cases 65 4 0 lt I lt 1 u amp 2 06x a6 U 1 x 1 5 0 lt lt 1 u amp 2 xy a 1 5 OIR 1 I Note that we may add the constraint u gt 1 to the constraint lt 1 to 5 since the correct value for these two cases would be zero if these extra constraints do not hold An alternative to adding the above two cases would be to simply ignore them i e assume a value
64. x1 1 66 Since we force the lower and upper bounds to be integers cases 4 and 5 do not occur while the conditions for cases I and 2 can be simplified to gt 0 and u amp lt 0 respectively If the variable x appears in any floor expression either because such an expression was present in the original weight function or because it was introduced when another variable with an affine bound in x was summed then the domain has to be splintered into D parts where D is the least common multiple of the denominators of the coeffi cients of x in any of the integer parts In particular the domain is split into x Dy i for each i in 0 D 1 Since D is proportional to the coefficients in the constraints it is exponential in the input size This splintering will therefore introduce exponential behavior even if the dimension is fixed Splintering is clearly the most expensive step in the algorithm so we want to avoid this step as much as possible already noted that summation should pro ceed over variables with integer bounds first This can be extended to choosing a vari able with the smallest coefficient in absolute value In this way we can avoid splinter ing on the largest denominator Sakellariou 1996 claims that splintering can be avoided altogether In particular Sakellariou 1996 Lemma 3 2 shows that 2 x x mod b x 0 with a and b integers is equal to 2 ifa lt b 1 bot T 5 35
65. 00 and is used by Rabl 2006 in his implementation The Cohen Hickey trian gulation is a much more efficient variation and uses one of the vertices instead of the barycenter The facets incident on the ver tex do not have to be considered in this case because the resulting subpolytopes would have zero volume Another possibility is to use a lifting triangulation De Loera 1995 In this triangulation each vertex is assigned a random height in an extra dimension The projection of the lower envelope of the resulting polytope onto the original space results in a subdivision which is a triangulation with very high probability A complication with the lifting triangulation is that the constraint system of the lifted polytope will in general not be linearly parameterized even if the original poly tope is It is however sufficient to perform the triangulation for a particular value of the 58 parameters inside the chamber since the parametric polytope has the same combinato rial structure throughout the chamber The triangulation obtained for the instantiated vertices can then be carried over to the corresponding parametric vertices We only need to be careful to select a value for the parameters that does not lie on any facet of the chambers On these chambers some of the vertices may coincide For linearly parametrized polytopes it is easy to find a parameter point in the interior of a chamber as explained in Section Note that
66. 005 e Volume Calculation and Estimation of Parameterized Integer Polytopes Rabl 2006 e Improved Derivation of Process Networks Verdoolaege Nikolov and Stefanov 2006 e Computing the Ehrhart quasi polynomial of a rational simplex Barvinok 2006 Memory Optimization by Counting Points in Integer Transformations of Para metric Polytopes Seghir and Loechner 2006 GRAPHITE Polyhedral Analyses and Optimizations for GCC Pop Silber Co hen Bastoul Girbal and Vasilache 2006 e Volume Computation for Polytopes and Partition Functions for Classical Root Systems Baldoni Silva Beck Cochet and Vergne 2006 e A primal Barvinok algorithm based on irrational decompositions K ppe 2007 e pn A Tool for Improved Derivation of Process Networks Verdoolaege Nikolov and Stefanov 2007a Theoretical and Computational Methods for Lattice Point Enumeration in Insid e Out Polytopes van Herick 2007 e On Ehrhart Polynomials and Probability Calculations in Voting Theory Lepel ley Louichi and Smaoui 2008 e Local Euler Maclaurin formula for polytopes Berline and Vergne 2006 e Exact algorithms and software in optimization and polyhedral computation Fukuda 2008 e Enumeration of 4 x 4 magic squares Beck and van Herick 2010 118 e Symbolic and Analytic Techniques for Resource Analysis of Java Bytecode As pinall Atkey MacKenzie and Sannella 2010 e Parametric Integer Programming Algorithm for B
67. 1 with d ni The lattice associated to the linear subspace orthogonal to the facet is the projection of into this space Since n is primitive a basis for this lattice can be identified with n d The coordinate of the whole facet in this space is therefore da p n m v p while the transverse cone is daz p R Similarly a linear functional projects onto a linear functional n d in the linear subspace Applying 5 37 with y D 2 D and t nivi p nzv2 p we therefore find b n 1 mi 5 n Di 7 Dic d 0 mio y pj J i j 1 ger 177 i 0 j 0 After applying this differential operator to the polynomial p x the resulting poly nomial h p x h p x needs to be integrated over the facet The measure 71 to be used is such that the integral of a lattice tile in the linear space parallel to the facet is 1 i e f 1 i 14 1 0 0 with z the coordinate along f Referring to and 5 41 all points of the facet have the form x p zf a p while the z coordinate of the vertices v p and v2 p are n2v1 1 2 4 and n1v22 d respectively That is the contribution of the facet is equal to n2v2 1 1 V2 2 d f h p zf a p n dz navi a n1vi2 d where again we need to ensure that the lower limit is small
68. 5 Symbolic Polynomial Maximization over Convex Sets and its Application to Memory Requirement Estimation Clauss Fern ndez Garbervetsky and Ver doolaege 2006 e Counting with rational generating functions Verdoolaege and Woods 2008 Counting integer points in parametric polytopes using Barvinok s rational func tions Verdoolaege Seghir Beyls Loechner and Bruynooghe 2007b Polynomial Approximations in the Polytope Model Bringing the Power of Quasi Poly nomials to the Masses Meister and Verdoolaege 2008 Bounds on Quasi Polynomials for Static Program Analysis Devos Verdoolaege Van Campenhout and Stroobandt 2007 Computing parametric rational generating functions with a primal Barvinok al gorithm K ppe and Verdoolaege 2008 An Implementation of the Barvinok Woods Integer Projection Algorithm K ppe Verdoolaege and Woods 2008 117 e Algorithms for Weighted Counting over Parametric Polytopes A Survey and a Practical Comparison Verdoolaege and Bruynooghe 2008 6 2 Publications Refering to the Library This is a list of some reports and publications refering to the barvinok library e Formulas of Brion Lawrence and Varchenko on rational generating functions for cones Beck Haase and Sottile 2009 Generating Cache Hints for Improved Program Efficiency Beyls and D Hollander 2005 e Analternative algorithm for counting lattice points in a convex polytope Lasserre and Zeron 2
69. 99 Also expanding the numerator in the first factor we find C1 xS Xia di i AT gt 0 60 with constant term 1 4 p DET K tdy i c1 4 t 15 Cj Vizo i CD GO XL To compute the first d 1 terms in the Taylor series 5 27 we write down the truncated Taylor series tdg_i c1 ca 5 28 el d d 1 d 1 d 1 pu 1012 aED where we 1 xD 1 d 1 i 1 d 1 Zir Computing the reciprocal as explained in Section 5 10 we find t t 1 ay EtG Dg D 5 29 1 e 1 t i 0 Gat Note that the constant term of the denominator is 1 d 1 The denominators of the quotient are therefore d 1 d 1 Also note that the b are independent of the generating function and can be computed in advance An alternative way of computing the b is to note that cupi 1 zo with B i b the Bernoulli numbers which can be computed using the recurrence 5 34 see Section 5 12 Substituting by c t in 5 29 we have c t i 0 Multiplication of these truncated Taylor series for each c results in the first d 1 terms of 5 27 d d B td 7c 2 Bm p 0 0 d D from which it is easy to compute the constant term 5 28 Note that this convolution can also be computed without the u
70. ASIP Journal on Embedded Systems spe cial issue on Embedded Digital Signal Processing Systems 2007 118 Verdoolaege S R Seghir K Beyls V Loechner and M Bruynooghe 2007b June Counting integer points in parametric polytopes using Barvinok s rational functions Algorithmica 48 1 37 66 11411117 Verdoolaege S and M Bruynooghe 2008 July Algorithms for weighted counting over parametric polytopes A survey and a practical comparison In M Beck and T Stoll Eds The 2008 International Conference on Information Theory and Statistical Learning Verdoolaege S and K M Woods 2008 Counting with rational generating func tions J Symb Comput 43 2 75 91 105 Verdoolaege S G Janssens and M Bruynooghe 2009 June Equivalence check ing of static affine programs using widening to handle recurrences In Computer Aided Verification 21 pp 599 613 Springer 11 Wilde D K 1993 A library for doing polyhedral operations Technical Report 785 IRISA Rennes France http www irisa fr EXTERNE bibli pi pi785 html 22 Woods M 2004 Rational Generating Functions and Lattice Point Sets Ph D thesis University of Michigan 102 Woods K M 2005 Computing the period of an Ehrhart quasi polynomial The Electronic Journal of Combinatorics 12 R34 104 Woods K M 2006 June personal communication 54 List of Acronyms DNF Disjunctive Normal Form GOD cer greatest c
71. Figure 5 17 Examples of decompositions in primal space 51 Figure 5 18 Possible locations of w with respect to u and uj projected onto a plane orthogonal to the other rays when aja lt 0 Figure 5 19 Possible locations of w with respect to u and uj projected onto a plane orthogonal to the other rays when aja gt 0 shown in Figure 5 18 If both no and no are closed then y lies in region 1 and therefore n as well as is closed too Similarly if both no and no are open then so is nj If only one of the facets is closed then as explained above we choose nj to be open i e we take y to lie in region 3 or 5 Figure 5 19 shows the possible configurations for the case that gt 0 If exactly one of no and no is closed then y lies in region 3 or region 5 and therefore n is closed iff no is closed Otherwise as explained above we choose nj to be closed if i j The algorithm is summarized in Algorithm I where we use the convention that in cone u refers to uo w Note that we do not need any of the rays or normals in this code The only information we need is the closedness of the facets in the original cone and the signs of the o 5 4 Triangulation in primal space As in the case for Barvinok s decomposition Section 5 3 we can transform a trian gulation of a closed cone into closed simple cones into a triangulation of half open simple cones that fully partitions the orig
72. I i e I K 0 87 and use TOPCOM s points2triangs to compute the regular triangulations of the points specified by the rows of K Each of the resulting triangulations corresponds to a cham ber in the chamber complex of the above vector partition problem Each simplex in a triangulation corresponds to a parametric vertex active on the corresponding chamber and each point in the simplex i e a row of corresponds to a variable x or sj that is set to zero to obtain this parametric vertex In the original formulation of the prob lem each such variable set to zero reflects the saturation of the corresponding constraint x 0 for x and aj x bj p for s 0 A description of the chamber can then be obtained by plugging in the parametric vertices in the remaining constraints Example 5 61 Consider the parametric polytope P pq r t pl0xizpA0zjz2i q A0xkszi ptr ApzO qzO0Arz0 The constraints involving the variables are 1 gt 0 1 j 0 gt 0 1 0 Oj i s p 1 0 1 lt 4 2 1 Ojikis p r We have 0 0 100102077254 1 1 0 1 O 1 Of 0 9 9 2 1 0 0 0 1 0 0 0 1 Computing the regular triangulations of the rows of K using TOPCOM we obtain gt cat e2 topcom 1 0 8 2 9 1 1 1 8 1 0 8 9 1 8 9 9 1 gt points2triangs regular lt e2 topcom T 1 0 1 2 1 2 3 0 1 4 1 3 4 0 2 5 2 3 5 0 4 5 3 4 5 T 2 1 2 3 1 3 4 2 3 5 3 3 4 5 1
73. Journal of Combinatorics 14 3 251 258 44 Tawbi N 1994 Estimation of nested loops execution time by integer arithmetic in convex polyhedra In Proceedings of the 8th International Parallel Processing Symposium pp 217 221 IEEE Computer Society Press 63 125 Van Engelen R A K Gallivan and B Walsh 2006 September Parametric tim ing estimation with the Newton Gregory formulae Journal of Concurrency and Computation Practice and Experience 18 10 1434 1464 van Herick A W 2007 July Theoretical and computational methods for lattice point enumeration in inside out polytopes Master s thesis San Francisco State University 118 Verdoolaege S 2005 April Incremental Loop Transformations and Enumeration of Parametric Sets Phd Department of Computer Science K U Leuven Leu ven Belgium Verdoolaege S Beyls M Bruynooghe and Catthoor 2004a October Experiences with enumeration of integer projections of parametric polytopes Report CW 395 K U Leuven Department of Computer Science URL http www cs kuleuven ac be publicaties rapporten cw CW395 abs html Verdoolaege S K Beyls M Bruynooghe R Seghir and V Loechner 2004b March Analytical computation of Ehrhart polynomials and its applications for embedded systems In 2nd Workshop on Optimization for DSP and Embedded Systems ODES 2 117 Verdoolaege S K Beyls M Bruynooghe R Seghir and V Loech ner 2004c jan
74. PolyLib to unions of polyhedra It checks whether every polyhedron in the union D2 is included in some polyhedron of D1 Polyhedron DomainConstraintSimplify Polyhedron P unsigned MaxRays The value returned by DomainConstraintSimplify is a pointer to a newly allocated Polyhedron that contains the same integer points as its first argument but possibly has simpler constraints Each constraint g a X gt c is replaced by a x gt where g is the greatest common divisor gcd of the coefficients in the original constraint The Polyhedron pointed to by P is destroyed Polyhedron Polyhedron_Project Polyhedron P int dim The function Polyhedron_Project projects P onto its last dim dimensions Matrix left inverse Matrix M Matrix Eq The left inverse function computes the left inverse of M as explained in Section 5 6 Matrix Polyhedron Reduced Basis Polyhedron P struct barvinok options options Polyhedron Reduced Basis computes a generalized reduced basis of P which is as sumed to be a polytope using the algorithm of Cook et al 1993 See subsection 5 20 for more information The basis vectors are stored in the rows of the matrix returned Vector Polyhedron Sample Polyhedron P struct barvinok options options Polyhedron Sample returns an integer point of P or NULL if P contains no integer points The integer point is found using the algorithm of Cook et al 1993 and uses Polyhedron Reduced Basis to com
75. To appear 87 Feautrier P 1988 Parametric integer programming Operationnelle Operations Research 22 3 243 268 41 Feautrier 2006 Solving systems of affine in equalities PIP s user s guide 41 Fukuda K 1993 cdd c C implementation of the double description method for computing all vertices and extremal rays of a convex polyhedron given by a sys tem of linear inequalities Technical report Department of Mathematics Swiss Federal Institute of Technology Lausanne Switzerland program available from http www ifor math ethz ch fukuda fukuda html 92 Fukuda K 2008 Exact algorithms and software in optimization and polyhedral computation In ISSAC 08 Proceedings of the twenty first international sym posium on Symbolic and algebraic computation New York NY USA pp 333 334 ACM 118 Gawrilow E and M Joswig 2000 polymake a framework for analyzing convex polytopes In G Kalai and G M Ziegler Eds Polytopes Combinatorics and Computation pp 43 74 Birkhauser 44 Gelfand I M M Kapranov and A V Zelevinsky 1994 Discriminants Resul tants and Multidimensional Determinants Birkhauser Boston 87 Gomory R E 1963 An algorithm for integer solutions to linear programming In R L Graves and P Wolfe Eds Recent Advances in Mathematical Program ming New York pp 269 302 McGraw Hill 41 Hartmann M E 1989 Cutting planes and the complexity of the
76. ailable information about the parametric polytopes and will therefore compute many redundant triangula tions chambers In particular any chamber that does not intersect with the parameter domain of the parametric polytope or only intersects in a face of this parameter do main is completely redundant Furthermore if the parametric polytope is not simple then many different combinations of the constraints will lead to the same parametric vertex Many triangulations will therefore correspond to one and the same chamber in the chamber complex of the parametric polytope For example for a dilated octa hedron TOPCOM will compute 150 triangulations chambers 104 of which are empty while the remaining 46 refer to the same single chamber 5 19 Computing the Hilbert basis of a cone To compute the Hilbert basis of a cone we use the zsolve library from 4ti2 Hem OWT WE fir E move all equalities from the cone through unimodular transformations and then apply the technique of to put the cone in standard form Note that for a non parametric cone the constant term b in is 0 The constraints Dx gt c 0 of 5 58 are therefore equivalent to x gt 0 91 5 20 Integer Feasibility For testing whether a polytope P Qf contains any integer points we use the technique of Cook et al 1993 based on generalized basis reduction The technique basically looks for a short vector c in the latt
77. ain set with at most one value of x associated to each value of 1 In particular let 21 2 be the generating function of the integer points in P Then g x zi f x z1 1 with f x z1 z2 f x z1 z2 21 22 zo f x z1 22 is the generating function of T To check whether we need to subtract any shifted copies of g x z1 from itself we compute zi go zi 1 2 The second argument of this Hadamard product is depicted in Figure 5 92 by its coef ficients The exponents in h x zi that have a non zero coefficient are therefore those marked by both a dot e and a number The total sum can be computed as h 1 1 26 which is finite but non zero We therefore need to subtract at least one shifted copy of gG zi Let h x Z1 g x zi gG z B x Z1 B x zi 218 zi 115 Computing z go z 1 z we would find that h 1 1 0 and so we do not need to shift and subtract any further However since we are dealing with a two dimensional problem we already know from Theorem 5 87 that we can stop after one shift and subtract so we do not even need to compute h x 21 here The function g x zi now only has non zero coefficients for at most one exponent of z for each exponent of x We may therefore project down to obtain the function g x 1 which is the generating function of the set in the lower left offFigure 592 For the chamber of the parameter domain the compu
78. ating function and g t c s 2 ke a Q NZ 0 24 ke a Q NZ a s k with o x x Q converges in a neighborhood of el while S dw ke a Q NZ 104 Since c s gt 0 we have d k gt and the above sum diverges iff any of the coefficients of the negative powers of t in the Laurent expansion of g t is non zero If the sum converges then the sum is simply the coefficient of the constant term in this expansion It only remains to show now that we can compute a suitable in polynomial time i e an such that rz a gt 0 for any extremal ray of Q and bij o 0 for any b in 5 26 By adding the rx to the list of bj if needed we can relax the first set of constraints to rj gt 0 Let Q be described by the constraints Ax c gt 0 and let B be d X d non singular submatrix of A obtained by removing some of the rows of A Such a B exists since Q does not contain any straight line Clearly Br gt 0 for any ray r of Q Let bi then since b 0 and B is non singular we have b z 0 We may therefore find in polynomial time a point o gt 0 on the moment curve such that bj a 0 Algorithm 5 2 Let a B a Then bij bij Bbij a b a 0 and rz rj Brg a gt 0 as required Note that in practice we would as usual first try a fixed number of random vectors a gt 0 before resorting to
79. barvinok User Guide Contents Contents 12 Calculatoy 2 2 6 Counting Function 2 7 Auxiliary Function 3 4 barvinok union 3 5 barvinok ehrhar 37 1 scan NETTE 1 2 7 X Comparison to other Interactive Polyhedral Tool PolyLib interface of the barvinok librar 2 1 Existing Data Structure 2 3 Data Structures for Quasi polynomial 2 4 Operations on Quasi polynomial 2 5 Generating Function 3 6 polyhedron sample Version barvinok 0 37 Sven Verdoolaege April 15 2014 obsolescent 3 3 GG eh ae eee w 3 12 polyhedron integer hull 3 13 polytope lattice widt NECEM AMT 5 5 Multivariate quasi polynomials as lists of polynomial ete TERT 5 7 Integral basis of the orthogonal complement of a linear subspace 5 8 Ensuring a polyhedron has only revlex positive ra HT Tu 5 14 Summation using local Euler Maclaurin formula 5 14 1 Reduction to the summation of a parametric polynomial over parametric polytope with a fixed combinatorial structure 5 15 Summation through exponential substitution and Laurent expansions Old ol nos od dp ay de
80. being considered as a pair of bounding constraints with opposite outer normals However in our implemen tation we have chosen to first compute a maximal set of affinely independent points by first taking any point from S and then adding points from S not on one of the equal ities satisfied by all points found so far This allows us to not have to worry about equalities in the main algorithm In the case of the computation of the integer hull finding these affinely independent points can be accomplished using the technique of subsection 5 20 Example 5 65 Assume we want to compute the integer hull of the polytope in the left part of Figure 5 64 We first compute a set of three affinely independent points shown in the same part of the figure Of the three facets of the corresponding convex hull optimization along the outer normal depicted by an arrow in the figure of only one facet will yield any additional points The other two are therefore facets of the integer hull Optimization along the above outer normal may yield any of the points marked by a o Assuming it is the bottom one we end up with the updated convex hull in the middle of the figure This convex hull has only one new facet Adding the point found 93 Figure 5 64 The integer hull of a polytope by optimizing over this facet s outer normal we obtain the convex hull on the right of the figure There
81. brand and Shmonin 2007 which improves upon the technique of Kannan 1992 Given a parametric polytope P p x Ax b p gt 0 the width along a direction c is defined in the same way as for non parametric polytopes see subsection 5 20 width P p max c x x e P p min c x x e P p 5 72 The lattice width is the minimum width over all non zero integer directions width P p min width P p cez4N0 We assume that the parameter domain Q of P p i e the set of parameter values for which 0 is full dimensional and that for each from the interior of P p is also full dimensional Clearly for any given direction c the minimum and maximum in are attained at different vertices of P p The idea of the algorithm is then to consider all pairs of parametric vertices of P p to compute all candidate integer directions for a given pair of vertices and then to compute the minimum width over all candidate integer directions found For any given parametric vertex v p the rational directions for which this ver tex is minimal can be found as follows Let v p C be the vertex cone of v p If v p is minimal for c then all other points in the vertex cone must yield a bigger or equal value i e y c gt 0 for all y C The set of directions is therefore the po lar cone C of C Note that in principle we should only do this for pairs of vertices that have common activity domain where th
82. ction can then be performed The standard way of dealing with variables of unrestricted sign is to replace a vari able x of unknown sign by the difference x x x of two non negative variables x gt 0 However some algorithms are somewhat sensitive with respect to the number of variables and so we would prefer to introduce as few extra variables as pos sible We will therefore apply a unimodular transformation on the variables such that all transformed variables are known to be non negative The first step is to compute the HNF of A i e a matrix H AU with U unimod ular in column echelon form such that the first entry in each column is positive and the other entries on the corresponding row are non negative and strictly smaller than this first entry By reordering the rows we may assume that the top square part of H is lower triangular By a further unimodular transformation the entries below the diago nal can be made non positive and strictly smaller in absolute value than the diagonal entry of the same row For each of the new variables we can take a positive combination of the corre sponding row and the previous rows to obtain a positive multiple of the corresponding unit vector implying that the variable has a lower bound A slack variable can then be 86 introduced for each of the rows in the top square part of H that is not already a positive multiple of a unit vector and for each of the rows below the to
83. da 2003 November A user s guide for latte v1 1 software package LattE is available at http www math ucdavis edu latte De Loera J A R Hemmecke J Tauzer and R Yoshida 2004 Effective lattice point counting in rational convex polytopes The Journal of Symbolic Computa tion 38 4 1273 1302 47 De Loera J A and M K ppe 2006 Experiments with an algebraic scheme for estimating the number of lattice points in polyhedra Manuscript in preparation Devos H S Verdoolaege J Van Campenhout and D Stroobandt 2007 Bounds on quasi polynomials for static program analysis manuscript in preparation 117 121 Edmonds J L Lov sz and W Pulleyblank 1982 Brick decompositions and the matching rank of graphs Combinatorica 2 3 247 274 93 Ehrhart E 1977 Polyn mes arithm tiques et M thode des Poly dres en Combi natoire Volume 35 of International Series of Numerical Mathematics Basel S tuttgart Birkhauser Verlag Eisenbrand F 2000 July Gomory Chv tal cutting planes and the elementary clo sure of polyhedra Ph D thesis Universitat des Saarlandes Eisenbrand and Shmonin 2007 Parametric integer programming in fixed dimension Eisenschmidt E and M K ppe 2007 Integrally indecomposable polytopes and the survivable network design problem In Electronic proceedings of the 6th International Workshop on the Design of Reliable Communication Networks DRCN 2007
84. ding on whether the constraint is an equality 0 or an inequality 1 The number of equalities is stored in NbEq If the constraint is a X c gt 0 then the next columns contain the coefficients a and the final column contains the constant c The first column of Ray is either 0 or 1 depending on whether the generator is a line 0 or a vertex or ray 1 The number of lines is stored in NbBid Let d be the least common multiple lcm of the denominators of the coordinates of a vertex v then the next columns contain dv and the final column contains d For a ray the final column contains 0 The field next points to the next polyhedron in the union of polyhedra It is 0 if this is the last or only polyhedron in the union For more information on this structure we refer to Wilde 1993 Quasi polynomials are represented using the evalue and enode structures 22 typedef enum polynomial periodic evector enode_type typedef struct _evalue Value d denominator union Value n numerator if denominator 8 struct _enode p pointer if denominator 0 x evalue typedef struct _enode enode_type type polynomial or periodic or evector int size number of attached pointers int pos parameter position evalue arr 1 array of rational pointer enode If the field d of an evalue is zero then the evalue is a placeholder for a pointer to an enode stored in x p Otherwise the
85. dings of the 15th IMACS World Congress on Scientific Computation Modelling and Applied Mathematics Volume 2 of Wis senschaft amp Technik Verlag pp 685 690 163 Scarf H E 1981 March Production sets with indivisibilities part II The case of two activities Econometrica 49 2 395 423 34 Scarf H E and K M Woods 2006 Neighborhood complexes and generating functions for affine semigroups Discrete amp Computational Geometry 35 3 385 403 134 Seghir 2002 June D nombrement des point entiers de l union et de l image des poly dres param tr s Master s thesis ICPS Universit Louis Pasteur de Strasbourg France 129 Seghir 5 Verdoolaege Beyls and Loechner 2004 February Analytical computation of Ehrhart polynomials and its application in compile time gen erated cache hints Technical Report 118 ICPS Universit Louis Pasteur de Strasbourg France 117 Seghir and V Loechner 2006 October Memory optimization by counting points in integer transformations of parametric polytopes In Proceedings of the International Conference on Compilers Architectures and Synthesis for Em bedded Systems CASES 2006 Seoul Korea 105 Shoup V 2004 NTL Available from http www shoup net nt1 30 Stanley R P 1986 Enumerative Combinatorics Volume 1 Cambridge University Press 53 Stanley R P 1993 A monotonicity property of h vectors and h vectors Euro pean
86. e activity domains have been partially opened using the technique of Theorem 5 12 to avoid multiple vertices that coincide on a lower dimensional chamber to all be considered on this intersection However this optimization has currently not been implemented Given a pair of vertices v p and v2 p we may assume that v p attains the mini mum and v2 p attains the maximum If v p and vo p are the corresponding vertex cones then the set of rational directions for this pair of vertices is Ci Cj N C5 0 The set of candidate integer directions are therefore the vertices of the integer hull of C45 which can be computed as explained in subsection 5 22 To see this note that by construction c v p lt c vo p and so we p width P p e 2 vi p gt 0 98 Figure 5 73 A polytope and its candidate width directions Any integer direction in C will therefore yield a width that is at least as large as that of one of the vertices of the integer hull Note that when using generalized basis reduction to compute the integer hull of these cones as in subsubsection 5 22 2 it can be helpful to use as vertices for the initial approximation not only the extremal rays of the cone but also those vertices of previously computed integer hulls that are elements of the current cone After computing a list of all possible candidate width directions
87. e introduced four new types fractional relation partition and flooring Also known as laziness 26 enode enode type polynomial type fractional size 3 size 3 pos 1 pos 1 2 0 arr 0 5 9 a d 1 d 1 arr 1 A 3 arr 1 1 0 d 1 arr 2 a arr 2 zm 2 enode type polynomial size 2 pos 1 d 1 arr 0 IA 9 d 2 arr 1 x 1 Figure 2 4 The quasi polynomial 1 2 a 3 typedef enum polynomial periodic evector fractional relation partition flooring enode type The field pos is not used in most of these additional types and is therefore set to 1 The types fractional and flooring represent polynomial expressions in a frac tional part or a floor respectively The generator is stored in arr 9 while the coeff cients are stored in the remaining array elements That is an enode oftype fractional represents 1 arr 2 arr 0 arr 3 arr 8 arr 1 1 arr 0 ja An enode of type flooring represents arr 1 arr 2 larr 9 arr 3 arr 9 arr l 1 Larr 9 77 Example 2 3 The internal representation of the quasi polynomial Pao 5 1 2 3p 2 205 2 2 is shown in Figure 2 4 The relation type is used to represent strides In particular if the value of size is 2 then the value of a relation is in pseudo code 27
88. edings of the eighth annual symposium on Computational geometry pp 161 170 ACM Press 45 Barvinok A I 1994 Computing the Ehrhart polynomial of a convex lattice poly tope Dicrete Comput Geom 12 35 48 47 Barvinok A 2006 Computing the Ehrhart quasi polynomial of a rational sim plex Math Comp 75 1449 1466 118 Barvinok A I and J Pommersheim 1999 An algorithmic theory of lattice points in polyhedra New Perspectives in Algebraic Combinatorics 38 91 147 105 Barvinok A I and K M Woods 2003 April Short rational generating functions for lattice point problems J Amer Math Soc 16 957 979 1105 1108 Beck M Haase Sottile 2009 January Formulas of Brion Lawrence and Varchenko on rational generating functions for cones The Mathematical Intelligencer 31 1 118 Beck M and A W van Herick 2010 Enumeration of 4 x 4 magic squares Math Comp 118 Berline N 2007 August personal communication 70 Berline N and M Vergne 2006 July Local Euler Maclaurin formula for poly topes http arXiv org abs math 0507256 120 Beyls K and E D Hollander 2005 4 Generating cache hints for improved pro gram efficiency Journal of Systems Architecture 51 4 223 250 118 Bik A J C 1996 Compiler Support for Sparse Matrix Computations Ph D thesis University of Leiden The Netherlands 56 Boulet P and X Redon 2001 Sppoc manipulation
89. ege 2005 Example 36 on page 129 cat projected 56 39 k i j p cst 1 0 1 9 6 1 1 80 1 0 9 8 1 0 0 1 6 1 1 0 1 1 9 9 1 6 7 E 2 P 1 gt barvinok_enumerate_e lt projected POLYHEDRON Dimension 4 Constraints 5 Equations 1 Rays 4 Lines 0 Constraints 5 6 Equality 1 6 9 8 7 Inequality 0 1 0 1 Inequality 9 1 0 8 Inequality 0 0 1 Inequality 89 1 1 0 Rays 4 6 Vertex 50 8 1 1 1 1 0 0 0 1 Ray L 9 0 1 1 Vertex 8 1 1 1 1 1 exist 2 nparam 1 3 gt 0 1 gt 0 C3 P 10 P 1 gt 0 P 2 gt 0 8 P 0 Options floor f convert fractionals to floorings convert convert fractionals to periodics omega o Omega as a preprocessor pip p call barvinok enumerate pip instead of barvinok enumerate e isl i call barvinok enumerate isl instead of barvinok enumerate e 3 4 barvinok union The program barvinok union enumerates a union of parametric polytopes It takes as input the number of parametric polytopes in the union the polytopes in combined data 40 and parameter space in PolyLib notation the context in parameter space in PolyLib notation and optionally a list of parameter names Options series s compute rational generating function instead of piecewise step polynomial 3 5 barvinok ehrhart The program barvinok ehrhart computes the Ehrhart quasi polynomial of a poly tope i e a quasi polynomial in
90. enrici P 59 723 Hermite Normal Form HNF 55 56 86 Hickey T 58 27 Hilbert basis Hubler S L 119 723 Huggins P 93 27 723 Hung M S 87 23 identity incremental specialization index 127 inner point input format constraints vertices integer hull integer point 35 integer projection 105 intersection iscc 6 isl DAW SIDOLI23 4042 isl basic 5 isl basic set 5 isl 5 isl_polyhedron_minimize isl pw qpolynomial 5 isl pw qpolynomial from evalue isl set 5 isl union 5 isl union pw qpolynomial 5 151 union set 5 iteration domain 7 Janssens G Jiang A X Joswig M 44 M 45 47 49 59 60 63 69 87 Kannan R 98 105 727 723 Kapranov M 22 Kelly W 13 25 Koeppe M 119 723 130 Lagarias J C 92 Lasserre J B last 7 LattE 36 60 63 LattE macchiato 601163 latte2polylib p1l 36 lattice width 43 92 98 98 LATTICE POINTS 44 lattice points 44 lattice width 16 laurent 42 77 Lebesgue measure 69 Lee C W 58 24 left inverse left inverse 35 55 Lenstra Jr W 24 Lenstra Lenstra and Lovasz basis reduc tion algorithm LLL 26 Lepelley D 118 724 lexmax 12 16 lexmin I 12 16 41 Leyton Brown K 119 25 lift 16 Lison k P Litvak A E 20 LLL a LLL b 26 local E
91. er than the upper limit using the usual method of plugging in a particular value of the parameters Finally we consider the contributions of the vertices The transverse cones are in this case simply the supporting cones Since p is a valuation we may apply Barvinok s decomposition and assume that the cone is unimodular For an affine cone K v p Rri Ry 1 2 R r2 with 1 r vp we have Berline Vergne 2006 Proposition 31 efiyithy2 1 u K amp Biy2 Ciy1 t2 BO 5 42 eny 1 e yiy2 with v b n 1 7 1 Dy BY 1 y Hei yi ri Ci V1 v2 vj Vi and t es Expanding 5 42 we find 60 0 SS bn Ln 60 22 SS bn Lo mel 5 2 Went Jl 2 CE s 1 5 75 x b n 1 t2 2 1 1 n 1 44 xD n b n Lt Yt p b n 1 5 Cay yq n 1 44 D y2 hi E 2 2 2 c C1 Co ti t2 72 with b n 1 1 b m 1 tr nj 1 0 1 b n 1 5 no 1 m 10 b n m ul 1 c Ci C2 tj n1 n2 jec 1 1 1 c ye For D we have yp GiiDi ri3D2 r23Di r22D2 ni n k k 1 pk pem 1 72 7 1 gt ratio Jo 195 Joi 1 0 and so Dpwp vp H K D 3 n 1 M72 ni nj Arnim rts ouo
92. ere exists a list of det L polynomials gy Q Ti Tn fori in the fundamental parallelepiped of L such that f s gils ifs i mod L To compute the period lattice from a fractional representation we compute the appropriate lattice for each fractional part and then take their intersection Recall that the argument of each fractional part is an affine expression in the parameters a p c m with a Z and c m Z Such a fractional part is translation invariant over any integer value of p such that a p mt 0 for some t Z Solving this homogeneous equation over the integers in our implementation we use PolyLib s SolveDiophantine gives the general solution PI U n for x 2 The matrix U Z then has the generators of the required lattice as columns The constituents are computed by plugging in each integer point in the fundamental par allelepiped of the lattice These points themselves are computed as explained in Sec tion 5 2 Note that for computing the constituents it is sufficient to take any represen tative of the residue class For example we could take w k W in the notations of Lemma 5 2 Example 5 24 Woods 2006 Consider the parametric polytope P i x 0xxz s c 0 2 The enumerator of Ps is 5 2 0 Sp i Z EUM im B 0 0 Js 21 2 1 titi 22 ol The corresponding output of barvinok enumerate is 54 s t gt 0 1 g
93. eriodic numbers for coefficients The period of a quasi polynomial is the lcm of the periods of its coefficients Other authors e g Stanley 1986 use the following definition of a quasi polynomial Definition 5 22 A function f Z Q is a univariate quasi polynomial of period q if there exists a list of q polynomials g Q T for 0 lt i lt q such that f s gils ifssi The functions g are called the constituents 53 In our implementation we use Definition 5 21 but whereas Ehrhart 1977 uses a list of q rational numbers enclosed in square brackets to represent periodic numbers our periodic numbers are polynomial expressions in fractional parts Section 2 3 These fractional parts naturally extend to multivariate quasi polynomials The bracketed plicit periodic numbers can be extended to multiple variables by nesting them e g Loechner 1999 Definition 5 22 could be extended in a similar way by having a constituent for each residue modulo a vector period q However as pointed out by Woods 2006 this may not result in the minimum number of constituents A vector period can be considered as a lattice with orthogonal generators and the number of constituents is equal to the index or determinant of that lattice By considering more general lattices we can potentially reduce the number of constituents Definition 5 23 A function f Z Q is a multivariate quasi polynomial of period L if th
94. espect since the polynomial substitution unlike the exponential substitution had not been optimized to take full advantage of the stopped Barvinok decomposition For comparison Table 2 shows running times for the same experiments of that paper but using barvinok version barvinok 9 23 47 gaa9024e on an Athlon MP 1500 with 512MiB internal memory This machine appears to be slightly slower than the machine used in the experiments of K ppe 2007 as computing hi ckerson 14 using the dual decomposition with polynomial substitution and maximal index 1 took 2768 seconds on this machine using LattE macchiato At this stage it is not clear yet why the number of cones in the dual decomposition of hickerson 13 differs from that of LattE and LattE macchiato K ppe 2006 We conclude from Table 2 that our implementation of the exponential substitution is always slightly faster than our implementation of the polynomial substitution The optimal maximal index for these examples is about 500 which agrees with the experiments of K ppe 2007 5 12 Approximate Enumeration using Nested Sums If P Q is a polyhedron and p x Q x is a polynomial and we want to sum p x over all integer values of a subset of the variables x then we can do this incre mentally by taking a variable x with lower bound L X and upper bound U X with amp Xq and computing U amp 5 31 x1 L amp Since P is a polytope the
95. f m to any element y in the domain of mz such that their images and overlap and such that lt may compute map that contains for any element in the domain of mz a mapping from the lexicographically last element x in the domain of m according to the schedule to y such that m x and m2 y overlap and such that m3 x lt m3 y Return a list contain ing this map and the set of elements in the domain of m for which there is no corresponding element in the domain of compute a map that contains for any element y in the domain of m3 a mapping from the lexicographically last element in the domain of m according to the schedule m4 to y such that and overlap and such that m4 x lt ma y as well as from any element z in the domain of m such that m z and ma y overlap m4 z lt m4 y and such that there is no x in the domain of m with 0 and ma z lt lt continued on next page 17 Syntax Meaning schedule respecting minimizing m 1 schedule forest respecting minimizing m i3 i th 3 1 S1 52 93 4 q2 h fi q h 4 fi S3 S1 5 5 204 8 5 S10 i5 ii i 3 S1 So m3 M 2 93 q1 42 13 ii in 3 S1 S2 M m q2 i q qQ q i q3 41 4 fi mM 5
96. formulated for a non parametric polytope and a non parametric polynomial However as we will see in each of the steps in the computation the parameters can be treated as symbolic constants without affecting the validity of the formula see also Berline and Vergne 2006 Section 6 The differential operator Dp y is obtained by plugging in the vector D D of first order differential operators i e D is the first order differential operator in the kth variable in the function up p rp This function is determined by the transverse cone of the polyhedron P p along its face F p which is the supporting cone of P p along F p projected into the linear subspace orthogonal to F p The lattice associated to this space is the projection of A into this space In particular for a zero dimensional affine cone in the zero dimensional space we have 1 Berline and Vergne 2006 Proposition 12 while for a one dimensional affine cone t R r in the one dimensional space where r is a primitive integer vector and f 0 1 we have Berline and Vergne 2006 13 b n 1 1 Ds 1 5 5 37 MK E 69 Dg with y amp and b n t the Bernoulli polynomials defined by the generating series o y _ v bnt Todd t y lt gt 2 mE 5 38 The constant terms of these Bernoulli polynomials are the Bernoulli numbers Applying 5 36 to a one dimensional parametric polytope P p vi
97. from Example 5 69 The initial approximation is C conv Q 3 3 4 pos 2 3 3 4 which is shown on the left of Figure 5 70 The only bounding constraint that does not correspond to a bounding constraint of C is 7x y 17 In the first step we will therefore look for a point minimizing 7x y with values in the interval 1 16 All values of this objective function in the given interval attained by points in C are shown in Figure 5 66 From Example 5 67 we know that the optimal value is 6 and this cor responds to the point 1 1 Adding this point to our hull we obtain the approximation in the middle of Figure 5 70 This approximation has two new facets The bounding constraint 2y gt 1 will not produce any new points since we would be looking for in the interval 1 0 The other new bounding constraint is 4x y gt 5 Minimiz ing 4x y with values in the interval 1 4 we find the minimal value 3 corresponding to the point 1 1 Adding this point we obtain the complete integer hull shown on 97 the right of Figure 5 70 Note that if in the first step we would have added not only the point corresponding to the optimal value but instead all points found in Example 5 67 then we would have obtained the complete integer hull directly 5 23 Computing the lattice width of a parametric polytope To compute the lattice width of a parametric polytope we essentially use the technique of Eisen
98. he same values for the iterators isl allows spaces to be named The name is placed in front of the enclosing the iterators Consider for example the program in Figure 1 5 The program contains three statements which have been labeled for convenience The iteration domain of the first statement T can be represented as n gt Ti 0 lt i lt n The union of all iteration domains can be represented as n gt T i 0 lt i lt F i j 0 lt i lt n and 0 lt lt n i B i 9 i n The semicolon is used to express a disjunction between spaces This should be con trasted with the or keyword which expresses a disjunction between conjunctions of constraints For example the result of the following iscc statement is True i i 3 ori 5 3 5 X 1 2 3 Maps and Access Relations A second important concept in the polyhedral model is that of an access relation An access relation connects iterations of a statement to the array elements accessed by those iterations Such a binary relation can be represented by a map in isl Maps are defined in similar way to sets except that the single space is replaced by a pair of spaces separated by gt As an example consider once more the program in Figure 1 5 In particular consider the access t j 1 in statement The index expression maps the pair of iterations i and j to 3 1 i e element j 1 of a space with name t Ignoring the loop bound con
99. he sum of x x over the integer points of P n 78 with n n n5 nj We observe that the sum of the monomial t over the integer points in P is equal to n times the coefficient of the y term in the Taylor expansion of fle As in the case of unweighted counting see subsection 5 11 we have to add the coefficients of these monomials in the Laurent expansions of the terms in 5 26 How ever unlike the case of unweighted counting we cannot transform this problem to a univariate problem and computing the coefficient of a monomial in the Laurent expan sions does not reduce to computing the coefficient of a single higher degree monomial in a Taylor expansion Consider now one of the terms g x f x in 5 26 si DXb y O amp T i 5 46 with w p ic s p b written in terms of the b which are assumed to form a basis and where we have made explicit the only place where the parameters p appear We rewrite this equation as d 1 5 ys auc cS Mes bj 5 leo 5 47 j l The second factor is analytic and is a product of generating functions Todd s p bj y of Bernoulli polynomials 5 38 Plugging in these expressions we find bj yyesiPXbj 5 2 y 2 0 P E k gt 0 ki b ki SP _ b y 2 ki d Ek Ek Qt iX Gul 2 i k p i 1 Todd s
100. heconeofdirections Co 100 5 76 Theconeofdirections C41 101 Ern 101 BEA 107 5 90 The parametric polytope from Example 5 91 for different values of the parametef OR e RR OR BR o de RD 113 5 92 The transformed parametric polytope from Example 1 151 interface 11 Library The barvinok library currently supports only a few functions that interface with the isl library In time this interface will grow and is set to replace the PolyLib interface For more information on the 151 data structures see the isl user manual __isl_give isl pw gpolynomial isl basic set card C __151 take isl basic set bset __isl_give isl pw qpolynomial isl set card isl take isl set set __isl_give isl union pw qpolynomial isl union set card __151 take isl union set uset Compute the number of elements in an isl basic set isl set or isl union set The resulting isl pw qpolynomialorisl union pw qpolynomial has purely para metric cells __isl_give isl pw qpolynomial isl basic map cardC __isl_take isl basic map bmap __isl_give isl pw qpolynomial isl map card isl take isl map map isl give isl union pw qpolynomial isl union map card isl take isl union map umap Compute a closed form expression for the number of image elements associated to any element in the domain of the given isl basic map isl or
101. ib PipLib and omega to perform its operations Although the first two of these libraries may have been state of the art at the time SPPoC was created they are no longer actively maintained and have been largely superseded by more recent libraries In particular PipLib effectively only supports a single oper ation which is now also available in both isl and PPL The operations on rational polyhedra in PolyLib are also available in PPL usually through a cleaner interface and with a more efficient implementation As to counting the elements in a paramet ric polytope Barvinok s algorithm implemented in the barvinok library is usually much faster than the algorithm implemented in PolyLib Verdoolaege et al 2007b Furthermore the ability to work with named and nested spaces and the ability of sets and maps to contain pairs of elements from different spaces are not available in the Omega calculator and SPPoC 14 Table 1 iscc operations The variables have the following types s set map q piecewise quasipolynomial f piecewise quasipolynomial fold list i integer boolean S string o object of any type Syntax Meaning S2 aff s affine hull of s m aff m affine hull of m q card s number of elements in the set s q card m number of elements in the image of a domain ele 59 coalesce s m coalesce m 2 coalesce q coalesce fi codegen s codegen m codegen s using m
102. ice 24 where short ness is measured in terms of the width of the polytope P along that direction width P max x x Pj min x x y E P The lattice width is the minimum width over all non zero integer directions width P min width P ceZ4 0 If the dimension d is fixed then the lattice width of any polytope P Q containing no integer points is bounded by a constant Lagarias et al 1990 Ba naszczyk et al 1999 If we slice the polytope using hyperplanes orthogonal to a short direction i e a direction where the width is small we will therefore only need to in spect few of them before either finding one with an integer point or running out of hyperplanes meaning that the polytope did not contain any integer points Each slice is checked for integer points by applying the above method recursively A nice feature of this technique is that it will not only tell you if there is any integer point in the given polytope but it will actually compute one if there is any The short vector is obtained as the first vector of a reduced basis of the lattice Z with respect to the polytope In particular the first vector b of this reduced basis will satisfy width P widthy P lt z7 with 0 lt lt 1 2 a fixed constant That is the width in direction b is no more than a constant factor bigger than the lattice width See for details In our
103. ilevel Mixed Integer Programs Koeppe Queyranne and Ryan 2010 e Symmetric games with piecewise linear utilities Ryan Jiang and Leyton Brown 2010 e Counting chemical compositions using Ehrhart quasi polynomials Hubler and Craciun 2012 119 References Aspinall D R Atkey K MacKenzie and D Sannella 2010 Symbolic and ana lytic techniques for resource analysis of java bytecode In Trustworthly Global Computing 5th International Symposium pp 1 22 Springer 119 Baldoni V N Berline and M Vergne 2008 March Sum over lattice points of a polygon with iterated Laurent series user s guide 77 Baldoni V N Berline and M Vergne 2009 May Summing a polynomial func tion over integral points of a polygon user s guide 82 Baldoni Silva M W M Beck C Cochet and M Vergne 2006 Volume compu tation for polytopes and partition functions for classical root systems Discrete amp Computational Geometry 35 4 551 595 118 Banaszczyk W A E Litvak A Pajor and S J Szarek 1999 August The flatness theorem for nonsymmetric convex bodies via the local theory of banach spaces Mathematics of Operations Research 24 3 728 750 19211108 Barvinok A 2002 A Course in Convexity Volume 54 of Graduate Studies in Mathematics Providence RI American Mathematical Society 9211108 Barvinok A I 1992 Computing the volume counting integral points and expo nential sums In Proce
104. implementation we use 1 4 When used in the above integer feasibility testing algorithm we will also terminate the reduced basis computation as soon as the width along the first basis vector is smaller than 2 This means that there will be at most 2 slices orthogonal to the chosen direction The computation of the above reduced basis requires the solution of many linear programs for which we use any of the following external solvers e GLPK Makhorin 2006 This solver is based on double precision floating point arithmetic and may there fore not be suitable if the coefficients of the constraints describing the polytope are large cdd Fukuda 1993 This solver is based on exact integer arithmetic Note that you need version cddlib 6 94e or newer Earlier versions 0 93 0 94d have a bug that may sometimes result in a polytope being reported as rationally empty even though it is not 92 The LP solver to use can be selected with the gbr option 5 21 Computing the integer hull of a polyhedron For computing the integer hull of a polyhedron we first describe how to compute the convex hull of a set given as an oracle for optimizing a linear objective function over the set and then we explain how to optimize a linear objective function over the integer points of a polyhedron Applying the first with the second as optimization oracle yields a method for computing the requested integer hull 5 211 Computing the convex hull based o
105. implementations of the algorithm of Loechner and Wilde 1997 for computing parametric vertices and the algorithm of Clauss and Loechner 1998 for computing chamber decompositions Initially our library was meant to be a replacement for the algorithm of Clauss and Loechner 1998 also implemented in PolyLib for computing quasi polynomials To ease the transition of application programs we tried to reuse the existing data structures as much as possible 21 Existing Data Structures Inside PolyLib integer values are represented by the Value data type Depending on a configure option the data type may either by a 32 bit integer a 64 bit integer or an arbitrary precision integer using GMP The barvinok library requires that PolyLib is compiled with support for arbitrary precision integers The basic structure for representing unions of polyhedra is a Polyhedron typedef struct polyhedron unsigned Dimension NbConstraints NbRays NbEq NbBid Value Constraint Value Ray Value p_Init int p_Init_size struct polyhedron next Polyhedron The attribute Dimension is the dimension of the ambient space i e the number of variables The attributes Constraint and Ray point to two dimensional arrays of con straints and generators respectively The number of rows is stored in NbConstraints and NbRays respectively The number of columns in both arrays is equal to 1 Dimension 1 The first column of Constraint is either 0 or 1 depen
106. inal cone i e such that the half open sim 52 Algorithm 1 Determine whether the facet opposite uj is closed in Kj if Qj 0 closed K u closed u else ifi j if Qj gt 0 closed K u closed K u else closed K u closed K u else if aja gt 0 if closed K uj closed K u closed K u i lt j else closed K u closed K u else closed K u closed K u and closed K u ple cones do not intersect at their facets Again we apply Proposition 5 12 with y an interior point of the cone Section 5 1 Note that the interior point y may still inter sect some of the internal facets so we may need to perturb it slightly In practice we apply a lexicographical rule for such internal facets which always appear in pairs we close the one with a lexico positive normal and open the one with a lexico negative normal 5 5 Multivariate quasi polynomials as lists of polynomials There are many definitions for a univariate quasi polynomial Ehrhart 1977 uses a definition based on periodic numbers Definition 5 20 A rational periodic number U p is a function Z Q such that there exists a period q such that U p U p whenever p mod q Definition 5 21 A univariate quasi polynomial f of degree d is a function f n ci n n where are rational periodic numbers I e it is a polynomial expression of degree d with rational p
107. ion to the stream pointed to by DST The argument pname points to an array of character strings repre senting the parameter names The array is assumed to be long enough int eequal const evalue el const evalue e2 29 The function eequal return true 1 if its two arguments are structurally identical Le it does not check whether the two piecewise quasi polynomial represent the same function void reduce evalue evalue e The function reduce evalue performs some simplifications on evalues Here we only describe the simplifications that are directly related to the internal representation Some other simplifications are explained in Verdoolaege 2005 Section 4 7 2 If the highest order coefficients of a polynomial fractional or flooring are zero pos sibly after some other simplifications then the size of the array is reduced If only the constant term remains 1 the size is reduced to 1 for polynomial or to 2 for the other types then the whole node is replaced by the constant term Additionally if the argument of a fractional has been reduced to a constant then the whole node is re placed by its partial evaluation A relation is similarly reduced if its second branch or both its branches are zero Chambers with zero associated quasi polynomials are discarded from a partition 2 5 Generating Functions The representation of rational generating functions uses some basic types from the NTL library Shoup 2004 for
108. j 0 k l i nj20 Oxkxn n2 gt 0 O lt I lt nz 8 T The contribution of this vertex is then h p v p with h p x Example 5 43 As simple example consider the non parametric triangle in Fig ure 5 44 and assume we want to compute xeTNZ Since T N Z 2 4 3 4 2 5 the result should be 2 44 3 44 2 5 30 Let us first consider the integral xix 5 7 du Integration along each of the edges of the triangle yields the following For the edge in the margin we have f 1 0 i e f 0 The contribution of this edge to the integral is therefore zero For this edge we have f 1 1 The contribution of this edge to the integral is EU therefore 3 AY 4 2937 T 24 73 2 5 2 4 3 4 Figure 5 44 Sum of polynomial over the integer points in a triangle For this edge we have f 0 1 The contribution of this edge to the integral is 4 B N oS therefore 1 42 2 5 A TOT A 0 2 The total integral is therefore mn 121 dx 9 or 2 24 24 Now let us consider the contributions of the edges We will need the following Bernoulli numbers in our computations 1 b 1 0 1 0 5 b 2 0 63 0 0 1 b 4 0 4 0 30 The normal to the facet F in the margin is n 0 1 The vector f 1 0 is parallel t
109. ler Maclaurin formula to compute the sum of a piecewise polynomial evaluated in all integer points of a two dimensional parametric polytope For the theory behind these formula and a discussion of the original implementation for non parametric simplices we refer to Berline and Vergne 2006 In particular consider a parametric piecewise polynomial in n parameters and m variables c 2 gt Z gt Q po c p with c p Z gt Q x gt c p x and ci p x ifxe Di p Cp x if x D p with the c polynomials c Q p x and the D disjoint linearly parametric poly topes We want to compute g p gt co 2 5 14 1 Reduction to the summation of a parametric polynomial over paramet ric polytope with a fixed combinatorial structure Since the D are disjoint we can consider each c D pair individually and compute r g p o i 1 i 1 xeD p nZ The second step is to compute the chamber decomposition Verdoolaege 2005 Section 4 2 3 of each parametric polytope D The result is a subdivision of the parameter space into chambers C such that D has a fixed combinatorial structure in particular a fixed set of parametric vertices on the interior of each C Applying Theorem 5 12 68 this subdivision can be transformed into a partition C j by making some of the facets of the chambers open Section 3 2 Since we are only interested in integer parameter
110. ll the points in the set The syntax for this weighted counting is to compute the sum of a piecewise quasipolynomial over its domain As in the case of the ub operator if the domain contains a pair of nested space the sum is computed over the range As an example the result of sum i gt j gt i j 0 lt lt 12 ist i gt 1 2 1 2 1 2 i73 i gt 9 After the computation of some sum or bound the result may have to be reformu lated in terms of other variables For example during inter procedural analysis a result computed in terms of the formal parameters may have to be reformulated in terms of the actual parameters iscc therefore allows maps and piecewise quasipolynomials or folds to be composed If the map is multi valued then the composition maps each domain element to the sum or upper bound of the values at its image elements Finally because of its high level representation iscc can provide a dependence analysis operation taking only three maps as input the sink accesses the potential source accesses and a schedule The result is a single dependence map 1 2 6 More Examples n m gt i j 0 lt i lt andi lt j lt m card P f n m gt i j gt i j n i i j i j gt 0 and 51 27j lt numm sum f s sum f s n m gt 0 lt n m lt 20 f n gt i gt 2 n i 3 n 1 2 i i 3 2 i 1 exists j 0 i
111. lower bound is a maximum of affine expressions in the remaining variables while the upper bound is a minimum of such expressions If the coefficients in these expressions are all integer then we can compute Q X exactly as a piecewise polynomial using formulas for sums of powers as proposed by e g Tawbi 1994 Sakellarion 1997 Nan Engelen e al 2006 IF some of the coeficients are not integer we can apply the same formulas to obtain an approximation which can is some cases be shown to be an overapproximation Van Engelen et al 2006 Note that if we take the initial polynomial to be the constant 1 then this gives us a method for computing an approximation of the number of integer points in a parametric polytope The first step is to compute the chamber decomposition of P when viewed as a 1 dimensional parametric polytope That is we need to partition the projection of P onto the remaining variables into polyhedral cells such that in each cell both the upper and the lower bound are described by a single affine expression Basically for each pair of lower and upper bound we compute the cell where the chosen lower bound is strictly smaller than all other lower bounds and similarly for the upper bound For any given pair of lower and upper bound u X the formula is computed for each monomial of p x separately For the constant term we have 63 Dual decomposition Primal decomposition Time s Time s Max index
112. ly defined by the other constraints Example 5 63 Consider the parametric polytope i p llzi 2is3j jzn The constraints are 1 O0 1 2 3 gt 0 o JU n The top 2 x 2 submatrix is already in HNF We have 3j 2 2i 2 2 so we can add a constraint of the form j 2 c n and obtain EE while K with A I K 0 is given by 1 2 3 90 The second row of K corresponds to the second variable which in turn corresponds to the newly added constraint Passing all rows of K to TOPCOM we would get gt points2triangs regular lt lt EOF gt 1 0 0 1 0 1 2 3 gt EOF TL1 0 1 0 2 1 3 2 3 T 2 0 2 2 3 0 33 TL3 The first vertex in the first chamber saturates the second row row 1 and therefore saturates both the first 0 and fourth 3 and it appears a second time as 1 3 Combining these two vertices into one as 8 33 results in the second identical chamber Removing the row corresponding to the new constraint from K we remove the duplicates points2triangs regular EOF gt 1 0 9 1 2 3 gt EOF T 1 0 1 1 2 0 23 T 2 Note that in this example we also could have interchanged the second and the third constraint and then have replaced j by j In practice this method of computing a chamber decomposition does not seem to perform very well mostly because TOPCOM can not exploit all av
113. ment in recent ver sions of PolyLib as shown below Enumeration barvinok enumerate Polyhedron P Polyhedron C unsigned MaxRays extern Enumeration Polyhedron Enumerate Polyhedron P Polyhedron C unsigned MAXRAYS char pname The argument MaxRays has the same meaning as the argument NbMaxCons above The argument P refers to the d n dimensional polyhedron defining the parametric polytope The argument C is an n dimensional polyhedron containing extra constraints on the parameter space Its primary use is to indicate how many of the dimensions in P refer to parameters as any constraint in C could equally well have been added to P itself Note that the dimensions referring to the parameters should appear last If either P or C is a union then only the first set in the union will be taken into account The result is a newly allocated Enumeration As an alternative we also provide a function barvinok enumerate ev or barvinok enumerate with options that returns an evalue evalue barvinok enumerate ev Polyhedron P Polyhedron C unsigned MaxRays evalue barvinok enumerate with options Polyhedron P Polyhedron C struct barvinok options options For enumerating parametric sets with existentially quantified variables we provide two functions barvinok enumerate e and barvinok enumerate isl evalue barvinok enumerate e Polyhedron P unsigned exist unsigned nparam unsigned MaxRays 33 evalue barvinok enumerate e
114. meter in the call to NTL s implemen tation of Lenstra Lenstra and Lovasz basis reduction algorithm LLL e barvinok specific options incremental specialization Selects the specialization algorithm to be used If set to O then a direct spe cialization is performed using a random vector Value 1 selects a depth first incremental specialization while value 2 selects a breadth first incremen tal specialization The default is selected by the enable incremental configure option For more information we refer to Verdoolaege 2005 Section 4 4 3 2 3 Data Structures for Quasi polynomials Internally we do not represent our quasi polynomials as step polynomials but instead as polynomials of fractional parts of degree 1 polynomials However we also allow our quasi polynomials to be represented as polynomials with periodic numbers for coeffi cients similarly to Loechner 1999 By default the current version of barvinok uses fractionals but this can be changed through the disable fractional config ure option When this option is specified the periodic numbers are represented as an explicit enumeration using the periodic type A quasi polynomial based on fractional parts can also be converted to an actual step polynomial using evalue frac2floor but this is not fully supported yet For reasons of compatibility we shoehorned our representations for piecewise quasi polynomials into the existing data structures To this effect w
115. n a generating function f x y for the set s t Ju Z s t u P 5 85 and then compute C x f x 1 There are however some complications Most no tably after applying the technique in one direction and projecting out the corresponding variable the resulting set i e S S t ui Um 1 dus Z stu e in general does not correspond to the integer points in some polytope For example as sume that the polytope in Figure 5 84 contains the values of u associated to a particular value of s t Since there are integer points in this polytope we should count this value of t but only once If we apply the above technique in the vertical direction u2 then we can compute a generating function for the set 5 shown on the bottom of the figure Note however that there are gaps in this set i e if we compute S V e 4 S then we will not end up with a single point for this value of s t Since the biggest gap is three wide we need to compute S eai S NQeiaqa S VOeiqa 5 to obtain a single point If we do the subtraction in the horizontal direction first then we end up with a set shown on the left with gaps at most two wide so afterwards we only need to subtract twice in the vertical direction In general there is no bound on the widths of the gaps we may encounter in any given direction However there are directions in which the gaps are known to be 106
116. n an optimization oracle The algorithm described below is presented by Cook et al 1992 Remark 2 5 as an extension of the algorithm by Edmonds et al 1982 Section 3 for computing the di mension of a polytope for which only an optimization oracle is available The algorithm is described in a bit more detail by Eisenbrand 2000 and reportedly stems from Hart mann 1989 Essentially the same algorithm has also been implemented by Huggins 2006 citing beneath beyond Preparata and Shamos 1985 as his inspiration The algorithm start out from an initial set of points from the set S After computing the convex hull of this set of points we take one of its bounding constraints and use the optimization oracle to compute an optimal point in 5 but on the other side of the bounding hyperplane along the outer normal of this bounding constraint If a new point is found it is added to the set of points and a new convex hull is computed or the old one is adapted in a beneath beyond fashion Otherwise the chosen bounding constraint is also a bounding constraint of S and need not be considered anymore The process continues until all bounding constraints in the description of the current convex hull have been considered In principle the initial set of points in the above algorithm may be empty with a convex hull described by a set of conflicting constraints and each equality in the description of any intermediate lower dimensional convex hull
117. n may access any element of the row The access relation corresponding to A 1 is therefore n gt Sli gt A i j 0 lt i j lt n This map associates n elements of A to each iteration of S 1 2 3 Nested Spaces Each space may contain a nested pair of spaces Such nested spaces are extremely useful in more advanced applications As an example suppose that during equivalence checking Verdoolaege et al 2009 of two programs the iterations of S1 in one program are supposed to produce the same results as the same iterations of SA in the other program which may be described using the following map n gt S1 i gt SA i lt i lt n 11 If the iterations of S1 depend on the same iterations of S2 i e S1 i gt S2 i while those of SA depend on the next iteration of B i e SA i gt SB i 1 then we can apply the cross product of these two dependence maps i e i S1 i gt SA i gt S2 i gt SB 1 i to the original map to find out which iterations of S2 should correspond to which iter ations of SB 1 2 4 Basic Operations Basic operations on sets and maps include intersection union difference cross product cross sampling sample affine hull aff lexicographic optimiza tion lexmin or lexmax subset lt and equality tests code generation codegen and the cardinality card Additional operations on maps include computing domain dom and range ran
118. nd bij Z4 0 After computing the rational generating function 5 26 of a polytope with k d for all i the number of lattice points in the polytope can be obtained by evaluating f 1 Since 1 is a pole of each term we need to compute the constant term in the Laurent expansions of each term in 5 26 about 1 Since it is easier to work with univariate series a substitution is usually applied either a polynomial substitution x42 as implemented in LattE De Loera et al 2003 or an exponential substitution see e g Barvinok and Pommersheim 1999 x e as implemented in LattE macchiato K ppe 2006 In each case A Zf is a vector that is not orthogonal to any of the Both substitutions also transform the problem of computing the constant term in the Laurent expansions about x 1 to that of com puting the constant term in the Laurent expansions about t 0 Here we discuss the exponential substitution Consider now one of the terms in 5 26 Q r S Ila 1 m with wi A and c bij We rewrite this equation as ic ean d a E il 1 ect The second factor is analytic in a neighborhood of the origin t c 0 and therefore has a Taylor series expansion H die y tds c1 5 27 I ice L AUS where td is a homogeneous polynomial of degree m called the m th Todd polyno mial Barvinok and Pommersheim 19
119. nstraints 2 Equations 9 Rays 2 Lines 0 Constraints 2 3 Inequality 1 8 Inequality 2 13 Rays 2 3 Vertex 0 1 1 Vertex 13 2 7 Note that if you use PolyLib version 5 22 0 or newer then the output may look slightly different as the computation of the Rays may have been postponed to a later stage The program latte2polyl ib pl can be used to convert a polytope from LattE De Loera et al 2003 notation to PolyLib notation As an alternative to the constraints based input the input polytope may also be spec ified by its Ray matrix The first line of the input contains the single word vertices The second line contains the number of rows and the number of columns in the Ray matrix The rest of the input is composed of the elements of the matrix Recall that 36 the number of columns is two more than the number of variables where the extra first columns is one or zero depending on whether the ray is a line or not The next columns contain the numerators of the coordinates and the final column contains the denomina tor of the vertex or 0 for a ray E g the above set can also be described as vertices 3 0 barvinok enumerate The program barvinok enumerate enumerates a parametric polytope as a piecewise step polynomial or rational generating function It takes two polytopes in PolyLib notation as input optionally followed by a list of parameter names The two polytopes refer to arguments P and C of the corresponding function
120. o the facet We have 4 1 Bet 74 Therefore t 4 0 y b 1 0 Drs OLO 4 I b 1 0 2 0 1 eee and 1 1 1 k X x1 x2 5 X032 33132 12 12 With xy and x 4 the contribution of this facet is 2 1 115 2z7 zdz I sE The normal to the facet F in the margin is n 1 0 The vector f 0 1 is x parallel to the facet We have 1 ol i Therefore t 2 0 y Dj b i 1 05 enr Gap bc 0 b 2 0 Ha EE 0 1 1 0 2 1 0 2 ls E and 2 12 2 12 With xy 2 and x z the contribution of this facet is 5 1 33 zdz 1224 7 The normal to the facet F4 in the margin is n 1 1 The vector f 1 1 is parallel to the facet We have Therefore t 7 0 3Di 1D b i 1 0 Deun Buddy 200 J 1 Dipi EET 1 279 b 1 1b 2 1 b 2 0 152 0 162 0 vw 2 1722 1 1 1 1 h x Drr xix z Di xi 2 5 75 and 1 1 1 1 D D 1 24 X1X2 9120 242 24 s d h x 5 24 With x z 1 and z 1 the contribution of this facet is Nie 47 P 2 y T 1 7 1 7 766 z 24 4 at 24 5 3 2 The total contribution of the edges is therefore 115 33 47 355
121. of f and combining the lists of quasipoly nomials over shared domains join of m and q taking the sum over all elements in the intersection of the range of m and the domain of 41 q qQ m q 42 41 join of and f computing bound over all el ements in the intersection of the range of m and the domain of f and returning a list containing the bound and a boolean that is true if the bound is known to be tight l m f l m f l m f universal map with domain s and range 52 evaluate the piecewise quasipolynomial q in each element of the set s and return a piecewise quasipolynomial mapping each of the individual el ements to the resulting constant evaluate the piecewise quasipolynomial fold f in each element of the set s and return a piecewise quasipolynomial mapping each of the individual el ements to the resulting constant simplify s in the context of s2 i e compute the gist of s given 52 simplify m in the context of m i e compute the gist of m given m simplify m in the context of the domain s i e com pute the gist of m given domain s simplify q in the context of the domain s i e com pute the gist of q given s simplify f in the context of the domain s i e com pute the gist of fj given s continued on next page 19 Syntax Meaning Il cu Il SS qo h SE m lt 5 lt nm lt 52 lt Hog Sus
122. ok library Some of these operations are extensions of the functions from PolyLib with the same name Most of these operation are also provided by islonisl pw qpolynomials which are set to replace evalues Use isl_pw_qpolynomial_from_evalue to convert from evalues to isl_pw_qpolynomials __isl_give isl_pw_qpolynomial isl_pw_qpolynomial_from_evalue __isl_take isl space dim const evalue e void eadd const evalue el evalue res void emul const evalue el evalue res The functions eadd and emul takes two pointers to evalues e1 and res and com putes their sum and product respectively The result is stored in res overwriting and deallocating the original value of res It is an error if exactly one of the arguments of 28 eadd is of type partition unless the other argument is 0 The addition and multipli cation operations are described in Verdoolaege 2005 Section 4 5 1 and Verdoolaege 2005 Section 4 5 2 respectively The function eadd is an extension of the function new_eadd from Seghir 2002 Apart from supporting the additional types from Section 2 3 the new version also additionally imposes an order on the nesting of different enodes Without such an ordering evalues could be constructed representing for example 059 0x9 1x5 y 39 0 1y 5x which is just a funny way of saying 0 void eor evalue 1 evalue res The function eor implements the union operation from Verdoolaege 2005 Section 4 5
123. ommon divisor FINE nee Hermite Normal Form LOM vex least common multiple LLE Ies Lenstra Lenstra and Lovasz basis reduction algorithm PIP Parametric Integer Programming SNF Smith Normal Form 127 Index x 8 18 18 12 12 convert 39 41 direct H1 direction 3 disable fractional 26 enable incremental 26 explicit 39 1 941 93 help 36 isl 40 lower 42 omega 40 pip 40 series 39 H1 summation 42 version 36 gt 19 36 42 V 36 39 41 d 41 e 39 BIHA i K0 0 40l p 40 s 39 12 Bia lt lt lt lt lt lt 12 20 20 gt 20 E gt 20 gt gt 20 gt gt 21 8 9 20 I2 20 2 4112 91 access relation 10 a I2 L5 affine embedding 55 after 17 Aspinall D 119 720 Atkey R 119 720 B eler B 58 727 Baldoni V 77 82 20 Baldoni Silva M W 118 Banaszczyk W barvinok 1 B 6 13 221 26 23 35 36 availability Barvinok s decomposition Barvinok A Barvinok A I barvinok_bound barvinok count barvinok count with options barvinok ehrhart barvinok enumerate barvinok enumerate e barvinok enumerate e series 34 barvinok enumerate ev 33 barvinok enumerate i
124. on trick then ensures that you do not need to worry about lower dimensional faces in the decomposition Another way of avoiding the lower dimensional faces in the primal space is to perturb the vertex of the cone such that none of the lower dimensional face encountered contain any integer points K ppe 2007 In this section we describe another technique that is based on allowing some of the facets of the cone to be open The basic step in Barvinok s decomposition is to replace a d dimensional simple cone K pos u 0 by a signed sum of at most d cones K with a smaller determinant in absolute value The cones are obtained by successively replacing each generator of K by an appropriately chosen w 1 K uj U w 5 6 To see that we can use these K to perform a decomposition rearrange the such that for all 1 lt i lt k we have a lt 0 and for all k 1 lt i d we have a gt 0 with d d the number of zero We may assume k lt d otherwise replace w B by w B We have k d wrt 2 i 1 i k 1 or Y u 5 7 i 0 i k 1 47 Figure 5 11 Possible locations of w with respect to the rays of a 3 dimensional cone The figure shows a section of the cones with up w Bo 1 and B a gt 0 for lt i lt k Any two uj and u on the same side of the equality are on opposite sides of the linear hull
125. ons in the parameters that may be valid only in parts chambers of the parameter domain Since the volume computation is based on the active vertices we perform the computation in each chamber separately Also note that since the vertices are affine expressions it is easy to check whether they belong to a facet The volume of a d simplex i e a d dimensional polytope with d 1 vertices is relatively easy to compute In particular if v p for O lt i lt d are the parametric vertices of the simplex P then vu p voi p vi2 ph Vo2 P via p voa p 1 Va3 p voi p v22 P vo p voa p voa p vol P di et 9 5 25 vai P 7010 vao2 p vo2 P voa p If P is not a simplex i e N gt d 1 with N the number of vertices of P then the standard way of computing the volume of P is to first triangulate P i e subdivide P into simplices and then to compute and sum the volumes of the resulting simplices One way of computing a triangulation is to compute the barycenter of and to perform a subdivision by computing the convex hulls of the barycenter with each of the facets of P If a given facet of P is itself a simplex then this convex hull is also a simplex Otherwise the facet is further subdivided This recursive process terminates as every 1 dimensional polytope is a simplex The triangulation described above is known as the boundary triangulation B eler et al 20
126. ons using oriented matroids In M Joswig and N Takayama Eds Algebra Geometry and Software Sys tems pp 49 75 Springer 187 Pop S G A Silber A Cohen C Bastoul S Girbal and N Vasilache 2006 GRAPHITE Polyhedral analyses and optimizations for GCC Technical Report A 378 CRI Centre de Recherche en Informatique cole des Mines de Paris Fontainebleau France Contribution to the GNU Compilers Collection Devel opers Summit 2006 GCC Summit 06 Ottawa Canada June 28 30 2006 118 Preparata F P and M I Shamos 1985 August Computational Geometry An Introduction Monographs in Computer Science Springer 93 124 Pugh W 1994 Counting solutions to Presburger formulas How and why In SIGPLAN Conference on Programming Language Design and Implementation PLDI 94 pp 121 134 167 Rabl T 2006 January Volume calculation and estimation of parameterized inte ger polytopes Master s thesis Universitat Passau Ryan C T A X Jiang and K Leyton Brown 2010 Symmetric games with piecewise linear utilities In 10 Proceedings of the Behavioral and Quantitative Game Theory New York NY USA pp 1 1 ACM 119 Sakellariou R 1996 October On the Quest for Perfect Load Balance in Loop Based Parallel Computations Ph D thesis University of Manchester 1671 Sakellariou 1997 August Symbolic evaluation of sums for parallelising com pilers In A Sydow Ed Procee
127. p vo p we find h p x Daorn AO P p xeP p nZ f ACP vi p p v2 p The transverse cone of a polytope along the whole polytope is a zero dimensional cone in a zero dimensional space and so Dpw p p LP p P p D 1 The transverse cone along vi p is vi p R and so Dpp H i p R D as in 5 37 with y D 1 D and t vi p vi p vi p Similarly we find Dp p u v p R D as in 5 37 with D 1 D and t v2 p lv p v2 p Summarizing we find va p npxo xeP p NZ v p bn 1 V oy 2 Gr COOP bin L W y Dyp COP Note that in order to apply this formula we need to verify first that v p is indeed smaller than or equal to v2 p Since the combinatorial structure of P p does not change throughout the interior of the chamber we only need to check the order of the two vertices for one value of the parameters from the interior of the chamber a point which we may compute as in subsection 5 1 5 14 3 Summation over a two dimensional parametric polytope For two dimensional polytope formula has three kinds of contributions the integral of the polynomial over the polytope contributions along edges and contribu tions along vertices As suggested by Berline 2007 the integral can be computed by applying the Green Stokes theorem M ji 22 f Ldx pp O
128. p square part of H Example 5 59 Consider the cone This cone is already situated in the first quadrant but this not be obvious from the constraints Furthermore directly adding slack variables would lead to a total of 4 variables whereas we can also represent this cone in standard form with only 3 variables We have JM end 53 peeTp 1 0 1331 2875 52 53 31 57 TEER Adding a slack variable for the second row of H we obtain the equivalent problem 1331 2875 1 0 x gt 0 with _ 6 13 0 x _ 31 57 ul A similar construction was used by Eisenbrand 2000 2000 Lemma 3 10 and Hung and Rom 1990 5 18 Using TOPCOM to compute Chamber Decompositions In this section we describe how to use the correspondence between the regular triangu lations of a point set and the chambers of the Gale transform of the point set Gelfand et al 1994 to compute the chamber decomposition of a parametric polytope This correspondence was also used by Pfeifle and Rambau 2003 and Rambau 2003 Eisenschmidt and K ppe 2007 Let us first assume that the parametric polytope can be written as an 5 60 Ax lt b p where the right hand side b p is arbitrary and may depend on the parameters The first step is to add slack variables s to obtain the vector partition problem 15 b p 5 gt 0 with the identity matrix Then we compute the right kernel K of the matrix A
129. points to lie on a non standard lattice Let the effective parameter domain be given as p Ap c gt 0 where A Z of row rank d otherwise the effective parameter domain would contain a line Let H be the left HNF of A i e A HQ with H lower triangular with positive diagonal elements and Q unimodular Let be the matrix obtained from Q by reversing its rows and similarly H from H by 56 reversing the columns After performing the transformation i e p the transformed parameter domain is given by 1401 e gt 0 or p Hp c gt 0 The first constraint of this domain is 2 0 A ray with non zero final co ordinate therefore has a positive final coordinate Similarly the second constraint is hap hoi Pm c2 gt 0 A ray with zero nth coordinate but non zero n 1st coor dinate will therefore have a positive n 1st coordinate Continuing this reasoning we see that all rays in the transformed domain are revlex positive If the parameter domain does contains lines but is not restricted to a non standard lattice then the number of points in the parametric polytope is invariant over a trans lation along the lines It is therefore sufficient to compute the number of points in the orthogonal complement of the linear subspace spanned by the lines That is we apply a prior transformation that maps a reduced parameter domain to this subspace 11 r
130. pute the reduced bases See subsection 5 20 for more information 35 3 Applications included in the barvinok distribution This section describes some application programs provided by the barvinok library available from http freshmeat net projects barvinok For compilation instructions we refer to the README file included in the distribution Common option to all programs version V print version help list available options 3 1 barvinok count The program barvinok count enumerates a non parametric polytope It takes one polytope in PolyLib notation as input and prints the number of integer points in the polytope The PolyLib notation corresponds to the internal representation of Polyhedrons as explained in Section 2 T The first line of the input contains the num ber of rows and the number of columns in the Constraint matrix The rest of the input is composed of the elements of the matrix Recall that the number of columns is two more than the number of variables where the extra first columns is one or zero depending on whether the constraint is an inequality gt 0 or an equality 0 The next columns contain the coefficients of the variables and the final column contains the constant in the constraint E g the set 5 5 52 0 25 lt 13 from Verdoolaege 2005 Example 38 on page 134 corresponds to the following input and output cat S 23 11989 1 2 13 barvinok count S POLYHEDRON Dimension 1 Co
131. re 4 5 and 3 4 The corresponding widths are 2 4 5 We 6 Co 3 4 Wa 4 We similarly find C31 pos 4 5 1 2 V 0 with integer hull pos 4 5 1 2 1 1 shown in Figure 5 76 yielding 4 5 We 6 c4 1 2 we 3 1 1 we 3 On the other hand C41 9 as shown in Figure 5 77 and so this combination does not yield any width direction candidates The other pairs of vertices further yield c6 1 2 we 3 3 5 We 5 3 4 we 4 3 5 We 5 2 3 Wep 3 Since the polytope under consideration is not parametric there is only one non empty O dimensional chamber and it corresponds to one of the directions say 1 2 with width 3 the other directions with the same width having been removed Each of the three directions that yield the minimal width of 3 is shown in Fig ure 5 78 Example 5 79 Consider the polytope 2x p 520 2x p 5 20 4 x 2x p 520 2x pt 52 gt 0 from Woods 2004 Example 2 1 7 The parametric vertices are _ 5 p 5 239 5 p vp o EE 5 v3 p 2 va p S JE Figure 5 78 A polytope and its lattice width directions We find two essentially different candidate width directions 0 1 wWe p 5 p 1 0 We p 5 p The first direction can
132. re 4 33 page 186 In the difference D there will be exactly one value of u for 105 each value of the remaining variables for which there was at least one value of u in P 5 0 du st u eP Clu s t u D The function c s can then be computed by counting the number of elements in D s These operations can be performed either in the space of unions of parametric poly topes or on generating functions In the first case D s can be written as a disjoint union of parametric polytopes that can be enumerated separately In the second case we first compute the generating function f x y of the set 5 s t Ju Z s t u P and then obtain the generating function C x of c s as C x f x 1 In the remainder of this section we will concentrate on the computation of the generating function of S To compute this generating function in the current case where there is only one existentially quantified variable we first compute the generating function g x y z of P s t u perform operations on the generating function equivalent to the set opera tions see e g Verdoolaege 2005 Section 4 5 3 resulting in a generating function g X z and then sum over all values at most one for each value of s and t of u i e f x y g e y 1 If there is more than one existentially quantified variable m gt 1 then we can in principle apply the above shifting and subtracting technique recursively to obtai
133. representing arbitrary precision integers ZZ as well as vectors vec ZZ and matrices mat ZZ of such integers We further introduces a type QQ for representing a rational number and use vectors 00 of such numbers struct QQ ZZ n ZZ NTL vector decl QQ vec Q0 Each term in a rational generating function is represented by a short rat structure struct short rat struct rows terms in numerator vec QQ coeff mat ZZ power ln struct rows factors in denominator mat ZZ power 1d The fields and represent the numerator and the denominator respectively Note that in our implementation we combine terms with the same denominator In the numerator 30 short_rat n coeff 3 2 2 1 n power 2 3 5 7 d power 1 3 0 2 Figure 2 6 Representation of 3 xix 2 xaT a 1 x each element of coeff and each row of power represents a single such term The vector coeff contains the rational coefficients of each term The columns of power correspond to the powers of the variables In the denominator each row of power corresponds to the power b of a factor in the denominator Example 2 5 Figure shows the internal representation of 3 2 3 5 7 5 2 xf The whole rational generating function is represented by a gen_fun structure typedef std set lt short_rat short_rat_lex_
134. rom the domain of m to the domain of of those elements such that their images live in the same space and such that the images of the elements of m are lexicographically smaller than those of map from s to s2 of those elements that live in the same space and such that the elements of s are lexicographically strictly greater than those of 52 a map from the domain of m to the domain of m of those elements such that their images live in the same space and such that the images of the el ements of m are lexicographically strictly greater than those of mp continued on next page 20 Syntax Meaning m 5 gt gt S2 gt gt map from s to s2 of those elements that live in the same space and such that the elements of s are lexicographically greater than those of 52 a map from the domain of m to the domain of of those elements such that their images live in the same space and such that the images of the elements of m are lexicographically greater than those of 21 2 PolyLib interface of the barvinok library obsoles cent Although barvinok currently still uses PolyLib internally this is likely to change in the not too distant future Consider using 151 based alternatives for the functions in this section as the latter are likely to be removed in future releases Our barvinok library is built on top of PolyLib Loechner 1999 In particular it reuses the
135. rwise we can apply the unimodular transformation bibis ub y m without changing the width in direction x If P contains x1 y1 it intersects x x in a segment a2 We may then similarly assume that a2 gt P will only cut x 2 in points with y coordinate smaller than 2 a2 The width in the y direction will therefore be smaller than 2 o lt 2 contradicting that x is a 108 Figure 5 88 Lattice point free polygon with lattice width 2 lattice width direction If P does not contain then it only intersects x x in points with y coordinate with 0 lt lt 1 Given any such point it is clear that P intersects x x 2 only in points with y coordinate strictly between y and y 1 a again showing that the width in the y direction is smaller than 2 and leading to the same contradiction The contradiction shows that there can be no gaps of width 2 in the lattice width direction of P 0 Note that the w 2 bound is too coarse to reach the above conclusion as w 2 gt 2 An example of a polygon with lattice with greater than 2 is the polygon with vertices 17 110 83 110 2 10 9 10 and 177 90 100 90 shown in Figure 5 88 which has width 221 110 The idea of the projection algorithm is now to first find a direction in which the gaps are expected to be small and to unimodularly transform the existentially
136. se of rational coefficients CY Sa pai 1X 4 d yi Ir aen 1 i 0 ja i i 0 with aj 61 Example 5 30 Consider the rational generating function 2 2 xi X 1 Q ax D xjbe xD0 xx x from Verdoolaege 2005 Example 39 Since this is a 2 dimensional problem we first compute the first Todd polynomials evaluated at 1 f T x 1 1 1 1 ET md pepe 6 3 1 t 2 6 6 and 1 3 l e e 1 II 6 36 where we represent each truncated power series by a vector of its coefficients The vector A 1 1 is not orthogonal to any of the rays so we can use the substitution x eo and obtain e 1 1 1 e7 1 e 1 e P 1 e e We have to 3 3 1 1 6 36 2 _ 6 l e 1 6 36 t 3 3 1 1 6 36 2t _ 6 1 e 1 6 36l The first term in the rational generating function evaluates to 1 2 41 4 3 IE 6 E jolt d 2 1 6 36 11 6 36 2 9 21 1 21 Lb 6 36 1 213 71 sc 2 6 ded 1 ELE zl 6 8 9 33 Due to symmetry the second term evaluates to the same value while for the third term we find 1 Z i 0 d 1 3 9e Oe ea The sum is um 24 24 12 62 Note that the run time complexities of polynomial and exponential substitution are basically the same The experiments of K ppe 2007 are somewhat misleading in this r
137. sl 33 40 barvinok_enumerate_pip 40 barvinok_enumerate_scarf 34 barvinok enumerate scarf 34 128 barvinok enumerate with options 33 barvinok_options 24 25 33 barvinok options new with defaults barvinok_series 34 barvinok series with options 34 56 barvinok_summate 2 41 barvinok_union barycenter 45 basis reduced see reduced basis Bastoul C 118 724 Beck M 118 720 before beneath beyond Berline N 68 69 69 70 72 118 720 bernoulli Bernoulli number Bernoulli polynomial Bernoulli polynomials Beyls K 117 118 727 251127 big parameter Bik A J C 56 727 Boulet P 13 27 box Brion s polarization trick 47 Brion M Bruynooghe M 127 card 12 L5 Catthoor F 117 26 cdd 92 chamber decomposition characteristic cone Chernikova Clauss P Z2 117 27 coalesce Cochet C 118 720 codegen 12 15 coeff BI coefficients Cohen A 118 24 Cohen J 58 27 compressed parameter 32 compute evalue compute poly 29 cone simple see simple cone configure constituent Constraint 2 36 context 37 Cook W 35 92 93 727 Cools R 117 726 Craciun G 119 723 cross 12 15 cutting plane 41 a 23 D Hollander E 18 727 De Loera J A 36 45 47 5860 63 2 7 deltas 2 deltas_map 15 Devos 117 Dimension 22 Disjunctive Normal Form DNF
138. smaller_denominator gt short_rat_list struct gen_fun short_rat_list term Polyhedron context void add const QQ amp c const vec ZZ amp num const mat ZZ amp den void add short rat r void add const QQ amp c const gen fun gf barvinok options options void substitute Matrix CP gen fun Hadamard product const gen fun gf barvinok options options void print std ostream amp os unsigned int nparam char param name const operator evalue const ZZ coefficient Value params barvinok options options const void coefficient Value params Value c const gen fun Polyhedron C gen fun Value c 31 gen fun const gen fun gf gen_fun 13 new gen fun can be constructed either as empty rational generating function possi bly with a given context C as a copy of an existing rational generating function gf or as constant rational generating function with value for the constant term specified by c The first gen fun add method adds a new term to the rational generating function described by the coefficient c the numerator num and the denominator den It makes all powers in the denominator lexico positive orders them in lexicographical order and inserts the new term in term according to the lexicographical order of the combined powers in the denominator The second gen fun add method adds c times gf to the rational generating function The method gen fun operator evalue performs
139. so among the bj this can be simplified to checking whether there exists a rational solution for a to J ajb witha gt Oand aj2 I iti iti Example 5 69 Consider the cone pos 2 3 3 4 shown in Figure 5 68 The Hilbert basis of this cone is 0 0 2 3 3 4 C1 D 0 1 1 0 We have 1 0 3 1 D 1 D while 1 1 and 1 1 can not be written as overconvex combinations of the other b 0 The vertices of the integer hull of C 0 are therefore Q2 3 3 4 1 1 CL 1 96 T T Figure 5 70 The integer hull of a truncated cone 5 22 2 Using generalized basis reduction Another way of computing the integer hull of a truncated cone is to apply the method of In this case the initial set of points will consist of the smallest integer representatives of the extremal rays of the cone together with the extremal rays themselves That is if C pos r with r Z4 then our initial approximation of the integer hull of C 0 is conv pos rj Furthermore we need never consider any of the bounding constraints that are also bounding constraints of the original cone When optimizing along the normal of any of the other facets we can take the lower bound to be 1 This will ensure that the origin is excluded without excluding any other integer points Example 5 71 Consider once more the cone pos 2 3 3 4
140. solution is 0 01s so very short execution times are not measured very accurately resulting in some anomalies in the graphs In the non parametric tri angle experiment 5 54 the new Laurent methods is by far the best For the parametric triangle Euler Maclaurin still beats Laurent but the new version is getting much closer than the old version 5 17 Some algorithms or tools expect a polyhedron to be specified in standard form i e Conversion to standard form Ax b 5 56 x 2 0 Given an arbitrary parametric polyhedron x Ax b p 0 5 57 a conversion to standard form requires the introduction of slack variables and a way of dealing with variables of unrestricted sign In this section we will be satisfied with a 85 28 45 ga380232 1000 T T T T T T x a z x 100 PR x o x a x 10 F J 2 x tom s 2 x ir a v d F 0 1 F 1 x K 0 01 reduction to unweighted Barvinok 7 Euler Maclaurin T Laurent expansion old ddoi r Laurent expansion new 0 2 4 6 8 10 12 14 16 18 Figure 5 55 Execution times for summing a monomial over a difficult parametric tri angle reduction to the form Ax b 5 58 Dx gt with D a diagonal matrix with positive entries That is we do not necessarily make all variables non negative but we do ensure that they have a lower bound If needed a subsequent redu
141. stic cone of P Adding a positive multiple of the sum of the extremal rays of C to the barycenter 1 em Y N i P of P where N is the number of vertices results in a point in the interior of P 5 2 The integer points in the fundamental parallelepiped of a sim ple cone This section is based on Barvinok 1992 Lemma 5 1 and De Loera and K ppe 2006 In this section we will deal exclusively with simple cones i e d dimensional cones with d extremal rays and d facets Some of the facets of these cones may be open Since we will mostly be dealing with cones in their explicit representation we will have occasion to speak of open rays by which we will mean that the facet not containing the ray is open There is only one such facet because the cone is simple Definition 5 1 Fundamental parallelepiped Let v pos u bea closed shifted cone then the fundamental parallelepiped II of K is Yaw lose 1 If some of the rays uj of are open then the constraints on the corresponding coeffi cient are such that lt a lt 1 Lemma 5 2 Integer points in the fundamental parallelepiped of a simple cone Let u bea closed simple cone and let A be the matrix with the generators uj of K as rows Furthermore let VAW S diags be the Smith Normal Form SNF of A Then the integer points in the fundamental parallelepiped of K are given by w 5 3 T
142. straints this access relation can be represented as t F i j gt t j 1 10 It is however customary to include the constraints on the iterators in the access relation resulting in n gt F i j t j 1 9 lt i lt n and 0 lt j lt n i The constraints can be added by intersecting the domains of the access relations with the iteration domains For example the following sequence constructs the access rela tion for all reads in the program D n gt T i 9 lt i lt n F i j 0 lt i lt n and 0 lt j lt n i B i 0 lt i lt n Read T i gt a i F i j gt t jl F i j gt t j 1 B i gt t i Read Read In this sequence the operator when applied to a and a set intersects the domain of the map with the set The notation RCS can be used to compute the image of the set S under the map For example given the sequence above 1 computes the array ele ments read in iteration 0 1 of statement F and is equal to n gt t 2 n gt 2 t 1 n 2 That is elements 1 and 2 of the t array are read provided n is at least 2 Maps need not be single valued As an example assume that A is a two dimensional array of size n in both dimensions Iterations i of a statement S may pass a pointer to an entire row of A to a function as in A i Without knowing anything about f we would have to assume that this functio
143. subtract once in the remaining horizontal direction In fact for two dimensional polytopes the gaps in the lattice width direction will always be one as shown by the following lemma Lemma 5 87 For any rational polygon the gaps in a lattice width direction are of width at most 1 Proof We may assume that x is the given lattice width direction of a given polygon P If there is a gap of width 2 then there is an integer value x of x such that PN x1 y 0 PA Gg 2 0 while Po x 1 N Z 0 Using Barvinok and Woods 2003 Lemma 4 2 we can put a scaled down P of P between x x and x x 2 and inside of P P meets the line x x 1 between two consecutive integer points y and y 1 Let P be the polygon bounded by x x and x x1 2 and two lines that separate P from these two integer points P will have the same width 2 in the x direction while P c P The x direction is therefore also a lattice width direction of P P cannot intersect both x x and x x 2 in a segment of length greater than or equal to 1 Otherwise it would also intersect x x lina segment of length greater than or equal to 1 We may therefore assume that the length of the intersection of P with x x is smaller than 1 If this line segment contains an integer point then call it y2 Otherwise let y be the greatest integer smaller than the points in the line segment We may assume that y y2 Othe
144. t We are now given a polyhedron P such that the lattice width direction of P is c We first extend c to an mxm unimodular matrix U using the technique of subsection 5 7 c i7 and then compute I 0 0 P 0 n O R 0 0 U We have T s t du Z st u e P i e we may have changed the values of the existentially quantified variables but we have not changed the set T Now consider the set T s t 305 u 2 situ e P This set has only m 1 existentially quantified variables so we may apply this projec tion algorithm recursively and obtain the generating function f x z for T The set T may no longer correspond to the integer points in a polytope but by construction the gaps in the final coordinate are small lt w m By now we have a generating function f x z for the set T with small gaps in the final coordinate and we have to compute the generating function f x y for T By computing Lw m yd fa y 2 CD f y 2 5 89 k 1 110 where represents the operation on generating functions that corresponds to set dif ference on the corresponding sets we obtain a generating for the set 7 where only the smallest value of uj is retained The total number of us associated to any s t is therefore either zero or one and so the multiset defined by taking as many copies of s t as there are associated values of is actually the set T That is
145. t 0 Lattice 1 1 2 8 0 0 1 2 5 1 2 t 1 1 0 1 2 s 1 2 t 1 2 5 6 Left inverse of an affine embedding We often map a polytope onto a lower dimensional space to remove possible equalities in the polytope These maps are typically represented by the inverse mapping the co ordinates x of the lower dimensional space to the coordinates x of an affine subspace of the original space 1 x T v x 1 0 1 where as usual in PolyLib we work with homogeneous coordinates To obtain the transformation that maps the coordinates of the original space to the coordinates of the lower dimensional space we need to compute the left inverse of the above affine embedding i e an A b and d such that x 1 To compute this left inverse we first compute the right Hermite Normal Form HNF of T JH 07 10 U dH v x 41 0 d b The left inverse is then simply 0 We often also want description of affine subspace that is range of the affine embedding and this is given by X lo U v 1 07 1 This computation is implemented left_inverse 55 5 7 Integral basis of the orthogonal complement of a linear sub space Let Z be a basis of a linear subspace We first extend with zero rows to obtain a square matrix M and then compute the left HNF of M _ H
146. t j lt i is equal to n gt a b 1 lt lt n and 1 lt b lt but it is not equal to n gt j i 1 lt i lt and 1 lt j lt because the order of the coordinates has changed or n gt i j 1 lt i lt m and 1 lt j lt because it involves a different parameter It is sometimes convenient to represent constraints that only involve parameters and that are not tied to any particular space To construct such a parameter domain the list of coordinates should simply be omitted Note that the colon is required in this case even if there are no constraints In particular icq represents the universal parameter domain which is very different from the empty set To plug in a particular value for a parameter the user should take the intersection with a parameter domain assigns a particular value to the parameter For example S n gt t i j 1 lt lt and 1 lt j lt i S n gt n 5 It should be noted though that the result is not the same as simply replacing n by 5 as the result of the above sequence will still have the global parameter n set to 5 To avoid this assignment the user should instead compute the gist of the original set in the context of setting n to 5 That is the result of the sequence below is True S1 i j 1 lt lt 5 and 1 lt j lt S2 n gt i j 1 lt lt and
147. tations are nearly identical and the final result is simply the sum of the two generating functions found for each of the two disjoint chambers h x z g G z x 116 6 Publications 6 1 Publications about the Library This is a list of some reports and publications explaining details of parts of the barvinok library Analytical computation of Ehrhart polynomials and its applications for embed ded systems Verdoolaege Beyls Bruynooghe Seghir and Loechner 2004b Analytical computation of Ehrhart polynomials and its applications for embed ded systems Verdoolaege Beyls Bruynooghe Seghir and Loechner 2004c Analytical Computation of Ehrhart Polynomials and its Application in Compile Time Generated Cache Hints Seghir Verdoolaege Beyls and Loechner 2004 Analytical computation of Ehrhart polynomials Enabling more compiler analy optimizations Verdoolaege Seghir Beyls Loechner and Bruynooghe 2004d Experiences with enumeration of integer projections of parametric polytopes Verdoolaege Beyls Bruynooghe and Catthoor 2004a Experiences with enumeration of integer projections of parametric polytopes Verdoolaege Beyls Bruynooghe and Catthoor 2005a Computation and Manipulation of Enumerators of Integer Projections of Para metric Polytopes Verdoolaege Woods Bruynooghe and Cools 2005b Incremental Loop Transformations and Enumeration of Parametric Sets Ver doolaege 200
148. the conversion from rational generating function to piecewise step polynomial explained in Verdoolaege 2005 Sec tion 4 5 5 The Polyhedron context is the superset of all points where the enumera tor is non zero used during this conversion i e it is the set from Verdoolaege 2005 Equation 4 31 If context is NULL the maximal allowed context is assumed i e the maximal region with lexico positive rays The method gen fun coefficient computes the coefficient of the term with power given by params and stores the result in c This method performs essentially the same computations as gen fun evalue except that it adds extra equality constraints based on the specified values for the power The method gen fun substitute performs the monomial substitution specified by the homogeneous matrix CP that maps a set of compressed parameters to the original set of parameters That is if we are given a rational generating func tion G z that encodes the explicit function g i where i are the coordinates of the transformed space and CP represents the map i Ai a back to the original space with coordinates i then this method transforms the rational generating function to F x encoding the same explicit function f W i e f f Ai a gi This means that the coefficient of the term x4i 4 in F x should be equal to the coefficient of the term z in G z In other words if
149. there is exactly one variable with a negative exponent where f is such that if is the first non zero coefficient of b and the sum of all exponents is 1 Let n contain the non negative exponents then the coefficient of such a term is e 5 51 poem J In the second kind of terms all exponents are non negative and the coefficient is given by bA Ying sp 5 52 n 7 n gt 0 Given as before a polynomial q p Bm P y that we wish to sum over the integer points of a polytope P we perform the following operations for each unimodular cone in the decomposition of each vertex cone e For each m with Bm p 0 83 Consider all decomposition nj such that each n corresponds to one of the two kinds of terms i e either lt 0 and njg 1 or nj gt 0 For each such decomposition compute the coefficient oq of y as the product of the corresponding coefficients from 5 51 and 5 52 e The contribution of the cone is the sum of m a m over all these m Example 5 53 As in Example let us compute 2 2 yeroz with T a triangle with generating function 2 Xx X 1 ag 0 259359 a fax Consider the first term As before we write the exponent of the numerator of this term as 2 0 2 1 0 O 1 1 and so we obtain 1 1 eC 9 0 1 1 evi ty2 1L
150. this point need not be integer A direct application of the above algorithm using any of the triangulations would yield for each chamber a volume expressed as the sum of the absolute values of poly nomials in the parameters To remove the absolute value we plug in a particular value of the parameters not necessarily integer belonging to the given chamber for which we know that the volume is non zero Again it is sufficient to take any point in the interior of the chamber The sign of the resulting value then determines the sign of the whole polynomial since polynomials are continuous functions and will not change sign without passing through zero 5 10 Maclaurin series division If P t and Q t are two Maclaurin series P t ao ait art Tee Q f bo bit ba then as outlined by Henrici 1974 241 247 we can compute the coefficients c in Po Q by applying the recurrence relation L 1 bo a gt i 1 opto To avoid dealing with denominators we also compute d pte instead as 1 bi bicis i l The coefficients c can then be directly read off as 1 1 by 5 11 Specialization through exponential substitution This section draws heavily from De Loera and K ppe 2006 59 We define a short rational generating function to be a function of the form fe Be 5 26 iel se xP with x C7 a Q wi Z a
151. total contribution of the vertices is then 1355 M 1067 E 253 6l 288 288 144 6 and the total sum is 121 355 61 24 24 6 Example 5 45 Consider the parametric polytope 30 2A3xi lt 9 4 lt x lt 5 Ifn 2 3 then the vertices of this polytope are 2 4 2 5 3 n 3 4 and 3 3 5 The contributions of the faces of P n to 2 m xeP n nZ2 for the chamber n gt 3 are shown in Table 5 The final result is 2 2026 45 ine320 5 15 Summation through exponential substitution and Laurent ex pansions old This section was inspired by Baldoni et al 2008 and presents the old implemena tion of the laurent summation method The implementation is described in subsec tion 5 16 Let f x be the generating function of a polytope P i e f x 2 xt tePNZ Substituting x we obtain a yy LE n J _ fe Au gt 2 t n tePNZ4 tePnzZ4 gt 0 gt 0 tePnZ 77 E N 3 n 3 E ILI A 3 n 3 5 s xd 3 n 3 4 n 36 3 3 AN 63 2 57 n 2 3 4 4 3 4 3 8 23 23 115 2167127724 31 5 30 4 155 216 337 24 31 31 217 589 n 36 13 72 24 3 24 13 144 341 144 253 144 23 3 23 23 2 161 2 437 3 72 24 3 24 13 144 Table 3 Contributions of the faces of P n to t
152. uler Maclaurin formula 68 Loechner V 221D4 D6 37 54 105 117 Louichi A 118 24 Lov sz L 22 lower bound 1 1 16 MacKenzie K 119 720 Makhorin A 92 724 Malkin 23 Maslov V 23 mat ZZ 30 MaxRays 26 33 McDiarmid C 27 Meister B 32 57 LLL7 24 minimizing I8 moment curve 105 monomial substitution 32 multinomial coefficient n NbBid 22 NbConstraints NbEq 22 NbRays 22 neighborhood complex new eadd 29 22 Nikolov 118 26 NTL 26 30 Omega open facet 45 open ray optimization oracle Pajor A 720 parallelepiped fundamental see fundamental paral lelepiped parameter compression Parametric Integer Programming PIP 34 parametric polytope 98 params 16 parse_file partition 26 28H30 period 53 periodic 23 26 BA periodic number 53 pet 16 Pfeifle J 87 724 piplib 41 points2triangs polar cone 98 poly 16 Polyhedron 22 32 35 36 Polyhedron_Enumerate polyhedron integer hull 2 42 Polyhedron_Polarize Polyhedron Project 35 Polyhedron Reduced Basis 35 41 Polyhedron Sample 35 41 polyhedron sample I 41 131 PolyhedronIncludes PolyLib I B 22 26 28 29 33 3518377 version 5 22 0 or newer polymake polynomial 23 30 polynomial substitution 60 polytope lattice width polytope minimize 2 42 1
153. values any of the resulting open facets a p c gt 0 with a Z and c Z can then be replaced by a p 1 gt 0 Again we have sie 2 800 2 gt j 2 After this reduction the technique of Berline and Vergne 2006 can be applied practically verbatim to the parametric polytope with a fixed combinatorial structure In principle we could also handle piecewise quasi polynomials using the technique of Section 4 5 4 except that we only need to create an extra variable for each distinct floor expression in a monomial rather than for each occurrence of a floor expression in a monomial However since we currently only support two dimensional polytopes this reduction has not been implemented yet 5 14 2 Summation over a one dimensional parametric polytope The basis for the summation technique is the local Euler Maclaurin formula Berline and Vergne 2006 Theorem 26 ww _ Dior 5 36 F p eF P p where is a parametric polytope A is a lattice F P p are the faces of P p Dp p r p 18 a specific differential operator associated to the face of polytope The Lebesgue measure used in the integral is such that the integral of the indicator function of a lattice element of the lattice A N aff F p F p is 1 i e the intersection of A with the linear subspace parallel to the affine hull of the face F p Note that the original theorem is
154. ven a piecewise step polynomial in isl format the program barvinok summate computes the sum of the piecewise quasi polynomial evaluated in all integer values of the variables The result is an expression in the parameters Note that barvinok enumerate 41 and barvinok_enumerate_e can produce piecewise step polynomials when given the I option but they will have only parameters and no variables For example gt cat square p3 pwgp n gt x y gt x y n gt 9 3x and x gt 2 and y gt 4 and y lt 5 gt barvinok summate square p3 pwqp n gt 45 63 2 3 9 2 EQD 3 12 n gt 3 3 Options summation specifies which summation method to use box refers to the method of ACE e bernoulli refers to the method of euler refers to the method of and laurent refers to the method of subsection 5 16 3 10 barvinok bound Given a piecewise step polynomial in 151 format the program barvinok bound com putes an upper bound or lower bound for the values attained by the piecewise quasi polynomial over all integer values of the variables The result is an expression in the parameters Note that barvinok enumerate and barvinok enumerate e can pro duce piecewise step polynomials when given the I option but they will have only parameters and no variables gt cat devos pwqp U gt V gt CCXG1 3 U 2 3 V WU 2V 3 2V gt 3 U and 2V lt U and U gt 8 and U lt 10
155. with options Polyhedron unsigned exist unsigned nparam struct barvinok options options evalue barvinok enumerate isl Polyhedron P unsigned exist unsigned nparam struct barvinok options options evalue barvinok enumerate scarf Polyhedron P unsigned exist unsigned nparam struct barvinok options options The first function tries the simplification rules from Verdoolaege 2005 Section 4 6 2 before resorting to the method based on Parametric Integer Programming PIP from Section 4 6 3 The second function immediately applies the tech Section 4 6 3 The argument exist refers to the num ber of existential variables whereas the argument nparam refers to the number of pa rameters The order of the dimensions in P is counted variables first then existential variables and finally the parameters The function barvinok_enumerate_scarf per forms the same computation as the function barvinok enumerate scarf series below but produces an explicit representation instead of a generating function gen fun barvinok series Polyhedron P Polyhedron C unsigned MaxRays gen fun barvinok series with options Polyhedron P Polyhedron C barvinok options options gen fun barvinok enumerate e series Polyhedron P unsigned exist unsigned nparam barvinok options options gen fun barvinok enumerate scarf series Polyhedron P unsigned exist unsigned nparam barvinok options options The function barvinok series
156. x 70 In particular if M p x is such that 2M y y then Ji T h p x y f MOEDA P p Care must be taken to integrate over the boundary the positive direction Assuming the vertices of the polygon are not given in a predetermined order we can check the correct orientation of the vertices of each edge individually Let be the inner normal of a facet and let v p and v2 p be the two vertices of the facet then the vertices are in the correct order if vai p via p mil vo p vio p Since these two vertices belong to the same edge their order will not change within a chamber and so we can again perform this check for a single value of the parameters To integrate M over an edge F let f be a primitive integer vector in the direction of the edge Then v2 p vi p p f and any point on the edge can be written as v p Af with 0 lt A lt k p That is k p M p G y dy M p i1 Afi vi2 p Af dA 5 39 For the edges we can again apply 5 37 but we must first project the supporting cone at the edge into the linear subspace orthogonal to the edge Let n be the primitive integer inner normal of this facet F p then f 75 71 is parallel to the facet and we can write one of the vertices v p as a linear combination of these two vectors vip n apy 77 MES 5 40 or E p 7 O ML 54
157. x z Note that we are talking about equal ity of rational functions here Any further application of the set difference operation will not change this rational function but it will typically produce a different more complex representation To check whether the current k is sufficient we are going to count how many times any element of still appears in a shifted copy of with shift greater than or equal to k 1 If this number is zero any further set difference will have no effect That is we want to compute P r N easi T gt l k 1 or in terms of generating functions hay d 2 2 2 l k 1 We should point out here that while the Hadamard product corresponds to intersec tion when applied to generator functions of indicator functions i e with coefficients one or zero in general it will produce a generating function with coefficients that are the product of the corresponding coefficients in the two operands We can therefore 111 rewrite the above equation as Key l k 1 r sa F draya l k 1 ZU f x z s f X y z _ Computing h x y 1 would give us a generating function with as coefficients how many times a point of T still appears in a shifted copy of T for any particular value of s and t However we want to know if this number is zero for all values of s and t so we compute 1 1 1 instead We have to be
158. y R N S Without loss of generality we can therefore assume for the purpose of showing that R N 5 is non empty that the u indeed form an orthogonal basis In the orthogonal basis we have b u and the corresponding inward normal is either u or u Furthermore each normal of a facet of S of the first type is of the form fi bu j with ax by gt O and ix lt jx while for the second type each normal is of the form fi a uj bu with aj by gt 0 If fij b uj is the normal of a facet of 5 then either N Nj Ur uj or Ng Nj u uj Otherwise the facet would not cut Similarly if fij j b uj is the normal of a facet of S then either uj uj or N Nj u uj Assume now that RNS is empty then there exist 44 gt 0 not all zero such that gt Y u N 0 Assume gt 0 for some facet of the first type If N uj then b can only be canceled by another facet k of the first type with j ip but then also uj Since the j are strictly increasing this sequence has to stop with a strictly positive coefficient for the largest u in this sequence If on the other hand uj then ag only be canceled by the normal of a facet k of the second kind with i jy but then u and we return to the first case Finally if 4 gt 0 only for normals of facets of the second type then
159. y matrix then our problem is of the form and we already know how to solve this problem Note that again saturating any of the transformed constraints in X is equivalent to saturating the corresponding constraint in x We therefore only need to compute A U for the computation of the kernel K To construct the parametric vertices in the original coordinate system we can simply use the original constraints The same reasoning holds if H is any diagonal matrix since we can virtually replace Hx by x without affecting the non negativity of the variables If H is not diagonal then we can introduce new constraints x 2 d p where d p is some symbolic constant These constraints do not remove any solutions since each row in H expresses that the corresponding variable is greater than or equal to a non negative combination of the previous variables plus some constant We can then proceed as before However to reduce unnecessary computations we may remove from K the rows that correspond to these new rows Any solution saturating the new constraint would also saturate the corresponding constraint h and all the constraints corresponding to the non zero entries in h If a chamber contains a vertex obtained by saturating such a new constraint it would appear multiple times in the same chamber each time combined with different constraints from the above set Furthermore there would also be another as it turns out identical chamber where the vertex is on

barvinok: User Guide

Contents

Download Pdf Manuals

Related Search

Related Contents