:

## The number of non-negative and positive solutions of linear Diophantine equations

This post is related to the this thread

The recursive procedure  PosIntSolve  finds the number of non-negative or positive solutions of any linear Diophantine equation

a1*x1+a2*x2+ ... +aN*xN = n  with positive coefficients a1, a2, ... aN .

Formal parameters: L is the list of coefficients of the left part, n  is the right part,  s (optional) is nonneg (by default) for nonnegint solutions and  pos  for positive solutions.

The basic ideas:

1) If you make shifts of the all unknowns by the formulas  x1'=x1-1,  x2'=x2-1, ... , xN'=xN-1  then  the number of positive solutions of the first equation equals the number of non-negative solutions of the second equation.

2) The recurrence formula (penultimate line of the procedure) can easily be proved by induction.

The code of the procedure:

restart;

PosIntSolve:=proc(L::list(posint), n::nonnegint, s::symbol:=nonneg)

local N, n0;

option remember;

if s=pos then n0:=n-`+`(op(L)) else n0:=n fi;

N:=nops(L);

if N=1 then if irem(n0,L[1])=0 then return 1 else return 0 fi; fi;

end proc:

Examples of use.

Finding of the all positive solutions of equation 30*a+75*b+110*c+85*d+255*e+160*f+15*g+12*h+120*i=8000:

st:=time():

PosIntSolve([30,75,110,85,255,160,15,12,120], 8000, pos);

time()-st;

13971409380

2.125

To test the procedure, solve (separately for non-negative and positive solutions) the simple equation  2*x1+7*x2+3*x3=2000  in two ways (by the  procedure and brute force method):

ts:=time():

PosIntSolve([2,7,3], 2000);

PosIntSolve([2,7,3], 2000, pos);

time()-ts;

47905

47334

0.281

ts:=time():

k:=0:

for x from 0 to 2000/2 do

for y from 0 to floor((2000-2*x)/7) do

for z from 0 to floor((2000-2*x-7*y)/3) do

if 2*x+7*y+3*z=2000 then k:=k+1 fi;

od: od: od:

k;

k:=0:

for x from 1 to 2000/2 do

for y from 1 to floor((2000-2*x)/7) do

for z from 1 to floor((2000-2*x-7*y)/3) do

if 2*x+7*y+3*z=2000 then k:=k+1 fi;

od: od: od:

k;

time()-ts;

47905

47334

50.063

Another example - the solution of the famous problem: how many ways can be exchanged \$ 1 using the coins of smaller denomination.

PosIntSolve([1,5,10,25,50],100);

292

Edit.  The code has been slightly edited

﻿