+++ /dev/null
-/* 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.
- ***************************************************************************
-
- giantIntegers.c - library for large-integer arithmetic.
-
- Revision History
- ----------------
- Fixed a==b bug in addg().
- 10/06/98 ap
- Changed to compile with C++.
- 13 Apr 98 Fixed shiftright(1) bug in modg_via_recip.
- 09 Apr 98 at Apple
- Major rewrite of core arithmetic routines to make this module
- independent of size of giantDigit.
- Removed idivg() and radixdiv().
- 20 Jan 98 at Apple
- Deleted FFT arithmetic; simplified mulg().
- 09 Jan 98 at Apple
- gshiftright() optimization.
- 08 Jan 98 at Apple
- newGiant() returns NULL on malloc failure
- 24 Dec 97 at Apple
- New grammarSquare(); optimized modg_via_recip()
- 11 Jun 97 at Apple
- Added modg_via_recip(), divg_via_recip(), make_recip()
- Added new multiple giant stack mechanism
- Fixed potential packing/alignment bug in copyGiant()
- Added profiling for borrowGiant(), returnGiant()
- Deleted obsolete ifdef'd code
- Deleted newgiant()
- All calls to borrowGiant() now specify required size (no more
- borrowGiant(0) calls)
- 08 May 97 at Apple
- Changed size of giantstruct.n to 1 for Mac build
- 05 Feb 97 at Apple
- newGiant() no longer modifies CurrentMaxShorts or giant stack
- Added modg profiling
- 01 Feb 97 at NeXT
- Added iszero() check in gcompg
- 17 Jan 97 at NeXT
- Fixed negation bug in gmersennemod()
- Fixed n[words-1] == 0 bug in extractbits()
- Cleaned up lots of static declarations
- 19 Sep 96 at NeXT
- Fixed --size underflow bug in normal_subg().
- 4 Sep 96 at NeXT
- Fixed (b<n), (sign<0) case in gmersennemod() to allow for arbitrary n.
- 9 Aug 96 at NeXT
- Fixed sign-extend bug in data_to_giant().
- 7 Aug 96 at NeXT
- Changed precision in newtondivide().
- Removed #ifdef UNUSED code.
- 24 Jul 96 at NeXT
- Added scompg().
- Created.
-
-*/
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include "platform.h"
-#include "giantIntegers.h"
-#include "feeDebug.h"
-#include "ckconfig.h"
-#include "ellipticMeasure.h"
-#include "falloc.h"
-#include "giantPortCommon.h"
-
-#ifdef FEE_DEBUG
-#if (GIANT_LOG2_BITS_PER_DIGIT == 4)
-#warning Compiling with two-byte giantDigits
-#endif
-#endif
-
-#if 0
-#if FEE_DEBUG
-char printbuf1[200];
-char printbuf2[200];
-char printbuf3[200];
-void printGiantBuf(giant x)
-{
- int i;
-
- sprintf(printbuf2, "sign=%d cap=%d n[]=", x->sign, x->capacity);
- for(i=0; i<abs(x->sign); i++) {
- sprintf(printbuf3 + 10*i, "%u:", x->n[i]);
- }
-}
-
-char printbuf4[200];
-char printbuf5[200];
-void printGiantBuf2(giant x)
-{
- int i;
-
- sprintf(printbuf4, "sign=%d cap=%d n[]=", x->sign, x->capacity);
- for(i=0; i<abs(x->sign); i++) {
- sprintf(printbuf5 + 10*i, "%u:", x->n[i]);
- }
-}
-#endif /* FEE_DEBUG */
-#endif /* 0 */
-
-/******** debugging flags *********/
-
-/*
- * Flag use of unoptimized divg, modg, binvg
- */
-//#define WARN_UNOPTIMIZE FEE_DEBUG
-#define WARN_UNOPTIMIZE 0
-
-/*
- * Log interesting giant stack events
- */
-#define LOG_GIANT_STACK 0
-
-/*
- * Log allocation of giant larger than stack size
- */
-#define LOG_GIANT_STACK_OVERFLOW 1
-
-/*
- * Flag newGiant(0) and borrowGiant(0) calls
- */
-#define WARN_ZERO_GIANT_SIZE FEE_DEBUG
-
-/* temp mac-only giant initialization debug */
-#define GIANT_MAC_DEBUG 0
-#if GIANT_MAC_DEBUG
-
-#include <string.h>
-#include <TextUtils.h>
-
-/* this one needs a writable string */
-static void logCom(unsigned char *str) {
- c2pstr((char *)str);
- DebugStr(str);
-}
-
-/* constant strings */
-void dblog0(const char *str) {
- Str255 outStr;
- strcpy((char *)outStr, str);
- logCom(outStr);
-}
-
-#else
-#define dblog0(s)
-
-#endif /* GIANT_MAC_DEBUG */
-
-#ifndef min
-#define min(a,b) ((a)<(b)? (a) : (b))
-#endif // min
-#ifndef max
-#define max(a,b) ((a)>(b)? (a) : (b))
-#endif // max
-
-#ifndef TRUE
-#define TRUE 1
-#endif // TRUE
-#ifndef FALSE
-#define FALSE 0
-#endif // FALSE
-
-static void absg(giant g); /* g := |g|. */
-
-/************** globals *******************/
-
-
-/* ------ giant stack package ------ */
-
-/*
- * The giant stack package is a local cache which allows us to avoid calls
- * to malloc() for borrowGiant(). On a 90 Mhz Pentium, enabling the
- * giant stack package shows about a 1.35 speedup factor over an identical
- * CryptKit without the giant stacks enabled.
- */
-
-#if GIANTS_VIA_STACK
-
-#if LOG_GIANT_STACK
-#define gstackDbg(x) printf x
-#else // LOG_GIANT_STACK
-#define gstackDbg(x)
-#endif // LOG_GIANT_STACK
-
-typedef struct {
- unsigned numDigits; // capacity of giants in this stack
- unsigned numFree; // number of free giants in stack
- unsigned totalGiants; // total number in *stack
- giant *stack;
-} gstack;
-
-static gstack *gstacks = NULL; // array of stacks
-static unsigned numGstacks = 0; // # of elements in gstacks
-static int gstackInitd = 0; // this module has been init'd
-
-#define INIT_NUM_GIANTS 16 /* initial # of giants / stack */
-#define MIN_GIANT_SIZE 4 /* numDigits for gstack[0] */
-#define GIANT_SIZE_INCR 2 /* in << bits */
-
-/*
- * Initialize giant stacks, with up to specified max giant size.
- */
-void initGiantStacks(unsigned maxDigits)
-{
- unsigned curSize = MIN_GIANT_SIZE;
- unsigned sz;
- unsigned i;
-
- dblog0("initGiantStacks\n");
-
- if(gstackInitd) {
- /*
- * Shouldn't be called more than once...
- */
- printf("multiple initGiantStacks calls\n");
- return;
- }
- gstackDbg(("initGiantStacks(%d)\n", maxDigits));
-
- /*
- * How many stacks?
- */
- numGstacks = 1;
- while(curSize<=maxDigits) {
- curSize <<= GIANT_SIZE_INCR;
- numGstacks++;
- }
-
- sz = sizeof(gstack) * numGstacks;
- gstacks = (gstack*) fmalloc(sz);
- bzero(gstacks, sz);
-
- curSize = MIN_GIANT_SIZE;
- for(i=0; i<numGstacks; i++) {
- gstacks[i].numDigits = curSize;
- curSize <<= GIANT_SIZE_INCR;
- }
-
- gstackInitd = 1;
-}
-
-/* called at shut down - free resources */
-void freeGiantStacks(void)
-{
- int i;
- int j;
- gstack *gs;
-
- if(!gstackInitd) {
- return;
- }
- for(i=0; i<numGstacks; i++) {
- gs = &gstacks[i];
- for(j=0; j<gs->numFree; j++) {
- freeGiant(gs->stack[j]);
- gs->stack[j] = NULL;
- }
- /* and the stack itself - may be null if this was never used */
- if(gs->stack != NULL) {
- ffree(gs->stack);
- gs->stack = NULL;
- }
- }
- ffree(gstacks);
- gstacks = NULL;
- gstackInitd = 0;
-}
-
-#endif // GIANTS_VIA_STACK
-
-giant borrowGiant(unsigned numDigits)
-{
- giant result;
-
- #if GIANTS_VIA_STACK
-
- unsigned stackNum;
- gstack *gs = gstacks;
-
- #if WARN_ZERO_GIANT_SIZE
- if(numDigits == 0) {
- printf("borrowGiant(0)\n");
- numDigits = gstacks[numGstacks-1].numDigits;
- }
- #endif // WARN_ZERO_GIANT_SIZE
-
- /*
- * Find appropriate stack
- */
- if(numDigits <= MIN_GIANT_SIZE)
- stackNum = 0;
- else if (numDigits <= (MIN_GIANT_SIZE << GIANT_SIZE_INCR))
- stackNum = 1;
- else if (numDigits <= (MIN_GIANT_SIZE << (2 * GIANT_SIZE_INCR)))
- stackNum = 2;
- else if (numDigits <= (MIN_GIANT_SIZE << (3 * GIANT_SIZE_INCR)))
- stackNum = 3;
- else if (numDigits <= (MIN_GIANT_SIZE << (4 * GIANT_SIZE_INCR)))
- stackNum = 4;
- else
- stackNum = numGstacks;
-
- if(stackNum >= numGstacks) {
- /*
- * out of bounds; just malloc
- */
- #if LOG_GIANT_STACK_OVERFLOW
- gstackDbg(("giantFromStack overflow; numDigits %d\n",
- numDigits));
- #endif // LOG_GIANT_STACK_OVERFLOW
- return newGiant(numDigits);
- }
- gs = &gstacks[stackNum];
-
- #if GIANT_MAC_DEBUG
- if((gs->numFree != 0) && (gs->stack == NULL)) {
- dblog0("borrowGiant: null stack!\n");
- }
- #endif
-
- if(gs->numFree != 0) {
- result = gs->stack[--gs->numFree];
- }
- else {
- /*
- * Stack empty; malloc
- */
- result = newGiant(gs->numDigits);
- }
-
- #else /* GIANTS_VIA_STACK */
-
- result = newGiant(numDigits);
-
- #endif /* GIANTS_VIA_STACK */
-
- PROF_INCR(numBorrows);
- return result;
-}
-
-void returnGiant(giant g)
-{
-
- #if GIANTS_VIA_STACK
-
- unsigned stackNum;
- gstack *gs;
- unsigned cap = g->capacity;
-
-
- #if FEE_DEBUG
- if(!gstackInitd) {
- CKRaise("returnGiant before stacks initialized!");
- }
- #endif // FEE_DEBUG
-
- #if GIANT_MAC_DEBUG
- if(g == NULL) {
- dblog0("returnGiant: null g!\n");
- }
- #endif
-
- /*
- * Find appropriate stack. Note we expect exact match of
- * capacity and stack's giant size.
- */
- /*
- * Optimized unrolled loop. Just make sure there are enough cases
- * to handle all of the stacks. Errors in this case will be flagged
- * via LOG_GIANT_STACK_OVERFLOW.
- */
- switch(cap) {
- case MIN_GIANT_SIZE:
- stackNum = 0;
- break;
- case MIN_GIANT_SIZE << GIANT_SIZE_INCR:
- stackNum = 1;
- break;
- case MIN_GIANT_SIZE << (2 * GIANT_SIZE_INCR):
- stackNum = 2;
- break;
- case MIN_GIANT_SIZE << (3 * GIANT_SIZE_INCR):
- stackNum = 3;
- break;
- case MIN_GIANT_SIZE << (4 * GIANT_SIZE_INCR):
- stackNum = 4;
- break;
- default:
- stackNum = numGstacks;
- break;
- }
-
- if(stackNum >= numGstacks) {
- /*
- * out of bounds; just free
- */
- #if LOG_GIANT_STACK_OVERFLOW
- gstackDbg(("giantToStack overflow; numDigits %d\n", cap));
- #endif // LOG_GIANT_STACK_OVERFLOW
- freeGiant(g);
- return;
- }
- gs = &gstacks[stackNum];
- if(gs->numFree == gs->totalGiants) {
- if(gs->totalGiants == 0) {
- gstackDbg(("Initial alloc of gstack(%d)\n",
- gs->numDigits));
- gs->totalGiants = INIT_NUM_GIANTS;
- }
- else {
- gs->totalGiants *= 2;
- gstackDbg(("Bumping gstack(%d) to %d\n",
- gs->numDigits, gs->totalGiants));
- }
- gs->stack = (giantstruct**) frealloc(gs->stack, gs->totalGiants*sizeof(giant));
- }
- g->sign = 0; // not sure this is important...
- gs->stack[gs->numFree++] = g;
-
- #if GIANT_MAC_DEBUG
- if((gs->numFree != 0) && (gs->stack == NULL)) {
- dblog0("borrowGiant: null stack!\n");
- }
- #endif
-
- #else /* GIANTS_VIA_STACK */
-
- freeGiant(g);
-
- #endif /* GIANTS_VIA_STACK */
-}
-
-void freeGiant(giant x) {
- ffree(x);
-}
-
-giant newGiant(unsigned numDigits) {
- // giant sufficient for 2^numbits+16 sized ops
- int size;
- giant result;
-
- #if WARN_ZERO_GIANT_SIZE
- if(numDigits == 0) {
- printf("newGiant(0)\n");
- #if GIANTS_VIA_STACK
- numDigits = gstacks[numGstacks-1].totalGiants;
- #else
- /* HACK */
- numDigits = 20;
- #endif
- }
- #endif // WARN_ZERO_GIANT_SIZE
-
- size = (numDigits-1) * GIANT_BYTES_PER_DIGIT + sizeof(giantstruct);
- result = (giant)fmalloc(size);
- if(result == NULL) {
- return NULL;
- }
- result->sign = 0;
- result->capacity = numDigits;
- return result;
-}
-
-giant copyGiant(giant x)
-{
- int bytes;
-
- giant result = newGiant(x->capacity);
-
- /*
- * 13 Jun 1997
- * NO! this assumes packed alignment
- */
- bytes = sizeof(giantstruct) +
- ((x->capacity - 1) * GIANT_BYTES_PER_DIGIT);
- bcopy(x, result, bytes);
- return result;
-}
-
-/* ------ initialization and utility routines ------ */
-
-
-unsigned bitlen(giant n) {
- unsigned b = GIANT_BITS_PER_DIGIT;
- giantDigit c = 1 << (GIANT_BITS_PER_DIGIT - 1);
- giantDigit w;
-
- if (isZero(n)) {
- return(0);
- }
- w = n->n[abs(n->sign) - 1];
- if (!w) {
- CKRaise("bitlen - no bit set!");
- }
- while((w&c) == 0) {
- b--;
- c >>= 1;
- }
- return(GIANT_BITS_PER_DIGIT * (abs(n->sign)-1) + b);
-}
-
-int bitval(giant n, int pos) {
- int i = abs(pos) >> GIANT_LOG2_BITS_PER_DIGIT;
- giantDigit c = 1 << (pos & (GIANT_BITS_PER_DIGIT - 1));
-
- return((n->n[i]) & c);
-}
-
-int gsign(giant g)
-/* returns the sign of g */
-{
- if (isZero(g)) return(0);
- if (g->sign > 0) return(1);
- return(-1);
-}
-
-/*
- * Adjust sign for possible leading (m.s.) zero digits
- */
-void gtrimSign(giant g)
-{
- int numDigits = abs(g->sign);
- int i;
-
- for(i=numDigits-1; i>=0; i--) {
- if(g->n[i] == 0) {
- numDigits--;
- }
- else {
- break;
- }
- }
- if(g->sign < 0) {
- g->sign = -numDigits;
- }
- else {
- g->sign = numDigits;
- }
-}
-
-
-int isone(giant g) {
- return((g->sign==1)&&(g->n[0]==1));
-}
-
-int isZero(giant thegiant) {
-/* Returns TRUE if thegiant == 0. */
- int count;
- int length = abs(thegiant->sign);
- giantDigit *numpointer;
-
- if (length) {
- numpointer = thegiant->n;
-
- for(count = 0; count<length; ++count,++numpointer)
- if (*numpointer != 0 ) return(FALSE);
- }
- return(TRUE);
-}
-
-int gcompg(giant a, giant b)
-/* returns -1,0,1 if a<b, a=b, a>b, respectively */
-{
- int sa = a->sign;
- int j;
- int sb = b->sign;
- giantDigit va;
- giantDigit vb;
- int sgn;
-
- if(isZero(a) && isZero(b)) return 0;
- if(sa > sb) return(1);
- if(sa < sb) return(-1);
- if(sa < 0) {
- sa = -sa; /* Take absolute value of sa */
- sgn = -1;
- } else sgn = 1;
- for(j = sa-1; j >= 0; j--) {
- va = a->n[j]; vb = b->n[j];
- if (va > vb) return(sgn);
- if (va < vb) return(-sgn);
- }
- return(0);
-}
-
-/* destgiant becomes equal to srcgiant */
-void gtog(giant srcgiant, giant destgiant) {
-
- int numbytes;
-
- CKASSERT(srcgiant != NULL);
- numbytes = abs(srcgiant->sign) * GIANT_BYTES_PER_DIGIT;
- if (destgiant->capacity < abs(srcgiant->sign))
- CKRaise("gtog overflow!!");
- memcpy((char *)destgiant->n, (char *)srcgiant->n, numbytes);
- destgiant->sign = srcgiant->sign;
-}
-
-void int_to_giant(int i, giant g) {
-/* The giant g becomes set to the integer value i. */
- int isneg = (i<0);
- unsigned int j = abs(i);
- unsigned dex;
-
- g->sign = 0;
- if (i==0) {
- g->n[0] = 0;
- return;
- }
-
- if(GIANT_BYTES_PER_DIGIT == sizeof(int)) {
- g->n[0] = j;
- g->sign = 1;
- }
- else {
- /* one loop per digit */
- unsigned scnt = GIANT_BITS_PER_DIGIT; // fool compiler
-
- for(dex=0; dex<sizeof(int); dex++) {
- g->n[dex] = j & GIANT_DIGIT_MASK;
- j >>= scnt;
- g->sign++;
- if(j == 0) {
- break;
- }
- }
- }
- if (isneg) {
- g->sign = -(g->sign);
- }
-}
-
-/*------------- Arithmetic --------------*/
-
-void negg(giant g) {
-/* g becomes -g */
- g->sign = -g->sign;
-}
-
-void iaddg(int i, giant g) { /* positive g becomes g + (int)i */
- int j;
- giantDigit carry;
- int size = abs(g->sign);
-
- if (isZero(g)) {
- int_to_giant(i,g);
- }
- else {
- carry = i;
- for(j=0; ((j<size) && (carry != 0)); j++) {
- g->n[j] = giantAddDigits(g->n[j], carry, &carry);
- }
- if(carry) {
- ++g->sign;
- // realloc
- if (g->sign > (int)g->capacity) CKRaise("iaddg overflow!");
- g->n[size] = carry;
- }
- }
-}
-
-/*
- * g *= (int n)
- *
- * FIXME - we can improve this...
- */
-void imulg(unsigned n, giant g)
-{
- giant tmp = borrowGiant(abs(g->sign) + sizeof(int));
-
- int_to_giant(n, tmp);
- mulg(tmp, g);
- returnGiant(tmp);
-}
-
-static void normal_addg(giant a, giant b)
-/* b := a + b, both a,b assumed non-negative. */
-{
- giantDigit carry1 = 0;
- giantDigit carry2 = 0;
- int asize = a->sign, bsize = b->sign;
- giantDigit *an = a->n;
- giantDigit *bn = b->n;
- giantDigit tmp;
- int j;
- int comSize;
- int maxSize;
-
- if(asize < bsize) {
- comSize = asize;
- maxSize = bsize;
- }
- else {
- comSize = bsize;
- maxSize = asize;
- }
-
- /* first handle the common digits */
- for(j=0; j<comSize; j++) {
- /*
- * first add the carry, then an[j] - either add could result
- * in another carry
- */
- if(carry1 || carry2) {
- tmp = giantAddDigits(bn[j], (giantDigit)1, &carry1);
- }
- else {
- carry1 = 0;
- tmp = bn[j];
- }
- bn[j] = giantAddDigits(tmp, an[j], &carry2);
- }
-
- if(asize < bsize) {
-
- /* now propagate remaining carry beyond asize */
- if(carry2) {
- carry1 = 1;
- }
- if(carry1) {
- for(; j<bsize; j++) {
- bn[j] = giantAddDigits(bn[j], (giantDigit)1, &carry1);
- if(carry1 == 0) {
- break;
- }
- }
- }
- } else {
- /* now propagate remaining an[] and carry beyond bsize */
- if(carry2) {
- carry1 = 1;
- }
- for(; j<asize; j++) {
- if(carry1) {
- bn[j] = giantAddDigits(an[j], (giantDigit)1, &carry1);
- }
- else {
- bn[j] = an[j];
- carry1 = 0;
- }
- }
- }
- b->sign = maxSize;
- if(carry1) {
- // realloc?
- bn[j] = 1;
- b->sign++;
- if (b->sign > (int)b->capacity) CKRaise("iaddg overflow!");
- }
-
-}
-
-static void normal_subg(giant a, giant b)
-/* b := b - a; requires b, a non-negative and b >= a. */
-{
- int j;
- int size = b->sign;
- giantDigit tmp;
- giantDigit borrow1 = 0;
- giantDigit borrow2 = 0;
- giantDigit *an = a->n;
- giantDigit *bn = b->n;
-
- if(a->sign == 0) {
- return;
- }
-
- for (j=0; j<a->sign; ++j) {
- if(borrow1 || borrow2) {
- tmp = giantSubDigits(bn[j], (giantDigit)1, &borrow1);
- }
- else {
- tmp = bn[j];
- borrow1 = 0;
- }
- bn[j] = giantSubDigits(tmp, an[j], &borrow2);
- }
- if(borrow1 || borrow2) {
- /* propagate borrow thru remainder of bn[] */
- borrow1 = 1;
- for (j=a->sign; j<size; ++j) {
- if(borrow1) {
- bn[j] = giantSubDigits(bn[j], (giantDigit)1, &borrow1);
- }
- else {
- break;
- }
- }
- }
-
- /* adjust sign for leading zero digits */
- while((size-- > 0) && (b->n[size] == 0))
- ;
- b->sign = (b->n[size] == 0)? 0 : size+1;
-}
-
-static void reverse_subg(giant a, giant b)
-/* b := a - b; requires b, a non-negative and a >= b. */
-{
- int j;
- int size = a->sign;
- giantDigit tmp;
- giantDigit borrow1 = 0;
- giantDigit borrow2 = 0;
- giantDigit *an = a->n;
- giantDigit *bn = b->n;
-
- if(b->sign == 0) {
- gtog(a, b);
- return;
- }
- for (j=0; j<b->sign; ++j) {
- if(borrow1 || borrow2) {
- tmp = giantSubDigits(an[j], (giantDigit)1, &borrow1);
- }
- else {
- tmp = an[j];
- borrow1 = 0;
- }
- bn[j] = giantSubDigits(tmp, bn[j], &borrow2);
- }
- if(borrow1 || borrow2) {
- /* propagate borrow thru remainder of bn[] */
- borrow1 = 1;
- }
- for (j=b->sign; j<size; ++j) {
- if(borrow1) {
- bn[j] = giantSubDigits(an[j], (giantDigit)1, &borrow1);
- }
- else {
- bn[j] = an[j];
- borrow1 = 0;
- }
- }
-
- b->sign = size; /* REC, 21 Apr 1996. */
- while(!b->n[--size]);
- b->sign = size+1;
-}
-
-
-void addg(giant a, giant b)
-/* b := b + a, any signs any result. */
-{ int asgn = a->sign, bsgn = b->sign;
- if(asgn == 0) return;
- if(bsgn == 0) {
- gtog(a,b);
- return;
- }
- if((asgn < 0) == (bsgn < 0)) {
- if(bsgn > 0) {
- normal_addg(a,b);
- return;
- }
- negg(a); if(a != b) negg(b); normal_addg(a,b); /* Fix REC 1 Dec 98. */
- negg(a); if(a != b) negg(b); return; /* Fix REC 1 Dec 98. */
- }
- if(bsgn > 0) {
- negg(a);
- if(gcompg(b,a) >= 0) {
- normal_subg(a,b);
- negg(a);
- return;
- }
- reverse_subg(a,b);
- negg(a);
- negg(b);
- return;
- }
- negg(b);
- if(gcompg(b,a) < 0) {
- reverse_subg(a,b);
- return;
- }
- normal_subg(a,b);
- negg(b);
- return;
-}
-
-void subg(giant a, giant b)
-/* b := b - a, any signs, any result. */
-{
- int asgn = a->sign, bsgn = b->sign;
- if(asgn == 0) return;
- if(bsgn == 0) {
- gtog(a,b);
- negg(b);
- return;
- }
- if((asgn < 0) != (bsgn < 0)) {
- if(bsgn > 0) {
- negg(a);
- normal_addg(a,b);
- negg(a);
- return;
- }
- negg(b);
- normal_addg(a,b);
- negg(b);
- return;
- }
- if(bsgn > 0) {
- if(gcompg(b,a) >= 0) {
- normal_subg(a,b);
- return;
- }
- reverse_subg(a,b);
- negg(b);
- return;
- }
- negg(a); negg(b);
- if(gcompg(b,a) >= 0) {
- normal_subg(a,b);
- negg(a);
- negg(b);
- return;
- }
- reverse_subg(a,b);
- negg(a);
- return;
-}
-
-static void bdivg(giant v, giant u)
-/* u becomes greatest power of two not exceeding u/v. */
-{
- int diff = bitlen(u) - bitlen(v);
- giant scratch7;
-
- if (diff<0) {
- int_to_giant(0,u);
- return;
- }
- scratch7 = borrowGiant(u->capacity);
- gtog(v, scratch7);
- gshiftleft(diff,scratch7);
- if(gcompg(u,scratch7) < 0) diff--;
- if(diff<0) {
- int_to_giant(0,u);
- returnGiant(scratch7);
- return;
- }
- int_to_giant(1,u);
- gshiftleft(diff,u);
- returnGiant(scratch7);
-}
-
-int binvaux(giant p, giant x)
-/* Binary inverse method.
- Returns zero if no inverse exists, in which case x becomes
- GCD(x,p). */
-{
- giant scratch7;
- giant u0;
- giant u1;
- giant v0;
- giant v1;
- int result = 1;
- int giantSize;
- PROF_START;
-
- if(isone(x)) return(result);
- giantSize = 4 * abs(p->sign);
- scratch7 = borrowGiant(giantSize);
- u0 = borrowGiant(giantSize);
- u1 = borrowGiant(giantSize);
- v0 = borrowGiant(giantSize);
- v1 = borrowGiant(giantSize);
- int_to_giant(1, v0); gtog(x, v1);
- int_to_giant(0,x); gtog(p, u1);
- while(!isZero(v1)) {
- gtog(u1, u0); bdivg(v1, u0);
- gtog(x, scratch7);
- gtog(v0, x);
- mulg(u0, v0);
- subg(v0,scratch7);
- gtog(scratch7, v0);
-
- gtog(u1, scratch7);
- gtog(v1, u1);
- mulg(u0, v1);
- subg(v1,scratch7);
- gtog(scratch7, v1);
- }
- if (!isone(u1)) {
- gtog(u1,x);
- if(x->sign<0) addg(p, x);
- result = 0;
- goto done;
- }
- if (x->sign<0) addg(p, x);
- done:
- returnGiant(scratch7);
- returnGiant(u0);
- returnGiant(u1);
- returnGiant(v0);
- returnGiant(v1);
- PROF_END(binvauxTime);
- return(result);
-}
-
-/*
- * Superceded by binvg_cp()
- */
-#if 0
-int binvg(giant p, giant x)
-{
- modg(p, x);
- return(binvaux(p,x));
-}
-#endif
-
-static void absg(giant g) {
-/* g becomes the absolute value of g */
- if (g->sign < 0) g->sign = -g->sign;
-}
-
-void gshiftleft(int bits, giant g) {
-/* shift g left bits bits. Equivalent to g = g*2^bits */
- int rem = bits & (GIANT_BITS_PER_DIGIT - 1);
- int crem = GIANT_BITS_PER_DIGIT - rem;
- int digits = 1 + (bits >> GIANT_LOG2_BITS_PER_DIGIT);
- int size = abs(g->sign);
- int j;
- int k;
- int sign = gsign(g);
- giantDigit carry;
- giantDigit dat;
-
- #if FEE_DEBUG
- if(bits < 0) {
- CKRaise("gshiftleft(-bits)\n");
- }
- #endif /* FEE_DEBUG */
-
- if(!bits) return;
- if(!size) return;
- if((size+digits) > (int)g->capacity) {
- CKRaise("gshiftleft overflow");
- return;
- }
- k = size - 1 + digits; // (MSD of result + 1)
- carry = 0;
-
- /* bug fix for 32-bit giantDigits; this is also an optimization for
- * other sizes. rem=0 means we're shifting strictly by digits, no
- * bit shifts. */
- if(rem == 0) {
- g->n[k] = 0; // XXX hack - for sign fixup
- for(j=size-1; j>=0; j--) {
- g->n[--k] = g->n[j];
- }
- do{
- g->n[--k] = 0;
- } while(k>0);
- }
- else {
- /*
- * normal unaligned case
- * FIXME - this writes past g->n[size-1] the first time thru!
- */
- for(j=size-1; j>=0; j--) {
- dat = g->n[j];
- g->n[k--] = (dat >> crem) | carry;
- carry = (dat << rem);
- }
- do{
- g->n[k--] = carry;
- carry = 0;
- } while(k>=0);
- }
- k = size - 1 + digits;
- if(g->n[k] == 0) --k;
- g->sign = sign * (k+1);
- if (abs(g->sign) > g->capacity) {
- CKRaise("gshiftleft overflow");
- }
-}
-
-void gshiftright(int bits, giant g) {
-/* shift g right bits bits. Equivalent to g = g/2^bits */
- int j;
- int size=abs(g->sign);
- giantDigit carry;
- int digits = bits >> GIANT_LOG2_BITS_PER_DIGIT;
- int remain = bits & (GIANT_BITS_PER_DIGIT - 1);
- int cremain = GIANT_BITS_PER_DIGIT - remain;
-
- #if FEE_DEBUG
- if(bits < 0) {
- CKRaise("gshiftright(-bits)\n");
- }
- #endif /* FEE_DEBUG */
- if(bits==0) return;
- if(isZero(g)) return;
- if (digits >= size) {
- g->sign = 0;
- return;
- }
-
- size -= digits;
-
-/* Begin OPT: 9 Jan 98 REC. */
- if(remain == 0) {
- if(g->sign > 0) {
- g->sign = size;
- }
- else {
- g->sign = -size;
- }
- for(j=0; j < size; j++) {
- g->n[j] = g->n[j+digits];
- }
- return;
- }
-/* End OPT: 9 Jan 98 REC. */
-
- for(j=0;j<size;++j) {
- if (j==size-1) {
- carry = 0;
- }
- else {
- carry = (g->n[j+digits+1]) << cremain;
- }
- g->n[j] = ((g->n[j+digits]) >> remain ) | carry;
- }
- if (g->n[size-1] == 0) {
- --size;
- }
- if(g->sign > 0) {
- g->sign = size;
- }
- else {
- g->sign = -size;
- }
- if (abs(g->sign) > g->capacity) {
- CKRaise("gshiftright overflow");
- }
-}
-
-
-void extractbits(unsigned n, giant src, giant dest) {
-/* dest becomes lowermost n bits of src. Equivalent to dest = src % 2^n */
- int digits = n >> GIANT_LOG2_BITS_PER_DIGIT;
- int numbytes = digits * GIANT_BYTES_PER_DIGIT;
- int bits = n & (GIANT_BITS_PER_DIGIT - 1);
-
- if (n <= 0) {
- return;
- }
- if (dest->capacity * 8 * GIANT_BYTES_PER_DIGIT < n) {
- CKRaise("extractbits - not enough room");
- }
- if (digits >= abs(src->sign)) {
- gtog(src,dest);
- }
- else {
- memcpy((char *)(dest->n), (char *)(src->n), numbytes);
- if (bits) {
- dest->n[digits] = src->n[digits] & ((1<<bits)-1);
- ++digits;
- }
- /* Next, fix by REC, 12 Jan 97. */
- // while((dest->n[words-1] == 0) && (words > 0)) --words;
- while((digits > 0) && (dest->n[digits-1] == 0)) {
- --digits;
- }
- if(src->sign < 0) {
- dest->sign = -digits;
- }
- else {
- dest->sign = digits;
- }
- }
- if (abs(dest->sign) > dest->capacity) {
- CKRaise("extractbits overflow");
- }
-}
-
-#define NEW_MERSENNE 0
-
-/*
- * New gmersennemod, 24 Dec 1997. This runs significantly slower than the
- * original.
- */
-#if NEW_MERSENNE
-
-void
-gmersennemod(
- int n,
- giant g
-)
-/* g := g (mod 2^n - 1) */
-{
- int the_sign;
- giant scratch3 = borrowGiant(g->capacity);
- giant scratch4 = borrowGiant(1);
-
- if ((the_sign = gsign(g)) < 0) absg(g);
- while (bitlen(g) > n) {
- gtog(g,scratch3);
- gshiftright(n,scratch3);
- addg(scratch3,g);
- gshiftleft(n,scratch3);
- subg(scratch3,g);
- }
- if(isZero(g)) goto out;
- int_to_giant(1,scratch3);
- gshiftleft(n,scratch3);
- int_to_giant(1,scratch4);
- subg(scratch4,scratch3);
- if(gcompg(g,scratch3) >= 0) subg(scratch3,g);
- if (the_sign < 0) {
- g->sign = -g->sign;
- addg(scratch3,g);
- }
-out:
- returnGiant(scratch3);
- returnGiant(scratch4);
-}
-
-#else /* NEW_MERSENNE */
-
-void gmersennemod(int n, giant g) {
-/* g becomes g mod ((2^n)-1)
- 31 Jul 96 modified REC.
- 17 Jan 97 modified REC.
-*/
- unsigned bits = n & (GIANT_BITS_PER_DIGIT - 1);
- unsigned digits = 1 + ((n-1) >> GIANT_LOG2_BITS_PER_DIGIT);
- int isPositive = (g->sign > 0);
- int j;
- int b;
- int size;
- int foundzero;
- giantDigit mask = (bits == 0) ? GIANT_DIGIT_MASK : (giantDigit)((1<<bits)-1);
- giant scratch1;
-
- b = bitlen(g);
- if(b < n) {
- unsigned numDigits = (n + GIANT_BITS_PER_DIGIT - 1) >>
- GIANT_LOG2_BITS_PER_DIGIT;
- giantDigit lastWord = 0;
- giantDigit bits = 1;
-
- if(g->sign >= 0) return;
-
- /*
- * Cons up ((2**n)-1), add to g.
- */
- scratch1 = borrowGiant(numDigits + 1);
- scratch1->sign = numDigits;
- for(j=0; j<(int)(numDigits-1); j++) {
- scratch1->n[j] = GIANT_DIGIT_MASK;
- }
-
- /*
- * Last word has lower (n & (GIANT_BITS_PER_DIGIT-1)) bits set.
- */
- for(j=0; j < (int)(n & (GIANT_BITS_PER_DIGIT-1)); j++) {
- lastWord |= bits;
- bits <<= 1;
- }
- scratch1->n[numDigits-1] = lastWord;
- addg(g, scratch1); /* One version. */
- gtog(scratch1, g);
- returnGiant(scratch1);
- return;
- }
- if(b == n) {
- for(foundzero=0, j=0; j<b; j++) {
- if(bitval(g, j)==0) {
- foundzero = 1;
- break;
- }
- }
- if (!foundzero) {
- int_to_giant(0, g);
- return;
- }
- }
-
- absg(g);
- scratch1 = borrowGiant(g->capacity);
- while ( ((unsigned)(g->sign) > digits) ||
- ( ((unsigned)(g->sign)==digits) && (g->n[digits-1] > mask))) {
- extractbits(n, g, scratch1);
- gshiftright(n, g);
- addg(scratch1, g);
- }
- size = g->sign;
-
-/* Commence new negation routine - REC 17 Jan 1997. */
- if (!isPositive) { /* Mersenne negation is just bitwise complement. */
- for(j = digits-1; j >= size; j--) {
- g->n[j] = GIANT_DIGIT_MASK;
- }
- for(j = size-1; j >= 0; j--) {
- g->n[j] = ~g->n[j];
- }
- g->n[digits-1] &= mask;
- j = digits-1;
- while((g->n[j] == 0) && (j > 0)) {
- --j;
- }
- size = j+1;
- }
-/* End new negation routine. */
-
- g->sign = size;
- if (abs(g->sign) > g->capacity) {
- CKRaise("gmersennemod overflow");
- }
- if (size < (int)digits) {
- goto bye;
- }
- if (g->n[size-1] != mask) {
- goto bye;
- }
- mask = GIANT_DIGIT_MASK;
- for(j=0; j<(size-1); j++) {
- if (g->n[j] != mask) {
- goto bye;
- }
- }
- g->sign = 0;
- bye:
- returnGiant(scratch1);
-}
-
-#endif /* NEW_MERSENNE */
-
-void mulg(giant a, giant b) { /* b becomes a*b. */
-
- int i;
- int asize, bsize;
- giantDigit *bptr = b->n;
- giantDigit mult;
- giant scratch1;
- giantDigit carry;
- giantDigit *scrPtr;
-
-
- if (isZero(b)) {
- return;
- }
- if (isZero(a)) {
- gtog(a, b);
- return;
- }
- if(a == b) {
- grammarSquare(b);
- return;
- }
-
- bsize = abs(b->sign);
- asize = abs(a->sign);
- scratch1 = borrowGiant((asize+bsize));
- scrPtr = scratch1->n;
-
- for(i=0; i<asize+bsize; ++i) {
- scrPtr[i]=0;
- }
-
- for(i=0; i<bsize; ++i, scrPtr++) {
- mult = bptr[i];
- if (mult != 0) {
- carry = VectorMultiply(mult,
- a->n,
- asize,
- scrPtr);
- /* handle MSD carry */
- scrPtr[asize] += carry;
- }
- }
- bsize+=asize;
- if(scratch1->n[bsize - 1] == 0) {
- --bsize;
- }
- scratch1->sign = gsign(a) * gsign(b) * bsize;
- if (abs(scratch1->sign) > scratch1->capacity) {
- CKRaise("GiantGrammarMul overflow");
- }
- gtog(scratch1,b);
- returnGiant(scratch1);
-
- #if FEE_DEBUG
- (void)bitlen(b); // Assertion....
- #endif /* FEE_DEBUG */
- PROF_INCR(numMulg); // for normal profiling
- INCR_MULGS; // for ellipticMeasure
-}
-
-void grammarSquare(giant a) {
- /*
- * For now, we're going to match the old implementation line for
- * line by maintaining prod, carry, and temp as double precision
- * giantDigits. There is probably a much better implementation....
- */
- giantDigit prodLo;
- giantDigit prodHi;
- giantDigit carryLo = 0;
- giantDigit carryHi = 0;
- giantDigit tempLo;
- giantDigit tempHi;
- unsigned int cur_term;
- unsigned asize;
- unsigned max;
- giantDigit *ptr = a->n;
- giantDigit *ptr1;
- giantDigit *ptr2;
- giant scratch;
-
- /* dmitch 11 Jan 1998 - special case for a == 0 */
- if(a->sign == 0) {
- goto end;
- }
- /* end a == 0 case */
- asize = abs(a->sign);
- max = asize * 2 - 1;
- scratch = borrowGiant(2 * asize);
- asize--;
-
- /*
- * temp = *ptr;
- * temp *= temp;
- * scratch->n[0] = temp;
- * carry = temp >> 16;
- */
- giantMulDigits(*ptr, *ptr, &tempLo, &tempHi);
- scratch->n[0] = tempLo;
- carryLo = tempHi;
- carryHi = 0;
-
- for (cur_term = 1; cur_term < max; cur_term++) {
- ptr1 = ptr2 = ptr;
- if (cur_term <= asize) {
- ptr2 += cur_term;
- } else {
- ptr1 += cur_term - asize;
- ptr2 += asize;
- }
-
- /*
- * prod = carry & 0xFFFF;
- * carry >>= 16;
- */
- prodLo = carryLo;
- prodHi = 0;
- carryLo = carryHi;
- carryHi = 0;
- while(ptr1 < ptr2) {
- /*
- * temp = *ptr1++ * *ptr2--;
- */
- giantMulDigits(*ptr1++, *ptr2--, &tempLo, &tempHi);
-
- /*
- * prod += (temp << 1) & 0xFFFF;
- */
- giantAddDouble(&prodLo, &prodHi, (tempLo << 1));
-
- /*
- * carry += (temp >> 15);
- * use bits from both product digits..
- */
- giantAddDouble(&carryLo, &carryHi,
- (tempLo >> (GIANT_BITS_PER_DIGIT - 1)));
- giantAddDouble(&carryLo, &carryHi, (tempHi << 1));
-
- /* snag the msb from that last shift */
- carryHi += (tempHi >> (GIANT_BITS_PER_DIGIT - 1));
- }
- if (ptr1 == ptr2) {
- /*
- * temp = *ptr1;
- * temp *= temp;
- */
- giantMulDigits(*ptr1, *ptr1, &tempLo, &tempHi);
-
- /*
- * prod += temp & 0xFFFF;
- */
- giantAddDouble(&prodLo, &prodHi, tempLo);
-
- /*
- * carry += (temp >> 16);
- */
- giantAddDouble(&carryLo, &carryHi, tempHi);
- }
-
- /*
- * carry += prod >> 16;
- */
- giantAddDouble(&carryLo, &carryHi, prodHi);
-
- scratch->n[cur_term] = prodLo;
- }
- if (carryLo) {
- scratch->n[cur_term] = carryLo;
- scratch->sign = cur_term+1;
- } else scratch->sign = cur_term;
-
- gtog(scratch,a);
- returnGiant(scratch);
-end:
- PROF_INCR(numGsquare);
-}
-
-/*
- * Clear all of a giant's data fields, for secure erasure of sensitive data.,
- */
-void clearGiant(giant g)
-{
- unsigned i;
-
- for(i=0; i<g->capacity; i++) {
- g->n[i] = 0;
- }
- g->sign = 0;
-}
-
-#if ENGINE_127_BITS
-/*
- * only used by engineNSA127.c, which is obsolete as of 16 Jan 1997
- */
-int
-scompg(int n, giant g) {
- if((g->sign == 1) && (g->n[0] == n)) return(1);
- return(0);
-}
-
-#endif // ENGINE_127_BITS
-
-/*
- */
-
-/*
- * Calculate the reciprocal of a demonimator used in divg_via_recip() and
- * modg_via_recip().
- */
-void
-make_recip(giant d, giant r)
-/* r becomes the steady-state reciprocal
- 2^(2b)/d, where b = bit-length of d-1. */
-{
- int b;
- int giantSize = 4 * abs(d->sign);
- giant tmp = borrowGiant(giantSize);
- giant tmp2 = borrowGiant(giantSize);
-
- if (isZero(d) || (d->sign < 0))
- {
- CKRaise("illegal argument to make_recip");
- }
- int_to_giant(1, r); subg(r, d); b = bitlen(d); addg(r, d);
- gshiftleft(b, r); gtog(r, tmp2);
- while(1) {
- gtog(r, tmp);
- gsquare(tmp);
- gshiftright(b, tmp);
- mulg(d, tmp);
- gshiftright(b, tmp);
- addg(r, r); subg(tmp, r);
- if(gcompg(r, tmp2) <= 0) break;
- gtog(r, tmp2);
- }
- int_to_giant(1, tmp);
- gshiftleft(2*b, tmp);
- gtog(r, tmp2); mulg(d, tmp2);
- subg(tmp2, tmp);
- int_to_giant(1, tmp2);
- while(tmp->sign < 0) {
- subg(tmp2, r);
- addg(d, tmp);
- }
-
- returnGiant(tmp);
- returnGiant(tmp2);
- return;
-}
-
-/*
- * Optimized divg, when reciprocal of denominator is known.
- */
-void
-divg_via_recip(giant d, giant r, giant n)
-/* n := n/d, where r is the precalculated
- steady-state reciprocal of d. */
-{
- int s = 2*(bitlen(r)-1), sign = gsign(n);
- int giantSize = (4 * abs(d->sign)) + abs(n->sign);
- giant tmp = borrowGiant(giantSize);
- giant tmp2 = borrowGiant(giantSize);
-
- if (isZero(d) || (d->sign < 0))
- {
- CKRaise("illegal argument to divg_via_recip");
- }
- n->sign = abs(n->sign);
- int_to_giant(0, tmp2);
- while(1) {
- gtog(n, tmp);
- mulg(r, tmp);
- gshiftright(s, tmp);
- addg(tmp, tmp2);
- mulg(d, tmp);
- subg(tmp, n);
- if(gcompg(n,d) >= 0) {
- subg(d,n);
- iaddg(1, tmp2);
- }
- if(gcompg(n,d) < 0) break;
- }
- gtog(tmp2, n);
- n->sign *= sign;
- returnGiant(tmp);
- returnGiant(tmp2);
- return;
-}
-
-/*
- * Optimized modg, when reciprocal of denominator is known.
- */
-
-/* New version, 24 Dec 1997. */
-
-void
-modg_via_recip(
- giant d,
- giant r,
- giant n
-)
-/* This is the fastest mod of the present collection.
- n := n % d, where r is the precalculated
- steady-state reciprocal of d. */
-
-{
- int s = (bitlen(r)-1), sign = n->sign;
- int giantSize = (4 * abs(d->sign)) + abs(n->sign);
- giant tmp, tmp2;
-
- tmp = borrowGiant(giantSize);
- tmp2 = borrowGiant(giantSize);
- if (isZero(d) || (d->sign < 0))
- {
- CKRaise("illegal argument to modg_via_recip");
- }
- n->sign = abs(n->sign);
- while (1)
- {
- gtog(n, tmp);
- /* bug fix 13 Apr 1998 */
- if(s == 0) {
- gshiftleft(1, tmp);
- }
- else {
- gshiftright(s-1, tmp);
- }
- /* end fix */
- mulg(r, tmp);
- gshiftright(s+1, tmp);
- mulg(d, tmp);
- subg(tmp, n);
- if (gcompg(n,d) >= 0)
- subg(d,n);
- if (gcompg(n,d) < 0)
- break;
- }
- if (sign >= 0)
- goto done;
- if (isZero(n))
- goto done;
- negg(n);
- addg(d,n);
-done:
- returnGiant(tmp);
- returnGiant(tmp2);
- return;
-}
-
-/*
- * Unoptimized, inefficient general modg, when reciprocal of denominator
- * is not known.
- */
-void
-modg(
- giant d,
- giant n
-)
-{
- /* n becomes n%d. n is arbitrary, but the denominator d must be
- * positive! */
-
- /*
- * 4/9/2001: seeing overflow on this recip. Alloc per
- * d->capacity, not d->sign.
- */
- //giant recip = borrowGiant(2 * abs(d->sign));
- giant recip = borrowGiant(2 * d->capacity);
-
- #if WARN_UNOPTIMIZE
- dbgLog(("Warning: unoptimized modg!\n"));
- #endif // WARN_UNOPTIMIZE
-
- make_recip(d, recip);
- modg_via_recip(d, recip, n);
- returnGiant(recip);
-}
-
-/*
- * Unoptimized, inefficient general divg, when reciprocal of denominator
- * is not known.
- */
-void
-divg(
- giant d,
- giant n
-)
-{
- /* n becomes n/d. n is arbitrary, but the denominator d must be
- * positive!
- */
-
- giant recip = borrowGiant(2 * abs(d->sign));
-
- #if WARN_UNOPTIMIZE
- dbgLog(("Warning: unoptimized divg!\n"));
- #endif // WARN_UNOPTIMIZE
-
- make_recip(d, recip);
- divg_via_recip(d, recip, n);
- returnGiant(recip);
-}
projects / apple / security.git / blobdiff