Re: linear equation with restrictions



could you show your sourcecode of that program?


Here we go, here is the main function, to be saved in one
m-file. The test program will follow.

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, his is Bruno's linprog, not MATLAB one in
optimization
% tool box, the result may be different.
sol = linprog([], [], a', b, cost, L, U);
if ~all(isinf(sol))
L(k) = max(L(k),sol(k)-epsilon);
end
cost = -basis(k,n);
sol = linprog([], [], a', b, cost, 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(1,n);
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
.

Relevant Pages

  • Re: Integer solutions for a_1*x_1+a_2*x_2+a_3*x_3+ ...+a_4*x_n = k
    ... This is recursion algo, but call for integer common division properties and linear programming to accelerate the computation. ... % x, a, lower, upper are n-dimensional vector ... set maxcount to finite value doesn't work as efficienly, ... elseif isinf% k has upper bound, but no lower bound ...
    (comp.soft-sys.matlab)
  • Re: Array issue has me stumped??
    ... ReDim Preserve sList ... The syntax for declaring an array dimension in VB is to specify the upper, ... By default the lower bound is 0 unless ...
    (microsoft.public.vb.general.discussion)
  • Re: Read previous line of textfile
    ... Then work from the upper to lower bound of the array. ... working backwards to the beginning etc etc. ...
    (microsoft.public.excel.programming)
  • Re: Cascading Windows
    ... Dim Resumes As String, FileArrayAs String, Num As Long, OpnFile As ... Dim Upper As Long 'Upper boundary of array ... Do While Len> 0 ' Loop to put files into an array ... and allowing me to use the bottom 50% for the other 2 apps). ...
    (microsoft.public.word.customization.menustoolbars)
  • Cascading Windows
    ... Dim Resumes As String, FileArrayAs String, Num As Long, OpnFile As ... Dim Upper As Long 'Upper boundary of array ... Do While Len> 0 ' Loop to put files into an array ... and allowing me to use the bottom 50% for the other 2 apps). ...
    (microsoft.public.word.customization.menustoolbars)
Loading