/*
 * giants.h 
 *
 * Header file for large-integer arithmetic library
 * giants.c. 
 *
 *
 */

/* 2^(16*MAX_SHORTS)-1 will fit into a giant. */
#define MAX_SHORTS ((1<<20)+(1<<8))
#define INFINITY (-1)
#define FALSE 0
#define TRUE 1
#define COLUMNWIDTH 64
#define DIVIDEBYZERO  1		/* error codes */
#define OVERFLOW      2
#define SIGN	      3
#define OVERRANGE     4
#define AUTO_MUL 0
#define GRAMMAR_MUL 1
#define FFT_MUL 2
#define USE_ASSEMBLER_MUL 1
#define TWOPI (2 * 3.14159265358979323846)
#define SQRT2 1.41421356237309504880
#define SQRTHALF 0.70710678118654752440
#define TWO16 65536.0
#define TWOM16 0.0000152587890625
#define MAX_DIGITS 10000	/* Decimal digit ceiling in I/O
				 * routines. */
#define BREAK_SHORTS 400	/* Number of shorts at which FFT
				 * breaks over. */
#define min(a,b) ((a)<(b)? (a) : (b))
#define max(a,b) ((a)>(b)? (a) : (b))
#define STEPS 32		/* Maximum number of recursive steps
				 * needed to calculate gcds of
				 * integers */
#define GCDLIMIT 5000		/* The limit below which hgcd is too
				 * ponderous */
#define INTLIMIT  31		/* The limit below which ordinary
				 * ints will be used */
#define gin(x)   gread(x,stdin)
#define gout(x)  gwriteln(x,stdout)

typedef struct _giant {
    int sign;			/* number of shorts = abs(sign) */
    unsigned short n[1];
} *giant;

typedef struct _gmatrix {
    giant ul, ur, ll, lr;	/* upper left, upper right, etc. */
} *gmatrix;

typedef struct {
    double re, im;
} complex;

/*
 * Creates a new giant, numshorts = INFINITY invokes the maximum
 * MAX_SHORTS. 
 */
giant newgiant(int numshorts);

/* Returns the bit-length n; e.g. n=7 returns 3. */
int bitlen(giant n);

/* Returns the value of the pos bit of n. */
int bitval(giant n, int pos);

int iszero(giant g);		/* Returns whether g is zero. */
int isone(giant g);		/* Returns whether g is one. */

/* Copies one giant to another. */
void gtog(giant src, giant dest);

/* Gives a giant an int value. */
void itog(int n, giant g);

/* Returns the sign of g: -1, 0, 1. */
int gsign(giant g);

/* Returns 1, 0, -1 as a>b, a=b, a<b. */
int gcompg(giant a, giant b);	

/* g := g/n, and (g mod n) is returned. */
int idivg(int n, giant g);

/*
 * Output the giant in decimal, with neither '\'-newline notation at
 * right terminal margin, nor a final newline.  
 */
void gwrite(giant g, FILE * fp, int newlines);	

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

void negg(giant g);		/* g := -g. */
void addg(giant a, giant b);	/* b += a. */
void subg(giant a, giant b);	/* b -= a. */

/* g += i, with i non-negative and < 2^16. */
void iaddg(int i, giant g);

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

/*
 * Stable-divide version of invg, init_divide() should be called
 * first if the denominator is new. 
 */
int stableinvg(giant n, giant x);
/* Same as stableinvg, except uses binary method. */
int stablebinvg(giant p, giant x);

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

void cgcdg(giant n, giant x);	/* Classical GCD, x:= GCD(n, x). */
void gcdg(giant n, giant x);	/* General GCD, x:= GCD(n, x). */
void bgcdg(giant n, giant x);	/* Binary GCD, x:= GCD(n, x). */

/* g := m^n, no mod is performed. */
void powerg(int a, int b, giant g);

/* num := num mod den, any positive den. */
void modg(giant den, giant num);

/* num := |num|/den, any positive den. */
void divg(giant den, giant num);

/* same as modg but  reciprocal storage is invoked for den. */
void stablemodg(giant den, giant num);

/* same as divg but reciprocal storage is  invoked for den. */
void stabledivg(giant den, giant num);

/* x := x^n (mod z). */
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 := g (mod 2^n-1). */
void mersennemod(int n, giant g);

void smulg(unsigned short s, giant g);	/* g *= s. */
void absg(giant g);			/* g := |g|. */
void squareg(giant g);			/* g *= g. */
void mulg(giant a, giant b);		/* b *= a. */

/* Shift g left by bits, introducing zeros on the right. */
void gshiftleft(int bits, giant g);
/* Shift g right by bits, losing bits on the right. */
void gshiftright(int bits, giant g);

/*
 * Set AUTO_MUL for automatic FFT crossover (this is the default),
 * set FFT_MUL for forced FFT multiply, set GRAMMAR_MUL for forced
 * grammar school multiply. 
 */
void setmulmode(int mode);