/****** eliptic.c *****/ /************************************************************************ * * * Elliptic curves over Galois Fields. These routines find points * * on a curve, add and double points for type 1 optimal normal bases. * * * * For succint explanaition, see Menezes, page 92 * * * * This file modified 6/23/97 to work with generalized ONB and * * tunable field sizes. mgr * * * ************************************************************************/ #include #include "field2n.h" #include "eliptic.h" /************************************************************************ * Note that the following is obvious to mathematicians. I thought it * * was pretty cool when I discovered it myself, . * * * * Routine to solve quadradic equation. Enter with coeficients * * a and b and it returns solutions y[2]: y^2 + ay + b = 0. * * If Tr(b/a^2) != 0, returns y=0 and error code 1. * * If Tr(b/a^2) == 0, returns y[2] and error code 0. * * If solution fails, returns y=0 and error code 2. * * * * Algorithm used based on normal basis GF math. Since (a+b)^2 = * * a^2 + b^2 it follows that (a+b)^.5 = a^.5 + b^.5. Note that squaring* * is a left shift and rooting is a right shift in a normal basis. * * Transforming the source equation with y = ax and dividing by a^2 * * gives: * * x^2 + x + b/a^2 = 0 * * * * or x = x^.5 + (b/a^2)^.5 * * * * Let k_i = the ith significant bit of (b/a^2)^.5 and * * x_i = the ith significant bit of x. * * The above equation is equivelent to a bitwise representation as * * * * x_i = x_(i+1) + k_i * * or * * x(i+1) = x_i + k_i. * * * * Since both x and x+1 are solutions, and 1 is represented by all * * bits set in a normal basis, we can start with x_0 = 0 or 1 at our * * pleasure and use the recursion relation to discover every bit of x. * * The answer is then ax and ax+a returned in y[0] and y[1] respectively* * If the sum of x_(n-1) + k_(n-1) != x_0, returns error code 2 and * * y = 0. * * * * error code returns * * 0 y[0] and y[1] values * * 1 y[0] = y[1] = 0 * * 2 mathematicly impossible !!!! * * * ************************************************************************/ int opt_quadratic(a, b, y) FIELD2N *a, *b, *y; { INDEX i, l, bits; FIELD2N x, k, a2; ELEMENT r, t, mask; /* test for a=0. Return y = square root of b. */ r = 0; SUMLOOP(i) r |= a->e[i]; if (!r) { copy( b, &y[0]); rot_right( &y[0]); copy( &y[0], &y[1]); return(0); } /* find a^-2 */ opt_inv( a, &a2); rot_left(&a2); /* find k=(b/a^2)^.5 */ opt_mul( b, &a2, &k); rot_right(&k); r = 0; /* check that Tr(k) is zero. Combine all words first. */ SUMLOOP(i) r ^= k.e[i]; /* take trace of word, combining half of all the bits each time */ mask = -1L; for (bits = WORDSIZE/2; bits > 0; bits >>= 1) { mask >>= bits; r = ((r & mask) ^ (r >> bits)); } /* if not zero, return error code 1. */ if (r) { null(&y[0]); null(&y[1]); return(1); } /* point is valid, proceed with solution. mask points to bit i, which is known, in x bits previously found and k (=b/a^2)^.5. */ null(&x); mask = 1; for (bits=0; bits < NUMBITS ; bits++) { /* source long word could be different than destination */ i = NUMWORD - bits/WORDSIZE; l = NUMWORD - (bits + 1)/WORDSIZE; /* use present bits to compute next one */ r = k.e[i] & mask; t = x.e[i] & mask; r ^= t; /* same word, so just shift result up */ if ( l == i ) { r <<= 1; x.e[l] |= r; mask <<= 1; } else { /* different word, reset mask and use a 1 */ mask = 1; if (r) x.e[l] = 1; } } /* test that last bit generates a zero */ r = k.e[0] & UPRBIT; t = x.e[0] & UPRBIT; if ( r^t ) { null(&y[0]); null(&y[1]); return(2); } /* convert solution back via y = ax */ opt_mul(a, &x, &y[0]); /* and create complementary (z+1) solution y = ax + a */ null (&y[1]); SUMLOOP(i) y[1].e[i] = y[0].e[i] ^ a->e[i]; /* no errors, bye! */ return(0); } /* compute R.H.S. f(x) = x^3 + a2*x^2 + a6 curv.form = 0 implies a2 = 0, so no extra multiply. curv.form = 1 is the "twist" curve. */ void fofx(x, curv, f) FIELD2N *x, *f; CURVE *curv; { FIELD2N x2,x3; INDEX i; copy(x, &x2); rot_left(&x2); opt_mul(x, &x2, &x3); if (curv->form) opt_mul(&x2, &curv->a2, f); else null(f); SUMLOOP(i) f->e[i] ^= (x3.e[i] ^ curv->a6.e[i]); } /**************************************************************************** * * * Implement elliptic curve point addition for optimal normal 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. * * * * This version is faster for inversion processes requiring fewer * * multiplies than projective math version. For NUMBITS = 148 or 226 this * * is the case because only 10 multiplies are required for inversion but * * 15 multiplies for projective math. I leave it as a paper to be written * * [HINT!!] to propagate TR-95-03 to normal basis inversion. In that case * * inversion will require order 2 multiplies and this method would be far * * superior to projective coordinates. * ****************************************************************************/ void esum (p1, p2, p3, curv) POINT *p1, *p2, *p3; CURVE *curv; { INDEX i; FIELD2N x1, y1, theta, onex, theta2; /* compute theta = (y_1 + y_2)/(x_1 + x_2) */ null(&x1); null(&y1); SUMLOOP(i) { x1.e[i] = p1->x.e[i] ^ p2->x.e[i]; y1.e[i] = p1->y.e[i] ^ p2->y.e[i]; } opt_inv( &x1, &onex); opt_mul( &onex, &y1, &theta); copy( &theta, &theta2); rot_left(&theta2); /* with theta and theta^2, compute x_3 */ if (curv->form) SUMLOOP (i) p3->x.e[i] = theta.e[i] ^ theta2.e[i] ^ p1->x.e[i] ^ p2->x.e[i] ^ curv->a2.e[i]; else SUMLOOP (i) p3->x.e[i] = theta.e[i] ^ theta2.e[i] ^ p1->x.e[i] ^ p2->x.e[i]; /* next find y_3 */ SUMLOOP (i) x1.e[i] = p1->x.e[i] ^ p3->x.e[i]; opt_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 normal basis. Enter with p1, p3 as source and destination as well as curv to operate on. Returns p3 = 2*p1. */ void edbl (p1, p3, curv) POINT *p1, *p3; CURVE *curv; { FIELD2N x1, y1, theta, theta2, t1; INDEX i; /* first compute theta = x + y/x */ opt_inv( &p1->x, &x1); opt_mul( &x1, &p1->y, &y1); SUMLOOP (i) theta.e[i] = p1->x.e[i] ^ y1.e[i]; /* next compute x_3 */ copy( &theta, &theta2); rot_left(&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 */ one( &y1); SUMLOOP (i) y1.e[i] ^= theta.e[i]; opt_mul( &y1, &p3->x, &t1); copy( &p1->x, &x1); rot_left( &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 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]; 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. This is for simplicity, larger or smaller fields can be independently implemented. Enter with: integer k, source point P, curve to compute over (curv) and Returns with: result point R. Reference: Koblitz, "CM-Curves with good Cryptografic Properties", Springer-Verlag LNCS #576, p279 (pg 284 really), 1992 */ void 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) { edbl(r, &temp, curv); bit_count--; switch (blncd[bit_count]) { case 1: esum (p, &temp, r, curv); break; case -1: esub (&temp, p, r, curv); break; case 0: copy_point (&temp, r); } } } /* One is not what it appears to be. In any normal basis, "1" is the sum of all powers of the generator. So this routine puts ones to fill the number size being used in the address of the FIELD2N supplied. */ void one (place) FIELD2N *place; { INDEX i; SUMLOOP(i) place->e[i] = -1L; place->e[0] &= UPRMASK; }