// ////////////////////////////////////////////////////////////////

// Quine-McCluskey Algorithm
// Derived from the public domain implementation by 
// Stefan Moebius (mail@stefanmoebius.de) [05/16/2012]
//
// License: Can be used freely (Public Domain)

#include <stdio.h>

#define TRUE 1
#define FALSE 0

// Count all set bits of the integer number

void explicitTerm(int bitfield, int mask, int dimension) {
	if (mask) {
		int z;
		int count = 0;

		for (z = 0; z < dimension; z++) {
			if (mask & (1 << z)) {
				if (bitfield & (1 << z))
				{
					//if (count>0) printf("*");
					printf("%c", 'z' - (dimension - 1) + z);
					count = 1;
				}
				else		
				{
					//if (count>0) printf("*");
					// printf("(1-%c)", 'z' - (dimension - 1) + z);
					printf("%c'", 'z' - (dimension - 1) + z);
					count = 1;
				}
			} 
		} 
	}
}

// Determines whether "value" contains "part"

int contains(value, mask, part, partmask)
{
	if ((value & partmask) == (part & partmask))
	{
		if ((mask & partmask) ==  partmask)
			return TRUE;
	}   
	return FALSE;
}

// Returns the number of products in the sum,
// and also the max of the degree of such products
// and also the sums of the degrees 
// and also the mask of the contained variables

void main_qmc(int dimension, int * table, int * num_products, int * max_degrees, int * sum_degrees, int * used_variablesp, int verbose)
{
	int MAXPOS = (1<<(dimension));	// 2 ^ dimension

	int minterm[MAXPOS][MAXPOS];
	int mask[MAXPOS][MAXPOS];	// mask of minterm
	int used[MAXPOS][MAXPOS];	// minterm used
	int primemask[MAXPOS];		// mask for prime implicants
	int prime[MAXPOS];			// prime implicant
	int isepi[MAXPOS];			// <IS> it an <E>ssential <P>rime <I>mplicant? (TRUE/FALSE)
	int rnepi[MAXPOS];			// <R>equired <N>on <E>ssential <P>rime <I>mplicant  (TRUE/FALSE)

	int cur = 0;
	int reduction = 0;			// Reduction step
	int reduction_done = FALSE;
	int prim_count = 0;
	int term = 0;
	int termmask = 0;
	int found = 0;
	int used_variables = 0;
	int x = 0;
	int y = 0;
	int z = 0;
	int count = 0;
	int max_degree = 0;
	int lastprime = 0; 
	int value = 0;				// Result of evaluation of an expression
	int sum_degs = 0;

	// ////////////////////////////////////////////////////////////////

	// Fill all arrays with default values

	for (x = 0; x < MAXPOS; x++)
	{
		primemask[x] = 0;
		prime[x] = 0;
		isepi[x] = FALSE;		// no monomial is a prime implicant at beginning
		rnepi[x] = TRUE;		// mark non-essential prime implicants as necessary
		for (y = 0; y < MAXPOS; y++)
		{
			mask[x][y] = 0;
			minterm[x][y] = 0;
			used[x][y] = FALSE;
		}
	}

	cur = 0;
	for ( x=0; x < MAXPOS; x++)
	{
		if (table[x])
		{
			mask[cur][0] = ((1 << dimension)- 1);   
			minterm[cur][0] = x;  
			cur++;  
		}
	}	
	
	// ////////////////////////////////////////////////////////////////

	for (reduction = 0; reduction < MAXPOS; reduction++)
	{
		cur = 0; 
		reduction_done = FALSE;
		for (y=0; y < MAXPOS; y++)
		{
			for (x=0; x < MAXPOS; x++)
			{   
				if ((mask[x][reduction]) && (mask[y][reduction]))
				{      
					if (popCount(mask[x][reduction]) > 1)
					// Do not allow complete removal (problem if all terms are 1)
					{
						if ((hammingDistance(minterm[x][reduction] & mask[x][reduction],
						                     minterm[y][reduction] & mask[y][reduction]) == 1)
						      && (mask[x][reduction] == mask[y][reduction]))
						{
						    // Simplification only possible if 1 bit differs
							term = minterm[x][reduction]; // could be mintern x or y
							// e.g.:
							// 1110
							// 1111
							// Should result in mask of 1110
							termmask = mask[x][reduction]  ^ (minterm[x][reduction] ^ minterm[y][reduction]); 
							term  &= termmask;

							found = FALSE;		

							for ( z=0; z<cur; z++)
							{
								if ((minterm[z][reduction+1] == term) && (mask[z][reduction+1] == termmask) )
									found = TRUE;
							}

							if (found == FALSE)
							{
								minterm[cur][reduction+1] = term;
								mask[cur][reduction+1] = termmask;
								cur++; 
								if (cur>=MAXPOS) { printf("\nWhat is happening? cur = %d\n",cur); fflush(stdout); goto breakout;}
							}
							used[x][reduction] = TRUE;
							used[y][reduction] = TRUE;  
							reduction_done = TRUE;
						}
					}
				} 
			}    
		}
breakout:
		if (reduction_done == FALSE)
			break; //exit loop early (speed optimisation)
	}

	prim_count = 0;
	for ( reduction = 0 ; reduction < MAXPOS; reduction++)
	{
		for ( x=0 ;x < MAXPOS; x++)
		{
			//Determine all not used minterms
			if ((used[x][reduction] == FALSE) && (mask[x][reduction]) )
			{
				//Check if the same prime implicant is already in the list
				found = FALSE;
				for ( z=0; z < prim_count; z++)
				{
					if (((prime[z] & primemask[z]) == (minterm[x][reduction] & mask[x][reduction])) &&  (primemask[z] == mask[x][reduction]) )					
						found = TRUE;
				} 
				if (found == FALSE)
				{
					//outputTerm(minterm[x][reduction], mask[x][reduction], dimension);
					//printf("\n");
					primemask[prim_count] = mask[x][reduction];
					prime[prim_count] = minterm[x][reduction];
					prim_count++;
				}     
			} 
		} 
	}  

	// Find essential and not essential prime implicants
	// All prime implicants are set to "not essential" so far

	for (y=0; y < MAXPOS; y++)  // for all minterms
	{
		count = 0;
		lastprime = 0;   
		if (mask[y][0])
		{
			for (x=0; x < prim_count; x++ )  // for all prime implicants
			{
				if (primemask[x]) {
					// Check if the minterm contains prime implicant
					if (contains(minterm[y][0], mask[y][0], prime[x], primemask[x])) {					
						count++;
						lastprime = x;          
					}  
				} 		
			}
			// If count = 1 then it is a essential prime implicant
			if (count == 1)
			{
				isepi[lastprime] = TRUE;
			}
		}
	}

	// successively testing if it is possible to remove prime implicants from the rest matrix

	for ( z=0; z < prim_count; z++) {
		if (primemask[z] ) {
			if (isepi[z] == FALSE) { // && (rnepi[z] == TRUE) ?)
				rnepi[z] = FALSE;    // mark as "not essential"
				for ( y=0; y < MAXPOS; y++) { // for all possibilities
					value = 0;
					for ( x=0; x < prim_count; x++) {
						if ( (isepi[x] == TRUE) || (rnepi[x] == TRUE)) {  // Essential prime implicant or marked as required
							if ((y & primemask[x]) == (prime[x] & primemask[x])) { // All bits must be 1
								value = 1; 
								break;
							}
						}
					}
					// printf(" %d\t%d\n", result, result[y]);
					if (value == table[y]) {  // compare calculated result with input value
						// printf("not needed\n"); // prime implicant not required
					}
					else {
						//printf("needed\n");
						rnepi[z] = TRUE; //prime implicant required  /  Primimplikant wird doch benoetigt
					}
				}
			}
		}
	}

	// ////////////////////////////////////////////////////////////////
	
	if (verbose)
	{
		for ( x = 0 ; x < MAXPOS; x++)
		{
			printf("%d", table[x]);
		}
	}
	
	if (verbose) printf(" -- Result: ");

	// Output of essential and required prime implicants
	count = 0;
	sum_degs = 0;
	used_variables = 0;
	for ( x = 0 ; x < prim_count; x++)
	{
		if (isepi[x] == TRUE)
		{
			if (count > 0) { if (verbose) printf(" + "); }
			if (verbose) explicitTerm(prime[x], primemask[x], dimension);
			used_variables |= primemask[x];
			count ++;
			y = popCount(primemask[x]);
			sum_degs += y;
			if (y > max_degree) max_degree = y;
		}
		else if ((isepi[x] == FALSE) && (rnepi[x] == TRUE))
		{
			if (count > 0) { if (verbose) printf(" + "); }
			if (verbose)
			{
				//printf(" _ ");
				explicitTerm(prime[x], primemask[x], dimension);
			}
			used_variables |= primemask[x];
			count ++;
			y = popCount(primemask[x]);
			sum_degs += y;
			if (y > max_degree) max_degree = y;
		}
	}

	//if (verbose) printf("\t\t(sum of degrees = %d)", sum_degs);
	if (verbose) printf("\n");

	*num_products = count;
	*max_degrees = max_degree;
	*sum_degrees = sum_degs;
	*used_variablesp = used_variables;
}

// returns the number of monomials in each of the degree many sums of products

void test_sbox_qmc (int * sbox, int dimension, int verbose, int * num_products, int * max_degrees, int * sums_degrees, int * used_variables)
{
	int size = 1<<dimension;
	int table[size];		// Truth table
	int weight;

	for (int i = 0; i < dimension; i ++) /* start with least significant bit, move up */
	{
		for (int j = 0; j < size; j ++)
		{
			table[j] = (sbox[j] >> i) & 1; /* start with least significant bit, move up */
		}
		main_qmc(dimension, table, num_products+i, max_degrees+i, sums_degrees+i, used_variables+i, verbose);
	}
}

#if !defined(QMC_NO_STANDALONE)
//int sbox[16] =
// {12, 14, 3, 2, 11, 6, 5, 10, 15, 13, 7, 4, 0, 9, 1, 8};
//	{5, 13, 4, 14, 2, 0, 7, 6, 15, 11, 12, 9, 10, 1, 3, 8}; // ***
// {10, 3, 11, 1, 8, 7, 13, 5, 4, 14, 0, 2, 15, 6, 9, 12};
// {12,10,13,3,14,11,15,7,8,9,1,5,0,2,4,6}; // midori
// {0xB,0xF,0x3,0x2,0xA,0xC,0x9,0x1,0x6,0x7,0x8,0x0,0xE,0x5,0xD,0x4}; // prince
// { 11, 14, 10,  4,  3, 15, 13,  9, 12,  7,  2,  0,  8,  6,  1,  5 };
// {   3, 7,11, 0, 5, 4  , 14, 1,10,13, 8, 2,15, 9, 6,12  };
// {1, 0, 7, 11, 14, 10, 9, 2, 12, 6, 5, 3, 8, 15, 4, 13};
// {9, 14, 3, 2, 12, 8, 11, 13, 5, 0, 15, 6, 4, 7, 1, 10};
// {4, 2, 1, 10, 0, 12, 14, 11, 9, 8, 3, 7, 5, 15, 6, 13};
// {3, 7, 6, 0, 12, 9, 2, 1, 14, 5, 11, 10, 4, 15, 8, 13};

int sbox[256] =
{ 0, 187,  52,  51,  34, 238,  61, 188, 170, 153, 152,  49, 103, 102, 245, 255, 
 74, 233, 217, 251, 234, 202, 219, 201,  90,  89,  73, 249, 250, 242, 218, 210, 
 68,  45,   4,  36,  35, 205, 228, 237, 113, 116,  84, 164, 161,  33, 132, 196, 
 66,  11,  67,   3,   2, 130, 131, 139,  82,  91,  75,  83,  86,   6,  70, 134, 
 64, 186,  48,  50,  32, 192,  62, 176, 112,  26,  16,  58,  96, 110, 208, 222, 
 80, 185,  56,  59,  42, 224,  60, 184, 122,  25,  24,  57, 106,  98, 240, 254, 
 76, 172,  93, 173, 171, 140,  13,  12, 123, 124,  92, 125, 115, 163,  77, 141, 
 72,  40, 248, 108,  41, 200, 252, 232, 121, 120,  88, 104, 105, 107, 216, 220, 
128, 178,  53,  54,  46, 206,  63, 190, 160, 154, 144,  55, 101, 111, 213, 223, 
138, 235, 209, 243, 226, 194, 211, 203,  10,   9, 137, 241, 247, 246, 215, 214, 
136,  44, 244, 109,  43, 204, 253, 236, 169, 168,   8, 100,  97,  99, 212, 221, 
 71, 227, 129, 225, 230, 198, 193, 195,  87,  81,  65,   1,   7, 231, 135, 199, 
 69, 182, 149, 183,  47, 207, 181, 191, 117,  23,  21, 151, 165,  37, 133, 197, 
 78, 146,  31, 150, 174, 142, 159, 158, 126,  18,  30,  22, 127, 175,  79, 143, 
 85, 179, 148, 177,  38, 239, 180, 189, 119,  17,  20, 145, 167,  39,   5, 229, 
 94, 155,  29, 147, 162,  14, 157, 156, 114,  27,  28,  19, 118, 166,  95,  15}; 

int main () {
	int dimension = 8; 		// Number of Variables
	int result[1<<dimension];		// Results
	int weight;
	int max_degree;
	int sum_degree;
	int used_variables;

	for (int i = 0; i < 8; i ++)
	{
		for (int j = 0; j < 256; j ++)
		{
			result[j] = (sbox[j] >> i) & 1;
		}
	
		main_qmc(dimension, result, &weight, &max_degree, &sum_degree, &used_variables, 1);
		printf("weight of bit %d is %d\n", i, weight);	
	}
}
#endif