elliptic.c   [plain text]

/* 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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#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));
}