/*
	Author:	James Pate Williams (c) 2000

	A naive PC algorithm from section 4.3.2 of
	_Foundations of Constraint Satisfaction_
	by E. P. K. Tsang page 92. See Figure 3.3
	page 65.
*/

#include <iostream.h>
#include <stdlib.h>

bool satisfies(int i, int a, int j, int b) {
// check to see if the compound label <i, a><j, b> satisfies the constraint
	if (i == 0) {
		if (j == 1) {
			if (a == 1 && b == 2) return true;
			if (a == 5 && b == 6) return true;
		}
		if (j == 3) {
			if (a == 1 && b == 4) return true;
			if (a == 5 && b == 8) return true;
		}
	}
	else if (i == 1) {
		if (j == 0) {
			if (a == 1 && b == 2) return true;
			if (a == 5 && b == 6) return true;
		}
		if (j == 2) {
			if (a == 2 && b == 3) return true;
			if (a == 6 && b == 7) return true;
		}
	}
	else if (i == 2) {
		if (j == 1) {
			if (a == 2 && b == 3) return true;
			if (a == 6 && b == 7) return true;
		}
		if (j == 3) {
			if (a == 3 && b == 4) return true;
			if (a == 7 && b == 8) return true;
		}
	}
	else if (i == 3) {
		if (j == 0) {
			if (a == 1 && b == 4) return true;
			if (a == 5 && b == 8) return true;
		}
		if (j == 2) {
			if (a == 3 && b == 4) return true;
			if (a == 7 && b == 8) return true;
		}
	}
	return false;
}

bool ****generateConstraint(int n, int *nD, int **D) {
// generate the n by n constraint matrices
	bool ****C = new bool***[n];
	int i, j, k, l;

	for (i = 0; i < n; i++) {
		C[i] = new bool**[n];
		for (j = 0; j < n; j++) {
			C[i][j] = new bool*[nD[i]];
			for (k = 0; k < nD[i]; k++) {
				C[i][j][k] = new bool[nD[j]];
				for (l = 0; l < nD[j]; l++)
					C[i][j][k][l] = satisfies(i, D[i][k], j, D[j][l]);
			}
		}
	}
	
	return C;
}

int **generateDomain(int n, int *nD) {
// generate the domains of the n variables
	int **domain = new int*[n];

	for (int i = 0; i < n; i++)
		domain[i] = new int[nD[i]];
	domain[0][0] = 1;
	domain[0][1] = 5;
	domain[1][0] = 2;
	domain[1][1] = 6;
	domain[2][0] = 3;
	domain[2][1] = 7;
	domain[3][0] = 4;
	domain[3][1] = 8;
	return domain;
}

void boolMultiply(bool **A, bool **B, bool **C, int n) {
	bool sum;
	int i, j, k;

	for (i = 0; i < n; i++) {
		for (j = 0; j < n; j++) {
			sum = false;
			for (k = 0; k < n; k++)
				sum = sum || A[i][k] && B[k][j];
			C[i][j] = sum;
		}
	}
}

void printMatrix(bool **matrix, int n) {
	int i, j;

	for (i = 0; i < n; i++) {
		for (j = 0; j < n; j++)
			cout << matrix[i][j] << ' ';
		cout << endl;
	}
	cout << endl;
}

bool ****PC1(int n) {
	bool equal;
	bool **temp0 = new bool*[n], **temp1 = new bool*[n];
	bool *****Y = new bool****[n], ****Z = new bool***[n];
	int i, j, k, l, m;
	int *nD = new int[n];

	for (i = 0; i < n; i++) nD[i] = 2;
	int **D = generateDomain(n, nD);
	bool ****C = generateConstraint(n, nD, D);
	for (i = 0; i < nD[0]; i++) {
		temp0[i] = new bool[nD[0]];
		temp1[i] = new bool[nD[0]];
	}
	for (i = 0; i < n; i++) {
		Y[i] = new bool***[n];
		Z[i] = new bool**[n];
		for (j = 0; j < n; j++) {
			Y[i][j] = new bool**[n];
			Z[i][j] = new bool*[nD[0]];
			for (k = 0; k < n; k++) {
				Y[i][j][k] = new bool*[nD[0]];
				for (l = 0; l < nD[0]; l++) 
					Y[i][j][k][l] = new bool[nD[0]];
			}
			for (k = 0; k < nD[0]; k++)
				Z[i][j][k] = new bool[nD[0]];
		}
	}
	for (i = 0; i < n; i++)
		for (j = 0; j < n; j++)
			for (k = 0; k < nD[0]; k++)
				for (l = 0; l < nD[0]; l++)
					Y[n - 1][i][j][k][l] = C[i][j][k][l];
	for (i = 0; i < n - 1; i++) {
		for (j = i + 1; j < n; j++) {
			cout << "i = " << i << " j = " << j << endl << endl;
			printMatrix(C[i][j], nD[0]);
		}
	}
	do {
		for (i = 0; i < n; i++)
			for (j = 0; j < n; j++)
				for (k = 0; k < nD[0]; k++)
					for (l = 0; l < nD[0]; l++)
						Z[i][j][k][l] = Y[n - 1][i][j][k][l];
		for (k = 0; k < n; k++) {
			for (i = 0; i < n; i++) {
				for (j = 0; j < n; j++) {
					if (k == 0) {
						boolMultiply(Z[i][k], Z[k][k], temp0, nD[0]);
						boolMultiply(temp0, Z[k][j], temp1, nD[0]);
						for (l = 0; l < nD[0]; l++)
							for (m = 0; m < nD[0]; m++)
								Y[0][i][j][l][m] = Z[i][j][l][m] && temp1[l][m];
					}
					else {
						boolMultiply(Y[k - 1][i][k], Y[k - 1][k][k], temp0, nD[0]);
						boolMultiply(temp0, Y[k - 1][k][j], temp1, nD[0]);
						for (l = 0; l < nD[0]; l++)
							for (m = 0; m < nD[0]; m++)
								Y[k][i][j][l][m] = Y[k - 1][i][j][l][m] && temp1[l][m];
					}
				}
			}
		}
		equal = true;
		for (i = 0; equal && i < n; i++)
			for (j = 0; equal && j < n; j++)
				for (k = 0; equal && k < nD[0]; k++)
					for (l = 0; equal && l < nD[0]; l++)
						equal = Z[i][j][k][l] == Y[n - 1][i][j][k][l];
	} while (!equal);
	for (i = 0; i < n; i++)
		for (j = 0; j < n; j++)
			for (k = 0; k < nD[0]; k++)
				for (l = 0; l < nD[0]; l++)
					C[i][j][k][l] = Y[n - 1][i][j][k][l];
	return C;
}
					
int main(int argc, char *argv[]) {
	int i, j, n = 4;
	bool ****C = PC1(n);
	for (i = 0; i < n - 1; i++) {
		for (j = i + 1; j < n; j++) {
			cout << "i = " << i << " j = " << j << endl << endl;
			printMatrix(C[i][j], 2);
		}
	}
	return 0;
}