// ////////////////////////////////////////////////////////////////////
//
// Search for an S-Box
//
// clang -O3 sbox-s1.c -o sbox-s1

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <fcntl.h>
#include <unistd.h>

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

#define SEPARATOR "------------------------------------------------------------------------------------------------"

void test_sbox_qmc (int * sbox, int dimension, int verbose, int * num_products, int * max_degrees, int * sums_degrees, int * used_variables);

#define OUTPUT_BFUNC_INPUT

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

#if 0
int popCount(unsigned x) { // Taken from book "Hackers Delight"
	x = x - ((x >> 1) & 0x55555555);
	x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
	x = (x + (x >> 4)) & 0x0F0F0F0F;
	x = x + (x >> 8);
	x = x + (x >> 16);
	return x & 0x0000003F;
}
#else
#define popCount __builtin_popcount
#endif

// Calculate hamming weight/distance of two integer numbers

int hammingDistance(int v1, int v2) {
	return popCount(v1 ^ v2);
} 

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

#define QMC_NO_STANDALONE
#include "./qmc.c"

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

#define limb_t uint32_t
#define BITS_IN_LIMB 32

#define DEGREE 4

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

int minimum(int a, int b)
{
	return (a<b) ? a : b;
}

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

// #define limb_t uint32_t
// #define BITS_IN_LIMB 32

#define MAXDEGREE 4
#define MAXSIZE (1<<MAXDEGREE)

int sbox[MAXSIZE];
int sinv[MAXSIZE];
int stmp[MAXSIZE];

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

// cryptographic desirable properties

int MAX_LINEARITY = 0;
int MAX_BITWISE_DIFF = 0;
int MAX_ARITH_DIFF = 0;

int MAX_FIXED_POINTS = 4;

// related to circuit complexity and depth

int MAX_WEIGHT_IN_SOP = 0;            /* how many products in a single sum  */
int MAX_SUM_WEIGHTS_IN_SOP = 0;       /* how many products in all four sums */
int MAX_DEGREE_IN_SOP = 0;            /* max degree of these products       */
int MAX_SUM_OF_DEGREES_IN_SOP = 0;    /* the sum of the degrees in a sum of products
                                       * the sum of the weights in an expression */
int MIN_USED_VARIABLES = 0;			  /* mask for used variables in all expressions */

#define COMPUTE_ANF

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

#include "walsh-anf.c"

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

void output_binary(int value, int dimension)
{
	int i, mask;
	for (i=dimension-1;i>=0;--i)
		putchar(value & (1<<i) ? '1' : '0');
}

void output_bfunc_input(int * sbox, int dimension)
{
	int size = 1 << dimension;
	int i;

	printf("\n");
	for (i=0;i<dimension;i++)
		//printf("%c", 'z' - (dimension - 1) + i);
		printf("%c", 'z' - i);

	printf(" ");
	for (i=dimension-1;i>=0;--i)
		printf("%d", i);
	printf("\n");

	for (i=0;i<size;i++)
	{
		output_binary(i,dimension);
		printf(" ");
		output_binary(sbox[i],dimension);
		printf("\n");
	}
}

int is_same(int * a, int * b, int size)
{
	for (int i=0;i<size;i++)
	{
		if (a[i] != b[i]) return 0;
	}
	
	return 1;
}

void generate_inversesbox(int * sinv, int * sbox, int size)
{
	int i;
	
	for (i=0;i<size;i++)
	{
		sinv[sbox[i]] = i;
	}
}

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

#define resetdiffs()											\
do {															\
    for (i=0;i<size;i++)										\
      for (j=0;j<size;j++)										\
         diff[i][j] = 0;										\
} while(0)

#define outputmaxdiffs()										\
do {															\
	for (i=0;i<size;i++)										\
	{															\
		int maxout = 0;											\
		for (j=0;j<size;j++)									\
			if (diff[i][j] > maxout) maxout = diff[i][j];		\
		printf("%3d ",maxout);									\
		if ((i % 16 == 15) || (i == size - 1)) printf("\n");	\
	}															\
}  while(0)

void evaluatesbox(int * sbox, int * numfix, int * maxdout_arith, int * maxdout_bitw, int dimension, int output)
{
	int i,j,k;
	int t,max;
	int size = 1 << dimension;

	int diff[size][size];

	* numfix = 0;
    for (i=0;i<size;i++)
		if (i == sbox[i])
		{
			if (output)
			{
				printf("fixed point: %d\n",i);
			}

			(*numfix)++;
		}

	resetdiffs();
	max = 0;
    for (i=0;i<size;i++)
		for (j=0;j<size;j++)
		    if (i != j)
		    {
		    	int idiff = (size+i-j) % size;
		    	int odiff = (size+sbox[i]-sbox[j]) % size;
		    	t = ++ diff[idiff][odiff];
		    	if (t>max) max=t;
			}
	*maxdout_arith = max;

	if (output)
	{
		printf("\nDOUTS (arithmetic)\n");
		outputmaxdiffs();
	}

	resetdiffs();
	max = 0;
    for (i=0;i<size;i++)
		for (j=0;j<size;j++)
		    if (i != j)
		    {
		    	int idiff = i ^ j;
		    	int odiff = sbox[i] ^ sbox[j];
		    	t = ++ diff[idiff][odiff];
		    	if (t>max) max=t;
			}
	*maxdout_bitw = max;

	if (output)
	{
		printf("\nDOUTS (bitwise)\n");
		outputmaxdiffs();
	}
}

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

void process_one_sbox(int * sbox, int * sinv, int dimension)
{
	int size = 1 << dimension;
    int i,j,lin, lin_inv;
	int mask;
	int numfix, maxdout_arith , maxdout_bitw;
	int numfix_inv, maxdout_arith_inv, maxdout_bitw_inv;

	// Store the NOT of an S-Box (also used for the NOT of the inverse)

	int sbox_not[size];

	// Stats for the S-Box

	int countz[dimension][dimension+1]; // count the bit flips

	int weights[dimension];		// weights of the SOPs, for each output bit
	int weights_not[dimension];	// weights of the NOT-SOP ---
	int weights_min[dimension];	// min of the two above

	int degrees[dimension];		// max degree of the products in the SOP, for each output bit
	int degrees_not[dimension];	// max degree of the products in the NOT-SOP, for each output bit
	int degrees_min[dimension]; // min of the two above

	int sums_degrees[dimension];	 // sum of the degrees of the products in the SOP (each output)
	int sums_degrees_not[dimension]; // same for the NOT-SOP
	int sums_degrees_min[dimension]; // min of the two abovey

	int used_variables[dimension];
	int used_variables_not[dimension];
	int used_variables_tot = 0;

	// Stats for the inverse S-Box (add "_inv", same descriptions)

	int countz_inv[dimension][dimension+1];

	int weights_inv[dimension];
	int weights_inv_not[dimension];
	int weights_inv_min[dimension];

	int degrees_inv[dimension];
	int degrees_inv_not[dimension];
	int degrees_inv_min[dimension];

	int sums_degrees_inv[dimension];
	int sums_degrees_inv_not[dimension];
	int sums_degrees_inv_min[dimension];

	int used_variables_inv[dimension];
	int used_variables_inv_not[dimension];
	/* int used_variables_tot_inv = 0; */

#ifdef COMPUTE_ANF
	int ANF[size];
	int ANF_inv[size];
#endif
	
	int is_involution;

	is_involution = is_same(sbox, sinv, size);

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

	evaluatesbox(sbox, &numfix, &maxdout_arith, &maxdout_bitw, dimension, 0);
	lin = linearity(sbox, dimension);
	
	for (i=0,mask=1; i<dimension; i++,mask<<=1)
	{
		for (j=0; j<=dimension; j++)
			countz[i][j] = 0; // i = flipped bit in input, j = count of flipped bits in output
		for (j=0;j<size;j++)
		{
			//if (!(j & mask)) /* each diff is counted twice */
			{
				int count = __builtin_popcount(sbox[j]^sbox[j^mask]);
				countz[i][count] ++;
			}
		}
	}

	if (!is_involution)
	{
		evaluatesbox(sinv, &numfix_inv, &maxdout_arith_inv, &maxdout_bitw_inv, dimension, 0);
		lin_inv = linearity(sinv, dimension);

		for (i=0,mask=1; i<dimension; i++,mask<<=1)
		{
			for (j=0; j<=dimension; j++)
				countz_inv[i][j] = 0;
			for (j=0;j<size;j++)
			{
				{
					int count = __builtin_popcount(sinv[j]^sinv[j^mask]);
					countz_inv[i][count] ++;
				}
			}	
		}
	}

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

	if ( ( (lin <= MAX_LINEARITY)
	        && (maxdout_bitw <= MAX_BITWISE_DIFF)
	        && (maxdout_arith <= MAX_ARITH_DIFF) )
	    && ( is_involution
	         || ( (lin <= MAX_LINEARITY)
	               && (maxdout_bitw_inv <= MAX_BITWISE_DIFF)
	               && (maxdout_arith_inv <= MAX_ARITH_DIFF) ) )
	        )
	{
		int max_weight = 0, max_weight_inv = 0;
		int max_degree = 0, max_degree_inv = 0;
		int max_sums_degrees = 0, max_sums_degrees_inv = 0;

		int max_weight_not = 0, max_weight_inv_not = 0;
		int max_degree_not = 0, max_degree_inv_not = 0;
		int max_sums_degrees_not = 0, max_sums_degrees_inv_not = 0;

		int max_weight_min = 0, max_weight_inv_min = 0;
		int max_degree_min = 0, max_degree_inv_min = 0;
		int max_sums_degrees_min = 0, max_sums_degrees_inv_min = 0;

		int tot_weight = 0, tot_weight_inv = 0;

		for (i=0; i<size; i++) sbox_not[i] = sbox[i] ^ (size-1);
		test_sbox_qmc(sbox, dimension, FALSE, weights, degrees, sums_degrees, used_variables);
		test_sbox_qmc(sbox_not, dimension, FALSE, weights_not, degrees_not, sums_degrees_not, used_variables_not);

		used_variables_tot = -1;
		for(i=0; i<dimension; i++)
		{
			weights_min[i] = minimum(weights[i], weights_not[i]);
			if (weights_min[i] > max_weight) max_weight = weights_min[i];

			degrees_min[i] = minimum(degrees[i], degrees_not[i]);
			if (degrees_min[i] > max_degree) max_degree = degrees_min[i];

			sums_degrees_min[i] = minimum(sums_degrees[i], sums_degrees_not[i]);
			if (sums_degrees_min[i] > max_sums_degrees) max_sums_degrees = sums_degrees_min[i];

			tot_weight += weights_min[i];

			used_variables_tot &= used_variables[i]; /* this same for function and not */
		}

		if (!is_involution)
		{
			for (i=0; i<size; i++) sbox_not[i] = sinv[i] ^ (size-1);
			test_sbox_qmc(sinv, dimension, FALSE, weights_inv, degrees_inv, sums_degrees_inv, used_variables_inv);
			test_sbox_qmc(sbox_not, dimension, FALSE, weights_inv_not, degrees_inv_not, sums_degrees_inv_not, used_variables_inv_not);

			/* used_variables_tot_inv = -1; */
			for(i=0; i<dimension; i++)
			{
				weights_inv_min[i] = minimum(weights_inv[i], weights_inv_not[i]);
				if (weights_inv_min[i] > max_weight_inv) max_weight_inv = weights_inv_min[i];

				degrees_inv_min[i] = minimum(degrees_inv[i], degrees_inv_not[i]);
				if (degrees_inv_min[i] > max_degree_inv) max_degree_inv = degrees_inv_min[i];

				sums_degrees_inv_min[i] = minimum(sums_degrees_inv[i], sums_degrees_inv_not[i]);
				if (sums_degrees_inv_min[i] > max_sums_degrees_inv) max_sums_degrees_inv = sums_degrees_inv_min[i];

				tot_weight_inv += weights_inv_min[i];

				/* used_variables_tot_inv &= used_variables_inv[i]; */
			}
		}

		if ( ( (max_weight <= MAX_WEIGHT_IN_SOP) && (tot_weight <= MAX_SUM_WEIGHTS_IN_SOP)
		         && (max_degree <= MAX_DEGREE_IN_SOP)
		         && (max_sums_degrees <= MAX_SUM_OF_DEGREES_IN_SOP)
		         && (used_variables_tot >= MIN_USED_VARIABLES)
		         && (numfix <= MAX_FIXED_POINTS) )
		     && ( is_involution
		          || ( (max_weight_inv <= MAX_WEIGHT_IN_SOP) && (tot_weight_inv <= MAX_SUM_WEIGHTS_IN_SOP)
		                 && (max_degree_inv <= MAX_DEGREE_IN_SOP)
		                 && (max_sums_degrees_inv <= MAX_SUM_OF_DEGREES_IN_SOP)
				         /* && (used_variables_tot_inv >= MIN_USED_VARIABLES) */
				         ) )
		    )
		{
			printf("SBOX = { ");
			for (i=0;i<size;i++)
			{
				printf("%d",sbox[i]);
				if (i < size - 1) printf(", "); else printf(" }");
				if ((i % 16 == 15) || (i == size - 1)) printf("\n");
			}

#if 0
			printf("SINV = { ");
			for (i=0;i<size;i++)
			{
				printf("%d",sinv[i]);
				if (i < size - 1) printf(", "); else printf(" }");
				if ((i % 16 == 15) || (i == size - 1)) printf("\n");
			}
#endif

			printf("\nfixed points %d, ",numfix);
			printf("max arith diff %d, ",maxdout_arith);
			printf("max bitwise diff %d, ",maxdout_bitw);
			printf("linearity %d\n",lin);

			for (i=0;i<dimension;i++)
			{
				printf("flipping bit %d in input, counts for 0..DEG flipped bits = ",i);
				for (j=0;j<=dimension;j++) printf("%d ", countz[i][j]/2);
				printf("\n");
			}

			for (i=0; i<size; i++) sbox_not[i] = sbox[i] ^ (size-1);
			printf("\nQMC of SBOX\n");
			test_sbox_qmc(sbox, dimension, 1, weights, degrees, sums_degrees, used_variables);

			printf("QMC of not SBOX\n");
			test_sbox_qmc(sbox_not, dimension, 1, weights_not, degrees_not, sums_degrees_not, used_variables_not);

			printf("\nweights: ");
			for (i=0;i<dimension;i++) printf("(%d,%d) ", weights[i], weights_not[i]);

			printf("    max weight = %d", max_weight);
			printf("    tot weight = %d", tot_weight);
			printf("    max degree = %d", max_degree);
			printf("    max sum of degrees = %d", max_sums_degrees);

			printf("\nused variables: ");
			for (i=0;i<dimension;i++) for (i=0;i<dimension;i++) printf("%d ", used_variables[i]);
			printf("[%d]\n", used_variables_tot);

#ifdef COMPUTE_ANF
			printf("ANF degrees = [");
			for (i=1;i<size;i++)
			{
				int extractor = i; // (1<<i);
				int w = anf(ANF, sbox, extractor, dimension, FALSE);
				printf("%d", w);
				if (i<size-1) printf(" ");
			}
			printf("]\n");
#endif

#ifdef OUTPUT_BFUNC_INPUT
				output_bfunc_input(sbox, dimension);
#endif

			if (!is_involution)
			{
				printf("\ntesting inverse sbox\n");
			
				printf("fixed points %d, ",numfix_inv);
				printf("max arith diff %d, ",maxdout_arith_inv);
				printf("max bitwise diff %d, ",maxdout_bitw_inv);
				printf("linearity %d\n",lin_inv);

				for (i=0;i<dimension;i++)
				{
					printf("flipping bit %d in input, counts for 0..DEG flipped bits = ",i);
					for (j=0;j<=dimension;j++) printf("%d ", countz_inv[i][j]/2);
					printf("\n");
				}

				for (i=0; i<size; i++) sbox_not[i] = sinv[i] ^ (size-1);
				printf("\nQMC of SBOX_inv\n");
				test_sbox_qmc(sinv, dimension, 1, weights_inv, degrees_inv, sums_degrees_inv, used_variables_inv);
				printf("QMC of not SBOX_inv\n");
				test_sbox_qmc(sbox_not, dimension, 1, weights_inv_not, degrees_inv_not, sums_degrees_inv_not, used_variables_inv_not);

				printf("\nweights:");
				for (i=0;i<dimension;i++) printf("(%d,%d) ", weights_inv[i], weights_inv_not[i]);

				printf("   max weight = %d", max_weight_inv);
				printf("   tot weight = %d", tot_weight_inv);
				printf("   max degree = %d", max_degree_inv);
				printf("   max sum of degrees = %d", max_sums_degrees_inv);

				//printf("\nused variables: ");
				//for (i=0;i<dimension;i++) for (i=0;i<dimension;i++) printf("%d ", used_variables_inv[i]);
				//printf("[%d]\n", used_variables_tot_inv);

#ifdef COMPUTE_ANF
				printf("ANF degrees = [");
				for (i=1;i<size;i++)
				{
					int extractor = i; //(1<<i);
					int w = anf(ANF, sinv, extractor, dimension, FALSE);
					printf("%d", w);
					if (i<size-1) printf(" ");
				}
				printf("]\n");
#endif			

#ifdef OUTPUT_BFUNC_INPUT
					output_bfunc_input(sinv, dimension);
#endif
			}
			printf("\n");
			printf(SEPARATOR);
			printf("\n");
			printf("\n");
		}
	}
}

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

#define SIZE 16

#include "./search.c"

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

int main ()
{
	printf(SEPARATOR"\n");

	int S[SIZE];

	MAX_WEIGHT_IN_SOP = 4;
	MAX_SUM_WEIGHTS_IN_SOP = 16;
	MAX_DEGREE_IN_SOP = 3;
	MAX_SUM_OF_DEGREES_IN_SOP = 10;
	MIN_USED_VARIABLES = 15;

	MAX_LINEARITY = 8;
	MAX_BITWISE_DIFF = 4;
	MAX_ARITH_DIFF = 3;

	num_fix_points = 0;

	create_empty_involution(SIZE-num_fix_points, S);
	fpfi(SIZE-num_fix_points,SIZE, S);

	return 0;
}

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