Satisfaction of polynomial constraints over finite domains using function values Federico Bergenti and Stefania Monica Dipartimento di Scienze Matematiche, Fisiche e Informatiche Università degli Studi di Parma, 43124 Parma, Italy {federico.bergenti,stefania.monica}@unipr.it Abstract. This paper shows how the solutions of constraint satisfac- tion problems that involve only polynomial constraints over finite do- mains can be enumerated by computing the values of related polynomial functions at appropriate points. The proposed algorithm first transforms constraints, which are expressed as equalities, inequalities, and disequal- ities of polynomials with integer coefficients and integer variables, into a canonical form that uses only inequalities. Then, starting from a bound- ing box, which is supposed to be known, the algorithm recursively sub- divides the box into disjoint boxes and it records boxes whose elements satisfy all constraints. The subdivision is driven by the study of the sign of polynomial functions over boxes, which is performed by means of a method that uses only the coefficients of polynomials and the values of functions at the corners of boxes. Keywords: polynomial constraints over finite domains, subdivision al- gorithms, Taylor’s theorem 1 Introduction Polynomial constraints are ubiquitous, and they find relevant applications in virtually all fields of science and engineering (see, e.g., [5, 13] and related liter- ature). In this paper we study a specific form of polynomial constraints, which we call polynomial constraints over finite domains [1], in the scope of constraint logic programming and constraint satisfaction problems. Polynomial constraints over finite domains are constraints expressed as equalities, inequalities, and dise- qualities of polynomials with integer coefficients whose variables take values from finite subsets of the integers. Normally, reasoning on polynomial constraints over finite domains assumes that an initial approximation of the domains of variables is available in terms of a suitable bounding box. We have already presented results on the use of the properties of Bernstein polynomials [4, 12] to reason on constraints satisfaction problems involving only polynomial constraints over finite domains [2, 3]. We have also proposed an al- gorithm to enforce hyper-arc consistency of this type of constraints using a sub- division algorithm [1], which also uses the properties of Bernstein polynomials. In this paper, the subdivision framework that we adopted in previous studies is used to enumerate the solutions of constraint satisfaction problems that in- volve only polynomial constraints over finite domains by computing the values of related polynomial functions at appropriate points, with no use of Bernstein polynomials. In detail, the proposed algorithm first rewrites constraints, which are expressed as equalities, inequalities and disequalities of polynomials with integer coefficients and integer variables, into a canonical form that uses only inequalities. Then, starting from an initial approximation of the domains of vari- ables, which is supposed to be known in term of a bounding box, the algorithm recursively subdivides the box and it records boxes whose elements satisfy all constraints. Thanks to the adopted canonical form of constraints, the subdivision is driven by the study of the sign of polynomial functions over boxes. Section 3 presents results that show how the computation of the values of a polynomial function at the corners of a box provides relevant information on the sign of the function over the box. In detail, Section 3 describes a sufficient condition for a polynomial function to be nonnegative over a box by using only the coefficients of the polynomial and the values of the function at the corners of the box, which is exactly what is needed to drive the subdivision. If the sufficient condition is met then the elements of the considered box satisfy all constraints, otherwise the box should be split into disjoint boxes and processed recursively. This paper is organised as follows. Section 2 recalls adopted canonical form of constraints and it outlines the subdivision algorithm, which is based on the possibility of identifying boxes where polynomial functions are nonnegative. Sec- tion 3 discusses how the computation of the value of a polynomial function at the corners of a box provides a sufficient condition for the function to be non- negative over the box. Section 4 presents illustrative examples meant to clarify the proposed algorithm. Finally, Section 5 concludes the paper and it outlines open questions and future research directions. 2 Reasoning on polynomial constraints over boxes Subdivision algorithms have been used extensively to reason on polynomial con- straints mostly because effective bounds on the range of polynomial functions can be computed using the properties of Bernstein polynomials (see, [6] for a survey of major results and a historical perspective). Starting from a given box, lower and upper bounds for a polynomial function over the box can be computed by using the properties of Bernstein polynomials. Such bounds give relevant in- formation on the sign of the function over the box, and they can be used to decide whether to recursively subdivide the box or not. This approach, which is infor- mally known as Bernstein algorithm [9], is very general, and it has been used, for example, to solve systems of polynomial equations [10, 14] and to support global optimisation [11, 15]. Besides applications in mixed-integer nonlinear optimisation [16], the use of the Bernstein algorithm has been mostly limited to reason on polynomials con- straints whose variables take values from the reals. Only recently, the properties of Bernstein polynomials have been used to reason on constraints whose vari- ables take values from the integers, as follows. First, constraints are rewritten into a canonical form, which turns a generic polynomial constraint with integer coefficient whose variables take values from the integers into an inequality. Given a polynomial constraint of the form p(x) q(x), (1) where p(x) and q(x) are polynomials with integer coefficients in n ∈ N variables x = (x1 , x2 , . . . , xn ), and ∈ {=, 6=, <, ≤, >, ≥}, a semantically equivalent con- straint of the form r(x) ≥ 0, (2) where r(x) is a polynomial with integer coefficients, can be written by applying simple rules detailed in [1]. The possibility of transforming a generic polynomial constraint (1) into a semantically equivalent constraint of the specific form (2) can be used to turn the problem of searching for the values of variables that satisfy (1) into the problem of searching for the values of variables that make the polynomial function r(x) nonnegative. This problem can be effectively solved using subdivision under the assumption that an initial approximation of the domains of variables is available. Actually, under this assumption, in [2, 3] we used subdivision to study the satisfiability of polynomial constraints over finite domains, while in [1] we used subdivision to make polynomial constraints over finite domains hyper-arc consistent. In both cases, we adopted the properties of Bernstein polynomials to study the signs of polynomial functions over boxes. In detail, the use of Bernstein polynomials to express a polynomial function over a given box allows the immediate identification of lower and upper bounds for the function over the box. Such bounds can be used to evaluate the sign of the function over the box, which can be recursively subdivided into disjoint boxes if the sign of the function changes over the box. This method is particularly appealing when variables are restricted to take values from finite subsets of the integers, because such a choice guarantees that the method always terminates, and also that results can be computed over the integers with no approximations, provided that arbitrary precision integer arithmetic is available. The subdivision algorithm that we have already adopted to reason on poly- nomial constraints over finite domains does not mandate the use of a specific algorithm to study the sign of polynomial functions over boxes, and it can be gen- eralised to a framework to reason on polynomial constraints over finite domains. In order to obtain working algorithms, the framework needs to be completed with a means to decide if a polynomial function is nonnegative over a box. This can be accomplished by means of any algorithm capable of computing suitable lower and upper bounds for a given polynomial function over a given box. If the lower bound is nonnegative, then the polynomial function is nonnegative over the box. Conversely, if the upper bound is negative, then the polynomial func- tion is negative over the box. Note that the trivial approach of evaluating the polynomial function at all points of the box, which is always finite, is surely pos- sible, but we expect that the use of lower and upper bounds would be beneficial especially for large boxes. Algorithm 1 Function enumerate(.,.) takes a set of constrains C in canonical form and a box B 6= ∅ over which the satisfiability of constraints in C is studied, and it returns a subset of B whose elements satisfy all constraints. 1: function enumerate(C, B) 2: input C : a set of constraints . current set of constraints in canonical form 3: input B : a box in Zn . nonempty box currently analysed 4: output S ⊆ B . subset of B in which all constraints are satisfied 5: if C = ∅ then 6: S←B 7: else 8: select c ∈ C with c = (p(x) ≥ 0) 9: (l, u) ← bounds(p, B) 10: if u < 0 then 11: S←∅ 12: else if l ≥ 0 then 13: S ← enumerate(C\{c}, B) 14: else 15: select xj with Ij = [xj ..xj ] for B = [x1 ..x1 ] × [x2 ..x2 ] × · · · × [xn ..xn ] 16: if xj = xj then 17: q(x) ← p(x) with xj replaced by xj 18: S ← enumerate(C\{c} ∪ {q(x) ≥ 0}, B) 19: else 20: select s ∈ Ij with s 6= xj 21: S 0 ← enumerate(C, Bj→[xj ..s] ) 22: S 00 ← enumerate(C, Bj→[s+1..xj ] ) 23: S ← S 0 ∪ S 00 24: end if 25: end if 26: end if 27: end function Algorithm 1 shows a pseudocode of the proposed framework to enumerate the solutions of constraint satisfaction problems that involve only polynomial constraints over finite domains. Function enumerate(.,.) takes a set of con- straints C in canonical form and a box B ⊂ Zn , with B 6= ∅, over which the signs of polynomials in C are studied. First, enumerate(.,.) selects a constraint and it uses the unspecified function bounds(.,.) to compute appropriate lower and upper bounds of the related polynomial function over B. If computed upper bound is negative, then the constraint is unsatisfiable in the box. Conversely, if computed lower bound is nonnegative, then the selected constraint is satisfiable for all elements of the box. If both conditions are not met, then a variable is selected and the current approximation of its domain is properly split into two disjoint intervals, which are processed recursively. The notation Bj→[a..b] is used to denote a box obtained from box B by replacing interval Ij , which corresponds to the current approximation of the domain of the j−th variable, with interval [a..b]. Note that the pseudocode in Algorithm 1 considers also the special case of variables whose domains are reduced to singleton sets. In this case, variables can be substituted in polynomials with appropriate values, so that the number of variables in polynomials is reduced. Finally, note that function enumerate(.,.) is first invoked with the initial set of constraints and the initial approximation of the domains of variables, which is expressed as a nonempty bounding box. 3 Study of the sign of polynomials using function values In this section we present a method to study the sign of a polynomial function over a given box that is based on a result originally proposed for univariate polynomials in [18], which was then extended to bivariate polynomials in [8]. The most general result, which holds for multivariate polynomials, is shown, but not proved for the sake of brevity, in Proposition 2. First, the case of univariate polynomials is discussed to explain the proposed method in a simple setting. Then, the proposed extension to multivariate polynomials is shown using the elegant notation of multi-indices. Such an extension can be used as a basis to implement function bounds(.,.) of the pseudocode in Algorithm 1. 3.1 Univariate polynomials Let p : R 7→ R be a polynomial function, and B ⊂ R be a closed interval. We are interested in computing a lower bound m ∈ R and an upper bound m ∈ R of p over B such that m ≤ min p(x) max p(x) ≤ m. (3) x∈B x∈B In order to compute m and m, we recall the Taylor’s theorem for functions in one variable, which is rephrased here to fix notation and because there is little consensus on what should be called Taylor’s theorem. The theorem is classic and it is shown without proof. Theorem 1 (Taylor’s theorem in one variable). Let f : R 7→ R be of class C (k+1) , with k ∈ N, on an open interval S ⊆ R. If a ∈ S, h ∈ R and a + h ∈ S, then k X D(i) f (a) i f (a + h) = h + Ra,k (h), (4) i=0 i! where D(i) f is the i−th derivative of f (with D(0) f = f ), and the remainder Ra,k (h) is given in Lagrange form by D(k+1) f (a + ch) (k+1) Ra,k (h) = h (5) (k + 1)! for some c ∈ (0, 1). The application of the Taylor’s theorem to a polynomial function on a closed interval easily proves the following corollary. Corollary 1. Let p : R 7→ R be a polynomial function and B = [b, b] ⊂ R be a closed interval. Then, for all x ∈ B, and for all k ∈ N k X D(i) p(b) p(x) = (x − b)i + Rb,k (x − b), (6) i=0 i! where the remainder Rb,k (h) is given in Lagrange form by D(k+1) p(b + ch) (k+1) Rb,k (h) = h (7) (k + 1)! for some c ∈ (0, 1). Proof. Recalling that any polynomial function is of class C ∞ on any open inter- val S ⊆ R with B ⊂ S, the proof of the corollary follows immediately from the Taylor’s theorem. t u In particular, if we truncate to k = 2, we have that D(2) p(b + c(x − b)) p(x) = p(b) + D(1) p(b)(x − b) + (x − b)2 (8) 2 for some c ∈ (0, 1). The truncated expansion (8) helps to prove the following proposition (see also [18]). Proposition 1. Let p : R 7→ R be a polynomial function of degree l ∈ N l X p(x) = ai xi (9) i=0 with coefficients {ai }li=0 ⊂ R, and let U = [0, 1] ⊂ R with m = min p(x) m = max p(x). (10) x∈U x∈U Then min p(x) − δ ≤ m m ≤ max p(x) + δ, (11) x∈{0,1} x∈{0,1} where l 1X δ= (i − 1)i|ai |. (12) 8 i=2 Proof. Let η ∈ U be such that p(η) = m, then if η = 0 or η = 1 the proposition is immediately proved because δ ≥ 0 and therefore min p(x) − δ ≤ m. (13) x∈{0,1} Conversely, if η 6= 0 and η 6= 1, then D(1) p(η) = 0. In this case, using the fact that the choice of interval U ensures that 0 ≤ xi ≤ 1, (14) with i ∈ N, and x ∈ U , the expansion of the second derivative in (8) is sufficient to prove the proposition. In detail, if η ≤ 21 , then D(2) p(cη) 2 m = p(0) + η (15) 2 for some c ∈ (0, 1). But, l X l X D(2) p(x) = i(i − 1)ai x(i−2) ≤ i(i − 1)|ai |, (16) i=2 i=2 for all x ∈ U , which proves the proposition for η ≤ 12 . For the case η > 21 , the proposition is proved by choosing the other endpoint of interval U to expand the function, which leads to an equality similar to (15). The case of m is treated with similar considerations. t u Proposition 1 ensures that, for the case of univariate polynomials, suitable bounds to support the pseudocode in Algorithm 1 can be obtained by computing the values of polynomial functions at the endpoints of considered intervals. Note that the assumption x ∈ U in Proposition 1 is not restrictive because, in order to consider a polynomial function p : R 7→ R in variable x which takes values from interval B = [b, b] ⊂ R, it is sufficient to change the variable of p to t ∈ U such that x = (b − b)t + b, (17) and to consider the coefficient of the new polynomial in variable t. In this case, we prefer to use the notation δB in (12), instead of simply δ, because the coefficients of the new polynomial make δB depend on B. 3.2 Multivariate polynomials In order to elegantly generalise Proposition 1 to the case of polynomials with n ∈ N variables, we need to introduce the notation of multi-indices. A multi- index I ∈ Nn is an n-tuple of the form I = (i1 , i2 , . . . , in ), (18) and we adopt the common notation of using an uppercase letter to denote a multi-index and the same letter, but lowercase and with an index, to denote its components. Given a multi-index I ∈ Nn , the following abbreviations n X n Y |I| = ij I! = ij ! (19) j=1 j=1 are often used to denote the order of I, and the factorial of I, respectively. Given two multi-indices I ∈ Nn and J ∈ Nn , they can be added I + J = (i1 + j1 , i2 + j2 , . . . , in + jn ), (20) and they can be compared I ≤ J ⇐⇒ i1 ≤ j2 ∧ i2 ≤ j2 ∧ · · · ∧ in ≤ jn , (21) which justify the following abbreviations, with k ∈ N, X j1 X X j2 jn X (·) = ··· (·) (22) I≤J i1 =0 i2 =0 in =0 X k X X k k X (·) = ··· (·) . (23) |I|≤k i1 =0 i2 =0 in =0 | {z } i1 +i2 +···+in ≤k Given a multi-index I ∈ Nn and v ∈ Rn , the following abbreviation is often used n i Y vI = vjj . (24) j=0 Finally, given a generic function f : Rn 7→ R in variables x = (x1 , x2 , . . . , xn ), the following notation based on multi-indices is used to denote elegantly the partial derivatives of f ∂ |I| f ∂ I f = ∂ i1 ∂ i2 · · · ∂ in f = . (25) ∂x1 i1 ∂x2 i2 · · · ∂xn in The notation of multi-indices can be used to state the following theorem, which is then used to study the sign of a polynomial function over a given box by computing the value of the function at the corners of the box. The theorem generalises the Taylor’s theorem recalled above, and it is shown without proof. Theorem 2 (Taylor’s theorem in several variables). Let f : Rn 7→ R be of class C (k+1) , with k ∈ N, on an open convex S ⊆ Rn . If a ∈ S, h ∈ Rn and a + h ∈ S, then X ∂ I f (a) f (a + h) = hI + Ra,k (h), (26) I! |I|≤k where the remainder Ra,k (h) is given in Lagrange form by X ∂ I f (a + ch) I Ra,k (h) = h (27) I! |I|=k+1 for some c ∈ (0, 1). Let p : Rn 7→ R be a multivariate polynomial function in n ∈ N variables x = (x1 , x2 , . . . , xn ), the adoption of multi-index notation allows expressing p as X p(x) = aI xI , (28) I≤L n where L ∈ N is the multi-degree of the polynomial function, and {aI }I≤L ⊂ R is the set of its coefficients. Given that we are interested in studying the sign of polynomial functions in boxes, we introduce the following notation to denote boxes. The box B built using v = (v 1 , v 2 , . . . , v n ) ∈ Rn and v = (v 1 , v 2 , . . . , v n ) ∈ Rn is denoted as B = [v, v] = [v 1 , v 1 ] × [v 2 , v 2 ] × · · · × [v n , v n ], (29) and it is the empty set if and only if for some 0 ≤ i ≤ n, v i > v i . Adopted notations can be used to state a generalisation of Corollary 1 for multivariate polynomials. Corollary 2. Let p : Rn 7→ R be a polynomial function in n ∈ N variables, and B = [b, b] ⊂ Rn be a box. Then, for all x ∈ B, and for all k ∈ N X ∂ I p(b) p(x) = (x − b)I + Rb,k (x − b) (30) I! |I|≤k and the remainder Rb,k (h) is given in Lagrange form by X ∂ I p(b + ch) I Rb,k (h) = h (31) I! |I|=k+1 for some c ∈ (0, 1). Proof. Recalling that any polynomial function is of class C ∞ for any open convex S ⊆ Rn with B ⊂ S, the proof of the corollary follows immediately from the Taylor’s theorem. t u Following the same path traced for the case of univariate polynomials, we can prove the following proposition, which is sufficient to compute appropriate lower and upper bounds to support the proposed algorithm. The proposition is shown without proof for the sake of brevity. Proposition 2. Let p : Rn 7→ R be a polynomial function with multi-degree L ∈ N in n ∈ N variables x = (x1 , x2 , . . . , xn ) X p(x) = aI xI (32) I≤L n n and U = [0, 1] ⊂ R , with m = min p(x) m = max p(x). (33) x∈U x∈U Then min p(x) − δ ≤ m m ≤ max p(x) + δ, (34) x∈C x∈C where C = {0, 1}n is the set of the corners of U and 1X δ= |I|(|I| − 1)|aI |. (35) 8 I≤L Note that the assumption x ∈ U in Proposition 2 is not restrictive. In order to consider a polynomial function p : Rn 7→ R in variable x which takes values from box B = [b, b] ⊂ Rn , it is sufficient to change the variable of p to t ∈ U such that x = (b − b)t + b, (36) and to consider the coefficient of the new polynomial in variable t. In this case, we prefer to use the notation δB in (35), instead of simply δ, because the coefficients of the new polynomial make δB depend on B. 4 Illustrative examples This section shows how the proposed algorithm can be applied to some illus- trative examples. Note that a systematic validation of the proposed algorithm with specific benchmarks is still in progress, and the following examples are only intended to clarify the algorithm. The first example is a quadratic constraint with only one variable p1 (x) = −x2 − 50 ≥ 0 x ∈ [1..25]. (37) The constraint is unsatisfiable in the given interval, and the proposed algorithm proves it by subdividing the initial interval into 4 intervals, the largest of which contains 12 values. Figure 1 shows the subdivision produced by the algorithm. The second example is an unsatisfiable constraint with two variables p2 (x, y) = −x2 − y 2 − 100 ≥ 0 x ∈ [1..25], y ∈ [1..25]. (38) The proposed algorithm needs to subdivide the initial box into 12 boxes, the largest of which contains 144 values, to prove that the constraint is unsatisfiable. Figure 2 shows the subdivision of the box that the algorithm produces. The following examples are relative to larger domains and the detailed sub- division of the initial approximations of domains are too complex to be shown. The following example is a satisfiable constraint with one variable p3 (x) = x2 − 16 ≥ 0 x ∈ [−100..100]. (39) The proposed algorithm builds a search tree made of 65 nodes to compute the following answer x ∈ [−100.. − 4] ∪ [4..100]. (40) 4 3 6 12 1 4 5 7 8 13 14 25 Fig. 1. Subdivision of interval I = [1..25] that the proposed algorithm produces to prove that p1 (x) = −x2 − 50 ≥ 0 is unsatisfiable in I. The cardinality of each interval is also shown. Given that each node requires the computation of the value of p3 at the two endpoints of the related interval, 130 evaluations of p3 are needed. Note that a brute force investigation requires 201 evaluations of the polynomial. The fourth example is still with only one variable p4 (x) = x2 − 50x + 1 ≥ 0 x ∈ [−100..100]. (41) In this case the answer x ∈ [−100..0] ∪ [50..100] (42) is returned after computing a search tree made of 61 nodes, which requires 122 evaluations of p4 , against the 201 evaluations needed by brute force investigation. The fifth illustrative example involves two variables p5 (x, y) = xy − 210 ≥ 0 x ∈ [−100..100], y ∈ [−100..100]. (43) The proposed algorithm needs to traverse a search tree made of 1, 909 nodes to compute the answer x ∈ [−100.. − 3] ∪ [3..100] y ∈ [−100.. − 3] ∪ [3..100]. (44) Note that the computation of the answer required 7, 636 evaluations of p5 , four for each considered box, while the brute force scan of the domains of variables requires 40, 401 evaluations of the polynomial. Finally, the last example involves two variables p6 (x, y) = x + xy − 1000 ≥ 0 x ∈ [−100..100], y ∈ [−100..100]. (45) In this case the answer x ∈ [−100.. − 11] ∪ [10..100] y ∈ [−100.. − 11] ∪ [9..100] (46) is computed after traversing a search tree made of 1, 347 nodes, which required 5, 388 evaluations of p6 , against the 40, 401 evaluations of brute force search. Note that discussed examples are processed using the current implementation of the algorithm, which intentionally does not use heuristics to select which variable to split, and where to split the domain of the selected variable. It is expected that the introduction of heuristics to drive the search would bring notable performance improvements. 25 42 20 72 144 19 42 14 13 42 36 72 8 7 21 5 4 42 42 42 28 1 7 8 13 14 19 20 25 Fig. 2. Subdivision of box B = [1..25] × [1..25] that the proposed algorithm produces to prove that p2 (x, y) = −x2 − y 2 − 100 ≥ 0 is unsatisfiable in B. The cardinality of each box is also shown. 5 Conclusions This paper presented an algorithm to enumerate the solutions of constraint sat- isfaction problems that involve only polynomial constraints over finite domains. Such constraints are expressed as equalities, inequalities, and disequalities of polynomials with integer coefficients whose variables take values from finite sub- sets of the integers. The proposed algorithm works under the assumption that an initial approximation of the domains of variables is available in terms of a bounding box. First, the algorithm rewrites all constraints into a canonical form in which only inequalities are used. Then, it recursively subdivides the initial bounding box into disjoint boxes, and it records boxes whose elements satisfy all constraints. The canonical form of constraints ensures that a box contains only elements that satisfy all constraints if and only if the polynomial functions re- lated to constraints in canonical form are all nonnegative over the box. In order to verify if a polynomial function is nonnegative over a given box, the algorithm uses lower and upper bounds described in Proposition 2. Such bounds are com- puted by considering the coefficients of the polynomial function and its values at the corners of the considered box. The algorithm was tested on illustrative examples and preliminary results are encouraging, but a systematic validation with specific benchmarks is still in progress. We expect that the proposed algorithm would outperform our previous proposals based on the properties of Bernstein polynomials because of the in- herent complexity of computing Bernstein coefficients, even using sophisticated techniques (see, e.g., [7, 17, 19, 20]). Moreover, it is evident that the performance of the algorithm on single problems severely depends on the details of the sub- division, which are deliberately not considered in this paper. In particular, it is expected that nontrivial heuristics to decide which domain to split and where to split it would be largely beneficial in practice. Heuristics could use the specific form of treated constraints to reason on the next variable to select for subdivi- sion, and on where to subdivide its domain, because the gradient of polynomial functions can be computed easily. The proposed algorithm is implemented as a reusable Java library that can be accessed as a Prolog module from SWI-Prolog [22]. The current implementation of the Java library and of the interface to SWI-Prolog can be downloaded in a single package (cmt.dmi.unipr.it/software/clppolyfd.zip). To ease experi- mentation, the Prolog module is fully compatible with the CLP(FD) module [21] of SWI-Prolog, and it can be configured to use the proposed algorithm or the algorithm based on Bernstein polynomials discussed in [1–3]. References 1. Bergenti, F., Monica, S.: Hyper-arc consistency of polynomial constraints over finite domains using the modified Bernstein form. Annals of Mathematics and Artificial Intelligence 80(2), 131–151 (2017) 2. Bergenti, F., Monica, S., Rossi, G.: Polynomial constraint solving over finite do- mains with the modified Bernstein form. In: Fiorentini, C., Momigliano, A. (eds.) 31st Italian Conference on Computational Logic. CEUR Workshop Proceedings, vol. 1645, pp. 118–131. RWTH Aachen (2016) 3. Bergenti, F., Monica, S., Rossi, G.: A subdivision approach to the solution of polynomial constraints over finite domains using the modified Bernstein form. In: Adorni, G., Cagnoni, S., Gori, M., Maratea, M. (eds.) AI*IA 2016 Advances in Artificial Intelligence. Lecture Notes in Computer Science, vol. 10037, pp. 179– 191. Springer International Publishing (2016) 4. Bernstein, S.N.: Démonstration du théorème de Weierstrass fondée sur le calcul des probabilités. Communications de la Société Mathématique de Kharkov 2:XIII(1), 1–2 (1912) 5. Borralleras, C., Lucas, S., Oliveras, A., Rodrı́guez-Carbonell, E., Rubio, A.: SAT modulo linear arithmetic for solving polynomial constraints. Journal of Automated Reasoning 48(1), 107–131 (2010) 6. Farouki, R.T.: The Bernstein polynomial basis: A centennial retrospective. Computer-Aided Geometric Design 29(6), 379–419 (2012) 7. Farouki, R.T., Rajan, V.T.: Algorithms for polynomials in Bernstein form. Computer-Aided Geometric Design 5(1), 1–26 (1988) 8. Garloff, J.: Convergent bounds for the range of multivariate polynomials. In: Nickel, K. (ed.) Interval Mathematics 1985. Lecture Notes in Computer Science, vol. 212, pp. 37–56. Springer International Publishing (1986) 9. Garloff, J.: The Bernstein algorithm. Interval Computations 2, 154–168 (1993) 10. Garloff, J., Smith, A.P.: Solution of systems of polynomial equations by using Bernstein expansion. In: Alefeld, G., Rohn, J., Rump, S., Yamamoto, T. (eds.) Symbolic Algebraic Methods and Verification Methods. pp. 87–97. Springer Inter- national Publishing (2001) 11. Grimstad, B., Sandnes, A.: Global optimization with spline constraints: A new branch-and-bound method based on B-splines. Journal of Global Optimization 65(3), 401–439 (2016) 12. Lorentz, G.G.: Bernstein Polynomials. University of Toronto Press, Toronto (1953) 13. Lucas, S., Navarro-Marset, R.: Comparing CSP and SAT solvers for polynomial constraints in termination provers. Electronic Notes in Theoretical Computer Sci- ence 206, 75–90 (2008) 14. Mourrain, B., Pavone, J.: Subdivision methods for solving polynomial equations. Journal of Symbolic Computation 44(3), 292–306 (2009) 15. Nataraj, P., Arounassalame, M.: A new subdivision algorithm for the Bernstein polynomial approach to global optimization. International Journal of Automation and Computing 4(4), 342–352 (2007) 16. Patil, B.V., Nataraj, P.S.V., Bhartiya, S.: Global optimization of mixed-integer nonlinear (polynomial) programming problems: The Bernstein polynomial ap- proach. Computing 94(2), 325–343 (2012) 17. Ray, S., Nataraj, P.: An efficient algorithm for range computation of polynomials using the Bernstein form. Journal of Global Optimization 45, 403–426 (2009) 18. Rivlin, T.J.: Bounds on a polynomial. Journal of Research of the National Bureau of Standards 74B(1), 47–54 (1970) 19. Sánchez-Reyes, J.: Algebraic manipulation in the Bernstein form made simple via convolutions. Computer-Aided Design 35, 959–967 (2003) 20. Smith, A.P.: Fast construction of constant bound functions for sparse polynomials. Journal of Global Optimization 43(2), 445–458 (2009) 21. Triska, M.: The finite domain constraint solver of SWI-Prolog. In: Schrijvers, T., Thiemann, P. (eds.) Functional and Logic Programming. Lecture Notes in Com- puter Science, vol. 7294, pp. 307–316. Springer International Publishing (2012) 22. Wielemaker, J., Schrijvers, T., Triska, M., Lager, T.: SWI-Prolog. Theory and Practice of Logic Programming 12(1–2), 67–96 (2012)