////////////////////////////////////////////////////////////////////////////////////////
//
//  Copyright (C) 2003 Jaap Spies, jaapspies@gmail.com
//
//  Distributed under the terms of the GNU General Public License (GPL):
//
//                  http://www.gnu.org/licenses/
//
// Quick and dirty implementation of NAW problem 29. Author: Jaap Spies
// Less dirty by adding more comments 20041010
////////////////////////////////////////////////////////////////////////////////////////

#include <stdio.h>

main()
{
	int n, h;
	printf("Enter n > 0: ");
	scanf("%d", &n);
	printf("Enter h >= 0: ");
	scanf("%d", &h);
	problem29(n, h);
}

problem29(int n, int h)
{
	int i, j;
	int *c;

	
	// this algorithm generates all possible n-subsets of the set {1,2,...n+h)
	// it is based on algorithm L from Knuth's taocp part 4: 7.2.1.3 p. 4
	// to be publiced (see his homepage)
	// Lexicographic combinations
	
	c = (int *) malloc((n+3)*sizeof(int));
	// L1. Initialize
	for(j=1; j<=n;j++)
		c[j] = j-1;
	c[n+1] = n+h;
	c[n+2] = 0;
	j = 1;

	while(j <= n)
	{
	// L2. Visit
		proces(c ,n ,h);
	// L3. Find j
		j = 1;
		while(c[j]+1 == c[j+1])
		{
			c[j] = j-1;
			j++;
		}
	// L5. Increase c[j]
		c[j] += 1;
	}
}

proces(int *c, int n, int h)
{
	int i;
	int *A, *B, *genB();

	// we adjust the input for genB(), the matrix generator

	A = (int *) malloc(n*sizeof(int));
	for(i=0; i<n; i++)
		*(A+i) = *(c+i+1)+1;

	// print the subset A
	printf("A = {");
	for(i = 0; i < n-1; i++)
		printf("%d, ", A[i]);
	printf("%d}   ", A[n-1]);
	
	// generate matrix B
	B = genB(n, h, A);
	// calculate the permanent of matrix B
	printf("per(B) = %d\n", per(n, B));
	free(A);
}

// not needed here

printmat(int n, int *B)
{
	int row, col;
	
	for (row = 0; row < n; row++)
	{
		for (col = 0; col < n ; col++)
			printf("%d ",*(B + n*row + col));
		printf("\n");
	}
}

Q(int *B, int n, int *v)
{
	int i, row, col, prod, sum;
	int *x;
	prod = 1;
	
	// we get (0,1) from per() and translate this to (-1,1)
	// and discard v[0]
	// for details see my note "The permanent of a matrix of order n"

	x = (int *) malloc(n*sizeof(int));
	for(i=0;i<n;i++)
	{
		*(x+i) = *(v+i+1)*2 - 1;
	}
	for (row=0; row<n; row++)
	{
		sum = 0;
		for (col=0; col<n; col++)
			sum += *(B+row*n+col)*(*(x+col));
		if(sum==0)
		{
			free(x);
			return 0;
		}
		prod *= sum;
	}
	for(i=0; i<n; i++)
		if (x[i] == -1)
			prod = -prod;
	free(x);
	return prod;
}

int *genB(int n, int h, int *a)
{
	int i, j;
	int *mat;
	mat = (int *) malloc(n*n*sizeof(int));
	
	for(j=0; j<n;j++)
		for(i=0;i<n;i++)
		{
	// is a[j] allowed in position i?
			if( i+1 <= *(a+j) && *(a+j) <= h + i+1)
				*(mat + i*n + j) = 1;
			else
				*(mat + i*n + j) = 0;
		}
	return mat;
}

per(int n, int *B)
{
	int i, j, sum;
	int *x;


	// this algorthm generates all n-tuples with entries 0 and 1
	// it is based on the algorithm M from Knuth's taocp part 4: 7.2.1.1
	// to be publiced (See his homepage).
	// However we need entries -1 and 1. we do the processing in Q(B, n, x)
	// See my note

	x = (int *) malloc((n+1)*sizeof(int));
	// M1. Initialize
	for (i=0; i < n+1; i++)
		*(x+i)=0;
	sum = 0;
	while(j > 0)
	{
	// M2. Visit
		sum += Q(B, n, x);
	// M3. Prepare to add 1
		j = n;
	// M4. Carry if necessary
		while( *(x+j) == 1)
		{
			*(x+j) = 0;
			j--;
		}
	// M5. Increase, unless done
		*(x+j) += 1;
	}
	// final division by 2^n see theorem 1
	for(i=0;i<n;i++)
		sum /= 2;
	return sum;
}