Security-59754.41.1.tar.gz

[apple/security.git] / OSX / libsecurity_cryptkit / lib / elliptic.c

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