Solving congruence equations using Bernstein forms
aa r X i v : . [ m a t h . O C ] A ug Solving congruence equations using Bernstein forms
C´esar Massri a,1, ∗ , Manuel Dubinsky b,2 a Department of Mathematics, University CAECE, Buenos Aires, Argentina b Computer Engineering Department, University of Avellaneda, Argentina
Abstract
We present a subdivision method to solve systems of congruence equations. This method isinspired in a subdivision method, based on Bernstein forms, to solve systems of polynomialinequalities in several variables and arbitrary degrees. The proposed method is exponential inthe number of variables.
Keywords:
Congruence equations, Bernstein form, Integer polynomial programming, Solver
Introduction
Overview of the problem
An important aspect within the study of algebraic varieties over characteristic p is to knowthe number of points that a variety has. This led to the statement of Weil’s conjectures proved byDeligne thanks to the analysis made by Grothendieck and others. In this article we are interestedin computing, not only the number of points in an algebraic variety modulo p but also, the co-ordinates of these points. We propose a computational method to find solutions to the followingsystem of congruence equations h ( x , . . . , x n ) ≡ p ) ... h r ( x , . . . , x n ) ≡ p r )where h , . . . , h r ∈ Z [ x , . . . , x n ] and p , . . . , p r ∈ Z ≥ . To achieve this goal we first analyze aknown method in nonlinear programming and then adapt it to our needs.Polynomial programming is the process of solving a system of equalities and inequalitiesbetween polynomial functions along with an objective function to be maximized (or minimized).Typically, such a system can be described as:maximize q ( x )subject to q i ( x ) ≥ , ≤ i ≤ r . ∗ Address for correspondence: Department of Mathematics, FCEN, University of Buenos Aires, Argentina. Postaladdress: 1428. Phone number: 54-11-4576-3335.
Email addresses: [email protected] (C´esar Massri), [email protected] (Manuel Dubinsky) The author was fully supported by CONICET, IMAS, Buenos Aires, Argentina The author was fully supported by UndavCyT 2014 project.
Preprint submitted to arXiv October 6, 2018 he polynomials q i are assumed to have real coe ffi cients but in this article we study integerpolynomial programming assuming q , . . . , q r ∈ Q [ x , . . . , x n ] and x ∈ Z n .Since the resolution of Hilbert’s Tenth problem by Matiyasevich, [19], several results ap-peared proving that there cannot exist any general algorithm to solve integer polynomial pro-gramming (it is non computable). For example, in [16], the author proves that the problem ofminimizing a linear form over quadratic constraints in integer variables is not computable by a re-cursive function. Hence, in order to avoid this phenomena, it is necessary to add some constrainsto the original problem. Given that we are interested in applications to congruence equations webound the variables to some region D = [ a , b ] × . . . × [ a n , b n ] ⊆ R n restricting the problem ofadapting an integer polynomial programmingmaximize q ( x )subject to q i ( x ) ≥ , ≤ i ≤ r and x ∈ D ∩ Z n . to a solver of congruence equations h ( x ) ≡ p ) ... h r ( x ) ≡ p r )where x ∈ D ∩ Z n , D = [0 , k ] n and k is the least common multiple of p , . . . , p r , k = lcm( p , . . . , p r ). Existing work
There are many local and global methods to solve a Diophantine equation. The first thingto do is to check whether the problem has a solution locally. If it has none, the problem issolved since our equation has no solution, but if it has, the local information that we obtain mayalready completely solve the problem, or may not so it is necessary to study the equation globally.The global methods include sophisticated tools from number theory, for example factorizationover a number field, computations of class groups and unit groups, Diophantine approximationtechniques, linear forms in logarithms, elliptic curves, modular forms and Galois representationsand more.Regarding congruence equations we can say that there are general methods for solving bothlinear and quadratic congruence equations however solutions to the general polynomial congru-ence is intractable. Standard algorithms can be found at [27] and also at [4, 5]. Several softwarepackages implement these techniques, see for example MagmaSoft, SageMath, Kant / Kash.Integer polynomial programming is a subclass of the more general class
Mixed integer non-linear programming (MINLP). Approaches to solve these problems may be classified as eitherstochastic or deterministic. Stochastic techniques employ some element of randomness in theirsearch for the global optimum and consequently rely on a statistical argument to prove con-vergence to the global optimum, see [26]. Deterministic methods have the advantage that theyprovide a rigorous guarantee of global optimality of any solution produced within a specifiedtolerance. Examples of deterministic techniques are cutting plane methods, interval methods,primal-dual methods, outer approximation, generating functions, etc. For instance using gener-ating functions of convex sets, it is proven in [1] that there exists a polynomial-time algorithm2o find the maximum value of the objective function without knowing where this maximumis achieved. Also, a very successful methodology is the spatial branch and bound algorithm[18, 28]. The aim of branch and bound techniques, is to divide the problem into subclasses tobe solved with convex or linear approximations, which are powerful and flexible for modelingand solving decision problems. Despite their widespread use, few available software packagesprovide any guarantee of correct answers or certification of results. Possible inaccuracy is causedby the use of floating-point numbers. Floating-point calculations require the use of built-in tol-erances for testing feasibility and optimality, and can lead to calculation errors in the solution oflinear programming relaxations and in the methods used for creating cutting planes to improvethese relaxations. See in [6] where the author gave a method to guarantee the result avoidingnumerical errors.There is a large variety of commercial and academic noncommercial software, but only afew of those software packages solve general non-convex problems to global optimality. Fora review see [3], [2], [17, §
1] and / or [8]. It is worth mention that there is a large literaturereviewing di ff erent methods (and software) developed to solve convex and non-convex MINLPproblems, furthermore, every two years, Robert Fourer publishes a list of currently availablecodes in the field of linear and integer programming, the latest edition can be found at the webpage of OR / MS-Today, Software Survey section, [12].In this article we use Bernstein forms. The advantage of using Bernstein forms is that theyproduce robust and reliable algorithms. They are used to solve systems of polynomial equationsand also integer polynomial programming, see for example [13, 22, 23]. We review an algorithmusing Bernstein forms to solve integer polynomial programming and then we adapt it to solvecongruence problems.Before recalling the next theorem, let us introduce some notations. Let ( d , . . . , d n ) be a multi-degree, s = ( d + . . . ( d n +
1) and let D = [ a , a + k ] × . . . × [ a n , a n + k n ], ( a , . . . , a n ) ∈ Z n ,( k , . . . , k n ) ∈ N n . Choose some index j , 1 ≤ j ≤ n , and then divide into two halves the j -side of D , D Lj = j − Y i = [ a i , a i + k i ] × [ a j , a j + k j − ] × n Y i = j + [ a i , a i + k i ] , D Rj = j − Y i = [ a i , a i + k i ] × [ a j + k j − , a j + k j ] × n Y i = j + [ a i , a i + k i ] . Theorem.
For every q , . . . , q r ∈ Z [ x , . . . , x n ] ≤ ( d ,..., d n ) there exists a rectangular matrix v ∈ Q s × r such that the values of q i over D Lj (resp. D Rj ) are bounded by the minimum and the maximum of { w Lki } sk = (resp. { w Rki } sk = ), where w L = M Lj v, w R = M Rj v, ≤ i ≤ r. The matrices { M Lj , M Rj } nj = depend on ( d , . . . , d n ) .Also, the value of q i at ( a , . . . , a n ) is equal to w L i and the value of q i at ( a , . . . , a j − , a j + k j − , a j + , . . . , a n ) is equal to w R i , ≤ i ≤ r. This theorem is related to the resolution of integer polynomial programming as follows. Thesubdivision process tests if a given region may have solutions or definitely not. If it may, thenthe region is divided and the process starts again, otherwise it is rejected. At the end, the processproduces candidates of the form Q ni = [ z i , z i +
1] and test if ( z , . . . , z n ) is a solution of the system.Both stages of this algorithm depend on the previous theorem.3 ain result In this article we present a subdivision method based on Bernstein forms to compute thesolutions of a system of congruence equations. The idea of using Bernstein forms to make a“solver” is not new, see for example [13, 22, 23] and [11, 20, 21]. One of our contribution is theapplication of Bernstein forms to congruence equations. The key ingredient is that we managedto make the subdivision process and the process of rewriting the system in Bernstein form as amatrix multiplication.An advantage of Bernstein forms is that the algorithm can test the existence of a solutionin a large region by testing bounds on the coordinates of a matrix. An important feature of thealgorithm is that it does not rely on numerical computations and as a result, it gives a certifiedanswer, see Theorem 4.2. The algorithm works by performing a matrix multiplication at each it-eration with triangular rational matrices { M Li , M Ri } ni = . These matrices depend on the multi-degreeof the system ( d , . . . , d n ) and we prove that it is convenient to have these 2 n matrices previouslycomputed for a large multi-degree ( d , . . . , d n ), see Theorem 4.4. It is worth mention that thealgorithm can work with any system of congruence equations, independently if the numbers p , . . . , p r are primes or not.The first contribution that we present is Theorem 3.3 where we give the expected numberof matrix multiplications that the algorithm for integer polynomial programming requires. Itis easily seen that this algorithm has an exponential complexity, but in practice, we prove thatthe complexity depends on the “size” of the real solution of the system, see Proposition 3.1.This “size” can be estimated using tools from algebraic geometry (dimension, degree) and inparticular, Theorem 3.3 implies that the algorithm is faster if we add more equations, that is, ifthere are less solutions to the system. We used concepts and results from the theory of Branchingprocesses to prove Theorem 3.3, Theorem.
Let λ be the complexity number (Definition 3.2) of { q , . . . , q r } over D = Q ni = [ a i , a i + k i ] and let K = k + . . . + k n . Then, the expected number of matrix multiplications produced bythe algorithm is K X i = i i − Y j = (cid:16) − (1 − λ ) K − j (cid:17) − . Regarding systems of congruence equations we first give Lemma 4.1 that characterizes ex-istence of solutions in a real interval. Then, we give the pseudocode of
SolveCE and proveTheorem 4.2
Theorem.
The output of
SolveCE ( ( p , . . . , p r ) ,v, ( k , . . . , k ) , (0 , . . . , ) is the set { x ∈ D ∩ Z n : p | h ( x ) , . . . , p r | h r ( x ) } , where D = [0 , k ] × . . . × [0 , k ] , k = ⌈ log ( lcm ( p , . . . , p r )) ⌉ and v is the matrix associated toh , . . . , h r . Finally, we compute the expected number of matrix multiplications produced by
SolveCE .This number depends on a parameter λ (similar to the complexity number),4 heorem. Let h , . . . , h r ∈ Z [ x , . . . , x n ] , p , . . . , p r ∈ Z ≥ , k = ⌈ log ( lcm ( p , . . . , p r )) ⌉ . Then,the expected number of matrix multiplications produced by SolveCE in D = Q ni = [0 , k ] is nk X i = i i − Y j = (cid:16) − (1 − λ ) nk − j (cid:17) − . This number is bounded between a polynomial and an exponential expression in λ .Summary In Section 1 we give some preliminary results; we recall the concept of Bernstein Basesand degree elevation, fixing a multi-degree ( d , . . . , d n ) we associate to every region D a squarematrix and to every system { q , . . . , q r } a rectangular matrix. Finally, with these results, weprove in Theorem 1.13 our main tool. In Section 2 we present the algorithm to solve integerpolynomial programming
SolveIPP and in
Section 3 we study its expected complexity andgive several examples computed using our implementation in SageMath [30]. Finally in
Section4 we present the pseudocode
SolveCE to solve congruence equations and study its expectedcomplexity. We compare in a table
SolveCE with the brute-force algorithm (Definition 4.5).
1. Preliminaries
In this and the next section we present a known algorithm to solve integer polynomial pro-gramming. Later we adapt it to solve congruence equations.We are interested in the integral points of semi-algebraic sets defined by a system of mul-tivariate rational polynomial equations in D = [ a , b ] × . . . × [ a n , b n ] ⊆ R n . Specifically, let q ∈ Q [ x , . . . , x n ] be some multivariate polynomial and let S ′ be the intersection between D ∩ Z n and some semi-algebraic set defined by some polynomials q , . . . , q r ∈ Q [ x , . . . , x n ]. We wantto find the maximum θ of q over S ′ and more general, to describe the set S , S = { z ∈ S ′ : q ( z ) = θ } . Notation 1.1.
First of all, we want to give a unifying description of the set S ′ . It is defined asthe intersection of equalities and inequalities. The first simplification to make is to assume thatthe polynomials have integer coe ffi cients and then, since we are working over the integers, wecan always assume that S ′ is defined only with inequalities of the form q ≥
0. Concretely, wecan make the following changes for q , q ′ ∈ Z [ x , . . . , x n ]: • Replace the inequality q ≤ q ′ with the inequality q ′ − q ≥ • Replace the equality q = q ′ with the inequalities q − q ′ ≥ q ′ − q ≥ • Replace the inequality q > q ′ with the inequality q − q ′ − ≥ • Replace the inequality q < q ′ with the inequality q ′ − q − ≥ S ′ is given as, S ′ = { z ∈ D ∩ Z n : q ( z ) ≥ , . . . , q r ( z ) ≥ } , where q , . . . , q r ∈ Z [ x , . . . , x n ] are some multivariate polynomials. Also, we can always assumethat q has integer coe ffi cients and that we are searching for the maximum of q over S ′ makingthe following changes, 5 Replace min( q ) with max( − q ). • If there is no q to maximize, then maximize q ≡ S , S = (cid:26) z ∈ S ′ : q ( z ) = max S ′ ( q ) (cid:27) . Now that we have a unifying way to present the system, let us define (or recall) Bernsteinbases. To simplify the notation, we use standard multi-index notation. Letters in boldface repre-sent multi-indexes. For example d = ( d , . . . , d n ). Definition 1.2 (Bernstein Basis) . According to [21, Proposition 2], any polynomial p ∈ R [ x , . . . , x n ]of multi-degree ( d , . . . , d n ) can be written in Bernstein form, p ( x ) = p ( x , . . . , x n ) = d X k = . . . d n X k n = β k ... k n B d ... d n , k ... k n ( x , . . . , x n ) = d X k = β k B d , k ( x ) , where β k = β k ... k n ∈ R and B d , k ( x ) = d k ! x k (1 − x ) d − k . . . d n k n ! x k n n (1 − x n ) d n − k n = dk ! x k ( − x ) d − k . The polynomials B d , k have rational coe ffi cients and multi-degree d for every k . The ordered set { B d , k : ≤ k ≤ d } is called the Bernstein basis of R [ x ] ≤ d or the d -Bernstein basis to emphasize the multi-degree.The order is the lexicographical order. The numbers β k are called Bernstein coe ffi cients . Remark 1.3.
If a polynomial is given in monomial form p ( x ) = P di = c i x i , it is easy to make achange of basis from monomial basis to Bernstein basis. Specifically, the Bernstein coe ffi cientsare computed as follows: β k = k X i = c i (cid:16) ki (cid:17)(cid:16) di (cid:17) Note that { β k } ⊆ Q if { c i } ⊆ Q . Recall from [21, Corollary 1] that if m (resp. m ) is the minimum(resp. maximum) of the Bernstein coe ffi cients β k , then we have the fundamental inequalitym ≤ p ( x ) ≤ m , ∀ x ∈ [0 , n . Lemma 1.4 (Degree elevation) . Let p ∈ Z [ x ] ≤ d be a polynomial defined over [0 , n and letm , m ∈ R be such that m < p ( x ) < m for all x ∈ [0 , n . Then, there exists d ′ such that forevery d ′′ ≥ d ′ the Bernstein coe ffi cients β k of p in d ′′ -Bernstein basis satisfy,m < β k < m , ∀ k ≤ d ′′ . roof. Given that m < p ( x ) < m , there exists δ > m + δ < p ( x ) < m − δ forall x ∈ [0 , n . Also, from [24, Eq. 3.1] there exists m ∈ N such that | β k − p ( k / m ) | ≤ δ for all k ≤ d ′ : = ( m , . . . , m ). Then, β k = ( β k − p ( k / m )) + p ( k / m ) < δ + ( m − δ ) = m and β k = ( β k − p ( k / m )) + p ( k / m ) > − δ + ( m + δ ) = m . Example 1.5.
Consider the following polynomial over [0 , p ( x ) = − x − ! − . Then max [0 , ( p ) = − /
10. The Bernstein coe ffi cients in degree 2 of p are ( − / , / , − / ffi cients in degree d ′ = − / , − / , − / , − / d ′′ ≥ − ≤ max k ≤ d ′′ ( β k ) ≤ − . Let us start analyzing the regions D . Our goal is to define subdivision matrices for [0 , n andthen, to treat any D in a uniform way. Definition 1.6 (Matrix associated to a box and to a multi-degree) . We say that D ⊆ R n is a box if there exist a , b ∈ R n such that D = [ a , b ] × . . . × [ a n , b n ] . Let ϕ : R n → R n be ϕ ( x ) = ( b − a ) x + a . Clearly ϕ maps bijectively [0 , n to D and defines, viapull-back, a linear map R [ x ] → R [ x ] given by q ϕ ∗ ( q ) : = q (( b − a ) x + a ).Fix a multi-degree d . We say that the matrix M of the map ϕ ∗ in the d -Bernstein basis is the matrix associated to the box D and to the multi-degree d . The matrix M is a s × s -square matrix,where s = ( d + . . . ( d n +
1) and also, if { a i , b i } ni = ⊆ Q , then M ∈ Q s × s . Notation 1.7.
Fix a multi-degree d . Let us denote M Li or M Li ( d ) (resp. M Ri or M Ri ( d )) to thematrix associated to the box [0 , i − × [0 , / × [0 , n − i (resp. [0 , i − × [1 / , × [0 , n − i )for 1 ≤ i ≤ n and to d (see Definition 1.6). Recall that if s = ( d + . . . ( d n + M Li , M Ri ∈ Q s × s . Proposition 1.8.
The matrix M Li ( d ) is lower triangular and M Ri ( d ) is upper triangular, ≤ i ≤ n.If i , j, M Li commutes with M Lj and with M Rj , but M Li M Ri , M Ri M Li .Proof. The shape of the matrices M Li and M Ri follows from [21, Proposition 5]. Let us prove thecommutativity.Let ϕ Li (resp. ϕ Ri ) be the unique a ffi ne isomorphism sending [0 , n to [0 , i − × [0 , / × [0 , n − i (resp. to [0 , i − × [1 / , × [0 , n − i ). The matrix representation of ϕ Li in d -Bernsteinbasis is M Li (same for ϕ Ri and M Ri ). Then, we need to prove that if i , j , then ϕ Li commutes7ith ϕ Lj and with ϕ Rj , but ϕ Li ϕ Ri , ϕ Ri ϕ Li . Clearly ϕ Li ϕ Lj and ϕ Lj ϕ Li (resp. ϕ Li ϕ Rj and ϕ Rj ϕ Li ) are a ffi neisomorphisms sending [0 , n to the same box. By uniqueness ϕ Li ϕ Lj = ϕ Lj ϕ Li (resp. ϕ Li ϕ Rj = ϕ Rj ϕ Li ).Given that the box associated to ϕ Li ϕ Ri is di ff erent from the box associated to ϕ Ri ϕ Li , we obtain ϕ Li ϕ Ri , ϕ Ri ϕ Li ,[0 , i − × [1 / , / × [0 , n − i , [0 , i − × [1 / , / × [0 , n − i . Proposition 1.9.
Fix a multi-degree d . Let k i , l i ∈ N be such that l i < k i , ≤ i ≤ n. Let D bethe box defined as D = " l k , l + k × . . . × " l n k n , l n + k n and let M be the matrix associated to D and d . Then, there exists a factorization of M into aproduct of matrices M Li ( d ) and M Ri ( d ) . The sum of the multiplicities of the matrices M Li ( d ) andM Ri ( d ) in this product is k i , ≤ i ≤ n.Proof. Without loss of generality, we may assume that n = D = [ l / k , ( l + / k ]. Let uscall M L = L and M R = R .Induction in k . If k =
1, then l = l =
1. Hence, M = L or M = R . Assume now the resultfor k − k . We have two possibilities, l is even or l + l = l ′ is even, then D is included in D ′ = [ l / k , ( l + / k ] = [ l ′ / k − , ( l ′ + / k − ]. By the inductivehypothesis, D ′ has associated a matrix M ′ that can be factorized as k − L and R , M ′ = L i R i . . . . Given that D is equal to the first half of D ′ , M = LM ′ . The case inwhich l + Definition 1.10.
Let F ( k ) be the set of functions from { , . . . , k } to { L , R } . The cardinal of F ( k )is 2 k .Fix a multi-degree d and k ∈ N n . The set of ( d , k ) -subdivision matrices is defined as n M f (1)1 . . . M f ( k )1 M f (1)2 . . . M f ( k )2 . . . M f n (1) n . . . M f n ( k n ) n : ( f , . . . , f n ) ∈ F ( k ) × . . . × F ( k n ) o . If k =
0, we use the convention F ( k ) = { } and M i is the identity matrix. According to Proposi-tion 1.9, the cardinal of this set is equal to the number of boxes obtained by subdividing k i timesin direction x i the box [0 , n . Hence, the cardinal is equal to 2 k + ... + k n .Now that we know how to treat di ff erent boxes D as matrices, let us give some results on howto treat polynomials as vectors. Proposition 1.11.
Let d ∈ N n , s = ( d + . . . ( d n + , D = [ a , b ] × . . . × [ a n , b n ] andlet q ∈ R [ x ] ≤ d be a multivariate polynomial of multi-degree ≤ d . Then, there exist a vectorv = v ( q ) ∈ R s such that, • The first coe ffi cient v is equal to q ( a , . . . , a n ) ,v = q ( a , . . . , a n ) . • The values of q over D are bounded by min( v ) and max( v ) , min( { v i } si = ) ≤ q ( x ) ≤ max( { v i } si = ) , ∀ x ∈ D . roof. Let ϕ : R n → R n be the unique a ffi ne isomorphism sending [0 , n to D , ϕ ( x ) = ( b − a ) x + a . Then, the pull-back of q over D , defines a polynomial ϕ ∗ ( q ) = q (( b − a ) x + a ) over [0 , n . Let v be the vector representation of ϕ ∗ ( q ) in d -Bernstein basis. Then,min( v ) ≤ q ( x ) ≤ max( v ) , ∀ x ∈ D and from [21, Proposition 4], v = b = ϕ ∗ ( q )( ) = q ( ϕ ( )) = q ( a , . . . , a n ) . Corollary 1.12.
Let d ∈ N n , s = ( d + . . . ( d n + , D = [ a , b ] × . . . × [ a n , b n ] and letq , . . . , q r ∈ R [ x ] ≤ d be polynomials of multi-degree ≤ d . Then, there exist a rectangular matrixv ∈ R s × r with columns v i = v i ( q i ) such that, • The coe ffi cient v i is equal to q i ( a , . . . , a n ) ,v i = q i ( a , . . . , a n ) , ≤ i ≤ r . • The values of q i over D are bounded by the minimum and the maximum of { v i , . . . , v si } , min( { v ji } sj = ) ≤ q i ( x ) ≤ max( { v ji } sj = ) , ∀ x ∈ D , ≤ i ≤ r . Proof.
Define the matrix v = ( v , . . . , v r ) ∈ R s × r , where v i = v i ( q i ) is defined as in Proposition1.11.Combining the previous results on boxes and polynomials, we obtain the following importantresult, Theorem 1.13.
Let d = ( d , . . . , d n ) be a multi-degree, s = ( d + . . . ( d n + and let D = [ a , a + k ] × . . . × [ a n , a n + k n ] be a box, where k = ( k , . . . , k n ) ∈ N n is fixed.For every q , . . . , q r ∈ Z [ x ] ≤ d there exists a rectangular matrix v ∈ Q s × r and for every boxD ′ ⊆ D of the form D ′ = [ z , z + k ′ ] × . . . × [ z n , z n + k ′ n ] , z i = a i + l i k ′ i , ≤ k ′ i ≤ k i , ≤ l i < k i − k ′ i , l i ∈ N , ≤ i ≤ n , there exists a ( d , k ) -subdivision matrix M ∈ Q s × s such that w = Mv satisfies • The coe ffi cient w i is equal to q i ( z , . . . , z n ) ,w i = q i ( z , . . . , z n ) , ≤ i ≤ r . • The values of q i over D ′ are bounded by the minimum and the maximum of { w i , . . . , w si } , min( { w ji } sj = ) ≤ q i ( x ) ≤ max( { w ji } sj = ) , ≤ i ≤ r , x ∈ D ′ . roof. Let ϕ ( x , . . . , x n ) = ( a + k x , . . . , a n + k n x n ) and consider the system ϕ ∗ ( q ) , . . . , ϕ ∗ ( q r )over [0 , n . Let v ∈ Q s × r be the rectangular matrix such that its i -th column is the vectorrepresenting ϕ ∗ ( q i ) in the Bernstein basis, 1 ≤ i ≤ r .The idea now is to compare the system { q i } ri = over D and the system { ϕ ∗ ( q i ) } ri = over [0 , n .First of all, it is easy to check that ϕ maps bijectively the box A ⊆ [0 , n to D ′ ⊆ D , A = n Y i = " z i − a i k i , z i − a i + k ′ i k i , D ′ = n Y i = h z i , z i + k ′ i i . Let φ A be the a ffi ne isomorphism sending [0 , n to A and let φ D ′ the one sending [0 , n to D ′ ,that is, φ A ( x , . . . , x n ) = (( z − a + x k ′ ) / k , . . . , ( z n − a n + x n k ′ n ) / k n ) and φ D ′ ( x , . . . , x n ) = ( z + x k ′ , . . . , z n + x n k ′ n ). Clearly φ D ′ = ϕφ A , hence the Bernstein coe ffi cients of φ ∗ D ′ ( q i ) arethe same as the Bernstein coe ffi cients of φ ∗ A ( ϕ ∗ ( q i ) ). Then, the values of { q i } ri = over D ′ are thesame as the values of { ϕ ∗ ( q i ) } ri = over A . The result follows from Proposition 1.9 by noting that A = Q ni = [ l i / k i − k ′ i , ( l i + / k i − k ′ i ] and letting M be the matrix associated to A and d . Definition 1.14.
Boxes appearing in the previous theorem can be parameterized by two multi-indexes ( l , k ′ ). Let a ∈ Z n , k ∈ N n and let D = Q ni = [ a i , a i + k i ]. For every ( l , k ′ ) such that ≤ k ′ ≤ k and 0 ≤ l i < k i − k ′ i , l i ∈ N , 1 ≤ i ≤ n define D l , k ′ as D l , k ′ = n Y i = [ a i + l i k ′ i , a i + ( l i + k ′ i ] . Otherwise, define D l , k ′ = ∅ . The vectors a and k are omitted from the notation. Notice that D , k = D and that the volume of D l , k ′ is 2 k ′ + ... + k ′ n (or 0). The index l represents the position andthe index k ′ the volume of the box D l , k ′ .
2. Subdivision algorithm
Let d ∈ N n be a multi-degree, s = ( d + . . . ( d n +
1) and let { M Li ( d ) , M Ri ( d ) } ni = ⊆ Q s × s .Given that these matrices depend only on d , it is possible to have them previously computed.For example, depending on the application, we can compute these matrices to solve systems in Z [ x ] ≤ d for some large d >> .Let a ∈ Z n , k ∈ N n , D = [ a , a + k ] × . . . × [ a n , a n + k n ] and let q , . . . , q r ∈ Z [ x ] ≤ d . Wewant to describe the following sets, S ′ = { z ∈ D ∩ Z n : q ( z ) ≥ , . . . , q r ( z ) ≥ } , S = (cid:26) z ∈ S ′ : q ( z ) = max S ′ ( q ) (cid:27) . Let v ∈ Q s × r denote the rectangular matrix associated to q , . . . , q r (see Theorem 1.13). Initialize10he global variables S : = ∅ and θ : = −∞ and run the function SolveIPP ( v , , k ). Function
SolveIPP( w , l , k ′ ) Data:
Global variable θ , a global set S and matrices { M Li ( d ) , M Ri ( d ) } ni = . if max( w ) ≥ θ and max( w ) ≥ and . . . and max( w r ) ≥ then k ′ j : = max( k ′ ) if k ′ j > then SolveIPP ( M Lj w , l + l j e j , k ′ − e j ) SolveIPP ( M Rj w , l + l j e j + e j , k ′ − e j ) else if w ≥ θ and w ≥ and . . . and w r ≥ then if w > θ then S : = {} θ : = w S : = S ∪ { a + l } Let us explain the subdivision process of
SolveIPP (Solve Integer Polynomial Program-ming). The first thing to say is that it is a recursive algorithm that makes two matrix multipli-cations in each call, M Lj v and M Rj v . The rationality of the matrices makes this algorithm robustand reliable and the resulting set S is equal to its mathematical counterpart S (see Theorem 2.1below).According to Theorem 1.13 the values in ( w , l , k ′ ) give information about the polynomials { q , . . . , q r } over the box D l , k ′ . The condition in Step 1 is given to test if there exist solutions on D l , k ′ , in other words, if we keep D l , k ′ or we reject it. Assuming that we do not reject D l , k ′ , wetest in Step 3 if D l , k ′ can be divided.If D l , k ′ can be divided, we define D L = D l + l j e j , k ′ − e j and D R = D l + l j e j + e j , k ′ − e j and we analyzethe values of the system over these boxes by defining w L = M Lj v and w R = M Rj v . The notation e j indicates the j -th cannonical vector,( e j ) i = i = j i , j , ≤ i ≤ n . If D l , k ′ can not be divided, that is k ′ = , we test in Step 6 if a + l is actually a solution of thesystem. If this is the case, we save z = a + l in S . If w > θ , we reset S and set θ = q ( z ). Thislast process guarantees that q ( z ) = θ for every z ∈ S . Theorem 2.1.
Let n ∈ N , d , k ∈ N n , a ∈ Z n , q , . . . , q r ∈ Z [ x , . . . , x n ] ≤ d , let D = Q ni = [ a i , a i + k i ] be a box and let v be the matrix associated to q , . . . , q r as in Theorem 1.13. Then, the output ofthe function SolveIPP (v, , k ) is the set S , S ′ = { z ∈ D ∩ Z n : q ( z ) ≥ , . . . , q r ( z ) ≥ } , S = (cid:26) z ∈ S ′ : q ( z ) = max S ′ ( q ) (cid:27) . Proof.
The compactness of D implies that the set S ′ is finite. If S ′ , ∅ , the set q ( S ′ ) has amaximum value θ and the space of solutions S is equal to q − ( θ ) ∩ S ′ .Let us prove that every z ∈ S is in the output of SolveIPP (the reciprocal is straightforward).Let D ′ be a box such that z ∈ D ′ ∩ S . Then max D ′ ( q ) ≥ θ and max D ′ ( q i ) ≥
0, 2 ≤ i ≤ r . Inparticular, the box D ′ is not rejected by the algorithm and is subdivided. Then, without loss ofgenerality, we may suppose that D ′ = Q ni = [ z i , z i +
1] and clearly, z is in the output of SolveIPP .Note that after some steps, the value of the global variable θ is equal to the maximum of q over S ′ . 11 . Expected complexity Given that the number of matrix multiplications can be computed as the number of times thefunction
SolveIPP is called (minus one), the complexity of the algorithm can be computed interms of boxes. The worst case scenario is when every box is subdivided (exponential case). Thebest case occurs when almost every box is unnecessary (linear case). Thus, the complexity of thealgorithm is related to the number of unnecessary boxes, and it can be described in probabilisticterms. The following proposition gives a geometric representation of the set of unnecessaryboxes.
Proposition 3.1.
Let n ∈ N , d , k ∈ N n , a ∈ Z n , q , . . . , q r ∈ Z [ x , . . . , x n ] ≤ d , D = Q ni = [ a i , a i + k i ] and let θ be the maximum of q over S ′ = { q i ≥ } ri = ∩ D ∩ Z n . Let R be the set of real solutionsof the system, R = { z ∈ D : q ( z ) ≥ θ, q ( z ) ≥ , . . . , q r ( z ) ≥ } . Then, there exists d ′ ∈ N n such that for every d ′′ ≥ d ′ the following sentences are equivalent forevery non-empty box D l , k ′ : The matrix associated to the system over D l , k ′ satisfies the condition in Step 1. D l , k ′ ∩ R , ∅ . There exists D l ′ , ⊆ D l , k ′ such that D l ′ , ∩ R , ∅ .Proof. Let us prove the existence of d ′ . First, replace q by q − θ , and write condition in Step 1as { max( w i ) ≥ } ri = . In each box D l , k ′ take some d l , k ′ satisfyingmax D l , k ′ ( q i ) < = ⇒ max( w i ) < , ≤ i ≤ r , where the rectangular matrix w = ( w , w , . . . , w r ) is the one associated to the system over thebox D l , k ′ and multi-degree d l , k ′ , see Lemma 1.4. Let d ′ be the maximum of { d l , k ′ } .Assume that we have a box D l , k ′ without real solutions in it. Let us see that the algorithmapplied in multi-degree d ′ rejects this box immediately. Given that D l , k ′ does not contain realsolutions, there exists some i , 1 ≤ i ≤ r , such that min( q i ) <
0. Then, min( w i ) < D l , k ′ is rejected.The equivalence between the last two sentences follows from the fact that it is possible tocover D l , k ′ with boxes of the form D l ′ , . Definition 3.2.
Let n ∈ N , d , k ∈ N n , a ∈ Z n , q , . . . , q r ∈ Z [ x , . . . , x n ] ≤ d and D = Q ni = [ a i , a i + k i ], let θ be the maximum of q over S ′ = { q i ≥ } ri = ∩ D ∩ Z n and let R be the set of realsolutions, R = { z ∈ D : q ( z ) ≥ θ, q ( z ) ≥ , . . . , q r ( z ) ≥ } .The complexity number of the system { q , . . . , q r } over D is defined as λ ( R ) = { l : D l , ∩ R , ∅} k + ... + k n , λ ( R ) ∈ Q , ≤ λ ( R ) ≤ . The complexity number can be described in probabilistic terms, it is the probability of a box D l , to have real solutions. Theorem 3.3.
Let λ be the complexity number of { q , . . . , q r } over D = Q ni = [ a i , a i + k i ] and letK = k + . . . + k n . Let d ′ be the multi-degree from Proposition 3.1. Then, the expected number of ultiplications between a square triangular matrix in Q s ′ × s ′ and a rectangular matrix in Q s ′ × r ,s ′ = ( d ′ + . . . ( d ′ n + produced by SolveIPP is K X i = i i − Y j = (cid:16) − (1 − λ ) K − j (cid:17) − . This number is bounded between ((2 λ ) K + − / (2 λ − − (or K if λ = / ) and K − .Proof. Each time a box is subdivided, the function
SolveIPP produces two matrix multiplica-tions. Hence, the expected number of matrix multiplications is equal to the expected number ofdivided boxes. Equivalently, it is the expected number of boxes processed by
SolveIPP minus 1( D is not a divided box). According to Proposition 3.1, the probability of a box D l , k ′ to be subdi-vided depends on the existence of some D l ′ , such that D l ′ , ⊆ D l , k ′ and D l ′ , ∩ R ,
0. There are atotal of 2 k ′ + ... + k ′ n possible boxes D l ′ , inside D l , k ′ hence the probability is equal to 1 − (1 − λ ) k ′ + ... + k ′ n .We can model the problem of finding the expected number of boxes as a Branching process , see[15]. We say that D l , k ′ is in generation i if k ′ + . . . + k ′ n = K − i . In other words, if the volume of D l , k ′ is equal to 2 K − i . Let { Z i } Ki = be the number of boxes in each generation. Boxes in generation i have probability λ i (resp. 1 − λ i ) to add two (resp. 0) boxes to the next generation, λ i = − (1 − λ ) K − i , ≤ i ≤ K . These boxes can be interpreted as independent identically distributed random variables { D i } withcommon generating function 1 − λ i + λ i x .Thus Z i + = D + . . . + D Z i and from [14, § E ( Z i + ) = λ i E ( Z i ), hence E ( Z ) = , E ( Z i + ) = λ i E ( Z i ) = i + λ i . . . λ , ≤ i ≤ K − . The number of expected boxes is E ( Z ) + . . . + E ( Z K ). Example.
Consider the equation y = x in [0 , × [0 , , × [0 ,
1] and [1 , × [1 , λ = / SolveIPP for d ′ as in Proposition 3.1 is 91 /
16, a number between 2 + = − =
7. The function applied to this example passes Step 1, 7 times.
Example.
Consider the system y − x = D = [0 , ] × [0 , ]. The system has three integralsolutions, (0 , ,
1) and (2 , D l , . Then, λ = / , ( λ , . . . , λ ) (cid:27) ( . , . , . , . , . , . , . . Hence, the expected number of boxes processed by
SolveIPP for d ′ as in Proposition 3.1 isapproximately 34. In this example Step 1 occurs 33 times. Remark.
In [23] the authors manage to to improve (a variant of)
SolveIPP to produce betterperformance. Also, it is possible to adapt
SolveIPP to treat systems of continuous functions byusing Stone-Weierstrass Theorem [29, Cor.1] and interpolation techniques [7, 9].Notice that Proposition 3.1 implies that the strategy of dividing a side of a box in half is notoptimal. For example, consider the equation x − , SolveIPP divides[0 ,
4] in two halves making two matrix multiplications, but given that
SolveIPP do not rejectthese two halves, it makes four multiplications more.Instead of dividing [0 ,
4] into [0 ,
2] and [2 ,
4] it is better (from a performance point of view)to divide it into [0 ,
2) and [2 ,
4] or even better into [0 ,
1] and [2 , emark. Given d ∈ N n , we need to compute the triangular matrices { M Li , M Ri } ni = ⊆ Q s × s , s = Q ni = ( d i +
1) and save them into the computer memory. If ns ( s +
1) is greater than thetotal amount of available memory it is possible to use techniques from computer science such asParallel computing or Map-Reduce [10, 25] in order to increase the available space of memory.An alternative (or complementary) solution is to use toric
Bernstein polynomials to minimize thenumber of monomials, [31, Cor.2].
4. Solving systems of congruence equations
In this section we apply the subdivision method presented in Section 2 to solve congruenceequations (
SolveCE ). Let h , . . . , h r ∈ Z [ x , . . . , x n ] ≤ d and let p , . . . , p r ∈ Z ≥ be natural num-bers (possibly coprimes). Here we propose a method to find solutions z ∈ Z n to the system h ( z ) ≡ , mod ( p ) ,... h r ( z ) ≡ , mod ( p r ) . In order to do so, we need to adapt some of the propositions presented so far. Notice that anysolution z ∈ Z n can be represented in the box [0 , k ] × . . . × [0 , k ], where k : = ⌈ log (lcm( p , . . . , p r )) ⌉ ∈ Z ≥ , lcm( x , y ) is the least common multiple of x and y and ⌈ x ⌉ is the least integer greater than or equalto x . The following simple Lemma is the key ingredient of our method. Lemma 4.1.
Let p ∈ Z ≥ , let m ≤ m be two integers and let r (resp. r ) be the remainder ofm (resp. m ) divided by p. The following sentences are equivalent, m − m < p − and < r ≤ r . There is no multiple of p between m and m .Proof. Let q , q ∈ Z be such that m = q p + r and m = q p + r . Given that m ≤ m < m + p −
1, we have the inequality, p − > m − m = ( q − q ) p + ( r − r ) ≥ . ⇒
2) Given that 0 ≤ r − r ≤ p −
1, we have ( q − q ) p ≥ p − > ( q − q ) p ≥ q = q . Then, there exists q ∈ Z such that m = qp + r and m = qp + r . Assumenow, that there exists a ∈ Z such that m ≤ ap ≤ m . Then, m − ap is also < p − a = q . Hence, m ≤ qp ≤ m ⇐⇒ qp + r ≤ qp ≤ qp + r ⇐⇒ r ≤ ≤ r . A contradiction, because r > ⇒
1) It is clear that if m − m ≥ p −
1, then there exists a multiple of p between m and m .Same if r = r =
0. Assume m − m < p − r > r and let us prove that there existsa multiple of p between m and m .Using the inequality p − > ( q − q ) p + ( r − r ) ≥ < r − r ≤ p − q = q +
1. Then, q p is in between m and m .14 direct application of Lemma 4.1 gives the following function, Function hasSol( m , m , p ) Data: m ≤ m ∈ Q and p ∈ Z ≥ . Output:
True if there is a multiple of p between m and m . False otherwise. ( m , m ) ← ( ⌈ m ⌉ , ⌊ m ⌋ ) ( r , r ) ← ( m % p , m % p ) if m − m < p − and < r ≤ r then return False return True
Step 1 in hasSol is based on the following fact. There exists an integer i such that m ≤ i ≤ m if only if ⌈ m ⌉ ≤ i ≤ ⌊ m ⌋ , where ⌈ m ⌉ is the least integer greater than or equal to m and ⌊ m ⌋ isthe greatest integer less than or equal to m . The notation m % p in Step 2 means the remainder of dividing m by p . Same for m % p .Now, we can adapt the function SolveIPP . Fix some d ∈ N n and let { M Li ( d ) , M Ri ( d ) } ni = bethe matrices defined in 1.7. Let h , . . . , h r ∈ Z [ x ] ≤ d and let p , . . . , p r ∈ Z ≥ be natural numbers.Let k : = ⌈ log (lcm( p , . . . , p r )) ⌉ ∈ Z ≥ and D = [0 , k ] × . . . × [0 , k ]. Denote p = ( p , . . . , p r ), = (0 , . . . ,
0) and k = ( k , . . . , k ).We want to describe the following set, S ′ = { z ∈ D ∩ Z n : h ( z ) ≡ p ) , . . . , h r ( z ) ≡ p r ) } . Let v be the rectangular matrix associated to h , . . . , h r as in Theorem 1.13. Initialize the globalset S ′ : = ∅ and run SolveCE ( p , v , k , ). Function
SolveCE( p , w , k ′ , l ) Data:
A global set S ′ and matrices { M Li ( d ) , M Ri ( d ) } ni = . Output: S ′ = { z ∈ D ∩ Z n : h ( z ) ≡ p ) , . . . , h r ( z ) ≡ p r ) } if hasSol ( min( w ) , max( w ) , p ) and . . . and hasSol ( min( w r ) , max( w r ) , p r ) then k ′ j : = max( k ′ ) if k ′ j > then SolveCE ( p , M Lj w , k ′ − e j , l + l j e j ) SolveCE ( p , M Rj w , k ′ − e j , l + l j e j + e j ) else if p | w and . . . and p r | w r then S ′ : = S ′ ∪ { l } Theorem 4.2.
The output of
SolveCE ( p ,v, ( k , . . . , k ) , (0 , . . . , ) is the set S ′ = { z ∈ D ∩ Z n : p | h ( z ) , . . . , p r | h r ( z ) } , where D = [0 , k ] × . . . × [0 , k ] , k = ⌈ log ( lcm ( p , . . . , p r )) ⌉ and v is the matrix associated toh , . . . , h r as in Theorem 1.13.Proof. If there is a solution in a box D ′ , then the process goes from Step 1 to Step 2 and D ′ getssubdivided unless the volume of D ′ is equal to 1. If the volume of D ′ is equal to 1 and thereexists a solution l ∈ D ′ , then the condition in Step 6 is satisfied and l is saved in S ′ .15s before, we can estimate the expected number of matrix multiplications that SolveCE do.Let h : R n → R r be the polynomial function h = ( h , . . . , h r ) and consider the discrete set X = { ( x , . . . , x r ) ∈ Z r : p | x , . . . , p r | x r } . Notice that the set h − ( x ) ⊆ R n is a real algebraic variety for each x ∈ X . Proposition 4.3.
Let n ∈ N , d ∈ N n , h , . . . , h r ∈ Z [ x , . . . , x n ] ≤ d , p , . . . , p r ∈ Z ≥ and letk = ⌈ log ( lcm ( p , . . . , p r )) ⌉ . Let X = { x ∈ Z r : p | x , . . . , p r | x r } and let h : R n → R r , h = ( h , . . . , h r ) . Then, there exists d ′ ∈ N n such that for every d ′′ ≥ d ′ and every non-empty boxD l , k ′ , the following sentences are equivalent: The matrix associated to ( h , . . . , h r ) over D l , k ′ satisfies the condition in Step 1. D l , k ′ ∩ h − ( X ) , ∅ . There exists D l ′ , ⊆ D l , k ′ such that D l ′ , ∩ h − ( X ) , ∅ .Proof. Clearly sentence 2 is equivalent to sentence 3. Let us prove the equivalence betweensentences 1 and 2. Fix some l and k ′ . Let us define m i and m i , 1 ≤ i ≤ r , m i : = ⌊ min D l , k ′ ( h i ) ⌋ if min D l , k ′ ( h i ) > ⌊ min D l , k ′ ( h i ) ⌋ min D l , k ′ ( h i ) − if not m i : = ⌈ max D l , k ′ ( h i ) ⌉ if max D l , k ′ ( h i ) < ⌈ max D l , k ′ ( h i ) ⌉ max D l , k ′ ( h i ) + if notChoose a multi-degree d l , k ′ as in Lemma 1.4 such that m i < min( w i ) ≤ max( w i ) < m i , ≤ i ≤ r , where w = ( w , . . . , w r ) is the rectangular matrix associated to ( h , . . . , h r ) over D l , k ′ . By con-struction, there exists a multiple of p i between min D l , k ′ ( h i ) and max D l , k ′ ( h i ) if and only if it isbetween min( w i ) and max( w i ). The result follows by taking d ′ as the maximum of { d l , k ′ } . Theorem 4.4.
Let h , . . . , h r ∈ Z [ x , . . . , x n ] , p , . . . , p r ∈ Z ≥ , k = ⌈ log ( lcm ( p , . . . , p r )) ⌉ andlet d ′ be as in Proposition 4.3. Let X = { x ∈ Z r : p | x , . . . , p r | x r } and let h : R n → R r , h = ( h , . . . , h r ) . Then, the expected number of matrix multiplications produced by SolveCE inD = Q ni = [0 , k ] is nk X i = i i − Y j = (cid:16) − (1 − λ ) nk − j (cid:17) − , λ : = { l : D l , ∩ h − ( X ) , ∅} nk . The number λ is the probability of a box D l , to intersect h − ( X ) .Proof. According to Proposition 4.3 the probability of a box D l , k ′ to be rejected in Step 1 is(1 − λ ) k ′ + ... + k ′ n . The proof continues as in Theorem 3.3. Definition 4.5.
We define the brute-force algorithm as the algorithm consisting in testing each z ∈ D ∩ Z n to be in S ′ , S ′ = { z ∈ D ∩ Z n : p | h ( z ) , . . . , p r | h r ( z ) } . r nk evaluations to characterize S ′ . Notice that there is no need for thesides of D to be a power of two and the same is true in SolveCE . We chose to have sides to be apower of two to make the presentation clearer and also to be able to compare this algorithm with
SolveCE .For example, if p = p = h = x + h = y we need to take the starting box to be D = [0 , . The brute-force algorithm with the box D = [0 , performs 25 × , SolveCE finds two (apparently di ff erent) solutions (4 , ,
5) in the box D = [0 , .In order to avoid benefiting either algorithm, we choose powers of two for the numbers p , . . . , p r . In the next table we compare SolveCE with the brute-force algorithm. All computa-tions were made using our implementation in SageMath. p h
Step 1 Step 6 BF512 2 x − x + x − x + x − , x + , y ) 40 4 262144(8 ,
8) ( xy + x + y , y + x ) 123 60 64(128 , y ( x − , y +
1) 29124 13675 16384(16 , ,
16) ( x + y − z , x + y − z , z − + x ) 491 84 4096(4 , ,
16) ( x + yz − z , y , x + z +
2) 1923 720 512The first two columns codify the system h ≡ mod ( p ). The column Step 1 is the number oftimes
SolveCE pass through the first step. The column
Step 6 is the number of boxes D l , testedto have a solution. Finally, the last column is the number of evaluations that the brute-forcealgorithm must do. We can infer from the table that SolveCE is faster for linear systems, but forlinear systems there are much better algorithms, [4, 27]. We have prioritized presentation overperformance.
References [1] Alexander Barvinok and Kevin Woods. Short rational generating functions for lattice point problems.
J. Amer.Math. Soc. , 16(4):957–979 (electronic), 2003.[2] Fani Boukouvala, Ruth Misener, and Christodoulos A. Floudas. Global optimization advances in Mixed-IntegerNonlinear Programming, MINLP, and Constrained Derivative-Free Optimization, CDFO.
European J. Oper. Res. ,252(3):701–727, 2016.[3] Michael R Bussieck and Stefan Vigerske. Minlp solver software.
Wiley encyclopedia of operations research andmanagement science , 2010.[4] Henri Cohen.
Number theory. Vol. I. Tools and Diophantine equations , volume 239 of
Graduate Texts in Mathe-matics . Springer, New York, 2007.[5] Henri Cohen.
Number theory. Vol. II. Analytic and modern tools , volume 240 of
Graduate Texts in Mathematics .Springer, New York, 2007.[6] William Cook, Thorsten Koch, Daniel E. Ste ff y, and Kati Wolter. A hybrid branch-and-bound approach for exactrational mixed-integer programming. Math. Program. Comput. , 5(3):305–344, 2013.[7] Wolfgang Dahmen. Subdivision algorithms converge quadratically.
J. Comput. Appl. Math. , 16(2):145–158, 1986.[8] Claudia D’Ambrosio and Andrea Lodi. Mixed integer nonlinear programming tools: a practical overview. ,9(4):329–349, 2011.[9] Carl de Boor. B -form basics. In Geometric modeling , pages 131–148. SIAM, Philadelphia, PA, 1987.[10] Je ff rey Dean and Sanjay Ghemawat. Mapreduce: Simplified data processing on large clusters. Commun. ACM ,51(1):107–113, January 2008.
11] Gershon Elber and Myung-Soo Kim. Geometric constraint solver using multivariate rational spline functions. In
Proceedings of the Sixth ACM Symposium on Solid Modeling and Applications , SMA ’01, pages 1–10, New York,NY, USA, 2001. ACM.[12] R. Fourer. Linear programming: Software survey. OR / MS TODAY , 42(3), 2015.[13] J¨urgen Garlo ff . The bernstein algorithm. Interval computations (now Reliable Computing) , 2(6):154–168, 1993.[14] Geo ff rey R. Grimmett and David R. Stirzaker. Probability and random processes . Oxford University Press, NewYork, third edition, 2001.[15] Theodore E. Harris.
The theory of branching processes . Die Grundlehren der Mathematischen Wissenschaften,Bd. 119. Springer-Verlag, Berlin; Prentice-Hall, Inc., Englewood Cli ff s, N.J., 1963.[16] R. C. Jeroslow. There cannot be any algorithm for integer programming with quadratic constraints. OperationsResearch , 21(1):221–224, 1973.[17] Sven Ley ff er (eds.) Jon Lee. Mixed Integer Nonlinear Programming . The IMA Volumes in Mathematics and itsApplications 154. Springer-Verlag New York, 1 edition, 2012.[18] Michael J¨unger, Thomas Liebling, Denis Naddef, George Nemhauser, William Pulleyblank, Gerhard Reinelt, Gio-vanni Rinaldi, and Laurence Wolsey, editors.
50 years of integer programming 1958–2008 . Springer-Verlag, Berlin,2010. From the early years to the state-of-the-art, Papers from the 12th Combinatorial Optimization Workshop(AUSSOIS 2008) held in Aussois, January 7–11, 2008.[19] Yuri V. Matiyasevich.
Hilbert’s tenth problem . Foundations of Computing Series. MIT Press, Cambridge, MA,1993. Translated from the 1993 Russian original by the author, With a foreword by Martin Davis.[20] B. Mourrain and J. P. Pavone. Subdivision methods for solving polynomial equations.
J. Symbolic Comput. ,44(3):292–306, 2009.[21] C´esar Mu˜noz and Anthony Narkawicz. Formalization of Bernstein polynomials and applications to global opti-mization.
J. Automat. Reason. , 51(2):151–196, 2013.[22] P SV Nataraj and M Arounassalame. A new subdivision algorithm for the bernstein polynomial approach to globaloptimization.
International journal of automation and computing , 4(4):342–352, 2007.[23] Bhagyesh V Patil and PSV Nataraj. An improved bernstein global optimization algorithm for minlp problems withapplication in process industry.
Mathematics in Computer Science , 8(3-4):357–377, 2014.[24] Hartmut Prautzsch and Leif Kobbelt. Convergence of subdivision and degree elevation.
Adv. Comput. Math. ,2(1):143–154, 1994.[25] Anand Rajaraman and Je ff rey David Ullman. Mining of Massive Datasets . Cambridge University Press, New York,NY, USA, 2011.[26] Fabio Schoen. Stochastic techniques for global optimization: a survey of recent advances.
J. Global Optim. ,1(3):207–228, 1991.[27] Nigel P. Smart.
The algorithmic resolution of Diophantine equations , volume 41 of
London Mathematical SocietyStudent Texts . Cambridge University Press, Cambridge, 1998.[28] E. M. B. Smith and C. C. Pantelides. Symbolic methods for global optimisation. In
Symbolic methods in controlsystem analysis and design , volume 56 of
IEE Control Eng. Ser. , pages 321–338. IEE, London, 1999.[29] M. H. Stone. The generalized Weierstrass approximation theorem.
Math. Mag. , 21:167–184, 237–254, 1948.[30] The Sage Developers.
SageMath, the Sage Mathematics Software System (Version 7.0) , 2016. .[31] Steve Zelditch. Bernstein polynomials, Bergman kernels and toric K¨ahler varieties.
J. Symplectic Geom. , 7(2):51–76, 2009., 7(2):51–76, 2009.