Minimal solutions of linear diophantine systems : bounds and algorithms

Loïc Pottier

January 9, 1991

S.A.F.I.R. Project Institut National de Recherche en Informatique et Automatique 2004 route des Lucioles, Sophia Antipolis, 06565 Valbonne CEDEX, FRANCE

                         Email : pottier@mirsa.inria.fr

Abstract : We give new bounds and algorithms for minimal solutions of linear diophantine systems. These bounds are simply exponential, while previous known bounds were, at least until recently, doubly exponential.

1 Introduction

A linear diophantine system Ax ≤ b is a set of inequations with integer coefficients (A is a matrix of integers with m rows and n columns, b a vector of Z^{m}and x a vector of n indeterminates), whose we search integer solutions. Recall that to decide if such a system has at least one integer solution is NP-complete (it is the NP-completeness of integer linear programming). We are interested here in describing and computing the set of solutions. Remark that these systems arise in pattern matching compilation theory. We will reduce our problem to the study of an equivalent problem, which is solving in non-negative integers the systems Ax = 0 (Frobenius problem). These systems arise in several sub-fields of equational rewriting theory, for instance in associative-commutative unification, or in Makanin algorithm. Remark that solving systems Ax = b, x ≥ 0, occurring in AC-unification, reduces also to Frobenius problem by adding a new variable to x. The non-negative integer solutions of Ax = 0 form a sub-monoid M of N^{n}, generated by its non zero minimal elements for the partial order (x_{1},…, x_{n}) ⪯ (y_{1},…, y_{n}) ⟺ ∀i, 1 ≤ i ≤ n, x_{i} ≤ y_{i}. They form a finite set. We will call this set “the Hilbert basis of M” (after [F.Giles and W.R.Pulleyblank 79]), and denote it by H(M). In the two next sections we will bound and compute the elements of H(M). The last section applies the previous results to the initial problem, i.e. the resolution of systems Ax ≤ b.

2 Bounds of H(M)

It is known since [J.Von zur Gathen and M.Sieveking 78] that if H(M) is non-empty, it contains an element with norm (for example ∥x∥_{∞}, where ∥x∥_{∞} = sup_{i}| x_{i} |) at most simply exponential in the size of A (for example n.m.(log∥A∥_{∞} + 2)). But we are interested here to uniformly bound the norms of the elements of H(M). Let

∥M∥_{∞} = sup_{x∈H(M)}∥x∥_{∞}

and

∥M∥_{1} = sup_{x∈H(M)}∥x∥_{1}

(with ∥x∥_{1} = ∑_{i} | x_{i} |). [I.Borosh and L.B.Treybig 76] have upper bounded ∥M∥_{∞} with an expression which is doubly exponential in the size of A. As far as we know, two simply exponential bounds exist to bound ∥M∥_{∞} or ∥M∥_{1}. We have given the first one in [L.Pottier 90]. The second one can be deduced from an rather unknown result of [J.L.Lambert 87], and has been found independently by [E.Domenjoud 90] in a better form. These two bounds are essentially different, because of their expressions and their proofs. We give here two new finest bounds, the first one inspired from [L.Pottier 90], the second one from [J.L.Lambert 87] and [A.Koscielski and L.Pacholski 90],. Recall that in the case of one equation (m = 1), [G.Huet 78] and [J.L.Lambert 87] have given bounds only depending of ∥A∥_{∞}. In the case of two equations, [J.F.Romeuf 89] gave a bound which is quadratic in the size of A.

2.1 First bound

This bound is inspirated by [L.Pottier 90]. Let ∥A∥_{1,∞} = sup_{i}{∑_{j} | a_{ij} |}, and let r be the rank of A.

Theorem 1

∥M∥_{1} ≤ (1 + ∥A∥_{1,∞})^{r} = B_{0}

Proof :

We can without restriction choose r independent equations of Ax = 0. Let x = (x_{1},…, x_{n}) be a non zero element of M, p = ∥x∥_{1}, and {e_{1},…, e_{n}} be the canonical basis of R^{n}. For every y in R^{n}, we note C_{y} the cube of volume 1 defined by

z ∈ C_{y} ⇔ z = y + ∑_{i=1}^{n} λ_{i}e_{i}, ∀i ∈ [1, n], λ_{i} ∈ [0, 1]

We will recursively define a sequence y^{0},…, y^{p} of N^{n} and a sequence z^{0},…, z^{p} of R^{n} verifying :

y^{0} = 0 ≺ y^{1} ≺ … ≺ y^{p} = x.
∀k ∈ [0, p − 1], ∃j, y^{k+1} = y^{k} + e_{j}
∀k ∈ [0, p], z^{k} ∈ C_{y^{k}} ∩ [0, x].

y^{0} = z^{0} = 0 are clearly convenient .

Suppose we have built y^{k}, and 0 ≤ k ≤ p − 1. [0, x] ∩ C_{y^{k}} is the set of all z which writes λ ∑_{i} x_{i}e_{i} with 0 ≤ λ ≤ 1 and

∀i ∈ [1, n], y_{i}^{k} ≤ λx_{i} ≤ y_{i}^{k} + 1

It is by hypothesis non empty, it is a segment. Take z^{k+1} = λ_{k} ∑_{i} x_{i}e_{i} its bound where λ_{k} is maximum. x is non zero, then there exists a j such that

λ_{k} = \frac{y_{j}^{k} + 1}{x_{j}} = inf_{i|x_{i}≠0}\{\frac{y_{i}^{k} + 1}{x_{i}}\}

Now, let y^{k+1} = y^{k} + ej. We have now

∀i ∈ [1, n], y_{i}^{k+1} ≤ λ_{k}x_{i} ≤ y_{i}^{k+1} + 1

and z^{k+1} belongs then to the cube C_{y^{k+1}}.

The points z^{k+1} and y^{k+1} are then correctly built. Finally, if k = p − 1, then y^{p} = x, because y^{p} ⪯ x and ∥y^{p}∥_{1} = p = ∥x∥_{1}, and we take z^{p} = x.

Now, let y^{′k} = z^{k} − y^{k}. We have now :

∀i, 0 ≤ y_{i}^{′k} ≤ 1

then, if (Ay^{k})_{i} is the i^{th} coordinate of Ay^{k} :

| (Ay^{k})_{i} | = | (Az^{k})_{i} − (Ay^{′k})_{i} | = | (Ay^{′k})_{i} |

As 0 ≤ y_{i}^{′k} ≤ 1, there is then at most ∑_{j} | a_{ij} | + 1 possible values for (Ay^{k})_{i} and then at most B_{0} distinct vectors Ay^{k}.

Now, suppose p > B_{0} . By the pigeon holes principle, it exists then i and j, i > j > 0 with Ay^{i} = Ay^{j}. Let z = y^{i} − y^{j}. We have now Az = 0. More we have 0 ≺ z ≺ x, and z ∈ M. Then x ∉ H(M).

□

2.2 Second bound

Let a_{ij} be the term of row i and of column j of the matrix A, and ∥A∥_{1} = ∑_{i,j} | a_{ij} |. Let D be the largest absolute value of the minors of A. [J.L.Lambert 87] gives the following result :

Theorem 2 (Lambert)

∥M∥_{∞} ≤ nD

From the proof of this result, and also from those of [A.Koscielski and L.Pacholski 90], also inspirated by [J.L.Lambert 87], we obtain the following improvements :

Theorem 3 Let D_{r} be the largest absolute value of the minors of order r of A .

∥M∥_{∞} ≤ (n − r)D_{r} = B_{1}

and then

∥M∥_{∞} ≤ (n − r) (\frac{∥A∥_{1}}{r} )^{r} = B_{2}

Remark : the first bound is the same as the bound of [E.Domenjoud 90], found independently. Proof :

Let C be the cone of R^{n} of non-negative real solutions of Ax = 0. Let C_{j} be its intersection with the hyperplane of equation x_{j} = 0. It is clear that C is the convex hull of the union of the cones C_{j}. We can recursively apply this decomposition of C to the C_{j}, while the dimension of built cones is largest than 1. C is then the convex hull of the union of these cones of dimension 1, called “edges” of C.

Every of these edges is then the set of non-negative solutions of a system of equations obtained by choosing r independent equations of Ax = 0, and by adding to them n − r − 1 equations of type x_{j} = 0 in order to keep the system of maximum rank, i.e. n − 1.

We can then obtain director vectors (with non-negative integer coefficients) of edges by computing the n minors of order n − 1 for every the preceding systems, that reduces to compute minors of order r of A.

Let g_{1},…, g_{k} be these vectors, which have then their coordinates upper bounded in absolute value by D_{r}.

M is then included in the non-negative cone that they generate (the linear combinations with non-negative real coefficients), which is exactly C, with dimension at most n − r. We have then, with the theorem of Carathéodory :

Ax = 0, x ≥ 0 ⟹ ∃j_{1},…, j_{n−r}, ∃α_{1},…, α_{n−r} ≥ 0, x = ∑_{l=1}^{n−r} α_{i}g_{j_{i}}

If now x is minimal, it is clear that the α_{i} are strictely smaller than 1. We obtain then the first part of the theorem.

The second part is a simple upper bound of the determinant of a square submatrix A′ of order r of A :

| det(A′) | ≤ ∏_{j} ∑_{i} | a′_{ij} | ≤ (\frac{∑_{i,j} | a′_{ij} |}{r})^{r} ≤ (\frac{∥A∥_{1}}{r} )^{r}

(The second inequality is an upper bound of geometric average with arithmetic average).

□

The bound B_{1} can be optimal , as we will see on examples, but it is not reasonably computable in practice: is it better to compute all the principal minors of A than to directly compute H(M), for example with the algorithm of [E.Contejean and H.Devie 89] which does not use a bound of H(M)?

2.3 Comparison of B_{0}, B_{1} and B_{2}

It is clear that these three bounds are simply exponential in the size of A. The following examples show that we can not compare in general the first and the last, the second being sometimes optimal, but being not computable in practice. We have the following inequalities :

∥x∥_{1} ≤ B_{0}
∥x∥_{1} ≤ nB_{2}
∥x∥_{∞} ≤ B_{0}
∥x∥_{∞} ≤ B_{2}

So we will study the behaviours of the ratios \frac{B_{2}}{B_{0}} (bounds of ∥x∥_{∞}) and \frac{nB_{2}}{B_{0}} (bounds of ∥x∥_{1}), when n or ∥A∥_{∞} increases to infinity.

2.3.1 Example 1

Let a be an integer greater or equal to 3 and A the matrix

( \begin{matrix}a & 1 − a \\ & ... & ... \\ & & a & 1 − a \end{matrix} )

where the non written coefficients are zero. We have r = m = n−1 and H(M) has only one element : ((a−1)^{n−1}, a(a−1)^{n−2},…, a^{n−1}). Then :

∥M∥_{∞} = B_{1} = a^{n−1}, ∥M∥_{1} = a^{n} − (a − 1)^{n}, B_{2} = (2a − 1)^{n−1}, B_{0} = (2a)^{n−1}

B_{1} is then optimal, B_{2} and B_{0} being very close. Asymptoticaly, we have finally :

lim_{n→∞} \frac{B_{2}}{B_{0}} = 0, lim_{n→∞} \frac{nB_{2}}{B_{0}} = 0, lim_{a→∞} \frac{B_{2}}{B_{0}} = 1, lim_{a→∞} \frac{nB_{2}}{B_{0}} = n

2.3.2 Example 2 : magic square matrices

A square matrix is called magic if the sums of its coefficients of a row and of a column are all equal. The magic square matrices of order k with non-negative coefficients are then the non-negative solutions of the system Ax = 0, where n = k^{2} + 1, m = 2k and :

A = ( \begin{matrix}1 & … & 1 & & & & & −1 \\ & & & 1 & … & 1 & & −1 \\ & & & & & & ... & ⋮ \\ 1 & & & 1 & & & & −1 \\ & ... & & & ... & & … & ⋮ \\ & & 1 & & & 1 & & −1 \end{matrix} )

where the non written coefficients are zero. We have r = 2k − 1. The Hilbert basis of M is the set of matrices of permutations of order k (cf [R.P.Stanley 83]). Then :

∥M∥_{∞} = 1, ∥M∥_{1} = k + 1, B_{1} ≥ k^{2} − 2k + 2

B_{2} = (k^{2} − 2k + 2) (\frac{2k(k + 1)}{2k − 1})^{2k−1} , B_{0} = (k + 2)^{2k−1}

and :

lim_{n→∞} \frac{B_{2}}{B_{0}} = ∞, lim_{n→∞} \frac{nB_{2}}{B_{0}} = ∞

which gives the inverse behaviour of the preceding example.

3 Algorithms

The subject of this section is the computation of all the elements of H(M). The first algorithms are based on the bounds of [G.Huet 78] and [J.L.Lambert 87] relative to one equation, extended to a system of equations, but giving then doubly exponential bounds. They are the followings :

Property 1 Let x = (x_{1},…, x_{p}, y_{1},…, y_{q}) be an element of the Hilbert basis of the equation

a_{1}x_{1} + … + a_{p}x_{p} + b_{1}y_{1} + … + b_{q}y_{q} = 0.

where the a_{i} are non-negative and the b_{j} are negative. Then :

∀i, | x_{i} | ≤ sup_{j}| b_{j} |

(Huet)

∑_{i} x_{i} ≤ sup_{j}| b_{j} |

(Lambert). (the part concerning the y_{j} is symmetric).

H(M) is then obtained by enumeration under the bound, and in the case of more than one equations, we iterate the method in injecting the solutions of treated equations in the followings (after eventualy having triangularized the matrix A). In the case of two equations [J.F.Romeuf 89] gives an original method for building a finite automaton enumerating H(M), and a quadratic bound in this case.

3.1 Algorithm of Contejean-Devie

[E.Contejean and H.Devie 89] have found a elegant algorithm which does not need any bound of H(M). The principle is the following. Let us order N^{n} by the order ⪯ defined before, and obtain a DAG (directed acyclic graph) of root 0. The algorithm enumerates a part of this DAG with the following principle : begin with 0, and if the current vertex is a non zero vector x such that for no one among its ancestors y we have A(x − y) = 0, visite its sons x + e_{j} verifying Ax.Ae_{j} ≤ 0 (the . denoting the scalar product of R^{n}). This algorithm suprizingly terminates and is complete. If we do not visite twice a vertex of the DAG, and keep only minimal solutions for ⪯, we then obtain H(M). Different refined versions of this algorithm exist, which eliminate early in the process some unusefull parts of the DAG. The only result of complexity about this algorithm is, to our knowledge, a consequence of [L.Baratchart and L.Pottier 89], which gives a doubly exponential bound on the number of visited vertices. This algorithm has good behaviour in practice, but is expensive if the elements of H(M) have large norms.

3.2 Algorithm of Domenjoud

In [E.Domenjoud 90] is described an algorithm which only builds solutions of Ax = 0 to compute minimal solutions (as the second algorithm that we present does). This recent algorithm would be interesting in pratice.

3.3 An algorithm inspirated by theorem 3

The analysis of the proof of the theorem 1 allows to modify the method of the algorithm of [E.Contejean and H.Devie 89] in only increment x by the e_{l} such that for every i, the i-th coordinate of A(x + e_{l}) is between − ∑_{j} a_{ij}^{+} and ∑_{j} a_{ij}^{−}. The generators are then all obtained as points of the sequences strictly increasing built similarly to the preceding algorithm.

3.4 Use of standard basis

We give here a new algorithm using the preceding bounds on ∥M∥_{∞} and ∥M∥_{1}, based on the theory of standard basis (or Gröbner basis). Let us recall basic notions of standard basis. For a polynomial P of the ring K[X_{1},…, X_{n}], we note in(f) the maximum monomial of f w.r.t a choosen admissible ordering on monomials (i.e. a total ordering stable by multiplication of monomials, and with the monomial 1 as minimal element). Then a familly F of polynomials of an ideal I is called a standard basis of I if and only if {in(f), f ∈ F } generates the ideal {in(f), f ∈ I }. A standard basis can be computed by completion algorithms (see [B.Buchberger 83], [A.Galligo 85]). In our problem, the idea is to see the columns of A as the exponents of monomials in m variables, and the solutions of Ax = 0 in Z^{n}as sysygies relative to these monomials. This idea has been introduced by [F.Ollivier 90] for computation of standard basis of sub-algebras. Then a computation of an appropriate standard basis gives a canonical rewriting system whose the inverse enumerates M by increasing norm. Finaly it suffices to only keep the minimal solutions for ⪯ and of norm smaller than inf{nB_{2}, B_{0}}. Let T, X_{1},…, X_{m}, Y_{1},…, Y_{n} be n + m + 1 variables, and k be an arbitrary field. We note a_{j} for the j^{th} column of A. For all α ∈ Z^{m} and β ∈ Z^{n} , we note X^{α} and Y^{β} the monomials X_{1}^{α_{1}} … X_{m}^{α_{m}} and Y_{1}^{β_{1}} … Y_{n}^{β_{n}}. α^{+} is the maximum of α and zero (for the partial order ⪯), and α^{−} is the maximum of −α and zero. Then α = α^{+} − α^{−}. For every j ∈ [1, n], we define a polynomial P_{j} in the ring R = k[T, X_{1},…, X_{m}, Y_{1},…, Y_{n}] :

P_{j} = X^{a_{j}^{+}} − Y_{j}X^{a_{j}^{−}}

Let I be the ideal of R generated by the P_{j} and the polynomial P_{0} = T Y_{1} … Y_{n} − 1, and J its trace (i.e. its intersection) on the ring R′ = k[Y_{1},…, Y_{n}]. Now, let BI be the reduced standard basis of I for the following ordering on the monomials of R : we compare first lexicographically the X_{i}, and in case of equality we use the degree order, and finally the lexicographic order. Let BJ be the set of polynomials of BI where the X_{i}’s and T do not appear. BJ is then a standard basis of the ideal J for the degree order (from a remark of D.Bayer and M.Stillman). More, its elements are differences of monomials (because those of BI are). Then let Y^{α_{k}} − Y^{β_{k}} be the elements of BJ , k ∈ [1, p] and Y^{α_{k}} being the leading monomials. Now, note ⟶ the rewriting relation corresponding to the division of polynomials by the standard basis BJ , and ⟶^{∗} its transitive reflexive closure. We write m1 ↓ m2 when two monomials m1 and m2 rewrite in the same monomial, or equivalently when m1 − m2 ⟶^{∗} 0. Then :

Property 2

∀x ∈ Z^{n} , Ax = 0 ⟺ Y^{x^{+}} − Y^{x^{−}} ∈ I ⟺ Y^{x^{+}} ↓ Y^{x^{−}}

Proof :

The first equivalence is easy to show by equational reasoning on the equations Y_{j} = X^{a_{j}^{+}}X^{−a_{j}^{−}} derived from polynomials P_{j} and with the equation P_{0} = 0 which only allows to eliminate monomials in factor in polynomials of I .

The second assertion is just the fact that a Gröbner basis is a canonical rewriting system equivalent to the relation P = Q ⇔ P − Q ∈ J . □

As a consequence :

Property 3

∀x ∈ N^{n} , x ∈ M ⟺ Y^{x} ⟶^{∗} 1

This last property allows to test if M is non reduced to {0} :

Theorem 4 The system Ax = 0 has a positive solution if and only if it exists in BI a polynomial of the form Y^{α} − 1.

More, we have an effective representation of M with of rewriting rules : Let SR_{M} the system of rewriting rules on monomials obtained in reversing the polynomials of BJ :

SR_{M} = {Y^{β_{1}}⟶Y^{α_{1}},…, Y^{β_{p}}⟶Y^{α_{p}}}

Note ⟶_{i} its rewriting relation (it is the symmetric of ⟶, and it is not noetherian). Then

x ∈ M ⟺ 1 ⟶_{i}^{∗} Y^{x}

We can then generate all the elements of M by exploration of the tree of rewritings of 1 by ⟶_{i}, and obtain H(M) in only keeping the minimal elements of degree smaller than the bounds nB_{2} and B_{0} (This method is complete because ⟶_{i} increases the degrees of monomials, and then the norms ∥.∥_{1} of the solutions). More precisely :

Theorem 5 The following algorithm stops and returns H(M) :

1. E := {1}

2. While ∃x ∈ E, y ∉ E, with x ⟶_{i} y, and deg(y) ≤ inf{nB_{2}, B_{0}}

Do E := E ∪ {y}

3. Return H(M) := minimal elements for ⪯ of vectors of exponents of monomials of E − {1}.

4 Application to Ax ≤ b

Now come back to the initial problem, i.e. the resolution of a system Ax ≤ b. Let C be the set of its solutions in Z^{n}. Then :

Corollary 1 It exists two finite parts C_{1} and C_{2} of Z^{n} such that :

x ∈ C ⇔ x = x_{1} + x_{2} + … + x_{k}, with x_{1} ∈ C_{1}, and x_{2},…, x_{k} ∈ C_{2}

and

∀x ∈ C_{1} ∪ C_{2}, ∥x∥_{1} ≤ (2 + ∥A∥_{1,∞} + ∥b∥_{∞})^{m}

Proof :

We will reduce the problem to solve in N a system of homogeneous equations.

Let ψ be an endomorphism of R^{n} which only change the signs of some coordinates of its argument, and ψ(A) the obtained matrix when changing the signs of the corresponding columns of A.

Let y = (y_{1},…, y_{m}) be a vector of m new variables, z a last variable, t the vector obtained in catenating x, y,and z, and let ϕ be the projection mapping t in x.

Let A′ be the matrix obtained in catenating ψ(A), the identity of order m and the opposite of b.

We have now clearly the equivalence :

Ax ≤ b ψ(x) ∈ N^{n} ⟺ ∃t ∈ N^{n+m+1}, A′t = 0, z = 1, x = ψ(ϕ(t))

More rank(A′) = m, and ∥A′∥_{1,∞} ≤ ∥A∥_{1,∞} + 1 + ∥b∥_{∞}.

Let H the Hilbert basis of A′t = 0, and C_{1}^{ψ} (resp. C_{2}^{ψ}) the image by ϕ of the elements of H such that z = 1 (resp. z = 0).

We take then C_{1} (resp. C_{2}) equal to the union of the C_{1}^{ψ} (resp. C_{2}^{ψ}) for the 2^{n} possible choices of ψ.

As ∥ψ(x)∥_{1} = ∥x∥_{1} and ∥ϕ(t)∥_{1} ≤ ∥t∥_{1}, we obtain the second part of the result. □

5 Acknowledgements

We would like to thank C.Traverso and A.Galligo for usefull discussions about first versions of algorithm 3.3., J.F.Romeuf who indicated us theorem 2, and the referees for their remarks on the first version of this paper.

References

[L.Baratchart and L.Pottier 89] “Un résultat sur les systèmes d’addition de vecteurs”, manuscript, INRIA Sophia Antipolis, France, feb. 1989.

[I.Borosh and L.B.Treybig 76] “Bounds of non-negative integral solutions of linear diophantine equations”, Proc. AMS v.55, n.2, march 1976.

[B.Buchberger 83] “Gröbner basis: an algorithmic method in polynomial ideal theory” Camp. Publ. Nr. 83-29. 0, nov. 1983.

[E.Contejean and H.Devie 89] “Solving systems of linear diophantine equations”, UNIF’89, proc. of the third international Workshop on unification, Lambrecht, RFA 89.

[E.Domenjoud 90] “Solving Systems fo Linear Diophantine Equations : An Algebraic Approch”, UNIF’90, International Workshop on Unification, Leeds UK, july 1990.

[A.Galligo 85] “Algorithmes de calcul de bases standard”, Preprint Université de Nice, France, 1985.

[F.Giles and W.R.Pulleyblank 79] “Total dual integrality and integer polyedra” Linear algebra and its applications, 25, pp191-196, 1979.

[G.Huet 78] “An algorithm to generate the basis of solutions to homogeneous linear diophantine equations”, Information Processing Letters, vol.3, No.7, 1978.

[A.Koscielski and L.Pacholski 90] “Exponent of periodicity of minimal solutions of word equations”, manuscript, University of Wroclaw, Poland, june 1990.

[J.L.Lambert 87] “Une borne pour les générateurs des solutions entières positives d’une équation diophantienne linéaire.” Comptes Rendus de l’Académie des Sciences de Paris, t.305, Série I, pp39-40, 1987.

[J.L.Lambert 87] “Un problème d’accessibilité dans les réseaux de Petri” Phd thesis, theorem I.5., p 18, University of Paris-Sud, Orsay, France, 1987.

[F.Ollivier 90] “Le problème de l’identifiabilité structurelle globale : approche théorique, méthodes effectives et bornes de complexité”, Phd Thesis, Ecole Polytechnique, France, june 1990.

[L.Pottier 90] “Bornes et algorithmes de calcul des générateurs des solutions de systèmes diophantiens linéaires”, internal report, INRIA, feb. 90, Comptes Rendus de l’Académie des Sciences de Paris, t.311, Série I, p813-816,1990.

[J.F.Romeuf 89] “Solutions of a linear diophantine system”, UNIF’89, proc. of the third international Workshop on unification, Lambrecht, RFA 89.

[R.P.Stanley 83] “Combinatorics and commutative algebra”, Progress in Mathematics, Birkäuser ed., 1983.

[L.B.Treybig 75] “Bounds in piecewise linear topology”, Trans.AMS, v.201, 1975. [J.Von zur Gathen and M.Sieveking 78] “A bound on solutions of linear integer equalities

and inequalities” Proc. AMS 72, pp155-158, 1978.