randomly + static.c

Elliptic Curve Cryptography full working C code for encryption and decryprion.

	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  */
		}
		k /= 2;
		bitmask ^= (bitmask >> k);
	}
	return( lg2);
}


INDEX degreeof(t, dim)
ELEMENT *t;
INDEX	dim;
{
	INDEX	degree, k;
	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);
}

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;
	
	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);
}

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);
}


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;
	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);
}



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];
}

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);

	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]);
	}

	null( &x2r);
	x2r.e[NUMWORD] = 2;
	for ( r=1; r <= deg_v/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);
}



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


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  */


	for ( row = 0; row < NUMBITS; row ++)
	{
		rowdex = NUMWORD - row/WORDSIZE;
		src_mask = 1L << (row % WORDSIZE);
		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   */


	
	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);
}

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];



/*  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]);

	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);
	}
	

	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 Tmatrix as diagonalized  */

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



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

		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   */



	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);
}


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);


	
		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 (sum)
		{
			null( &y[0]);
			null( &y[1]);
			return(1);
		}


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


		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
	{

		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);
	}
}

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] );
}


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);
}

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



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)
	{
		printf("summing point at infinity\n");
		copy_point( p2, p3);
		return;
	}
	check = 0;
	SUMLOOP(i) check |= p2->x.e[i] | p2->y.e[i];
	if (!check) 
	{
		printf("summing point at infinity\n");
		copy_point( p1, p3);
		return;
	}
	
/*  compute theta = (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/x_1 = theta */
    poly_mul(&theta, &theta, &theta2);   /* then theta^2  */

/*  with theta and theta^2, compute x_3  */