[apple/security.git] / libsecurity_cryptkit / lib / giantIntegers.c

/* Copyright (c) 1998 Apple Computer, 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 COMPUTER, INC. AND THE
 * ORIGINAL LICENSEE THAT OBTAINED THESE MATERIALS FROM APPLE COMPUTER,
 * 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
 ----------------
 12/04/98       dmitch/R. Crandall
        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      Doug Mitchell at Apple
        Major rewrite of core arithmetic routines to make this module
                independent of size of giantDigit.
        Removed idivg() and radixdiv().
 20 Jan 98      Doug Mitchell at Apple
        Deleted FFT arithmetic; simplified mulg().
 09 Jan 98      Doug Mitchell and Richard Crandall at Apple
        gshiftright() optimization.
 08 Jan 98      Doug Mitchell at Apple
        newGiant() returns NULL on malloc failure
 24 Dec 97      Doug Mitchell and Richard Crandall at Apple
        New grammarSquare(); optimized modg_via_recip()
 11 Jun 97      Doug Mitchell and Richard Crandall 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      Doug Mitchell at Apple
        Changed size of giantstruct.n to 1 for Mac build
 05 Feb 97      Doug Mitchell at Apple
        newGiant() no longer modifies CurrentMaxShorts or giant stack
        Added modg profiling
 01 Feb 97      Doug Mitchell at NeXT
        Added iszero() check in gcompg
 17 Jan 97      Richard Crandall, Doug Mitchell at NeXT
        Fixed negation bug in gmersennemod()
        Fixed n[words-1] == 0 bug in extractbits()
        Cleaned up lots of static declarations
 19 Sep 96      Doug Mitchell at NeXT
        Fixed --size underflow bug in normal_subg().
  4 Sep 96      Doug Mitchell at NeXT
        Fixed (b<n), (sign<0) case in gmersennemod() to allow for arbitrary n.
  9 Aug 96      Doug Mitchell at NeXT
        Fixed sign-extend bug in data_to_giant().
  7 Aug 96      Doug Mitchell at NeXT
        Changed precision in newtondivide().
        Removed #ifdef UNUSED code.
 24 Jul 96      Doug Mitchell at NeXT
        Added scompg().
      1991      Richard Crandall at NeXT
        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     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 */

/******** 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 <CoreServices/../Frameworks/CarbonCore.framework/Headers/MacTypes.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()
{
        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);
}

/* new addg, negg, and gcompg from Crandall 6/95 */
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

/*
 * New modg() and divg() arithmetic, from R. Crandall, 11 June 1997.
 */

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