eliptic_poly.c

多项式的计算

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