Fastest way to solve non-negative linear diophantine equations
Let $A$ be a matrix in $M_{n \times m}(\mathbb{Z}_{\ge 0})$ without zero column.
Let $V$ be a vector in $\mathbb{Z}_{> 0}^m$.
What is the fastest way to find all the solutions $X \in \mathbb{Z}_{\ge 0}^n$ of the following equation? $$AX=V$$ By non-negativity and "without zero column" assumption, there always are finitely many solutions only.
The usual functions on SageMath for solving such a system may underuse the non-negativity. To prove so, I wrote a naive code finding all the solutions of such a system (see below), but I'm still looking for something significantly faster than it (see this crosspost on mathoverflow).
Here is the kind of system I need to solve (where $x_i \in \mathbb{Z}_{\ge 0}$), but with possibly more equations and variables.
[
5*x0 + 5*x1 + 5*x2 + 6*x3 + 7*x4 + 7*x5 == 24,
5*x1 + 7*x10 + 5*x6 + 5*x7 + 6*x8 + 7*x9 == 25,
5*x11 + 6*x12 + 7*x13 + 7*x14 + 5*x2 + 5*x7 == 25,
5*x12 + 6*x15 + 7*x16 + 7*x17 + 5*x3 + 5*x8 == 30,
5*x13 + 6*x16 + 7*x18 + 7*x19 + 5*x4 + 5*x9 == 35,
5*x10 + 5*x14 + 6*x17 + 7*x19 + 7*x20 + 5*x5 == 35,
5*x1 + 7*x10 + 5*x6 + 5*x7 + 6*x8 + 7*x9 == 25,
5*x21 + 5*x22 + 6*x23 + 7*x24 + 7*x25 + 5*x6 == 24,
5*x22 + 5*x26 + 6*x27 + 7*x28 + 7*x29 + 5*x7 == 25,
5*x23 + 5*x27 + 6*x30 + 7*x31 + 7*x32 + 5*x8 == 30,
5*x24 + 5*x28 + 6*x31 + 7*x33 + 7*x34 + 5*x9 == 35,
5*x10 + 5*x25 + 5*x29 + 6*x32 + 7*x34 + 7*x35 == 35,
5*x11 + 6*x12 + 7*x13 + 7*x14 + 5*x2 + 5*x7 == 25,
5*x22 + 5*x26 + 6*x27 + 7*x28 + 7*x29 + 5*x7 == 25,
5*x11 + 5*x26 + 5*x36 + 6*x37 + 7*x38 + 7*x39 == 24,
5*x12 + 5*x27 + 5*x37 + 6*x40 + 7*x41 + 7*x42 == 30,
5*x13 + 5*x28 + 5*x38 + 6*x41 + 7*x43 + 7*x44 == 35,
5*x14 + 5*x29 + 5*x39 + 6*x42 + 7*x44 + 7*x45 == 35,
5*x12 + 6*x15 + 7*x16 + 7*x17 + 5*x3 + 5*x8 == 30,
5*x23 + 5*x27 + 6*x30 + 7*x31 + 7*x32 + 5*x8 == 30,
5*x12 + 5*x27 + 5*x37 + 6*x40 + 7*x41 + 7*x42 == 30,
5*x15 + 5*x30 + 5*x40 + 6*x46 + 7*x47 + 7*x48 == 35,
5*x16 + 5*x31 + 5*x41 + 6*x47 + 7*x49 + 7*x50 == 42,
5*x17 + 5*x32 + 5*x42 + 6*x48 + 7*x50 + 7*x51 == 42,
5*x13 + 6*x16 + 7*x18 + 7*x19 + 5*x4 + 5*x9 == 35,
5*x24 + 5*x28 + 6*x31 + 7*x33 + 7*x34 + 5*x9 == 35,
5*x13 + 5*x28 + 5*x38 + 6*x41 + 7*x43 + 7*x44 == 35,
5*x16 + 5*x31 + 5*x41 + 6*x47 + 7*x49 + 7*x50 == 42,
5*x18 + 5*x33 + 5*x43 + 6*x49 + 7*x52 + 7*x53 == 48,
5*x19 + 5*x34 + 5*x44 + 6*x50 + 7*x53 + 7*x54 == 49,
5*x10 + 5*x14 + 6*x17 + 7*x19 + 7*x20 + 5*x5 == 35,
5*x10 + 5*x25 + 5*x29 + 6*x32 + 7*x34 + 7*x35 == 35,
5*x14 + 5*x29 + 5*x39 + 6*x42 + 7*x44 + 7*x45 == 35,
5*x17 + 5*x32 + 5*x42 + 6*x48 + 7*x50 + 7*x51 == 42,
5*x19 + 5*x34 + 5*x44 + 6*x50 + 7*x53 + 7*x54 == 49,
5*x20 + 5*x35 + 5*x45 + 6*x51 + 7*x54 + 7*x55 == 48
]
Here are explicit $A$ and $V$ from above system (in list form):
A=[ [5,5,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,5,0,0,0,0,5,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,5,0,0,0,0,5,0,0,0,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,5,0,0,0,0,5,0,0,0,5,0,0,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,5,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,5,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,5,0,0,0,0,5,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,0,0,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,5,0,0,0,0,5,0,0,0,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0], [0,0,0,5,0,0,0,0,5,0,0,0,5,0,0,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,0,0,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,6,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,7,7,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,7,7,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,7,7,0,0,0,0], [0,0,0,0,5,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,7,7,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,0,7,7,0,0], [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,0,7,7,0], [0,0,0,0,0,5,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,6,0,7,7,0,0,0,0,0,0,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,7,7,0,0,0,0], [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,0,7,7,0], [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,6,0,0,7,7] ]
V=[24,25,25,30,35,35,25,24,25,30,35,35,25,25,24,30,35,35,30,30,30,35,42,42,35,35,35,42,48,49,35,35,35,42,49,48]
I wrote a naive code finding all the solutions of such a system, then applied it to A, V above. It took 41 seconds to find all 5499 solutions. I'm looking for something significantly faster than this.
sage: %time LX=NonnegativeSolverPartition(A,V)
CPU times: user 40.7 s, sys: 0 ns, total: 40.7 s
Wall time: 40.9 s
sage: len(LX)
5499
Here is my code:
import copy
def NonnegativeSolverPartition(A,V):
global LX
WB=WeakBasis(A)
VB=VarBound(A,V)
m=len(A[0])
PP=[]
LX=[]
Or=[]
c=0
for l in WB:
[i,ll]=l
L=[A[i][j] for j in ll]
B=[VB[j] for j in ll]
n=V[i]
P=WeightedPartitionSolver(L,B,n)
Or.append([len(P),c])
c+=1
PP.append(P)
Or.sort()
for i in range(len(Or)):
if Or[i][0]!=1:
break
ii=Or[i][1]
for p in PP[ii]:
X=[-1 for i in range(m)]
PPP=copy.deepcopy(PP)
PPP[ii]=[p]
Fi=Filter(PPP,X,WB)
if Fi!=[]:
[PPPP,XX]=Fi
NonnegativeSolverPartitionInter(A,V,PPPP,WB,VB,XX)
return LX
def NonnegativeSolverPartitionInter(A,V,PP,WB,VB,X):
global LX
Or=[[len(PP[i]),i] for i in range(len(PP))]
Or.sort()
c=0
for i in range(len(Or)):
if Or[i][0]!=1:
c=1
break
ii=Or[i][1]
if c==0:
if -1 in X:
print('problem: -1 in X')
else:
if V==[sum([A[i][j]*X[j] for j in range(len(A[0]))]) for i in range(len(A))]:
LX.append(X)
else:
for p in PP[ii]:
PPP=copy.deepcopy(PP)
PPP[ii]=[p]
Fi=Filter(PPP,X,WB)
if Fi!=[]:
[PPPP,XX]=Fi
NonnegativeSolverPartitionInter(A,V,PPPP,WB,VB,XX)
def WeakBasis(A):
n=len(A)
m=len(A[0])
L=[]; F=[]
for i in range(n):
FF=[]
for j in range(m):
if A[i][j]!=0:
FF.append(j)
F.append(FF)
for j in range(m):
LL=[]
for i in range(n):
if A[i][j]!=0:
LL.append([len(F[i]),i])
else:
LL.append([0,i])
LL.sort(reverse=true)
L.append(LL[0])
return [[i,F[i]] for i in range(n)]
def VarBound(A,V):
M=[]
n=len(A)
m=len(A[0])
for j in range(m):
N=[]
for i in range(n):
if A[i][j]!=0:
N.append(V[i]//A[i][j])
M.append(min(N)+1)
return M
def WeightedPartitionSolver(L,B,n): # the entries of L must be non-negative
global P #B for upperbound coming from other equations (see VarBound)
P=[]
l=len(L)
S=[0 for i in range(l)]
m=n//L[0]+1
for i in range(min(m,B[0])):
S[0]=i
WeightedPartitionSolverInter(L,B,n-L[0]*i,1,S)
return P
def WeightedPartitionSolverInter(L,B,n,j,S):
global P
l=len(L)
if j==l:
if n==0:
P.append(S)
else:
if n==0:
P.append(S)
else:
m=n//L[j]+1
for i in range(min(m,B[j])):
SS=copy.deepcopy(S)
SS[j]=i
WeightedPartitionSolverInter(L,B,n-L[j]*i,j+1,SS)
def Filter(PP,X,W):
if [] in PP:
return []
XX=copy.deepcopy(X)
LL=[]
for i in range(len(PP)):
P=PP[i]
F=FixedPoints(P)
Wi=W[i]
for j in F:
P0j=P[0][j]
Wij=Wi[1][j]
if XX[Wij]==-1:
XX[Wij]=P0j
elif XX[Wij]!=P0j:
return []
for i in range(len(PP)):
L=[]
P=PP[i]
Wi=W[i]
lP0=len(P[0])
for p in P:
c=0
for j in range(lP0):
Wij=Wi[1][j]
if XX[Wij]!=-1:
if p[j]!=XX[Wij]:
c=1
break
if c==0:
L.append(p)
if L==[]:
return []
LL.append(L)
if [PP,X]==[LL,XX]:
return [LL,XX]
else:
return Filter(LL,XX,W)
def FixedPoints(P):
m=len(P)
n=len(P[0])
if m==1:
return [i for i in range(n)]
F=[]
for i in range(n):
c=0
for j in range(m):
if P[j][i]!=P[0][i]:
c=1
break
if c==0:
F.append(i)
return F
Use
Polyhedra
I suggest to employ
normaliz
backend and use.integral_points()
function. See https://doc.sagemath.org/html/en/refe...@Max Alekseyev Can you show me on above example? I added the explicit matrix $A$ and vector $V$ (in list form).