The fastest combination generator

[dirt_bagz]

Dear all,

I am looking for an expert opinion on what the fastest combination generator is. I am currently facing a problem of generating a huge combinations (k elements out of n) and as it stands now, the running time of my program will take more than three years. My current implementation uses the Algorithm R (see Knuth's book pre-fascicle 3a), a revolving-door algorithm.

If anyone knows a better/faster way of enumerating combinations, please let me know. Many thanks.

[iamthwee]

Dear all,

I am looking for an expert opinion on what the fastest combination generator is. I am currently facing a problem of generating a huge combinations (k elements out of n) and as it stands now, the running time of my program will take more than three years. My current implementation uses the Algorithm R (see Knuth's book pre-fascicle 3a), a revolving-door algorithm.

If anyone knows a better/faster way of enumerating combinations, please let me know. Many thanks.

Well unless you have access to a quantum computer I doubt you will be able to cut down on run time. Knuth's work is generally accepted as being the best for algorithms and reducing time complexity.

However, you may be looking at your problem in the wrong way. Take for example the classic travelling salesman problem.

Applying a brute force algorithm to calculate the shortest route for a journey with 100+ pit-stops will take an unfathomable amount of time.

However, if you adopt a heuristic approach you can find a high probabilistic solution,within a much smaller time scale.

Perhaps you need to think about approaching your problem from a different angle.

See this:

http://en.wikipedia.org/wiki/Dynamic_programming

Incase anyone wanted to know here is Knuth's algorithm for producing combinations in C.

/* Algorithm by Donald Knuth. */
 
#include <stdio.h>
#include <stdlib.h> 
void main( void)
{
	int i, j=1, k, n, *c, x;
	printf( "Enter n,k: ");
	scanf( "%d,%d", &n, &k);
	c = malloc( (k+3) * sizeof(int));
 
	for (i=1; i <= k; i++) c[i] = i;
	c[k+1] = n+1;
	c[k+2] = 0;
	j = k;
visit:
	for (i=k; i >= 1; i--) printf( "%3d", c[i]);
	printf( "\n");
	if (j > 0) {x = j+1; goto incr;}
	if (c[1] + 1 < c[2])
	 {
	 c[1] += 1;
	 goto visit;
	 }
	j = 2;
do_more:
	c[j-1] = j-1;
	x = c[j] + 1;
	if (x == c[j+1]) {j++; goto do_more;}
	if (j > k) exit(0);
incr:
	c[j] = x;
	j--;
	goto visit;
}

And java

//CombinationGenerator.java
//--------------------------------------
// Systematically generate combinations.
//--------------------------------------
import java.math.BigInteger;
public class CombinationGenerator {
private int[] a;
private int n;
private int r;
private BigInteger numLeft;
private BigInteger total;
//------------
// Constructor
//------------
public CombinationGenerator (int n, int r) {
	if (r > n) {
	 throw new IllegalArgumentException ();
	}
	if (n < 1) {
	 throw new IllegalArgumentException ();
	}
	this.n = n;
	this.r = r;
	a = new int[r];
	BigInteger nFact = getFactorial (n);
	BigInteger rFact = getFactorial (r);
	BigInteger nminusrFact = getFactorial (n - r);
	total = nFact.divide (rFact.multiply (nminusrFact));
	reset ();
}
//------
// Reset
//------
public void reset () {
	for (int i = 0; i < a.length; i++) {
	 a[i] = i;
	}
	numLeft = new BigInteger (total.toString ());
}
//------------------------------------------------
// Return number of combinations not yet generated
//------------------------------------------------
public BigInteger getNumLeft () {
	return numLeft;
}
//-----------------------------
// Are there more combinations?
//-----------------------------
public boolean hasMore () {
	return numLeft.compareTo (BigInteger.ZERO) == 1;
}
//------------------------------------
// Return total number of combinations
//------------------------------------
public BigInteger getTotal () {
	return total;
}
//------------------
// Compute factorial
//------------------
private static BigInteger getFactorial (int n) {
	BigInteger fact = BigInteger.ONE;
	for (int i = n; i > 1; i--) {
	 fact = fact.multiply (new BigInteger (Integer.toString (i)));
	}
	return fact;
}
//--------------------------------------------------------
// Generate next combination (algorithm from Rosen p. 286)
//--------------------------------------------------------
public int[] getNext () {
	if (numLeft.equals (total)) {
	 numLeft = numLeft.subtract (BigInteger.ONE);
	 return a;
	}
	int i = r - 1;
	while (a[i] == n - r + i) {
	 i--;
	}
	a[i] = a[i] + 1;
	for (int j = i + 1; j < r; j++) {
	 a[j] = a[i] + j - i;
	}
	numLeft = numLeft.subtract (BigInteger.ONE);
	return a;
}
}

/*
* test.java
*
* Created on 14 January 2006, 16:29
*/
/**
*
* @author iamthwee
*/
public class test {
 
	/** Creates a new instance of test */
	public test() {
	}
 
	/**
	 * @param args the command line arguments
	 */
	public static void main(String[] args) {
		String[] elements = {"a", "b","c","d"};
int[] indices;
CombinationGenerator x = new CombinationGenerator (elements.length, 2);
StringBuffer combination;
while (x.hasMore ()) {
combination = new StringBuffer ();
indices = x.getNext ();
for (int i = 0; i < indices.length; i++) {
	combination.append (elements[indices[i]]);
}
System.out.println (combination.toString ());
}
 
	}
 
}

[dirt_bagz]

Well unless you have access to a quantum computer I doubt you will be able to cut down on run time. Knuth's work is generally accepted as being the best for algorithms and reducing time complexity.

However, you may be looking at your problem in the wrong way. Take for example the classic travelling salesman problem.

Applying a brute force algorithm to calculate the shortest route for a journey with 100+ pit-stops will take an unfathomable amount of time.

However, if you adopt a heuristic approach you can find a high probabilistic solution,within a much smaller time scale.

Perhaps you need to think about approaching your problem from a different angle.

See this:

http://en.wikipedia.org/wiki/Dynamic_programming

Incase anyone wanted to know here is Knuth's algorithm for producing combinations in C.

/* Algorithm by Donald Knuth. */
 
#include <stdio.h>
#include <stdlib.h> 
void main( void)
{
	int i, j=1, k, n, *c, x;
	printf( "Enter n,k: ");
	scanf( "%d,%d", &n, &k);
	c = malloc( (k+3) * sizeof(int));
 
	for (i=1; i <= k; i++) c[i] = i;
	c[k+1] = n+1;
	c[k+2] = 0;
	j = k;
visit:
	for (i=k; i >= 1; i--) printf( "%3d", c[i]);
	printf( "\n");
	if (j > 0) {x = j+1; goto incr;}
	if (c[1] + 1 < c[2])
	 {
	 c[1] += 1;
	 goto visit;
	 }
	j = 2;
do_more:
	c[j-1] = j-1;
	x = c[j] + 1;
	if (x == c[j+1]) {j++; goto do_more;}
	if (j > k) exit(0);
incr:
	c[j] = x;
	j--;
	goto visit;
}

And java

//CombinationGenerator.java
//--------------------------------------
// Systematically generate combinations.
//--------------------------------------
import java.math.BigInteger;
public class CombinationGenerator {
private int[] a;
private int n;
private int r;
private BigInteger numLeft;
private BigInteger total;
//------------
// Constructor
//------------
public CombinationGenerator (int n, int r) {
	if (r > n) {
	 throw new IllegalArgumentException ();
	}
	if (n < 1) {
	 throw new IllegalArgumentException ();
	}
	this.n = n;
	this.r = r;
	a = new int[r];
	BigInteger nFact = getFactorial (n);
	BigInteger rFact = getFactorial (r);
	BigInteger nminusrFact = getFactorial (n - r);
	total = nFact.divide (rFact.multiply (nminusrFact));
	reset ();
}
//------
// Reset
//------
public void reset () {
	for (int i = 0; i < a.length; i++) {
	 a[i] = i;
	}
	numLeft = new BigInteger (total.toString ());
}
//------------------------------------------------
// Return number of combinations not yet generated
//------------------------------------------------
public BigInteger getNumLeft () {
	return numLeft;
}
//-----------------------------
// Are there more combinations?
//-----------------------------
public boolean hasMore () {
	return numLeft.compareTo (BigInteger.ZERO) == 1;
}
//------------------------------------
// Return total number of combinations
//------------------------------------
public BigInteger getTotal () {
	return total;
}
//------------------
// Compute factorial
//------------------
private static BigInteger getFactorial (int n) {
	BigInteger fact = BigInteger.ONE;
	for (int i = n; i > 1; i--) {
	 fact = fact.multiply (new BigInteger (Integer.toString (i)));
	}
	return fact;
}
//--------------------------------------------------------
// Generate next combination (algorithm from Rosen p. 286)
//--------------------------------------------------------
public int[] getNext () {
	if (numLeft.equals (total)) {
	 numLeft = numLeft.subtract (BigInteger.ONE);
	 return a;
	}
	int i = r - 1;
	while (a[i] == n - r + i) {
	 i--;
	}
	a[i] = a[i] + 1;
	for (int j = i + 1; j < r; j++) {
	 a[j] = a[i] + j - i;
	}
	numLeft = numLeft.subtract (BigInteger.ONE);
	return a;
}
}

/*
* test.java
*
* Created on 14 January 2006, 16:29
*/
/**
*
* @author iamthwee
*/
public class test {
 
	/** Creates a new instance of test */
	public test() {
	}
 
	/**
	 * @param args the command line arguments
	 */
	public static void main(String[] args) {
		String[] elements = {"a", "b","c","d"};
int[] indices;
CombinationGenerator x = new CombinationGenerator (elements.length, 2);
StringBuffer combination;
while (x.hasMore ()) {
combination = new StringBuffer ();
indices = x.getNext ();
for (int i = 0; i < indices.length; i++) {
	combination.append (elements[indices[i]]);
}
System.out.println (combination.toString ());
}
 
	}
 
}

Dear iamthwee, thanks for your response. I understand that there is a probabilistic way of looking at this problem. However, for my this particular problem, probabilistic solution is not acceptable. This problem is similar to the search for optimum Golomb Ruler of certain marks or minimum distance computation in linear error correcting codes. As you can see here, in order to prove something, a non probabilistic solution is required. Anyway, many thanks.

[iamthwee]

>This problem is similar to the search for optimum Golomb Ruler Well finding a solution AND proving it is quite beyond me. Oh dear brain lock...:cry: Sorrie.

[Lord Soth]

Hi,

Do you need a solution which is in a smaller time complexity class than Knuth's or have a bound for n,k and just need a faster code to get it done in less than 3 years ? +what is the maximum acceptable time span for you ? because I have some ideas for optimization which can reduce your 3 year but not for n limit to infinity.

Loren Soth

[dirt_bagz]

Hi Loren Soth,

Thanks for your response. I just need a faster code to get my computation done. The value of n does not go to infinity, in fact n <= 512 and 1 <= k <= n/2. If possible, I would like to have to computation done in less than 1 year. Well, 1 year may seem rather long, but I probably can run this computation in parallel using many computers. Hope to hear your optimisation ideas soon.

dirt_bagz.

Hi,

Do you need a solution which is in a smaller time complexity class than Knuth's or have a bound for n,k and just need a faster code to get it done in less than 3 years ? +what is the maximum acceptable time span for you ? because I have some ideas for optimization which can reduce your 3 year but not for n limit to infinity.

Loren Soth

[sandysat]

Hi there!...I'd like to point you to http://www.merriampark.com/comb.htm
I don't know if it's any faster than your current considerations but hope it'll open up some more avenues for thinking...
rgds,
Sandeep

[piagina]

I doubt it is as fast as the Knuth code but here is a slightly more modern version of a combination generator using some features of C++ so I guess its C++

nb. if you copy the code and run it you may like to change the lines that set j and k as this currently give all possible combinations of the letters of the English alphabet.

#include <iostream>
#include <vector>
#include <string>

using namespace std;

// just to check the numbers are right

long totalcombs(int n, int r){
	long c=1;
	if (r > n) return 0;
	for (long d=1; d <= r; d++) {
		c *= n--;
		c /= d;
	}
	return c;
}

// example code to get all combination 

int main(){
	vector<int> indx;
	string alpha("ABCDEFGHIJKLMNOPQRSTUVWXYZ");
	int n=26;
	int j=1;
	int k=26;
	for(int twk=j;twk<=k;twk++){
		int r=twk;
		int total=totalcombs(n,r);
		int ccount=1;
		bool done=true;
		for(int iwk=0;iwk<r;iwk++)indx.push_back(iwk);
		while(done){
			done=false;
			cout << ccount++ << " of " << total << " ";
			for(int owk=0;owk<r;owk++){
				cout << alpha[indx[owk]] << " ";
			}
			cout << endl;
			for(int iwk=r-1;iwk>=0;iwk--){
				if(indx[iwk]<=(n-1)-(r-iwk)){
					indx[iwk]++;
					for(int swk=iwk+1;swk<r;swk++){
						indx[swk]=indx[swk-1]+1;
					}
					iwk=-1;
					done=true;
				}	
			}	
		}
		cout << " --------------------------- " << endl;
		indx.clear();
	}
	return 0;
}

[wwx]

Does the problem resolved ? Can you describe what you're doing ? Generate C(N,1..N/2) ?

[ithelp]

It is an NP complete problem , so you can only get some approximation algorithm to get it faster.

[androidandrew]

Hi there,

I have done a bit of reading on this topic and I know there are some research papers on the topic of multiprocessor enumeration optimisation.

for example one is:
'A Parallel Algorithm for Enumerating Combinations', By Martha Torres, Alfredo Goldman,
Junior Barrera

good luck.

[androidandrew]

oh sorry i forgot to mention that was in

'Proceedings of the 2003 International Conference on Parallel Processing'

i think IEEE Computer Society. I found it using my uni journal sources so if you know someone with access to a similar source then you should be able to find it.

[wwx]

Here is a link for the source if someone interesting in paralel processing

Link

[shaunwilliams]

Here's an impressive source of many different combination algorithms. They're documented in a book called "Algorithms for Programmers" by Jörg Arndt. It will be published later this year, but as of now, you can download the PDF.
http://www.jjj.de/fxt/#fxtbook

I will certainly buy this book when it's published!

[sarehu]

It is an NP complete problem , so you can only get some approximation algorithm to get it faster.

No it's not. Learn what NP complete means. Learn what NP means. This problem isn't even in NP.

[bhavanaindia]

the code is not correct. If u enter n= 5 and k=5 it should generate only one combination
but generating lot many.
nCr = !n/(!r*!(n-r))

Well unless you have access to a quantum computer I doubt you will be able to cut down on run time. Knuth's work is generally accepted as being the best for algorithms and reducing time complexity.

However, you may be looking at your problem in the wrong way. Take for example the classic travelling salesman problem.

Applying a brute force algorithm to calculate the shortest route for a journey with 100+ pit-stops will take an unfathomable amount of time.

However, if you adopt a heuristic approach you can find a high probabilistic solution,within a much smaller time scale.

Perhaps you need to think about approaching your problem from a different angle.

See this:

http://en.wikipedia.org/wiki/Dynamic_programming

Incase anyone wanted to know here is Knuth's algorithm for producing combinations in C.

/* Algorithm by Donald Knuth. */
 
#include <stdio.h>
#include <stdlib.h> 
void main( void)
{
	int i, j=1, k, n, *c, x;
	printf( "Enter n,k: ");
	scanf( "%d,%d", &n, &k);
	c = malloc( (k+3) * sizeof(int));
 
	for (i=1; i <= k; i++) c[i] = i;
	c[k+1] = n+1;
	c[k+2] = 0;
	j = k;
visit:
	for (i=k; i >= 1; i--) printf( "%3d", c[i]);
	printf( "\n");
	if (j > 0) {x = j+1; goto incr;}
	if (c[1] + 1 < c[2])
	 {
	 c[1] += 1;
	 goto visit;
	 }
	j = 2;
do_more:
	c[j-1] = j-1;
	x = c[j] + 1;
	if (x == c[j+1]) {j++; goto do_more;}
	if (j > k) exit(0);
incr:
	c[j] = x;
	j--;
	goto visit;
}

And java

//CombinationGenerator.java
//--------------------------------------
// Systematically generate combinations.
//--------------------------------------
import java.math.BigInteger;
public class CombinationGenerator {
private int[] a;
private int n;
private int r;
private BigInteger numLeft;
private BigInteger total;
//------------
// Constructor
//------------
public CombinationGenerator (int n, int r) {
	if (r > n) {
	 throw new IllegalArgumentException ();
	}
	if (n < 1) {
	 throw new IllegalArgumentException ();
	}
	this.n = n;
	this.r = r;
	a = new int[r];
	BigInteger nFact = getFactorial (n);
	BigInteger rFact = getFactorial (r);
	BigInteger nminusrFact = getFactorial (n - r);
	total = nFact.divide (rFact.multiply (nminusrFact));
	reset ();
}
//------
// Reset
//------
public void reset () {
	for (int i = 0; i < a.length; i++) {
	 a[i] = i;
	}
	numLeft = new BigInteger (total.toString ());
}
//------------------------------------------------
// Return number of combinations not yet generated
//------------------------------------------------
public BigInteger getNumLeft () {
	return numLeft;
}
//-----------------------------
// Are there more combinations?
//-----------------------------
public boolean hasMore () {
	return numLeft.compareTo (BigInteger.ZERO) == 1;
}
//------------------------------------
// Return total number of combinations
//------------------------------------
public BigInteger getTotal () {
	return total;
}
//------------------
// Compute factorial
//------------------
private static BigInteger getFactorial (int n) {
	BigInteger fact = BigInteger.ONE;
	for (int i = n; i > 1; i--) {
	 fact = fact.multiply (new BigInteger (Integer.toString (i)));
	}
	return fact;
}
//--------------------------------------------------------
// Generate next combination (algorithm from Rosen p. 286)
//--------------------------------------------------------
public int[] getNext () {
	if (numLeft.equals (total)) {
	 numLeft = numLeft.subtract (BigInteger.ONE);
	 return a;
	}
	int i = r - 1;
	while (a[i] == n - r + i) {
	 i--;
	}
	a[i] = a[i] + 1;
	for (int j = i + 1; j < r; j++) {
	 a[j] = a[i] + j - i;
	}
	numLeft = numLeft.subtract (BigInteger.ONE);
	return a;
}
}

/*
* test.java
*
* Created on 14 January 2006, 16:29
*/
/**
*
* @author iamthwee
*/
public class test {
 
	/** Creates a new instance of test */
	public test() {
	}
 
	/**
	 * @param args the command line arguments
	 */
	public static void main(String[] args) {
		String[] elements = {"a", "b","c","d"};
int[] indices;
CombinationGenerator x = new CombinationGenerator (elements.length, 2);
StringBuffer combination;
while (x.hasMore ()) {
combination = new StringBuffer ();
indices = x.getNext ();
for (int i = 0; i < indices.length; i++) {
	combination.append (elements[indices[i]]);
}
System.out.println (combination.toString ());
}
 
	}
 
}

[henrik14]

Honored sir,
I am intersted in Your problem :-)) Long ago I have solved >>Traveling salesman<< problem when I realized that contrary to overall opinion there is only one optimal path between cities, and difference is only in starting and ending city. Implementing one kind of fractall geometry (Sierpinski Curve) as basis for adressing cities, I was able to get 1000 cities solution under one hour on ATARI TT which worked on 60 Mhz.
Now, my question is, do You just need all combinations in order and eventually listed on file?
But, I suppose You have some good use for this combinations, so if You explain exactly what it is supposed to be, perhaps I could find solution in a way You dont expect.
Last of all, if my solution satisfies Your requirements and You can get results fast, how much it would be worth to You?
Sorry that I have to ask this, but I am poor man from poor country and I have to fight for survival just to stay alive (I am retired)....
I am also married to one woman from India, who would like to visit her sick father who need heart operation and pacemaker, and I cannot buy even airoplane tickets..........
Regards from Zagreb, the capitol of Croatia, Europe!
Marijan Pollak, IT SA/SE 1st. Class, Instructor & Team Leader

[overtomanu]

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

int nCk(int n,int loopno,int ini,int *a,int k)
{
	static int count=0;
	int i;
	loopno--;
	if(loopno<0)
	{
		a[k-1]=ini;
		for(i=0;i<k;i++)
		{
			printf("%d,",a[i]);
		}
		printf("\n");
		count++;
		return 0;
	}
	for(i=ini;i<=n-loopno-1;i++)
	{
		a[k-1-loopno]=i+1;
		nCk(n,loopno,i+1,a,k);
	}
	if(ini==0)
	return count;
	else
	return 0;
}

void main()
{
	int n,k,*a,count;
	printf("Enter the value of n and k\n");
	scanf("%d %d",&n,&k);
	a=(int*)malloc(k*sizeof(int));
	count=nCk(n,k,0,a,k);
	printf("No of combinations=%d\n",count);
}

[overtomanu]

Here is a recursive c program thats generates all of the combinations nCk in a fast possible way.
This code works for nC0 also.

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

int nCk(int n,int loopno,int ini,int *a,int k)
{
	if(k==0)
	return 1;
	static int count=0;
	int i;
	loopno--;
	if(loopno<0)
	{
		a[k-1]=ini;
		for(i=0;i<k;i++)
		{
			printf("%d,",a[i]);
		}
		printf("\n");
		count++;
		return 0;
	}
	for(i=ini;i<=n-loopno-1;i++)
	{
		a[k-1-loopno]=i+1;
		nCk(n,loopno,i+1,a,k);
	}
	if(ini==0)
	return count;
	else
	return 0;
}

void main()
{
	int n,k,*a,count;
	printf("Enter the value of n and k\n");
	scanf("%d %d",&n,&k);
	a=(int*)malloc(k*sizeof(int));
	count=nCk(n,k,0,a,k);
	printf("No of combinations=%d\n",count);
	getch();
}

[babak_behravesh]

Jörg Arndt the author of Matters Computational, published by Springer in 2011, addressed an approach using Co-lexicographical order to generate combinations. In this method a coposition set is generated, then based on composition set the algorithm generates all combinatins.
1) Generating Composition Set:
Given n and k, the algorithm generates composition set, say n=5 and k =3.
The length of each set in composition set is calculated as, n -k +1. So for the given n and k, we have 3 which is the size of each composition. The first item of first composition is k as bellow,

3 . .

The next round the algorithm creates the next composition from current one by subtracting first non-zero item from left side, placing the result in left most place; and increasing the next place by 1. The process will be stopped when r locates in left most possition.
So, composition set for (n k) is presented as bellow,

3 . .
2 1 .
1 2 .
. 3 .
2 . 1
1 1 1
. 2 1
. 1 2
. . 3

Then, each row of composition set transforms into its respective combination. In first row 3 means, 3 conseuitive numbers followed by nothing, so the first combination is, 1 2 3
The second comopsition generates two consecutive numbers following a gap in numbers and then a number which creates, 1 2 4. The entire set is as bellow,
1 2 3
1 2 4
2 3 4
1 2 5
1 3 5
2 3 5
2 4 5
3 4 5

#ifndef KNUTH_H_
#define KNUTH_H_
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
using namespace std;

class composition_colex
{
    public:
    unsigned long n_, k_; // composition of n into k parts
    unsigned long nk1_;
    unsigned long *x_;
    unsigned long *b_;

    // data (k elements)

    void first()
    {
        x_[0] = k_; // all in first position
        for (unsigned long k=1; k<nk1_; ++k) x_[k] = 0;

    }

    void last()
    {
        for (unsigned long k=0; k<nk1_; ++k) x_[k] = 0;
        x_[k_-1] = k_; // all in last position
    }

    unsigned long next()
    // Return position of rightmost change, return k with last composition.
    {
        unsigned long j = 0;
        while ( 0==x_[j] ) ++j; // find first nonzero

        if ( j==k_-1 && x_[nk1_-1]==k_)            // current composition is last
            return k_;

        unsigned long v = x_[j];// value of first nonzero
        x_[j] = 0;                // set to zero
        x_[0] = v - 1;            // value-1 to first position
        ++j;
        ++x_[j];                // increment next position
        return j;
    }

    unsigned long prev()
    // Return position of rightmost change, return k with last composition.
    {
        const unsigned long v = x_[0];
        // value at first position

        if ( n_==v ) return k_;   // current composition is first
        x_[0] = 0;                  // set first position to zero
        unsigned long j = 1;              // set first position to zero

        while ( 0==x_[j] ) ++j;   // find next nonzero

        --x_[j];                  // decrement value
        x_[j-1] = 1 + v;          // set previous position
        return     j;

    }

    composition_colex(unsigned long n, unsigned long k)
    {
    // Must have n>=k
        n_ = n;
        k_ = k;
        nk1_ = n - k + 1;                 // must be >= 1
        if ( (long)nk1_ < 1 ) nk1_ = 1; // avoid hang with invalid pair n,k
        x_ = new unsigned long[nk1_];
        b_ = new unsigned long[k_];
        x_[nk1_-1] = 0;                         // not one
        first();
    }

    void comp2comb(const unsigned long *p, unsigned long nk1_, unsigned long *b)
    // Convert composition P(*, k) in p[] to combination in b[]
    {
        for (unsigned long j=0,i=0,z=0; j< nk1_; ++j)
        {
            unsigned long pj = p[j];
            for (unsigned long w=0; w<pj; ++w)
                b[i++] = z++;
            ++z;
        }
    }

    void print(unsigned long *b_)
    {
        for (unsigned long i=0; i<k_;i++)
            cout << b_[i];
        cout << "\n";
    }

};

class composition_colex2
{
    public:
    unsigned long n_, k_; // composition of n into k parts
    unsigned long *x_;
    unsigned long p0_;

    unsigned long next()
    // Return position of rightmost change, return k with last composition.
    {
        unsigned long j = p0_; // position of first nonzero
        if ( j==k_-1 )
            return k_;
        unsigned long v = x_[j];
        x_[j] = 0;
        --v;
        x_[0] = v;
        ++p0_;
        if ( 0!=v )  p0_ = 0;
            ++j;
        ++x_[j];
        return j;
        // current composition is last
        // value of first nonzero
        // set to zero
        // value-1 to first position
        // first nonzero one more right except ...

        // ... if value v was not one
        // increment next position

    }

};

class composition_ex_colex
{
public:
    unsigned long n_, k_; // composition of n into exactly k parts
    unsigned long *x_;
    // data (k elements)
    unsigned long nk1_;
    // ==n-k+1

public:
    composition_ex_colex(unsigned long n, unsigned long k)
    // Must have n>=k
    {
        n_ = n;
        k_ = k;
        nk1_ = n - k + 1;                 // must be >= 1
        if ( (long)nk1_ < 1 ) nk1_ = 1; // avoid hang with invalid pair n,k
        x_ = new unsigned long[k_ + 1];
        x_[k] = 0;                         // not one
        first();
    }

    void first()
    {
        x_[0] = nk1_; // all in first position
        for (unsigned long k=1; k<k_; ++k) x_[k] = 1;
    }

    void last()
    {
        for (unsigned long k=0; k<k_; ++k) x_[k] = 1;

        x_[k_-1] = nk1_; // all in last position
    }


    unsigned long next()
    // Return position of rightmost change, return k with last composition.
    {
        unsigned long j = 0;
        while ( 1==x_[j] ) ++j;    // find first greater than one

        if ( j==k_ )               // current composition is last
            return k_;
        unsigned long v = x_[j];   // value of first greater one
        x_[j] = 1;                   // set to 1
        x_[0] = v - 1;               // value-1 to first position
        ++j;
        ++x_[j];                   // increment next position

        return j;
    }

    unsigned long  prev()
    // Return position of rightmost change, return k with last composition.
    {
        const unsigned long  v = x_[0];     // value at first position
        if ( nk1_==v )
            return k_;             // current composition is first
        x_[0] = 1;                // set first position to 1
        unsigned long  j = 1;
        while ( 1==x_[j] )    ++j;// find next greater than 1

        --x_[j];
        x_[j-1] = 1 + v;   // set previous position
        return     j;
    }


    inline void comp2comb(const unsigned long *p, unsigned long k, unsigned long *b)
    // Convert composition P(*, k) in p[] to combination in b[]
    {
        for (unsigned long j=0,i=0,z=0; j<k; ++j)
        {
            unsigned long pj = p[j];
            for (unsigned long w=0; w<pj; ++w)
                b[i++] = z++;
            ++z;
        }
    }

};