Algorithm to solve a system of linear inequalities over the natural numbers

232 Views Asked by At

I am looking for an algorithm to solve a system of linear inequalities of the form

$$A\,\vec n + \vec c \;\geq\; 0$$

for $\vec n$ where $A$ is a matrix of integers, $\vec c$ consists of integers and $\vec n$ of natural numbers (i.e. an algorithm to find a parametric representation of all integer points inside an in general unbounded convex polytope).

The solution should look like the output of Mathematica's Reduce command, e.g.

In[1]:= Reduce[1 - n2 + n3 >= 0 && -1 + n2 >= 0 && -n1 + n2 + n3 >= 0 && n1 - n2 >= 0 && n1 - n3 >= 0 && 1 - n1 + n2 >= 0 && 3 + n2 - n3 >= 0 && n1 >= 0 && n2 >= 0 && n3 >= 0,Integers]
Out[1]= C[1] ∈ Integers && C[1] >= 0 && ((n1 == 3 + C[1] && n2 == 2 + C[1] && n3 == 1 + C[1]) || (n1 == 2 + C[1] && n2 == 1 + C[1] && n3 == 2 + C[1]) || (n1 == 2 + C[1] && n2 == 1 + C[1] && n3 == 1 + C[1]) || (n1 == 1 + C[1] && n2 == 1 + C[1] && n3 == 1 + C[1]) || (n1 == 1 + C[1] && n2 == 1 + C[1] && n3 == C[1]))

(here the result is a union of lines)

or

In[1]:= Reduce[n1 >= 0 && n2 >= 0 && n1 - n2 + 1 >= 0, Integers]
Out[1]= (C[1] | C[2]) ∈ Integers && C[1] >= 0 && C[2] >= 0 && ((n1 == C[1] + C[2] && n2 == 1 + C[2]) || (n1 == C[1] + C[2] && n2 == C[2]))

(a union of two planes)

Note that the solutions describe infinite numbers of integer points.

Unfortunately, the Mathematica documentation doesn't mention which algorithm is used.

I am porting an existing Mathematica package (used in particle physics) for performance reasons to C++. The only missing ingredient is an equivalent to the Reduce command. The algorithm is needed to rewrite multi-fold sums of the form $$ \sum_{n_1=0}^\infty \cdots \sum_{n_N=0}^\infty f(\vec n)\, \theta(\vec a_1\cdot\vec n + \vec c_1) \cdots \theta(\vec a_J\cdot\vec n +\vec c_J) \,,\quad \theta(x)=\begin{cases} 1\;;\; x\geq 0 \\ 0\;;\; x<0 \end{cases} $$ to sums without $\theta$-functions.

Any hints are appreciated.

1

There are 1 best solutions below

1
On

Unfortunately, this problem is a subset of Integer Programming, which is NP-Complete, so you won't find anything necessarily "efficient".

Fortunately, now that you know the name, are plenty of libraries that have fancy heuristics for finding solutions. I know about lp_solve, but you can choose one that's more appropriate for your project.