[Search for users] [Overall Top Noters] [List of all Conferences] [Download this site]

Conference rusure::math

Title:Mathematics at DEC
Moderator:RUSURE::EDP
Created:Mon Feb 03 1986
Last Modified:Fri Jun 06 1997
Last Successful Update:Fri Jun 06 1997
Number of topics:2083
Total number of notes:14613

926.0. "Nearest product" by LOCLE::RATCLIFF (What does "curiosity" mean?) Mon Sep 05 1988 12:21

I have the following problem and would like to find an algorithm to solve it:

  let K be an integer, and c_1, ..., c_n integers between 1 and K exclusive.
  Which subset of the c_is gives the largest product under K, when multiplying
  its elements?, e.g.:
   {2,3,5}, K=11; best subset is {2,5}.

  Typically, K ~ 3000, 5< n <20, 1< c_i < 300, but these limits are
  not absolute.

  An obvious algorithm is to try all 2^n-1 combinations, but I'm looking
  for a more elegant solution. Has anyone ideas/already done this?

  Thanks, John.
T.RTitleUserPersonal
Name
DateLines
926.1CLT::GILBERTspeciation sometimes convergesMon Sep 05 1988 15:2359
    By "largest product under K", I'm unsure whether you mean "the largest
    product <= K", or "the largest value when taken modulo K".  I'll assume
    the first of these.

>  An obvious algorithm is to try all 2^n-1 combinations, but I'm looking
>  for a more elegant solution. Has anyone ideas/already done this?

    I think this approach sounds quite workable, as follows:

1.  Set m <- 0.

2.  Set T <- { }; the empty set.

3.  Now T = { p | p is the product of some subset of { c_1 .. c_m },
	and p <= K }.

4.  Set m <- m + 1.
    If m > n, then we are done, and the answer is the largest value in T.

5.  Update T.  The new value of T is given by:
    newT <- oldT union { c_m * p | p is an element of oldT, and c_m * p <= K }

6.  Go to step 3.

The difficult part is providing an efficient implementation of step 5, and that 
requires a good representation for T.

Note that T never has more than K elements in it (assuming the c_i are all > 0).
Suppose T is kept as a linked list, with the elements in order of ascending
values (a priority queue would provide a better structure, but never mind that).

Suppose that the c_i are in sorted order, so that if i < j then c_i <= c_j.
(this seems useful, but is not strictly necessary).

Now step 5 can be performed by:

5a. Initialize.  Set U <- { }; the set U is also implemented as a linked list.
    Let x point to the first element of T.  If there are no elements of T,
    then we are done with step 5.

5b. Let v be the value pointed to by x.
    If c_m*v <= K, then insert c_m*v at the end of U.

5c. If U has no elements, goto step 5d.
    Compare the first element of U with v.
    If the first element of U < v, then:
	remove the first element of U; insert it before x; goto step 5c.
    If the first element of U = v, then:
	remove the first element of U; (discard it); goto step 5c.
    If the first element of U > v, then goto step 5d.

5d. If c_m*v <= K *and* there are more elements in T, then advance x to
    the next element, and goto step 5b.  Otherwise, proceed to step 5e.

5e. If U has no elements, we are done with step 5.
    Otherwise, remove the first element of U; insert it at the end of T,
    and goto step 5e.

I think this will give you very good performance.
926.2knapsack note is 377LISP::DERAMODaniel V. {AITG,LISP,ZFC}:: D&#039;EramoTue Sep 06 1988 01:2513
     You could start with the logarithms of the numbers and
     do additions instead of multiplications, although for
     the sizes you gave in .0 it's not really necessary. 
     And unless you are doing many of these at once the brute
     force method isn't too bad, either.
     
     You can split the n numbers into two lists of n/2, and store
     the 2^(n/2) products per list for all subsets of each list.
     Then sort one list ascending and the other descending, and
     run through them as I just described in "Solution 2" in
     reply 377.1 (knapsack problems).
     
     Dan
926.3CLT::GILBERTspeciation sometimes convergesTue Sep 06 1988 13:206
>    You can split the n numbers into two lists of n/2, and store
>    the 2^(n/2) products per list for all subsets of each list.

     Actually, you needn't store all 2 * 2^(n/2) products -- just
     those <= K.  In fact, because K is so small, it may be better
     to avoid the complication of splitting the list into two parts.
926.4N Log NDWOVAX::YOUNGfeet of clay, too.Tue Sep 06 1988 21:093
    I believe that this can be solved in O(n) or at worst O(n*Log(n)).
    
    --  Barry
926.5How about "close enough for government work"?POOL::HALLYBThe smart money was on GoliathMon Sep 12 1988 17:105
    Is it necessary to find the "best" subset, or merely a "good" one?
    (Assuming it speeds things up by an order of magnitude to find just
    a "good" solution).
    
      John
926.6LOCLE::RATCLIFFWhat does &quot;curiosity&quot; mean?Fri Sep 16 1988 08:494
    Re .5: ideally, the "best" would be better, but a "good" one would
    still be useful, e.g. 0.8*K < product < K.
    
    John.
926.7brute force is fast enough for this problemCLT::GILBERT$8,000,000,000 in damagesSat Sep 17 1988 21:00123
program x926 ( input, output ) ;
 
const
    maxk = 10000;	{ Largest K value we are concerned with handling }
    maxn = 100;		{ Maximum cardinality of C }

type
    c_array = array [ 0..maxn ] of integer;	{ element 0 holds the count }



function maximize_product (
    k : integer;	{ Final product must be <= k }
    var c : [ readonly ] c_array
    ) : integer;
type
    bitarray = packed array [ 1..maxk ] of boolean;
var
    i,j : integer;
    b : bitarray;
    b0: [ static, readonly ] bitarray := (repeat false);
    d : array [ 0..maxk ] of integer;
    highwater, newhigh : integer;
label
    10;
begin

{ We'll process all the elements of c in their cardinal order }

b := b0;
b[1] := true; highwater := 1;	{ only one product so far }
newhigh := 1;	{ tricky -- initialize this outside the following loop }

if k >= maxk then begin maximize_product := 0; goto 10; end;

for i := 1 to c[0] do
    begin
    d[0] := 0;	{ nothing new found for c[i] yet }
    { we only need to go up to the highest value that was set }
    { and we need never go beyond k/c[i], since j*c[i] would be > k }
    for j := 1 to min ( highwater, k div c[i] ) do
        begin
	if b[j] then if j*c[i] <= k then
	    begin
	    newhigh := j*c[i];	    { perhaps }
	    if not b[newhigh] then
		begin
		d[0] := d[0] + 1;
		d[d[0]] := newhigh;	{ save this product for now }
		end;
	    end;
	end;
    { now we include the new products from c[i] }
    for j := 1 to d[0] do
	begin
	b[d[j]] := true;
	end;
    highwater := max(highwater,newhigh);
    end;
maximize_product := highwater;
10:
end;
 


procedure write_factors (
    pp : integer;
    var c : c_array
    );
var
    c0 : integer;
    f,p,q : integer;
begin
c0 := c[0];	{ save this before we trash it }
p := pp;
write (p:0, ' = 1');

while c[0] > 0 do
    begin
    f := c[c[0]];	{ grab the potential factor we're removing }
    c[0] := c[0] - 1;
    if p <> maximize_product (p,c) then
	begin
	{ we just removed an important factor }
	write ('*', f:0);
	p := p div f;
	end;
    end;
writeln;
c[0] := c0;
end;

procedure main;
var
    k,p : integer;
    c : c_array;
begin
writeln ('Enter K followed by the C values (separated by spaces)');
while not eof(input) do
    begin
    write ('> ');
    read (k);
    c[0] := 0;
    while not eoln(input) and (c[0] < maxn) do
	begin
	c[0] := c[0] + 1;
	read (c[c[0]]);
	end;
    readln;
    p := maximize_product (k, c);
    writeln ('Maximized product is ', p:0);

    { having done that, let's write out the factors }
    { we do this in what may be the worst possible way }

    write_factors (p, c);
    
    end;
end;

begin
main;
end.
926.8brute force, huh?LISP::DERAMODaniel V. {AITG,LISP,ZFC}:: D&#039;EramoSat Sep 17 1988 23:528
     re .-1
     
>>     < Note 926.7 by CLT::GILBERT "$8,000,000,000 in damages" >
>>                -< brute force is fast enough for this problem >-

     Are you going to change your node::name to HURRICANE::GILBERT :-) ?
     
     Dan
926.9Thank you GilbertNSDC::GASCHENWed Sep 21 1988 12:51170
	Re 926.7  
    
    	Thank you for your x926 program. You are right it is fast enough
    for our needs. I have just written a second output procedure to
    display all possible combinations of c-i leading to the maximum
    product. 
    
    	(John has put note 926 in this conference for me during my leaves.)
    
program x926b ( input, output ) ;
 
const
    maxk = 10000;	{ Largest K value we are concerned with handling }
    maxn = 100;		{ Maximum cardinality of C }

type
    c_array = array [ 0..maxn ] of integer;	{ element 0 holds the count }

var
    d : array [ 0..maxk ] of integer;


function maximize_product (
    k : integer;	{ Final product must be <= k }
    var c : [ readonly ] c_array
    ) : integer;
type
    bitarray = packed array [ 1..maxk ] of boolean;
var
    i,j : integer;
    b : bitarray;
    b0: [ static, readonly ] bitarray := (repeat false);
    highwater, newhigh : integer;
label
    10;
begin

{ We'll process all the elements of c in their cardinal order }

b := b0;
b[1] := true; highwater := 1;	{ only one product so far }
newhigh := 1;	{ tricky -- initialize this outside the following loop }

if k >= maxk then begin maximize_product := 0; goto 10; end;

for i := 1 to c[0] do
    begin
    d[0] := 0;	{ nothing new found for c[i] yet }
    { we only need to go up to the highest value that was set }
    { and we need never go beyond k/c[i], since j*c[i] would be > k }
    for j := 1 to min ( highwater, k div c[i] ) do
        begin
	if b[j] then if j*c[i] <= k then
	    begin
	    newhigh := j*c[i];	    { perhaps }
	    if not b[newhigh] then
		begin
		d[0] := d[0] + 1;
		d[d[0]] := newhigh;	{ save this product for now }
		end;
	    end;
	end;
    { now we include the new products from c[i] }
    for j := 1 to d[0] do
	begin
	b[d[j]] := true;
	end;
    highwater := max(highwater,newhigh);
    end;
maximize_product := highwater;
10:
end;
 


procedure write_factors (
    pp : integer;
    var c : c_array
    );
var
    c0 : integer;
    f,p,q : integer;
begin
c0 := c[0];	{ save this before we trash it }
p := pp;
write (p:0, ' = 1');

while c[0] > 0 do
    begin
    f := c[c[0]];	{ grab the potential factor we're removing }
    c[0] := c[0] - 1;
    if p <> maximize_product (p,c) then
	begin
	{ we just removed an important factor }
	write ('*', f:0);
	p := p div f;
	end;
    end;
writeln;
c[0] := c0;
end;

procedure find_factors (
    p : integer;
    c : c_array;
    istart : integer
    );
var
    i,j,newstart,newp : integer;
    star : char;
begin
   for i := istart to c[0] do begin
      newp := p div c[i];
      if newp*c[i] = p then begin
         if newp > 1 then begin
            d[0] := d[0] + 1;
            d[d[0]] := i;
            newstart := i + 1;
            if newstart <= c[0] then begin
               find_factors(newp,c,newstart);
            end;
            d[0] := d[0] - 1;
         end else begin
            write('new solution :');
            star := ' ';
            d[0] := d[0] + 1;
            d[d[0]] := i;
            for j:=1 to d[0] do begin
               write(star,c[d[j]]:0,'(',d[j]:0,')');
               star := '*';
            end;
            writeln;
            d[0] := d[0] - 1;
         end;
      end;
   end;
end;

procedure main;
var
    k,p : integer;
    c : c_array;
begin
writeln ('Enter K followed by the C values (separated by spaces)');
while not eof(input) do
    begin
    write ('> ');
    read (k);
    c[0] := 0;
    while not eoln(input) and (c[0] < maxn) do
	begin
	c[0] := c[0] + 1;
	read (c[c[0]]);
	end;
    readln;
    p := maximize_product (k, c);
    writeln ('Maximized product is ', p:0);

    { having done that, let's write out the factors }
    { we do this in what may be the worst possible way }
    d[0] := 0;
    find_factors(p,c,1);
    write_factors (p, c);
    
    end;
end;

begin
main;
end.