eliptic_poly.c

多项式的计算

/*  eliptic_poly.c  */
/************************************************************************************
*																					*
*		Polynomial version of elliptic curve functions.  essentially the same as	*
*  optimal normal basis but uses polynomial functions.  Biggest difference is that	*
*  once poly_prime is picked (see polymain.c) you need to compute the "T matrix" to	*
*  help embed data onto a curve by solving a quadratic equation.					*
*																					*
*							Author = mike rosing									*
*							 date  = June 27, 1997									*
*																					*
************************************************************************************/

#include <stdio.h>
#include "field2n.h"
#include "poly.h"
#include "eliptic.h"

/*  Global parameters for polynomial math.  poly_prime is the irreducible polynomial
	for the Galois field arithmetic.  Tmatrix is used to solve quadratic equations
	and Smatrix is used to compute square roots.
*/

extern	FIELD2N poly_prime;
FIELD2N	Tmatrix[NUMBITS], Smatrix[NUMBITS], Trace_Vector;
INDEX	Null_Row;

/*  Subroutine to invert a binary matrix.  The method is brute force but simple.
	Start with a diagonal matrix in I, and for each row operation in mat_in that
	eliminates non diagonal terms, do the same in I.  The result is that mat_in 
	= 1 and I = 1/mat_in (transpose).  
	If the input transpose is set to any value, we then have to transpose I back
	into mat_out to get the correct result.  If transpose is zero, I is copied to
	mat_out and the value of error is returned to indicate any null row.
	
	Returns 0 if matrix invertible, row number of zero column if not invertible.
*/

INDEX poly_matrix_invert( mat_in, mat_out)
FIELD2N mat_in[NUMBITS], mat_out[NUMBITS];
{
	INDEX	row, col, i, j, rowdex, found, error;
	ELEMENT	src_mask, dst_mask;
	FIELD2N	dummy, Imatrix[NUMBITS];

/*  Create Imatrix as diagonalized  */

	for ( row=0; row<NUMBITS; row++)
	{
		null( &Imatrix[row]);
		Imatrix[row].e[NUMWORD - row/WORDSIZE] = 1L << (row%WORDSIZE);
	}
	error = 0;  /* hope this is return value  */
	
/*  Diagonalize input matrix.  Eliminate all other bits in each column. 
	First find a column bit that is set, and swap with diagonal
	row if needed.
*/

	for ( row = 0; row < NUMBITS; row ++)
	{
		rowdex = NUMWORD - row/WORDSIZE;
		src_mask = 1L << (row % WORDSIZE);

/*  find a row of input matrix which has col = row bit set.
	First see if we get lucky, then do search.
*/
		found = 0;
		if ( !(mat_in[row].e[rowdex] & src_mask))
		{
			for (j = row+1; j<NUMBITS; j++)
			{
				if ( mat_in[j].e[rowdex] & src_mask)  /*  found one, swap rows  */
				{
					copy( &Imatrix[j], &dummy);
					copy( &Imatrix[row], &Imatrix[j]);
					copy( &dummy, &Imatrix[row]);
					copy( &mat_in[j], &dummy);
					copy( &mat_in[row], &mat_in[j]);
					copy( &dummy, &mat_in[row]);
					found = 1;
					break;
				}
			}
		}
		else found = 1;
		
/*  eliminate all other terms in this column  */

		if (found)
		{
			for ( i=0; i<NUMBITS; i++)
			{
				if ( i == row) continue;
				if ( mat_in[i].e[rowdex] & src_mask)
				{
					SUMLOOP(j)
					{
						Imatrix[i].e[j] ^= Imatrix[row].e[j];
						mat_in[i].e[j] ^= mat_in[row].e[j];
					}
				}
			}
		}			/*  end column eliminate  */
		else error = row;
	}				/*  end diagonalization   */


/*  Next step is to transpose diagonalized matrix.  This completes the 
	inversion process.  Clear Tmatrix to begin with. */
	
	for (row=0; row<NUMBITS; row++) null(&mat_out[row]);
	for (col = 0; col<NUMBITS; col++)
	{
		j = NUMWORD - col/WORDSIZE;
		src_mask = 1L << (col % WORDSIZE);
		for (row = 0; row<NUMBITS; row++)
		{
			if (Imatrix[row].e[j] & src_mask)
			{
				i = NUMWORD - row/WORDSIZE;
				dst_mask = 1L << (row  % WORDSIZE);
				mat_out[col].e[i] |= dst_mask;
			}
		}
	}
	return(error);
}

/*  initialize polynomial math for elliptic curve calculations.  
	
	Creates Tmatrix and Smatrix to help solve quadratic equations as
	well as Trace_Vector.  The Smatrix is only invoked to compute the
	square root of a polynomial modulo poly_prime.  The Tmatrix is a
	set of basis vectors which are summed assuming a solution to the
	quadratic exists.  The Trace_Vector is used to determine if the
	solution exists.  If an error occurs in creating the Smatrix the
	routine returns the row it could not diagonalize, otherwise the
	routine returns 0 for no errors.
*/

INDEX init_poly_math()
{
	INDEX	i, j, error, k, sum, row, rowdex, found, nulldex;
	ELEMENT	src_mask;
	FIELD2N	x, c, dummy;
	FIELD2N Trace[2*NUMBITS], Tz[NUMBITS], T2[NUMBITS];

/*  Create Trace_Vector for computing the trace in polynomial basis.  
	Given any input c = c_m*x^m + ... + c_1*x + c_0  the Trace_Vector
	will mask off and sum the correct coefficients of a_i.  The
	following method was suggested by Richard Pinch ([e-mail])
*/

/*  step 1 Build a list of powers of x modulo poly_prime from 0 to 2n-2 */

	null( &Trace[0]);
	Trace[0].e[NUMWORD] = 1L;
	for (i=1; i<2*NUMBITS-1; i++)
		mul_x_mod( &Trace[i-1], &poly_prime, &Trace[i]);

/*  step 2: Sum diagonals of nxn matrix using each row as a starting
	point.  Each sum amounts to the trace vector coefficient for that
	power of x.
*/

	null( &Trace_Vector);
	for (i=0; i<NUMBITS; i++)
	{
		sum = 0;
		j = NUMWORD;
		src_mask = 1;
		for ( k=i; k<i+NUMBITS; k++)
		{
			if ( Trace[k].e[j] & src_mask) sum ^= 1;
			src_mask <<= 1;
			if (!src_mask)
			{
				src_mask = 1;
				j--;
			}
		}
		if (sum) Trace_Vector.e[NUMWORD - i/WORDSIZE] |= 1L << (i % WORDSIZE);
	}
	
/*  Next compute Tz = z^2 + z matrix.  We already have every power of x in
	the Trace matrix, so each row of Tz is just two rows from Trace.  Do
	partial inversion to get Tmatrix.
*/
	for ( i=0; i<NUMBITS; i++)
	{
		SUMLOOP(j) 
		{
			Tz[i].e[j] = Trace[i].e[j] ^ Trace[i<<1].e[j];
			T2[i].e[j] = Trace[i<<1].e[j];
		}
	}
	
/*  Invert T2 matrix to get special square root matrix (called Smatrix)  */

	if (error = poly_matrix_invert( T2, Smatrix))
	{
		printf("Can not invert square root matrix.  Null row = %d\n", error);
		return(error);
	}

/*  Create a set of basis vectors for use in finding roots of z^2 + z = c
	This is similar to the matrix inversion, but there is no transpose and
	the row order is different.  Exactly how this works is not clear to me,
	but Prof. Pinch gets it and it works.
*/

/*  Create Tmatrix as diagonalized  */

	for ( row=0; row<NUMBITS; row++)
	{
		null( &Tmatrix[row]);
		Tmatrix[row].e[NUMWORD - row/WORDSIZE] = 1L << ( row % WORDSIZE);
	}

/*  Semi diagonalize input matrix.  
	First find a column bit that is set, and swap with diagonal
	row if needed.  Then eliminate all bits below it.
*/

	Null_Row = 0;
	nulldex = 0;
	for ( row = 0; row < NUMBITS; row++)
	{
		rowdex = NUMWORD - row/WORDSIZE;
		src_mask = 1L << (row % WORDSIZE);
		
/*  find a row of input matrix which has col = row bit set.
	First see if we get lucky, then do search.
*/
		found = 0;
		if ( !(Tz[row].e[rowdex] & src_mask))
		{
			for (j = row+1; j<NUMBITS; j++)
			{
				if ( Tz[j].e[rowdex] & src_mask)  /*  found one, swap rows  */
				{
					copy( &Tmatrix[j], &dummy);
					copy( &Tmatrix[row], &Tmatrix[j]);
					copy( &dummy, &Tmatrix[row]);
					copy( &Tz[j], &dummy);
					copy( &Tz[row], &Tz[j]);
					copy( &dummy, &Tz[row]);
					found = 1;
			/*  keep track of original null row */
					if (row == nulldex) nulldex = j;
					break;
				}
			}
		}
		else found = 1;
		
/*  eliminate all other terms below diagonal in this column  */

		if (found)
		{
			for ( i=row+1; i<NUMBITS; i++)
			{
				if ( Tz[i].e[rowdex] & src_mask)
				{
					SUMLOOP(j)
					{
						Tmatrix[i].e[j] ^= Tmatrix[row].e[j];
						Tz[i].e[j] ^= Tz[row].e[j];
					}
				}
			}
		}			/*  end column eliminate  */
		else
		{		/*  mark null row and swap position with original  */
			Null_Row = row;
			copy( &Tmatrix[nulldex], &dummy);
			copy( &Tmatrix[row], &Tmatrix[nulldex]);
			copy( &dummy, &Tmatrix[row]);
			copy( &Tz[nulldex], &dummy);
			copy( &Tz[row], &Tz[nulldex]);
			copy( &dummy, &Tz[row]);
		}
	}				/*  end diagonalization   */

/*  finally eliminate all other terms above diagonal except for
	Null_Row.  Result is a set of basis vectors which converts
	c to z and z^2 + z = c.
*/

	for ( row = NUMBITS-1; row > 0; row--)
	{
		if (row == Null_Row) continue;
		
		rowdex = NUMWORD - row/WORDSIZE;
		src_mask = 1L << (row % WORDSIZE);

		for (i = row-1; i>=0; i--)
		{
			if (Tz[i].e[rowdex] & src_mask)
			{
				SUMLOOP(j)
				{
					Tmatrix[i].e[j] ^= Tmatrix[row].e[j];
					Tz[i].e[j] ^= Tz[row].e[j];
				}
			}
		}
	}			/*  end column eliminate  */
	return(0);
}

/*********************************************************************
*																	 *
*	solve a quadratic equation over an irreducible polynomial field. *
*   Enter with parameters for the equation y^2 + xy + f = 0, where y *
*  is the unknown, x is	usually the x coordinate of a point and 	 *
*  f = f(x) of a curve.  											 *
*																	 *
*	Assumes init_poly_math already run with Trace_Vector, Null_Row   *
*   and Tmatrix filled out.											 *
*																	 *
*	computes c = f/x^2 and tests if the trace condition is met.  If  *
*   not	returns an error code of 1.  Otherwise it means we can solve *
*   the reduced	equation z^2 + z + c = 0.  The change of variable is *
*   y = xz.															 *
*																	 *
*	returns error code zero, y[0] = xz, y[1] = y[0] + x				 *
*********************************************************************/

INDEX poly_quadratic( x, f, y)
FIELD2N *x, *f, y[2];
{
	FIELD2N	c, z, dummy;
	FIELD2N test1, test2;
	INDEX	i, j, k;
	ELEMENT	sum, mask;

/*  first check to see if x is zero  */

	sum = 0;
	SUMLOOP(i) sum |= x->e[i];
	
	if (sum)
	{
/*  compute c = x^-2 * f  */

		copy (x, &c);
		poly_mul(&c, x, &dummy);		/*  get x^2  */
		poly_inv(&dummy, &z);
		poly_mul( f, &z, &c);

/*  verify that trace of c = 0.  If not, no solution
	is possible.  */
	
		sum = 0;
		SUMLOOP(i) sum ^= c.e[i] & Trace_Vector.e[i];
        mask = ~0;
        for (i = WORDSIZE/2; i > 0; i >>= 1)
	    {
    	    mask >>= i;
        	sum = ((sum & mask) ^ (sum >> i));
        }
 
   	/* if last bit is set, there is no solution to equation.  
	This eliminates half the points in the field, which might
	make sense to a mathematician.
	*/
		if (sum)
		{
			null( &y[0]);
			null( &y[1]);
			return(1);
		}

/*  clear out null row bit.  that part of matrix will not work.  */

		j = NUMWORD - Null_Row / WORDSIZE;
		c.e[j] &= ~( 1L << (Null_Row % WORDSIZE));

/*  for every bit set in c, add that row of Tmatrix to solution.  */

		null( &z);
		mask = 1;
		j = NUMWORD;
		for (i=0; i<NUMBITS; i++)
		{
			if ( c.e[j] & mask)
				SUMLOOP(k) z.e[k] ^= Tmatrix[i].e[k];
			mask <<= 1;
			if ( !mask)
			{
				mask = 1;
				j--;
			}
		}
		
/*  compute final solution using input parameters.  */

		poly_mul( x, &z, &y[0]);	
		SUMLOOP(i) y[1].e[i] = y[0].e[i] ^ x->e[i];
		return(0);
	}
	else
	{
/*  x input was zero.  Return y = square root of f. Process
	involves ANDing each row of Smatrix with f and summing
	all bits to find coefficient for that power of x */

		null( &z);
		for (j = 0; j<NUMBITS; j++)
		{
			sum = 0;
			SUMLOOP(i) sum ^= Smatrix[j].e[i] & f->e[i];
			if (sum)
			{
        		mask = -1L;
        		for (i = WORDSIZE/2; i > 0; i >>= 1)
	       	 	{
    	    		mask >>= i;
        	    	sum = ((sum & mask) ^ (sum >> i));
            	}
	        } 
	        if (sum) z.e[NUMWORD - j/WORDSIZE] |= (1L << j%WORDSIZE);
		}
		copy( &z, &y[0]);
		copy( &z, &y[1]);
		return(0);
	}
}

/*  print field matrix.  an array of FIELD2N of length NUMBITS is assumed.  The
	bits are printed out as a 2D matrix with an extra space every 5 characters.
*/

void matrix_print( name, matrix, file)
FIELD2N matrix[NUMBITS];
char *name;
FILE *file;
{
	INDEX i,j;
	
	fprintf(file,"%s\n",name);
	for (i = NUMBITS-1; i>=0; i--)
	{
		if (i%5==0) fprintf(file,"\n");	/* extra line every 5 rows  */
		for (j = NUMBITS-1; j>=0; j--)
		{
			if ( matrix[i].e[NUMWORD - j/WORDSIZE] & (1L<<(j%WORDSIZE)))
				fprintf(file,"1");
			else fprintf(file,"0");
			if (j%5==0) fprintf(file," ");  /* extra space every 5 characters  */
		}
		fprintf(file,"\n"); 
	}
	fprintf(file,"\n");
}

/*  compute f(x) = x^3 + a_2*x^2 + a_6 for non-supersingular elliptic curves */

void poly_fofx( x, curv, f)
FIELD2N *x, *f;
CURVE *curv;
{
	FIELD2N 	x2, x3;
	INDEX		i;
	
	copy ( x, &x3);