vfactor source/vgiants.c

/**************************************************************
 *
 *  vgiants.c
 *
 *  AltiVec-based giant integer package
 *  Derived from original giants.c project at 
 *  Next Software, Inc.
 *
 *  This package is part of ongoing research in the
 *  Advanced Computation Group, Apple Computer.
 *  
 *  c. 1999 Apple Computer, Inc.
 *  All Rights Reserved.
 *
 *
 *************************************************************/
 
/*
    Disclaimer: IMPORTANT:  This Apple software is supplied to you by Apple Computer, Inc.
                ("Apple") in consideration of your agreement to the following terms, and your
                use, installation, modification or redistribution of this Apple software
                constitutes acceptance of these terms.  If you do not agree with these terms,
                please do not use, install, modify or redistribute this Apple software.
 
                In consideration of your agreement to abide by the following terms, and subject
                to these terms, Apple grants you a personal, non-exclusive license, under AppleÕs
                copyrights in this original Apple software (the "Apple Software"), to use,
                reproduce, modify and redistribute the Apple Software, with or without
                modifications, in source and/or binary forms; provided that if you redistribute
                the Apple Software in its entirety and without modifications, you must retain
                this notice and the following text and disclaimers in all such redistributions of
                the Apple Software.  Neither the name, trademarks, service marks or logos of
                Apple Computer, Inc. may be used to endorse or promote products derived from the
                Apple Software without specific prior written permission from Apple.  Except as
                expressly stated in this notice, no other rights or licenses, express or implied,
                are granted by Apple herein, including but not limited to any patent rights that
                may be infringed by your derivative works or by other works in which the Apple
                Software may be incorporated.
 
                The Apple Software is provided by Apple on an "AS IS" basis.  APPLE MAKES NO
                WARRANTIES, EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION THE IMPLIED
                WARRANTIES OF NON-INFRINGEMENT, MERCHANTABILITY AND FITNESS FOR A PARTICULAR
                PURPOSE, REGARDING THE APPLE SOFTWARE OR ITS USE AND OPERATION ALONE OR IN
                COMBINATION WITH YOUR PRODUCTS.
 
                IN NO EVENT SHALL APPLE BE LIABLE FOR ANY SPECIAL, INDIRECT, INCIDENTAL OR
                CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
                GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
                ARISING IN ANY WAY OUT OF THE USE, REPRODUCTION, MODIFICATION AND/OR DISTRIBUTION
                OF THE APPLE SOFTWARE, HOWEVER CAUSED AND WHETHER UNDER THEORY OF CONTRACT, TORT
                (INCLUDING NEGLIGENCE), STRICT LIABILITY OR OTHERWISE, EVEN IF APPLE HAS BEEN
                ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
 
#include "giantsdebug.h"
#if __MWERKS__
#include <altivec.h>
#endif
#include "vgiants.h"
#include "vecarith.h"
#include <Memory.h>
#include <Quickdraw.h>      // for random
#include <stddef.h>
#include <stdlib.h>
#include <math.h>
#include "giantsstack.h"
#include <segload.h>        // for ExitToShell
 
//////////////////////////////////////////////
// globals
//////////////////////////////////////////////
 
static giant    cur_recip = NULL;
static giant    cur_den = NULL;
 
 
 
//////////////////////////////////////////////
 
#if 0
#define BORROW_GIANT()  newgiantbits(DEFAULT_GIANT_BITS)
#define RETURN_GIANT(a) disposegiant(a)
#else
#define BORROW_GIANT()  popg()
#define RETURN_GIANT(a) pushg(1)
#endif
 
#define ADDVECS addvecs
#define SUBVECS subvecs
 
#define max(a, b) (a > b ? a : b)
 
#if GDEBUG 
 
void ASSERT_GVALID(giant g)
{
    GAssert((g->giantSign == -1) || (g->giantSign == 1));
    
    if (((long)g->vectors & 0x0f) != 0) {
        DebugStr("\pvectors not 16-byte aligned.");
    }
    
    if (g->giantDigits > g->giantCapacity) {
        DebugStr("\pgiant digits overflow!");       
    }
    
    if (g->giantDigits != 0) {
        vector unsigned long    veczero = (vector unsigned long)(0);
        
        if (vec_all_eq(g->vectors[g->giantDigits-1], veczero)) {
            DebugStr("\pMSV is zero");
        }
    }
}
#endif
 
static void FATAL_ERROR(long err)
{
#pragma unused(err)
    DebugStr("\pfatal error in vgiants");
    //ExitToShell();
}
 
 
 
 
 
gmatrix
newgmatrix(
    void
)
/* Allocates space for a gmatrix struct, but does not
 * allocate space for the giants. */
{
    return((gmatrix) malloc (4*sizeof(giant)));
}
 
 
GiantStructPtr  newgiantbits(unsigned long capacityBits)
{
    unsigned long   giantCapacity;
    unsigned long   allocSize;
    GiantStructPtr  pg;
    
    if (!capacityBits) capacityBits = 1;
    
    giantCapacity = (capacityBits + BITS_PER_VECTOR-1) / BITS_PER_VECTOR;
    allocSize = giantCapacity * BYTES_PER_VECTOR;
    pg = (GiantStructPtr)NewPtr(sizeof(GiantStruct));
 
    GAssert(pg != nil); 
    
    pg->giantSign = 1;
    pg->giantDigits = 0;
    pg->giantCapacity = giantCapacity;  
    
    pg->vectors = (vector unsigned long *)NewPtr(allocSize);
 
    GAssert(pg->vectors != nil);
    
    return pg;
}
 
GiantStructPtr  newgiant(unsigned long capacityShorts)
{
    return newgiantbits(capacityShorts << 4);
}
 
 
void disposegiant(GiantStructPtr g)
{
    ASSERT_GVALID(g);
    
    DisposePtr((Ptr)g->vectors);    
    DisposePtr((Ptr)g);
}
 
static void reallocgdigits(giant g, unsigned long vectorcount, Boolean save)
{
    long                i;
            
    vector unsigned long *pNewVectors = (vector unsigned long *)NewPtr( vectorcount * BYTES_PER_VECTOR );
 
    GAssert(vectorcount >= g->giantDigits);
 
    GAssert(pNewVectors != nil);
    
    if (pNewVectors != nil) {
        if (save) {
            for (i=0; i<g->giantDigits; i++) {
                pNewVectors[i] = g->vectors[i];
            }
        }
                
        DisposePtr((Ptr)g->vectors);
        g->vectors = pNewVectors;
        g->giantCapacity = vectorcount;
    }
}
 
 
 
static void reallocg(giant g, unsigned long newbits, Boolean save)
{
    unsigned long   vectorcount =  (newbits + 127) / 128;
    
    reallocgdigits(g, vectorcount, save);
}
 
static void EnsureDigits(giant g, unsigned long newdigits, Boolean save)
{
    if (g->giantCapacity < newdigits) {
        reallocgdigits(g, newdigits, save);
    }
}
 
static void EnsureBits(giant g, unsigned long newbits, Boolean save)
{
    if (g->giantCapacity * BITS_PER_VECTOR < newbits) {
        reallocg(g, newbits, save);
    }
}
 
static void EnsureBytes(giant g, unsigned long newbytes, Boolean save)
{
    EnsureBits(g, newbytes << 3, save);
}
 
static void AddMSDLong(giant g, unsigned long digit)
{
    unsigned long *topvector;
 
    if (g->giantCapacity == g->giantDigits) {
        reallocgdigits(g, g->giantDigits+1, true);
    }   
    
    topvector = (unsigned long*)&g->vectors[g->giantDigits];
    
    *topvector++ = 0;
    *topvector++ = 0;
    *topvector++ = 0;
    *topvector = digit;
    
    g->giantDigits++;
}
 
#if 1
static void justg(giant g)
{
    if (g->giantDigits) {
        unsigned long       counter = g->giantDigits;
        unsigned long       *pLongs = (unsigned long*)&g->vectors[g->giantDigits];  // point to one past last
 
        while (counter--) {
            if (*--pLongs) break;
            if (*--pLongs) break;
            if (*--pLongs) break;
            if (*--pLongs) break;
            
            g->giantDigits--;
        }
    
    }
 
}
#else
static void justg(giant g)
{
    if (g->giantDigits) {
        vector unsigned long vecZero = (vector unsigned long)(0);
        
        while (vec_all_eq(vecZero, g->vectors[g->giantDigits-1]) && g->giantDigits) {
            // do nothing
            g->giantDigits--;
        }
    }
}
#endif
 
static unsigned long gtoul(giant g)
{
    GAssert(bitlen(g) <= 32);
    
    if (!g->giantSign) return 0;
    
    return ((unsigned long *)g->vectors)[LONGS_PER_VECTOR-1];
}
 
 
 
 
 
int             cur_run = 0;
double *        sinCos=NULL;
 
void
init_sinCos(
    int         n
)
{
    int         j;
    double  e = TWOPI/n;
 
    if (n<=cur_run)
        return;
    cur_run = n;
    if (sinCos)
        free(sinCos);
    sinCos = (double *)malloc(sizeof(double)*(1+(n>>2)));
    for (j=0;j<=(n>>2);j++)
    {
        sinCos[j] = sin(e*j);
    }
}
 
 
double
s_sin(
    int     n
)
{
    int     seg = n/(cur_run>>2);
 
    switch (seg)
    {
        case 0: return(sinCos[n]);
        case 1: return(sinCos[(cur_run>>1)-n]);
        case 2: return(-sinCos[n-(cur_run>>1)]);
        case 3: return(-sinCos[cur_run-n]);
    }
    return 0;
}
 
 
double
s_cos(
    int     n
)
{
    int     quart = (cur_run>>2);
 
    if (n < quart)
        return(s_sin(n+quart));
    return(-s_sin(n-quart));
}
 
 
 
 
unsigned char GetNthByte(giant g, unsigned long pos)
{
    unsigned long   byteindex;
    
    //GAssert(pos < g->giantDigits * BYTES_PER_VECTOR);
    
    if (pos >= g->giantDigits * BYTES_PER_VECTOR) {
        return 0;
    }
    
    byteindex = pos & ~0xf;
    
    byteindex += BYTES_PER_VECTOR - (pos & 0xf)-1;
    
    return ((unsigned char*)g->vectors)[byteindex];
 
}
 
static unsigned short * GetNthShortPtr(giant g, unsigned long pos)
{
    unsigned long   shortindex;
        
    if (pos >= g->giantDigits * SHORTS_PER_VECTOR) {
        return 0;
    }
    
    shortindex = pos & ~0x7;
    
    shortindex += SHORTS_PER_VECTOR - (pos & 0x7)-1;
    
    return &((unsigned short*)g->vectors)[shortindex];
 
}
 
 
static unsigned long * GetNthLongPtr(giant g, unsigned long pos)
{
    unsigned long   longindex;
        
    if (pos >= g->giantDigits * LONGS_PER_VECTOR) {
        return 0;
    }
    
    longindex = pos & ~0x3;
    
    longindex += LONGS_PER_VECTOR - (pos & 0x3)-1;
    
    return &((unsigned long*)g->vectors)[longindex];
 
}
 
static unsigned long NumLongs(giant g)
{
    unsigned long result = 0;
    unsigned long *longPtr;
    unsigned long i;
    
    if (g->giantDigits) {
        result = g->giantDigits * LONGS_PER_VECTOR;
        
        longPtr = (unsigned long*)&g->vectors[g->giantDigits-1];
        
        for (i=0; i<LONGS_PER_VECTOR; i++) {
            if (*longPtr++ != 0)  break;
            result--;   
        }
    }
    
    return result;
}
 
 
static unsigned long NumBytes(giant g)
{
    unsigned long result = 0;
    unsigned char *bytePtr;
    unsigned long i;
    
    if (g->giantDigits) {
        result = g->giantDigits * BYTES_PER_VECTOR;
        
        bytePtr = (unsigned char*)&g->vectors[g->giantDigits-1];
        
        for (i=0; i<BYTES_PER_VECTOR; i++) {
            if (*bytePtr++ != 0)  break;
            result--;   
        }
    }
    
    return result;
}
 
static unsigned long BytesToDigits(unsigned long bytes)
{
    return (bytes+(BYTES_PER_VECTOR-1))/BYTES_PER_VECTOR;
}
 
 
 
 
static Boolean gshiftleftsubvector(unsigned long bits, giant g)
{
    vector unsigned char    temp;   
    vector unsigned char    shiftFactor, negShiftFactor;
    vector unsigned char    vecZero;
    vector unsigned long    t1, t2;
    long                        i;
    Boolean                 shiftedOut = false;
 
    ASSERT_GVALID(g);
            
    if (bits) {
        vecZero = (vector unsigned char)(0);
 
        *((unsigned char*)&temp) = bits;    
        
        shiftFactor = vec_splat(temp, 0);               // fill vector with shift amount
        negShiftFactor = vec_sub(vecZero, shiftFactor);
        
        // figure out if we need to expand
        t2 = vec_sro(g->vectors[g->giantDigits-1], negShiftFactor);
        t2 = vec_srl(t2, negShiftFactor);
        if (!vec_all_eq((vector unsigned long)vecZero, t2)) {
            shiftedOut = true;
            g->vectors[g->giantDigits] = t2;
        }
 
        for (i=g->giantDigits-1; i>0; i--) {
            t1 = vec_slo(g->vectors[i], shiftFactor);
            t1 = vec_sll(t1, shiftFactor);
            t2 = vec_sro(g->vectors[i-1], negShiftFactor);
            t2 = vec_srl(t2, negShiftFactor);
            g->vectors[i] = vec_or(t1, t2);
        }
        
        t1 = vec_slo(g->vectors[i], shiftFactor);
        g->vectors[i] = vec_sll(t1, shiftFactor);
    }
        
    return shiftedOut;
}
 
void        gshiftleft(unsigned long bits, giant g)
{
    unsigned long   bitlength;
    long            multiples;
 
    if (!bits) return;
    if (!g->giantDigits) return;
 
    bitlength = bitlen(g);
 
    if (bitlength + bits > (g->giantCapacity*BITS_PER_VECTOR)) {
        reallocg(g, bitlength + bits, true);
    }
    
    multiples = bits >> 7;
    
    if (multiples) {
        vector unsigned long    *pSrc;
        vector unsigned long    *pDst;
        unsigned long           digitCount = g->giantDigits;
        
        pSrc = &g->vectors[digitCount-1];
        pDst = &g->vectors[digitCount+multiples-1];
 
        g->giantDigits += multiples;
        
        // move up digits
        while (digitCount--) {
            *pDst-- = *pSrc--;
        }
        
        // fill remaining shifting with zeros
        while (multiples--) {
            *pDst-- = (vector unsigned long)(0);
        }
    }
    
    if (gshiftleftsubvector(bits & 0x7f, g)) {
        g->giantDigits++;
    }
    
    ASSERT_GVALID(g);
}
 
 
static void gshiftrightsubvector(unsigned long bits, giant g)
{
    vector unsigned char    temp;   
    vector unsigned char    shiftFactor, negShiftFactor;
    vector unsigned char    vecZero;
    vector unsigned long    t1, t2;
    long                        i;
    
    if (bits) {
        vecZero = (vector unsigned char)(0);
        
        // set 0th element of shiftFactor to be bits
        *((unsigned char*)&temp) = bits;
            
        shiftFactor = vec_splat(temp, 0);   // fill vector with shift amount
        negShiftFactor = vec_sub(vecZero, shiftFactor);
        
        for (i=0; i<g->giantDigits-1; i++) {
            t1 = vec_sro(g->vectors[i], shiftFactor);
            t1 = vec_srl(t1, shiftFactor);
            t2 = vec_slo(g->vectors[i+1], negShiftFactor);
            t2 = vec_sll(t2, negShiftFactor);
            g->vectors[i] = vec_or(t1, t2);
        }
        
        t1 = vec_sro(g->vectors[i], shiftFactor);
        g->vectors[i] = vec_srl(t1, shiftFactor);
    }
}
 
void        gshiftright(unsigned long bits, giant g)
{
    long        multiples;
 
    if (!bits) return;
    if (!g->giantDigits) return;
 
    multiples = bits >> 7;
    
    if (multiples) {
        if (multiples >= g->giantDigits) {
            g->giantDigits = 0; 
        } else {        
            vector unsigned long    *pSrc;
            vector unsigned long    *pDst;
            unsigned long           newdigits = g->giantDigits-multiples;
            
            g->giantDigits = newdigits;
            
            pSrc = &g->vectors[multiples];
            pDst = g->vectors;
            
            while (newdigits--) {
                *pDst++ = *pSrc++;
            }
        }
    }
 
        
    if (g->giantDigits) {
        gshiftrightsubvector(bits& 0x7f, g);
    }
    
    justg(g);
    
    ASSERT_GVALID(g);   
}
 
 
 
 
#if 1
unsigned long           bitlen(giant g)
{
    unsigned long           length;
    unsigned long           *longPtr;
    unsigned long           i;
    register unsigned long  longdigit;
    ASSERT_GVALID(g);
 
    if (!g->giantDigits) return 0;
    
    length = length = BITS_PER_VECTOR * g->giantDigits;
 
    longPtr = (unsigned long*)&g->vectors[g->giantDigits-1];
    
    for (i=0; i<4; i++) {
        longdigit = longPtr[i];
 
        if (!longdigit) {
            length -= 32;   // bits per long
        } else {
#if __MWERKS__
            register unsigned long leadingzeros=0;  // initialized to eliminate compiler warning
            asm {
                cntlzw  leadingzeros, longdigit
            }
            length -= leadingzeros;         
#else               
            while (!(longdigit & 0x80000000)) {
                longdigit = longdigit << 1;
                length -= 1;
            }
 
#endif      
            break;  
        }
    }
    
    return length;
 
}
#else
unsigned long           bitlen(giant g)
{
    unsigned long           length;
    unsigned long           count;
#if __MWERKS__
    vector unsigned long    vecones = (vector signed long)(-1);
#else   
    vector unsigned long    vecones = (vector unsigned long)(-1);
#endif  
    vector unsigned long    vecmask = (vector unsigned long)(0);
    vector unsigned long    veczero = (vector unsigned long)(0);
    vector unsigned long    digit;
    unsigned long           longindex;
    unsigned long           thelong;
    
    if (!g->giantDigits) return 0;
    
    count = 3;
        
    digit = g->vectors[g->giantDigits-1];
    
    longindex = 0;
 
    length = BITS_PER_VECTOR * g->giantDigits;
 
    vecmask = vec_sld(vecones, vecmask, 12);
    
    while (count-- && vec_all_eq(veczero, vec_and(vecmask, digit))) {
        vecmask = vec_sld(vecones, vecmask, 12);
        length -= 32;
        longindex++;
    }   
        
    thelong = ((unsigned long*)(&g->vectors[g->giantDigits-1]))[longindex];
    
    while (!(thelong & 0x80000000)) {
        thelong = thelong << 1;
        length -= 1;
    }
        
    return length;  
}
#endif
 
long    bitval(giant g, long pos)
{
#if 1
    return (GetNthByte(g, pos >> 3) & (1 << (pos & 0x07)));
#else
    unsigned char   mask;
    unsigned long   byteindex;
    
    byteindex = (pos >> 7) << 4;
    
    byteindex += BYTES_PER_VECTOR - ((pos & 0x7f) >> 3)-1;
    
    mask = 1 << (pos & 0x07);
    
    return ((unsigned char*)g->vectors)[byteindex] & mask;
#endif
}
 
 
 
void gtog(giant a, giant b)
{
    long i;
 
    ASSERT_GVALID(a);
    ASSERT_GVALID(b);
 
    if (b->giantCapacity < a->giantDigits) {
        reallocgdigits(b, a->giantDigits, false);
    }
 
    GAssert(b->giantCapacity >= a->giantDigits)
    
    
    for (i=0; i<a->giantDigits; i++) {
        b->vectors[i] = a->vectors[i];
    }
    
    b->giantDigits = a->giantDigits;
    b->giantSign = a->giantSign;
}
 
 
long CompareMagnitude(giant a, giant b)
{
    long                vectorindex, longindex;
    unsigned long       *aPtr, *bPtr;
 
    ASSERT_GVALID(a);
    ASSERT_GVALID(b);
    
    if (a->giantDigits > b->giantDigits) return 1;
    if (a->giantDigits < b->giantDigits) return -1;
 
 
    for (vectorindex=a->giantDigits-1; vectorindex >= 0; vectorindex--) {
        aPtr = (unsigned long*)&a->vectors[vectorindex];
        bPtr = (unsigned long*)&b->vectors[vectorindex];
        
        for (longindex = 0; longindex< LONGS_PER_VECTOR; longindex++) {
            if (*aPtr > *bPtr) return 1;
            if (*aPtr < *bPtr) return -1;
            aPtr++;
            bPtr++;     
        }
    }
        
    return 0;
}
 
 
 
static void subg_normal(GiantStructPtr a, GiantStructPtr b, GiantStructPtr result)
{
    unsigned long endDigits = b->giantDigits;
    
    EnsureDigits(result, b->giantDigits, true);
        
    SUBVECS(a->vectors, a->giantDigits, b->vectors, b->giantDigits, result->vectors);
 
    result->giantDigits = endDigits;
            
    justg(result);
}
 
 
void addg(giant a, giant b)
{
    ASSERT_GVALID(a);
    ASSERT_GVALID(b);
 
    if (isZero(a)) {
        return;
    }
    
    if (isZero(b)) {
        gtog(a, b);
        return;
    }
 
    if (a->giantSign == b->giantSign) {
        unsigned long   newdigits;
        
        newdigits = max(a->giantDigits, b->giantDigits)+1;
        EnsureDigits(b, newdigits, true);
        if (!ADDVECS(a->vectors, a->giantDigits, b->vectors, b->giantDigits, b->vectors)) {
            newdigits--;
        }
        b->giantDigits = newdigits;
 
    } else {
        if (b->giantSign > 0) {
            // a is negative, b is positive
            if (CompareMagnitude(b, a) >= 0) {
                subg_normal(a, b, b);
            } else {
                // a neg, b pos, b<a
                subg_normal(b, a, b);
                negg(b);
            }
        } else {
            // a is positive, b is negative
            if (CompareMagnitude(b, a) < 0) {
                subg_normal(b, a, b);
                negg(b);
            } else {
                // a neg, b pos, b>a
                subg_normal(a, b, b);
            }           
        }
    }   
    ASSERT_GVALID(a);
    ASSERT_GVALID(b);
}
 
 
 
// result = b-a
void subg(giant a, giant b)
{
    ASSERT_GVALID(a);
    ASSERT_GVALID(b);
 
    if (isZero(a)) {
        return;
    }
    
    if (isZero(b)) {
        gtog(a, b);
        negg(b);
        return;
    }
 
    if (a->giantSign != b->giantSign) {
        unsigned long   newdigits;
        
        newdigits = max(a->giantDigits, b->giantDigits)+1;
        EnsureDigits(b, newdigits, true);
        if (!ADDVECS(a->vectors, a->giantDigits, b->vectors, b->giantDigits, b->vectors)) {
            newdigits--;
        }
        b->giantDigits = newdigits;
    } else {
        if (b->giantSign > 0) {
            // a is negative, b is positive
            if (CompareMagnitude(b, a) >= 0) {
                subg_normal(a, b, b);
            } else {
                // a neg, b pos, b<a
                subg_normal(b, a, b);
                negg(b);
            }
        } else {
            // a is positive, b is negative
            if (CompareMagnitude(b, a) < 0) {
                subg_normal(b, a, b);
                negg(b);
            } else {
                // a neg, b pos, b>a
                subg_normal(a, b, b);
            }           
        }
    }
    
    ASSERT_GVALID(a);
    ASSERT_GVALID(b);
}
 
 
 
 
 
void grammarmulg(giant a, giant b)
{
    unsigned long   resultdigits = a->giantDigits + b->giantDigits;
    
    if (isZero(b)) return;
 
    if (isZero(a)) {
        gtog(a, b);
        return;
    }
 
    if (a->giantDigits == 1) {
        if (NumLongs(a) == 1) {
            
            smulg(gtoul(a), b);
            b->giantSign *= a->giantSign;
        
        } else {
            EnsureDigits(b, b->giantDigits+1, true);
            
            VecMult1(   a->vectors,
                        b->giantDigits,
                        b->vectors,
                        b->vectors);
                        
            b->giantDigits++;
            
            b->giantSign *= a->giantSign;
            
            justg(b);
        }
    } else  if (b->giantDigits == 1) {
        if (NumLongs(b) == 1) {
            unsigned long   blong = gtoul(b);
            long            bsign = b->giantSign;
            
            gtog(a, b);
            smulg(blong, b);
            b->giantSign *= bsign;
 
        } else {
            EnsureDigits(b, a->giantDigits+1, true);
            
            VecMult1(   b->vectors,
                        a->giantDigits,
                        a->vectors,
                        b->vectors);
                        
            b->giantDigits = a->giantDigits+1;
            
            b->giantSign *= a->giantSign;
            
            justg(b);
        }
    } else {
        GiantStructPtr  temp = BORROW_GIANT();
        
        if (a->giantDigits == 2 && b->giantDigits == 2) {
            EnsureDigits(temp, 4, false);   
 
            mult2by2(   &a->vectors[1],
                        &a->vectors[0],
                        &b->vectors[1],
                        &b->vectors[0],
                        temp->vectors);
                    
            temp->giantDigits = 4;
            
            temp->giantSign = a->giantSign*b->giantSign;
 
            justg(temp);
            
            gtog(temp, b);
 
            RETURN_GIANT(temp);
        
        } else {    
            EnsureDigits(temp, resultdigits, false);    
 
            VecMult(a->vectors,
                    a->giantDigits,
                    b->vectors,
                    b->giantDigits,
                    temp->vectors);
                    
            temp->giantDigits = resultdigits;
            
            temp->giantSign = a->giantSign*b->giantSign;
            justg(temp);
            
            gtog(temp, b);
 
            RETURN_GIANT(temp);
        }
    }
    
    ASSERT_GVALID(b);
    
}
 
void
karatmulg(
        giant   x,
        giant   y)
/* y becomes x*y. */
{
        int     s = x->giantDigits, t = y->giantDigits, w, bits,
                sg = gsign(x)*gsign(y);
        giant   a, b, c, e, f;
 
        if((s <= KARAT_BREAK) || (t <= KARAT_BREAK)) {
                grammarmulg(x,y);
                return;
        }
        w = (s + t + 2)/4; bits = 16*w;
        
        a = BORROW_GIANT(); b = BORROW_GIANT(); c = BORROW_GIANT(); e = BORROW_GIANT(); f = BORROW_GIANT();
        gtog(x,a); if (w <= s) {a->giantDigits = w; justg(a);}
        absg(a);
        gtog(x,b); absg(b);
        gshiftright(bits, b);
 
        gtog(y,c); if (w <= t) {c->giantDigits = w; justg(c);}
        absg(c);
        absg(y);
        gshiftright(bits,y);
        gtog(a,e); addg(b,e);
        gtog(c,f); addg(y,f);
        karatmulg(e,f);  /* f := (a + b)(c + d). */
        karatmulg(c,a); /* a := a c. */
        karatmulg(b,y); /* d := b d. */
        subg(a,f); subg(y,f);
        gshiftleft(bits, y);
        addg(f, y);
        gshiftleft(bits, y);
 
        addg(a, y);
        setsign(y, sg);
        
        RETURN_GIANT(f);
        RETURN_GIANT(e);
        RETURN_GIANT(c);
        RETURN_GIANT(b);
        RETURN_GIANT(a);
 
        return;
}
 
void mulg(giant a, giant b)
{
    if((a->giantDigits <= KARAT_BREAK) || (b->giantDigits <= KARAT_BREAK)) {
        grammarmulg(a,b);
    } else {
        karatmulg(a, b);
    }
}
 
 
#if 1
/* g += i, with i non-negative. */
void        iaddg(unsigned long i,giant g)
{
    if (isZero(g)) {
        itog(i, g);
    } else {    
        ASSERT_GVALID(g);
 
        if (AddULongToVecs(g->vectors, g->giantDigits, i)) {
            AddMSDLong(g, 1);
        }
        
        justg(g);
 
        ASSERT_GVALID(g);
    }
}
#else
/* g += i, with i non-negative. */
void        iaddg(unsigned long i,giant g)
{
    giant       temp=BORROW_GIANT();
    
    ASSERT_GVALID(g);
 
    itog(i, temp);
    addg(temp, g);
 
    ASSERT_GVALID(g);
    
    RETURN_GIANT(temp);
}
#endif
 
 
 
/* Integer <-> giant. */
void itog(long n, giant g)
{
    unsigned long   absn = abs(n);
    
    if (absn) {
        unsigned long   *pLong = (unsigned long *)g->vectors;
        
        *pLong++ = 0;
        *pLong++ = 0;
        *pLong++ = 0;
        *pLong = absn;
 
        g->giantSign = (n < 0 ? -1 : 1);
        g->giantDigits = 1;
        
    } else {
        g->giantDigits = 0;
    }
        
    ASSERT_GVALID(g);   
}
 
 
#if 1
void        squareg(giant g)
{
    GiantStructPtr  temp;
    unsigned long   digits;
 
    if (isZero(g)) return;
 
    temp = BORROW_GIANT();
        
    digits = g->giantDigits * 2;
    
    EnsureDigits(temp, digits, false);  
        
    VecSquare(g->vectors,
            g->giantDigits,
            temp->vectors);
            
    temp->giantDigits = digits;
    
    temp->giantSign = 1;
    
    justg(temp);
    
    gtog(temp, g);
 
    RETURN_GIANT(temp);
    
    ASSERT_GVALID(g);
}   
 
#else
/* g *= g. */
void        squareg(giant g)
{
    ASSERT_GVALID(g);
    mulg(g, g);
 
    ASSERT_GVALID(g);
 
}
#endif
 
void        randomg(giant g, unsigned long numvecs)
{
    long i;
    
    if (numvecs > g->giantCapacity) numvecs = g->giantCapacity;
    
    for (i=0; i<numvecs*SHORTS_PER_VECTOR; i++) {
        ((unsigned short*)g->vectors)[i] = Random();
    }
    
    g->giantSign = (Random() & 1) ? -1 : 1;
    g->giantDigits = numvecs;
    
    justg(g);
 
    ASSERT_GVALID(g);
}
 
 
 
 
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 = gsign(n);
    giant   tmp, tmp2;
 
    if (isZero(d) || (gsign(d) < 0))
    {
        FATAL_ERROR(SIGN);
    }
 
    tmp = BORROW_GIANT();
    tmp2 = BORROW_GIANT();
    
    if (isone(d)) {
        itog(0, n);
        goto done;
    }
    
    absg(n);
    while (1) 
    {
        gtog(n, tmp); gshiftright(s-1, tmp);    
        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:
    RETURN_GIANT(tmp);
    RETURN_GIANT(tmp2);
    return;
}
 
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;
    giant   tmp, tmp2;
 
    if (isZero(d) || (gsign(d) < 0))
    {
        FATAL_ERROR(SIGN);
    }
    tmp = BORROW_GIANT();
    tmp2 = BORROW_GIANT();
    itog(1, r); 
    subg(r, d); 
    b = bitlen(d); 
    addg(r, d);
    gshiftleft(b, r); 
    gtog(r, tmp2);
    while (1) 
    {
        gtog(r, tmp);
        squareg(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);
    }
    itog(1, tmp);
    gshiftleft(2*b, tmp);
    gtog(r, tmp2); 
    mulg(d, tmp2);
    subg(tmp2, tmp);
    itog(1, tmp2);
    while (gsign(tmp) < 0) 
    {
        subg(tmp2, r);
        addg(d, tmp);
    }
    RETURN_GIANT(tmp);
    RETURN_GIANT(tmp2);
}
 
 
void
modg(
    giant   d,
    giant   n
)
/* n becomes n%d. n is arbitrary, but the denominator d must be positive! */
{
    if (cur_recip == NULL) {
        cur_recip = newgiantbits(1);
        cur_den = newgiantbits(1);
        gtog(d, cur_den);
        make_recip(d, cur_recip);
    } else if (gcompg(d, cur_den)) {
        gtog(d, cur_den);
        make_recip(d, cur_recip);
    }
    modg_via_recip(d, cur_recip, n);
}
 
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);
    giant   tmp, tmp2;
 
    if (isZero(d) || (gsign(d) < 0))
    {
        FATAL_ERROR(SIGN);
    }
    
    tmp = BORROW_GIANT();
    tmp2 = BORROW_GIANT();
    
    if (isone(d)) {
        // dividing by one!
        goto done;
    }
    
    absg(n);
    itog(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);
    setsign(n, gsign(n) * sign);
done:
    RETURN_GIANT(tmp);
    RETURN_GIANT(tmp2);
}
 
 
 
void
divg(
    giant   d,
    giant   n
)
/* n becomes n/d. n is arbitrary, but the denominator d must be positive! */
{
    if (cur_recip == NULL) {
        cur_recip = newgiantbits(1);
        cur_den = newgiantbits(1);
        gtog(d, cur_den);
        make_recip(d, cur_recip);
    } else if (gcompg(d, cur_den)) {
        gtog(d, cur_den);
        make_recip(d, cur_recip);
    }
    divg_via_recip(d, cur_recip, n);
}
 
static long
radixdiv(
    long                base,
    long                divisor,
    giant           thegiant
)
/* Divides giant of arbitrary base by divisor.
 * Returns remainder. Used by idivg and gread. */
{
    long                j = (thegiant->giantDigits*SHORTS_PER_VECTOR)-1;
    unsigned int        num,rem=0;
    unsigned short      *digitpointer;
 
    if (divisor == 0)
    {
        FATAL_ERROR(DIVIDEBYZERO);
    }
 
    while (j>=0)
    {
        digitpointer = GetNthShortPtr(thegiant, j);
        num=rem*base + *digitpointer;       
        *digitpointer = (unsigned short)(num/divisor);
 
        rem = num % divisor;
        --j;
    }
 
    if (divisor < 0) negg(thegiant);
 
    justg(thegiant);
 
    return(rem);
}
 
#if 1
long
idivg(
    long    divisor,
    giant   theg
)
{
    /* theg becomes theg/divisor. Returns remainder. */
    int     n;
    int     base = 1<<(8*sizeof(short));
 
    n = radixdiv(base,divisor,theg);
    return(n);
}
#else
long
idivg(
    long        divisor,
    giant       theg
)
{
    giant   originalg = BORROW_GIANT();
    giant   divisorg = BORROW_GIANT();
    long    remainder;
        
    gtog(theg, originalg);
    
    itog(divisor, divisorg);
    
    divg(divisorg, theg);
 
    mulg(theg, divisorg);
    
    subg(divisorg, originalg);
        
    remainder = gtoi(originalg);
    
    RETURN_GIANT(divisorg);
    RETURN_GIANT(originalg);
    
    return remainder;
}
#endif
 
 
 
 
 
static void
gswap(
    giant   *p,
    giant   *q
)
{
    giant   t;
 
    t = *p;
    *p = *q;
    *q = t;
}
 
 
static void
fix(
    giant   *p,
    giant   *q
)
/* Insure that x > y >= 0. */
{
    if( gsign(*p) < 0 )
        negg(*p);
    if( gsign(*q) < 0 )
        negg(*q);
    if( gcompg(*p,*q) < 0 )
        gswap(p,q);
}
 
 
unsigned long   numtrailzeros(giant g)
{   
    unsigned long trailingbits = 0;
 
    ASSERT_GVALID(g);
    
    if (g->giantDigits != 0) {
        vector unsigned long    veczero = (vector unsigned long)(0);
        long                    vecindex;
        long                    byteindex;
        long                    bitindex;
        unsigned char           lastbyte;
        
        for (vecindex = 0; vecindex < g->giantDigits; vecindex++) {
            if (!vec_all_eq(veczero, g->vectors[vecindex])) {
                break;
            }   
            
            trailingbits += BITS_PER_VECTOR;
        }
    
        // now go by bytes;
        for (byteindex = 0; byteindex < BYTES_PER_VECTOR; byteindex++) {
            lastbyte = GetNthByte(g, byteindex + (vecindex * BYTES_PER_VECTOR));
            if (lastbyte != 0) {
                break;
            } 
            
            trailingbits += 8;
        }
        
        
        for (bitindex=0; bitindex < 8; bitindex++) {
            if (lastbyte & 1) {
                break;
            }
            
            lastbyte = lastbyte >> 1;
            trailingbits++;
        }   
    }
    
    return trailingbits;
}
 
static void
dotproduct(
    giant   a,
    giant   b,
    giant   c,
    giant   d
)
/* Replace last argument with the dot product of two 2-vectors. */
{
    giant s4 = BORROW_GIANT();
 
    gtog(c,s4);
    mulg(a, s4);
    mulg(b,d);
    addg(s4,d);
    
    RETURN_GIANT(s4);
}
 
 
static   void
mulvM(
    gmatrix     A,
    giant       x,
    giant       y
)
/* Multiply vector by Matrix; changes x,y. */
{
    giant s0 = BORROW_GIANT();
    giant s1 = BORROW_GIANT();
    
    gtog(A->ur,s0);
    gtog( A->lr,s1);
    dotproduct(x,y,A->ul,s0);
    dotproduct(x,y,A->ll,s1);
    gtog(s0,x);
    gtog(s1,y);
    
    RETURN_GIANT(s0);
    RETURN_GIANT(s1);
}
 
 
static void
bgcdg(
    giant   a,
    giant   b
)
/* Binary form of the gcd. b becomes the gcd of a,b. */
{
    int     k = isZero(b), m = isZero(a);
    giant   u, v, t;
 
    if (k || m)
    {
        if (m)
        {
            if (k)
                itog(1,b);
            return;
        }
        if (k)
        {
            if (m)
                itog(1,b);
            else
                gtog(a,b);
            return;
        }
    }
 
    u = BORROW_GIANT();
    v = BORROW_GIANT();
    t = BORROW_GIANT();
 
    /* Now neither a nor b is zero. */
    gtog(a, u);
    absg(u);
    gtog(b, v);
    absg(v);
    k = numtrailzeros(u);
    m = numtrailzeros(v);
    if (k>m)
        k = m;
    gshiftright(k,u);
    gshiftright(k,v);
    if (isOdd(u)) {
        gtog(v, t);
        negg(t);
    }
    else
    {
        gtog(u,t);
    }
 
    while (!isZero(t))
    {       
        m = numtrailzeros(t);
        gshiftright(m, t);
        if (gsign(t) > 0)
        {
            gtog(t, u);
            subg(v,t);
        }
        else
        {
            gtog(t, v);
            negg(v);
            addg(u,t);
        }
    }
    gtog(u,b);
    gshiftleft(k, b);
    RETURN_GIANT(v);
    RETURN_GIANT(u);
    RETURN_GIANT(t);
}
 
static void
shgcd(
    register int    x,
    register int    y,
    gmatrix         A
)
/*
 * Do a half gcd on the integers a and b, putting the result in A
 * It is fairly easy to use the 2 by 2 matrix description of the
 * extended Euclidean algorithm to prove that the quantity q*t
 * never overflows.
 */
{
    register int    q, t, start = x;
    int             Aul = 1, Aur = 0, All = 0, Alr = 1;
 
    while   (y != 0 && y > start/y)
    {
        q = x/y;
        t = y;
        y = x%y;
        x = t;
        t = All;
        All = Aul-q*t;
        Aul = t;
        t = Alr;
        Alr = Aur-q*t;
        Aur = t;
    }
    itog(Aul,A->ul);
    itog(Aur,A->ur);
    itog(All,A->ll);
    itog(Alr,A->lr);
}
 
static void
punch(
    giant       q,
    gmatrix     A
)
/* Multiply the matrix A on the left by [0,1,1,-q]. */
{
    giant s0 = BORROW_GIANT();
    
    gtog(A->ll,s0);
    mulg(q,A->ll);
    gswap(&A->ul,&A->ll);
    subg(A->ul,A->ll);
    gtog(s0,A->ul);
    gtog(A->lr,s0);
    mulg(q,A->lr);
    gswap(&A->ur,&A->lr);
    subg(A->ur,A->lr);
    gtog(s0,A->ur);
 
    RETURN_GIANT(s0);
}
 
 
static void
onestep(
    giant       x,
    giant       y,
    gmatrix     A
)
/* Do one step of the euclidean algorithm and modify
 * the matrix A accordingly. */
{
    giant s4 = BORROW_GIANT();
 
    gtog(x,s4);
    gtog(y,x);
    gtog(s4,y);
    divg(x,s4);
    punch(s4,A);
    mulg(x,s4);
    subg(s4,y);
    
    RETURN_GIANT(s4);
}
 
static void
mulmM(
    gmatrix     A,
    gmatrix     B
)
/* Multiply matrix by Matrix; changes second matrix. */
{
    giant s0 = BORROW_GIANT();
    giant s1 = BORROW_GIANT();
    giant s2 = BORROW_GIANT();
    giant s3 = BORROW_GIANT();
 
    gtog(B->ul,s0);
    gtog(B->ur,s1);
    gtog(B->ll,s2);
    gtog(B->lr,s3);
 
    dotproduct(A->ur,A->ul,B->ll,s0);
    dotproduct(A->ur,A->ul,B->lr,s1);
    dotproduct(A->ll,A->lr,B->ul,s2);
    dotproduct(A->ll,A->lr,B->ur,s3);
 
    gtog(s0,B->ul);
    gtog(s1,B->ur);
    gtog(s2,B->ll);
    gtog(s3,B->lr);
 
    RETURN_GIANT(s0);
    RETURN_GIANT(s1);
    RETURN_GIANT(s2);
    RETURN_GIANT(s3);
}
 
 
 
static void
hgcd(
    int         n,
    giant       xx,
    giant       yy,
    gmatrix     A
)
/* hgcd(n,x,y,A) chops n bits off x and y and computes th
 * 2 by 2 matrix A such that A[x y] is the pair of terms
 * in the remainder sequence starting with x,y that is
 * half the size of x. Note that the argument A is modified
 * but that the arguments xx and yy are left unchanged.
 */
{
    giant       x, y;
 
    if (isZero(yy))
        return;
    
    x = BORROW_GIANT();
    y = BORROW_GIANT();
    gtog(xx,x);
    gtog(yy,y);
    gshiftright(n,x);
    gshiftright(n,y);
    if (bitlen(x) <= INTLIMIT )
    {
        shgcd(gtoi(x),gtoi(y),A);
    }
    else
    {
        GMatrixStruct   B;
        int             m = bitlen(x)/2;
 
        hgcd(m,x,y,A);
        mulvM(A,x,y);
        if (gsign(x) < 0)
        {
            negg(x); negg(A->ul); negg(A->ur);
        }
        if (gsign(y) < 0)
        {
            negg(y); negg(A->ll); negg(A->lr);
        }
        if (gcompg(x,y) < 0)
        {
            gswap(&x,&y);
            gswap(&A->ul,&A->ll);
            gswap(&A->ur,&A->lr);
        }
        if (!isZero(y))
        {
            onestep(x,y,A);
            m /= 2;
            B.ul = BORROW_GIANT();
            B.ur = BORROW_GIANT();
            B.ll = BORROW_GIANT();
            B.lr = BORROW_GIANT();
            itog(1,B.ul);
            itog(0,B.ur);
            itog(0,B.ll);
            itog(1,B.lr);
            hgcd(m,x,y,&B);
            mulmM(&B,A);
            
            RETURN_GIANT(B.ul);
            RETURN_GIANT(B.ur);
            RETURN_GIANT(B.ll);
            RETURN_GIANT(B.lr);
        }
    }
    RETURN_GIANT(x);
    RETURN_GIANT(y);
}
 
 
 
 
static void
ggcd(
    giant   xx,
    giant   yy
)
/* A giant gcd.  Modifies its arguments. */
{
    giant           x = BORROW_GIANT();
    giant           y = BORROW_GIANT();
    GMatrixStruct   A;
 
    gtog(xx,x); gtog(yy,y);
    for(;;)
    {
        fix(&x,&y);
        if (bitlen(y) <= GCDLIMIT )
            break;
        A.ul = BORROW_GIANT();
        A.ur = BORROW_GIANT();
        A.ll = BORROW_GIANT();
        A.lr = BORROW_GIANT();
        itog(1,A.ul);
        itog(0,A.ur);
        itog(0,A.ll);
        itog(1,A.lr);
        hgcd(0,x,y,&A);
        mulvM(&A,x,y);
        RETURN_GIANT(A.ul);
        RETURN_GIANT(A.ur);
        RETURN_GIANT(A.ll);
        RETURN_GIANT(A.lr);
        fix(&x,&y);
        if (bitlen(y) <= GCDLIMIT )
            break;
        modg(y,x);
        gswap(&x,&y);
    }
    bgcdg(x,y);
    gtog(y,yy);
    RETURN_GIANT(x);
    RETURN_GIANT(y);
}
 
void
gcdg(
    giant       x,
    giant       y
)
{
    if (bitlen(y)<= GCDLIMIT)
        bgcdg(x,y);
    else
        ggcd(x,y);
}
 
 
void
cgcdg(
    giant   a,
    giant   v
)
/* Classical Euclid GCD. v becomes gcd(a, v). */
{
    giant   u, r;
 
    absg(v);
    if (isZero(a))
        return;
    
    u = BORROW_GIANT();
    r = BORROW_GIANT();
    
    gtog(a, u);
    absg(u);
    while (!isZero(v))
    {
        gtog(u, r);
        modg(v, r);
        gtog(v, u);
        gtog(r, v);
    }
    gtog(u,v);
    RETURN_GIANT(u);
    RETURN_GIANT(r);
}
 
 
 
 
static void
columnwrite(
    FILE    *filepointer,
    short   *column,
    char    *format,
    short   arg,
    int     newlines
)
/* Used by gwriteln. */
{
    char    outstring[10];
    short   i;
 
    sprintf(outstring,format,arg);
    for (i=0; outstring[i]!=0; ++i)
    {
        if (newlines)
        {
            if (*column >= COLUMNWIDTH)
            {
                fputs("\\\n",filepointer);
                *column = 0;
            }
        }
        fputc(outstring[i],filepointer);
        ++*column;
    }
}
 
 
void
gwrite(
    giant           thegiant,
    FILE            *filepointer,
    int             newlines
)
/* Outputs thegiant to filepointer. Output is terminated by a newline. */
{
    short           column;
    static giant    garbagegiant = NULL;
    unsigned long   *pLongStorage;
    unsigned long   numlongs = (thegiant->giantDigits * 40);    
    unsigned long   *numpointer;
    Boolean         firstdigits=true;
    
    pLongStorage = (unsigned long *)NewPtr(numlongs*sizeof(long));
    
    if (garbagegiant == NULL)
        garbagegiant = newgiantbits(1);
 
    if (isZero(thegiant))
    {
        fputs("0",filepointer);
    }
    else
    {
        numpointer = pLongStorage;
        gtog(thegiant,garbagegiant);
        absg(garbagegiant);
 
        do
        {
            *numpointer = (unsigned short)idivg(10000,garbagegiant);
            ++numpointer;
        }   while (!isZero(garbagegiant));
 
        column = 0;
        numpointer--;
 
        if (gsign(thegiant)<0) {
            columnwrite(filepointer,&column,"-",0, newlines);
        }
        
        do {
            if (firstdigits) {
                columnwrite(filepointer,&column,"%d",*numpointer--,newlines);       
                firstdigits = false;
            } else {
                columnwrite(filepointer,&column,"%04d",*numpointer--,newlines);                 
            }
        } while (numpointer >= pLongStorage);
        
    }
    
    DisposePtr((Ptr)pLongStorage);
}
 
void
gwriteln(
    giant       theg,
    FILE        *filepointer
)
{
    gwrite(theg, filepointer, 1);
    fputc('\n',filepointer);
}
 
 
void
gread(
    giant           theg,
    FILE            *filepointer
)
{
    char            currentchar;
    int             isneg,size,backslash=false,numdigits=0;
    static giant    basetenthousand = NULL;
    static char     *inbuf = NULL;
    static giant    currentmult;
    
    if (basetenthousand == NULL)
        basetenthousand = newgiantbits(1);
    if (currentmult == NULL)
        currentmult = newgiantbits(1);
    if (inbuf == NULL)
        inbuf = (char*)malloc(MAX_DIGITS);
 
    itog(1, currentmult);
    itog(0, theg);
    
    currentchar = (char)fgetc(filepointer);
    if (currentchar=='-')
    {
        isneg=true;
    }
    else
    {
        isneg=false;
        if (currentchar!='+')
            ungetc(currentchar,filepointer);
    }
 
    do
    {
        currentchar = (char)fgetc(filepointer);
        if ((currentchar>='0') && (currentchar<='9'))
        {
            inbuf[numdigits]=currentchar;
            if(++numdigits==MAX_DIGITS)
                break;
            backslash=false;
        }
        else
        {
            if (currentchar=='\\')
                backslash=true;
        }
    } while(((currentchar!=' ') && (currentchar!='\n') &&
                (currentchar!='\t')) || (backslash) );
    if (numdigits)
    {
        size = 0;
        do
        {
            inbuf[numdigits] = 0;
            numdigits-=4;
            if (numdigits<0)
                numdigits=0;
                
            itog((unsigned short)strtol(&inbuf[numdigits],NULL,10), basetenthousand);
            mulg(currentmult, basetenthousand);
            addg(basetenthousand, theg);
            itog(10000, basetenthousand);
            mulg(basetenthousand, currentmult);
 
        } while (numdigits>0);
 
 
        if (isneg) {
            setsign(theg, -1);
        } else {
            setsign(theg, 1);
        }
    }
}
 
void        gdumphex(giant g)
{
    long            i;
    unsigned char   b;
    
    if (g->giantDigits) {       
        if (g->giantSign == -1) {
            printf("-");
        }
 
        for (i=BYTES_PER_VECTOR*(g->giantDigits)-1; i>=0; i--) {
            b = GetNthByte(g, i);
            if (b == 0) {
                printf("00");
            } else if (b<0x10) {
                printf("0%x", b);
            } else {
                printf("%x", b);
            }
        }
    } else {
        printf("0");
    }
}
 
 
 
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)
    {
        itog(0,u);
        return;
    }
    scratch7 = BORROW_GIANT();
    gtog(v, scratch7);
    gshiftleft(diff,scratch7);
    if (gcompg(u,scratch7) < 0)
        diff--;
    if (diff<0)
    {
        itog(0,u);
        RETURN_GIANT(scratch7);
        return;
    }
    itog(1,u);
    gshiftleft(diff,u);
 
    RETURN_GIANT(scratch7);
}
 
 
 
static giant    u0=NULL, u1=NULL, v0=NULL, v1=NULL;
 
 
static 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;
 
    if (isone(x))
        return(1);
    if (!v0)
        v0 = newgiantbits(1);
    if (!v1)
        v1 = newgiantbits(1);
    if (!u1)
        u1 = newgiantbits(1);
    if (!u0)
        u0 = newgiantbits(1);
    itog(1, v0);
    gtog(x, v1);
    itog(0,x);
    gtog(p, u1);
 
    scratch7 = BORROW_GIANT();
    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);
    }
    
    RETURN_GIANT(scratch7);
 
    if (!isone(u1))
    {
        gtog(u1,x);
        if(gsign(x)<0) addg(p, x);
            return(0);
    }
    if(gsign(x)<0)
        addg(p, x);
    return(1);
}
 
int
binvg(
    giant   p,
    giant   x
)
{
    modg(p, x);
    return(binvaux(p,x));
}
 
static int
invaux(
    giant   p,
    giant   x
)
/* Returns zero if no inverse exists, in which case x becomes
 * GCD(x,p). */
{
 
    giant scratch7;
    
    if (isone(x))
        return(1);
    if (!v0)
        v0 = newgiantbits(1);
    if (!v1)
        v1 = newgiantbits(1);
    if (!u1)
        u1 = newgiantbits(1);
    if (!u0)
        u0 = newgiantbits(1);
    itog(1,u1);
    gtog(p, v0);
    gtog(x, v1);
    itog(0,x);
 
    scratch7 = BORROW_GIANT();
    while (!isZero(v1))
    {
        gtog(v0, u0);
        divg(v1, u0);
        gtog(u0, scratch7);
        mulg(v1, scratch7);
        subg(v0, scratch7);
        negg(scratch7);
        gtog(v1, v0);
        gtog(scratch7, v1);
        gtog(u1, scratch7);
        mulg(u0, scratch7);
        subg(x, scratch7);
        negg(scratch7);
        gtog(u1,x);
        gtog(scratch7, u1);
    }
 
    RETURN_GIANT(scratch7);
    
    if (!isone(v0))
    {
        gtog(v0,x);
        return(0);
    }
    if(gsign(x)<0)
        addg(p, x);
    return(1);
}
 
 
int
invg(
    giant   p,
    giant   x
)
{
    modg(p, x);
    return(invaux(p,x));
}
 
void
mersennemod (
    int n,
    giant g
)
/* g := g (mod 2^n - 1) */
{
    int the_sign, s;
    giant scratch3 = BORROW_GIANT();
    giant scratch4 = BORROW_GIANT();
    
    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)) {
        if ((s = gsign(g)) < 0) absg(g);
        itog(1,scratch3);
        gshiftleft(n,scratch3);
        itog(1,scratch4);
        subg(scratch4,scratch3);
        if(gcompg(g,scratch3) >= 0) subg(scratch3,g);
        if (s < 0) {
            negg(g);
            addg(scratch3,g);
        }
        if (the_sign < 0) {
            negg(g);
            addg(scratch3,g);
        }
    }
    
    RETURN_GIANT(scratch3);
    RETURN_GIANT(scratch4);
}
 
 
 
int
mersenneinvg(
    int     q,
    giant   x
)
{
    int     k;
 
    if (!u0)
        u0 = newgiantbits(1);
    if (!u1)
        u1 = newgiantbits(1);
    if (!v1)
        v1 = newgiantbits(1);
 
    itog(1, u0);
    itog(0, u1);
    itog(1, v1);
    gshiftleft(q, v1);
    subg(u0, v1);
    mersennemod(q, x);
    while (1)
    {
        k = -1;
        if (isZero(x))
        {
            gtog(v1, x);
            return(0);
        }
        while (bitval(x, ++k) == 0)
            ;
 
        gshiftright(k, x);
        if (k)
        {
            gshiftleft(q-k, u0);
            mersennemod(q, u0);
        }
        if (isone(x))
            break;
        addg(u1, u0);
        mersennemod(q, u0);
        negg(u1);
        addg(u0, u1);
        mersennemod(q, u1);
        if (!gcompg(v1,x))
            return(0);
        addg(v1, x);
        negg(v1);
        addg(x, v1);
        mersennemod(q, v1);
    }
    gtog(u0, x);
    mersennemod(q,x);
    return(1);
}
 
void
powermod(
    giant       x,
    int         n,
    giant       g
)
/* x becomes x^n (mod g). */
{
    giant scratch2 = BORROW_GIANT();
    gtog(x, scratch2);
    itog(1, x);
    while (n)
    {
        if (n & 1)
        {
            mulg(scratch2, x);
            modg(g, x);
        }
        n >>= 1;
        if (n)
        {
            squareg(scratch2);
            modg(g, scratch2);
        }
    }
    
    RETURN_GIANT(scratch2);
}
 
 
void
powermodg(
    giant       x,
    giant       n,
    giant       g
)
/* x becomes x^n (mod g). */
{
    int         len, pos;
    giant       scratch2 = BORROW_GIANT();
 
    gtog(x, scratch2);
    itog(1, x);
    len = bitlen(n);
    pos = 0;
    while (1)
    {
        if (bitval(n, pos++))
        {
            mulg(scratch2, x);
            modg(g, x);
        }
        if (pos>=len)
            break;
        squareg(scratch2);
        modg(g, scratch2);
    }
    RETURN_GIANT(scratch2);
}
 
 
void
fermatpowermod(
    giant   x,
    int     n,
    int     q
)
/* x becomes x^n (mod 2^q+1). */
{
    giant scratch2 = BORROW_GIANT();
    
    gtog(x, scratch2);
    itog(1, x);
    while (n)
    {
        if (n & 1)
        {
            mulg(scratch2, x);
            fermatmod(q, x);
        }
        n >>= 1;
        if (n)
        {
            squareg(scratch2);
            fermatmod(q, scratch2);
        }
    }
    RETURN_GIANT(scratch2);
}
 
 
void
fermatpowermodg(
    giant   x,
    giant   n,
    int     q
)
/* x becomes x^n (mod 2^q+1). */
{
    int     len, pos;
    giant   scratch2 = BORROW_GIANT();
 
    gtog(x, scratch2);
    itog(1, x);
    len = bitlen(n);
    pos = 0;
    while (1)
    {
        if (bitval(n, pos++))
        {
            mulg(scratch2, x);
            fermatmod(q, x);
        }
        if (pos>=len)
            break;
        squareg(scratch2);
        fermatmod(q, scratch2);
    }
    RETURN_GIANT(scratch2);
}
 
 
void
mersennepowermod(
    giant   x,
    int     n,
    int     q
)
/* x becomes x^n (mod 2^q-1). */
{
    giant scratch2 = BORROW_GIANT();
 
    gtog(x, scratch2);
    itog(1, x);
    while (n)
    {
        if (n & 1)
        {
            mulg(scratch2, x);
            mersennemod(q, x);
        }
        n >>= 1;
        if (n)
        {
            squareg(scratch2);
            mersennemod(q, scratch2);
        }
    }
    RETURN_GIANT(scratch2);
}
 
 
void
mersennepowermodg(
    giant   x,
    giant   n,
    int     q
)
/* x becomes x^n (mod 2^q-1). */
{
    int     len, pos;
    giant   scratch2 = BORROW_GIANT();
 
    gtog(x, scratch2);
    itog(1, x);
    len = bitlen(n);
    pos = 0;
    while (1)
    {
        if (bitval(n, pos++))
        {
            mulg(scratch2, x);
            mersennemod(q, x);
        }
        if (pos>=len)
            break;
        squareg(scratch2);
        mersennemod(q, scratch2);
    }
    RETURN_GIANT(scratch2);
}
 
 
 
void
fermatmod (
    int             n,
    giant           g
)
/* g := g (mod 2^n + 1). */
{
    int the_sign, s;
    giant scratch3 = BORROW_GIANT();
    
    if ((the_sign = gsign(g)) < 0) absg(g);
    while (bitlen(g) > n) {
        gtog(g,scratch3);
        gshiftright(n,scratch3);
        subg(scratch3,g);
        gshiftleft(n,scratch3);
        subg(scratch3,g);
    }
    if(!isZero(g)) {
        if ((s = gsign(g)) < 0) absg(g);
        itog(1,scratch3);
        gshiftleft(n,scratch3);
        iaddg(1,scratch3);
        if(gcompg(g,scratch3) >= 0) subg(scratch3,g);
        if (s < 0) {
            negg(g);
            addg(scratch3,g);
        }
        if (the_sign < 0) {
            negg(g);
        }
    }
    
    RETURN_GIANT(scratch3);
}
 
 
#if 1
void
smulg(
    unsigned long   i,
    giant           g
)
{
    if (!i) {
        g->giantDigits = 0;
        return;
    }
    
    if (!isZero(g)) {   
        unsigned long topresult = MultVecsByULong(  g->vectors,
                                                    g->giantDigits,
                                                    i);
 
        if (topresult) {
            AddMSDLong(g, topresult);
        }
    }
    
    ASSERT_GVALID(g);
}
#else
void
smulg(
    unsigned short  i,
    giant           g
)
/* g becomes g * i. */
{
    giant   temp = BORROW_GIANT();
    
    itog(i, temp);
    
    mulg(temp, g);
    
    RETURN_GIANT(temp);
 
}
#endif
 
long            icompg(unsigned long i, giant g)
{
    long    result;
 
    giant   temp = BORROW_GIANT();
    
    itog(i, temp);
    result = gcompg(temp, g);   
    
    RETURN_GIANT(temp);
    
    return result;
}
 
 
void        trunctobitlen(long len, giant g)
{
    unsigned long   neededDigits;
    unsigned long   trimbits;
    unsigned char   *pDigit;
    unsigned char   maskByte;
    long            i;
    
    if (bitlen(g) < len) return;
 
    neededDigits = (len + BITS_PER_VECTOR - 1)/BITS_PER_VECTOR;
    
    g->giantDigits = neededDigits;
    
    trimbits = len & 0x7f;
    
    trimbits = BITS_PER_VECTOR - trimbits;
    
    if (trimbits) {
        pDigit = (unsigned char*)&g->vectors[neededDigits-1];
 
        for (i=0; i<trimbits >> 3; i++) {
            pDigit[i] = 0;
        }   
        
        trimbits &= 0x07;
        maskByte = 0xff;
        
        maskByte = maskByte >> trimbits;
        
        pDigit[i] &= maskByte;
        
    }
}
 
 
#if !INLINE_GIANTS
 
/* Returns 1, 0, -1 as a>b, a=b, a<b. */
long gcompg(giant a, giant b)
{
    ASSERT_GVALID(a);
    ASSERT_GVALID(b);
 
    if (a->giantSign > b->giantSign) return 1;
    if (a->giantSign < b->giantSign) return -1;
    
    return a->giantSign * CompareMagnitude(a, b);
}
 
signed long     gtoi (giant g)
{
    GAssert(bitlen(g) <= 32);
    
    if (!g->giantSign) return 0;
    
    return g->giantSign * ((unsigned long *)g->vectors)[LONGS_PER_VECTOR-1];
}
 
/* Returns the giantSign of g: -1, 0, 1. */
long            gsign(giant g)
{
    ASSERT_GVALID(g);
    
    if (!g->giantDigits) return 0;
    
    return g->giantSign;
}
 
void        negg(giant g)
{
    g->giantSign = -g->giantSign;
}
 
void        setsign(giant g, long sign)
{
    if (sign > 0) {
        g->giantSign = 1;
    } else {
        g->giantSign = -1;
    }
}
 
void        absg(giant g)
{
    if (g->giantSign < 0) g->giantSign = 1;
}
 
long isZero(giant g)
{
    ASSERT_GVALID(g);
 
    return (g->giantDigits == 0);
}
 
Boolean isone(giant g)
{
    Boolean                 result = false;
    unsigned long           *pLongPtr;
    unsigned long           accumulator;
    long                    counter;
    
    if (g->giantSign == 1 && g->giantDigits == 1) {
        pLongPtr = (unsigned long*)g->vectors;
        accumulator = 0;
        for (counter=0; counter < LONGS_PER_VECTOR -1 ; counter++) {
            accumulator |= pLongPtr[counter];
        }
        
        // if first three longs are zero
        if (!accumulator) {
            // and last one is one, then it's one
            result = (pLongPtr[counter] == 1);
        }
    }
    
    return result;
 
}
 
Boolean         isOdd(giant g)
{
    return (GetNthByte(g, 0) & 0x01);
}
 
 
#endif //#if INLINE_GIANTS


Copyright © 2003 Apple Computer, Inc. All Rights Reserved. Terms of Use | Privacy Policy | Updated: 2003-01-14