imath.h (11559B)
1 /* 2 ** imath.h 3 ** Arbitrary precision integer arithmetic routines. 4 ** See Copyright Notice in klisp.h 5 */ 6 7 /* 8 ** SOURCE NOTE: This is mostly from the IMath library, written by 9 ** M.J. Fromberger. It is adapted to klisp, mainly in the use of the 10 ** klisp allocator and fixing of digit size to 32 bits. 11 ** Imported from version (1.15) updated 01-Feb-2011 at 03:10 PM. 12 */ 13 14 #ifndef IMATH_H_ 15 #define IMATH_H_ 16 17 #include <limits.h> 18 19 /* Andres Navarro: use c99 constants */ 20 #ifndef USE_C99 21 #define USE_C99 1 22 #endif 23 24 /* Andres Navarro: klisp includes */ 25 #include "kobject.h" 26 #include "kstate.h" 27 28 #ifdef USE_C99 29 #include <stdint.h> 30 #endif 31 32 #ifdef __cplusplus 33 extern "C" { 34 #endif 35 36 #if USE_C99 37 typedef unsigned char mp_sign; 38 typedef uint32_t mp_size; 39 typedef int mp_result; 40 typedef int32_t mp_small; /* must be a signed type */ 41 typedef uint32_t mp_usmall; /* must be an unsigned type */ 42 typedef uint32_t mp_digit; 43 typedef uint64_t mp_word; 44 #else /* USE_C99 */ 45 typedef unsigned char mp_sign; 46 typedef unsigned int mp_size; 47 typedef int mp_result; 48 typedef long mp_small; /* must be a signed type */ 49 typedef unsigned long mp_usmall; /* must be an unsigned type */ 50 #ifdef USE_LONG_LONG 51 typedef unsigned int mp_digit; 52 typedef unsigned long long mp_word; 53 #else /* USE_LONG_LONG */ 54 typedef unsigned short mp_digit; 55 typedef unsigned int mp_word; 56 #endif /* USE_LONG_LONG */ 57 #endif /* USE_C99 */ 58 59 /* Andres Navarro: Use kobject type instead */ 60 typedef Bigint mpz_t, *mp_int; 61 62 #if 0 63 typedef struct mpz { 64 mp_digit single; 65 mp_digit *digits; 66 mp_size alloc; 67 mp_size used; 68 mp_sign sign; 69 } mpz_t, *mp_int; 70 #endif 71 72 #define MP_SINGLE(Z) ((Z)->single) /* added to correct check in mp_int_clear */ 73 #define MP_DIGITS(Z) ((Z)->digits) 74 #define MP_ALLOC(Z) ((Z)->alloc) 75 #define MP_USED(Z) ((Z)->used) 76 #define MP_SIGN(Z) ((Z)->sign) 77 78 extern const mp_result MP_OK; 79 extern const mp_result MP_FALSE; 80 extern const mp_result MP_TRUE; 81 extern const mp_result MP_MEMORY; 82 extern const mp_result MP_RANGE; 83 extern const mp_result MP_UNDEF; 84 extern const mp_result MP_TRUNC; 85 extern const mp_result MP_BADARG; 86 extern const mp_result MP_MINERR; 87 88 #define MP_DIGIT_BIT (sizeof(mp_digit) * CHAR_BIT) 89 #define MP_WORD_BIT (sizeof(mp_word) * CHAR_BIT) 90 /* Andres Navarro: USE_C99 */ 91 #ifdef USE_C99 92 #define MP_SMALL_MIN INT32_MIN 93 #define MP_SMALL_MAX INT32_MAX 94 #define MP_USMALL_MIN UINT32_MIN 95 #define MP_USMALL_MAX UINT32_MAX 96 #define MP_DIGIT_MAX ((uint64_t) UINT32_MAX) 97 #define MP_WORD_MAX UINT64_MAX 98 #else /* USE_C99 */ 99 #define MP_SMALL_MIN LONG_MIN 100 #define MP_SMALL_MAX LONG_MAX 101 #define MP_USMALL_MIN ULONG_MIN 102 #define MP_USMALL_MAX ULONG_MAX 103 #ifdef USE_LONG_LONG 104 # ifndef ULONG_LONG_MAX 105 # ifdef ULLONG_MAX 106 # define ULONG_LONG_MAX ULLONG_MAX 107 # else 108 # error "Maximum value of unsigned long long not defined!" 109 # endif 110 # endif 111 # define MP_DIGIT_MAX (UINT_MAX * 1ULL) 112 # define MP_WORD_MAX ULONG_LONG_MAX 113 #else /* USE_LONG_LONG */ 114 # define MP_DIGIT_MAX (USHRT_MAX * 1UL) 115 # define MP_WORD_MAX (UINT_MAX * 1UL) 116 #endif /* USE_LONG_LONG */ 117 #endif /* USE_C99 */ 118 119 #define MP_MIN_RADIX 2 120 #define MP_MAX_RADIX 36 121 122 /* Values with fewer than this many significant digits use the 123 standard multiplication algorithm; otherwise, a recursive algorithm 124 is used. Choose a value to suit your platform. 125 */ 126 #define MP_MULT_THRESH 22 127 128 #define MP_DEFAULT_PREC 8 /* default memory allocation, in digits */ 129 130 extern const mp_sign MP_NEG; 131 extern const mp_sign MP_ZPOS; 132 133 #define mp_int_is_odd(Z) ((Z)->digits[0] & 1) 134 #define mp_int_is_even(Z) !((Z)->digits[0] & 1) 135 136 /* NOTE: this doesn't use the allocator */ 137 mp_result mp_int_init(mp_int z); 138 mp_int mp_int_alloc(klisp_State *K); 139 mp_result mp_int_init_size(klisp_State *K, mp_int z, mp_size prec); 140 mp_result mp_int_init_copy(klisp_State *K, mp_int z, mp_int old); 141 mp_result mp_int_init_value(klisp_State *K, mp_int z, mp_small value); 142 mp_result mp_int_set_value(klisp_State *K, mp_int z, mp_small value); 143 void mp_int_clear(klisp_State *K, mp_int z); 144 void mp_int_free(klisp_State *K, mp_int z); 145 146 mp_result mp_int_copy(klisp_State *K, mp_int a, mp_int c); /* c = a */ 147 /* NOTE: this doesn't use the allocator */ 148 void mp_int_swap(mp_int a, mp_int c); /* swap a, c */ 149 /* NOTE: this doesn't use the allocator */ 150 void mp_int_zero(mp_int z); /* z = 0 */ 151 mp_result mp_int_abs(klisp_State *K, mp_int a, mp_int c); /* c = |a| */ 152 mp_result mp_int_neg(klisp_State *K, mp_int a, mp_int c); /* c = -a */ 153 /* c = a + b */ 154 mp_result mp_int_add(klisp_State *K, mp_int a, mp_int b, mp_int c); 155 mp_result mp_int_add_value(klisp_State *K, mp_int a, mp_small value, 156 mp_int c); 157 /* c = a - b */ 158 mp_result mp_int_sub(klisp_State *K, mp_int a, mp_int b, mp_int c); 159 mp_result mp_int_sub_value(klisp_State *K, mp_int a, mp_small value, 160 mp_int c); 161 /* c = a * b */ 162 mp_result mp_int_mul(klisp_State *K, mp_int a, mp_int b, mp_int c); 163 mp_result mp_int_mul_value(klisp_State *K, mp_int a, mp_small value, 164 mp_int c); 165 mp_result mp_int_mul_pow2(klisp_State *K, mp_int a, mp_small p2, mp_int c); 166 mp_result mp_int_sqr(klisp_State *K, mp_int a, mp_int c); /* c = a * a */ 167 /* q = a / b */ 168 /* r = a % b */ 169 mp_result mp_int_div(klisp_State *K, mp_int a, mp_int b, mp_int q, 170 mp_int r); 171 /* q = a / value */ 172 /* r = a % value */ 173 mp_result mp_int_div_value(klisp_State *K, mp_int a, mp_small value, 174 mp_int q, mp_small *r); 175 /* q = a / 2^p2 */ 176 /* r = q % 2^p2 */ 177 mp_result mp_int_div_pow2(klisp_State *K, mp_int a, mp_small p2, 178 mp_int q, mp_int r); 179 /* c = a % m */ 180 mp_result mp_int_mod(klisp_State *K, mp_int a, mp_int m, mp_int c); 181 #define mp_int_mod_value(K, A, V, R) \ 182 mp_int_div_value((K), (A), (V), 0, (R)) 183 /* c = a^b */ 184 mp_result mp_int_expt(klisp_State *K, mp_int a, mp_small b, mp_int c); 185 /* c = a^b */ 186 mp_result mp_int_expt_value(klisp_State *K, mp_small a, mp_small b, mp_int c); 187 /* c = a^b */ 188 mp_result mp_int_expt_full(klisp_State *K, mp_int a, mp_int b, mp_int c); 189 190 /* NOTE: this doesn't use the allocator */ 191 int mp_int_compare(mp_int a, mp_int b); /* a <=> b */ 192 /* NOTE: this doesn't use the allocator */ 193 int mp_int_compare_unsigned(mp_int a, mp_int b); /* |a| <=> |b| */ 194 /* NOTE: this doesn't use the allocator */ 195 int mp_int_compare_zero(mp_int z); /* a <=> 0 */ 196 /* NOTE: this doesn't use the allocator */ 197 int mp_int_compare_value(mp_int z, mp_small value); /* a <=> v */ 198 199 /* Returns true if v|a, false otherwise (including errors) */ 200 int mp_int_divisible_value(klisp_State *K, mp_int a, mp_small v); 201 202 /* NOTE: this doesn't use the allocator */ 203 /* Returns k >= 0 such that z = 2^k, if one exists; otherwise < 0 */ 204 int mp_int_is_pow2(mp_int z); 205 206 mp_result mp_int_exptmod(klisp_State *K, mp_int a, mp_int b, mp_int m, 207 mp_int c); /* c = a^b (mod m) */ 208 mp_result mp_int_exptmod_evalue(klisp_State *K, mp_int a, mp_small value, 209 mp_int m, mp_int c); /* c = a^v (mod m) */ 210 mp_result mp_int_exptmod_bvalue(klisp_State *K, mp_small value, mp_int b, 211 mp_int m, mp_int c); /* c = v^b (mod m) */ 212 mp_result mp_int_exptmod_known(klisp_State *K, mp_int a, mp_int b, 213 mp_int m, mp_int mu, 214 mp_int c); /* c = a^b (mod m) */ 215 mp_result mp_int_redux_const(klisp_State *K, mp_int m, mp_int c); 216 217 /* c = 1/a (mod m) */ 218 mp_result mp_int_invmod(klisp_State *K, mp_int a, mp_int m, mp_int c); 219 220 /* c = gcd(a, b) */ 221 mp_result mp_int_gcd(klisp_State *K, mp_int a, mp_int b, mp_int c); 222 223 /* c = gcd(a, b) */ 224 /* c = ax + by */ 225 mp_result mp_int_egcd(klisp_State *K, mp_int a, mp_int b, mp_int c, 226 mp_int x, mp_int y); 227 228 /* c = lcm(a, b) */ 229 mp_result mp_int_lcm(klisp_State *K, mp_int a, mp_int b, mp_int c); 230 231 /* c = floor(a^{1/b}) */ 232 mp_result mp_int_root(klisp_State *K, mp_int a, mp_small b, mp_int c); 233 /* c = floor(sqrt(a)) */ 234 #define mp_int_sqrt(K, a, c) mp_int_root((K), a, 2, c) 235 236 /* Convert to a small int, if representable; else MP_RANGE */ 237 /* NOTE: this doesn't use the allocator */ 238 mp_result mp_int_to_int(mp_int z, mp_small *out); 239 /* NOTE: this doesn't use the allocator */ 240 mp_result mp_int_to_uint(mp_int z, mp_usmall *out); 241 242 /* Convert to nul-terminated string with the specified radix, writing at 243 most limit characters including the nul terminator */ 244 mp_result mp_int_to_string(klisp_State *K, mp_int z, mp_size radix, 245 char *str, int limit); 246 247 /* Return the number of characters required to represent 248 z in the given radix. May over-estimate. */ 249 /* NOTE: this doesn't use the allocator */ 250 mp_result mp_int_string_len(mp_int z, mp_size radix); 251 252 /* Read zero-terminated string into z */ 253 mp_result mp_int_read_string(klisp_State *K, mp_int z, mp_size radix, 254 const char *str); 255 mp_result mp_int_read_cstring(klisp_State *K, mp_int z, mp_size radix, 256 const char *str, char **end); 257 258 /* Return the number of significant bits in z */ 259 /* NOTE: this doesn't use the allocator */ 260 mp_result mp_int_count_bits(mp_int z); 261 262 /* Convert z to two's complement binary, writing at most limit bytes */ 263 mp_result mp_int_to_binary(klisp_State *K, mp_int z, unsigned char *buf, 264 int limit); 265 266 /* Read a two's complement binary value into z from the given buffer */ 267 mp_result mp_int_read_binary(klisp_State *K, mp_int z, unsigned char *buf, 268 int len); 269 270 /* Return the number of bytes required to represent z in binary. */ 271 /* NOTE: this doesn't use the allocator */ 272 mp_result mp_int_binary_len(mp_int z); 273 274 /* Convert z to unsigned binary, writing at most limit bytes */ 275 mp_result mp_int_to_unsigned(klisp_State *K, mp_int z, unsigned char *buf, 276 int limit); 277 278 /* Read an unsigned binary value into z from the given buffer */ 279 mp_result mp_int_read_unsigned(klisp_State *K, mp_int z, unsigned char *buf, 280 int len); 281 282 /* Return the number of bytes required to represent z as unsigned output */ 283 /* NOTE: this doesn't use the allocator */ 284 mp_result mp_int_unsigned_len(mp_int z); 285 286 /* Return a statically allocated string describing error code res */ 287 /* NOTE: this doesn't use the allocator */ 288 const char *mp_error_string(mp_result res); 289 290 #if DEBUG 291 void s_print(char *tag, mp_int z); 292 void s_print_buf(char *tag, mp_digit *buf, mp_size num); 293 #endif 294 295 #ifdef __cplusplus 296 } 297 #endif 298 299 #endif /* end IMATH_H_ */