/*  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 
#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 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; ie[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; je[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);
	poly_mul( x, &x3, &x2);       /*  get x^2  */
	if (curv->form) poly_mul ( &x2, &curv->a2, f);
	else null( f);
	poly_mul( x, &x2, &x3);       /*  get x^3  */
	SUMLOOP (i) f->e[i] ^= ( x3.e[i] ^ curv->a6.e[i] );
}

/*  embed data onto a curve.
	Enter with data, curve, ELEMENT offset to be used as increment, and
	which root (0 or 1).
	Returns with point having data as x and correct y value for curve.
	Will use y[0] for last bit of root clear, y[1] for last bit of root set.
	if ELEMENT offset is out of range, default is 0.
*/

void poly_embed( data, curv, incrmt, root, pnt)
FIELD2N	*data;
CURVE	*curv;
INDEX	incrmt, root;
POINT	*pnt;
{
	FIELD2N		f, y[2];
	INDEX		inc = incrmt;
	INDEX		i;
	
	if ( (inc < 0) || (inc > NUMWORD) ) inc = 0;
	copy( data, &pnt->x);
	poly_fofx( &pnt->x, curv, &f);
	while (poly_quadratic( &pnt->x, &f, y))
	{
		pnt->x.e[inc]++;
		poly_fofx( &pnt->x, curv, &f);
	}
	copy ( &y[root&1], &pnt->y);
}

/****************************************************************************
*                                                                           *
*   Implement elliptic curve point addition for polynomial basis form.  	*
*  This follows R. Schroeppel, H. Orman, S. O'Mally, "Fast Key Exchange with*
*  Elliptic Curve Systems", CRYPTO '95, TR-95-03, Univ. of Arizona, Comp.   *
*  Science Dept.                                                            *
****************************************************************************/

void poly_esum (p1, p2, p3, curv)
POINT   *p1, *p2, *p3;
CURVE   *curv;
{
    INDEX   i;
    FIELD2N  x1, y1, theta, onex, theta2;
    ELEMENT  check;
	
/*  check if p1 or p2 is point at infinity  */

	check = 0;
	SUMLOOP(i) check |= p1->x.e[i] | p1->y.e[i];
	if (!check)
	{
		copy_point( p2, p3);
		return;
	}
	check = 0;
	SUMLOOP(i) check |= p2->x.e[i] | p2->y.e[i];
	if (!check) 
	{
		copy_point( p1, p3);
		return;
	}
	
/*  compute lambda = (y_1 + y_2)/(x_1 + x_2)  */

    null(&x1);
    null(&y1);
    check = 0;
    SUMLOOP(i) 
    {
		x1.e[i] = p1->x.e[i] ^ p2->x.e[i];
		y1.e[i] = p1->y.e[i] ^ p2->y.e[i];
		check |= x1.e[i];
    }
    if (!check)   /* return point at infinity  */
    {
/*    	printf("input to elliptic sum has duplicate x values\n"); */
    	null(&p3->x);
    	null(&p3->y);
    	return;
    }
    poly_inv( &x1, &onex);
    poly_mul( &onex, &y1, &theta);  /*  compute y_1/lambda */
    poly_mul(&theta, &theta, &theta2);   /* then lambda^2  */

/*  with lambda and lamda^2, compute x_3  */

    if (curv->form)
		SUMLOOP (i)
	    	p3->x.e[i] = theta.e[i] ^ theta2.e[i] ^ x1.e[i] ^ curv->a2.e[i];
    else
		SUMLOOP (i)
	    	p3->x.e[i] = theta.e[i] ^ theta2.e[i] ^ x1.e[i];

/*  next find y_3  */

    SUMLOOP (i) x1.e[i] = p1->x.e[i] ^ p3->x.e[i];
    poly_mul( &x1, &theta, &theta2);
    SUMLOOP (i) p3->y.e[i] = theta2.e[i] ^ p3->x.e[i] ^ p1->y.e[i];
}

/*  elliptic curve doubling routine for Schroeppel's algorithm over polymomial
    basis.  Enter with p1, p3 as source and destination as well as curv
    to operate on.  Returns p3 = 2*p1.
*/

void poly_edbl (p1, p3, curv)
POINT *p1, *p3;
CURVE *curv;
{
    FIELD2N  x1, y1, theta, theta2, t1;
    INDEX   i;
    ELEMENT  check;
    
    check = 0;
    SUMLOOP (i) check |= p1->x.e[i];
    if (!check)
    {
    	null(&p3->x);
    	null(&p3->y);
    	return;
    }

/*  first compute theta = x + y/x  */

    poly_inv( &p1->x, &x1);
    poly_mul( &x1, &p1->y, &y1);
    SUMLOOP (i) theta.e[i] = p1->x.e[i] ^ y1.e[i];

/*  next compute x_3  */

    poly_mul( &theta, &theta, &theta2);
    if(curv->form)
		SUMLOOP (i) p3->x.e[i] = theta.e[i] ^ theta2.e[i] ^ curv->a2.e[i];
    else
		SUMLOOP (i) p3->x.e[i] = theta.e[i] ^ theta2.e[i];

/*  and lastly y_3  */

    theta.e[NUMWORD] ^= 1;			/*  theta + 1 */
    poly_mul( &theta, &p3->x, &t1);
    poly_mul( &p1->x, &p1->x, &x1);
    SUMLOOP (i) p3->y.e[i] = x1.e[i] ^ t1.e[i];
}

/*  subtract two points on a curve.  just negates p2 and does a sum.
    Returns p3 = p1 - p2 over curv.
*/

void poly_esub (p1, p2, p3, curv)
POINT   *p1, *p2, *p3;
CURVE   *curv;
{
    POINT   negp;
    INDEX   i;

    copy ( &p2->x, &negp.x);
    null (&negp.y);
    SUMLOOP(i) negp.y.e[i] = p2->x.e[i] ^ p2->y.e[i];
    poly_esum (p1, &negp, p3, curv);
}

/*  need to move points around, not just values.  Optimize later.  */

void copy_point (p1, p2)
POINT *p1, *p2;
{
	copy (&p1->x, &p2->x);
	copy (&p1->y, &p2->y);
}

/*  Routine to compute kP where k is an integer (base 2, not normal basis)
	and P is a point on an elliptic curve.  This routine assumes that K
	is representable in the same bit field as x, y or z values of P.  Since
	the field size determines the largest possible order, this makes sense.
    Enter with: integer k, source point P, curve to compute over (curv) 
    Returns with: result point R.

  Reference: Koblitz, "CM-Curves with good Cryptografic Properties", 
	Springer-Verlag LNCS #576, p279 (pg 284 really), 1992
*/

void  poly_elptic_mul(k, p, r, curv)
FIELD2N	*k;
POINT	*p, *r;
CURVE	*curv;
{
	char		blncd[NUMBITS+1];
	INDEX		bit_count, i;
	ELEMENT		notzero;
	FIELD2N		number;
	POINT		temp;

/*  make sure input multiplier k is not zero.
	Return point at infinity if it is.
*/
	copy( k, &number);
	notzero = 0;
	SUMLOOP (i) notzero |= number.e[i];
	if (!notzero)
	{
		null (&r->x);
		null (&r->y);
		return;
	}

/*  convert integer k (number) to balanced representation.
	Called non-adjacent form in "An Improved Algorithm for
	Arithmetic on a Family of Elliptic Curves", J. Solinas
	CRYPTO '97. This follows algorithm 2 in that paper.
*/
	bit_count = 0;
	while (notzero)
	{
	/*  if number odd, create 1 or -1 from last 2 bits  */
	
		if ( number.e[NUMWORD] & 1 )
		{
			blncd[bit_count] = 2 - (number.e[NUMWORD] & 3);
			
	/*  if -1, then add 1 and propagate carry if needed  */
			
			if ( blncd[bit_count] < 0 )
			{
				for (i=NUMWORD; i>=0; i--)
				{
					number.e[i]++;
					if (number.e[i]) break;
				}
			}
		}
		else
			blncd[bit_count] = 0;
	
	/*  divide number by 2, increment bit counter, and see if done  */
	
		number.e[NUMWORD] &= ~0 << 1;
		rot_right( &number);
		bit_count++;
		notzero = 0;
		SUMLOOP (i) notzero |= number.e[i];
	}
		
/*  now follow balanced representation and compute kP  */

	bit_count--;
	copy_point(p,r);		/* first bit always set */
	while (bit_count > 0) 
	{
	  poly_edbl(r, &temp, curv);
	  bit_count--;
	  switch (blncd[bit_count]) 
	  {
	     case 1: poly_esum (p, &temp, r, curv);
				 break;
	     case -1: poly_esub (&temp, p, r, curv);
				  break;
	     case 0: copy_point (&temp, r);
	   }
	}
}

/*void print_field( string, x)
char *string;
FIELD2N *x;
{
	INDEX i;
	
	printf("%s\n",string);
	SUMLOOP(i) printf("%8x ",x->e[i]);
	printf("\n");
}

void print_point( string, point)
char *string;
POINT *point;
{
	INDEX i;
	
	printf("%s\n",string);
	printf("x: ");
/*	printf("%s  ",string);
	SUMLOOP(i) printf("%8x ",point->x.e[i]);
	printf("\n");
	printf("y: ");
	SUMLOOP(i) printf("%8x ",point->y.e[i]);
	printf("\n");
}

void print_curve( string, curv)
char *string;
CURVE *curv;
{
	INDEX i;
	
	printf("%s\n",string);
	printf("form: %d\n",curv->form);
	if (curv->form)
	{
		printf("a2: ");
		SUMLOOP(i) printf("%lx ",curv->a2.e[i]);
		printf("\n");
	}
	printf("a6: ");
	SUMLOOP(i) printf("%lx ",curv->a6.e[i]);
	printf("\n\n");
}


main()
{
	FIELD2N	t1, t2, test;
	FIELD2N q, r, y, x, y2, xy, g[3];
	INDEX	i, error, j, order, k, m, n;
	ELEMENT index, check;
	FILE *del;
	CURVE  crv;
	POINT   p2, p3, p4, p5, p6, p7;
	char curve_string[80];
	
/*	del = fopen("curves_points_5","w");
	if (!del) return(0);
*/
/*	if (!irreducible(&poly_prime)) return(0);
	print_field("poly_prime = ", &poly_prime);
	
	if (error = init_poly_math())
	{
		printf("Can't initialize S matrix, row = %d\n", error);
		return(-1);
	}

	poly_gf8( &poly_prime, g);
	print_field(" g_0", &g[0]);
	print_field("g_1", &g[1]);
	print_field("g_2", &g[2]);
	
/*  check solutions actually work  */

/*	copy (&g[0], &t1);
	print_field("g[0]", &t1);
	
	poly_mul( &g[0], &t1, &t2);
	print_field("g[0]^2", &t2);
	poly_mul( &g[0], &t2, &t1);
	print_field("g[0]^3", &t1);
	SUMLOOP (i) test.e[i] = g[0].e[i] ^ t1.e[i];
	print_field(" one = g[0] + g[0]^3", &test);
	return;
	crv.form = 0;
	null(&crv.a2);
	null(&crv.a6);
/*	crv.a2.e[NUMWORD] = 5;*/
/*	crv.a6.e[NUMWORD] = 9;
	null(&test);
	test.e[NUMWORD] = 5;
	poly_embed( &test, &crv, NUMWORD, 0, &p2);
	
/*  check that point is in fact on curve  */

/*	copy(&p2.y, &y);
	copy(&p2.x, &x);
	print_point("for point", &p2);
	poly_mul( &p2.y, &y, &y2);
	poly_mul( &y, &x, &xy);
	SUMLOOP(i) r.e[i] = y2.e[i] ^ xy.e[i];
	poly_fofx( &x, &crv, &q);
	SUMLOOP(i) test.e[i] = r.e[i] ^ q.e[i];
	print_field("rhs + lhs =",&test);

/*  check that elliptic multiply works  */

/*	print_point(" base point P", &p2);
	poly_edbl( &p2, &p3, &crv);
	print_point(" 2P ", &p3);
	poly_esum( &p2, &p3, &p4, &crv);
	print_point(" 3P ", &p4);
	poly_esum( &p2, &p4, &p5, &crv);
	print_point(" 4P ", &p5);
	poly_esum( &p2, &p5, &p6, &crv);
	print_point(" 5P ", &p6);
	
	poly_esum( &p2, &p6, &p4, &crv);
	print_point(" 6P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point(" 7P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point(" 8P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point(" 9P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point("10P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point("11P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point("12P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point("13P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point("14P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point("15P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point("16P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point("17P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point("18P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point("19P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point("20P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point("21P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point("22P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point("23P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point("24P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point("25P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point("26P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point("27P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point("28P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point("29P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point("30P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point("31P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point("32P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point("33P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point("34P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point("35P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point("36P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point("37P ", &p6);
	poly_esum( &p2, &p6, &p4, &crv);
	print_point("38P ", &p4);
	poly_esum( &p2, &p4, &p6, &crv);
	print_point("39P ", &p6);
	
	for (i=1; i<40; i++)
	{
		null( &t1);
		t1.e[NUMWORD] = i;
		poly_elptic_mul( &t1, &p2, &p7, &crv);
		sprintf( curve_string, "%d*P", i);
		print_point( curve_string, &p7);
	}
	scanf("%c",&i);
}
*/