@Ravindra Pawar Replace eqt 2 to eqt 1 to remove one variable, you can solve it with gcd or more generally see this oldthread
% Test program:
% This example would take more or less a minute to solve
a = [30 12 4]
b = 1440;
b2 = 100;
bb = b - a(3)*b2;
aa = a(1:2) - a(3)
%
% solving
% c1*a1 + ... cn*an = b
% c >=0
sol = intlin(aa, bb);
% Recompute c from f and f
m = sol(:,1);
f = sol(:,2);
c = b2-f-m;
% Only keep solution that give integer c
keep = c >= 0;
f = f(keep);
m = m(keep);
c = c(keep);
T = table(m,f,c)
% Bruno
disp(30*m+12*f+4*c)
disp(m+f+c)
function x = intlin(a, b)
% function x = intlin(a, b)
%
% PURPOSE: solving integer equation:
%
% x(1)*a(1) + x(2)*a(2) + ... + x(n)*a(n) = b
% x >=0
%
% Input: a (1 x n) integer (could be negative)
% b (1 x 1) integer (could be negative)
% Output: x (n x 1) integer such that
% a*x = b, x>=0
%
% Author: Bruno Luong
% Last update: 24/July/2008
% Default value, empty solution
x = zeros(0,length(a));
if ~isscalar(b)
fprintf('b must be scalar\n');
return
end
% cast to double, and reshape as long thin array
a = double(a(:));
b = double(b);
if ~all(mod(a,1)==0) || mod(b,1)~=0
fprintf('a and b must be integers\n');
return
end
% d*a' = g
[g d] = gcdn(a);
if mod(b,g)==0
an = a/g;
bn = b/g;
% (d*a' = b) OR (d*an' = bn)
d = bn*d;
%
% General gcd final solution would be
% x = d + k;
% with k such that: k*an' = 0, and
% k>=-d (<=> x>=0)
kmin = ceil(-d);
kmax = +inf(size(a));
% Find all k such that k.an = 0
k = allintlin0(an, kmin, kmax);
x = bsxfun(@plus, d, k);
else % mode(b,g)~=0
fprintf('WARNING: there is no solution\n');
end
end
function x = allintlin0(a, lower, upper)
%
% x, a, lower, upper are n-dimensional vector
% List all interger x such that
% x.a = 0
% lower <= x <= upper
%
a = a(:);
n = length(a);
lower = lower(:);
upper = upper(:);
% Adjust upper and lower bounds by LP
L = lower;
U = upper;
b = 0;
epsilon = 1e-6;
for k=1:n
cost = basis(k,n);
% Beware, this is Bruno's linprog, not MATLAB one in optimization
% tool box, the result may be different.
sol = linprog(cost, [], [], a', b, L, U);
if ~all(isinf(sol))
L(k) = max(L(k),sol(k)-epsilon);
end
cost = -basis(k,n);
sol = linprog(cost, [], [], a', b, L, U);
if ~all(isinf(sol))
U(k) = min(U(k),sol(k)+epsilon);
end
end
L = ceil(L);
U = floor(U);
if all(~isinf(L)) && all(~isinf(U))
maxcount = Inf;
else % Limit the number of solutions that will be listed
% NOTE: set maxcount to finite value doesn't work as efficienly,
% as expected because the recursive engine might spend much of
% CPU time to look for unvalid solutions
% maxcount = 100;
fprintf('There is infinity of solutions\n');
x = NaN;
return
end
% Call the engine
count = 0;
x = ilinengine(a, L, U, b, count, maxcount);
if maxcount<inf
fprintf('WARNING: only %d solutions will be provided\n', maxcount);
x = x(1:min(maxcount,end),:); % clipping
end
end
function v = basis(k,n) % generate a k-th basis vector of dimension n
v = zeros(n,1);
v(k) = 1;
end
% Solver engine for integer x
% a*x = b
% lower <= x <= upper
% RESTRICTION:
% a must be primary array, i.e., they greatest common divisor is one
function x = ilinengine(a, lower, upper, b, count, maxcount)
% default value, empty result
x = zeros(0,length(a));
if count>maxcount
return
end
% Trivial solution with one variable
if length(a)==1
if mod(b,a(1))==0
xtmp = b / a(1);
if (xtmp >= lower) && (xtmp <=upper)
x = xtmp;
end
end
return
end
% preliminary check where as the sum b is possible
% This check is to speed up (?), and deos not affect the result
as = bsxfun(@times,a,[lower upper]);
as = sum(sort(as,2),1);
if b<as(1) || b>as(2)
return
end
% Greatest common divisor for the tail
g = gcdn(a(2:end));
% find r1, the inverse of a(1) in g-modulo group
[uno r1] = gcd(a(1),g);
clear uno; % should be 1
r=r1*mod(b,g);
s = (b - a(1)*r)/g;
a(2:end) = a(2:end)/g; % the tail is primary among them
kmin = ceil((lower(1)-r)/g);
kmax = floor((upper(1)-r)/g);
% perform a basic step
function newcount = kstep(k)
% Recursive call
xk = ilinengine(a(2:end), lower(2:end), upper(2:end), ...
s-k*a(1), count, maxcount);
nxk = size(xk,1);
x = [x; ...
(r+g*k)+zeros(nxk,1) xk]; % append the new solutions
newcount = count + nxk; % adjust the counter
end
if isinf(kmin) && isinf(kmax) % k is unbounded
for absk=0:inf
for k=unique([-absk absk])
if kstep(k)>maxcount
break
end
end
if count>maxcount
break
end
end
elseif isinf(kmin) % k has upper bound, but no lower bound
for k=kmax:-1:-inf
if kstep(k)>maxcount
break
end
end
else % k has lower bound
for k=kmin:kmax
if kstep(k)>maxcount
break
end
end
end
end
function [g varargout] = gcdn(varargin)
% function g = gcdn(a1, a2, ..., an);
%
% Return g, Greatest common divisor of a1, ... an
%
% [g c1 c2 ... cn]=gcdn(a1, a2, ..., an)
% return also c1, ..., cn
% So that a1*c1 + ... an*cn = g
%
% Compact calling form:
% [g c]=gcdn(a1, a2, ..., an) or [g c]=gcdn(a)
% assumes a and c are array
%
if nargin<2
a = varargin{1};
if length(a)<2
if a
g = abs(a);
if nargout>=2
varargout{1}=sign(a);
end
return
else
error('gcdn cannot compute for a = 0');
end
end
a = reshape(a,1,[]);
else
% Put all numbers in array
a=cell2mat(varargin);
end
g=a(1);
c=zeros(size(a));
c(1)=1;
for k=2:length(a)
[g cg c(k)]= gcd(g, a(k));
c(1:k-1) = c(1:k-1)*cg;
end
if nargout>=2
switch (nargout-1)
case 1,
varargout{1}=c;
case length(c)
varargout=num2cell(c);
otherwise
error('number of output is incompatible with input');
end
end
end
