X-Git-Url: https://git.saurik.com/apple/security.git/blobdiff_plain/80e2389990082500d76eb566d4946be3e786c3ef..d8f41ccd20de16f8ebe2ccc84d47bf1cb2b26bbb:/Security/libsecurity_cryptkit/lib/elliptic.c diff --git a/Security/libsecurity_cryptkit/lib/elliptic.c b/Security/libsecurity_cryptkit/lib/elliptic.c new file mode 100644 index 00000000..3f45d199 --- /dev/null +++ b/Security/libsecurity_cryptkit/lib/elliptic.c @@ -0,0 +1,1437 @@ +/* Copyright (c) 1998,2011,2014 Apple Inc. All Rights Reserved. + * + * NOTICE: USE OF THE MATERIALS ACCOMPANYING THIS NOTICE IS SUBJECT + * TO THE TERMS OF THE SIGNED "FAST ELLIPTIC ENCRYPTION (FEE) REFERENCE + * SOURCE CODE EVALUATION AGREEMENT" BETWEEN APPLE, INC. AND THE + * ORIGINAL LICENSEE THAT OBTAINED THESE MATERIALS FROM APPLE, + * INC. ANY USE OF THESE MATERIALS NOT PERMITTED BY SUCH AGREEMENT WILL + * EXPOSE YOU TO LIABILITY. + *************************************************************************** + + elliptic.c - Library for Apple-proprietary Fast Elliptic + Encryption. The algebra in this module ignores elliptic point's + y-coordinates. + + Patent information: + + FEE is patented, U.S. Patents #5159632 (1992), #5271061 (1993), + #5463690 (1994). These patents are implemented + in various elliptic algebra functions such as + numer/denom_plus/times(), and in the fact of special + forms for primes: p = 2^q-k. + + Digital signature using fast elliptic addition, in + U. S. Patent #5581616 (1996), is implemented in the + signature_compare() function. + + FEED (Fast Elliptic Encryption) is patent pending (as of Jan 1998). + Various functions such as elliptic_add() implement the patent ideas. + + + Modification history since the U.S. Patent: + ------------------------------------------- + 10/06/98 ap + Changed to compile with C++. + 9 Sep 98 at Apple + cp->curveType optimizations. + Removed code which handled "unknown" curve orders. + elliptic() now exported for timing measurements. + 21 Apr 98 at Apple + Used inline platform-dependent giant arithmetic. + 20 Jan 98 at Apple + feemod now handles PT_MERSENNE, PT_FEE, PT_GENERAL. + Added calcGiantSizes(), rewrote giantMinBytes(), giantMaxShorts(). + Updated heading comments on FEE curve algebra. + 11 Jan 98 at Apple + Microscopic feemod optimization. + 10 Jan 98 at Apple + ell_odd, ell_even() Montgomery optimization. + 09 Jan 98 at Apple + ell_odd, ell_even() Atkin3 optimization. + 08 Jan 97 at Apple + Cleaned up some debugging code; added gsquareTime + 11 Jun 97 at Apple + Mods for modg_via_recip(), divg_via_recip() math + Deleted a lot of obsolete code (previously ifdef'd out) + Added lesserX1OrderJustify(), x1OrderPlusJustify() + Added binvg_cp(), avoiding general modg in favor of feemod + 05 Feb 97 at Apple + New optimized numer_plus(), denom_double(), and numer_times() + All calls to borrowGiant() and newGiant have explicit giant size + 08 Jan 97 at NeXT + Major mods to accomodate IEEE-style curve parameters. + New functions feepowermodg() and curveOrderJustify(); + elliptic_add(), elliptic(), signature_compare(), and + which_curve() completely rewritten. + 19 Dec 96 at NeXT + Added mersennePrimes[24..26] + 08 Aug 96 at NeXT + Fixed giant leak in elliptic_add() + 05 Aug 96 at NeXT + Removed dead code + 24 Jul 96 at NeXT + Added ENGINE_127_BITS dependency for use of security engine + 24 Oct 92 at NeXT + Modified new_public_from_private() + Created. + + + FEE curve algebra, Jan 1997. + + Curves are: + + y^2 = x^3 + c x^2 + a x + b + + where useful parameterizations for practical purposes are: + + Montgomery: a = 1, b = 0. (The original 1991 FEE system.) + Weierstrass: c = 0. (The basic IEEE form.) + Atkin3: c = a = 0. + Atkin4: c = b = 0. + + However, the code is set up to work with any a, b, c. + The underlying fields F_{p^m} are of odd characteristic, + with all operations are (mod p). The special FEE-class + primes p are of the form: + + p = 2^q - k = 3 (mod 4) + + where k is single-precision. For such p, the mod operations + are especially fast (asymptotically vanishing time with respect + to a multiply). Note that the whole system + works equally well (except for slower execution) for arbitrary + primes p = 3 (mod 4) of the same bit length (q or q+1). + + The elliptic arithmetic now proceeds on the basis of + five fundamental operations that calculate various + numerator/denominator parts of the elliptic terms: + + numer_double(giant x, giant z, giant res, curveParams *par) + res := (x^2 - a z^2)^2 - 4 b (2 x + c z) z^3. + + numer_plus(giant x1, giant x2, giant res, curveParams *par) + res = (x1 x2 + a)(x1 + x2) + 2(c x1 x2 + b). + + denom_double(giant x, giant z, giant res, curveParams *par) + res = 4 z (x^3 + c x^2 z + a x z^2 + b z^3). + + denom_times(giant x1, giant z1, giant x2, giant z2, giant res, + curveParams *par) + res := (x1 z2 - x2 z1)^2 + + numer_times(giant x1, giant z1, giant x2, giant z2, giant res, + curveParams *par) + res := (x1 x2 - a z1 z2)^2 - 4 b(x1 z2 + x2 z1 + c z1 z2) z1 z2 + + If x(+-) represent the sum and difference x-coordinates + respectively, then, in pseudocode, + + For unequal starting coords: + x(+) + x(-) = U = 2 numer_plus/denom_times + x(+) x(-) = V = numer_times/denom_times + + and for equal starting coords: + x(+) = numer_double/denom_double + + The elliptic_add() function uses the fact that + + x(+) = U/2 + s*Sqrt[U^2/4 - V] + + where the sign s = +-1. + +*/ + +#include "ckconfig.h" +#include +#include +#include +#include "platform.h" + +#include "giantIntegers.h" +#include "elliptic.h" +#include "ellipticProj.h" +#include "ckutilities.h" +#include "curveParams.h" +#include "feeDebug.h" +#include "ellipticMeasure.h" +#include "falloc.h" +#include "giantPortCommon.h" + +#if FEE_PROFILE + +unsigned numerDoubleTime; +unsigned numerPlusTime; +unsigned numerTimesTime; +unsigned denomDoubleTime; +unsigned denomTimesTime; +unsigned ellipticTime; +unsigned sigCompTime; +unsigned powerModTime; +unsigned ellAddTime; +unsigned whichCurveTime; +unsigned modgTime; +unsigned mulgTime; +unsigned binvauxTime; +unsigned gsquareTime; + +unsigned numMulg; +unsigned numFeemod; +unsigned numGsquare; +unsigned numBorrows; + +void clearProfile() +{ + numerDoubleTime = 0; + numerPlusTime = 0; + numerTimesTime = 0; + denomDoubleTime = 0; + denomTimesTime = 0; + ellipticTime = 0; + sigCompTime = 0; + powerModTime = 0; + ellAddTime = 0; + whichCurveTime = 0; + modgTime = 0; + mulgTime = 0; + binvauxTime = 0; + gsquareTime = 0; + numMulg = 0; + numFeemod = 0; + numGsquare = 0; + numBorrows = 0; +} + +#endif // FEE_PROFILE + +#if ELL_PROFILE +unsigned ellOddTime; +unsigned ellEvenTime; +unsigned numEllOdds; +unsigned numEllEvens; + +void clearEllProfile() +{ + ellOddTime = 0; + ellEvenTime = 0; + numEllOdds = 0; + numEllEvens = 0; +} + +#endif /* ELL_PROFILE */ +#if ELLIPTIC_MEASURE + +int doEllMeasure; // gather stats on/off */ +int bitsInN; +int numFeeMods; +int numMulgs; + +void dumpEll() +{ + printf("\nbitlen(n) : %d\n", bitsInN); + printf("feemods : %d\n", numFeeMods); + printf("mulgs : %d\n", numMulgs); +} + +#endif // ELLIPTIC_MEASURE + +/********** Globals ********************************/ + +static void make_base(curveParams *par, giant result); // result = with 2^q-k +static int keys_inconsistent(key pub1, key pub2); +/* Return non-zero if pub1, pub2 have inconsistent parameters. + */ + + +static void ell_even(giant x1, giant z1, giant x2, giant z2, curveParams *par); +static void ell_odd(giant x1, giant z1, giant x2, giant z2, giant xxor, + giant zor, curveParams *par); +static void numer_double(giant x, giant z, giant res, curveParams *par); +static void numer_plus(giant x1, giant x2, giant res, curveParams *par); +static void denom_double(giant x, giant z, giant res, curveParams *par); +static void denom_times(giant x1, giant z1, giant x2, giant z2, giant res, + curveParams *par); +static void numer_times(giant x1, giant z1, giant x2, giant z2, giant res, + curveParams *par); +static void feepowermodg(curveParams *par, giant x, giant n); +static void curveOrderJustifyWithRecip(giant g, giant curveOrder, giant recip); + +/* + * Completely rewritten in CryptKit-18, 13 Jan 1997, for new IEEE-style + * curveParameters. + */ +int which_curve(giant x, curveParams *par) + /* Returns (+-1) depending on whether x is on curve + (+-)y^2 = x^3 + c x^2 + a x + b. + */ +{ + giant t1; + giant t2; + giant t3; + int result; + PROF_START; + + t1 = borrowGiant(par->maxDigits); + t2 = borrowGiant(par->maxDigits); + t3 = borrowGiant(par->maxDigits); + + /* First, set t2:= x^3 + c x^2 + a x + b. */ + gtog(x, t2); addg(par->c, t2); + mulg(x, t2); addg(par->a, t2); /* t2 := x^2 + c x + a. */ + feemod(par, t2); + mulg(x, t2); addg(par->b, t2); + feemod(par, t2); + /* Next, test whether t2 is a square. */ + gtog(t2, t1); + make_base(par, t3); iaddg(1, t3); gshiftright(1, t3); /* t3 = (p+1)/2. */ + feepowermodg(par, t1, t3); /* t1 := t2^((p+1)/2) (mod p). */ + if(gcompg(t1, t2) == 0) + result = CURVE_PLUS; + else + result = CURVE_MINUS; + returnGiant(t1); + returnGiant(t2); + returnGiant(t3); + PROF_END(whichCurveTime); + return result; +} + +key new_public(curveParams *cp, int twist) { + key k; + + k = (key) fmalloc(sizeof(keystruct)); + k->cp = cp; + k->twist = twist; + + k->x = newGiant(cp->maxDigits); + if((twist == CURVE_PLUS) && (cp->curveType == FCT_Weierstrass)) { + k->y = newGiant(cp->maxDigits); + } + else { + /* + * no projective algebra. We could optimize and save a few bytes + * here by setting y to NULL, but that really complicates things + * in may other places. Best to have a real giant. + */ + k->y = newGiant(1); + } + return(k); +} + +key new_public_with_key(key old_key, curveParams *cp) +{ + key result; + + result = new_public(cp, old_key->twist); + CKASSERT((old_key->x != NULL) && (old_key->y != NULL)); + CKASSERT((result->x != NULL) && (result->y != NULL)); + gtog(old_key->x, result->x); + gtog(old_key->y, result->y); + return result; +} + +void free_key(key pub) { + if(!pub) { + return; + } + if (pub->x) { + freeGiant(pub->x); + } + if (pub->y) { + freeGiant(pub->y); + } + ffree(pub); +} + +/* + * Specify private data for key created by new_public(). + * Generates k->x. + */ +void set_priv_key_giant(key k, giant privGiant) +{ + curveParams *cp = k->cp; + + /* elliptiy multiply of initial public point times private key */ + #if CRYPTKIT_ELL_PROJ_ENABLE + if((k->twist == CURVE_PLUS) && (cp->curveType == FCT_Weierstrass)) { + /* projective */ + + pointProj pt1 = newPointProj(cp->maxDigits); + + CKASSERT((cp->y1Plus != NULL) && (!isZero(cp->y1Plus))); + CKASSERT(k->y != NULL); + + /* pt1 := {x1Plus, y1Plus, 1} */ + gtog(cp->x1Plus, pt1->x); + gtog(cp->y1Plus, pt1->y); + int_to_giant(1, pt1->z); + + /* pt1 := pt1 * privateKey */ + ellMulProjSimple(pt1, privGiant, cp); + + /* result back to {k->x, k->y} */ + gtog(pt1->x, k->x); + gtog(pt1->y, k->y); + freePointProj(pt1); // FIXME - clear the giants + } + else { + #else + { + #endif /* CRYPTKIT_ELL_PROJ_ENABLE */ + /* FEE */ + if(k->twist == CURVE_PLUS) { + gtog(cp->x1Plus, k->x); + } + else { + gtog(cp->x1Minus, k->x); + } + elliptic_simple(k->x, privGiant, k->cp); + } +} + +int key_equal(key one, key two) { + if (keys_inconsistent(one, two)) return 0; + return !gcompg(one->x, two->x); +} + +static void make_base(curveParams *par, giant result) +/* Jams result with 2^q-k. */ +{ + gtog(par->basePrime, result); +} + +void make_base_prim(curveParams *cp) +/* Jams cp->basePrime with 2^q-k. Assumes valid maxDigits, q, k. */ +{ + giant tmp = borrowGiant(cp->maxDigits); + + CKASSERT(cp->primeType != FPT_General); + int_to_giant(1, cp->basePrime); + gshiftleft((int)cp->q, cp->basePrime); + int_to_giant(cp->k, tmp); + subg(tmp, cp->basePrime); + returnGiant(tmp); +} + +static int sequalg(int n, giant g) { + if((g->sign == 1) && (g->n[0] == n)) return(1); + return(0); +} + + +/* + * Elliptic multiply: x := n * {x, 1} + */ +void elliptic_simple(giant x, giant n, curveParams *par) { + giant ztmp = borrowGiant(par->maxDigits); + giant cur_n = borrowGiant(par->maxDigits); + + START_ELL_MEASURE(n); + int_to_giant(1, ztmp); + elliptic(x, ztmp, n, par); + binvg_cp(par, ztmp); + mulg(ztmp, x); + feemod(par, x); + END_ELL_MEASURE; + + returnGiant(cur_n); + returnGiant(ztmp); +} + +/* + * General elliptic multiply. + * + * {xx, zz} := k * {xx, zz} + */ +void elliptic(giant xx, giant zz, giant k, curveParams *par) { + int len = bitlen(k); + int pos = len - 2; + giant xs; + giant zs; + giant xorg; + giant zorg; + + PROF_START; + if(sequalg(1,k)) return; + if(sequalg(2,k)) { + ell_even(xx, zz, xx, zz, par); + goto out; + } + zs = borrowGiant(par->maxDigits); + xs = borrowGiant(par->maxDigits); + zorg = borrowGiant(par->maxDigits); + xorg = borrowGiant(par->maxDigits); + gtog(xx, xorg); gtog(zz, zorg); + ell_even(xx, zz, xs, zs, par); + do { + if(bitval(k, pos--)) { + ell_odd(xs, zs, xx, zz, xorg, zorg, par); + ell_even(xs, zs, xs, zs, par); + } else { + ell_odd(xx, zz, xs, zs, xorg, zorg, par); + ell_even(xx, zz, xx, zz, par); + } + } while (pos >= 0); // REC fix 9/23/94 + returnGiant(xs); + returnGiant(zs); + returnGiant(xorg); + returnGiant(zorg); +out: + PROF_END(ellipticTime); +} + +/* + * Completely rewritten in CryptKit-18, 13 Jan 1997, for new IEEE-style + * curveParameters. + */ +void elliptic_add(giant x1, giant x2, giant x3, curveParams *par, int s) { + + /* Addition algorithm for x3 = x1 + x2 on the curve, with sign ambiguity s. + From theory, we know that if {x1,1} and {x2,1} are on a curve, then + their elliptic sum (x1,1} + {x2,1} = {x3,1} must have x3 as one of two + values: + + x3 = U/2 + s*Sqrt[U^2/4 - V] + + where sign s = +-1, and U,V are functions of x1,x2. Tho present function + is called a maximum of twice, to settle which of +- is s. When a call + is made, it is guaranteed already that x1, x2 both lie on the same curve + (+- curve); i.e., which curve (+-) is not connected at all with sign s of + the x3 relation. + */ + + giant cur_n; + giant t1; + giant t2; + giant t3; + giant t4; + giant t5; + + PROF_START; + cur_n = borrowGiant(par->maxDigits); + t1 = borrowGiant(par->maxDigits); + t2 = borrowGiant(par->maxDigits); + t3 = borrowGiant(par->maxDigits); + t4 = borrowGiant(par->maxDigits); + t5 = borrowGiant(par->maxDigits); + + if(gcompg(x1, x2)==0) { + int_to_giant(1, t1); + numer_double(x1, t1, x3, par); + denom_double(x1, t1, t2, par); + binvg_cp(par, t2); + mulg(t2, x3); feemod(par, x3); + goto out; + } + numer_plus(x1, x2, t1, par); + int_to_giant(1, t3); + numer_times(x1, t3, x2, t3, t2, par); + int_to_giant(1, t4); int_to_giant(1, t5); + denom_times(x1, t4, x2, t5, t3, par); + binvg_cp(par, t3); + mulg(t3, t1); feemod(par, t1); /* t1 := U/2. */ + mulg(t3, t2); feemod(par, t2); /* t2 := V. */ + /* Now x3 will be t1 +- Sqrt[t1^2 - t2]. */ + gtog(t1, t4); gsquare(t4); feemod(par, t4); + subg(t2, t4); + make_base(par, cur_n); iaddg(1, cur_n); gshiftright(2, cur_n); + /* cur_n := (p+1)/4. */ + feepowermodg(par, t4, cur_n); /* t4 := t2^((p+1)/4) (mod p). */ + gtog(t1, x3); + if(s != SIGN_PLUS) negg(t4); + addg(t4, x3); + feemod(par, x3); + +out: + returnGiant(cur_n); + returnGiant(t1); + returnGiant(t2); + returnGiant(t3); + returnGiant(t4); + returnGiant(t5); + + PROF_END(ellAddTime); +} + +/* + * Key exchange atom. + */ +giant make_pad(giant privGiant, key publicKey) { + curveParams *par = publicKey->cp; + giant result = newGiant(par->maxDigits); + + gtog(publicKey->x, result); + elliptic_simple(result, privGiant, par); + return result; +} + +static void ell_even(giant x1, giant z1, giant x2, giant z2, curveParams *par) { + giant t1; + giant t2; + giant t3; + + EPROF_START; + + t1 = borrowGiant(par->maxDigits); + t2 = borrowGiant(par->maxDigits); + t3 = borrowGiant(par->maxDigits); + + if(par->curveType == FCT_Montgomery) { + /* Begin Montgomery OPT: 10 Jan 98 REC. */ + gtog(x1, t1); gsquare(t1); feemod(par, t1); /* t1 := x1^2. */ + gtog(z1, t2); gsquare(t2); feemod(par, t2); /* t2 := z1^2. */ + + gtog(x1, t3); mulg(z1, t3); feemod(par, t3); + gtog(t3, z2); mulg(par->c, z2); feemod(par, z2); + addg(t1, z2); addg(t2, z2); mulg(t3, z2); gshiftleft(2, z2); + feemod(par, z2); /* z2 := 4 x1 z1 (x1^2 + c x1 z1 + z1^2). */ + gtog(t1, x2); subg(t2, x2); gsquare(x2); feemod(par, x2); + /* x2 := (x1^2 - z1^2)^2. */ + /* End OPT: 10 Jan 98 REC. */ + } + else if((par->curveType == FCT_Weierstrass) && isZero(par->a)) { + /* Begin Atkin3 OPT: 9 Jan 98 REC. */ + gtog(x1, t1); + gsquare(t1); feemod(par, t1); + mulg(x1, t1); feemod(par, t1); /* t1 := x^3. */ + gtog(z1, t2); + gsquare(t2); feemod(par, t2); + mulg(z1, t2); feemod(par, t2); /* t2 := z1^3 */ + mulg(par->b, t2); feemod(par, t2); /* t2 := b z1^3. */ + gtog(t1, t3); addg(t2, t3); /* t3 := x^3 + b z1^3 */ + mulg(z1, t3); feemod(par, t3); /* t3 *= z1 + * = z1 ( x^3 + b z1^3 ) */ + gshiftleft(2, t3); feemod(par, t3); /* t3 = 4 z1 (x1^3 + b z1^3) */ + + gshiftleft(3, t2); /* t2 = 8 b z1^3 */ + subg(t2, t1); /* t1 = x^3 - 8 b z1^3 */ + mulg(x1, t1); feemod(par, t1); /* t1 = x1 (x1^3 - 8 b z1^3) */ + + gtog(t3, z2); + gtog(t1, x2); + /* End OPT: 9 Jan 98 REC. */ + } + else { + numer_double(x1, z1, t1, par); + denom_double(x1, z1, t2, par); + gtog(t1, x2); gtog(t2, z2); + } + returnGiant(t1); + returnGiant(t2); + returnGiant(t3); + + EPROF_END(ellEvenTime); + EPROF_INCR(numEllEvens); + + /* + printf("ell_even end\n"); + printf(" x1 : "); printGiant(x1); + printf(" z1 : "); printGiant(z1); + printf(" x2 : "); printGiant(x2); + printf(" z2 : "); printGiant(z2); + */ +} + +static void ell_odd(giant x1, giant z1, giant x2, giant z2, giant xxor, + giant zor, curveParams *par) +{ + + giant t1; + giant t2; + + EPROF_START; + t1 = borrowGiant(par->maxDigits); + t2 = borrowGiant(par->maxDigits); + + if(par->curveType == FCT_Montgomery) { + /* Begin Montgomery OPT: 10 Jan 98 REC. */ + giant t3 = borrowGiant(par->maxDigits); + giant t4 = borrowGiant(par->maxDigits); + + gtog(x1, t1); addg(z1, t1); /* t1 := x1 + z1. */ + gtog(x2, t2); subg(z2, t2); /* t2 := x2 - z2. */ + gtog(x1, t3); subg(z1, t3); /* t3 := x1 - z1. */ + gtog(x2, t4); addg(z2, t4); /* t4 := x2 + z2. */ + mulg(t2, t1); feemod(par, t1); /* t1 := (x1 + z1)(x2 - z2) */ + mulg(t4, t3); feemod(par, t3); /* t4 := (x2 + z2)(x1 - z1) */ + gtog(t1, z2); subg(t3, z2); /*???gshiftright(1, z2);*/ + /* z2 := ((x1 + z1)(x2 - z2) - x2)/2 */ + gsquare(z2); feemod(par, z2); + mulg(xxor, z2); feemod(par, z2); + gtog(t1, x2); addg(t3, x2); /*????gshiftright(1, x2);*/ + gsquare(x2); feemod(par, x2); + mulg(zor, x2); feemod(par, x2); + + returnGiant(t3); + returnGiant(t4); + } + else if((par->curveType == FCT_Weierstrass) && isZero(par->a)) { + /* Begin Atkin3 OPT: 9 Jan 98 REC. */ + + giant t3 = borrowGiant(par->maxDigits); + giant t4 = borrowGiant(par->maxDigits); + + gtog(x1, t1); mulg(x2, t1); feemod(par, t1); /* t1 := x1 x2. */ + gtog(z1, t2); mulg(z2, t2); feemod(par, t2); /* t2 := z1 z2. */ + gtog(x1, t3); mulg(z2, t3); feemod(par, t3); /* t3 := x1 z2. */ + gtog(z1, t4); mulg(x2, t4); feemod(par, t4); /* t4 := x2 z1. */ + gtog(t3, z2); subg(t4, z2); gsquare(z2); feemod(par, z2); + mulg(xxor, z2); feemod(par, z2); + gtog(t1, x2); gsquare(x2); feemod(par, x2); + addg(t4, t3); mulg(t2, t3); feemod(par, t3); + mulg(par->b, t3); feemod(par, t3); + addg(t3, t3); addg(t3, t3); + subg(t3, x2); mulg(zor, x2); feemod(par, x2); + + returnGiant(t3); + returnGiant(t4); + /* End OPT: 9 Jan 98 REC. */ + } + else { + numer_times(x1, z1, x2, z2, t1, par); + mulg(zor, t1); feemod(par, t1); + denom_times(x1, z1, x2, z2, t2, par); + mulg(xxor, t2); feemod(par, t2); + + gtog(t1, x2); gtog(t2, z2); + } + + returnGiant(t1); + returnGiant(t2); + + EPROF_END(ellOddTime); + EPROF_INCR(numEllOdds); + + /* + printf("ell_odd end\n"); + printf(" x2 : "); printGiant(x2); + printf(" z2 : "); printGiant(z2); + */ +} + +/* + * return codes from keys_inconsistent() and signature_compare(). The actual + * values are not public; they are defined here for debugging. + */ +#define CURVE_PARAM_MISMATCH 1 +#define TWIST_PARAM_MISMATCH 2 +#define SIGNATURE_INVALID 3 + + +/* + * Determine whether two keystructs have compatible parameters (i.e., same + * twist and curveParams). Return 0 if compatible, else non-zero. + */ +static int keys_inconsistent(key pub1, key pub2){ + if(!curveParamsEquivalent(pub1->cp, pub2->cp)) { + return CURVE_PARAM_MISMATCH; + } + if(pub1->twist != pub2->twist) { + return TWIST_PARAM_MISMATCH; + } + return 0; +} + +int signature_compare(giant p0x, giant p1x, giant p2x, curveParams *par) +/* Returns non-zero iff p0x cannot be the x-coordinate of the sum of two points whose respective x-coordinates are p1x, p2x. */ +{ + int ret = 0; + giant t1; + giant t2; + giant t3; + giant t4; + giant t5; + + PROF_START; + + t1 = borrowGiant(par->maxDigits); + t2 = borrowGiant(par->maxDigits); + t3 = borrowGiant(par->maxDigits); + t4 = borrowGiant(par->maxDigits); + t5 = borrowGiant(par->maxDigits); + + if(gcompg(p1x, p2x) == 0) { + int_to_giant(1, t1); + numer_double(p1x, t1, t2, par); + denom_double(p1x, t1, t3, par); + mulg(p0x, t3); subg(t3, t2); + feemod(par, t2); + } else { + numer_plus(p1x, p2x, t1, par); + gshiftleft(1, t1); feemod(par, t1); + int_to_giant(1, t3); + numer_times(p1x, t3, p2x, t3, t2, par); + int_to_giant(1, t4); int_to_giant(1, t5); + denom_times(p1x, t4 , p2x, t5, t3, par); + /* Now we require t3 x0^2 - t1 x0 + t2 == 0. */ + mulg(p0x, t3); feemod(par, t3); + subg(t1, t3); mulg(p0x, t3); + feemod(par, t3); + addg(t3, t2); + feemod(par, t2); + } + + if(!isZero(t2)) ret = SIGNATURE_INVALID; + returnGiant(t1); + returnGiant(t2); + returnGiant(t3); + returnGiant(t4); + returnGiant(t5); + PROF_END(sigCompTime); + return(ret); +} + + +static void numer_double(giant x, giant z, giant res, curveParams *par) +/* Numerator algebra. + res := (x^2 - a z^2)^2 - 4 b (2 x + c z) z^3. + */ +{ + giant t1; + giant t2; + + PROF_START; + t1 = borrowGiant(par->maxDigits); + t2 = borrowGiant(par->maxDigits); + + gtog(x, t1); gsquare(t1); feemod(par, t1); + gtog(z, res); gsquare(res); feemod(par, res); + gtog(res, t2); + if(!isZero(par->a) ) { + if(!isone(par->a)) { /* Speedup - REC 17 Jan 1997. */ + mulg(par->a, res); feemod(par, res); + } + subg(res, t1); feemod(par, t1); + } + gsquare(t1); feemod(par, t1); + /* t1 := (x^2 - a z^2)^2. */ + if(isZero(par->b)) { /* Speedup - REC 17 Jan 1997. */ + gtog(t1, res); + goto done; + } + if(par->curveType != FCT_Weierstrass) { // i.e., !isZero(par->c) + // Speedup - REC 17 Jan 1997. + gtog(z, res); mulg(par->c, res); feemod(par, res); + } else { + int_to_giant(0, res); + } + addg(x, res); addg(x, res); mulg(par->b, res); + feemod(par, res); + gshiftleft(2, res); mulg(z, res); feemod(par, res); + mulg(t2, res); feemod(par, res); + negg(res); addg(t1, res); + feemod(par, res); + +done: + returnGiant(t1); + returnGiant(t2); + PROF_END(numerDoubleTime); +} + +static void numer_plus(giant x1, giant x2, giant res, curveParams *par) +/* Numerator algebra. + res = (x1 x2 + a)(x1 + x2) + 2(c x1 x2 + b). + */ +{ + giant t1; + giant t2; + + PROF_START; + t1 = borrowGiant(par->maxDigits); + t2 = borrowGiant(par->maxDigits); + + gtog(x1, t1); mulg(x2, t1); feemod(par, t1); + gtog(x2, t2); addg(x1, t2); feemod(par, t2); + gtog(t1, res); + if(!isZero(par->a)) + addg(par->a, res); + mulg(t2, res); feemod(par, res); + if(par->curveType == FCT_Weierstrass) { // i.e., isZero(par->c) + int_to_giant(0, t1); + } + else { + mulg(par->c, t1); feemod(par, t1); + } + if(!isZero(par->b)) + addg(par->b, t1); + gshiftleft(1, t1); + addg(t1, res); feemod(par, res); + + returnGiant(t1); + returnGiant(t2); + PROF_END(numerPlusTime); +} + +static void denom_double(giant x, giant z, giant res, curveParams *par) +/* Denominator algebra. + res = 4 z (x^3 + c x^2 z + a x z^2 + b z^3). */ +{ + giant t1; + giant t2; + + PROF_START; + t1 = borrowGiant(par->maxDigits); + t2 = borrowGiant(par->maxDigits); + + gtog(x, res); gtog(z, t1); + if(par->curveType != FCT_Weierstrass) { // i.e., !isZero(par->c) + gtog(par->c, t2); mulg(t1, t2); feemod(par, t2); + addg(t2, res); + } + mulg(x, res); feemod(par, res); + gsquare(t1); feemod(par, t1); + if(!isZero(par->a)) { + gtog(t1, t2); + mulg(par->a, t2); feemod(par, t2); + addg(t2, res); + } + mulg(x, res); feemod(par, res); + if(!isZero(par->b)) { + mulg(z, t1); feemod(par, t1); + mulg(par->b, t1); feemod(par, t1); + addg(t1, res); + } + mulg(z, res); gshiftleft(2, res); + feemod(par, res); + + returnGiant(t1); + returnGiant(t2); + PROF_END(denomDoubleTime); +} + + + +static void denom_times(giant x1, giant z1, giant x2, giant z2, giant res, + curveParams *par) +/* Denominator algebra. + res := (x1 z2 - x2 z1)^2 + */ +{ + giant t1; + + PROF_START; + t1 = borrowGiant(par->maxDigits); + + gtog(x1, res); mulg(z2, res); feemod(par, res); + gtog(z1, t1); mulg(x2, t1); feemod(par, t1); + subg(t1, res); gsquare(res); feemod(par, res); + + returnGiant(t1); + PROF_END(denomTimesTime); +} + +static void numer_times(giant x1, giant z1, giant x2, giant z2, giant res, + curveParams *par) +/* Numerator algebra. + res := (x1 x2 - a z1 z2)^2 - + 4 b(x1 z2 + x2 z1 + c z1 z2) z1 z2 + */ +{ + giant t1; + giant t2; + giant t3; + giant t4; + + PROF_START; + t1 = borrowGiant(par->maxDigits); + t2 = borrowGiant(par->maxDigits); + t3 = borrowGiant(par->maxDigits); + t4 = borrowGiant(par->maxDigits); + + gtog(x1, t1); mulg(x2, t1); feemod(par, t1); + gtog(z1, t2); mulg(z2, t2); feemod(par, t2); + gtog(t1, res); + if(!isZero(par->a)) { + gtog(par->a, t3); + mulg(t2, t3); feemod(par, t3); + subg(t3, res); + } + gsquare(res); feemod(par, res); + if(isZero(par->b)) + goto done; + if(par->curveType != FCT_Weierstrass) { // i.e., !isZero(par->c) + gtog(par->c, t3); + mulg(t2, t3); feemod(par, t3); + } else int_to_giant(0, t3); + gtog(z1, t4); mulg(x2, t4); feemod(par, t4); + addg(t4, t3); + gtog(x1, t4); mulg(z2, t4); feemod(par, t4); + addg(t4, t3); mulg(par->b, t3); feemod(par, t3); + mulg(t2, t3); gshiftleft(2, t3); feemod(par, t3); + subg(t3, res); + feemod(par, res); + +done: + returnGiant(t1); + returnGiant(t2); + returnGiant(t3); + returnGiant(t4); + PROF_END(numerTimesTime); +} + +/* + * New, 13 Jan 1997. + */ +static void feepowermodg(curveParams *par, giant x, giant n) +/* Power ladder. + x := x^n (mod 2^q-k) + */ +{ + int len, pos; + giant t1; + + PROF_START; + t1 = borrowGiant(par->maxDigits); + gtog(x, t1); + int_to_giant(1, x); + len = bitlen(n); + pos = 0; + while(1) { + if(bitval(n, pos++)) { + mulg(t1, x); + feemod(par, x); + } + if(pos>=len) break; + gsquare(t1); + feemod(par, t1); + } + returnGiant(t1); + PROF_END(powerModTime); +} + +/* + * Set g := g mod curveOrder; + * force g to be between 2 and (curveOrder-2), inclusive. + * + * Tolerates zero curve orders (indicating "not known"). + */ + +/* + * This version is not normally used; it's for when the reciprocal of + * curveOrder is not known and won't be used again. + */ +void curveOrderJustify(giant g, giant curveOrder) +{ + giant recip; + + CKASSERT(!isZero(curveOrder)); + + recip = borrowGiant(2 * abs(g->sign)); + make_recip(curveOrder, recip); + curveOrderJustifyWithRecip(g, curveOrder, recip); + returnGiant(recip); +} + +/* + * New optimzation of curveOrderJustify using known reciprocal, 11 June 1997. + * g is set to be within [2, curveOrder-2]. + */ +static void curveOrderJustifyWithRecip(giant g, giant curveOrder, giant recip) +{ + giant tmp; + + CKASSERT(!isZero(curveOrder)); + + modg_via_recip(curveOrder, recip, g); // g now in [0, curveOrder-1] + + if(isZero(g)) { + /* + * First degenerate case - (g == 0) : set g := 2 + */ + dbgLog(("curveOrderJustify: case 1\n")); + int_to_giant(2, g); + return; + } + if(isone(g)) { + /* + * Second case - (g == 1) : set g := 2 + */ + dbgLog(("curveOrderJustify: case 2\n")); + int_to_giant(2, g); + return; + } + tmp = borrowGiant(g->capacity); + gtog(g, tmp); + iaddg(1, tmp); + if(gcompg(tmp, curveOrder) == 0) { + /* + * Third degenerate case - (g == (curveOrder-1)) : set g -= 1 + */ + dbgLog(("curveOrderJustify: case 3\n")); + int_to_giant(1, tmp); + subg(tmp, g); + } + returnGiant(tmp); + return; +} + +#define RECIP_DEBUG 0 +#if RECIP_DEBUG +#define recipLog(x) printf x +#else // RECIP_DEBUG +#define recipLog(x) +#endif // RECIP_DEBUG + +/* + * curveOrderJustify routines with specific orders, using (and possibly + * generating) associated reciprocals. + */ +void lesserX1OrderJustify(giant g, curveParams *cp) +{ + /* + * Note this is a pointer to a giant in *cp, not a newly + * malloc'd giant! + */ + giant lesserX1Ord = lesserX1Order(cp); + + if(isZero(lesserX1Ord)) { + return; + } + + /* + * Calculate reciprocal if we don't have it + */ + if(isZero(cp->lesserX1OrderRecip)) { + if((lesserX1Ord == cp->x1OrderPlus) && + (!isZero(cp->x1OrderPlusRecip))) { + /* + * lesserX1Ord happens to be x1OrderPlus, and we + * have a reciprocal for x1OrderPlus. Copy it over. + */ + recipLog(( + "x1OrderPlusRecip --> lesserX1OrderRecip\n")); + gtog(cp->x1OrderPlusRecip, cp->lesserX1OrderRecip); + } + else { + /* Calculate the reciprocal. */ + recipLog(("calc lesserX1OrderRecip\n")); + make_recip(lesserX1Ord, cp->lesserX1OrderRecip); + } + } + else { + recipLog(("using existing lesserX1OrderRecip\n")); + } + curveOrderJustifyWithRecip(g, lesserX1Ord, cp->lesserX1OrderRecip); +} + +/* + * Common code used by x1OrderPlusJustify() and x1OrderPlusMod() to generate + * reciprocal of x1OrderPlus. + * 8 Sep 1998 - also used by feeSigSign(). + */ +void calcX1OrderPlusRecip(curveParams *cp) +{ + if(isZero(cp->x1OrderPlusRecip)) { + if((cp->x1OrderPlus == lesserX1Order(cp)) && + (!isZero(cp->lesserX1OrderRecip))) { + /* + * lesserX1Order happens to be x1OrderPlus, and we + * have a reciprocal for lesserX1Order. Copy it over. + */ + recipLog(( + "lesserX1OrderRecip --> x1OrderPlusRecip\n")); + gtog(cp->lesserX1OrderRecip, cp->x1OrderPlusRecip); + } + else { + /* Calculate the reciprocal. */ + recipLog(("calc x1OrderPlusRecip\n")); + make_recip(cp->x1OrderPlus, cp->x1OrderPlusRecip); + } + } + else { + recipLog(("using existing x1OrderPlusRecip\n")); + } +} + +void x1OrderPlusJustify(giant g, curveParams *cp) +{ + CKASSERT(!isZero(cp->x1OrderPlus)); + + /* + * Calculate reciprocal if we don't have it + */ + calcX1OrderPlusRecip(cp); + curveOrderJustifyWithRecip(g, cp->x1OrderPlus, cp->x1OrderPlusRecip); +} + +/* + * g := g mod x1OrderPlus. Result may be zero. + */ +void x1OrderPlusMod(giant g, curveParams *cp) +{ + CKASSERT(!isZero(cp->x1OrderPlus)); + + /* + * Calculate reciprocal if we don't have it + */ + calcX1OrderPlusRecip(cp); + modg_via_recip(cp->x1OrderPlus, cp->x1OrderPlusRecip, g); +} + +/* + * New general-purpose giant mod routine, 8 Jan 97. + * x := x mod basePrime. + */ + +/* + * This stuff is used to analyze the critical loop behavior inside feemod(). + */ +#define FEEMOD_LOOP_TEST 0 +#if FEEMOD_LOOP_TEST +/* + * these two are only examined via debugger + */ +unsigned feemodCalls = 0; // calls to feemod() +unsigned feemodMultLoops = 0; // times while() loop runs > once +#define FEEMOD_LOOP_INCR feemodLoops++ +#define FEEMOD_CALL_INCR feemodCalls++ +#else // FEEMOD_LOOP_TEST +#define FEEMOD_LOOP_INCR +#define FEEMOD_CALL_INCR +#endif // FEEMOD_LOOP_TEST + + +void +feemod(curveParams *par, giant x) +{ + int sign, sign2; + giant t1; + giant t3; + giant t4; + giant t5; + #if FEEMOD_LOOP_TEST + unsigned feemodLoops = 0; + #endif // FEEMOD_LOOP_TEST + + FEEMOD_CALL_INCR; // for FEEMOD_LOOP_TEST + INCR_FEEMODS; // for ellipticMeasure + PROF_INCR(numFeemod); // for general profiling + + switch(par->primeType) { + case FPT_Mersenne: + /* + * Super-optimized Mersenne prime modulus case + */ + gmersennemod(par->q, x); + break; + + case FPT_FEE: + /* + * General 2**q-k case + */ + sign = (x->sign < 0) ? -1 : 1; + sign2 = 1; + x->sign = abs(x->sign); + if(gcompg(par->basePrime, x) >= 0) { + goto outFee; + } + t1 = borrowGiant(par->maxDigits); + t3 = borrowGiant(par->maxDigits); + t4 = borrowGiant(par->maxDigits); + t5 = borrowGiant(par->maxDigits); + + /* Begin OPT: 11 Jan 98 REC. */ + if( ((par->q & (GIANT_BITS_PER_DIGIT - 1)) == 0) && + (par->k >= 0) && + (par->k < GIANT_DIGIT_MASK)) { + + /* + * Microscopic mod for certain regions of {q,k} + * parameter space. + */ + int j, digits, excess, max; + giantDigit carry; + giantDigit termHi; // double precision + giantDigit termLo; + giantDigit *p1, *p2; + + digits = par->q >> GIANT_LOG2_BITS_PER_DIGIT; + while(bitlen(x) > par->q) { + excess = (x->sign) - digits; + max = (excess > digits) ? excess : digits; + carry = 0; + + p1 = &x->n[0]; + p2 = &x->n[digits]; + + if(excess <= digits) { + carry = VectorMultiply(par->k, + p2, + excess, + p1); + + /* propagate final carry */ + p1 += excess; + for(j = excess; j < digits; j++) { + + /* + * term = *p1 + carry; + * *p1++ = term & 65535; + * carry = term >> 16; + */ + termLo = giantAddDigits(*p1, carry, &carry); + *p1++ = termLo; + } + } else { + carry = VectorMultiply(par->k, + p2, + digits, + p1); + p1 += digits; + p2 += digits; + for(j = digits; j < excess; j++) { + /* + * term = (par->k)*(*p2++) + carry; + */ + giantMulDigits(par->k, + *p2++, + &termLo, + &termHi); + giantAddDouble(&termLo, &termHi, carry); + + /* + * *p1++ = term & 65535; + * carry = term >> 16; + */ + *p1++ = termLo; + carry = termHi; + } + } + + if(carry > 0) { + x->n[max] = carry; + } else { + while(--max){ + if(x->n[max] != 0) break; + } + } + x->sign = max + 1; + FEEMOD_LOOP_INCR; + } + } else { /* Macroscopic mod for general PT_FEE case. */ + int_to_giant(par->k, t4); + while(bitlen(x) > par->q) { + /* Enter fast loop, as in FEE patent. */ + int_to_giant(1, t5); + gtog(x, t3); + extractbits(par->q, t3, x); + while(bitlen(t3) > par->q) { + gshiftright(par->q, t3); + extractbits(par->q, t3, t1); + PAUSE_ELL_MEASURE; + mulg(t4, t5); + mulg(t5, t1); + RESUME_ELL_MEASURE; + addg(t1, x); + } + FEEMOD_LOOP_INCR; + } + } + /* End OPT: 11 Jan 98 REC. */ + + sign2 = (x->sign < 0)? -1: 1; + x->sign = abs(x->sign); + returnGiant(t1); + returnGiant(t3); + returnGiant(t4); + returnGiant(t5); + outFee: + if(gcompg(x, par->basePrime) >=0) subg(par->basePrime, x); + if(sign != sign2) { + giant t2 = borrowGiant(par->maxDigits); + gtog(par->basePrime, t2); + subg(x, t2); + gtog(t2, x); + returnGiant(t2); + } + break; + /* end case PT_FEE */ + + case FPT_General: + /* + * Use brute-force modg. + */ + #if FEE_DEBUG + if(par->basePrimeRecip == NULL) { + CKRaise("feemod(PT_GENERAL): No basePrimeRecip!\n"); + } + #endif /* FEE_DEBUG */ + modg_via_recip(par->basePrime, par->basePrimeRecip, x); + break; + + case FPT_Default: + /* never appears here */ + CKASSERT(0); + break; + } /* switch primeType */ + + #if FEEMOD_LOOP_TEST + if(feemodLoops > 1) { + feemodMultLoops++; + if(feemodLoops > 2) { + printf("feemod while loop executed %d times\n", feemodLoops); + } + } + #endif // FEEMOD_LOOP_TEST + + return; +} + +/* + * For a given curveParams, calculate minBytes and maxDigits. + * Assumes valid primeType, and also a valid basePrime for PT_GENERAL. + */ +void calcGiantSizes(curveParams *cp) +{ + + if(cp->primeType == FPT_General) { + cp->minBytes = (bitlen(cp->basePrime) + 7) / 8; + } + else { + cp->minBytes = giantMinBytes(cp->q, cp->k); + } + cp->maxDigits = giantMaxDigits(cp->minBytes); +} + +unsigned giantMinBytes(unsigned q, int k) +{ + unsigned kIsNeg = (k < 0) ? 1 : 0; + + return (q + 7 + kIsNeg) / 8; +} + +/* + * min value for "extra" bytes. Derived from the fact that during sig verify, + * we have to multiply a giant representing a digest - which may be + * 20 bytes for SHA1 - by a giant of minBytes. + */ +#define MIN_EXTRA_BYTES 20 + +unsigned giantMaxDigits(unsigned minBytes) +{ + unsigned maxBytes = 4 * minBytes; + + if((maxBytes - minBytes) < MIN_EXTRA_BYTES) { + maxBytes = minBytes + MIN_EXTRA_BYTES; + } + return BYTES_TO_GIANT_DIGITS(maxBytes); +} + + +/* + * Optimized binvg(basePrime, x). Avoids the general modg() in favor of + * feemod. + */ +int binvg_cp(curveParams *cp, giant x) +{ + feemod(cp, x); + return(binvaux(cp->basePrime, x)); +} + +/* + * Optimized binvg(x1OrderPlus, x). Uses x1OrderPlusMod(). + */ +int binvg_x1OrderPlus(curveParams *cp, giant x) +{ + x1OrderPlusMod(x, cp); + return(binvaux(cp->x1OrderPlus, x)); +}