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

// Derived from Cyril Prissette's algorithm
// described in http://arxiv.org/pdf/1006.3993v1.pdf

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

// values 0,2 or 4
int num_fix_points = 0;
// non defined fix points must be -1.
// they must be in increasing orderi, i.e.
// fix_point_1 < fix_point_2 < fix_point_3 < fix_point_4
// where the relation applies only to defined fix_point_X values

int fix_point_1 = -1;
int fix_point_2 = -1;
int fix_point_3 = -1;
int fix_point_4 = -1;

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

// Take an involution on {0,1,...,n-1} with size=n+num_fix_points
// and turn it into an involution on {0,1,...,size-1} with fix points
// at values fix_point_X (for X=1,...,num_fix_points

void reconstruct_and_analyse(int n, int size, int * S)
{
	int i;

	if (num_fix_points)
	{
		for (i=0;i<n;i++)
		{
			if (fix_point_1 >= 0) if (S[i] >= fix_point_1) S[i]++;
			if (fix_point_2 >= 0) if (S[i] >= fix_point_2) S[i]++;
			if (fix_point_3 >= 0) if (S[i] >= fix_point_3) S[i]++;
			if (fix_point_4 >= 0) if (S[i] >= fix_point_4) S[i]++;
		}
	}
	
	for (i=0;i<n;i+=2)
	{
		sbox[S[i]  ] = S[i+1];
		sbox[S[i+1]] = S[i  ];
		sinv[S[i]  ] = S[i+1];
		sinv[S[i+1]] = S[i  ];
	}

	if (num_fix_points)
	{
		if (fix_point_1 >= 0) sbox[fix_point_1] = fix_point_1;
		if (fix_point_2 >= 0) sbox[fix_point_2] = fix_point_2;
		if (fix_point_3 >= 0) sbox[fix_point_3] = fix_point_3;
		if (fix_point_4 >= 0) sbox[fix_point_4] = fix_point_4;
	}

#if 0
	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");
	}

	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

	generate_inversesbox(sinv, sbox, 16);
	process_one_sbox(sbox, sinv, 4);
}

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

void create_empty_involution(int n, int * s)
{
	for (int i=0;i<n;i++)
		s[i] = -1;
}

void output_involution(int n, int * s)
{
	int i;
	for (i=0;i<n;i+=2)
		printf("(%d,%d)", s[i], s[i+1]);
	printf("\n");
}

// A set of bijections
//   from S = {0, 1, .., 2k-1} to T = {1,..,2k+1} \ i , with i in {1,2,..,2k+1}.
// Just add 1 or 2 to map S to T while skipping the value i

int B(int i, int x)
{
	return (x+1 >= i) ? x+2 : x+1;
}

// n : cardinality of the final set
// S : current set (initial value = empty if (|S| = n) then

// Array S[] fatto così:
// swap S[0] with S[1]
// swap S[2] with S[3] (if S[2] and S[3] != -1)
// swap S[4] with S[5] (if S[4] and S[5] != -1)
// and so on

// n = max depth = size-0,2 or 4.

void fpfi(int n, int size, int * S)
{
	int Sp[n]; // create new one.
	int i,j,k;

	if (S[n-2] != -1) // finito
	{
		reconstruct_and_analyse(n, size, S);
		return;
	}

	for (i=1; i<n; i++)
	{
		// S′ = {{0,i}}
		create_empty_involution(n, Sp);
		Sp[0] = 0;
		Sp[1] = i;

		// forall {e, e′} in S, S′ = S′ U {{Bi(e),Bi(e′)}} 
		for (j=0;j+4<=n;j+=2)
		{
			// the pairs {e,e'} are {S[j],S[j+1]} in S with j even
			if (S[j] == -1) break;
			if ((Sp[j+2] = B(i,S[j  ])) >= n) return;
			if ((Sp[j+3] = B(i,S[j+1])) >= n) return;
		}
		fpfi(n,size,Sp);
	}
}