This paper describes an algorithm that can be used to effectively solve polynomial constraints over finite domains. Such constraints are expressed in terms of inequalities of polynomials with integer coefficients whose variables are assumed to be defined over proper finite domains. The proposed algorithm first reduces each constraint to a canonical form, i.e., a specific form of inequality, then it uses the modified Bernstein form of resulting polynomials to incrementally restrict the domains of variables. No approximation is involved in the solving process because the coefficients of the modified Bernstein form of the considered type of polynomials are always integer numbers.
In this paper we study constraints written in terms of polynomials whose coefficients are integer numbers and whose variables are defined over proper finite subsets of integer numbers. Such constraints, that we call polynomial constraints over finite domains, are ubiquitous in many areas of system analysis and verification (as discussed, e.g., in [1]) and they deserve specific treatment because common finite-domain techniques cannot adequately process them, even in simple cases. The language of CLP(FD) is suitable to express polynomial constraints over finite domains, but the fact that FD solvers are incomplete sometimes does not allow to treat such constraints adequately. For example, if constraints are expressed as non-linear polynomials, the FD solver that ships with SWI-Prolog (version 6.6.6) may fail in completely reducing the domains of variables [6]. Consider the constraint
X in -10..10, X^2 #>= 9 which involves only one variable over a finite interval, namely I = [−10..10]. It is easy to observe that such a constraint is satisfied only for values in the interval IR = [−10.. − 3] ∪ [3..10]. However, the FD solver exhibits its well-known incompleteness and it rewrites the constraint as which introduces an unneeded variable and counts {−2, −1, 1, 2} as admissible values. Concerning non-linear constraints that involve more than one variable, consider where both variables are over the finite domain I. Again, such a constraint is satisfied only for values in IR for both variables, but the FD solver introduces an unneeded variable and it also allows inadmissible values in the domains of variables. As a matter of fact, the FD solver rewrites the constraint as X in -10..-1\/1..10, Y in -10..-1\/1..10, Z in 21..100, X*Y #= Z.
In this paper we describe an algorithm to solve polynomial constraints over finite domains that overcomes the well-known incompleteness of FD solvers. The proposed algorithm ensures that the domains of variables after the solution process contain only admissible values by a proper use of the specific form of allowed constraints, namely polynomial inequalities. Polynomials with real coefficients whose variables admit values in known intervals can always be rewritten in so called Bernstein Form (BF ) (see, e.g., [4] or [5]) to effectively compute a lower bound and an upper bound of their range. Such bounds can be used to effectively verify polynomial inequalities, and they have been used extensively for optimization problems. Unfortunately, even if a polynomial has only integer coefficients, no guarantees can be provided that its equivalent in BF has integer coefficients. Therefore, the BF has been traditionally used only with polynomials with real or rational variables. In this paper we show how to use a variant of the BF, that we call Modified Bernstein Form (MBF ), to process polynomial constraints over finite domains with no involvement of rational or real numbers. Note that the MBF is also known in the literature as scaled BF [2], but we adopted the modified adjective since the introduced multiplicative factor is not the same for all coefficients of the BF, and the resulting patch is not geometrically scaled. As a matter of fact, if a polynomial has integer coefficients, its equivalent in MBF always has integer coefficients. This property allows effective reasoning on polynomial inequalities with integer coefficients. In details, even if the MBF does not provide accurate bounds for the range of polynomials, it provides information on the sign of suitable bounds, which is sufficient to reason on polynomial constraints over finite domains. The algorithm that we propose uses only such signs, and the actual values of bounds are not used. The current implementation of the proposed algorithm was used to validate all examples presented in this paper and it is available for download (cmt.dmi.unipr.it/software/clppolyfd.zip).
The paper is organized as follows. Section 2 presents the proposed approach and the underlying constraint language. Section 3 introduces the MBF and it proves notable properties that allow using the coefficients of the MBF to compute the signs of suitable lower and upper bounds of polynomials, as needed by the proposed algorithm. Section 4 illustrates in details the behaviour of the proposed algorithm by means of specific examples. Finally, Section 5 concludes the paper and presents planned developments.
The interesting characteristics of the constraint solving approach that we propose derive from the specific form of constraints that we address. Polynomial constraints over finite domains are expressed in terms of inequalities of polynomials with integer coefficients whose variables are defined over finite subsets of integer numbers. The relative constraint language is described as follows. 2.1
The syntax of the language of polynomial constraints that we consider is defined, as usual, by its signature Σ, a triple Σ =< V , F , Π > where V is a denumerable set of variable symbols, F is the set of constant and function symbols, and Π is the finite set of constraint predicate symbols allowed in the language. The set F of constant and function symbols is F = O ∪ Z, where O = {+, ∗} is the set of function symbols representing binary operations over integer numbers, and Z = {0, 1, −1, 2, −2, . . .} is the denumerable set of constants representing integer numbers. The set Π of constraint predicate symbols is Π = {=, 6=, <, ≤, >, ≥}, and it contains only binary predicate symbols.
A primitive constraint is any atomic predicate built using the symbols from signature Σ. A (non-primitive) constraint is a conjunction of primitive constraints. The semantics of the language defined over signature Σ is trivial and it allows expressing systems of polynomial equalities, inequalities and disequalities with integer coefficients.
Example 1 The following are primitive constraints expressed using signature Σ with V = {x, y, z}:
x ∗ x ≥ 9 3 ∗ x ∗ x + 9 ≤ 5 ∗ y ∗ y + z
x 6= y + z.
Note that other common function symbols and predicate symbols are supported in terms of syntactic abbreviations, as follows
v in a..b → (v + (−1) ∗ a) ∗ (v + (−1) ∗ b) ≤ 0 v nin a..b → (v + (−1) ∗ a) ∗ (v + (−1) ∗ b) > 0
−t → (−1) ∗ t
t1 − t2 → t1 + (−1) ∗ t2
t0 → 1
tn → t ∗ t ∗ · · · ∗ t
| n{>z0 }
(
Note that the extended language introduced by the aforementioned syntactic abbreviations does not increase the set of expressible constraints, and that we allow operator precedence and parenthesized expressions in this language.
The following example is meant to clarify why v in a..b and v nin a..b (read v not in a..b) are rewritten using the mentioned syntactic abbreviations. Example 2 Syntactic rules rewrite the constraint x in 1..3 as
(x + (−1) ∗ 1) ∗ (x + (−1) ∗ 3) ≤ 0.
As a matter of fact, p(x) = (x − 1)(3 − x) represents a downward parabola whose
vertex is x = 2 (which corresponds to p(
Analogously, syntactic rules rewrite the constraint x nin 1..3 as
(x + (−1) ∗ 1) ∗ (x + (−1) ∗ 3) > 0.
As a matter of fact, p(x) = (x − 1)(x − 3) represents an upward parabola whose
vertex is x = 2 (which corresponds to p(
If a total order on variable symbols is fixed, any primitive constraint can be rewritten to a canonical form. The chosen canonical form expresses any primitive constraint as t ≥ 0, where t is a term in canonical form. We say that a term t is in canonical form if it is a summation term (i.e., t = t1 + . . . + th) involving only terms with the following structure ti = k∗x1 ∗ x1 ∗ · · · ∗ x1 ∗ x2 ∗ x2 ∗ · · · ∗ x2 ∗ · · ·∗xm ∗ xm ∗ · · · ∗ xm | {nz1 } | {nz2 } | n{mz } i ∈ [1..h] where k is a constant symbol and the {xi}im=1 are distinct variable symbols arranged using the fixed total order. The fixed total order induces a lexicographic h order of product terms that is used to arrange the {ti}i=1 in t. We can define a function deg(v, t) on product terms by returning nv if variable symbol v is repeated nv times in t. Function deg(v, t) allows writing product terms in a compact form as
deg({xz1,ti) ti = k ∗ x1 ∗ x1 ∗ · · · ∗ x1 ∗ x2 ∗ x2 ∗ · · · ∗ x2 ∗ · · · ∗ xm ∗ xm ∗ · · · ∗ xm | } | } | } deg({xz2,ti) deg({xzm,ti) = k ∗ x1ˆdeg(x1, ti) ∗ x2ˆdeg(x2, ti) ∗ · · · ∗ xmˆdeg(xm, ti) i ∈ [1..h]
Besides ordinary algebraic manipulations and the total order of variable symbols, the conversion to the canonical form uses the following ideas. First, observe that any constraint of the form f (x) ⊙ g(x), where f (x) and g(x) are polynomials with integer coefficients, x is the vector of variables of polynomials, and ⊙ ∈ Π, can be expressed in the form p(x) ⊙ 0, where p(x) = f (x) − g(x) is still a polynomial with integer coefficients. Second, in order to convert constraints in canonical form we use the following lemma, which holds for every polynomial p(x) with integer coefficients.
Lemma 1 The following co-implications hold for every polynomial p(x) with integer coefficients and variables x p(x) ≤ 0 ⇐⇒ −p(x) ≥ 0 p(x) > 0 ⇐⇒ p(x) ≥ 1
⇐⇒ p(x) − 1 ≥ 0 p(x) < 0 ⇐⇒ p(x) ≤ −1 ⇐⇒ −p(x) − 1 ≥ 0 p(x) = 0 ⇐⇒ p2(x) ≤ 0 ⇐⇒ −p2(x) ≥ 0 p(x) 6= 0 ⇐⇒ p2(x) > 0 ⇐⇒ p2(x) − 1 ≥ 0 (7) (8) (9) (10) (11)
The following illustrative examples are provided to exemplify how generic constraints can be written in canonical form.
Example 3 Consider the following primitive constraints in the form (x1+x2)⊙0 where ⊙ ∈ Π. 1. The constraint x1 + x2 ≤ 0 can be written as −x1 − x2 ≥ 0 by co-implication (7) of lemma 1, which is in canonical form. 2. The constraint x1 +x2 > 0 can be written as x1 +x2 −1 ≥ 0 by co-implication (8) of lemma 1, which is in canonical form. 3. The constraint x1 + x2 < 0 can be written as x1 + x2 ≤ −1 by co-implication (9) of lemma 1, which corresponds to −x1 − x2 − 1 ≥ 0 in canonical form. 4. The constraint x1 + x2 = 0 can be written as (x1 + x2)2 ≤ 0 by co-implication (10) of lemma 1, which is −x21 − x22 − 2x1x2 ≥ 0 in canonical form. 5. The constraint x1 + x2 6= 0 can be written as (x1 + x2)2 > 0 by co-implication (11) of lemma 1, and it can be written as x21 +2x1x2 +x22 −1 ≥ 0 in canonical form. 2.3
The canonical form of primitive constraints has the interesting property of reducing constraint solving to the study of the sign of polynomials with integer coefficients whose variables are defined over finite subsets of integer numbers. In the rest of this paper we assume that constraints are always written in canonical form since the reformulation of constraints to their respective canonical forms can be considered a simple preprocessing step.
As discussed in Section 1, we assume that the user provides initial bounds for all variables, i.e., we can assume that before the solving process begins xi ∈ [xi..xi] xi ∈ Z xi ∈ Z i ∈ [1..m] (12) where {xi}im=1 are all the variables involved in the constraints, {xi}im=1 are their lower bounds, and {xi}im=1 are their upper bounds. The lower bounds and the upper bounds identify a box in Zm that can be used as an initial approximation of the domain of variables. The constraint solving algorithm refines the initial box and possibly breaks it into disjoint boxes until boxes are sufficiently refined and the consistency of constraints can be certified in each of them. In details, the algorithm considers all primitive constraints and it verifies if each constraint is consistent in the current box. If no consistency can be certified, a variable is selected and the current box is split into two disjoint boxes by splitting the current domain of the selected variable. The obtained disjoint boxes are processed recursively. The algorithms always terminates because the consistency of constraints can always be verified in a box formed by intervals that contain only one element. The performances of the algorithm depend significantly on the details of the subdivision, i.e., on the techniques adopted to select a variable and to split a box into two disjoint boxes, and a large number of alternatives have been studied in the literature (see, e.g., [3]). In the rest of this paper, we assume that variables are selected using a fixed lexicographic order and that once a variable xi is selected, its domain [xi..xi] is split at si = (xi + xi) 2 (13) into two disjoint intervals [xi..si] and [(si + 1)..xi].
For each primitive constraint p(x) ≥ 0, the algorithm computes suitable lower bound l and upper bound u, with l ≤ u, for the polynomial p(x) in the current box, and it uses such values to verify if p(x) ≥ 0 in the current box. The following procedure summarizes the core of the constraint solving algorithm: 1. If u < 0 the constraint is not consistent in the current box; 2. If l ≥ 0 the constraint is consistent in the current box; 3. Otherwise, the current box should be split and analysed recursively. The following example shows how the proposed algorithm works for a simple constraint under the assumption that suitable lower bounds and upper bounds can be computed.
Example 4 Let us consider the constraint x − 3 ≥ 0 involving only variable x,
which is assumed to be defined in I = [0..5]. Using the notation introduced in
the algorithm described above, suitable values of the lower and upper bounds are
l = −3 and u = 2, respectively, so that none of the first two conclusive cases of
the algorithm applies. The interval I is then split into two subintervals, namely
IL = [0..2] and IR = [3..5]. First, consider IL with suitable values of the lower
and upper bounds as l = −3 and u = −1, respectively, so that case (
It is worth noting that we are only interested in the signs of l and u, and not in their actual value. Any algorithm that allows the effective computation of the signs of suitable l and u can be used to support the proposed algorithm. In the following section we present a technique to effectively compute such signs using the modified Bernstein form of polynomials. 3
The Modified Bernstein Form of Polynomials In this section we introduce the Modified Bernstein Form (MBF ) of polynomials and we explain how to use it to solve polynomial constraints over finite domains. As an illustrative special case, we start by considering univariate polynomials. Problems involving only polynomials in one variable are not relevant for typical applications, however, they allow simplifying the notation and, therefore, they allow the proposed approach to be explained in a simpler way. 3.1
We consider a single constraint expressed in canonical form as (14) (15) (16) is a generic polynomial of degree n, x is the (unique) variable, {ai}in=0 are the coefficients relative to the canonical basis B given by the following set of monomials
B = {1, x, x2, . . . , xn} and we assume that x is defined in a given interval I = [x..x]. In order to write the explicit expression of the polynomial p(x) in the Bernstein basis, we first need to consider a new variable t defined in [0, 1], which is related to x according to the following affine transformation
x = x + (x − x)t.
By substituting (17) in (15), one obtains
Xn ai(x + (x − x)t)i = Xn Xi ai ki xi−k(x − x)ktk = Xn citi i=0 i=0 k=0 i=0 where the coefficients {ci}in=0 are defined as
n ci = X k=i k i akxk−i(x − x)i i ∈ [0..n].
Observe that {ci}in=0 are related to the coefficients {ai}in=0 relative to the canonical basis (16) and also to the extremes x and x of the interval I where p(x) is defined. Moreover, since the coefficients {ai}in=0 are integer numbers, it can be concluded that also {ci}in=0 are integer numbers. The coefficients {ci}in=0 defined in (19) are necessary to derive the MBF of p(x).
The Bernstein basis is given by the set of n + 1 polynomials {Bin(x)}in=0, the i-th of which is given by
The coefficients {bi}in=0 are useful to find a lower bound and an upper bound of the polynomial p(x). As a matter of fact the following inequalities hold i∈m[0i.n.n]{bi} ≤ p(x) ≤ i∈m[0a..xn]{bi} x ∈ I.
Observe that, even if the coefficients {ci}in=0 are integer numbers, the coefficients {bi}in=0 are not guaranteed to be integer because of the divisions by binomial coefficients. For this reason, the treatment of the {bi}in=0 requires representing rational numbers, which often involves approximations. However, it is possible to consider a modified Bernstein basis whose i-th element is given by B˜in(x) = (x − x)i(x − x)n−i (x − x)n i ∈ [0..n].
(24) n p(x) = X biBin(x)
i=0 i i bi = X nj cj j=0 j i ∈ [0..n]. where the coefficients {bi}in=0 are related to {ci}in=0, and, hence, to {ai}in=0, according to Each element B˜in(x) of the modified Bernstein basis is obtained from the corresponding element Bin(x) of the Bernstein basis using the following equality where the coefficients {˜bi}in=0 are ˜bi = n i bi
The right-hand side of (26) represents the MBF of polynomial p(x). The advantage of considering the MBF of p(x) instead of its Bernstein form is that the coefficients {˜bi}in=0 are integer numbers and they can therefore be represented with no approximation, provided that integer numbers can be represented with arbitrary precision. As a matter of fact, since the coefficients {˜bi}in=0 can be written as n i i nj = j n − j i − j i ˜bi = X j=0 n − j i − j cj i ∈ [0..n] from which it can be easily deduced that they are integer numbers. It is worth noting that binomial coefficients are integer numbers and they can be computed without divisions. As a matter of fact, the following recursive formula holds for binomial coefficients, when n ≥ j > 0 n j = n − 1 j + n − 1 j − 1 .
The recursion terminates because the binomial coefficient equals 1 for j = 0.
Observe that (23) holds for the coefficients {bi}in=0 and not for {˜bi}in=0. However, as previously discussed, we are only interested in the signs of the coefficients {bi}in=0, which are the same as that of corresponding {˜bi}in=0 since they only differ by a positive multiplicative factor. For this reason, it is possible to focus only on the integer coefficients {˜bi}in=0. (25) (26) (27) (28) (29) (30) x = (x1, . . . , xk) where N = (n1, . . . , nk) is a multi-index, and operations on multi-indices are defined component-wise. Observe that the inequality I ≤ N among multi-indices is also component-wise, namely it corresponds to As in the univariate case, we start by considering a single constraint, namely and we assume that each variable is defined over a finite interval. More precisely, we denote as x = (x1, . . . , xk) and x = (x1, . . . , xk) the extremes of the box I = [x, x] where the polynomial p(x) is defined. In order to find the explicit expression in the Bernstein basis of the polynomial p(x) defined in (32), it is necessary to consider an affine change of variable between x and a new variable t defined over the unit box [0, 1]k. Using the affine transformation A more sophisticated notation needs to be introduced in the multivariate case.
k We denote as {xj }j=1 the considered variables where k ≥ 1 and j is an index that we associate with each variable. We define the vector x of all variables and we define a multi-index I as I = (i1, . . . , ik), so that xI = xi11 · . . . · xikk .
A k-variate polynomial p(x) with real coefficients can be written as it is possible to show that
xi = xi + (xi − xi)ti where the coefficients cI are related to aI according to and Then, it is possible to write the polynomial p(x) defined in (32) as cI =
N X aI L=I
L xL−I (x − x)I
I L I l1 i1 = · . . . · lk .
ik and Binjj (xj) = nj (xj − xj)ij (xj − xj )nj−ij ij
(xj − xj)nj The coefficients bI are related to cI (and, hence, to aI ) according to ij ∈ {0, . . . , nj} j ∈ {1, . . . , k}. where the polynomials {BIN (x)}I≤N represent the Bernstein basis of the considered polynomial space with the following explicit expression
BIN (x) = Bin11(x1) · Bin22(x2) · . . . · Binkk (xk) Notably, the range of the polynomial p(x) satisfies
I bI = X NJ cJ
J≤I J
I ≤ N. range(p) ⊆ [min{bI}, max{bI }].
I≤N I≤N As in the univariate case, such a well-known property of the Bernstein form is useful to study polynomial constraints because it provides effective lower and upper bounds of a polynomial, and it can be used to study its sign in a box.
We now aim at deriving the MBF of the polynomial p(x) starting from its Bernstein form in (39). In order to do so, let us consider the modified Bernstein basis given by the following set of polynomials
B˜IN (x) = B˜in11(x1) · B˜in22(x2) · . . . · B˜inkk (xk) where
B˜injj (xj ) = (xj − xj)ij (xj − xj)nj−ij
(xj − xj)nj Using such a basis, the polynomial p(x) can be written as ij ∈ {0, . . . , nj} j ∈ {1, . . . , k}.
n p(x) = X ˜bI B˜IN (x)
i=0 where the coefficients {˜bI }I≤N are ˜bI = X
J≤I
N I
I NJ cJ = X J J≤I
N − J I − J cJ
I ≤ N.
Observe that, since all factors involved in (45) are integer numbers, the coefficients {˜bI}I≤N of the polynomial p(x) in the MBF are integer numbers. Moreover, the binomial coefficients in (45) can be computed with no divisions using (30). This allows evaluating the coefficients {˜bI}I≤N with no approximation. As in the univariate case, even if condition (42) holds for only the coefficients {bI}I≤N , we can focus on the coefficients {˜bI}I≤N since we are only interested in their signs. (40) (41) (42) (43) (44) (45) and such a constraint can be rewritten in canonical form and interpreted as p(x) = x2 − 9 ≥ 0 x ∈ [−10..10].
(46) The trace of the solving process for such a constraint is shown in Figure 2 where light blue nodes are non-terminal nodes, green nodes are terminal nodes with consistent domains and red nodes are terminal nodes with inconsistent domains. In non-terminal nodes, we denote as ˜land u˜ a negative value and a positive value of {˜bi}i=0, respectively. The figure shows that after 7 branches the algorithm 2 computes the correct answer x ∈ [−10.. − 3] ∪ [3..10]. (47) The second motivating example shown in Section 1 is
X in -10..10, Y in -10..10, X*Y #>= 21 and such a constraint can be rewritten in canonical form and interpreted as p(x) = xy − 21 ≥ 0 x ∈ [−10..10] y ∈ [−10..10].
(48) The trace of the solving process is too complex to be shown because it counts 41 branches. After all such branches, the algorithm computes the correct answer
Finally, the following is a well-known example used to show that CLP(FD) solver available in SWI-Prolog (version 6.6.6) sometimes does not capture inconsistent constraints
X in 1..2, Y in 1..2, Z in 1..2, X #\=Y, X #\=Z, Z #\=Y. The expected result is false, but the solver does not change the constraint, even if all different is used instead of the three inequalities. Such a constraint can be written in canonical form and interpreted as x2 − 2xy + y2 − 1 ≥ 0 x2 − 2xz + z2 ≥ 0 y2 − 2yz + z2 ≥ 0.
(50) Figure 3 shows the solving process of the proposed algorithm, where, after 5 branches, all terminal nodes are red, i.e., the constraint is inconsistent. In nonterminal nodes, we denote as l˜ and u˜ a negative value and a positive value of {˜bI}I≤N , respectively, relative to a constraint.
In this paper we presented a solver for polynomial constraints over finite domains. Such constraints are expressed as polynomials with integer coefficients whose variables are defined over finite subsets of integer numbers. The proposed algorithm uses the MBF of polynomials to effectively compute the signs of suitable lower and upper bounds of polynomials. It uses such signs to reason on constraints and to remove inadmissible values from the domains of variables. It is worth noting that the algorithm only involves algebraic manipulations over integer numbers because the coefficients of the MBF of a polynomial with integer coefficients are provably integer. The availability of arbitrary precision integer arithmetic is therefore sufficient to ensure that no approximation is involved in the solution process.
The proposed algorithm is open for further developments on possible optimizations, at least, for two important choices that impact significantly on the performances of the solution process. First, the choice of the variable for splitting the current box is crucial. Second, once a variable is selected, the choice of simply halving the current box can be too simplistic. For both choices a number of static and dynamic techniques can be found in the literature (see, e.g., [3]) and we think that both choices could benefit from the peculiar form of constraints that we address. In particular, we plan to investigate how the MFB could help in the computation of the gradient of polynomials, which is a valuable information for both such choices.
The current implementation of the proposed algorithm is available for download (cmt.dmi.unipr.it/software/clppolyfd.zip) in terms of a reusable Java library that is also usable as a drop-in replacement of the CLP(FD) solver that ships with SWI-Prolog.