\Encryption\polynominal
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
#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;
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--)
{
{
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				 *
*********************************************************************/

FIELD2N *x, *f, y[2];
{
FIELD2N	c, z, dummy;
FIELD2N test1, test2;
INDEX	i, j, k;

/*  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];
for (i = WORDSIZE/2; i > 0; i >>= 1)
{
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);
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)
{
for (i = WORDSIZE/2; i > 0; i >>= 1)
{
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);
{
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);
}
*/```