1 \section{Internal Representation of the
\protect\ai[\tt]{barvinok
} library
}
3 Our
\barvinok/ library is built on top of
\PolyLib/
4 \shortcite{Wilde1993,Loechner1999
}.
5 In particular, it reuses the implementations
7 \shortciteN{Loechner97parameterized
}
8 for computing parametric vertices
10 \shortciteN{Clauss1998parametric
}
11 for computing chamber decompositions.
12 Initially, our library was meant to be a replacement
13 for the algorithm of
\shortciteN{Clauss1998parametric
},
14 also implemented in
\PolyLib/, for computing quasi-polynomials.
15 To ease the transition of application programs we
16 tried to reuse the existing data structures as much as possible.
18 \subsection{Existing Data Structures
}
21 Inside
\PolyLib/ integer values are represented by the
22 \ai[\tt]{Value
} data type.
23 Depending on a configure option, the data type may
24 either by a
32-bit integer, a
64-bit integer
25 or an arbitrary precision integer using
\ai[\tt]{GMP
}.
26 The
\barvinok/ library requires that
\PolyLib/ is compiled
27 with support for arbitrary precision integers.
29 The basic structure for representing (unions of) polyhedra is a
32 typedef struct polyhedron
{
33 unsigned Dimension, NbConstraints, NbRays, NbEq, NbBid;
38 struct polyhedron *next;
41 The attribute
\ai[\tt]{Dimension
} is the dimension
42 of the ambient space, i.e., the number of variables.
43 The attributes
\ai[\tt]{Constraint
}
44 and
\ai[\tt]{Ray
} point to two-dimensional arrays
45 of constraints and generators, respectively.
46 The number of rows is stored in
47 \ai[\tt]{NbConstraints
} and
48 \ai[\tt]{NbRays
}, respectively.
49 The number of columns in both arrays is equal
50 to
\verb!
1+Dimension+
1!.
51 The first column of
\ai[\tt]{Constraint
} is either
52 $
0$ or $
1$ depending on whether the constraint
53 is an equality ($
0$) or an inequality ($
1$).
54 The number of equalities is stored in
\ai[\tt]{NbEq
}.
55 If the constraint is $
\sp a x + c
\ge 0$, then
56 the next columns contain the coefficients $a_i$
57 and the final column contains the constant $c$.
58 The first column of
\ai[\tt]{Ray
} is either
59 $
0$ or $
1$ depending on whether the generator
60 is a line ($
0$) or a vertex or ray ($
1$).
61 The number of lines is stored in
\ai[\tt]{NbBid
}.
62 Let $d$ be the
\ac{lcm
} of the denominators of the coordinates
63 of a vertex $
\vec v$, then the next columns contain
64 $d v_i$ and the final column contains $d$.
65 For a ray, the final column contains $
0$.
66 The field
\ai[\tt]{next
} points to the next polyhedron in
67 the union of polyhedra.
68 It is
\verb+
0+ if this is the last (or only) polyhedron in the union.
69 For more information on this structure, we refer to
\shortciteN{Wilde1993
}.
71 Quasi-polynomials are represented using the
72 \ai[\tt]{evalue
} and
\ai[\tt]{enode
} structures.
74 typedef enum
{ polynomial, periodic, evector
} enode_type;
76 typedef struct _evalue
{
77 Value d; /* denominator */
79 Value n; /* numerator (if denominator !=
0) */
80 struct _enode *p; /* pointer (if denominator ==
0) */
84 typedef struct _enode
{
85 enode_type type; /* polynomial or periodic or evector */
86 int size; /* number of attached pointers */
87 int pos; /* parameter position */
88 evalue arr
[1]; /* array of rational/pointer */
91 If the field
\ai[\tt]{d
} of an
\ai[\tt]{evalue
} is zero, then
92 the
\ai[\tt]{evalue
} is a placeholder for a pointer to
93 an
\ai[\tt]{enode
}, stored in
\ai[\tt]{x.p
}.
94 Otherwise, the
\ai[\tt]{evalue
} is a rational number with
95 numerator
\ai[\tt]{x.n
} and denominator
\ai[\tt]{d
}.
96 An
\ai[\tt]{enode
} is either a
\ai[\tt]{polynomial
}
97 or a
\ai[\tt]{periodic
}, depending on the value
99 The length of the array
\ai[\tt]{arr
} is stored in
\ai[\tt]{size
}.
100 For a
\ai[\tt]{polynomial
},
\ai[\tt]{arr
} contains the coefficients.
101 For a
\ai[\tt]{periodic
}, it contains the values for the different
102 residue classes modulo the parameter indicated by
\ai[\tt]{pos
}.
103 For a polynomial,
\ai[\tt]{pos
} refers to the variable
105 The value of
\ai[\tt]{pos
} is
\verb+
1+ for the first parameter.
106 That is, if the value of
\ai[\tt]{pos
} is
\verb+
1+ and the first
107 parameter is $p$, and if the length of the array is $l$,
108 then in case it is a polynomial, the
109 \ai[\tt]{enode
} represents
111 \verb+arr
[0]+ +
\verb+arr
[1]+ p +
\verb+arr
[2]+ p^
2 +
\cdots +
112 \verb+arr
[l-
1]+ p^
{l-
1}
115 If it is a periodic, then it represents
118 \verb+arr
[0]+,
\verb+arr
[1]+,
\verb+arr
[2]+,
\ldots,
123 Note that the elements of a
\ai[\tt]{periodic
} may themselves
124 be other
\ai[\tt]{periodic
}s or even
\ai[\tt]{polynomial
}s.
125 In our library, we only allow the elements of a
\ai[\tt]{periodic
}
126 to be other
\ai[\tt]{periodic
}s or rational numbers.
127 The chambers and their corresponding quasi-polynomial are
128 stored in
\ai[\tt]{Enumeration
} structures.
130 typedef struct _enumeration
{
131 Polyhedron *ValidityDomain; /* constraints on the parameters */
132 evalue EP; /* dimension = combined space */
133 struct _enumeration *next; /* Ehrhart Polynomial,
134 corresponding to parameter
135 values inside the domain
136 ValidityDomain above */
139 For more information on these structures, we refer to
\shortciteN{Loechner1999
}.
142 Figure~
\ref{f:Loechner
} is a skillful reconstruction
143 of Figure~
2 from
\shortciteN{Loechner1999
}.
144 It shows the contents of the
\ai[\tt]{enode
} structures
145 representing the quasi-polynomial
147 [1,
2]_p p^
2 +
3 p +
\frac 5 2
154 \begin{tabular
}{|c|c|c|
}
156 \multicolumn{2}{|c|
}{type
} & polynomial \\
158 \multicolumn{2}{|c|
}{size
} &
3 \\
160 \multicolumn{2}{|c|
}{pos
} &
1 \\
162 \smash{\lower 6.25pt
\hbox{arr
[0]}} & d &
2 \\
166 \smash{\lower 6.25pt
\hbox{arr
[1]}} & d &
1 \\
170 \smash{\lower 6.25pt
\hbox{arr
[2]}} & d &
0 \\
177 +DR*!DR
\hbox{\strut\hskip 1.5\tabcolsep\phantom{\tt polynomial
}\hskip 1.5\tabcolsep}+C="a"
178 \POS(
60,-
15)*!UL
{\hbox{
180 \begin{tabular
}{|c|c|c|
}
182 \multicolumn{2}{|c|
}{type
} & periodic \\
184 \multicolumn{2}{|c|
}{size
} &
2 \\
186 \multicolumn{2}{|c|
}{pos
} &
1 \\
188 \smash{\lower 6.25pt
\hbox{arr
[0]}} & d &
1 \\
192 \smash{\lower 6.25pt
\hbox{arr
[1]}} & d &
1 \\
199 +UL+<
0.5\tabcolsep,
0pt>*!UL
\hbox{\strut}+CL="b"
201 \POS"box1"+UC*++!D
\hbox{\tt enode
}
202 \POS"box2"+UC*++!D
\hbox{\tt enode
}
204 \caption{The quasi-polynomial $
[1,
2]_p p^
2 +
3 p +
\frac 5 2$.
}
212 The
\ai[\tt]{barvinok
\_options} structure contains various
213 options that influence the behavior of the library.
216 struct barvinok_options
{
217 struct barvinok_stats *stats;
219 /* PolyLib options */
223 /* LLL reduction parameter delta=LLL_a/LLL_b */
227 /* barvinok options */
228 #define BV_SPECIALIZATION_BF
2
229 #define BV_SPECIALIZATION_DF
1
230 #define BV_SPECIALIZATION_RANDOM
0
231 #define BV_SPECIALIZATION_TODD
3
232 int incremental_specialization;
234 unsigned long max_index;
237 int count_sample_infinite;
239 int try_Delaunay_triangulation;
241 #define BV_APPROX_SIGN_NONE
0
242 #define BV_APPROX_SIGN_APPROX
1
243 #define BV_APPROX_SIGN_LOWER
2
244 #define BV_APPROX_SIGN_UPPER
3
245 int polynomial_approximation;
246 #define BV_APPROX_NONE
0
247 #define BV_APPROX_DROP
1
248 #define BV_APPROX_SCALE
2
249 #define BV_APPROX_VOLUME
3
250 #define BV_APPROX_BERNOULLI
4
251 int approximation_method;
252 #define BV_APPROX_SCALE_FAST (
1 <<
0)
253 #define BV_APPROX_SCALE_NARROW (
1 <<
1)
254 #define BV_APPROX_SCALE_NARROW2 (
1 <<
2)
255 #define BV_APPROX_SCALE_CHAMBER (
1 <<
3)
257 #define BV_VOL_LIFT
0
258 #define BV_VOL_VERTEX
1
259 #define BV_VOL_BARYCENTER
2
260 int volume_triangulate;
262 /* basis reduction options */
263 #define BV_GBR_NONE
0
264 #define BV_GBR_GLPK
1
268 /* bernstein options */
269 #define BV_BERNSTEIN_NONE
0
270 #define BV_BERNSTEIN_MAX
1
271 #define BV_BERNSTEIN_MIN -
1
272 int bernstein_optimize;
274 #define BV_BERNSTEIN_FACTORS
1
275 #define BV_BERNSTEIN_INTERVALS
2
276 int bernstein_recurse;
278 #define BV_LP_POLYLIB
0
285 struct barvinok_options *barvinok_options_new_with_defaults();
288 The function
\ai[\tt]{barvinok
\_options\_new\_with\_defaults}
289 can be used to create a
\ai[\tt]{barvinok
\_options} structure
293 \item \PolyLib/ options
297 \item \ai[\tt]{MaxRays
}
299 The value of
\ai[\tt]{MaxRays
} is passed to various
\PolyLib/
300 functions and defines the
301 maximum size of a table used in the
\ai{double description
} computation
302 in the
\PolyLib/ function
\ai[\tt]{Chernikova
}.
303 In earlier versions of
\PolyLib/,
304 this parameter had to be conservatively set
305 to a high number to ensure successful operation,
306 resulting in significant memory overhead.
307 Our change to allow this table to grow
308 dynamically is available in recent versions of
\PolyLib/.
309 In these versions, the value no longer indicates the maximal
310 table size, but rather the size of the initial allocation.
311 This value may be set to
\verb+
0+ or left as set
312 by
\ai[\tt]{barvinok
\_options\_new\_with\_defaults}.
316 \item \ai[\tt]{NTL
} options
320 \item \ai[\tt]{LLL
\_a}
321 \item \ai[\tt]{LLL
\_b}
323 The values used for the
\ai{reduction parameter
}
324 in the call to
\ai[\tt]{NTL
}'s implementation of
\indac{LLL
}.
328 \item \ai[\tt]{barvinok
} specific options
332 \item \ai[\tt]{incremental
\_specialization}
334 Selects the
\ai{specialization
} algorithm to be used.
335 If set to
{\tt 0} then a direct specialization is performed
336 using a random vector.
337 Value
{\tt 1} selects a depth first incremental specialization,
338 while value
{\tt 2} selects a breadth first incremental specialization.
339 The default is selected by the
\ai[\tt]{--enable-incremental
}
340 \ai[\tt]{configure
} option.
341 For more information we refer to~
\citeN[Section~
4.4.3]{Verdoolaege2005PhD
}.
347 \subsection{Data Structures for Quasi-polynomials
}
350 Internally, we do not represent our
\ai{quasi-polynomial
}s
351 as step-polynomials, but, similarly to
\shortciteN{Loechner1999
},
352 as polynomials with periodic numbers for coefficients.
353 However, we also allow our periodic numbers to be represented by
354 fractional parts of degree-$
1$ polynomials rather than
355 an explicit enumeration using the
\ai[\tt]{periodic
} type.
356 By default, the current version of
\barvinok/ uses
357 \ai[\tt]{periodic
}s, but this can be changed through
358 the
\ai[\tt]{--enable-fractional
} configure option.
359 In the latter case, the quasi-polynomial using fractional
360 parts can also be converted to an actual step-polynomial
361 using
\ai[\tt]{evalue
\_frac2floor}, but this is not fully
364 For reasons of compatibility,
%
365 \footnote{Also known as laziness.
}
366 we shoehorned our representations for piecewise quasi-polynomials
367 into the existing data structures.
368 To this effect, we introduced four new types,
369 \ai[\tt]{fractional
},
\ai[\tt]{relation
},
370 \ai[\tt]{partition
} and
\ai[\tt]{flooring
}.
372 typedef enum
{ polynomial, periodic, evector, fractional,
373 relation, partition, flooring
} enode_type;
375 The field
\ai[\tt]{pos
} is not used in most of these
376 additional types and is therefore set to
\verb+-
1+.
378 The types
\ai[\tt]{fractional
} and
\ai[\tt]{flooring
}
379 represent polynomial expressions in a fractional part or a floor respectively.
380 The generator is stored in
\verb+arr
[0]+, while the
381 coefficients are stored in the remaining array elements.
382 That is, an
\ai[\tt]{enode
} of type
\ai[\tt]{fractional
}
385 \verb+arr
[1]+ +
\verb+arr
[2]+ \
{\verb+arr
[0]+\
} +
386 \verb+arr
[3]+ \
{\verb+arr
[0]+\
}^
2 +
\cdots +
387 \verb+arr
[l-
1]+ \
{\verb+arr
[0]+\
}^
{l-
2}
390 An
\ai[\tt]{enode
} of type
\ai[\tt]{flooring
}
393 \verb+arr
[1]+ +
\verb+arr
[2]+
\lfloor\verb+arr
[0]+
\rfloor +
394 \verb+arr
[3]+
\lfloor\verb+arr
[0]+
\rfloor^
2 +
\cdots +
395 \verb+arr
[l-
1]+
\lfloor\verb+arr
[0]+
\rfloor^
{l-
2}
400 The internal representation of the quasi-polynomial
401 $$
\left(
1+
2 \left\
{\frac p
2\right\
}\right) p^
2 +
3 p +
\frac 5 2$$
402 is shown in Figure~
\ref{f:fractional
}.
408 \begin{tabular
}{|c|c|c|
}
410 \multicolumn{2}{|c|
}{type
} & polynomial \\
412 \multicolumn{2}{|c|
}{size
} &
3 \\
414 \multicolumn{2}{|c|
}{pos
} &
1 \\
416 \smash{\lower 6.25pt
\hbox{arr
[0]}} & d &
2 \\
420 \smash{\lower 6.25pt
\hbox{arr
[1]}} & d &
1 \\
424 \smash{\lower 6.25pt
\hbox{arr
[2]}} & d &
0 \\
431 +DR*!DR
\hbox{\strut\hskip 1.5\tabcolsep\phantom{\tt polynomial
}\hskip 1.5\tabcolsep}+C="a"
432 \POS(
60,
0)*!UL
{\hbox{
434 \begin{tabular
}{|c|c|c|
}
436 \multicolumn{2}{|c|
}{type
} & fractional \\
438 \multicolumn{2}{|c|
}{size
} &
3 \\
440 \multicolumn{2}{|c|
}{pos
} & -
1 \\
442 \smash{\lower 6.25pt
\hbox{arr
[0]}} & d &
0 \\
446 \smash{\lower 6.25pt
\hbox{arr
[1]}} & d &
1 \\
450 \smash{\lower 6.25pt
\hbox{arr
[2]}} & d &
1 \\
457 +UL+<
0.5\tabcolsep,
0pt>*!UL
\hbox{\strut}+CL="b"
459 \POS"box2"+UR*!UR
{\hbox{
473 }+CD*!U
{\strut}+C="c"
474 \POS(
60,-
50)*!UL
{\hbox{
476 \begin{tabular
}{|c|c|c|
}
478 \multicolumn{2}{|c|
}{type
} & polynomial \\
480 \multicolumn{2}{|c|
}{size
} &
2 \\
482 \multicolumn{2}{|c|
}{pos
} &
1 \\
484 \smash{\lower 6.25pt
\hbox{arr
[0]}} & d &
1 \\
488 \smash{\lower 6.25pt
\hbox{arr
[1]}} & d &
2 \\
495 +UR-<
0.8\tabcolsep,
0pt>*!UR
\hbox{\strut}+CR="d"
497 \POS"box1"+UC*++!D
\hbox{\tt enode
}
498 \POS"box2"+UC*++!D
\hbox{\tt enode
}
499 \POS"box3"+UC*++!D
\hbox{\tt enode
}
501 \caption{The quasi-polynomial
502 $
\left(
1+
2 \left\
{\frac p
2\right\
}\right) p^
2 +
3 p +
\frac 5 2$.
}
508 The
\ai[\tt]{relation
} type is used to represent
\ai{stride
}s.
509 In particular, if the value of
\ai[\tt]{size
} is
2, then
510 the value of a
\ai[\tt]{relation
} is (in pseudo-code):
512 (value(arr
[0]) ==
0) ? value(arr
[1]) :
0
514 If the size is
3, then the value is:
516 (value(arr
[0]) ==
0) ? value(arr
[1]) : value(arr
[2])
518 The type of
\verb+arr
[0]+ is typically
\ai[\tt]{fractional
}.
520 Finally, the
\ai[\tt]{partition
} type is used to
521 represent piecewise quasi-polynomials.
522 We prefer to encode this information inside
\ai[\tt]{evalue
}s
524 rather than using
\ai[\tt]{Enumeration
}s since we want
525 to perform the same kinds of operations on both quasi-polynomials
526 and piecewise quasi-polynomials.
527 An
\ai[\tt]{enode
} of type
\ai[\tt]{partition
} may not be nested
528 inside another
\ai[\tt]{enode
}.
529 The size of the array is twice the number of ``chambers''.
530 Pointers to chambers are stored in the even slots,
531 whereas pointer to the associated quasi-polynomials
532 are stored in the odd slots.
533 To be able to store pointers to chambers, the
534 definition of
\ai[\tt]{evalue
} was changed as follows.
536 typedef struct _evalue
{
537 Value d; /* denominator */
539 Value n; /* numerator (if denominator >
0) */
540 struct _enode *p; /* pointer (if denominator ==
0) */
541 Polyhedron *D; /* domain (if denominator == -
1) */
545 Note that we allow a ``chamber'' to be a union of polyhedra
546 as discussed in
\citeN[Section~
4.5.1]{Verdoolaege2005PhD
}.
547 Chambers with extra variables, i.e., those of
548 \citeN[Section~
4.6.5]{Verdoolaege2005PhD
},
549 are only partially supported.
550 The field
\ai[\tt]{pos
} is set to the actual dimension,
551 i.e., the number of parameters.
553 \subsection{Operations on Quasi-polynomials
}
556 In this section we discuss some of the more important
557 operations on
\ai[\tt]{evalue
}s provided by the
559 Some of these operations are extensions
560 of the functions from
\PolyLib/ with the same name.
563 void eadd(const evalue *e1,evalue *res);
564 void emul(const evalue *e1, evalue *res);
566 The functions
\ai[\tt]{eadd
} and
\ai[\tt]{emul
} takes
567 two (pointers to)
\ai[\tt]{evalue
}s
\verb+e1+ and
\verb+res+
568 and computes their sum and product respectively.
569 The result is stored in
\verb+res+, overwriting (and deallocating)
570 the original value of
\verb+res+.
571 It is an error if exactly one of
572 the arguments of
\ai[\tt]{eadd
} is of type
\ai[\tt]{partition
}
573 (unless the other argument is
\verb+
0+).
574 The addition and multiplication operations are described
575 in
\citeN[Section~
4.5.1]{Verdoolaege2005PhD
}
576 and~
\citeN[Section~
4.5.2]{Verdoolaege2005PhD
}
579 The function
\ai[\tt]{eadd
} is an extension of the function
580 \ai[\tt]{new
\_eadd} from
\shortciteN{Seghir2002
}.
581 Apart from supporting the additional types from Section~
\ref{a:data
},
582 the new version also additionally imposes an order on the nesting of
583 different
\ai[\tt]{enode
}s.
584 Without such an ordering,
\ai[\tt]{evalue
}s could be constructed
585 representing for example
587 (
0 y^
0 + (
0 x^
0 +
1 x^
1 ) y^
1 ) x^
0 + (
0 y^
0 -
1 y^
1) x^
1
590 which is just a funny way of saying $
0$.
593 void eor(evalue *e1, evalue *res);
595 The function
\ai[\tt]{eor
} implements the
\ai{union
}
596 operation from
\citeN[Section~
4.5.3]{Verdoolaege2005PhD
}. Both arguments
597 are assumed to correspond to indicator functions.
600 evalue *esum(evalue *E, int nvar);
601 evalue *evalue_sum(evalue *E, int nvar, unsigned MaxRays);
603 The function
\ai[\tt]{esum
} has been superseded by
604 \ai[\tt]{evalue
\_sum}.
605 The function
\ai[\tt]{evalue
\_sum} performs the summation
606 operation from
\citeN[Section~
4.5.4]{Verdoolaege2005PhD
}.
607 The piecewise step-polynomial represented by
\verb+E+ is summated
608 over its first
\verb+nvar+ variables.
609 Note that
\verb+E+ must be zero or of type
\ai[\tt]{partition
}.
610 The function returns the result in a newly allocated
612 Note also that
\verb+E+ needs to have been converted
613 from
\ai[\tt]{fractional
}s to
\ai[\tt]{flooring
}s using
614 the function
\ai[\tt]{evalue
\_frac2floor}.
616 void evalue_frac2floor(evalue *e);
618 This function also ensures that the arguments of the
619 \ai[\tt]{flooring
}s are positive in the relevant chambers.
620 It currently assumes that the argument of each
621 \ai[\tt]{fractional
} in the original
\ai[\tt]{evalue
}
622 has a minimum in the corresponding chamber.
625 double compute_evalue(const evalue *e, Value *list_args);
626 Value *compute_poly(Enumeration *en,Value *list_args);
627 evalue *evalue_eval(const evalue *e, Value *values);
629 The functions
\ai[\tt]{compute
\_evalue},
630 \ai[\tt]{compute
\_poly} and
631 \ai[\tt]{evalue
\_eval}
632 evaluate a (piecewise) quasi-polynomial
633 at a certain point. The argument
\verb+list_args+
634 points to an array of
\ai[\tt]{Value
}s that is assumed
636 The
\verb+double+ return value of
\ai[\tt]{compute
\_evalue}
637 is inherited from
\PolyLib/.
640 void print_evalue(FILE *DST, const evalue *e, char **pname);
642 The function
\ai[\tt]{print
\_evalue} dumps a human-readable
643 representation to the stream pointed to by
\verb+DST+.
644 The argument
\verb+pname+ points
645 to an array of character strings representing the parameter names.
646 The array is assumed to be long enough.
649 int eequal(const evalue *e1, const evalue *e2);
651 The function
\ai[\tt]{eequal
} return true (
\verb+
1+) if its
652 two arguments are structurally identical.
653 I.e., it does
{\em not\/
} check whether the two
654 (piecewise) quasi-polynomial represent the same function.
657 void reduce_evalue (evalue *e);
659 The function
\ai[\tt]{reduce
\_evalue} performs some
660 simplifications on
\ai[\tt]{evalue
}s.
661 Here, we only describe the simplifications that are directly
662 related to the internal representation.
663 Some other simplifications are explained in
664 \citeN[Section~
4.7.2]{Verdoolaege2005PhD
}.
665 If the highest order coefficients of a
\ai[\tt]{polynomial
},
666 \ai[\tt]{fractional
} or
\ai[\tt]{flooring
} are zero (possibly
667 after some other simplifications), then the size of the array
668 is reduced. If only the constant term remains, i.e.,
669 the size is reduced to $
1$ for
\ai[\tt]{polynomial
} or to $
2$
670 for the other types, then the whole node is replaced by
672 Additionally, if the argument of a
\ai[\tt]{fractional
}
673 has been reduced to a constant, then the whole node
674 is replaced by its partial evaluation.
675 A
\ai[\tt]{relation
} is similarly reduced if its second
676 branch or both its branches are zero.
677 Chambers with zero associated quasi-polynomials are
678 discarded from a
\ai[\tt]{partition
}.
680 \subsection{Generating Functions
}
682 The representation of
\rgf/s uses
683 some basic types from the
\ai[\tt]{NTL
} library
\shortcite{NTL
}
684 for representing arbitrary precision integers
686 as well as vectors (
\ai[\tt]{vec
\_ZZ}) and matrices (
\ai[\tt]{mat
\_ZZ})
688 We further introduces a type
\ai[\tt]{QQ
} for representing a rational
689 number and use vectors (
\ai[\tt]{vec
\_QQ}) of such numbers.
696 NTL_vector_decl(QQ,vec_QQ);
699 Each term in a
\rgf/ is represented by a
\ai[\tt]{short
\_rat}
704 /* rows: terms in numerator */
709 /* rows: factors in denominator */
714 The fields
\ai[\tt]{n
} and
\ai[\tt]{d
} represent the
715 numerator and the denominator respectively.
716 Note that in our implementation we combine terms
717 with the same denominator.
718 In the numerator, each element of
\ai[\tt]{coeff
} and each row of
\ai[\tt]{power
}
719 represents a single such term.
720 The vector
\ai[\tt]{coeff
} contains the rational coefficients
721 $
\alpha_i$ of each term.
722 The columns of
\ai[\tt]{power
} correspond to the powers
724 In the denominator, each row of
\ai[\tt]{power
}
725 corresponds to the power $
\vec b_
{ij
}$ of a
726 factor in the denominator.
730 shows the internal representation of
732 \frac{\frac 3 2 \, x_0^
2 x_1^
3 +
2 \, x_0^
5 x_1^
{-
7}}
733 { (
1 - x_0 x_1^
{-
3}) (
1 - x_1^
2)
}
739 \begin{minipage
}{0cm
}
743 \begin{tabular
}{|c|c|c|
}
758 }+UC*++!D
\hbox{\tt short
\_rat}
762 \caption{Representation of
764 \left(
\frac 3 2 \, x_0^
2 x_1^
3 +
2 \, x_0^
5 x_1^
{-
7}\right)
765 /
\left( (
1 - x_0 x_1^
{-
3}) (
1 - x_1^
2)
\right)
772 The whole
\rgf/ is represented by a
\ai[\tt]{gen
\_fun}
775 typedef std::set<short_rat *,
776 short_rat_lex_smaller_denominator > short_rat_list;
782 void add(const QQ& c, const vec_ZZ& num, const mat_ZZ& den);
783 void add(short_rat *r);
784 void add(const QQ& c, const gen_fun *gf);
785 void substitute(Matrix *CP);
786 gen_fun *Hadamard_product(const gen_fun *gf,
787 barvinok_options *options);
788 void print(std::ostream& os,
789 unsigned int nparam, char **param_name) const;
790 operator evalue *() const;
791 ZZ coefficient(Value* params, barvinok_options *options) const;
792 void coefficient(Value* params, Value* c) const;
794 gen_fun(Polyhedron *C = NULL);
796 gen_fun(const gen_fun *gf);
800 A new
\ai[\tt]{gen
\_fun} can be constructed either as empty
\rgf/ (possibly
801 with a given context
\verb+C+), as a copy of an existing
\rgf/
\verb+gf+, or as
802 constant
\rgf/ with value for the constant term specified by
\verb+c+.
804 The first
\ai[\tt]{gen
\_fun::add
} method adds a new term to the
\rgf/,
805 described by the coefficient
\verb+c+, the numerator
\verb+num+ and the
806 denominator
\verb+den+.
807 It makes all powers in the denominator lexico-positive,
808 orders them in lexicographical order and inserts the new
809 term in
\ai[\tt]{term
} according to the lexicographical
810 order of the combined powers in the denominator.
811 The second
\ai[\tt]{gen
\_fun::add
} method adds
\verb+c+ times
\verb+gf+
814 The method
\ai[\tt]{gen
\_fun::operator evalue *
} performs
815 the conversion from
\rgf/ to
\psp/ explained in
816 \citeN[Section~
4.5.5]{Verdoolaege2005PhD
}.
817 The
\ai[\tt]{Polyhedron
} \ai[\tt]{context
} is the superset
818 of all points where the enumerator is non-zero used during this conversion,
819 i.e., it is the set $Q$ from
\citeN[Equation~
4.31]{Verdoolaege2005PhD
}.
820 If
\ai[\tt]{context
} is
\verb+NULL+ the maximal
821 allowed context is assumed, i.e., the maximal
822 region with lexico-positive rays.
824 The method
\ai[\tt]{gen
\_fun::coefficient
} computes the coefficient
825 of the term with power given by
\verb+params+ and stores the result
827 This method performs essentially the same computations as
828 \ai[\tt]{gen
\_fun::operator evalue *
}, except that it adds extra
829 equality constraints based on the specified values for the power.
831 The method
\ai[\tt]{gen
\_fun::substitute
} performs the
832 \ai{monomial substitution
} specified by the homogeneous matrix
\verb+CP+
833 that maps a set of ``
\ai{compressed parameter
}s''
\shortcite{Meister2004PhD
}
834 to the original set of parameters.
835 That is, if we are given a
\rgf/ $G(
\vec z)$ that encodes the
836 explicit function $g(
\vec i')$, where $
\vec i'$ are the coordinates of
837 the transformed space, and
\verb+CP+ represents the map
838 $
\vec i = A
\vec i' +
\vec a$ back to the original space with coordinates $
\vec i$,
839 then this method transforms the
\rgf/ to $F(
\vec x)$ encoding the
840 same explicit function $f(
\vec i)$, i.e.,
841 $$f(
\vec i) = f(A
\vec i' +
\vec a) = g(
\vec i ').$$
842 This means that the coefficient of the term
843 $
\vec x^
{\vec i
} =
\vec x^
{A
\vec i' +
\vec a
}$ in $F(
\vec x)$ should be equal to the
844 coefficient of the term $
\vec z^
{\vec i'
}$ in $G(
\vec z)$.
848 \sum_i \epsilon_i \frac{\vec z^
{\vec v_i
}}{\prod_j (
1-
\vec z^
{\vec b_
{ij
}})
}
853 \sum_i \epsilon_i \frac{\vec x^
{A
\vec v_i +
\vec a
}}
854 {\prod_j (
1-
\vec x^
{A
\vec b_
{ij
}})
}
858 The method
\ai[\tt]{gen
\_fun::Hadamard
\_product} computes the
859 \ai{Hadamard product
} of the current
\rgf/ with the
\rgf/
\verb+gf+,
860 as explained in
\citeN[Section~
4.5.2]{Verdoolaege2005PhD
}.
862 \subsection{Counting Functions
}
863 \label{a:counting:functions
}
865 Our library provides essentially three different counting functions:
866 one for non-parametric polytopes, one for parametric polytopes
867 and one for parametric sets with existential variables.
868 The old versions of these functions have a ``
\ai[\tt]{MaxRays
}''
869 argument, while the new versions have a more general
870 \ai[\tt]{barvinok
\_options} argument.
871 For more information on
\ai[\tt]{barvinok
\_options}, see Section~
\ref{a:options
}.
874 void barvinok_count(Polyhedron *P, Value* result,
876 void barvinok_count_with_options(Polyhedron *P, Value* result,
877 struct barvinok_options *options);
879 The function
\ai[\tt]{barvinok
\_count} or
880 \ai[\tt]{barvinok
\_count\_with\_options} enumerates the non-parametric
881 polytope
\verb+P+ and returns the result in the
\ai[\tt]{Value
}
882 pointed to by
\verb+result+, which needs to have been allocated
884 If
\verb+P+ is a union, then only the first set in the union will
885 be taken into account.
886 For the meaning of the argument
\verb+NbMaxCons+, see
887 the discussion on
\ai[\tt]{MaxRays
} in Section~
\ref{a:options
}.
889 The function
\ai[\tt]{barvinok
\_enumerate} for enumerating
890 parametric polytopes was meant to be
891 a drop-in replacement of
\PolyLib/'s
\ai[\tt]{Polyhedron
\_Enumerate}
893 Unfortunately, the latter has been changed to
894 accept an extra argument in recent versions of
\PolyLib/ as shown below.
896 Enumeration* barvinok_enumerate(Polyhedron *P, Polyhedron* C,
898 extern Enumeration *Polyhedron_Enumerate(Polyhedron *P,
899 Polyhedron *C, unsigned MAXRAYS, char **pname);
901 The argument
\verb+MaxRays+ has the same meaning as the argument
902 \verb+NbMaxCons+ above.
903 The argument
\verb+P+ refers to the $(d+n)$-dimensional
904 polyhedron defining the parametric polytope.
905 The argument
\verb+C+ is an $n$-dimensional polyhedron containing
906 extra constraints on the parameter space.
907 Its primary use is to indicate how many of the dimensions
908 in
\verb+P+ refer to parameters as any constraint in
\verb+C+
909 could equally well have been added to
\verb+P+ itself.
910 Note that the dimensions referring to the parameters should
912 If either
\verb+P+ or
\verb+C+ is a union,
913 then only the first set in the union will be taken into account.
914 The result is a newly allocated
\ai[\tt]{Enumeration
}.
915 As an alternative we also provide a function
916 (
\ai[\tt]{barvinok
\_enumerate\_ev} or
917 \ai[\tt]{barvinok
\_enumerate\_with\_options}) that returns
920 evalue* barvinok_enumerate_ev(Polyhedron *P, Polyhedron* C,
922 evalue* barvinok_enumerate_with_options(Polyhedron *P,
923 Polyhedron* C, struct barvinok_options *options);
926 For enumerating parametric sets with existentially quantified variables,
927 we provide two functions:
928 \ai[\tt]{barvinok
\_enumerate\_e}
930 \ai[\tt]{barvinok
\_enumerate\_pip}.
932 evalue* barvinok_enumerate_e(Polyhedron *P,
933 unsigned exist, unsigned nparam, unsigned MaxRays);
934 evalue* barvinok_enumerate_e_with_options(Polyhedron *P,
935 unsigned exist, unsigned nparam,
936 struct barvinok_options *options);
937 evalue *barvinok_enumerate_pip(Polyhedron *P,
938 unsigned exist, unsigned nparam, unsigned MaxRays);
939 evalue* barvinok_enumerate_pip_with_options(Polyhedron *P,
940 unsigned exist, unsigned nparam,
941 struct barvinok_options *options);
942 evalue *barvinok_enumerate_scarf(Polyhedron *P,
943 unsigned exist, unsigned nparam,
944 struct barvinok_options *options);
946 The first function tries the simplification rules from
947 \citeN[Section~
4.6.2]{Verdoolaege2005PhD
} before resorting to the method
948 based on
\indac{PIP
} from
\citeN[Section~
4.6.3]{Verdoolaege2005PhD
}.
949 The second function immediately applies the technique from
950 \citeN[Section~
4.6.3]{Verdoolaege2005PhD
}.
951 The argument
\verb+exist+ refers to the number of existential variables,
953 the argument
\verb+nparam+ refers to the number of parameters.
954 The order of the dimensions in
\verb+P+ is:
955 counted variables first, then existential variables and finally
957 The function
\ai[\tt]{barvinok
\_enumerate\_scarf} performs the same
958 computation as the function
\ai[\tt]{barvinok
\_enumerate\_scarf\_series}
959 below, but produces an explicit representation instead of a generating function.
962 gen_fun * barvinok_series(Polyhedron *P, Polyhedron* C,
964 gen_fun * barvinok_series_with_options(Polyhedron *P,
965 Polyhedron* C, barvinok_options *options);
966 gen_fun *barvinok_enumerate_scarf_series(Polyhedron *P,
967 unsigned exist, unsigned nparam,
968 barvinok_options *options);
971 \ai[\tt]{barvinok
\_series} or
972 \ai[\tt]{barvinok
\_series\_with\_options} enumerates parametric polytopes
973 in the form of a
\rgf/.
974 The polyhedron
\verb+P+ is assumed to have only
975 revlex-positive rays.
977 The function
\ai[\tt]{barvinok
\_enumerate\_scarf\_series} computes a
978 generating function for the number of point in the parametric set
979 defined by
\verb+P+ with
\verb+exist+ existentially quantified
980 variables, which is assumed to be
2.
981 This function implements the technique of
982 \shortciteN{Scarf2006Neighborhood
} using the
\ai{neighborhood complex
}
983 description of
\shortciteN{Scarf1981indivisibilities:II
}.
984 It is currently restricted to problems with
3 or
4 constraints involving
985 the existentially quantified variables.
987 \subsection{Auxiliary Functions
}
989 In this section we briefly mention some auxiliary functions
990 available in the
\barvinok/ library.
993 void Polyhedron_Polarize(Polyhedron *P);
995 The function
\ai[\tt]{Polyhedron
\_Polarize}
996 polarizes its argument and is explained
997 in
\citeN[Section~
4.4.2]{Verdoolaege2005PhD
}.
1000 int unimodular_complete(Matrix *M, int row);
1002 The function
\ai[\tt]{unimodular
\_complete} extends
1003 the first
\verb+row+ rows of
1004 \verb+M+ with an integral basis of the orthogonal complement
1005 as explained in Section~
\ref{s:completion
}.
1007 if the resulting matrix is unimodular
\index{unimodular matrix
}.
1010 int DomainIncludes(Polyhedron *D1, Polyhedron *D2);
1012 The function
\ai[\tt]{DomainIncludes
} extends
1013 the function
\ai[\tt]{PolyhedronIncludes
}
1014 provided by
\PolyLib/
1015 to unions of polyhedra.
1016 It checks whether every polyhedron in the union
{\tt D2
}
1017 is included in some polyhedron of
{\tt D1
}.
1020 Polyhedron *DomainConstraintSimplify(Polyhedron *P,
1023 The value returned by
1024 \ai[\tt]{DomainConstraintSimplify
} is a pointer to
1025 a newly allocated
\ai[\tt]{Polyhedron
} that contains the
1026 same integer points as its first argument but possibly
1027 has simpler constraints.
1028 Each constraint $ g
\sp a x
\ge c $
1029 is replaced by $
\sp a x
\ge \ceil{ \frac c g
} $,
1030 where $g$ is the
\ac{gcd
} of the coefficients in the original
1032 The
\ai[\tt]{Polyhedron
} pointed to by
\verb+P+ is destroyed.
1035 Polyhedron* Polyhedron_Project(Polyhedron *P, int dim);
1037 The function
\ai[\tt]{Polyhedron
\_Project} projects
1038 \verb+P+ onto its last
\verb+dim+ dimensions.
1041 Matrix *left_inverse(Matrix *M, Matrix **Eq);
1043 The
\ai[\tt]{left
\_inverse} function computes the left inverse
1044 of
\verb+M+ as explained in Section~
\ref{s:inverse
}.
1046 \sindex{reduced
}{basis
}
1047 \sindex{generalized
}{reduced basis
}
1049 Matrix *Polyhedron_Reduced_Basis(Polyhedron *P,
1050 struct barvinok_options *options);
1052 \ai[\tt]{Polyhedron
\_Reduced\_Basis} computes
1053 a
\ai{generalized reduced basis
} of
{\tt P
}, which
1054 is assumed to be a polytope, using the algorithm
1055 of~
\shortciteN{Cook1993implementation
}.
1056 The basis vectors are stored in the rows of the matrix returned.
1057 This function currently uses
\ai[\tt]{GLPK
}~
\shortcite{GLPK
}
1058 to perform the linear optimizations and so is only available
1059 if you have
\ai[\tt]{GLPK
}.
1062 Vector *Polyhedron_Sample(Polyhedron *P,
1063 struct barvinok_options *options);
1065 \ai[\tt]{Polyhedron
\_Sample} returns an
\ai{integer point
} of
{\tt P
}
1066 or
{\tt NULL
} if
{\tt P
} contains no integer points.
1067 The integer point is found using the algorithm
1068 of~
\shortciteN{Cook1993implementation
} and uses
1069 \ai[\tt]{Polyhedron
\_Reduced\_Basis} to compute the reduced bases
1070 and therefore also requires
\ai[\tt]{GLPK
}.
1072 \subsection{\protect\ai[\tt]{bernstein
} Data Structures and Functions
}
1074 The
\bernstein/ library used
\ai[\tt]{GiNaC
} data structures to
1075 represent the data it manipulates.
1076 In particular, a polynomial is stored in a
\ai[\tt]{GiNaC::ex
},
1077 a list of variable or parameter names is stored in a
\ai[\tt]{GiNaC::exvector
},
1078 while the parametric vertices or generators are stored in a
\ai[\tt]{GiNaC::matrix
},
1079 where the rows refer to the generators and the columns to the coordinates
1083 namespace bernstein
{
1084 GiNaC::exvector constructParameterVector(
1085 const char * const *param_names, unsigned nbParams);
1086 GiNaC::exvector constructVariableVector(unsigned nbVariables,
1087 const char *prefix);
1090 The functions
\ai[\tt]{constructParameterVector
}
1091 and
\ai[\tt]{constructVariableVector
} construct a list of variable
1092 names either from a list of
{\tt char *
}s or
1093 by suffixing
{\tt prefix
} with a number starting from
0.
1094 Such lists are needed for the functions
1095 \ai[\tt]{domainVertices
},
\ai[\tt]{bernsteinExpansion
}
1096 and
\ai[\tt]{evalue
\_bernstein\_coefficients}.
1099 namespace bernstein
{
1100 GiNaC::matrix domainVertices(Param_Polyhedron *PP, Param_Domain *Q,
1101 const GiNaC::exvector& params);
1104 The function
\ai[\tt]{domainVertices
} constructs a matrix representing the
1105 generators (in this case vertices) of the
\ai[\tt]{Param
\_Polyhedron} {\tt PP
}
1106 for the
\ai[\tt]{Param
\_Domain} {\tt Q
}, to be used
1107 in a call to
\ai[\tt]{bernsteinExpansion
}.
1108 The elements of
{\tt params
} are used in the resulting matrix
1109 to refer to the parameters.
1112 namespace bernstein
{
1113 GiNaC::lst bernsteinExpansion(const GiNaC::matrix& vert,
1114 const GiNaC::ex& poly,
1115 const GiNaC::exvector& vars,
1116 const GiNaC::exvector& params);
1119 The function
\ai[\tt]{bernsteinExpansion
} computes the
1120 \ai{Bernstein coefficient
}s of the polynomial
\verb+poly+
1121 over the
\ai{parametric polytope
} that is the
\ai{convex hull
}
1122 of the rows in
\verb+vert+. The vectors
\verb+vars+
1123 and
\verb+params+ identify the variables (i.e., the coordinates
1124 of the space in which the parametric polytope lives) and
1125 the parameters, respectively.
1128 namespace bernstein
{
1130 typedef std::pair< Polyhedron *, GiNaC::lst > guarded_lst;
1132 struct piecewise_lst
{
1133 const GiNaC::exvector vars;
1134 std::vector<guarded_lst> list;
1135 /*
0: just collect terms
1136 *
1: remove obviously smaller terms (maximize)
1137 * -
1: remove obviously bigger terms (minimize)
1141 piecewise_lst(const GiNaC::exvector& vars);
1142 piecewise_lst& combine(const piecewise_lst& other);
1144 void simplify_domains(Polyhedron *ctx, unsigned MaxRays);
1145 GiNaC::numeric evaluate(const GiNaC::exvector& values);
1146 void add(const GiNaC::ex& poly);
1151 A
\ai[\tt]{piecewise
\_list} structure represents a list of (disjoint)
1152 polyhedral domains, each with an associated
\ai[\tt]{GiNaC::lst
}
1154 The
\ai[\tt]{vars
} member contains the variable names of the
1155 dimensions of the polyhedral domains.
1157 \ai[\tt]{piecewise
\_lst::combine
} computes the
\ai{common refinement
}
1158 of the polyhedral domains in
\verb+this+ and
\verb+other+ and associates
1159 to each of the resulting subdomains the union of the sets of polynomials
1160 associated to the domains from
\verb+this+ and
\verb+other+ that contain
1162 If the
\verb+sign+s of the
\ai[\tt]{piecewise
\_list}s are not zero,
1163 then the (obviously) redundant elements of these sets are removed
1165 The result is stored in
\verb+this+.
1167 \ai[\tt]{piecewise
\_lst::maximize
} removes polynomials from domains that evaluate
1168 to a value that is smaller than or equal to the value of some
1169 other polynomial associated to the same domain for each point in the domain.
1171 \ai[\tt]{piecewise
\_lst::evaluate
} ``evaluates'' the
\ai[\tt]{piecewise
\_list}
1172 by looking for the domain (if any) that contains the point given by
1173 \verb+values+ and computing the maximal value attained by any of the
1174 associated polynomials evaluated at that point.
1176 \ai[\tt]{piecewise
\_lst::add
} adds the polynomial
\verb+poly+
1177 to each of the polynomial associated to each of the domains.
1179 \ai[\tt]{piecewise
\_lst::simplify
\_domains} ``simplifies'' the domains
1180 by removing the constraints that are implied by the constraints
1181 in
\verb+ctx+, basically by calling
\PolyLib/'s
1182 \ai[\tt]{DomainSimplify
}. Note that you should only do this
1183 at the end of your computation. In particular, you do not
1184 want to call this method before calling
1185 \ai[\tt]{piecewise
\_lst::maximize
}, since this method will then
1186 have less information on the domains to exploit.
1190 namespace barvinok
{
1191 bernstein::piecewise_lst *evalue_bernstein_coefficients(
1192 bernstein::piecewise_lst *pl_all, evalue *e,
1193 Polyhedron *ctx, const GiNaC::exvector& params);
1194 bernstein::piecewise_lst *evalue_bernstein_coefficients(
1195 bernstein::piecewise_lst *pl_all, evalue *e,
1196 Polyhedron *ctx, const GiNaC::exvector& params,
1197 barvinok_options *options);
1200 The
\ai[\tt]{evalue
\_bernstein\_coefficients} function will compute the
1201 \ai{Bernstein coefficient
}s of the piecewise parametric polynomial stored in the
1202 \ai[\tt]{evalue
} \verb+e+.
1203 The
\verb+params+ vector specifies the names to be used for the parameters,
1204 while the context
\ai[\tt]{Polyhedron
} \verb+ctx+ specifies extra constraints
1206 The dimension of
\verb+ctx+ needs to be the same as the length of
\verb+params+.
1207 The
\ai[\tt]{evalue
} \verb+e+ is assumed to be of type
\ai[\tt]{partition
}
1208 and each of the domains in this
\ai[\tt]{partition
} is interpreted
1209 as a parametric polytope in the given parameters. The procedure
1210 will compute the
\ai{Bernstein coefficient
}s of the associated polynomial
1211 over each such parametric polytope.
1212 The resulting
\ai[\tt]{bernstein::piecewise
\_lst} collects the
1213 Bernstein coefficients over all parametric polytopes in
\verb+e+.
1214 If
\verb+pl_all+ is not
\verb+NULL+ then this list will be combined
1215 with the list computed by calling
\ai[\tt]{piecewise
\_lst::combine
}.
1216 If
\ai[\tt]{bernstein
\_optimize} is set to
\ai[\tt]{BV
\_BERNSTEIN\_MAX}
1217 in
\verb+options+, then this combination will remove obviously
1218 redundant Bernstein coefficients with respect to upper bound computation
1219 and similarly for
\ai[\tt]{BV
\_BERNSTEIN\_MIN}.
1220 The default (
\ai[\tt]{BV
\_BERNSTEIN\_NONE}) is to only remove duplicate
1221 Bernstein coefficients.