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