Path: news.mathworks.com!not-for-mail
From: "John D'Errico" <woodchips@rochester.rr.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Integer solutions for a_1*x_1+a_2*x_2+a_3*x_3+ ...+a_4*x_n = k
Date: Tue, 18 Nov 2008 18:33:03 +0000 (UTC)
Organization: John D'Errico (1-3LEW5R)
Lines: 140
Message-ID: <gfv1ov$rrh$1@fred.mathworks.com>
References: <gfun9a$41t$1@fred.mathworks.com> <gfuq8v$n6o$1@fred.mathworks.com>
Reply-To: "John D'Errico" <woodchips@rochester.rr.com>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1227033183 28529 172.30.248.35 (18 Nov 2008 18:33:03 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Tue, 18 Nov 2008 18:33:03 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 869215
Xref: news.mathworks.com comp.soft-sys.matlab:501528


"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gfuq8v$n6o$1@fred.mathworks.com>...
> "Oriole " <oriole_ni@hotmail.com> wrote in message <gfun9a$41t$1@fred.mathworks.com>...
> > Hello, 
> > 
> > I am wondering how to write a function for solving
> > 
> > a_1*x_1+a_2*x_2+a_3*x_3+ ...+a_4*x_n = k
> > 
> > where
> > Inputs: 
> > n, k, and vector of non-negaive integer coefficients a=[a_1, a_2, ...a_n].     
> > k is some small positive integer, for example, k= 3, or k= 6, n is rather large positve integer, for example,, n = 20, or n= 60 or in some rare cases even n=100.
> > 
> > Output:
> > The unknowns x_1, x_2, ...x_n are non-negative integers.
> > 
> > I have asked a simplier version of this problem here, where a=[1, 1, ...1], and got lots of help from Walter, Bruno, Roger, John (thank you). 
> > 
> > But this problem... seems really difficult,  is it even possible to solve it in Matlab? Any ideas?
> > 
> > Thank you for your help!
> > Oriole
> 
> This is a Diophantine equation, as is the partitions
> problem, which is just a special case of the
> general Diophantine equation.
> 
> For n very large, this is a nasty problem. You can still
> solve it recursively as we have shown you, but that
> may be computationally expensive.
> 
> John

For kicks, I wrote a simple recursive code to compute
the integer solutions. For a problem with 15 coefficients,
it runs reasonably quickly, for small k (k = 5 here.)

A = rem(1:6,3)+1
A =
     2     3     1     2     3     1


Each row of x is a solution.

>> x = diophantine(A,5,0:5)
x =
     1     1     0     0     0     0
     2     0     1     0     0     0
     0     1     2     0     0     0
     1     0     3     0     0     0
     0     0     5     0     0     0
     0     1     0     1     0     0
     1     0     1     1     0     0
     0     0     3     1     0     0
     0     0     1     2     0     0
     1     0     0     0     1     0
     0     0     2     0     1     0
     0     0     0     1     1     0
     2     0     0     0     0     1
     0     1     1     0     0     1
     1     0     2     0     0     1
     0     0     4     0     0     1
     1     0     0     1     0     1
     0     0     2     1     0     1
     0     0     0     2     0     1
     0     0     1     0     1     1
     0     1     0     0     0     2
     1     0     1     0     0     2
     0     0     3     0     0     2
     0     0     1     1     0     2
     0     0     0     0     1     2
     1     0     0     0     0     3
     0     0     2     0     0     3
     0     0     0     1     0     3
     0     0     1     0     0     4
     0     0     0     0     0     5

>> tic,x = diophantine(rem(1:15,3)+1,5,0:5);toc
Elapsed time is 1.112639 seconds.
>> size(x)
ans =
   476    15

tic,x = diophantine(rem(1:20,3)+1,5,0:5);toc
Elapsed time is 2.848407 seconds.
>> size(x)
ans =
        1008          20

>> tic,x = diophantine(rem(1:25,3)+1,5,0:5);toc
Elapsed time is 8.448537 seconds.
>> size(x)
ans =
        2592          25

>> tic,x = diophantine(rem(1:30,3)+1,5,0:5);toc
Elapsed time is 19.793219 seconds.
>> size(x)
ans =
        5402          30


Expect this to be slow for most serious problems. 


function xsol = diophantine(A,c,xset)
% solves the diophantine equation dot(A,x) = c, where x is a member of xset
% 
% arguments
%  A    - row vector of length p, integer coefficients
%  c    - scalar (integer) constant
%  xset - set to sample from, sampling with replacement
%  
%  xsol - mxp array of solutions found

% check sizes, shapes
A = A(:)';
p = length(A);
xset = xset(:);
nx = length(xset);

xsol = zeros(0,p);
if length(A) == 1
  ind = find((A*xset - c) == 0);
  if ~isempty(ind)
    xsol = xset(ind);
  end
else
  % recursive calls
  for i = 1:nx
    if (c-A(end)*xset(i))>=0
      xsub = diophantine(A(1:(end-1)),c-A(end)*xset(i),xset);
      nsub = size(xsub,1);
      xsol = [xsol;[xsub,repmat(xset(i),nsub,1)]];
    end
  end
end

Have fun,
John