/ OSX / libsecurity_cryptkit / lib / elliptic.c
elliptic.c
   1  /* Copyright (c) 1998,2011,2014 Apple Inc.  All Rights Reserved.
   2   *
   3   * NOTICE: USE OF THE MATERIALS ACCOMPANYING THIS NOTICE IS SUBJECT
   4   * TO THE TERMS OF THE SIGNED "FAST ELLIPTIC ENCRYPTION (FEE) REFERENCE
   5   * SOURCE CODE EVALUATION AGREEMENT" BETWEEN APPLE, INC. AND THE
   6   * ORIGINAL LICENSEE THAT OBTAINED THESE MATERIALS FROM APPLE,
   7   * INC.  ANY USE OF THESE MATERIALS NOT PERMITTED BY SUCH AGREEMENT WILL
   8   * EXPOSE YOU TO LIABILITY.
   9   ***************************************************************************
  10  
  11     elliptic.c - Library for Apple-proprietary Fast Elliptic
  12     Encryption. The algebra in this module ignores elliptic point's
  13     y-coordinates.
  14  
  15     Patent information:
  16  
  17     FEE is patented, U.S. Patents #5159632 (1992), #5271061 (1993),
  18  	#5463690 (1994).  These patents are implemented
  19          in various elliptic algebra functions such as
  20          numer/denom_plus/times(), and in the fact of special
  21          forms for primes: p = 2^q-k.
  22  
  23     Digital signature using fast elliptic addition, in
  24          U. S. Patent #5581616 (1996), is implemented in the
  25          signature_compare() function.
  26  
  27     FEED (Fast Elliptic Encryption) is patent pending (as of Jan 1998).
  28     	Various functions such as elliptic_add() implement the patent ideas.
  29  
  30  
  31     Modification history since the U.S. Patent:
  32     -------------------------------------------
  33  	10/06/98		ap
  34    		Changed to compile with C++.
  35      9 Sep 98 at Apple
  36      	cp->curveType optimizations.
  37  	Removed code which handled "unknown" curve orders.
  38  	elliptic() now exported for timing measurements.
  39     21 Apr 98 at Apple
  40     	Used inline platform-dependent giant arithmetic.
  41     20 Jan 98 at Apple
  42     	feemod now handles PT_MERSENNE, PT_FEE, PT_GENERAL.
  43  	Added calcGiantSizes(), rewrote giantMinBytes(), giantMaxShorts().
  44  	Updated heading comments on FEE curve algebra.
  45     11 Jan 98 at Apple
  46     	Microscopic feemod optimization.
  47     10 Jan 98 at Apple
  48   	ell_odd, ell_even() Montgomery optimization.
  49     09 Jan 98 at Apple
  50   	ell_odd, ell_even() Atkin3 optimization.
  51     08 Jan 97 at Apple
  52     	Cleaned up some debugging code; added gsquareTime
  53     11 Jun 97 at Apple
  54   	Mods for modg_via_recip(), divg_via_recip() math
  55  	Deleted a lot of obsolete code (previously ifdef'd out)
  56  	Added lesserX1OrderJustify(), x1OrderPlusJustify()
  57  	Added binvg_cp(), avoiding general modg in favor of feemod
  58     05 Feb 97 at Apple
  59     	New optimized numer_plus(), denom_double(), and numer_times()
  60  	All calls to borrowGiant() and newGiant have explicit giant size
  61     08 Jan 97 at NeXT
  62     	Major mods to accomodate IEEE-style curve parameters.
  63  	New functions feepowermodg() and curveOrderJustify();
  64  	elliptic_add(), elliptic(), signature_compare(), and
  65  	which_curve() completely rewritten.
  66     19 Dec 96 at NeXT
  67     	Added mersennePrimes[24..26]
  68     08 Aug 96 at NeXT
  69     	Fixed giant leak in elliptic_add()
  70     05 Aug 96 at NeXT
  71     	Removed dead code
  72     24 Jul 96 at NeXT
  73   	Added ENGINE_127_BITS dependency for use of security engine
  74     24 Oct 92 at NeXT
  75     	Modified new_public_from_private()
  76  	Created.
  77  
  78  
  79     FEE curve algebra, Jan 1997.
  80  
  81     Curves are:
  82  
  83        y^2 = x^3 + c x^2 + a x + b
  84  
  85     where useful parameterizations for practical purposes are:
  86  
  87        Montgomery: a = 1, b = 0. (The original 1991 FEE system.)
  88        Weierstrass: c = 0.  (The basic IEEE form.)
  89        Atkin3: c = a = 0.
  90        Atkin4: c = b = 0.
  91  
  92     However, the code is set up to work with any a, b, c.
  93     The underlying fields F_{p^m} are of odd characteristic,
  94     with all operations are (mod p).  The special FEE-class
  95     primes p are of the form:
  96  
  97        p = 2^q - k = 3 (mod 4)
  98  
  99     where k is single-precision.  For such p, the mod operations
 100     are especially fast (asymptotically vanishing time with respect
 101     to a multiply).  Note that the whole system
 102     works equally well (except for slower execution) for arbitrary
 103     primes p = 3 (mod 4) of the same bit length (q or q+1).
 104  
 105     The elliptic arithmetic now proceeds on the basis of
 106     five fundamental operations that calculate various
 107     numerator/denominator parts of the elliptic terms:
 108  
 109     numer_double(giant x, giant z, giant res, curveParams *par)
 110     res := (x^2 - a z^2)^2 - 4 b (2 x + c z) z^3.
 111  
 112     numer_plus(giant x1, giant x2, giant res, curveParams *par)
 113     res = (x1 x2 + a)(x1 + x2) + 2(c x1 x2 + b).
 114  
 115     denom_double(giant x, giant z, giant res, curveParams *par)
 116     res = 4 z (x^3 + c x^2 z + a x z^2 + b z^3).
 117  
 118     denom_times(giant x1, giant z1, giant x2, giant z2, giant res,
 119       curveParams *par)
 120     res := (x1 z2 - x2 z1)^2
 121  
 122     numer_times(giant x1, giant z1, giant x2, giant z2, giant res,
 123        curveParams *par)
 124     res := (x1 x2 - a z1 z2)^2 - 4 b(x1 z2 + x2 z1 + c z1 z2) z1 z2
 125  
 126     If x(+-) represent the sum and difference x-coordinates
 127     respectively, then, in pseudocode,
 128  
 129     For unequal starting coords:
 130      x(+) + x(-) = U = 2 numer_plus/denom_times
 131       x(+) x(-)  = V = numer_times/denom_times
 132  
 133     and for equal starting coords:
 134       x(+) = numer_double/denom_double
 135  
 136     The elliptic_add() function uses the fact that
 137  
 138       x(+) = U/2 + s*Sqrt[U^2/4 - V]
 139  
 140     where the sign s = +-1.
 141  
 142  */
 143  
 144  #include "ckconfig.h"
 145  #include <stdio.h>
 146  #include <stdlib.h>
 147  #include <string.h>
 148  #include "platform.h"
 149  
 150  #include "giantIntegers.h"
 151  #include "elliptic.h"
 152  #include "ellipticProj.h"
 153  #include "ckutilities.h"
 154  #include "curveParams.h"
 155  #include "feeDebug.h"
 156  #include "ellipticMeasure.h"
 157  #include "falloc.h"
 158  #include "giantPortCommon.h"
 159  
 160  #if	FEE_PROFILE
 161  
 162  unsigned numerDoubleTime;
 163  unsigned numerPlusTime;
 164  unsigned numerTimesTime;
 165  unsigned denomDoubleTime;
 166  unsigned denomTimesTime;
 167  unsigned ellipticTime;
 168  unsigned sigCompTime;
 169  unsigned powerModTime;
 170  unsigned ellAddTime;
 171  unsigned whichCurveTime;
 172  unsigned modgTime;
 173  unsigned mulgTime;
 174  unsigned binvauxTime;
 175  unsigned gsquareTime;
 176  
 177  unsigned numMulg;
 178  unsigned numFeemod;
 179  unsigned numGsquare;
 180  unsigned numBorrows;
 181  
 182  void clearProfile()
 183  {
 184  	numerDoubleTime = 0;
 185  	numerPlusTime = 0;
 186  	numerTimesTime = 0;
 187  	denomDoubleTime = 0;
 188  	denomTimesTime = 0;
 189  	ellipticTime = 0;
 190  	sigCompTime = 0;
 191  	powerModTime = 0;
 192  	ellAddTime = 0;
 193  	whichCurveTime = 0;
 194  	modgTime = 0;
 195  	mulgTime = 0;
 196  	binvauxTime = 0;
 197  	gsquareTime = 0;
 198  	numMulg = 0;
 199  	numFeemod = 0;
 200  	numGsquare = 0;
 201  	numBorrows = 0;
 202  }
 203  
 204  #endif	// FEE_PROFILE
 205  
 206  #if	ELL_PROFILE
 207  unsigned ellOddTime;
 208  unsigned ellEvenTime;
 209  unsigned numEllOdds;
 210  unsigned numEllEvens;
 211  
 212  void clearEllProfile()
 213  {
 214  	ellOddTime = 0;
 215  	ellEvenTime = 0;
 216  	numEllOdds = 0;
 217  	numEllEvens = 0;
 218  }
 219  
 220  #endif	/* ELL_PROFILE */
 221  #if	ELLIPTIC_MEASURE
 222  
 223  int doEllMeasure;	// gather stats on/off */
 224  int bitsInN;
 225  int numFeeMods;
 226  int numMulgs;
 227  
 228  void dumpEll()
 229  {
 230  	printf("\nbitlen(n) : %d\n", bitsInN);
 231  	printf("feemods   : %d\n", numFeeMods);
 232  	printf("mulgs     : %d\n", numMulgs);
 233  }
 234  
 235  #endif	// ELLIPTIC_MEASURE
 236  
 237  /********** Globals ********************************/
 238  
 239  static void make_base(curveParams *par, giant result); // result = with 2^q-k
 240  static int keys_inconsistent(key pub1, key pub2);
 241  /* Return non-zero if pub1, pub2 have inconsistent parameters.
 242   */
 243  
 244  
 245  static void ell_even(giant x1, giant z1, giant x2, giant z2, curveParams *par);
 246  static void ell_odd(giant x1, giant z1, giant x2, giant z2, giant xxor,
 247  	giant zor, curveParams *par);
 248  static void numer_double(giant x, giant z, giant res, curveParams *par);
 249  static void numer_plus(giant x1, giant x2, giant res, curveParams *par);
 250  static void denom_double(giant x, giant z, giant res, curveParams *par);
 251  static void denom_times(giant x1, giant z1, giant x2, giant z2, giant res,
 252  	curveParams *par);
 253  static void numer_times(giant x1, giant z1, giant x2, giant z2, giant res,
 254  	curveParams *par);
 255  static void feepowermodg(curveParams *par, giant x, giant n);
 256  static void curveOrderJustifyWithRecip(giant g, giant curveOrder, giant recip);
 257  
 258  /*
 259   * Completely rewritten in CryptKit-18, 13 Jan 1997, for new IEEE-style
 260   * curveParameters.
 261   */
 262  int which_curve(giant x, curveParams *par)
 263   /* Returns (+-1) depending on whether x is on curve
 264     (+-)y^2 = x^3 + c x^2 + a x + b.
 265   */
 266  {
 267      giant t1;
 268      giant t2;
 269      giant t3;
 270      int result;
 271      PROF_START;
 272  
 273      t1 = borrowGiant(par->maxDigits);
 274      t2 = borrowGiant(par->maxDigits);
 275      t3 = borrowGiant(par->maxDigits);
 276  
 277     /* First, set t2:= x^3 + c x^2 + a x + b. */
 278      gtog(x, t2); addg(par->c, t2);
 279      mulg(x, t2); addg(par->a, t2);  /* t2 := x^2 + c x + a. */
 280      feemod(par, t2);
 281      mulg(x, t2); addg(par->b, t2);
 282      feemod(par, t2);
 283      /* Next, test whether t2 is a square. */
 284      gtog(t2, t1);
 285      make_base(par, t3); iaddg(1, t3); gshiftright(1, t3); /* t3 = (p+1)/2. */
 286      feepowermodg(par, t1, t3); 		      /* t1 := t2^((p+1)/2) (mod p). */
 287      if(gcompg(t1, t2) == 0)
 288              result = CURVE_PLUS;
 289      else
 290              result = CURVE_MINUS;
 291      returnGiant(t1);
 292      returnGiant(t2);
 293      returnGiant(t3);
 294      PROF_END(whichCurveTime);
 295      return result;
 296  }
 297  
 298  key new_public(curveParams *cp, int twist) {
 299      key k;
 300  
 301      k = (key) fmalloc(sizeof(keystruct));
 302      k->cp = cp;
 303      k->twist = twist;
 304  
 305      k->x = newGiant(cp->maxDigits);
 306      if((twist == CURVE_PLUS) && (cp->curveType == FCT_Weierstrass)) {
 307  		k->y = newGiant(cp->maxDigits);
 308      }
 309      else {
 310      	/*
 311  		 * no projective algebra. We could optimize and save a few bytes
 312  		 * here by setting y to NULL, but that really complicates things
 313  		 * in may other places. Best to have a real giant.
 314  		 */
 315  		k->y = newGiant(1);
 316      }
 317      return(k);
 318  }
 319  
 320  key new_public_with_key(key old_key, curveParams *cp)
 321  {
 322  	key result;
 323  
 324  	result = new_public(cp, old_key->twist);
 325  	CKASSERT((old_key->x != NULL) && (old_key->y != NULL));
 326  	CKASSERT((result->x != NULL) && (result->y != NULL));
 327  	gtog(old_key->x, result->x);
 328  	gtog(old_key->y, result->y);
 329  	return result;
 330  }
 331  
 332  void free_key(key pub) {
 333  	if(!pub) {
 334  		return;
 335  	}
 336  	if (pub->x) {
 337  		freeGiant(pub->x);
 338  	}
 339  	if (pub->y) {
 340  		freeGiant(pub->y);
 341  	}
 342  	ffree(pub);
 343  }
 344  
 345  /*
 346   * Specify private data for key created by new_public().
 347   * Generates k->x.
 348   */
 349  void set_priv_key_giant(key k, giant privGiant)
 350  {
 351  	curveParams *cp = k->cp;
 352  
 353  	/* elliptiy multiply of initial public point times private key */
 354  	if((k->twist == CURVE_PLUS) && (cp->curveType == FCT_Weierstrass)) {
 355  		/* projective */
 356  
 357  		pointProj pt1 = newPointProj(cp->maxDigits);
 358  
 359  		CKASSERT((cp->y1Plus != NULL) && (!isZero(cp->y1Plus)));
 360  		CKASSERT(k->y != NULL);
 361  
 362  		/* pt1 := {x1Plus, y1Plus, 1} */
 363  		gtog(cp->x1Plus, pt1->x);
 364  		gtog(cp->y1Plus, pt1->y);
 365  		int_to_giant(1, pt1->z);
 366  
 367  		/* pt1 := pt1 * privateKey */
 368  		ellMulProjSimple(pt1, privGiant, cp);
 369  
 370  		/* result back to {k->x, k->y} */
 371  		gtog(pt1->x, k->x);
 372  		gtog(pt1->y, k->y);
 373  		freePointProj(pt1);	// FIXME - clear the giants
 374  	}
 375  	else {
 376          
 377  		/* FEE */
 378  		if(k->twist == CURVE_PLUS) {
 379  			gtog(cp->x1Plus, k->x);
 380  		}
 381  		else {
 382  			gtog(cp->x1Minus, k->x);
 383  		}
 384  		elliptic_simple(k->x, privGiant, k->cp);
 385  	}
 386  }
 387  
 388  int key_equal(key one, key two) {
 389      if (keys_inconsistent(one, two)) return 0;
 390      return !gcompg(one->x, two->x);
 391  }
 392  
 393  static void make_base(curveParams *par, giant result)
 394  /* Jams result with 2^q-k. */
 395  {
 396      gtog(par->basePrime, result);
 397  }
 398  
 399  void make_base_prim(curveParams *cp)
 400  /* Jams cp->basePrime with 2^q-k. Assumes valid maxDigits, q, k. */
 401  {
 402      giant tmp = borrowGiant(cp->maxDigits);
 403  
 404      CKASSERT(cp->primeType != FPT_General);
 405      int_to_giant(1, cp->basePrime);
 406      gshiftleft((int)cp->q, cp->basePrime);
 407      int_to_giant(cp->k, tmp);
 408      subg(tmp, cp->basePrime);
 409      returnGiant(tmp);
 410  }
 411  
 412  static int sequalg(int n, giant g) {
 413  	if((g->sign == 1) && (g->n[0] == n)) return(1);
 414  	return(0);
 415  }
 416  
 417  
 418  /*
 419   * Elliptic multiply: x := n * {x, 1}
 420   */
 421  void elliptic_simple(giant x, giant n, curveParams *par) {
 422      giant ztmp = borrowGiant(par->maxDigits);
 423      giant cur_n = borrowGiant(par->maxDigits);
 424  
 425      START_ELL_MEASURE(n);
 426      int_to_giant(1, ztmp);
 427      elliptic(x, ztmp, n, par);
 428      binvg_cp(par, ztmp);
 429      mulg(ztmp, x);
 430      feemod(par, x);
 431      END_ELL_MEASURE;
 432  
 433      returnGiant(cur_n);
 434      returnGiant(ztmp);
 435  }
 436  
 437  /*
 438   * General elliptic multiply.
 439   *
 440   * {xx, zz} := k * {xx, zz}
 441   */
 442  void elliptic(giant xx, giant zz, giant k, curveParams *par) {
 443  	int len = bitlen(k);
 444          int pos = len - 2;
 445          giant xs;
 446          giant zs;
 447          giant xorg;
 448          giant zorg;
 449  
 450  	PROF_START;
 451  	if(sequalg(1,k)) return;
 452  	if(sequalg(2,k)) {
 453  		ell_even(xx, zz, xx, zz, par);
 454  		goto out;
 455  	}
 456          zs = borrowGiant(par->maxDigits);
 457          xs = borrowGiant(par->maxDigits);
 458          zorg = borrowGiant(par->maxDigits);
 459          xorg = borrowGiant(par->maxDigits);
 460  	gtog(xx, xorg); gtog(zz, zorg);
 461  	ell_even(xx, zz, xs, zs, par);
 462  	do {
 463  	   if(bitval(k, pos--)) {
 464  	   	ell_odd(xs, zs, xx, zz, xorg, zorg, par);
 465  		ell_even(xs, zs, xs, zs, par);
 466  	   } else {
 467  	   	ell_odd(xx, zz, xs, zs, xorg, zorg, par);
 468  		ell_even(xx, zz, xx, zz, par);
 469  	   }
 470          } while (pos >= 0);	// REC fix 9/23/94
 471          returnGiant(xs);
 472          returnGiant(zs);
 473          returnGiant(xorg);
 474          returnGiant(zorg);
 475  out:
 476  	PROF_END(ellipticTime);
 477  }
 478  
 479  /*
 480   * Completely rewritten in CryptKit-18, 13 Jan 1997, for new IEEE-style
 481   * curveParameters.
 482   */
 483  void elliptic_add(giant x1, giant x2, giant x3, curveParams *par, int s) {
 484  
 485   /* Addition algorithm for x3 = x1 + x2 on the curve, with sign ambiguity s.
 486      From theory, we know that if {x1,1} and {x2,1} are on a curve, then
 487      their elliptic sum (x1,1} + {x2,1} = {x3,1} must have x3 as one of two
 488      values:
 489  
 490         x3 = U/2 + s*Sqrt[U^2/4 - V]
 491  
 492      where sign s = +-1, and U,V are functions of x1,x2.  Tho present function
 493      is called a maximum of twice, to settle which of +- is s.  When a call
 494      is made, it is guaranteed already that x1, x2 both lie on the same curve
 495      (+- curve); i.e., which curve (+-) is not connected at all with sign s of
 496      the x3 relation.
 497    */
 498  
 499      giant cur_n;
 500      giant t1;
 501      giant t2;
 502      giant t3;
 503      giant t4;
 504      giant t5;
 505  
 506      PROF_START;
 507      cur_n = borrowGiant(par->maxDigits);
 508      t1 = borrowGiant(par->maxDigits);
 509      t2 = borrowGiant(par->maxDigits);
 510      t3 = borrowGiant(par->maxDigits);
 511      t4 = borrowGiant(par->maxDigits);
 512      t5 = borrowGiant(par->maxDigits);
 513  
 514      if(gcompg(x1, x2)==0) {
 515  	int_to_giant(1, t1);
 516  	numer_double(x1, t1, x3, par);
 517  	denom_double(x1, t1, t2, par);
 518  	binvg_cp(par, t2);
 519  	mulg(t2, x3); feemod(par, x3);
 520  	goto out;
 521      }
 522      numer_plus(x1, x2, t1, par);
 523      int_to_giant(1, t3);
 524      numer_times(x1, t3, x2, t3, t2, par);
 525      int_to_giant(1, t4); int_to_giant(1, t5);
 526      denom_times(x1, t4, x2, t5, t3, par);
 527      binvg_cp(par, t3);
 528      mulg(t3, t1); feemod(par, t1); /* t1 := U/2. */
 529      mulg(t3, t2); feemod(par, t2); /* t2 := V. */
 530      /* Now x3 will be t1 +- Sqrt[t1^2 - t2]. */
 531      gtog(t1, t4); gsquare(t4); feemod(par, t4);
 532      subg(t2, t4);
 533      make_base(par, cur_n); iaddg(1, cur_n); gshiftright(2, cur_n);
 534      	/* cur_n := (p+1)/4. */
 535      feepowermodg(par, t4, cur_n);      /* t4 := t2^((p+1)/4) (mod p). */
 536      gtog(t1, x3);
 537      if(s != SIGN_PLUS) negg(t4);
 538      addg(t4, x3);
 539      feemod(par, x3);
 540  
 541  out:
 542      returnGiant(cur_n);
 543      returnGiant(t1);
 544      returnGiant(t2);
 545      returnGiant(t3);
 546      returnGiant(t4);
 547      returnGiant(t5);
 548  
 549      PROF_END(ellAddTime);
 550  }
 551  
 552  /*
 553   * Key exchange atom.
 554   */
 555  giant make_pad(giant privGiant, key publicKey) {
 556      curveParams *par = publicKey->cp;
 557      giant result = newGiant(par->maxDigits);
 558  
 559      gtog(publicKey->x, result);
 560      elliptic_simple(result, privGiant, par);
 561      return result;
 562  }
 563  
 564  static void ell_even(giant x1, giant z1, giant x2, giant z2, curveParams *par) {
 565      giant t1;
 566      giant t2;
 567      giant t3;
 568  
 569      EPROF_START;
 570  
 571      t1 = borrowGiant(par->maxDigits);
 572      t2 = borrowGiant(par->maxDigits);
 573      t3 = borrowGiant(par->maxDigits);
 574  
 575      if(par->curveType == FCT_Montgomery) {
 576  	/* Begin Montgomery OPT: 10 Jan 98 REC. */
 577  	gtog(x1, t1); gsquare(t1); feemod(par, t1); /* t1 := x1^2. */
 578  	gtog(z1, t2); gsquare(t2); feemod(par, t2); /* t2 := z1^2. */
 579  
 580  	gtog(x1, t3); mulg(z1, t3); feemod(par, t3);
 581  	gtog(t3, z2); mulg(par->c, z2); feemod(par, z2);
 582  	addg(t1, z2); addg(t2, z2); mulg(t3, z2); gshiftleft(2, z2);
 583  	feemod(par, z2);  /* z2 := 4 x1 z1 (x1^2 + c x1 z1 + z1^2). */
 584  	gtog(t1, x2); subg(t2, x2); gsquare(x2); feemod(par, x2);
 585  						/* x2 := (x1^2 - z1^2)^2. */
 586  	/* End OPT: 10 Jan 98 REC. */
 587      }
 588      else if((par->curveType == FCT_Weierstrass) && isZero(par->a)) {
 589  	/* Begin Atkin3 OPT: 9 Jan 98 REC. */
 590  	gtog(x1, t1);
 591  	gsquare(t1); feemod(par, t1);
 592  	mulg(x1, t1); feemod(par, t1);   	/* t1 := x^3. */
 593  	gtog(z1, t2);
 594  	gsquare(t2); feemod(par, t2);
 595  	mulg(z1, t2); feemod(par, t2);	/* t2 := z1^3 */
 596  	mulg(par->b, t2); feemod(par, t2); 	/* t2 := b z1^3. */
 597  	gtog(t1, t3); addg(t2, t3);		/* t3 := x^3 + b z1^3 */
 598  	mulg(z1, t3); feemod(par, t3);	/* t3 *= z1
 599  						*     = z1 ( x^3 + b z1^3 ) */
 600  	gshiftleft(2, t3); feemod(par, t3);	/* t3 = 4 z1 (x1^3 + b z1^3) */
 601  
 602  	gshiftleft(3, t2);			/* t2 = 8 b z1^3 */
 603  	subg(t2, t1);			/* t1 = x^3 - 8 b z1^3 */
 604  	mulg(x1, t1); feemod(par, t1);	/* t1 = x1 (x1^3 - 8 b z1^3) */
 605  
 606  	gtog(t3, z2);
 607  	gtog(t1, x2);
 608  	/* End OPT: 9 Jan 98 REC. */
 609      }
 610      else {
 611  	numer_double(x1, z1, t1, par);
 612  	denom_double(x1, z1, t2, par);
 613  	gtog(t1, x2); gtog(t2, z2);
 614      }
 615      returnGiant(t1);
 616      returnGiant(t2);
 617      returnGiant(t3);
 618  
 619      EPROF_END(ellEvenTime);
 620      EPROF_INCR(numEllEvens);
 621  
 622      /*
 623      printf("ell_even end\n");
 624      printf(" x1 : "); printGiant(x1);
 625      printf(" z1 : "); printGiant(z1);
 626      printf(" x2 : "); printGiant(x2);
 627      printf(" z2 : "); printGiant(z2);
 628      */
 629  }
 630  
 631  static void ell_odd(giant x1, giant z1, giant x2, giant z2, giant xxor,
 632  	giant zor, curveParams *par)
 633  {
 634  
 635      giant t1;
 636      giant t2;
 637  
 638      EPROF_START;
 639      t1 = borrowGiant(par->maxDigits);
 640      t2 = borrowGiant(par->maxDigits);
 641  
 642      if(par->curveType == FCT_Montgomery) {
 643  	/* Begin Montgomery OPT: 10 Jan 98 REC. */
 644  	giant t3 = borrowGiant(par->maxDigits);
 645  	giant t4 = borrowGiant(par->maxDigits);
 646  
 647  	gtog(x1, t1); addg(z1, t1);  	/* t1 := x1 + z1. */
 648  	gtog(x2, t2); subg(z2, t2);  	/* t2 := x2 - z2. */
 649  	gtog(x1, t3); subg(z1, t3);  	/* t3 := x1 - z1. */
 650  	gtog(x2, t4); addg(z2, t4);  	/* t4 := x2 + z2. */
 651  	mulg(t2, t1); feemod(par, t1);      /* t1 := (x1 + z1)(x2 - z2) */
 652  	mulg(t4, t3); feemod(par, t3);      /* t4 := (x2 + z2)(x1 - z1) */
 653  	gtog(t1, z2); subg(t3, z2); 	/*???gshiftright(1, z2);*/
 654  			    /* z2 := ((x1 + z1)(x2 - z2) - x2)/2 */
 655  	gsquare(z2); feemod(par,  z2);
 656  	mulg(xxor, z2); feemod(par, z2);
 657  	gtog(t1, x2); addg(t3, x2); 	/*????gshiftright(1, x2);*/
 658  	gsquare(x2); feemod(par, x2);
 659  	mulg(zor, x2); feemod(par, x2);
 660  
 661  	returnGiant(t3);
 662  	returnGiant(t4);
 663      }
 664      else if((par->curveType == FCT_Weierstrass) && isZero(par->a)) {
 665  	/* Begin Atkin3 OPT: 9 Jan 98 REC. */
 666  
 667  	giant t3 = borrowGiant(par->maxDigits);
 668  	giant t4 = borrowGiant(par->maxDigits);
 669  
 670  	gtog(x1, t1); mulg(x2, t1);  feemod(par, t1);  /* t1 := x1 x2. */
 671  	gtog(z1, t2); mulg(z2, t2);  feemod(par, t2);  /* t2 := z1 z2. */
 672  	gtog(x1, t3); mulg(z2, t3);  feemod(par, t3);  /* t3 := x1 z2. */
 673  	gtog(z1, t4); mulg(x2, t4);  feemod(par, t4);  /* t4 := x2 z1. */
 674  	gtog(t3, z2); subg(t4, z2); gsquare(z2); feemod(par, z2);
 675  	mulg(xxor, z2); feemod(par, z2);
 676  	gtog(t1, x2); gsquare(x2); feemod(par, x2);
 677  	addg(t4, t3); mulg(t2, t3); feemod(par, t3);
 678  	mulg(par->b, t3); feemod(par, t3);
 679  	addg(t3, t3); addg(t3, t3);
 680  	subg(t3, x2); mulg(zor, x2); feemod(par, x2);
 681  
 682  	returnGiant(t3);
 683  	returnGiant(t4);
 684  	/* End OPT: 9 Jan 98 REC. */
 685      }
 686      else {
 687  	numer_times(x1, z1, x2, z2, t1, par);
 688  	mulg(zor, t1); feemod(par, t1);
 689  	denom_times(x1, z1, x2, z2, t2, par);
 690  	mulg(xxor, t2); feemod(par, t2);
 691  
 692  	gtog(t1, x2); gtog(t2, z2);
 693      }
 694  
 695      returnGiant(t1);
 696      returnGiant(t2);
 697  
 698      EPROF_END(ellOddTime);
 699      EPROF_INCR(numEllOdds);
 700  
 701      /*
 702      printf("ell_odd end\n");
 703      printf(" x2 : "); printGiant(x2);
 704      printf(" z2 : "); printGiant(z2);
 705      */
 706  }
 707  
 708  /*
 709   * return codes from keys_inconsistent() and signature_compare(). The actual
 710   * values are not public; they are defined here for debugging.
 711   */
 712  #define CURVE_PARAM_MISMATCH	1
 713  #define TWIST_PARAM_MISMATCH 	2
 714  #define SIGNATURE_INVALID 	3
 715  
 716  
 717  /*
 718   * Determine whether two keystructs have compatible parameters (i.e., same
 719   * twist and curveParams). Return 0 if compatible, else non-zero.
 720   */
 721  static int keys_inconsistent(key pub1, key pub2){
 722  	if(!curveParamsEquivalent(pub1->cp, pub2->cp)) {
 723  		return CURVE_PARAM_MISMATCH;
 724  	}
 725  	if(pub1->twist != pub2->twist) {
 726  		return TWIST_PARAM_MISMATCH;
 727  	}
 728  	return 0;
 729  }
 730  
 731  int signature_compare(giant p0x, giant p1x, giant p2x, curveParams *par)
 732  /* Returns non-zero iff p0x cannot be the x-coordinate of the sum of two points whose respective x-coordinates are p1x, p2x. */
 733  {
 734          int ret = 0;
 735          giant t1;
 736  	giant t2;
 737          giant t3;
 738          giant t4;
 739          giant t5;
 740  
 741  	PROF_START;
 742  
 743          t1 = borrowGiant(par->maxDigits);
 744  	t2 = borrowGiant(par->maxDigits);
 745          t3 = borrowGiant(par->maxDigits);
 746          t4 = borrowGiant(par->maxDigits);
 747          t5 = borrowGiant(par->maxDigits);
 748  
 749          if(gcompg(p1x, p2x) == 0) {
 750  		int_to_giant(1, t1);
 751  		numer_double(p1x, t1, t2, par);
 752  		denom_double(p1x, t1, t3, par);
 753  		mulg(p0x, t3); subg(t3, t2);
 754  		feemod(par, t2);
 755          } else {
 756  		numer_plus(p1x, p2x, t1, par);
 757  		gshiftleft(1, t1); feemod(par, t1);
 758  		int_to_giant(1, t3);
 759  		numer_times(p1x, t3, p2x, t3, t2, par);
 760  		int_to_giant(1, t4); int_to_giant(1, t5);
 761  		denom_times(p1x, t4 , p2x, t5, t3, par);
 762  		/* Now we require t3 x0^2 - t1 x0 + t2 == 0. */
 763  		mulg(p0x, t3); feemod(par, t3);
 764  		subg(t1, t3); mulg(p0x, t3);
 765  		feemod(par, t3);
 766  		addg(t3, t2);
 767  		feemod(par, t2);
 768          }
 769  
 770  	if(!isZero(t2)) ret = SIGNATURE_INVALID;
 771          returnGiant(t1);
 772          returnGiant(t2);
 773          returnGiant(t3);
 774          returnGiant(t4);
 775          returnGiant(t5);
 776  	PROF_END(sigCompTime);
 777  	return(ret);
 778  }
 779  
 780  
 781  static void numer_double(giant x, giant z, giant res, curveParams *par)
 782  /* Numerator algebra.
 783     res := (x^2 - a z^2)^2 - 4 b (2 x + c z) z^3.
 784   */
 785  {
 786      giant t1;
 787      giant t2;
 788  
 789      PROF_START;
 790      t1 = borrowGiant(par->maxDigits);
 791      t2 = borrowGiant(par->maxDigits);
 792  
 793      gtog(x, t1); gsquare(t1); feemod(par, t1);
 794      gtog(z, res); gsquare(res); feemod(par, res);
 795      gtog(res, t2);
 796      if(!isZero(par->a) ) {
 797          if(!isone(par->a)) { /* Speedup - REC 17 Jan 1997. */
 798  	    mulg(par->a, res); feemod(par, res);
 799          }
 800          subg(res, t1); feemod(par, t1);
 801      }
 802      gsquare(t1); feemod(par, t1);
 803      /* t1 := (x^2 - a z^2)^2. */
 804      if(isZero(par->b))  {   /* Speedup - REC 17 Jan 1997. */
 805  	gtog(t1, res);
 806          goto done;
 807      }
 808      if(par->curveType != FCT_Weierstrass) {	// i.e., !isZero(par->c)
 809      						// Speedup - REC 17 Jan 1997.
 810  	gtog(z, res); mulg(par->c, res); feemod(par, res);
 811      } else {
 812          int_to_giant(0, res);
 813      }
 814      addg(x, res); addg(x, res); mulg(par->b, res);
 815      feemod(par, res);
 816      gshiftleft(2, res); mulg(z, res); feemod(par, res);
 817      mulg(t2, res); feemod(par, res);
 818      negg(res); addg(t1, res);
 819      feemod(par, res);
 820  
 821  done:
 822      returnGiant(t1);
 823      returnGiant(t2);
 824      PROF_END(numerDoubleTime);
 825  }
 826  
 827  static void numer_plus(giant x1, giant x2, giant res, curveParams *par)
 828  /* Numerator algebra.
 829     res = (x1 x2 + a)(x1 + x2) + 2(c x1 x2 + b).
 830   */
 831  {
 832      giant t1;
 833      giant t2;
 834  
 835      PROF_START;
 836      t1 = borrowGiant(par->maxDigits);
 837      t2 = borrowGiant(par->maxDigits);
 838  
 839      gtog(x1, t1); mulg(x2, t1); feemod(par, t1);
 840      gtog(x2, t2); addg(x1, t2); feemod(par, t2);
 841      gtog(t1, res);
 842      if(!isZero(par->a))
 843      	addg(par->a, res);
 844      mulg(t2, res); feemod(par, res);
 845      if(par->curveType == FCT_Weierstrass) {	// i.e., isZero(par->c)
 846      	int_to_giant(0, t1);
 847      }
 848      else {
 849          mulg(par->c, t1); feemod(par, t1);
 850      }
 851      if(!isZero(par->b))
 852      	addg(par->b, t1);
 853      gshiftleft(1, t1);
 854      addg(t1, res); feemod(par, res);
 855  
 856      returnGiant(t1);
 857      returnGiant(t2);
 858      PROF_END(numerPlusTime);
 859  }
 860  
 861  static void denom_double(giant x, giant z, giant res, curveParams *par)
 862  /* Denominator algebra.
 863      res = 4 z (x^3 + c x^2 z + a x z^2 + b z^3). */
 864  {
 865      giant t1;
 866      giant t2;
 867  
 868      PROF_START;
 869      t1 = borrowGiant(par->maxDigits);
 870      t2 = borrowGiant(par->maxDigits);
 871  
 872      gtog(x, res); gtog(z, t1);
 873      if(par->curveType != FCT_Weierstrass) {	// i.e., !isZero(par->c)
 874  	gtog(par->c, t2); mulg(t1, t2); feemod(par, t2);
 875  	addg(t2, res);
 876      }
 877      mulg(x, res); feemod(par, res);
 878      gsquare(t1); feemod(par, t1);
 879      if(!isZero(par->a)) {
 880  	gtog(t1, t2);
 881      	mulg(par->a, t2); feemod(par, t2);
 882      	addg(t2, res);
 883      }
 884      mulg(x, res); feemod(par, res);
 885      if(!isZero(par->b)) {
 886  	mulg(z, t1); feemod(par, t1);
 887      	mulg(par->b, t1); feemod(par, t1);
 888      	addg(t1, res);
 889      }
 890      mulg(z, res); gshiftleft(2, res);
 891      feemod(par, res);
 892  
 893      returnGiant(t1);
 894      returnGiant(t2);
 895      PROF_END(denomDoubleTime);
 896  }
 897  
 898  
 899  
 900  static void denom_times(giant x1, giant z1, giant x2, giant z2, giant res,
 901  	curveParams *par)
 902  /* Denominator algebra.
 903      res := (x1 z2 - x2 z1)^2
 904   */
 905  {
 906      giant t1;
 907  
 908      PROF_START;
 909      t1 = borrowGiant(par->maxDigits);
 910  
 911      gtog(x1, res); mulg(z2, res); feemod(par, res);
 912      gtog(z1, t1); mulg(x2, t1); feemod(par, t1);
 913      subg(t1, res); gsquare(res); feemod(par, res);
 914  
 915      returnGiant(t1);
 916      PROF_END(denomTimesTime);
 917  }
 918  
 919  static void numer_times(giant x1, giant z1, giant x2, giant z2, giant res,
 920  	curveParams *par)
 921  /* Numerator algebra.
 922      res := (x1 x2 - a z1 z2)^2 -
 923    	          4 b(x1 z2 + x2 z1 + c z1 z2) z1 z2
 924   */
 925  {
 926      giant t1;
 927      giant t2;
 928      giant t3;
 929      giant t4;
 930  
 931      PROF_START;
 932      t1 = borrowGiant(par->maxDigits);
 933      t2 = borrowGiant(par->maxDigits);
 934      t3 = borrowGiant(par->maxDigits);
 935      t4 = borrowGiant(par->maxDigits);
 936  
 937      gtog(x1, t1); mulg(x2, t1); feemod(par, t1);
 938      gtog(z1, t2); mulg(z2, t2); feemod(par, t2);
 939      gtog(t1, res);
 940      if(!isZero(par->a)) {
 941  	gtog(par->a, t3);
 942        	mulg(t2, t3); feemod(par, t3);
 943        	subg(t3, res);
 944      }
 945      gsquare(res); feemod(par, res);
 946      if(isZero(par->b))
 947          goto done;
 948      if(par->curveType != FCT_Weierstrass) {	// i.e., !isZero(par->c)
 949          gtog(par->c, t3);
 950      	mulg(t2, t3); feemod(par, t3);
 951      } else int_to_giant(0, t3);
 952      gtog(z1, t4); mulg(x2, t4); feemod(par, t4);
 953      addg(t4, t3);
 954      gtog(x1, t4); mulg(z2, t4); feemod(par, t4);
 955      addg(t4, t3); mulg(par->b, t3); feemod(par, t3);
 956      mulg(t2, t3); gshiftleft(2, t3); feemod(par, t3);
 957      subg(t3, res);
 958      feemod(par, res);
 959  
 960  done:
 961      returnGiant(t1);
 962      returnGiant(t2);
 963      returnGiant(t3);
 964      returnGiant(t4);
 965      PROF_END(numerTimesTime);
 966  }
 967  
 968  /*
 969   * New, 13 Jan 1997.
 970   */
 971  static void feepowermodg(curveParams *par, giant x, giant n)
 972  /* Power ladder.
 973     x := x^n  (mod 2^q-k)
 974   */
 975  {
 976      int len, pos;
 977      giant t1;
 978  
 979      PROF_START;
 980      t1 = borrowGiant(par->maxDigits);
 981      gtog(x, t1);
 982      int_to_giant(1, x);
 983      len = bitlen(n);
 984      pos = 0;
 985      while(1) {
 986  	if(bitval(n, pos++)) {
 987  	    mulg(t1, x);
 988  	    feemod(par, x);
 989  	}
 990  	if(pos>=len) break;
 991  	gsquare(t1);
 992  	feemod(par, t1);
 993      }
 994      returnGiant(t1);
 995      PROF_END(powerModTime);
 996  }
 997  
 998  /*
 999   * Set g := g mod curveOrder;
1000   * force g to be between 2 and (curveOrder-2), inclusive.
1001   *
1002   * Tolerates zero curve orders (indicating "not known").
1003   */
1004  
1005  /*
1006   * This version is not normally used; it's for when the reciprocal of
1007   * curveOrder is not known and won't be used again.
1008   */
1009  void curveOrderJustify(giant g, giant curveOrder)
1010  {
1011      giant recip;
1012  
1013      CKASSERT(!isZero(curveOrder));
1014  
1015      recip = borrowGiant(2 * abs(g->sign));
1016      make_recip(curveOrder, recip);
1017      curveOrderJustifyWithRecip(g, curveOrder, recip);
1018      returnGiant(recip);
1019  }
1020  
1021  /*
1022   * New optimzation of curveOrderJustify using known reciprocal, 11 June 1997.
1023   * g is set to be within [2, curveOrder-2].
1024   */
1025  static void curveOrderJustifyWithRecip(giant g, giant curveOrder, giant recip)
1026  {
1027      giant tmp;
1028  
1029      CKASSERT(!isZero(curveOrder));
1030  
1031      modg_via_recip(curveOrder, recip, g);	// g now in [0, curveOrder-1]
1032  
1033      if(isZero(g)) {
1034      	/*
1035  	 * First degenerate case - (g == 0) : set g := 2
1036  	 */
1037  	dbgLog(("curveOrderJustify: case 1\n"));
1038     	int_to_giant(2, g);
1039  	return;
1040      }
1041      if(isone(g)) {
1042      	/*
1043  	 * Second case - (g == 1) : set g := 2
1044  	 */
1045   	dbgLog(("curveOrderJustify: case 2\n"));
1046     	int_to_giant(2, g);
1047  	return;
1048      }
1049      tmp = borrowGiant(g->capacity);
1050      gtog(g, tmp);
1051      iaddg(1, tmp);
1052      if(gcompg(tmp, curveOrder) == 0) {
1053      	/*
1054  	 * Third degenerate case - (g == (curveOrder-1)) : set g -= 1
1055  	 */
1056  	dbgLog(("curveOrderJustify: case 3\n"));
1057  	int_to_giant(1, tmp);
1058  	subg(tmp, g);
1059      }
1060      returnGiant(tmp);
1061      return;
1062  }
1063  
1064  #define RECIP_DEBUG	0
1065  #if	RECIP_DEBUG
1066  #define recipLog(x)	printf x
1067  #else	// RECIP_DEBUG
1068  #define recipLog(x)
1069  #endif	// RECIP_DEBUG
1070  
1071  /*
1072   * curveOrderJustify routines with specific orders, using (and possibly
1073   * generating) associated reciprocals.
1074   */
1075  void lesserX1OrderJustify(giant g, curveParams *cp)
1076  {
1077  	/*
1078  	 * Note this is a pointer to a giant in *cp, not a newly
1079  	 * malloc'd giant!
1080  	 */
1081  	giant lesserX1Ord = lesserX1Order(cp);
1082  
1083  	if(isZero(lesserX1Ord)) {
1084  		return;
1085  	}
1086  
1087  	/*
1088  	 * Calculate reciprocal if we don't have it
1089  	 */
1090  	if(isZero(cp->lesserX1OrderRecip)) {
1091  		if((lesserX1Ord == cp->x1OrderPlus) &&
1092  		   (!isZero(cp->x1OrderPlusRecip))) {
1093  		   	/*
1094  			 * lesserX1Ord happens to be x1OrderPlus, and we
1095  			 * have a reciprocal for x1OrderPlus. Copy it over.
1096  			 */
1097  			recipLog((
1098  				"x1OrderPlusRecip --> lesserX1OrderRecip\n"));
1099  		 	gtog(cp->x1OrderPlusRecip, cp->lesserX1OrderRecip);
1100  		}
1101  		else {
1102  			/* Calculate the reciprocal. */
1103  			recipLog(("calc lesserX1OrderRecip\n"));
1104  			make_recip(lesserX1Ord, cp->lesserX1OrderRecip);
1105  		}
1106  	}
1107  	else {
1108  		recipLog(("using existing lesserX1OrderRecip\n"));
1109  	}
1110  	curveOrderJustifyWithRecip(g, lesserX1Ord, cp->lesserX1OrderRecip);
1111  }
1112  
1113  /*
1114   * Common code used by x1OrderPlusJustify() and x1OrderPlusMod() to generate
1115   * reciprocal of x1OrderPlus.
1116   * 8 Sep 1998 - also used by feeSigSign().
1117   */
1118  void calcX1OrderPlusRecip(curveParams *cp)
1119  {
1120  	if(isZero(cp->x1OrderPlusRecip)) {
1121  		if((cp->x1OrderPlus == lesserX1Order(cp)) &&
1122  		   (!isZero(cp->lesserX1OrderRecip))) {
1123  		   	/*
1124  			 * lesserX1Order happens to be x1OrderPlus, and we
1125  			 * have a reciprocal for lesserX1Order. Copy it over.
1126  			 */
1127  			recipLog((
1128  				"lesserX1OrderRecip --> x1OrderPlusRecip\n"));
1129  		 	gtog(cp->lesserX1OrderRecip, cp->x1OrderPlusRecip);
1130  		}
1131  		else {
1132  			/* Calculate the reciprocal. */
1133  			recipLog(("calc x1OrderPlusRecip\n"));
1134  			make_recip(cp->x1OrderPlus, cp->x1OrderPlusRecip);
1135  		}
1136  	}
1137  	else {
1138  		recipLog(("using existing x1OrderPlusRecip\n"));
1139  	}
1140  }
1141  
1142  void x1OrderPlusJustify(giant g, curveParams *cp)
1143  {
1144  	CKASSERT(!isZero(cp->x1OrderPlus));
1145  
1146  	/*
1147  	 * Calculate reciprocal if we don't have it
1148  	 */
1149  	calcX1OrderPlusRecip(cp);
1150  	curveOrderJustifyWithRecip(g, cp->x1OrderPlus, cp->x1OrderPlusRecip);
1151  }
1152  
1153  /*
1154   * g := g mod x1OrderPlus. Result may be zero.
1155   */
1156  void x1OrderPlusMod(giant g, curveParams *cp)
1157  {
1158  	CKASSERT(!isZero(cp->x1OrderPlus));
1159  
1160  	/*
1161  	 * Calculate reciprocal if we don't have it
1162  	 */
1163  	calcX1OrderPlusRecip(cp);
1164  	modg_via_recip(cp->x1OrderPlus, cp->x1OrderPlusRecip, g);
1165  }
1166  
1167  /*
1168   * New general-purpose giant mod routine, 8 Jan 97.
1169   * x := x mod basePrime.
1170   */
1171  
1172  /*
1173   * This stuff is used to analyze the critical loop behavior inside feemod().
1174   */
1175  #define FEEMOD_LOOP_TEST	0
1176  #if	FEEMOD_LOOP_TEST
1177  /*
1178   * these two are only examined via debugger
1179   */
1180  unsigned feemodCalls = 0;		// calls to feemod()
1181  unsigned feemodMultLoops = 0;		// times while() loop runs > once
1182  #define FEEMOD_LOOP_INCR	feemodLoops++
1183  #define FEEMOD_CALL_INCR	feemodCalls++
1184  #else	// FEEMOD_LOOP_TEST
1185  #define FEEMOD_LOOP_INCR
1186  #define FEEMOD_CALL_INCR
1187  #endif	// FEEMOD_LOOP_TEST
1188  
1189  
1190  void
1191  feemod(curveParams *par, giant x)
1192  {
1193      int sign, sign2;
1194      giant t1;
1195      giant t3;
1196      giant t4;
1197      giant t5;
1198      #if		FEEMOD_LOOP_TEST
1199      unsigned feemodLoops = 0;
1200      #endif	// FEEMOD_LOOP_TEST
1201  
1202      FEEMOD_CALL_INCR;		// for FEEMOD_LOOP_TEST
1203      INCR_FEEMODS;		// for ellipticMeasure
1204      PROF_INCR(numFeemod);	// for general profiling
1205  
1206      switch(par->primeType) {
1207          case FPT_Mersenne:
1208  	    /*
1209  	     * Super-optimized Mersenne prime modulus case
1210  	     */
1211  	    gmersennemod(par->q, x);
1212  	    break;
1213  
1214          case FPT_FEE:
1215  	    /*
1216  	     * General 2**q-k case
1217  	     */
1218  	    sign = (x->sign < 0) ? -1 : 1;
1219  	    sign2 = 1;
1220  	    x->sign = abs(x->sign);
1221  	    if(gcompg(par->basePrime, x) >= 0) {
1222  	    	goto outFee;
1223  	    }
1224  	    t1 = borrowGiant(par->maxDigits);
1225  	    t3 = borrowGiant(par->maxDigits);
1226  	    t4 = borrowGiant(par->maxDigits);
1227  	    t5 = borrowGiant(par->maxDigits);
1228  
1229  	    /* Begin OPT: 11 Jan 98 REC. */
1230  	    if( ((par->q & (GIANT_BITS_PER_DIGIT - 1)) == 0) &&
1231  	        (par->k >= 0) &&
1232  		(par->k < GIANT_DIGIT_MASK)) {
1233  
1234  		/*
1235  		 * Microscopic mod for certain regions of {q,k}
1236  		 * parameter space.
1237  		 */
1238  		int j, digits, excess, max;
1239  		giantDigit carry;
1240  		giantDigit termHi;	// double precision
1241  		giantDigit termLo;
1242  		giantDigit *p1, *p2;
1243  
1244  		digits = par->q >> GIANT_LOG2_BITS_PER_DIGIT;
1245  		while(bitlen(x) > par->q) {
1246  		    excess = (x->sign) - digits;
1247  		    max = (excess > digits) ? excess : digits;
1248  		    carry = 0;
1249  
1250  		    p1 = &x->n[0];
1251  		    p2 = &x->n[digits];
1252  
1253  		    if(excess <= digits) {
1254  			carry = VectorMultiply(par->k,
1255  				p2,
1256  				excess,
1257  				p1);
1258  
1259  			/* propagate final carry */
1260  			p1 += excess;
1261  			for(j = excess; j < digits; j++) {
1262  
1263  			    /*
1264  			     * term = *p1 + carry;
1265  			     * *p1++ = term & 65535;
1266  			     * carry = term >> 16;
1267  			     */
1268  			    termLo = giantAddDigits(*p1, carry, &carry);
1269  			    *p1++ = termLo;
1270  			}
1271  		    } else {
1272  			carry = VectorMultiply(par->k,
1273  				p2,
1274  				digits,
1275  				p1);
1276  			p1 += digits;
1277  			p2 += digits;
1278  			for(j = digits; j < excess; j++) {
1279  			    /*
1280  			     * term = (par->k)*(*p2++) + carry;
1281  			     */
1282  			    giantMulDigits(par->k,
1283  			    	*p2++,
1284  				&termLo,
1285  				&termHi);
1286  			    giantAddDouble(&termLo, &termHi, carry);
1287  
1288  			    /*
1289  			     * *p1++ = term & 65535;
1290  			     * carry = term >> 16;
1291  			     */
1292  			    *p1++ = termLo;
1293  			    carry = termHi;
1294  			}
1295  		    }
1296  
1297  		    if(carry > 0) {
1298  			x->n[max] = carry;
1299  		    } else {
1300  			while(--max){
1301  			    if(x->n[max] != 0) break;
1302  			}
1303  		    }
1304  		    x->sign = max + 1;
1305  		    FEEMOD_LOOP_INCR;
1306  		}
1307  	    } else { /* Macroscopic mod for general PT_FEE case. */
1308  		int_to_giant(par->k, t4);
1309  		while(bitlen(x) > par->q) {
1310  		    /* Enter fast loop, as in FEE patent. */
1311  		    int_to_giant(1, t5);
1312  		    gtog(x, t3);
1313  		    extractbits(par->q, t3, x);
1314  		    while(bitlen(t3) > par->q) {
1315  			gshiftright(par->q, t3);
1316  			extractbits(par->q, t3, t1);
1317  			PAUSE_ELL_MEASURE;
1318  			mulg(t4, t5);
1319  			mulg(t5, t1);
1320  			RESUME_ELL_MEASURE;
1321  			addg(t1, x);
1322  		    }
1323  		    FEEMOD_LOOP_INCR;
1324  		}
1325  	    }
1326  	    /* End OPT: 11 Jan 98 REC. */
1327  
1328  	    sign2 = (x->sign < 0)? -1: 1;
1329  	    x->sign = abs(x->sign);
1330  	    returnGiant(t1);
1331  	    returnGiant(t3);
1332  	    returnGiant(t4);
1333  	    returnGiant(t5);
1334  	 outFee:
1335  	    if(gcompg(x, par->basePrime) >=0) subg(par->basePrime, x);
1336  	    if(sign != sign2) {
1337  		    giant t2 = borrowGiant(par->maxDigits);
1338  		    gtog(par->basePrime, t2);
1339  		    subg(x, t2);
1340  		    gtog(t2, x);
1341  		    returnGiant(t2);
1342  	    }
1343  	    break;
1344  	    /* end case PT_FEE */
1345  
1346          case FPT_General:
1347  	    /*
1348  	     * Use brute-force modg.
1349  	     */
1350  	    #if		FEE_DEBUG
1351  	    if(par->basePrimeRecip == NULL) {
1352  	    	CKRaise("feemod(PT_GENERAL): No basePrimeRecip!\n");
1353  	    }
1354  	    #endif	/* FEE_DEBUG */
1355  	    modg_via_recip(par->basePrime, par->basePrimeRecip, x);
1356  	    break;
1357  	    
1358  	case FPT_Default:
1359  	    /* never appears here */
1360  	    CKASSERT(0);
1361  	    break;
1362      } /* switch primeType */
1363  
1364      #if		FEEMOD_LOOP_TEST
1365      if(feemodLoops > 1) {
1366      	feemodMultLoops++;
1367  		if(feemodLoops > 2) {
1368  			printf("feemod while loop executed %d times\n", feemodLoops);
1369  		}
1370      }
1371      #endif	// FEEMOD_LOOP_TEST
1372  
1373      return;
1374  }
1375  
1376  /*
1377   * For a given curveParams, calculate minBytes and maxDigits.
1378   * Assumes valid primeType, and also a valid basePrime for PT_GENERAL.
1379   */
1380  void calcGiantSizes(curveParams *cp)
1381  {
1382  
1383  	if(cp->primeType == FPT_General) {
1384  	    cp->minBytes = (bitlen(cp->basePrime) + 7) / 8;
1385  	}
1386  	else {
1387  	    cp->minBytes = giantMinBytes(cp->q, cp->k);
1388  	}
1389  	cp->maxDigits = giantMaxDigits(cp->minBytes);
1390  }
1391  
1392  unsigned giantMinBytes(unsigned q, int k)
1393  {
1394  	unsigned kIsNeg = (k < 0) ? 1 : 0;
1395  
1396  	return (q + 7 + kIsNeg) / 8;
1397  }
1398  
1399  /*
1400   * min value for "extra" bytes. Derived from the fact that during sig verify,
1401   * we have to multiply a giant representing a digest - which may be
1402   * 20 bytes for SHA1 - by a giant of minBytes.
1403   */
1404  #define MIN_EXTRA_BYTES		20
1405  
1406  unsigned giantMaxDigits(unsigned minBytes)
1407  {
1408  	unsigned maxBytes = 4 * minBytes;
1409  
1410  	if((maxBytes - minBytes) < MIN_EXTRA_BYTES) {
1411  		maxBytes = minBytes + MIN_EXTRA_BYTES;
1412  	}
1413  	return BYTES_TO_GIANT_DIGITS(maxBytes);
1414  }
1415  
1416  
1417  /*
1418   * Optimized binvg(basePrime, x). Avoids the general modg() in favor of
1419   * feemod.
1420   */
1421  int binvg_cp(curveParams *cp, giant x)
1422  {
1423  	feemod(cp, x);
1424  	return(binvaux(cp->basePrime, x));
1425  }
1426  
1427  /*
1428   * Optimized binvg(x1OrderPlus, x). Uses x1OrderPlusMod().
1429   */
1430  int binvg_x1OrderPlus(curveParams *cp, giant x)
1431  {
1432  	x1OrderPlusMod(x, cp);
1433  	return(binvaux(cp->x1OrderPlus, x));
1434  }