1/************************************************************** 2 * 3 * giants.h 4 * 5 * Header file for large-integer arithmetic library giants.c. 6 * 7 * Updates: 8 * 18 Jul 99 REC Added fer_mod(). 9 * 30 Apr 98 JF USE_ASSEMBLER_MUL removed 10 * 29 Apr 98 JF Function prototypes cleaned up 11 * 20 Apr 97 RDW 12 * 13 * c. 1997 Perfectly Scientific, Inc. 14 * All Rights Reserved. 15 * 16 **************************************************************/ 17 18 19/************************************************************** 20 * 21 * Error Codes 22 * 23 **************************************************************/ 24 25#define DIVIDEBYZERO 1 26#define OVFLOW 2 27#define SIGN 3 28#define OVERRANGE 4 29#define AUTO_MUL 0 30#define GRAMMAR_MUL 1 31#define FFT_MUL 2 32#define KARAT_MUL 3 33 34/************************************************************** 35 * 36 * Preprocessor definitions 37 * 38 **************************************************************/ 39 40/* 2^(16*MAX_SHORTS)-1 will fit into a giant, but take care: 41 * one usually has squares, etc. of giants involved, and 42 * every intermediate giant in a calculation must fit into 43 * this many shorts. Thus, if you want systematically to effect 44 * arithmetic on B-bit operands, you need MAX_SHORTS > B/8, 45 * perferably a tad larger than this; e.g. MAX_SHORTS > B/7. 46 */ 47#define MAX_SHORTS (1<<19) 48 49#define INFINITY (-1) 50#define FA 0 51#define TR 1 52#define COLUMNWIDTH 64 53 54#define TWOPI (double)(2*3.1415926535897932384626433) 55#define SQRT2 (double)(1.414213562373095048801688724209) 56#define SQRTHALF (double)(0.707106781186547524400844362104) 57#define TWO16 (double)(65536.0) 58#define TWOM16 (double)(0.0000152587890625) 59 60/* Decimal digit ceiling in digit-input routines. */ 61#define MAX_DIGITS 10000 62 63/* Next, mumber of shorts per operand 64 at which Karatsuba breaks over. */ 65#define KARAT_BREAK 40 66 67/* Next, mumber of shorts per operand 68 at which FFT breaks over. */ 69#define FFT_BREAK 200 70 71#define newmin(a,b) ((a)<(b)? (a) : (b)) 72#define newmax(a,b) ((a)>(b)? (a) : (b)) 73 74/* Maximum number of recursive steps needed to calculate 75 * gcds of integers. */ 76#define STEPS 32 77 78/* The limit below which hgcd is too ponderous */ 79#define GCDLIMIT 5000 80 81/* The limit below which ordinary ints will be used */ 82#define INTLIMIT 31 83 84/* Size by which to increment the stack used in pushg() and popg(). */ 85#define STACK_GROW 16 86 87#define gin(x) gread(x,stdin) 88#define gout(x) gwriteln(x,stdout) 89 90 91/************************************************************** 92 * 93 * Structure definitions 94 * 95 **************************************************************/ 96 97typedef struct 98{ 99 int sign; 100 unsigned short n[1]; /* number of shorts = abs(sign) */ 101} giantstruct; 102 103typedef giantstruct *giant; 104 105typedef struct _matrix 106{ 107 giant ul; /* upper left */ 108 giant ur; /* upper right */ 109 giant ll; /* lower left */ 110 giant lr; /* lower right */ 111} *gmatrix; 112 113typedef struct 114{ 115 double re; 116 double im; 117} complex; 118 119 120/************************************************************** 121 * 122 * Function Prototypes 123 * 124 **************************************************************/ 125 126/************************************************************** 127 * 128 * Initialization and utility functions 129 * 130 **************************************************************/ 131 132/* trig lookups. */ 133void init_sinCos(int); 134double s_sin(int); 135double s_cos(int); 136 137 138/* Creates a new giant, numshorts = INFINITY invokes the 139 * maximum MAX_SHORTS. */ 140giant newgiant(int numshorts); 141 142/* Creates a new giant matrix, but does not malloc 143 * the component giants. */ 144gmatrix newgmatrix(void); 145 146/* Returns the bit-length n; e.g. n=7 returns 3. */ 147int bitlen(giant n); 148 149/* Returns the value of the pos bit of n. */ 150int bitval(giant n, int pos); 151 152/* Returns whether g is one. */ 153int isone(giant g); 154 155/* Returns whether g is zero. */ 156int isZero(giant g); 157 158/* Copies one giant to another. */ 159void gtog(giant src, giant dest); 160 161/* Integer <-> giant. */ 162void itog(int n, giant g); 163signed int gtoi (giant); 164 165/* Returns the sign of g: -1, 0, 1. */ 166int gsign(giant g); 167 168/* Returns 1, 0, -1 as a>b, a=b, a<b. */ 169int gcompg(giant a, giant b); 170 171/* Set AUTO_MUL for automatic FFT crossover (this is the 172 * default), set FFT_MUL for forced FFT multiply, set 173 * GRAMMAR_MUL for forced grammar school multiply. */ 174void setmulmode(int mode); 175 176/************************************************************** 177 * 178 * I/O Routines 179 * 180 **************************************************************/ 181 182/* Output the giant in decimal, with optional newlines. */ 183void gwrite(giant g, FILE *fp, int newlines); 184 185/* Output the giant in decimal, with both '\'-newline 186 * notation and a final newline. */ 187void gwriteln(giant g, FILE *fp); 188 189/* Input the giant in decimal, assuming the formatting of 190 * 'gwriteln'. */ 191void gread(giant g, FILE *fp); 192 193/************************************************************** 194 * 195 * Math Functions 196 * 197 **************************************************************/ 198 199/* g := -g. */ 200void negg(giant g); 201 202/* g := |g|. */ 203void absg(giant g); 204 205/* g += i, with i non-negative and < 2^16. */ 206void iaddg(int i,giant g); 207 208/* b += a. */ 209void addg(giant a, giant b); 210 211/* b -= a. */ 212void subg(giant a, giant b); 213 214/* Returns the number of trailing zero bits in g. */ 215int numtrailzeros(giant g); 216 217/* u becomes greatest power of two not exceeding u/v. */ 218void bdivg(giant v, giant u); 219 220/* Same as invg, but uses bdivg. */ 221int binvg(giant n, giant x); 222 223/* If 1/x exists (mod n), 1 is returned and x := 1/x. If 224 * inverse does not exist, 0 is returned and x := GCD(n, x). */ 225int invg(giant n, giant x); 226 227int mersenneinvg(int q, giant x); 228 229/* Classical GCD, x:= GCD(n, x). */ 230void cgcdg(giant n, giant x); 231 232/* General GCD, x:= GCD(n, x). */ 233void gcdg(giant n, giant x); 234 235/* Binary GCD, x:= GCD(n, x). */ 236void bgcdg(giant n, giant x); 237 238/* g := m^n, no mod is performed. */ 239void powerg(int a, int b, giant g); 240 241/* r becomes the steady-state reciprocal 2^(2b)/d, where 242 * b = bit-length of d-1. */ 243void make_recip(giant d, giant r); 244 245/* n := [n/d], d positive, using stored reciprocal directly. */ 246void divg_via_recip(giant d, giant r, giant n); 247 248/* n := n % d, d positive, using stored reciprocal directly. */ 249void modg_via_recip(giant d, giant r, giant n); 250 251/* num := num % den, any positive den. */ 252void modg(giant den, giant num); 253 254/* a becomes (a*b) (mod 2^q-k) where q % 16 == 0 and k is "small" 255 * (0 < k < 65535). Returns 0 if unsuccessful, otherwise 1. */ 256int feemulmod(giant x, giant y, int q, int k); 257 258/* g := g/n, and (g mod n) is returned. */ 259int idivg(int n, giant g); 260 261/* num := [num/den], any positive den. */ 262void divg(giant den, giant num); 263 264/* num := [num/den], any positive den. */ 265void powermod(giant x, int n, giant z); 266 267/* x := x^n (mod z). */ 268void powermodg(giant x, giant n, giant z); 269 270/* x := x^n (mod 2^q+1). */ 271void fermatpowermod(giant x, int n, int q); 272 273/* x := x^n (mod 2^q+1). */ 274void fermatpowermodg(giant x, giant n, int q); 275 276/* x := x^n (mod 2^q-1). */ 277void mersennepowermod(giant x, int n, int q); 278 279/* x := x^n (mod 2^q-1). */ 280void mersennepowermodg(giant x, giant n, int q); 281 282/* Shift g left by bits, introducing zeros on the right. */ 283void gshiftleft(int bits, giant g); 284 285/* Shift g right by bits, losing bits on the right. */ 286void gshiftright(int bits, giant g); 287 288/* dest becomes lowermost n bits of src. 289 * Equivalent to dest = src % 2^n. */ 290void extractbits(int n, giant src, giant dest); 291 292/* negate g. g is mod 2^n+1. */ 293void fermatnegate(int n, giant g); 294 295/* g := g (mod 2^n-1). */ 296void mersennemod(int n, giant g); 297 298/* g := g (mod 2^n+1). */ 299void fermatmod(int n, giant g); 300 301/* g := g (mod 2^n+1). */ 302void fer_mod(int n, giant g, giant modulus); 303 304/* g *= s. */ 305void smulg(unsigned short s, giant g); 306 307/* g *= g. */ 308void squareg(giant g); 309 310/* b *= a. */ 311void mulg(giant a, giant b); 312 313/* A giant gcd. Modifies its arguments. */ 314void ggcd(giant xx, giant yy); 315