polymain.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,0x0000000,0x0000000,0x0000000,0x000000B1};  /*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_counte[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);
	}
	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; ke[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 1) return (0);
		for (i=0; i