giants.h

/**************************************************************
 *
 *	giants.h
 *
 *	Header file for large-integer arithmetic library giants.c.
 *
 *	Updates:
 *          18 Jul 99  REC  Added fer_mod().
 *          30 Apr 98  JF   USE_ASSEMBLER_MUL removed
 *          29 Apr 98  JF   Function prototypes cleaned up
 *			20 Apr 97  RDW
 *
 *	c. 1997 Perfectly Scientific, Inc.
 *	All Rights Reserved.
 *
 **************************************************************/


/**************************************************************
 *
 * Error Codes
 *
 **************************************************************/

#define DIVIDEBYZERO  1
#define OVFLOW      2
#define SIGN	      3
#define OVERRANGE     4
#define AUTO_MUL 0
#define GRAMMAR_MUL 1
#define FFT_MUL 2
#define KARAT_MUL 3

/**************************************************************
 *
 * Preprocessor definitions
 *
 **************************************************************/

/* 2^(16*MAX_SHORTS)-1 will fit into a giant, but take care:
 * one usually has squares, etc. of giants involved, and
 * every intermediate giant in a calculation must fit into
 * this many shorts. Thus, if you want systematically to effect
 * arithmetic on B-bit operands, you need MAX_SHORTS > B/8,
 * perferably a tad larger than this; e.g. MAX_SHORTS > B/7.
 */
#define MAX_SHORTS (1<<19)	

#define GIANT_INFINITY (-1)
#define FA 0
#define TR 1
#define COLUMNWIDTH 64

#define TWOPI (double)(2*3.1415926535897932384626433)
#define SQRT2 (double)(1.414213562373095048801688724209)
#define SQRTHALF (double)(0.707106781186547524400844362104)
#define TWO16 (double)(65536.0)
#define TWOM16 (double)(0.0000152587890625)

/* Decimal digit ceiling in digit-input routines. */
#define MAX_DIGITS 10000

/* Next, mumber of shorts per operand 
   at which Karatsuba breaks over. */
#define KARAT_BREAK 40

/* Next, mumber of shorts per operand 
   at which FFT breaks over. */
#define FFT_BREAK 200

#define newmin(a,b) ((a)<(b)? (a) : (b))
#define newmax(a,b) ((a)>(b)? (a) : (b))

/* Maximum number of recursive steps needed to calculate
 * gcds of integers. */
#define STEPS 32

/* The limit below which hgcd is too ponderous */
#define GCDLIMIT 5000

/* The limit below which ordinary ints will be used */
#define INTLIMIT  31

/* Size by which to increment the stack used in pushg() and popg(). */
#define	STACK_GROW	16

#define gin(x)   gread(x,stdin)
#define gout(x)  gwriteln(x,stdout)


/**************************************************************
 *
 * Structure definitions
 *
 **************************************************************/

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

typedef giantstruct *giant;

typedef struct _matrix
{
	 giant 				ul;			/* upper left */
	 giant				ur;         /* upper right */
	 giant				ll;         /* lower left */
	 giant				lr;       	/* lower right */
} *gmatrix;

typedef struct
{
	double				re;
	double				im;
} complex;


/**************************************************************
 *
 * Function Prototypes
 *
 **************************************************************/

/**************************************************************
 *
 * Initialization and utility functions
 *
 **************************************************************/

/* trig lookups. */
void			init_sinCos(int);
double			s_sin(int);
double			s_cos(int);


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

/* Creates a new giant matrix, but does not malloc
 * the component giants. */
gmatrix			newgmatrix(void);

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

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

/* Returns whether g is zero. */
int 			isZero(giant g);

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

/* Integer <-> giant. */
void 			itog(int n, giant g);
signed int		gtoi (giant);

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

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

/**************************************************************
 *
 * I/O Routines
 *
 **************************************************************/

/* Output the giant in decimal, with optional newlines.  */
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);

/**************************************************************
 *
 * Math Functions
 *
 **************************************************************/

/* g := -g. */
void 		negg(giant g);

/* g := |g|. */
void 		absg(giant g);

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

/* b += a. */
void 		addg(giant a, giant b);

/* b -= a. */
void 		subg(giant a, giant b);

/* Returns the number of trailing zero bits in g. */
int		numtrailzeros(giant g);

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

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

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

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

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

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

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

/* a becomes (a*b) (mod 2^q-k) where q % 16 == 0 and k is "small"
 * (0 < k < 65535). Returns 0 if unsuccessful, otherwise 1. */
int			feemulmod(giant x, giant y, int q, int k);

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

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

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

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

/* dest becomes lowermost n bits of src.
 * Equivalent to dest = src % 2^n. */
void 		extractbits(int n, giant src, giant dest);

/* negate g. g is mod 2^n+1. */
void 		fermatnegate(int n, giant g);

/* g := g (mod 2^n-1). */
void 		mersennemod(int n, giant g);

/* g := g (mod 2^n+1). */
void 		fermatmod(int n, giant g);

/* g := g (mod 2^n+1). */
void 		fer_mod(int n, giant g, giant modulus);

/* g *= s. */
void 		smulg(unsigned short s, giant g);

/* g *= g. */
void 		squareg(giant g);

/* b *= a. */
void 		mulg(giant a, giant b);

/* A giant gcd.  Modifies its arguments. */
void		ggcd(giant xx, giant yy);