Thread Subject: linear equation with restrictions

Hi,
i would carry out the following equation
(figures are as examples):

A1 * c1 + A2 * c2 = b
A = [2.30 4.56]
b = 57.4

the core problem is I'am looking for c1 and c2 which
can/should assume only integer values and are declared in an
array c = int32[(d e)]
where d = int32(1:1000) AND e = int32(1:1000)

Now the term to get computing started:

                   c = A\b

to solve the equation where c is restricted to an array of
given integers (1 to 10000)does not involve the restrictions
I had given to refer to by making computations and gives me
solving like this:

c =

         0
   12.5877

I'would be glad to have tips to solve this!

"Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
message <g5sgla$b7r$1@fred.mathworks.com>...
> Hi,
> i would carry out the following equation
> (figures are as examples):
>
> A1 * c1 + A2 * c2 = b
> A = [2.30 4.56]
> b = 57.4

OK, in order to avoid work with float, I'll transform your
fisrt equation by multiplying it by 100

Ax100_1 * c1 + Ax100_2 * c2 = bx100 (eqt 1)

Ax100 = [230 456]
bx100 = 5740

Now search for greater common divisor of Ax1000 by MATLAB
function GCD:

[g D1 D2]=gcd(Ax100(1),Ax100(2))
% it is gcd(230,456) % g ==2, D1=-113, D2=57

Which mean
All the numbers generated by the lhs of (eqt 1) is
+ multiple of 2 (=g), and furthermore
+ g = 2 = -113*230 + 57*456 =
          -113*Ax100(1) + 47*Ax100(2) (eqt 2).

Now all you need is to see whereas bx100 is divisible by g.

if mod(bx100 ,g)==0 % yes
  C1 = D1 * bx100/g % -324310
  C2 = D2 * bx100/g % +163590
else
  error('No integer solution exists')
end

Verification:

C1*2.30+C2*4.56 % 57.4000

Bruno

Hi Bruno,
thanks indeed for your proposal to solve this!
I'm beginner of MatLab und therofore very thankfull for
initial support.
Now i have a question to this: g so far as i grasped is the
number of factors which are in a equation (C1, C2). Why is
'all the numbers generated by the ihs of (eqt1) must be
multiple of 2 (=g)'? And what is ihs?
Question2: c1 and c2 must be onle positiv, c1, c2 > 0.
How can this be secured?
Thanks!
Leon


"Bruno Luong" <b.luong@fogale.fr> wrote in message
<g5sjfa$fbg$1@fred.mathworks.com>...
> "Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
> message <g5sgla$b7r$1@fred.mathworks.com>...
> > Hi,
> > i would carry out the following equation
> > (figures are as examples):
> >
> > A1 * c1 + A2 * c2 = b
> > A = [2.30 4.56]
> > b = 57.4
>
> OK, in order to avoid work with float, I'll transform your
> fisrt equation by multiplying it by 100
>
> Ax100_1 * c1 + Ax100_2 * c2 = bx100 (eqt 1)
>
> Ax100 = [230 456]
> bx100 = 5740
>
> Now search for greater common divisor of Ax1000 by MATLAB
> function GCD:
>
> [g D1 D2]=gcd(Ax100(1),Ax100(2))
> % it is gcd(230,456) % g ==2, D1=-113, D2=57
>
> Which mean
> All the numbers generated by the lhs of (eqt 1) is
> + multiple of 2 (=g), and furthermore
> + g = 2 = -113*230 + 57*456 =
> -113*Ax100(1) + 47*Ax100(2) (eqt 2).
>
> Now all you need is to see whereas bx100 is divisible by g.
>
> if mod(bx100 ,g)==0 % yes
> C1 = D1 * bx100/g % -324310
> C2 = D2 * bx100/g % +163590
> else
> error('No integer solution exists')
> end
>
> Verification:
>
> C1*2.30+C2*4.56 % 57.4000
>
> Bruno
>

"Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
message <g5td8a$s9q$1@fred.mathworks.com>...
> Why is
> 'all the numbers generated by the ihs of (eqt1) must be
> multiple of 2 (=g)'? And what is ihs?

This is a well known theorem in algebra, it saids that all
the ring defined by {c1*a1 + c2*a2; c1, c2 in N} is an
"principal ideal" and the generator is the gcd of a1, a2.


> Question2: c1 and c2 must be onle positiv, c1, c2 > 0.
> How can this be secured?

In your case it cannot, or in other word there is no
solution that satisfies this constraint.

Here is why: If {d1,d2} is one couple integer solution (as
provided by gcd function) of (eqt2):

d1*a1 + d2*a2 = g, (eqt2)

then all solutions of (eqt 2) are (e1,e2) with

e1 = d1 + k*a2/g
e2 = d2 - k*a1/g
    for k=..., -1, 0, 1, 2, ....

In your example, none of the solution has both (e1,e2)
positive. So (c1,c2) must have alternated sign.

Other example this might be possible, but not in your example.

Bruno

"Bruno Luong" <b.luong@fogale.fr> wrote in message
<g5tep3$g7r$1@fred.mathworks.com>...
> "Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
> message <g5td8a$s9q$1@fred.mathworks.com>...
> > Why is
> > 'all the numbers generated by the ihs of (eqt1) must be
> > multiple of 2 (=g)'? And what is ihs?
>
> This is a well known theorem in algebra, it saids that all
> the ring defined by {c1*a1 + c2*a2; c1, c2 in N} is an
> "principal ideal" and the generator is the gcd of a1, a2.
>
>
> > Question2: c1 and c2 must be onle positiv, c1, c2 > 0.
> > How can this be secured?
>
> In your case it cannot, or in other word there is no
> solution that satisfies this constraint.
>
> Here is why: If {d1,d2} is one couple integer solution (as
> provided by gcd function) of (eqt2):
>
> d1*a1 + d2*a2 = g, (eqt2)
>
> then all solutions of (eqt 2) are (e1,e2) with
>
> e1 = d1 + k*a2/g
> e2 = d2 - k*a1/g
> for k=..., -1, 0, 1, 2, ....
>
> In your example, none of the solution has both (e1,e2)
> positive. So (c1,c2) must have alternated sign.
>
> Other example this might be possible, but not in your example.
>
> Bruno

Well, i might be bold and turn my qeustion more precisely,
indeed sorry for having failed this previosly!

The example with given figures means only to have made
problem clearly.

The assignment for MatLab would be to get only solvings
outputed which are positiv. If integer and positiv solutions
cant be found, then should be equally given a message like
you did in the code before: else
> error('No (positive)integer solution exists')

Probably there is a function for that or it can be done
through conditional looping, anyway you might have a right idea.
Leon

"Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
message <g5thqi$h7d$1@fred.mathworks.com>...
> then should be equally given a message like
> you did in the code before: else
> > error('No (positive)integer solution exists')
>

It's not that difficult, all we have to do is find a bracket
of an interval for 'k' so that both e1 and e2 are positives.
If the bracket contains at least one integer k, then there
is a solution (pick the k).

In the above example the bracket interval is
[0.495614035087719 0.495652173913044], and of course no
integer k belongs to it.

I let you figure out how to program this part in MATLAB.

Bruno

"Bruno Luong" <b.luong@fogale.fr> wrote in message
<g5tj02$rj4$1@fred.mathworks.com>...
> "Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
> message <g5thqi$h7d$1@fred.mathworks.com>...
> > then should be equally given a message like
> > you did in the code before: else
> > > error('No (positive)integer solution exists')
> >
>
> It's not that difficult, all we have to do is find a bracket
> of an interval for 'k' so that both e1 and e2 are positives.
> If the bracket contains at least one integer k, then there
> is a solution (pick the k).
>
> In the above example the bracket interval is
> [0.495614035087719 0.495652173913044], and of course no
> integer k belongs to it.
>
> I let you figure out how to program this part in MATLAB.
>
> Bruno

Well done!
my really last question to:
in case there are more variables than 2 in the equation
(a1 * b2 + a2 * b2 + a3 * b3 = c) the solving with refering
to GCD - Function of MATLAB falls down becaus gcd operates
only with 2 elements?
Thanks!
Leon

"Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
message <g5tleq$siu$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.fr> wrote in message
> <g5tj02$rj4$1@fred.mathworks.com>...
> > "Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
> > message <g5thqi$h7d$1@fred.mathworks.com>...
> > > then should be equally given a message like
> > > you did in the code before: else
> > > > error('No (positive)integer solution exists')
> > >
> >
> > It's not that difficult, all we have to do is find a bracket
> > of an interval for 'k' so that both e1 and e2 are positives.
> > If the bracket contains at least one integer k, then there
> > is a solution (pick the k).
> >
> > In the above example the bracket interval is
> > [0.495614035087719 0.495652173913044], and of course no
> > integer k belongs to it.
> >
> > I let you figure out how to program this part in MATLAB.
> >
> > Bruno
>
> Well done!
> my really last question to:
> in case there are more variables than 2 in the equation
> (a1 * b2 + a2 * b2 + a3 * b3 = c) the solving with refering
> to GCD - Function of MATLAB falls down becaus gcd operates
> only with 2 elements?

No big deal, first find the gcd (g12) of the first two
coefs, then find the gcd (g123) of the g12 and a3, etc...

Bruno

"Bruno Luong" <b.luong@fogale.fr> wrote in message
<g5tlrt$44m$1@fred.mathworks.com>...
> "Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
> message <g5tleq$siu$1@fred.mathworks.com>...
> > "Bruno Luong" <b.luong@fogale.fr> wrote in message
> > <g5tj02$rj4$1@fred.mathworks.com>...
> > > "Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
> > > message <g5thqi$h7d$1@fred.mathworks.com>...
> > > > then should be equally given a message like
> > > > you did in the code before: else
> > > > > error('No (positive)integer solution exists')
> > > >
> > >
> > > It's not that difficult, all we have to do is find a
bracket
> > > of an interval for 'k' so that both e1 and e2 are
positives.
> > > If the bracket contains at least one integer k, then there
> > > is a solution (pick the k).
> > >
> > > In the above example the bracket interval is
> > > [0.495614035087719 0.495652173913044], and of course no
> > > integer k belongs to it.
> > >
> > > I let you figure out how to program this part in MATLAB.
> > >
> > > Bruno
> >
> > Well done!
> > my really last question to:
> > in case there are more variables than 2 in the equation
> > (a1 * b2 + a2 * b2 + a3 * b3 = c) the solving with refering
> > to GCD - Function of MATLAB falls down becaus gcd operates
> > only with 2 elements?
>
> No big deal, first find the gcd (g12) of the first two
> coefs, then find the gcd (g123) of the g12 and a3, etc...
>
> Bruno
>
>
Pardon me for being stupid, i need some clarification:
if i find gcd for a1 and a2 by this (refering to your code):
[g D1 D2]=gcd(Ax100(1),Ax100(2))
> % it is gcd(230,456) % g ==2, D1=-113, D2=57
then how should i proceed with coefs for a3 etc?
[g ? D3]=gcd(Ax100(?),Ax100(3))

Try this:

[g12 F1 F2]=gcd(Ax100(1),Ax100(2))
[g F12 F3]=gcd(g12,Ax100(3))

I little calculation is required to compute D1, D2, D3 from
F.

Bruno

"Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in message
<g5tleq$siu$1@fred.mathworks.com>...
> Well done!
> my really last question to:
> in case there are more variables than 2 in the equation
> (a1 * b2 + a2 * b2 + a3 * b3 = c) the solving with refering
> to GCD - Function of MATLAB falls down becaus gcd operates
> only with 2 elements?
> Thanks!
> Leon

  Leon, I understand that you wish to solve the equation

 a1*b1 + a2*b2 + a3*b3 = c

where a1, a2, a3, and c are known and without loss of generality can be
assumed to be positive integers, and where b1, b2, and b3 are unknown but
are restricted to also being positive integers.

  If

 [g12,-,-] = gcd(a1,a2)

is done, you already know how to find b1 and b2 which can make

 a1*b1 + a2*b2

equal to any desired multiple of g12. If

 [g123,-,-] = gcd(g12,a3)

is done, you can then find from this a set of b1, b2, and b3 values which can
make

 a1*b1 + a2*b2 + a3*b3

equal to any desired multiple of g123. If c is not a multiple of g123, then the
original equation above has no solution. Otherwise b1, b2, and b3 solutions
can always be found, but one or two of them will typically be negative.

  It is not hard to also show from the above reasoning that b3 can be adjusted
up or down by any desired multiple of g12/g123 if corresponding changes
are made to b1 and b2 that cause the left side of the equation to remain
unchanged. That is, if the multiple of g12 that a1*b1+a2*b2 is equal to is
changed in the opposite direction by the same multiple of a3/g123.

  This should then be done in such a way that b3 assumes its minimum
possible positive value, and that will give c-a3*b3 its maximum possible
value. From what Bruno has told you, you already know how to then look for
possible solutions to the two variable problem

 a1*b1 + a2*b2 = c3 = c-a3*b3

where b3 has been so determined. This will tell you if any valid positive
solutions exist and if so, will give you one of them.

  If b3 were to be increased from this minimum value, then c-a3*b3 is
decreased and the number of possible positive b1 and b2 solutions can only
decrease. If you also have an upper bound for your b's, as in your first
article, the latter process of increasing b3 from its minimum might
conceivably be necessary to keep all b's within upper bounds simultaneously.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
wrote in message <g5ujb2$7c1$1@fred.mathworks.com>...
> "Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
message
> <g5tleq$siu$1@fred.mathworks.com>...
> > Well done!
> > my really last question to:
> > in case there are more variables than 2 in the equation
> > (a1 * b2 + a2 * b2 + a3 * b3 = c) the solving with refering
> > to GCD - Function of MATLAB falls down becaus gcd operates
> > only with 2 elements?
> > Thanks!
> > Leon
>
> Leon, I understand that you wish to solve the equation
>
> a1*b1 + a2*b2 + a3*b3 = c
>
> where a1, a2, a3, and c are known and without loss of
generality can be
> assumed to be positive integers, and where b1, b2, and b3
are unknown but
> are restricted to also being positive integers.
>
> If
>
> [g12,-,-] = gcd(a1,a2)
>
> is done, you already know how to find b1 and b2 which can make
>
> a1*b1 + a2*b2
>
> equal to any desired multiple of g12. If
>
> [g123,-,-] = gcd(g12,a3)
>
> is done, you can then find from this a set of b1, b2, and
b3 values which can
> make
>
> a1*b1 + a2*b2 + a3*b3
>
> equal to any desired multiple of g123. If c is not a
multiple of g123, then the
> original equation above has no solution. Otherwise b1,
b2, and b3 solutions
> can always be found, but one or two of them will typically
be negative.
>
> It is not hard to also show from the above reasoning
that b3 can be adjusted
> up or down by any desired multiple of g12/g123 if
corresponding changes
> are made to b1 and b2 that cause the left side of the
equation to remain
> unchanged. That is, if the multiple of g12 that
a1*b1+a2*b2 is equal to is
> changed in the opposite direction by the same multiple of
a3/g123.
>
> This should then be done in such a way that b3 assumes
its minimum
> possible positive value, and that will give c-a3*b3 its
maximum possible
> value. From what Bruno has told you, you already know how
to then look for
> possible solutions to the two variable problem
>
> a1*b1 + a2*b2 = c3 = c-a3*b3
>
> where b3 has been so determined. This will tell you if
any valid positive
> solutions exist and if so, will give you one of them.
>
> If b3 were to be increased from this minimum value, then
c-a3*b3 is
> decreased and the number of possible positive b1 and b2
solutions can only
> decrease. If you also have an upper bound for your b's,
as in your first
> article, the latter process of increasing b3 from its
minimum might
> conceivably be necessary to keep all b's within upper
bounds simultaneously.
>
> Roger Stafford
>
Roger,
thanks to your answer. I'm afraid i have but to clarify some
detail which Bruno almost tryed to satisfy.
(It's going about the same equation
(a1 * b2 + a2 * b2 + a3 * b3 = c) you understand rightly)
Say, if g1 and g2 through gcd(a1, a2) have been got
(output would be D1 and D2), then
to get g3 it's to state: gcd(g12,a3), where i understand
g12 to be Greatest Common Deviser for A1 and A1 (g12).

Where is have to be deduced g12 from?
Is it correct:

            g12 = D1*A1 + D2*A2 ?
Thanks
Leon

For you Leon,

Bruno


function [g varargout] = gcdn(varargin)
% function g = gcdn(a1, a2, ..., an); OR
%
% 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
% [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
        error('gcdn requires two or more inputs');
    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

Example:

a = [75 180 35];
[g c]=gcdn(a)

sum(c.*a) % is g

Bruno

"Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in message
<g5v5rm$iav$1@fred.mathworks.com>...
> Roger,
> thanks to your answer. I'm afraid i have but to clarify some
> detail which Bruno almost tryed to satisfy.
> (It's going about the same equation
> (a1 * b2 + a2 * b2 + a3 * b3 = c) you understand rightly)
> Say, if g1 and g2 through gcd(a1, a2) have been got
> (output would be D1 and D2), then
> to get g3 it's to state: gcd(g12,a3), where i understand
> g12 to be Greatest Common Deviser for A1 and A1 (g12).
>
> Where is have to be deduced g12 from?
> Is it correct:
>
> g12 = D1*A1 + D2*A2 ?
> Thanks
> Leon

  I'll give a specific example, Leon, of what I was describing in the earlier
article to clarify things.

  Suppose we have a1 = 770, a2 = 105, a3 = 60, and c = 2215, so we are
trying to solve

 770*b1 + 105*b2 + 60*b3 = 2215.

First we do

 [g12,d1,d2] = gcd(a1,a2)

which yields g12 = 35, d1 = 1, and d2 = -7. Next we do

 [g123,e12,e3] = gcd(g12,a3)

which produces g123 = 5, e12 = -5, and e3 = 3. We notice that c is indeed a
multiple of g123: c = 2215 = 443*5, so we proceed. From above we know
that

 e12*g12+e3*a3 = g123 [-5*35+3*60=5]

Therefore

 (443*e12)*g12+(443*e3)*a3 = 443*g123, or
 -2215*g12 + 1329*a3 = 2215

We can subtract any multiple of g12/g123=7 from 1329 if we add the same
multiple of a3/g123=12 to -2215, so we use the multiple 189 which gives

 53*g12 + 6*a3 = 2215

This is the least positive value b3 could possibly have.

  Hence we now wish to solve

 a1*b1 + a2*b2 = c-6*a3 = 1855 = 53*g12

From above we know that

 d1*a1 + d2*a2 = g12 [1*770-7*105=35]

so

 (53*d1)*a1 + (53*d2)*a2 = 53*g12, or
 53*a1 - 371*a2 = 1855

We can add any multiple of a1/g12=22 to -371 provided we subtract the
same multiple of a2/g12=3 from 53, so we choose 17 which gives

 2*a1+3*a2 = 1855

  Therefore we have a solution with all positive b's:

b1 = 2, b2 = 3, and b3 = 6.

  Now you need to translate all the above reasoning into a general form that
can handle any situation.

  I hope this example has been of assistance to you.

Roger Stafford

I think finding all solutions (if they are finite) of this
natural-integer linear equations with three or more
variables is a fun problem. Personally I would use a
recursive algorithm for it.

I'm not sure if there is a convenient way to parametrize the
set of infinite solutions. It would be nice if someone who
works with in number group/ring theory or
integer-programming can shed a light.

Bruno

For fun, I finally made a program to find all the solutions
of the integer equations:

given a: (1 x n), b : (1 x 1)

solve x integer, (1 x n) such that
  a*x' = b
  x >= 0

Example:

a = [770 105 60] % Like Roger's example
c = 7755

x =

     0 3 124
     0 7 117
     0 11 110
     0 15 103
     0 19 96
     0 23 89
     0 27 82
     0 31 75
     0 35 68
     0 39 61
     0 43 54
     0 47 47
     0 51 40
     0 55 33
     0 59 26
     0 63 19
     0 67 12
     0 71 5
     3 1 89
     3 5 82
     3 9 75
     3 13 68
     3 17 61
     3 21 54
     3 25 47
     3 29 40
     3 33 33
     3 37 26
     3 41 19
     3 45 12
     3 49 5
     6 3 47
     6 7 40
     6 11 33
     6 15 26
     6 19 19
     6 23 12
     6 27 5
     9 1 12
     9 5 5

Not 100% sure whereas my program is BUG free and would not
miss any solution.

Bruno

"Bruno Luong" <b.luong@fogale.fr> wrote in message
<g62uv6$r2k$1@fred.mathworks.com>...
> I think finding all solutions (if they are finite) of this
> natural-integer linear equations with three or more
> variables is a fun problem.

No one else wants to tackle the integer problem "Help me to
spend all my money")? Roger?

I would like to see good ideas from you in solving this
problem (mine is rather brute-force).

I propose this a test example for us to compare notes:

a = [770 105 60 85 50]
b = 9435

Find *all* x (5x1), integer arrays such that
  a*x= b and x>=0.

Bruno

"Bruno Luong" <b.luong@fogale.fr> wrote in message
<g68au9$b8u$1@fred.mathworks.com>...
> For fun, I finally made a program to find all the solutions
> of the integer equations:
>
> given a: (1 x n), b : (1 x 1)
>
> solve x integer, (1 x n) such that
> a*x' = b
> x >= 0
>
> Example:
>
> a = [770 105 60] % Like Roger's example
> c = 7755
>
> x =
>
> 0 3 124
> 0 7 117
> 0 11 110
> 0 15 103
> 0 19 96
> 0 23 89
> 0 27 82
> 0 31 75
> 0 35 68
> 0 39 61
> 0 43 54
> 0 47 47
> 0 51 40
> 0 55 33
> 0 59 26
> 0 63 19
> 0 67 12
> 0 71 5
> 3 1 89
> 3 5 82
> 3 9 75
> 3 13 68
> 3 17 61
> 3 21 54
> 3 25 47
> 3 29 40
> 3 33 33
> 3 37 26
> 3 41 19
> 3 45 12
> 3 49 5
> 6 3 47
> 6 7 40
> 6 11 33
> 6 15 26
> 6 19 19
> 6 23 12
> 6 27 5
> 9 1 12
> 9 5 5
>
> Not 100% sure whereas my program is BUG free and would not
> miss any solution.
>
> Bruno

Bruno,
you said you made a program for outputting a set of all
possible combinations of solvings.

could you show your sourcecode of that program?

Leon

>
> 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

Test program:

a=[770 105 60 85 50];
b = 9435;

%
% solving
% c1*a1 + ... cn*an = b
% c >=0

c = intlin(a, b);

if isnan(c) % unbounded, infinity solutions
elseif ~isempty(c)

    fprintf('\nSolutions c:\n\n');
    disp(c);
    fprintf('\nNumber of solutions = %d\n\n', size(c,1));
   
    bverif = c*a';
    bverifu = unique(bverif);
    
    if length(bverifu)~=1 || bverifu~=b
        fprintf('Wrong solution(s)\n');
        disp(c(bverif~=b,:));
        save test_inlin_debug.mat a cc b
    end
    
else % mode(b,g)==0
    fprintf('There is no solution\n');
end

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in
message <g6ic71$jqt$1@fred.mathworks.com>...
> Test program:
>
> a=[770 105 60 85 50];
> b = 9435;
>
> %
> % solving
> % c1*a1 + ... cn*an = b
> % c >=0
>
> c = intlin(a, b);
>
> if isnan(c) % unbounded, infinity solutions
> elseif ~isempty(c)
>
> fprintf('\nSolutions c:\n\n');
> disp(c);
> fprintf('\nNumber of solutions = %d\n\n', size(c,1));
>
> bverif = c*a';
> bverifu = unique(bverif);
>
> if length(bverifu)~=1 || bverifu~=b
> fprintf('Wrong solution(s)\n');
> disp(c(bverif~=b,:));
> save test_inlin_debug.mat a cc b
> end
>
> else % mode(b,g)==0
> fprintf('There is no solution\n');
> end
>
Hi Bruno,
Thank you for that!
I'am just trying to test it, now it comes to the following
error message:
??? Error using ==> linprog at 175
The number of rows in A must be the same as the length of b.

In your sourcecode I found a comment:

'% 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);'

Wouldn't run it with some Stadard- MatLab of mine?
Leon

"Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
message <g6l6lu$40r$1@fred.mathworks.com>...

>
> '% 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);'
>
> Wouldn't run it with some Stadard- MatLab of mine?
> Leon

I don't know. I do not have MATLAB linprog and I try to
simulate it with my own LP solver obviously without succeed.

What I want to solve is this continue pb
  minimize cost * x with constraints
     L <= x <= U and a'*x = b

I'm not sure what is the correct syntax of MATLAB linprog to
ignore linear-inequality constraint, I put empty matrix but
MATLAB obviously doesn't like. Hope someone can help.

Otherwise could you try to replace with the following line
with a dummy inequality matrix:

sol = linprog(zeros(size(a')), 1, a', b, cost, L, U);

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in
message <g6l7nl$hpr$1@fred.mathworks.com>...

>
> sol = linprog(zeros(size(a')), 1, a', b, cost, L, U);
>

My bad, the cost function should be the first argument:

sol = linprog(cost', [], [], a', b, L, U);

OR dummy syntax

sol = linprog(cost', zeros(size(a')), 1, a', b, L, U);

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in
message <g6lh5j$laa$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in
> message <g6l7nl$hpr$1@fred.mathworks.com>...
>
> >
> > sol = linprog(zeros(size(a')), 1, a', b, cost, L, U);
> >
>
> My bad, the cost function should be the first argument:
>
> sol = linprog(cost', [], [], a', b, L, U);
>
> OR dummy syntax
>
> sol = linprog(cost', zeros(size(a')), 1, a', b, L, U);
>
> Bruno

Bruno,
with your last syntax it did!

Now would you know how to get alter the syntax to get
results where (a1,a2, etc) > 0
Leon

"Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
message <g6qj0s$mut$1@fred.mathworks.com>...

>
> Now would you know how to get alter the syntax to get
> results where (a1,a2, etc) > 0
> Leon
>
>

Sorry, I don't understand your question. a1, a2 are the
known inputs, what you want to do with them?

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in
message <g6qm7f$5mc$1@fred.mathworks.com>...
> "Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
> message <g6qj0s$mut$1@fred.mathworks.com>...
>
> >
> > Now would you know how to get alter the syntax to get
> > results where (a1,a2, etc) > 0
> > Leon
> >
> >
>
> Sorry, I don't understand your question. a1, a2 are the
> known inputs, what you want to do with them?
>
> Bruno

Sorry, ich meant c>0 not a.

"Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
message <g6sgu6$f5s$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in
> message <g6qm7f$5mc$1@fred.mathworks.com>...

>
> Sorry, ich meant c>0 not a.
>

I see. To ways to do it.

1. Brute force: find all solutions with c>=0, then filter
out by find.

2. Neater way: c>0 means cm1:=c-1 >=0.

So solving
(A) a*c = b such that c > 0

is equivalent to solving

(B) a*cm1 = b-sum(a) such that cm1>=0,
Then take c=cm1+1.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in
message <g6st08$o4i$1@fred.mathworks.com>...
> "Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
> message <g6sgu6$f5s$1@fred.mathworks.com>...
> > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in
> > message <g6qm7f$5mc$1@fred.mathworks.com>...
>
> >
> > Sorry, ich meant c>0 not a.
> >
>
> I see. To ways to do it.
>
> 1. Brute force: find all solutions with c>=0, then filter
> out by find.
>
> 2. Neater way: c>0 means cm1:=c-1 >=0.
>
> So solving
> (A) a*c = b such that c > 0
>
> is equivalent to solving
>
> (B) a*cm1 = b-sum(a) such that cm1>=0,
> Then take c=cm1+1.
>
> Bruno

Sorry for pestering you with that but i'm unclear
where would [a*cm1 = b-sum(a) and c=cm1+1] be belonging?
If this refers to MainProg in m-file, thenn to what
syntax sequence have to be that(i thought just you know your
own code anyway better than me bein unexperienced in MathLab).

Thank you
Leon

"Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
message <g6ubj2$6mb$1@fred.mathworks.com>...

>
> Sorry for pestering you with that but i'm unclear
> where would [a*cm1 = b-sum(a) and c=cm1+1] be belonging?
> If this refers to MainProg in m-file, thenn to what
> syntax sequence have to be that

I let it to you to figure it out Leon.

> (i thought just you know your
> own code anyway better than me bein unexperienced in MathLab).
>

I'm sure you can do it on your own. Don't be hesitate to
play with the code.

Start with simple example, for example:

a = [770 105];
b = 4830;

There is three solutions c>=0, But only two c>0.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in
message <g6uj4u$n2m$1@fred.mathworks.com>...
> "Leon galushko" <leonid.galushko@rwth-aachen.de> wrote in
> message <g6ubj2$6mb$1@fred.mathworks.com>...
>
> >
> > Sorry for pestering you with that but i'm unclear
> > where would [a*cm1 = b-sum(a) and c=cm1+1] be belonging?
> > If this refers to MainProg in m-file, thenn to what
> > syntax sequence have to be that
>
> I let it to you to figure it out Leon.
>
> > (i thought just you know your
> > own code anyway better than me bein unexperienced in
MathLab).
> >
>
> I'm sure you can do it on your own. Don't be hesitate to
> play with the code.
>
> Start with simple example, for example:
>
> a = [770 105];
> b = 4830;
>
> There is three solutions c>=0, But only two c>0.
>
> Bruno

O.k I'll do it.
I'm trying to store the solvings in the separated txt file
and alter your syntax:

fid = fopen('C:.../file.txt', 'wt');
disp(c);
fprintf(fid, '\nSolutions c:\n\n');
But it doesn't work!

Leon

Bruno,
i tryed to see through for myself your function
x = intlin(a, b), but have had difficulties by understanding
it. (i'm not trained mathematican, unfortunately).
I intended to adjust it from x>=0 to x>0, but must say that
my knowledge in linear algebra yet lacking for it.
I would like to use your function for some purposes in my
fee work, but need it really with restrictions there x are
only not null. I would be lucky when you could adjust it for.
Kindly,
Leon