/****************************************************************
*																*
*		periferal functions useful for setting up polynomial	*
*	math but not called routinely.								*
*																*
****************************************************************/

#include "field2n.h"
#include "poly.h"

/*  for given irreducible polynomial solve g^3 = g + 1.
	Useful for GF(8^n) Koblitz curves.
	Trick is to lineraize equation into form g^4 + g^2 + g = 0.
	This is linear because (a + b)^(2^n) = a^(2^n) + b^(2^n).
	There are four solutions: one is zero, two are found as
	independent vectors and the third is their sum.
	
	Input:  pointers to irreducble polynomial and g[3]
	Output: g[0], g[1], and g[3] = g[0] + g[1] filled in
*/

void poly_gf8( prime, g)
FIELD2N  *prime, g[3];
{
	FIELD2N		gamma_matrix[NUMBITS], solve[NUMBITS];
	FIELD2N		power_table[4*NUMBITS], temp;
	INDEX		row, column, i, j, found, vector;
	ELEMENT		tobit, frombit, bit, null_check;
	
/*  step 1:  compute all powers of x modulo prime polynomial
	with simple function */

	null( &power_table[0]);
	power_table[0].e[NUMWORD] = 1L;
	for (row = 1; row < 4*NUMBITS; row++)
		mul_x_mod( &power_table[row-1], prime, &power_table[row]);
	
/*  step 2:  sum powers of rows to create g^4 + g^2 + g 
	coefficients matrix.
*/
	for( row=0; row < NUMBITS; row++)
	{
		copy( &power_table[row], &gamma_matrix[row]);
		SUMLOOP (i) gamma_matrix[row].e[i] ^=
				power_table[row<<1].e[i] ^ power_table[row<<2].e[i];
	}

/*  step 3:  transpose matrix and work with single powers of x */

	for ( row=0; row < NUMBITS; row++)  null( &solve[row]);
	for ( row=0; row < NUMBITS; row++)
	{
		bit = 1L << (row % WORDSIZE);
		i = NUMWORD - (row/WORDSIZE);
		frombit = 1;
		j = NUMWORD;
		for ( column = 0; column < NUMBITS; column++)
		{
			if (gamma_matrix[row].e[j] & frombit)
						solve[column].e[i] |= bit;
			frombit <<= 1;
			if ( !frombit)
			{
				frombit = 1;
				j--;
			}
		}
	}

/*  step 4: lower diagonalize solve matrix.
	2 columns will go null, search for null rows instead
	of diagnalizing. */

	vector = 0;		/*  set up solution space  */
	null( &g[0]);
	null( &g[1]);

	for (column = NUMBITS - 1; column > 0; column--)
	{
		bit = 1L << (column % WORDSIZE);
		i = NUMWORD - column/WORDSIZE;
		if ( !(bit & solve[column].e[i]) )
		{
		/*  go look for a bit set in this column  */
		
			found = 0;
			for (row = column - 1; row >= 0; row--)
			{
				if ( bit & solve[row].e[i])
				{
					copy ( &solve[row], &temp);
					copy ( &solve[column], &solve[row]);
					copy ( &temp, &solve[column]);
					found = 1;
					break;
				}
			}
		}
		else found = 1;
	/*  and eliminate any set bits below it  */
	
		if (found)
		{
			for ( row = column-1; row >=0; row--)
			{
				if (solve[row].e[i] & bit)
					SUMLOOP (j) solve[row].e[j] ^= solve[column].e[j];
			}
		}
		else
		{
	/*  if null column, see if we have null row to match it.  */
			
			null_check = 0;
			SUMLOOP (j) null_check |= solve[column].e[j];
			if ( null_check)
			{
			/*  search for a null row to put here  */
			
				for ( row = column-1; row >=0; row--)
				{
					null_check = 0;
					SUMLOOP (j) null_check |= solve[row].e[j];
					if ( !null_check) 
					{
						copy ( &solve[row], &temp);
						copy ( &solve[column], &solve[row]);
						copy ( &temp, &solve[column]);
						break;
					}
				}
			}
			
	/*  mark this vector with column bit, and use next vector */
	
			g[vector].e[i] |= bit;
			vector++;
		}
	}

/*  last row may be null and not marked.  if INDEX vector < 2,
	set g[1] = 1
*/
	if (vector < 2) g[1].e[NUMWORD] = 1;
	
/*  find two solution vectors by back solving matrix.  use bits
	previoiusly solved to find bits not yet solved.  Initial vectors
	come from assuming g0 and g1 are variables.  
*/

	for ( row=1; row= 0; column--)
			{
				frombit = 1L << (column % WORDSIZE);
				i = NUMWORD - column/WORDSIZE;
				if ( solve[row].e[i] & frombit)
				{
					if (g[0].e[i] & frombit) g[0].e[j] ^= tobit;
					if (g[1].e[i] & frombit) g[1].e[j] ^= tobit;
				}
			}
		}
	}

/*  last step:  g[2] = g[1] + g[0]  */

	SUMLOOP (i) g[2].e[i] = g[0].e[i] ^ g[1].e[i];
}