/**************************************************************

 *

 *  vgiants.h

 *

 *  AltiVec-based giant integer package

 *  Derived from original giants.c project at 

 *  Next Software, Inc.

 *

 *  Updates:

 *

 *  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.

*/

#ifndef __VECTORGIANTS_H__

#define __VECTORGIANTS_H__

#if __MWERKS__

#include <altivec.h>

#endif

#define vector __vector

#include <Memory.h>

#include <stdio.h>

#include "giantsdebug.h"

#ifdef __cplusplus

extern "C" {

#endif

#ifdef __cplusplus

#define INLINE_GIANTS 1

#else

#define INLINE_GIANTS 0     // if not c++, don't inline

#endif

/**************************************************************

 *

 * Error Codes

 *

 **************************************************************/

#define DIVIDEBYZERO  1

#define OVFLOW      2

#define SIGN          3

/**************************************************************

**************************************************************/

    /* vectors beyond which we resort to Karatsuba mul. */

    #define KARAT_BREAK 24

    /* The limit below which hgcd is too ponderous */

    #define GCDLIMIT 5000

    /* The limit below which ordinary ints will be used */

    #define INTLIMIT  31

    /* columns for gwrite output */

    #define COLUMNWIDTH 64

    /* Decimal digit ceiling in I/O routines. */

    #define MAX_DIGITS 10000

#define LONGS_PER_VECTOR    4

#define SHORTS_PER_VECTOR   8

#define BYTES_PER_VECTOR    16

#define BITS_PER_VECTOR     128

#define SIGN_NEGATIVE       (-1)

#define SIGN_POSITIVE       (1)

#define DEFAULT_GIANT_BITS  128

#define TWOPI (double)(2*3.1415926535897932384626433)

#define gin(x)   gread(x,stdin)

#define gout(x)  gwriteln(x,stdout)

typedef struct giantStruct {

    signed long             giantSign;      // +/- 1

    unsigned long           giantDigits;

    unsigned long           giantCapacity;

    vector unsigned long    *vectors;

} GiantStruct, *GiantStructPtr;

// for compatability with old code:

typedef GiantStructPtr  giant;

typedef GiantStruct     giantStruct;

typedef struct _matrix

{

     giant              ul;         /* upper left */

     giant              ur;         /* upper right */

     giant              ll;         /* lower left */

     giant              lr;         /* lower right */

} GMatrixStruct, *GMatrixPtr;

typedef GMatrixPtr gmatrix; // for compatability with old code

//-------------------------------------------------

#if GDEBUG 

void            ASSERT_GVALID(giant g);

#else

#define         ASSERT_GVALID(g)

#endif

//-------------------------------------------------

/* trig lookups. */

void            init_sinCos(int);

double          s_sin(int);

double          s_cos(int);

void mulg(giant a, giant b);

void karatmulg(giant x, giant y);

void grammarmulg(giant a, giant b);

unsigned char GetNthByte(giant a, unsigned long byte);

GiantStructPtr  newgiantbits(unsigned long capacityBits);

GiantStructPtr  newgiant(unsigned long capacityShorts);

void            disposegiant(GiantStructPtr g);         

/* Creates a new giant matrix, but does not malloc

 * the component giants. */

gmatrix         newgmatrix(void);

/* Returns 1, 0, -1 bitlen(a) > bitlen(b), or =, or < */

long CompareMagnitude(giant a, giant b);

/* g += i, with i non-negative. */

void            iaddg(unsigned long i,giant g);

/* Shift g left by bits, introducing zeros on the right. */

void            gshiftleft(unsigned long bits, giant g);

/* Shift g right by bits, losing bits on the right. */

void            gshiftright(unsigned long bits, giant g);

/* Integer <-> giant. */

void            itog(long n, giant g);

long            isZero(giant g);

/* n := n % d, d positive, using stored reciprocal directly. */

void            modg_via_recip(giant d, giant r, giant n);

/* num := num % den, any positive den. */

void            modg(giant den, giant num);

/* r becomes the steady-state reciprocal 2^(2b)/d, where

 * b = bit-length of d-1. */

void            make_recip(giant d, giant r);

/* n := [n/d], d positive, using stored reciprocal directly. */

void        divg_via_recip(giant d, giant r, giant n);

/* num := [num/den], any positive den. */

void        divg(giant den, giant num);

/* g := g/n, and (g mod n) is returned. */

long        idivg(long n, giant g);

void        setsign(giant g, long sign);

/* Returns the number of trailing zero bits in g. */

unsigned long   numtrailzeros(giant g);

/* General GCD, x:= GCD(n, x). */

void        gcdg(giant n, giant x);

/* Classical GCD, x:= GCD(n, x). */

void        cgcdg(giant n, giant x);

/* g *= g. */

void        squareg(giant g);

void        randomg(giant g, unsigned long numvecs);

void            addg(giant a, giant b);

void            addglongs(giant a, giant b);

void            subg(giant a, giant b);

void            negg(giant g);

void            absg(giant g);

void            gtog(giant a, giant b);

long            gsign(giant g);

long            isZero(giant g);

unsigned long   bitlen(giant g);

long            bitval(giant g, long pos);

/* Output the giant in decimal, with optional newlines.  */

void        gwrite(giant g, FILE *fp, int newlines);

void        gdumphex(giant g);

/* Output the giant in decimal, with both '\'-newline

 * notation and a final newline. */

void        gwriteln(giant g, FILE *fp);

/* Input the giant in decimal, assuming the formatting of

 * 'gwriteln'. */

void        gread(giant g, FILE *fp);

/* u becomes greatest power of two not exceeding u/v. */

void        bdivg(giant v, giant u);

/* Same as invg, but uses bdivg. */

int         binvg(giant n, giant x);

/* If 1/x exists (mod n), 1 is returned and x := 1/x.  If

 * inverse does not exist, 0 is returned and x := GCD(n, x). */

int         invg(giant n, giant x);

int         mersenneinvg(int q, giant x);

/* g := g (mod 2^n-1). */

void        mersennemod(int n, giant g);

/* num := [num/den], any positive den. */

void        powermod(giant x, int n, giant z);

/* x := x^n (mod z). */

void        powermodg(giant x, giant n, giant z);

/* x := x^n (mod 2^q+1). */

void        fermatpowermod(giant x, int n, int q);

/* x := x^n (mod 2^q+1). */

void        fermatpowermodg(giant x, giant n, int q);

/* x := x^n (mod 2^q-1). */

void        mersennepowermod(giant x, int n, int q);

/* x := x^n (mod 2^q-1). */

void        mersennepowermodg(giant x, giant n, int q);

/* g := g (mod 2^n+1). */

void        fermatmod(int n, giant g);

/* g *= s. */

void        smulg(unsigned long s, giant g);

void        trunctobitlen(long len, giant g);

long            icompg(unsigned long i, giant g);

#if INLINE_GIANTS

/* Returns 1, 0, -1 as a>b, a=b, a<b. */

inline 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);

}

inline 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. */

inline long             gsign(giant g)

{

    ASSERT_GVALID(g);

    if (!g->giantDigits) return 0;

    return g->giantSign;

}

inline void         negg(giant g)

{

    g->giantSign = -g->giantSign;

}

inline void     setsign(giant g, long sign)

{

    if (sign > 0) {

        g->giantSign = 1;

    } else {

        g->giantSign = -1;

    }

}

inline void     absg(giant g)

{

    if (g->giantSign < 0) g->giantSign = 1;

}

inline Boolean          isOdd(giant g)

{

    return (GetNthByte(g, 0) & 0x01);

}

inline 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;

}

inline long isZero(giant g)

{

    ASSERT_GVALID(g);

    return (g->giantDigits == 0);

}

#else

/* Returns 1, 0, -1 as a>b, a=b, a<b. */

long            gcompg(giant a, giant b);

signed long     gtoi (giant g);

long            gsign(giant g);

void            negg(giant g);

void            setsign(giant g, long sign);

void            absg(giant g);

long            isZero(giant g);

Boolean         isone(giant g);

Boolean         isOdd(giant g);

#endif

#ifdef __cplusplus

}

#endif

#endif //#ifndef __VECTORGIANTS_H__