]> git.saurik.com Git - apple/security.git/blobdiff - Security/libsecurity_cryptkit/lib/elliptic.c
Security-57336.1.9.tar.gz
[apple/security.git] / Security / libsecurity_cryptkit / lib / elliptic.c
diff --git a/Security/libsecurity_cryptkit/lib/elliptic.c b/Security/libsecurity_cryptkit/lib/elliptic.c
deleted file mode 100644 (file)
index 3f45d19..0000000
+++ /dev/null
@@ -1,1437 +0,0 @@
-/* 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));
-}