--- /dev/null
+/**************************************************************
+ *
+ * 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 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 = 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);