eliptic_poly.c

ECC的C++源码

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

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

    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)
    {
    	printf("doubling point at infinity\n");
    	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: ");
	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");
}