Combinations

Basic Combination Generator

Generating combinations is easy. You just start at some value, then put the remaining values one by one. After that, you change the starting value and do the same thing over again. In my problem, I needed to generate combinations of bits, so this is the initial code:

void gencom(int N,int value, int start, int remaining)
{
  int i;
  if (remaining==0)
  {
    for(i=0;i<N;i++) if (value & (1<<i)) putc('a'+i,stdout);
    puts("");
    return ;
  }

  for(i=start;i+remaining<=N;i++)
  {
    value= value | (1<<i);
    gencom(N, value, i+1, remaining-1);
    value= value & ~(1<<i);
  }
}

void generate_combinations(int N,int M)
{
  gencom(N, 0, 0, M);
}

In my problem, a non-recursive function was more useful. I can figure out the non-recursive version of the above code but I want to resolve the recursion in a systematic way, just to develop the technique.

First of all, I'm going to jump into and out of the loop, so I better expand it into some labels and jumps:

void gencom(int N,int value, int start, int remaining)
{
  int i;
  if (remaining==0)
  {
    for(i=0;i<N;i++) if (value & (1<<i)) putc('a'+i,stdout);
    puts("");
    return ;
  }

  i= start;
loop_next:
  if (i+remaining>N) goto loop_end;

  value= value | (1<<i);
  gencom(N, value, i+1, remaining-1);
  value= value & ~(1<<i);
  
  i++;
  goto loop_next;
loop_end:
  return ;
}

Now, here is my algorithm to convert a recursive function to a non-recursive function:

1. Make new variable, depth, to represent the depth of the call stack for the recursive function.
2. Store all local variables (including arguments) in arrays, indexed by depth.
3. When there is a call to the recursive function: store the values in the next entries for the argument arrays, increase depth and restart the function (except for the initial setup code).
4. When there is a return from the function, if the return value is meaningful, store it in the variable ret to be used later. If the depth is zero, then make a real return. Otherwise, decrease depth and jump to the statement which follows the call.

Basicly, I'm doing the function calls manually, arranging both the stack and flow control by imitating the compiler. Naturally, the stack will be limited to a pre-determined depth. Fortunately for us, combinatorics problems don't need too much space since numbers used are usually quite small. Here is my first non-recursive version:

#define MAXDEPTH 64

void gencom(int N0,int value0, int start0, int remaining0)
{
  int N[MAXDEPTH], value[MAXDEPTH], 
      start[MAXDEPTH], remaining[MAXDEPTH];
  int i[MAXDEPTH];
  int depth;

  depth= 0;
  N[0]= N0;
  value[0]= value0;
  start[0]= start0;
  remaining[0]= remaining0;

restart:

  if (remaining[depth]==0)
  {
    int j;
    for(j=0;j<N[depth];j++) 
       if (value[depth] & (1<<j)) 
         putc('a'+j,stdout);
    puts("");
    goto ret;
  }

  i[depth]=start[depth];
loop_next:
  if (i[depth]+remaining[depth]>N[depth]) goto loop_end;

  value[depth]= value[depth] | (1<<i[depth]);
  goto call;
loop_cont:
  value[depth]= value[depth] & ~(1<<i[depth]);
  
  i[depth]++;
  goto loop_next;

ret:
  if (depth==0) return ;
  depth--;
  goto loop_cont;

loop_end:
  goto ret;

call:
  depth++;
  N[depth]= N[depth-1]; 
  value[depth]= value[depth-1];
  start[depth]= i[depth-1]+1;
  remaining[depth]= remaining[depth-1]-1;
  goto restart;
}

Now it's time to make some optimizations. First of all, N[x] has only one value: N0. N[0] is initialized to N0, and other entries are just copied from there at the call. So we can replace N[x] simply by the integer argument N, as it was in the original function.

Second, start[depth] is only used to initialize i[depth]. There is only one reference to it. Also, there is only one assignment to it (during the call). We can therefore skip this temporary variable and initialize i[depth] at the time of the call and the entry initialization.

void gencom(int N,int value0, int start0, int remaining0)
{
  int value[MAXDEPTH], remaining[MAXDEPTH];
  int i[MAXDEPTH];
  int depth;

  depth= 0;
  value[0]= value0;
  i[0]= start0;
  remaining[0]= remaining0;

restart:

  if (remaining[depth]==0)
  {
    int j;
    for(j=0;j<N;j++) if (value[depth] & (1<<j)) putc('a'+j,stdout);
    puts("");
    goto ret;
  }

loop_next:
  if (i[depth]+remaining[depth]>N) goto loop_end;

  value[depth]= value[depth] | (1<<i[depth]);
  goto call;
loop_cont:
  value[depth]= value[depth] & ~(1<<i[depth]);
  
  i[depth]++;
  goto loop_next;

ret:
  if (depth==0) return ;
  depth--;
  goto loop_cont;
loop_end:
  goto ret;

call:
  depth++;
  value[depth]= value[depth-1];
  i[depth]= i[depth-1]+1;
  remaining[depth]= remaining[depth-1]-1;
  goto restart;
}

Now we see that remaining[] and depth are interrelated. Whenever depth gets incremented, remaining[depth] is initialized to one less remaining[depth-1]. This establishes the following pattern for these two variables:

    depth     remaining[depth]
    --------------------------
      0       remaining0-0
      1       remaining0-1
      2       remaining0-2
      3       remaining0-3
      ...

Since there is no other assigment to remaining[depth], we can replace this expression with its value: remaining0-depth. This also eliminates the remaining[] array.

void gencom(int N,int value0, int start0, int remaining0)
{
  int value[MAXDEPTH];
  int i[MAXDEPTH];
  int depth;

  depth= 0;
  value[0]= value0;
  i[0]= start0;

restart:

  if ((remaining0-depth)==0)
  {
    int j;
    for(j=0;j<N;j++) if (value[depth] & (1<<j)) putc('a'+j,stdout);
    puts("");
    goto ret;
  }
loop_next:
  if (i[depth]+(remaining0-depth)>N) goto loop_end;

  value[depth]= value[depth] | (1<<i[depth]);
  goto call;
loop_cont:
  value[depth]= value[depth] & ~(1<<i[depth]);
  
  i[depth]++;
  goto loop_next;

ret:
  if (depth==0) return ;
  depth--;
  goto loop_cont;
loop_end:
  goto ret;

call:
  depth++;
  value[depth]= value[depth-1];
  i[depth]= i[depth-1]+1;
 goto restart;
}

Now I will make another subtle optimization. Each time the loop finishes, value[depth] has the same value as value[depth-1]. This is because the loop body always clears the bits it has set. This means that we could use just one variable instead of the value array. Sure, the sub-function will modify the variable of its parent, but it's going to restore it back. No harm done.

void gencom(int N,int value0, int start0, int remaining0)
{
  int value;
  int i[MAXDEPTH];
  int depth;

  depth= 0;
  value= value0;
  i[0]= start0;

restart:

  if ((remaining0-depth)==0)
  {
    int j;
    for(j=0;j<N;j++) if (value & (1<<j)) putc('a'+j,stdout);
    puts("");
    goto ret;
  }

loop_next:
  if (i[depth]+(remaining0-depth)>N) goto loop_end;
  value= value | (1<<i[depth]);
  goto call;
loop_cont:
  value= value & ~(1<<i[depth]);
  
  i[depth]++;
  goto loop_next;

ret:
  if (depth==0) return ;
  depth--;
  goto loop_cont;
loop_end:
  goto ret;

call:
  depth++;
  i[depth]= i[depth-1]+1;
  goto restart;
}

I don't see how I can get rid of the array i[]. That has to stay there, I guess. Anyway, we can inline loop_end and call labels since they are quite simple now and are only called from one place each. While we're doing that, we can also inline the whole function into generate_combinations.

Recall that the original call from generate_combinations was:

 gencom(N, 0, 0, M);

This means that:

• N(gencom)= N(generate_combinations)
• value0 = 0
• start0 = 0
• remaining0 = M

So the final thing looks like this. I looks quite complicated, but it was all generated from a simple recursive function:

void generate_combinations(int N,int M)
{
  int value;
  int i[MAXDEPTH];
  int depth;

  depth= 0;
  value= 0;
  i[0]= 0;

restart:
  if ((M-depth)==0)
  {
    int j;
    for(j=0;j<N;j++) if (value & (1<<j)) putc('a'+j,stdout);
    puts("");
    goto ret;
  }

loop_next:
  if (i[depth]+(M-depth)>N) goto ret;
  value= value | (1<<i[depth]);
  depth++;
  i[depth]= i[depth-1]+1;
  goto restart;

loop_cont:
  value= value & ~(1<<i[depth]);
  i[depth]++;
  goto loop_next;

ret:
  if (depth==0) return ;
  depth--;
  goto loop_cont;
}

Generating Some Subsets

A problem related to combinations is the generation of subsets. Here, I'm interested in all subsets with a maximum element count.

The idea is simple. I always try to add one element to the current subset. In order to avoid duplicates due to ordering, I never add an element which is smaller than the last one in the current subset. Here, A < B means that A occurs before B in the input array.

This operation will eventually fail, due to two reasons. First, we might run out of elements. Second, we might have filled all the slots in the subset. In both cases we do the same thing: remove the last element and try to add a new one.

#include <stdio.h>
#include <stdlib.h>

void prt_all_com(int N, int M)
{
  char R[200];  // auxilliary array for printing
  int idx[200]; // the subset 
  int P;        // current number of elements
  int S;        // next element to be put

  P= 0;
  R[0]= 0;
found:
  R[P]= 0;
  printf("C: %s\n", R);

  if (P==M) goto remove;
  if (P==0) S= 0;
       else S= idx[P-1]+1;

add:
  if (S==N) goto remove;
  P++;
  idx[P-1]= S;
  R[P-1]= 'A' + idx[P-1]; 
  goto found;

remove:
  S= idx[P-1] + 1;
  P--;
  if (P==-1) return ;
  goto add;
}

int main(int argc,char **argv)
{
  prt_all_com(atoi(argv[1]), atoi(argv[2]));
  return 0;
}

A lower bound for the number of elements can be implemented very easily. Just check that the number of elements is acceptable before printing or processing the array. Here is the order traversed:

$ ./a.out 4 3 
C: 
C: A
C: AB
C: ABC
C: ABD
C: AC
C: ACD
C: AD
C: B
C: BC
C: BCD
C: BD
C: C
C: CD
C: D

Now that I look at it again, this is the same thing as the combination generator. The only difference is that, we now print the array even when the number of elements is less. In fact, both functions are variations on an all-subsets-generator. We simply cut off the search when we exceed the number of elements required. Here is the generator which does everything:

#include <stdio.h>
#include <stdlib.h>

                        // max number of elements
                              // min number of elements
void prt_all_com(int N, int M,int Z)
{
  char R[200];  // auxilliary array for printing
  int idx[200]; // the subset 
  int P;        // current number of elements
  int S;        // next element to be put

  P= 0;
  R[0]= 0;
found:
  R[P]= 0;
  if (P>=Z)
    printf("C: %s\n", R);

  if (P==M) goto remove;
  if (P==0) S= 0;
       else S= idx[P-1]+1;

add:
  if (S==N) goto remove;
  P++;
  idx[P-1]= S;
  R[P-1]= 'A' + idx[P-1]; 
  goto found;

remove:
  S= idx[P-1] + 1;
  P--;
  if (P==-1) return ;
  goto add;
}

int main(int argc,char **argv)
{
  prt_all_com(atoi(argv[1]), atoi(argv[2]),atoi(argv[3]));
  return 0;
}

For generating combinations with an exact number of elements, the above function takes more time than the purpose-built function when the prefix isn't suitable for further progress. For example if the problem is N=5 and M=3, then element 3 can't be put into position 0: we wouldn't be able to put anything to position 2.

The purpose-built function discovers this by simply checking the numbers. The subset generator has to proceed with that and fail later even if M=Z=3. This is because the subset function is designed to process shorter sequences as well.

I can do the same check here as well. I just need to modify the condition S==N. Instead of checking that we have reached the end of the available elements, we just check the numbers and see whether starting the subsequence at this element could yield a valid subset or not. Just like the combination generator.

When we have P elements, we need to add (Z-P) elements to make a subset of size Z or greater. If we add this required number of elements to the starting element S, we mustn't exceed N.

If you add the following under S==N check, this function will run as efficient as the purpose built function for M=Z.

  if (P<Z && S+Z-P>N) goto remove; 

We can safely remove condition P<Z from above. S starts out from value 0 and increases by one now and then, until it reaches N. When that happens, the previous condition S==N prevents the code from increasing S further and cuts off the search. So, we know that S<N at this point. If P≥Z, then the expression S+Z-P can never be greater than N since S is smaller than N and Z-P is non-positive. Therefore checking for P<Z is useless at this point. When combined with the previous check, the cut-off condition becomes:

  if (S==N || S+Z-P>N) goto remove;

Here is the final version:

#include <stdio.h>
#include <stdlib.h>

                        // max number of elements
                               // min number of elements
void prt_all_com(int N, int M, int Z)
{
  char R[200];  // auxilliary array for printing
  int idx[200]; // the subset 
  int P;        // current number of elements
  int S;        // next element to be put

  P= 0;
found:
  R[P]= 0;
  if (P>=Z)
    printf("C: %s\n", R);

  if (P==M) goto remove;
  if (P==0) S= 0;
       else S= idx[P-1]+1;

add:
  if (S==N || S+Z-P>N) goto remove;
  P++;
  idx[P-1]= S;
  R[P-1]= 'A' + idx[P-1]; 
  goto found;

remove:
  if (P==0) return ;
  S= idx[P-1] + 1;
  P--;
  goto add;
}

int main(int argc,char **argv)
{
  prt_all_com(atoi(argv[1]), atoi(argv[2]),atoi(argv[3]));
  return 0;
}

This is so nice. The variables idx and P are necessary for storing output. R is only there to print out the thing. The only extra variable is S. This function could be further optimized, but I think it's quite easy to understand in this form so I don't want to mess with it.

Binary counting can generate all subsets as well. It's also much easier to understand. However, it can't limit the size of the output. That takes extra work and is not worth pursue.

OK, here's a version with no gotos. It looks nicer but is harder to understand. The output is in idx[1..P].

#include <stdio.h>
#include <stdlib.h>

                        // max number of elements
                               // min number of elements
void prt_all_com(int N, int M, int Z)
{
  char R[200];  // auxilliary array for printing
  int idx[200]; // the subset 
  int P;        // current number of elements
  int S;        // next element to be put

  P= 0;
  idx[0]= -1;
  while(1)
  {
    if (P>=Z) { R[P]= 0; printf("C: %s\n", R); }
    S= idx[P]+1;
    while(P==M || S==N || S+Z-P>N)
    {
      if (P==0) return ;
      S= idx[P] + 1;
      P--;
    }
    P++;
    idx[P]= S;
    R[P-1]= 'A' + idx[P];
  }
}

int main(int argc,char **argv)
{
  prt_all_com(atoi(argv[1]), atoi(argv[2]),atoi(argv[3]));
  return 0;
}

Just modify the hilighted areas to suit your data structures.