/*
 * gencombo.c  --  Functions to generate combinations.
 * Copyright (c) 2006 by James Dow Allen
 *
 *  This program can be distributed or modified subject to
 *  the terms of the GNU General Public License (version 2
 *  or, at your option, any later version) as published by
 *  the Free Software Foundation.
 *
 *
 * Two different beautiful algorithms from Knuth are provided.
 * (See Fund. Comp. Prog. 7.2.1.3 preprint by Donald Knuth.)
 *	*  A revolving-door algorithm
 *	*  Random-access to combinations
 * Coding by James Dow Allen (jamesdowallen at gmail)
 *
 * In either case, the routine creates m-sized combinations
 *	of {0, 1, 2, ..., n-1}.  (Often you will use these
 *	as indices into an array whose elements are those
 *	the application wishes to combinate.)
 * Note: `m' must be in range:   0 < m <= n
 *
 * Random-access assigns a canonical number to each combination
 *   and allows non-successive retrieval; it uses Pascal's
 *   Triangle.
 *
 * Revolving-door is used to retrieve combinations in sequence,
 *   with successive results differing in only two positions.
 *   The order in which Revolving-door presents combinations
 *   is NOT the canonical order used by the Random-access method.
 *
 * To demonstrate the routines, compile via
 * 	gcc -Wall -DTESTING -O -o combo_test gencombo.c
 * and get usage instructions via
 *	combo_test
 * The 'B' option uses 'R' (revolving-door) algorithm, but
 * with alternate formatting support.
 *
 * Modify this module to support your needs (e.g. remove Pascal
 *   Triangle generation if you don't need it), but do not
 *   remove license notice or attribution at top.
 */

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

/*
 * Number types we use working with C(n,m) :
 *   Dindex admits {0, 1, 2, ..., n}
 *   Combnum admits {0, 1, 2, ..., C(n,m)}, but also some
 *	negative numbers in processing internal to kth_combo().
 */
typedef unsigned short Dindex;
typedef long Combnum;



/* ========== Code assoc. with Random-access Combination =========  */

/*
 * There are two ways to store Pascal's Triangle:
 *   (1) as a Triangle
 *   (2) as a rectangle truncated by MaxM
 * Either works, but the latter saves space when N is much larger than M.
 * Turn on FULL_PASCT to get the former.
 *
 * The variable MaxM shouldn't really exist (should just be a
 *	macro arg) but helps keep these *Variable* macros simple.
 */

#if	1
#define	FULL_PASCT
#endif

/* The size of a Triangle: */
#ifdef	FULL_PASCT
#define	Size_Triangle(n)	((((n)+2)*((n)+1))/2)
#else
static Dindex	MaxM;
#define	Size_Triangle(n)	(((n)+1) * (MaxM + 1))
#endif

/* Triangular addressing: */
#define	PasTri(ptri, n, m)	((ptri)[(m) + Size_Triangle(n-1)])

/*
 *  Initialize pasctri[] to the first maxn rows of Pascal's Triangle
 */
void mk_pasc_tri(Dindex maxn, Combnum *pasctri)
{
	Dindex	n, m;
	Combnum	*tp = pasctri;

	for (n = 0; n <= maxn; n++) {
		*tp++ = 1;
#ifdef	FULL_PASCT
		for (m = 1; m < n; m++)
#else
		for (m = 1; m < n && m <= MaxM; m++)
#endif
			*tp++ = PasTri(pasctri, n-1, m)
				+ PasTri(pasctri, n-1, m-1);
#ifdef	FULL_PASCT
		if (m == n)
			*tp++ = 1;
#else
		for (; m <= MaxM; m++)
			*tp++ = m == n;
#endif
	}
}

/* Next is just for "paranoid" checking */
int check_table(Dindex n, Dindex m, Combnum *pasctri)
{
	Dindex	k;
	Combnum	x = 1, y;

	for (k = m + 1; k <= n; k++) {
		y = PasTri(pasctri, k, m);
		if (y <= x)
			return 0;
		x = y;
	}
	return 1;
}

/*
 * Prepare the k'th combo of numbers {0,1,2,...,n-1} taken m at a time,
 *  and store it in result[].  Assumes pasctri is Pascal's triangle.
 * Return 1 on success, 0 if k is out of bounds.
 * (Note: argument `n' has no purpose except bounds check.)
 */
int kth_combo(Combnum k, Dindex n, Dindex m, Combnum *pasctri, Dindex *result)
{
	int	i, j;

	if (k < 0)
		return 0;
	for (i = m; i; i--) {
		for (j = i; PasTri(pasctri, j, i) <= k; j++) {
			if (j == n)
				return 0;
		}
		k -= PasTri(pasctri, j - 1, i);
		*result++ = j - 1;
	}
	return 1;
}

/* ========== End of Random-access Combination Code =========  */




/* ========== Code assoc. with Revolving-Door Combination =========  */

/*
 * Revolving-Door (Payne's algorithm,
 *	cited as Algorithm R in Knuth 7.2.1.3)
 * Three entry_points:
 *  revdoor_init(n, m) -- malloc's array for combo and sets the
 *	first combo.  For the algorithm's own convenience,
 *	arr[-1] is set to m, and arr[m] set to n.
 *	Returns pointer to the array, or NULL on failure.
 *  revdoor_next(arr, in, out) -- modifies the array passed to it
 *	and returns 1 on success, or 0 if all combinations have
 *	been presented. It also returns, by writing argument pointers
 *	`in' and `out', the values changed (this is suppressed
 *	if `in' is NULL).
 *  revdoor_free() -- free the malloc'ed array
 */

Dindex *revdoor_init(Dindex n, Dindex m)
{
	Dindex	*arr, j;

	if (m < 1 || m > n)
		return 0;
	arr = malloc((m + 2) * sizeof(Dindex));
	if (arr == 0)
		return 0;
	*arr++ = m;
	for (j = 0; j < m; j++)
		*arr++ = j;
	*arr = n;
	return arr - m;
}

void revdoor_free(Dindex *arr)
{
	free((void *)(arr - 1));
}

int revdoor_next(Dindex *arr, Dindex *in, Dindex *out)
{
	Dindex	j, k, a0, a1;

	k = arr[-1];
	a0 = arr[0];
	a1 = arr[1];
	if (k & 1) {
		if (a0 + 1 < a1) {
			if (in)
				*in = a0 + 1, *out = a0;
			arr[0] = a0 + 1;
			return 1;
		}
		j = 1;
	} else if (a0 > 0) {
		if (in)
			*in = a0 - 1, *out = a0;
		arr[0] = a0 - 1;
		return 1;
	} else if (a1 + 1 < arr[2]) {
		if (in)
			*in = a1 + 1, *out = a0;
		arr[0] = a1;
		arr[1] = a1 + 1;
		return 1;
	} else {
		j = 2;
		a0 = a1;
		a1 = arr[2];
	}
	while (j < k) {
		if (a1 > j) {
			if (in)
				*in = j - 1, *out = a1;
			arr[j] = a0;
			arr[j - 1] = j - 1;
			return 1;
		}
		a0 = a1;
		a1 = arr[++j];
		if (a1 + 1 < arr[j + 1]) {
			if (in)
				*in = a1 + 1, *out = a0;
			arr[j - 1] = a1;
			arr[j] = a1 + 1;
			return 1;
		}
		a0 = a1;
		a1 = arr[++j];
	}
	return 0;
}

/* ========== End of Revolving-Door Combination Code =========  */



/* ========== Simple support for Bit maps =========  */
#define	Bit_set(b, x)	(   (b)[(x) >> 3] |=   128u >> ((x) & 7))
#define	Bit_clr(b, x)	(   (b)[(x) >> 3] &= ~(128u >> ((x) & 7)))
#define	Bit_tst(b, x)	(!!((b)[(x) >> 3] &    128u >> ((x) & 7)))

void bitm_dump(unsigned char *from, int nw)
{
	nw = (nw + 7) / 8;
	while (nw--)
		printf("%02x", *from++);
	printf("\n");
}

void bitm_clear(unsigned char *from, int nw)
{
	nw = (nw + 7) / 8;
	while (nw--)
		*from++ = 0;
}

int bit_nextcombo(Dindex  *arr, unsigned char *bits)
{
	Dindex	in, out;

	if (revdoor_next(arr, &in, &out)) {
		Bit_set(bits, in);
		Bit_clr(bits, out);
		return 1;
	}
	return 0;
}

void bit_freecombo(Dindex  *arr)
{
	revdoor_free(arr);
}

Dindex  *bit_initcombo(unsigned char *bits, Dindex n, Dindex m)
{
	int	i;

	bitm_clear(bits, n);
	for (i = 0; i < m; i++)
		Bit_set(bits, i);
	return revdoor_init(n, m);
}
/* ========== End of Simple Bit map support  =========  */




/* ========== Test program for afore-going routines =========  */
#ifdef	TESTING
/* gcc -Wall -DTESTING -O -o combo_test gencombo.c */

int main(int argc, char *argv[])
{
	Dindex	*arr, in, out;
	Combnum	*pascal, ixcomb;
	unsigned char *Bits;
	Dindex	N, M, i;

	if (argc != 4) {
	usage:
		fprintf(stderr, "Generate the combinations of #N-choose-#M\n");
		fprintf(stderr, "  where #N >= #M >= 1.\n");
		fprintf(stderr,
	"Usage: combo_test #N #M #K -- for K'th combo\n");
		fprintf(stderr,
	"Usage: combo_test #N #M B -- for all combos, as bits\n");
		fprintf(stderr,
	"Usage: combo_test #N #M R -- for all combos, as revolvers\n");
		exit(1);
	}
	N = atoi(argv[1]);
	M = atoi(argv[2]);
	if (1 > M || M > N)
		goto usage;

	if (argv[3][0] == 'B') {
		printf("Demonstrating bitfields.\n");
		Bits = malloc((N + 7) / 8);
		if (Bits == 0)
			goto failure;
		arr = bit_initcombo(Bits, N, M);
		if (arr == 0)
			goto failure;
		printf("Combinations of %d taken %d at a time:\n", N, M);
		do {
			bitm_dump(Bits, N);
		} while (bit_nextcombo(arr, Bits));
		bit_freecombo(arr);
	} else if (argv[3][0] == 'R') {
		printf("Demonstrating revolving door.\n");
		arr = revdoor_init(N, M);
		if (arr == 0)
			goto failure;
		printf("Combinations of %d taken %d at a time:\n", N, M);
		in = out = 0;
		do {
			printf("In %d Out %d : ", in, out);
			for (i = 0; i < M; i++)
				printf(" %d", arr[i]);
			printf("\n");
		} while (revdoor_next(arr, &in, &out));
		revdoor_free(arr);
	} else {
		printf("Demonstrating random access.\n");
		ixcomb = atol(argv[3]);
#ifndef	FULL_PASCT
		MaxM = M;
#endif
		pascal = malloc(sizeof(Combnum) * Size_Triangle(N));
		if (pascal == 0)
			goto failure;
		mk_pasc_tri(N, pascal);
		if (! check_table(N, M, pascal)) {
			printf("These N, M too large for implementation\n"); 
		}
		arr = malloc(sizeof(Dindex) * N);
		if (arr == 0)
			goto failure;
		if (! kth_combo(ixcomb, N, M, pascal, arr)) {
			printf("Index %ld is out of range 0 <= x < %ld\n",
				(long)ixcomb, (long)PasTri(pascal, N, M));
			exit(1);
		}
		printf("The %ld'th combination of %d taken %d at a time:\n",
				ixcomb, N, M);
		for (i = 0; i < M; i++)
			printf(" %d", arr[i]);
		printf("\n");
	}
	exit(0);
    failure:
	printf("Failed: probably excessive N\n");
	exit(1);
}
#endif