polymain.c

ECC的C++源码

/*  This is the main test driver for polynomial math routines*/

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

/*  Define irreducible polynomial for field here.  Use table, or evaluate your
	own, but the math won't work if it' isn't really a prime polynomial.  
	NOTE: the value you pick must be consistent with NUMBITS.
*/

/*FIELD2N  poly_prime = {0x08000000,0x00000000,0x00000000,0x00000000,0x000000B1};  /*155*/
/*FIELD2N  poly_prime = {0x08000000,0x00000000,0x00000000,0x40000000,0x00000001};  /*155*/
/*FIELD2N poly_prime = {0x00000040,0x00000000,0x00000000,0x02000000,0x00000001}; /*134*/
/*FIELD2N poly_prime = {0x00000008,0x00000000,0x00000000,0x00000000,0x0000010D}; /*131*/
/*FIELD2N poly_prime = {0x00000004,0x00000000,0x00000000,0x00000000,0x00000009}; /*130*/
/*FIELD2N poly_prime = {0x00000002,0x00000000,0x00000000,0x00000000,0x00000021}; /*129*/
/*FIELD2N poly_prime = {0x00000001,0x00000000,0x00000000,0x00000000,0x00000087}; /*128*/
FIELD2N poly_prime = {0x00008000,0x00000000,0x00000000,0x00000401};  /* 111*/
/*FIELD2N  poly_prime = {0x20000000,0x00000000,0x00000005};	/*93*/
/*FIELD2N  poly_prime = {0x0000800,0x00000000,0x0000004b};*/	/*75*/
/*FIELD2N  poly_prime = {0x0000020,0x00000000,0x00000065};*/	/*69*/
/*FIELD2N  poly_prime = {0x00000002,0x00000000,0x00000001b};	/*65*/
/*FIELD2N  poly_prime = {0x00000001,0x00000000,0x00000001b};*/	/*64*/
/*FIELD2N poly_prime = {0x80000000,0x00000003};*/ /*63*/
/*FIELD2N poly_prime = {0x00000002,0x00002001}; /*33*/
/*FIELD2N poly_prime = {0x00000001,0x000000c5}; /*32*/
/*FIELD2N poly_prime = {0x80000009}; /*31*/
/*FIELD2N poly_prime = {0x800021};*/ /* 23*/
/*FIELD2N poly_prime = {0x100009};*/ /* 20*/
/*FIELD2N poly_prime = {0x20021};*/ /* 17*/
/*FIELD2N poly_prime = {0x1002B};*/ /* 16*/
/*FIELD2N poly_prime = {0x8003}; */ /* 15*/
/*FIELD2N poly_prime = {0x89};*/ /* 7*/
/*FIELD2N poly_prime = {0x43}; /*6*/
/*FIELD2N poly_prime = {0x29}; /*5*/
/*FIELD2N poly_prime = {0x13};*/ /*4*/

/*  rotation routines copied from normal basis.  Useful
	in general for high level operations.
*/

void rot_left(a)
FIELD2N *a;
{
        INDEX i;
        ELEMENT bit,temp;

        bit = (a->e[0] & UPRBIT) ? 1L : 0L;
        for (i=NUMWORD; i>=0; i--) {
           temp = (a->e[i] & MSB) ? 1L : 0L;
           a->e[i] = ( a->e[i] << 1) | bit;
           bit = temp;
        }
        a->e[0] &= UPRMASK;
}

void rot_right(a)
FIELD2N *a;
{
        INDEX i;
        ELEMENT bit,temp;

        bit = (a->e[NUMWORD] & 1) ? UPRBIT : 0L;
        SUMLOOP(i) {
           temp = ( a->e[i] >> 1)  | bit;
           bit = (a->e[i] & 1) ? MSB : 0L;
           a->e[i] = temp;
        }
        a->e[0] &= UPRMASK;
}

/*  Shift left one bit used by multiply routine.  Make inline for speed.  */

void mul_shift(a)
DBLFIELD *a;
{
	ELEMENT	*eptr, temp, bit;	/* eptr points to one ELEMENT at a time  */
	INDEX	i;
	
	eptr = &a->e[DBLWORD];		/* point to end, note: bigendian processing  */
	bit = 0;					/*  initial carry bit is clear */
	DBLLOOP (i)
	{
		temp = (*eptr << 1) | bit;		/* compute result as temporary  */
		bit = (*eptr & MSB) ? 1L : 0L;	/* get carry bit from shift  */
		*eptr-- = temp;					/* save new result  */
	}
}

/*  null out a FIELD2N variable.  Make inline for speed.  */

void null(a)
FIELD2N *a;
{
	INDEX i;
	
	SUMLOOP(i) a->e[i] = 0L;
}

void dblnull(a)
DBLFIELD *a;
{
	INDEX i;
	
	DBLLOOP(i) a->e[i] = 0L;
}

/*  copy one FIELD2N variable to another.  Make inline for speed.  */

void copy(from, to)
FIELD2N *from, *to;
{
	INDEX	i;
	
	SUMLOOP(i) to->e[i] = from->e[i];
}

/*  copy a FIELD2N variable into a DBLFIELD variable */

void sngltodbl(from, to)
FIELD2N *from;
DBLFIELD *to;
{
	INDEX i;
	
	dblnull (to);
	SUMLOOP(i) to->e[DBLWORD - NUMWORD + i] = from->e[i];
}

/*  copy the bottom portion of DBLFIELD to FIELD2N variable  */

void dbltosngl(from, to)
FIELD2N *to;
DBLFIELD *from;
{
	INDEX i;
	
	SUMLOOP(i) to->e[i] = from->e[DBLWORD - NUMWORD + i];
}

/*  Simple polynomial multiply. 
	Enter with 2 single FIELD2N variables to multiply.
	Result is double length DBLFIELD variable.
	user supplys all storage space, only pointers used here.
*/

void poly_mul_partial(a, b, c)
FIELD2N *a, *b;
DBLFIELD *c;
{
	INDEX i, bit_count, word;
	ELEMENT mask;
	DBLFIELD B;
	
/*  clear all bits in result  */

	dblnull(c);
	
/*  initialize local copy of b so we can shift it  */

	sngltodbl(b, &B);
	
/*  for every bit in 'a' that is set, add present B to result  */

	mask = 1;
	for (bit_count=0; bit_count<NUMBITS; bit_count++)
	{
		word = NUMWORD - bit_count/WORDSIZE;
		if (mask & a->e[word])
		{
			DBLLOOP(i) c->e[i] ^= B.e[i];
		}
		mul_shift( &B);			/* multiply copy of b by x  */
		mask <<= 1;				/* shift mask bit up  */
		if (!mask) mask = 1;	/* when it goes to zero, reset to 1  */
	}
}

/*  Shift right routine for polynomial divide.  Convert to inline for speed.
	Enter with pointer to DBLFIELD variable.
	shifts whole thing right one bit;
*/

void div_shift(a)
DBLFIELD *a;
{
	ELEMENT	*eptr, temp, bit;
	INDEX	i;
	
	eptr = (ELEMENT*) &a->e[0];
	bit = 0;
	DBLLOOP (i)
	{
		temp = (*eptr>>1) | bit;		/*  same as shift left but  */
		bit = (*eptr & 1) ? MSB : 0L;	/*  carry bit goes off other end  */
		*eptr++ = temp;
	}
}

/*  binary search for most significant bit within word */

INDEX log_2 (x)
ELEMENT x;
{
	ELEMENT ebit, bitsave, bitmask;
	INDEX	k, lg2;
	
	lg2 = 0;
	bitsave = x;					/* grab bits we're interested in.  */
	k = WORDSIZE/2;					/* first see if msb is in top half */
	bitmask = -1L<<k;				/* of all bits  */
	while (k)
	{
		ebit = bitsave & bitmask;	/* did we hit a bit?  */
		if (ebit)					/* yes  */
		{
			lg2 += k;				/* increment degree by minimum possible offset */
			bitsave = ebit;			/* and zero out non useful bits  */
		}
	/*  The following two lines are slick and tricky.  We are binary searching,
		so k is divided by 2.  It's also the base amount we can add at any given
		search point.  The mask continuously searches for upper blocks being set.
		Once found, the remaining lower blocks are zeroed, to make sure we don't
		use bits which aren't the most significant.
	*/
		k /= 2;
		bitmask ^= (bitmask >> k);
	}
	return( lg2);
}

/*  find most significant bit in multiple ELEMENT array.
	Enter with pointer to array (t) and number of elements (dim).
	Returns degree of polynomial == position count of most signicant bit set
*/

INDEX degreeof(t, dim)
ELEMENT *t;
INDEX	dim;
{
	INDEX	degree, k;

/*  This is generic routine for arbitrary length arrays. use ELEMENT
	pointer to find first non-zero ELEMENT.  We will add degree from some base,
	so initial base is at least one WORDSIZE smaller than maximum possible.
*/	

	degree = dim * WORDSIZE;
	for (k=0; k<dim; k++)					/* search for a non-zero ELEMENT  */
	{
		if (*t) break;
		degree -= WORDSIZE;
		t++;
	}

/* for inversion routine, need to know when all bits zero.  */

	if (!*t) return(-1);		

/*  find most significant bit in this element */
	
	degree += log_2( *t);
	return(degree);
}

/*  division of two polynomials.  Major use is reduction modulo 
	irreducible polynomials.
	Enter with pointers to top, bottom, quotient and remainder.
	  Assumes top is DBLFIELD and all other arrays are 
	  FIELD2N variables.
	Returns with top destroyed, and the following satisfied:
	  top = quotient * bottom + remainder.
*/

void poly_div(top, bottom, quotient, remainder)
DBLFIELD *top;
FIELD2N  *bottom, *quotient, *remainder;
{
	INDEX	deg_top;
	INDEX	deg_bot;
	INDEX	deg_quot;
	INDEX	bit_count;
	INDEX	i;
	INDEX	equot;
	ELEMENT	topbit, *tptr;
	DBLFIELD	shift;
	
/*  Step 1: find degree of top and bottom polynomials.  */

	deg_top = degreeof( top, DBLWORD);
	deg_bot = degreeof( bottom, NUMWORD);
	
/*  prepare for null return and check for null quotient */

	null (quotient);
	if (deg_top < deg_bot)
	{
		dbltosngl(top, remainder);
		return;
	}

/*  Step 2: shift bottom to align with top.  Note that there 
		are much more efficient ways to do this. */
	
	deg_quot = deg_top - deg_bot;
	bit_count = deg_quot + 1;
	sngltodbl(bottom, &shift);
	for (i = 0; i<deg_quot; i++)
		 mul_shift(&shift);

/* Step 3: create bit mask to check if msb of top is set  */

	topbit = 1L << (deg_top % WORDSIZE);
	tptr = (ELEMENT*) top + DBLWORD - deg_top/WORDSIZE;
	
/*  for each possible quotient bit, see if bottom can be subtracted
	(added, it's modulo 2!) from top. If it can, set that bit in
	quotient. If it can't, clear that bit in quotient.  Shift bottom
	right once and keep going for total bit_count.
*/

	while (bit_count)
	{
	
/*  Step 4: determine one bit of the quotient.  */

		if (*tptr & topbit)		/* is bit set in top?  */
		{
			DBLLOOP (i)			/* yes, subtract shift from top  */
				top->e[i] ^= shift.e[i];
				
	/* find word and bit in quotient to be set  */
	
			equot = 	NUMWORD - deg_quot/WORDSIZE;	
			quotient->e[equot] |= 1L << (deg_quot % WORDSIZE);	
		}

/*  Step 5: advance to the next quotient bit and perform the necessary
			shifts.  */
			
		bit_count--;			/*  number of bits is one more  */
		deg_quot--;				/*  than polynomial degree  */
		div_shift(&shift);		/*  divide by 2  */
		topbit >>= 1;			/*  move mask over one bit also  */
		if (!topbit)
		{
			topbit = MSB;		/*  reset mask bit to next word  */
			tptr++;				/*  when it goes to zero  */
		}
	}

/*  Step 6: return the remainder in FIELD2N size  */

	dbltosngl (top, remainder);
}

/*  Polynomial multiplication modulo poly_prime.  */

void poly_mul(a, b, c)
FIELD2N *a, *b, *c;
{
	DBLFIELD temp;
	FIELD2N  dummy;
	
	poly_mul_partial(a, b, &temp);
	poly_div(&temp, &poly_prime, &dummy, c);
}

/*  Polynomial inversion routine.  Computes inverse of polynomial field element
	assuming irreducible polynomial "poly_prime" defines the field. 
*/

void poly_inv(a, inverse)
FIELD2N *a, *inverse;
{
	FIELD2N	pk, pk1, pk2;
	FIELD2N rk, rk1;
	FIELD2N qk, qk1, qk2;
	INDEX	i;
	DBLFIELD rk2;
	
/*  initialize remainder, quotient and product terms  */

	sngltodbl (&poly_prime, &rk2);	
	copy( a, &rk1);
	null( &pk2);
	null( &pk1);
	pk1.e[NUMWORD] = 1L;
	null( &qk2);
	null( &qk1);
	qk1.e[NUMWORD] = 1L;

/*  compute quotient and remainder for Euclid's algorithm.
    when degree of remainder is < 0, there is no remainder, and we're done.
    At that point, pk is the answer.
*/
	null( &pk);
	pk.e[NUMWORD] = 1L;
	poly_div(&rk2, &rk1, &qk, &rk);
		
	while (degreeof(&rk, NUMWORD) >= 0)
	{
		poly_mul_partial( &qk, &pk1, &rk2);
		SUMLOOP(i) pk.e[i] = rk2.e[i+DBLWORD-NUMWORD] ^ pk2.e[i];
		
/*  set up variables for next loop */

		sngltodbl(&rk1, &rk2);
		copy( &rk, &rk1);
		copy( &qk1, &qk2);
		copy( &qk, &qk1);
		copy( &pk1, &pk2);
		copy( &pk, &pk1);
		poly_div( &rk2, &rk1, &qk, &rk);
	}
	copy( &pk, inverse);		/* copy answer to output  */
}

/*  polynomial greatest	common divisor routine.  Same as Euclid's algorithm.  */

void poly_gcd( u, v, gcd)
FIELD2N *u, *v, *gcd;
{
	DBLFIELD	top;
	FIELD2N		r, dummy, temp;
	
	sngltodbl( u, &top);
	copy( v, &r);
	
	while( degreeof( &r, NUMWORD) >= 0)
	{
		poly_div( &top, &r, &dummy, &temp);
		sngltodbl( &r, &top);
		copy( &temp, &r);
	}
	dbltosngl( &top, gcd);
}

/*  multiply a polynomial u by x modulo polynomial v.  Useful in
	several places.  Enter with u and v, returns polynomial w.
	w can equal to u, so this will work in place.
*/

void mul_x_mod( u, v, w)
FIELD2N *u, *v, *w;
{
	DBLFIELD mulx;
	INDEX    i, deg_v;
	
	deg_v = degreeof(v, NUMWORD);
	sngltodbl( u, &mulx);
	mul_shift( &mulx);  /*  multiply u by x  */
	dbltosngl( &mulx, w);
	if (w->e[ NUMWORD - (deg_v/WORDSIZE)] & ( 1L << (deg_v % WORDSIZE)))
		SUMLOOP (i) w->e[i] ^= v->e[i];
}

/*  check to see if input polynomial is irreducible.
	Returns 1 if irreducible, 0 if not.

	This method explained by Richard Pinch.  The idea is to use
	gcd algorithm to test if x^2^r - x has input v as a factor
	for all r <= degree(v)/2.  If there is any common factor,
	then v is not irreducible.  This works because x^2^r - x
	contains all possible irreducible factors of degree r, and
	because we only need to find the smallest degree such factor
	if one exists.
*/

INDEX irreducible( v)
FIELD2N *v;
{
	FIELD2N		vprm, gcd, x2r, x2rx, temp;
	FIELD2N		sqr_x[NUMBITS+1];
	INDEX		i, r, deg_v, k;
		
/*  check that gcd(v, v') = 1.  If not, then v not irreducible.  */ 

	SUMLOOP(i) vprm.e[i] = (v->e[i] >> 1) & DERIVMASK;
	poly_gcd( v, &vprm, &gcd);
	if (gcd.e[NUMWORD] > 1) return(0);
	for (i=0; i<NUMWORD; i++)  if (gcd.e[i]) return(0);

/*  find maximum power we have to deal with  */

	deg_v = degreeof(v, NUMWORD);

/*  create a vector table of powers of x^2k mod v.
	this will be used to compute square of x^2^r
*/

	null (&sqr_x[0]);
	sqr_x[0].e[NUMWORD] = 1;
	for (i=1; i<= deg_v; i++)
	{
		mul_x_mod( &sqr_x[i-1], v, &temp);
		mul_x_mod( &temp, v, &sqr_x[i]);
	}

/*  check that gcd( x^2^r - x, v) == 1 for all r <= degreeof(v)/2.
	set x^2^r = x for r = 0 to initialize.
*/
	null( &x2r);
	x2r.e[NUMWORD] = 2;
	for ( r=1; r <= deg_v/2; r++)
	{
/*  square x^2^r mod v.  We do this by seeing that s^2(x) = s(x^2).
	for each bit j set in x^2^(r-1) add in x^2j mod v to find x^2^r.
*/
		null( &x2rx);
		for (i=0; i <= deg_v; i++)
		{
			if ( x2r.e[NUMWORD - (i/WORDSIZE)] & ( 1L << (i%WORDSIZE) ))
				SUMLOOP(k) x2rx.e[k] ^= sqr_x[i].e[k];
		}
		
/*  save new value of x^2^r mod v and compute x^2^r - x mod v */

		copy( &x2rx, &x2r);
		x2rx.e[NUMWORD] ^= 2;

/*  is gcd( x^2^r - x, v) == 1?  if not, exit with negative result  */

		poly_gcd( &x2rx, v, &gcd);
		if (gcd.e[NUMWORD] > 1) return (0);
		for (i=0; i<NUMWORD; i++) if ( gcd.e[i]) return (0);
	}

/*  passed all tests, provably irreducible  */

	return (1);
}