1/* 2 Name: imath.h 3 Purpose: Arbitrary precision integer arithmetic routines. 4 Author: M. J. Fromberger 5 6 Copyright (C) 2002-2007 Michael J. Fromberger, All Rights Reserved. 7 8 Permission is hereby granted, free of charge, to any person obtaining a copy 9 of this software and associated documentation files (the "Software"), to deal 10 in the Software without restriction, including without limitation the rights 11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 12 copies of the Software, and to permit persons to whom the Software is 13 furnished to do so, subject to the following conditions: 14 15 The above copyright notice and this permission notice shall be included in 16 all copies or substantial portions of the Software. 17 18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 21 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 23 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 24 SOFTWARE. 25 */ 26 27#ifndef IMATH_H_ 28#define IMATH_H_ 29 30#include <limits.h> 31#include <stdbool.h> 32#include <stdint.h> 33 34#ifdef __cplusplus 35extern "C" { 36#endif 37 38typedef unsigned char mp_sign; 39typedef unsigned int mp_size; 40typedef int mp_result; 41typedef long mp_small; /* must be a signed type */ 42typedef unsigned long mp_usmall; /* must be an unsigned type */ 43 44 45/* Build with words as uint64_t by default. */ 46#ifdef USE_32BIT_WORDS 47typedef uint16_t mp_digit; 48typedef uint32_t mp_word; 49# define MP_DIGIT_MAX (UINT16_MAX * 1UL) 50# define MP_WORD_MAX (UINT32_MAX * 1UL) 51#else 52typedef uint32_t mp_digit; 53typedef uint64_t mp_word; 54# define MP_DIGIT_MAX (UINT32_MAX * UINT64_C(1)) 55# define MP_WORD_MAX (UINT64_MAX) 56#endif 57 58typedef struct { 59 mp_digit single; 60 mp_digit* digits; 61 mp_size alloc; 62 mp_size used; 63 mp_sign sign; 64} mpz_t, *mp_int; 65 66static inline mp_digit* MP_DIGITS(mp_int Z) { return Z->digits; } 67static inline mp_size MP_ALLOC(mp_int Z) { return Z->alloc; } 68static inline mp_size MP_USED(mp_int Z) { return Z->used; } 69static inline mp_sign MP_SIGN(mp_int Z) { return Z->sign; } 70 71extern const mp_result MP_OK; 72extern const mp_result MP_FALSE; 73extern const mp_result MP_TRUE; 74extern const mp_result MP_MEMORY; 75extern const mp_result MP_RANGE; 76extern const mp_result MP_UNDEF; 77extern const mp_result MP_TRUNC; 78extern const mp_result MP_BADARG; 79extern const mp_result MP_MINERR; 80 81#define MP_DIGIT_BIT (sizeof(mp_digit) * CHAR_BIT) 82#define MP_WORD_BIT (sizeof(mp_word) * CHAR_BIT) 83#define MP_SMALL_MIN LONG_MIN 84#define MP_SMALL_MAX LONG_MAX 85#define MP_USMALL_MAX ULONG_MAX 86 87#define MP_MIN_RADIX 2 88#define MP_MAX_RADIX 36 89 90/** Sets the default number of digits allocated to an `mp_int` constructed by 91 `mp_int_init_size()` with `prec == 0`. Allocations are rounded up to 92 multiples of this value. `MP_DEFAULT_PREC` is the default value. Requires 93 `ndigits > 0`. */ 94void mp_int_default_precision(mp_size ndigits); 95 96/** Sets the number of digits below which multiplication will use the standard 97 quadratic "schoolbook" multiplication algorithm rather than Karatsuba-Ofman. 98 Requires `ndigits >= sizeof(mp_word)`. */ 99void mp_int_multiply_threshold(mp_size ndigits); 100 101/** A sign indicating a (strictly) negative value. */ 102extern const mp_sign MP_NEG; 103 104/** A sign indicating a zero or positive value. */ 105extern const mp_sign MP_ZPOS; 106 107/** Reports whether `z` is odd, having remainder 1 when divided by 2. */ 108static inline bool mp_int_is_odd(mp_int z) { return (z->digits[0] & 1) != 0; } 109 110/** Reports whether `z` is even, having remainder 0 when divided by 2. */ 111static inline bool mp_int_is_even(mp_int z) { return (z->digits[0] & 1) == 0; } 112 113/** Initializes `z` with 1-digit precision and sets it to zero. This function 114 cannot fail unless `z == NULL`. */ 115mp_result mp_int_init(mp_int z); 116 117/** Allocates a fresh zero-valued `mpz_t` on the heap, returning NULL in case 118 of error. The only possible error is out-of-memory. */ 119mp_int mp_int_alloc(void); 120 121/** Initializes `z` with at least `prec` digits of storage, and sets it to 122 zero. If `prec` is zero, the default precision is used. In either case the 123 size is rounded up to the nearest multiple of the word size. */ 124mp_result mp_int_init_size(mp_int z, mp_size prec); 125 126/** Initializes `z` to be a copy of an already-initialized value in `old`. The 127 new copy does not share storage with the original. */ 128mp_result mp_int_init_copy(mp_int z, mp_int old); 129 130/** Initializes `z` to the specified signed `value` at default precision. */ 131mp_result mp_int_init_value(mp_int z, mp_small value); 132 133/** Initializes `z` to the specified unsigned `value` at default precision. */ 134mp_result mp_int_init_uvalue(mp_int z, mp_usmall uvalue); 135 136/** Sets `z` to the value of the specified signed `value`. */ 137mp_result mp_int_set_value(mp_int z, mp_small value); 138 139/** Sets `z` to the value of the specified unsigned `value`. */ 140mp_result mp_int_set_uvalue(mp_int z, mp_usmall uvalue); 141 142/** Releases the storage used by `z`. */ 143void mp_int_clear(mp_int z); 144 145/** Releases the storage used by `z` and also `z` itself. 146 This should only be used for `z` allocated by `mp_int_alloc()`. */ 147void mp_int_free(mp_int z); 148 149/** Replaces the value of `c` with a copy of the value of `a`. No new memory is 150 allocated unless `a` has more significant digits than `c` has allocated. */ 151mp_result mp_int_copy(mp_int a, mp_int c); 152 153/** Swaps the values and storage between `a` and `c`. */ 154void mp_int_swap(mp_int a, mp_int c); 155 156/** Sets `z` to zero. The allocated storage of `z` is not changed. */ 157void mp_int_zero(mp_int z); 158 159/** Sets `c` to the absolute value of `a`. */ 160mp_result mp_int_abs(mp_int a, mp_int c); 161 162/** Sets `c` to the additive inverse (negation) of `a`. */ 163mp_result mp_int_neg(mp_int a, mp_int c); 164 165/** Sets `c` to the sum of `a` and `b`. */ 166mp_result mp_int_add(mp_int a, mp_int b, mp_int c); 167 168/** Sets `c` to the sum of `a` and `value`. */ 169mp_result mp_int_add_value(mp_int a, mp_small value, mp_int c); 170 171/** Sets `c` to the difference of `a` less `b`. */ 172mp_result mp_int_sub(mp_int a, mp_int b, mp_int c); 173 174/** Sets `c` to the difference of `a` less `value`. */ 175mp_result mp_int_sub_value(mp_int a, mp_small value, mp_int c); 176 177/** Sets `c` to the product of `a` and `b`. */ 178mp_result mp_int_mul(mp_int a, mp_int b, mp_int c); 179 180/** Sets `c` to the product of `a` and `value`. */ 181mp_result mp_int_mul_value(mp_int a, mp_small value, mp_int c); 182 183/** Sets `c` to the product of `a` and `2^p2`. Requires `p2 >= 0`. */ 184mp_result mp_int_mul_pow2(mp_int a, mp_small p2, mp_int c); 185 186/** Sets `c` to the square of `a`. */ 187mp_result mp_int_sqr(mp_int a, mp_int c); 188 189/** Sets `q` and `r` to the quotent and remainder of `a / b`. Division by 190 powers of 2 is detected and handled efficiently. The remainder is pinned 191 to `0 <= r < b`. 192 193 Either of `q` or `r` may be NULL, but not both, and `q` and `r` may not 194 point to the same value. */ 195mp_result mp_int_div(mp_int a, mp_int b, mp_int q, mp_int r); 196 197/** Sets `q` and `*r` to the quotent and remainder of `a / value`. Division by 198 powers of 2 is detected and handled efficiently. The remainder is pinned to 199 `0 <= *r < b`. Either of `q` or `r` may be NULL. */ 200mp_result mp_int_div_value(mp_int a, mp_small value, mp_int q, mp_small *r); 201 202/** Sets `q` and `r` to the quotient and remainder of `a / 2^p2`. This is a 203 special case for division by powers of two that is more efficient than 204 using ordinary division. Note that `mp_int_div()` will automatically handle 205 this case, this function is for cases where you have only the exponent. */ 206mp_result mp_int_div_pow2(mp_int a, mp_small p2, mp_int q, mp_int r); 207 208/** Sets `c` to the remainder of `a / m`. 209 The remainder is pinned to `0 <= c < m`. */ 210mp_result mp_int_mod(mp_int a, mp_int m, mp_int c); 211 212/** Sets `c` to the value of `a` raised to the `b` power. 213 It returns `MP_RANGE` if `b < 0`. */ 214mp_result mp_int_expt(mp_int a, mp_small b, mp_int c); 215 216/** Sets `c` to the value of `a` raised to the `b` power. 217 It returns `MP_RANGE` if `b < 0`. */ 218mp_result mp_int_expt_value(mp_small a, mp_small b, mp_int c); 219 220/** Sets `c` to the value of `a` raised to the `b` power. 221 It returns `MP_RANGE`) if `b < 0`. */ 222mp_result mp_int_expt_full(mp_int a, mp_int b, mp_int c); 223 224/** Sets `*r` to the remainder of `a / value`. 225 The remainder is pinned to `0 <= r < value`. */ 226static inline 227mp_result mp_int_mod_value(mp_int a, mp_small value, mp_small* r) { 228 return mp_int_div_value(a, value, 0, r); 229} 230 231/** Returns the comparator of `a` and `b`. */ 232int mp_int_compare(mp_int a, mp_int b); 233 234/** Returns the comparator of the magnitudes of `a` and `b`, disregarding their 235 signs. Neither `a` nor `b` is modified by the comparison. */ 236int mp_int_compare_unsigned(mp_int a, mp_int b); 237 238/** Returns the comparator of `z` and zero. */ 239int mp_int_compare_zero(mp_int z); 240 241/** Returns the comparator of `z` and the signed value `v`. */ 242int mp_int_compare_value(mp_int z, mp_small v); 243 244/** Returns the comparator of `z` and the unsigned value `uv`. */ 245int mp_int_compare_uvalue(mp_int z, mp_usmall uv); 246 247/** Reports whether `a` is divisible by `v`. */ 248bool mp_int_divisible_value(mp_int a, mp_small v); 249 250/** Returns `k >= 0` such that `z` is `2^k`, if such a `k` exists. If no such 251 `k` exists, the function returns -1. */ 252int mp_int_is_pow2(mp_int z); 253 254/** Sets `c` to the value of `a` raised to the `b` power, reduced modulo `m`. 255 It returns `MP_RANGE` if `b < 0` or `MP_UNDEF` if `m == 0`. */ 256mp_result mp_int_exptmod(mp_int a, mp_int b, mp_int m, mp_int c); 257 258/** Sets `c` to the value of `a` raised to the `value` power, modulo `m`. 259 It returns `MP_RANGE` if `value < 0` or `MP_UNDEF` if `m == 0`. */ 260mp_result mp_int_exptmod_evalue(mp_int a, mp_small value, mp_int m, mp_int c); 261 262/** Sets `c` to the value of `value` raised to the `b` power, modulo `m`. 263 It returns `MP_RANGE` if `b < 0` or `MP_UNDEF` if `m == 0`. */ 264mp_result mp_int_exptmod_bvalue(mp_small value, mp_int b, mp_int m, mp_int c); 265 266/** Sets `c` to the value of `a` raised to the `b` power, reduced modulo `m`, 267 given a precomputed reduction constant `mu` defined for Barrett's modular 268 reduction algorithm. 269 270 It returns `MP_RANGE` if `b < 0` or `MP_UNDEF` if `m == 0`. */ 271mp_result mp_int_exptmod_known(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c); 272 273/** Sets `c` to the reduction constant for Barrett reduction by modulus `m`. 274 Requires that `c` and `m` point to distinct locations. */ 275mp_result mp_int_redux_const(mp_int m, mp_int c); 276 277/** Sets `c` to the multiplicative inverse of `a` modulo `m`, if it exists. 278 The least non-negative representative of the congruence class is computed. 279 280 It returns `MP_UNDEF` if the inverse does not exist, or `MP_RANGE` if `a == 281 0` or `m <= 0`. */ 282mp_result mp_int_invmod(mp_int a, mp_int m, mp_int c); 283 284/** Sets `c` to the greatest common divisor of `a` and `b`. 285 286 It returns `MP_UNDEF` if the GCD is undefined, such as for example if `a` 287 and `b` are both zero. */ 288mp_result mp_int_gcd(mp_int a, mp_int b, mp_int c); 289 290/** Sets `c` to the greatest common divisor of `a` and `b`, and sets `x` and 291 `y` to values satisfying Bezout's identity `gcd(a, b) = ax + by`. 292 293 It returns `MP_UNDEF` if the GCD is undefined, such as for example if `a` 294 and `b` are both zero. */ 295mp_result mp_int_egcd(mp_int a, mp_int b, mp_int c, mp_int x, mp_int y); 296 297/** Sets `c` to the least common multiple of `a` and `b`. 298 299 It returns `MP_UNDEF` if the LCM is undefined, such as for example if `a` 300 and `b` are both zero. */ 301mp_result mp_int_lcm(mp_int a, mp_int b, mp_int c); 302 303/** Sets `c` to the greatest integer not less than the `b`th root of `a`, 304 using Newton's root-finding algorithm. 305 It returns `MP_UNDEF` if `a < 0` and `b` is even. */ 306mp_result mp_int_root(mp_int a, mp_small b, mp_int c); 307 308/** Sets `c` to the greatest integer not less than the square root of `a`. 309 This is a special case of `mp_int_root()`. */ 310static inline 311mp_result mp_int_sqrt(mp_int a, mp_int c) { return mp_int_root(a, 2, c); } 312 313/** Returns `MP_OK` if `z` is representable as `mp_small`, else `MP_RANGE`. 314 If `out` is not NULL, `*out` is set to the value of `z` when `MP_OK`. */ 315mp_result mp_int_to_int(mp_int z, mp_small *out); 316 317/** Returns `MP_OK` if `z` is representable as `mp_usmall`, or `MP_RANGE`. 318 If `out` is not NULL, `*out` is set to the value of `z` when `MP_OK`. */ 319mp_result mp_int_to_uint(mp_int z, mp_usmall *out); 320 321/** Converts `z` to a zero-terminated string of characters in the specified 322 `radix`, writing at most `limit` characters to `str` including the 323 terminating NUL value. A leading `-` is used to indicate a negative value. 324 325 Returns `MP_TRUNC` if `limit` was to small to write all of `z`. 326 Requires `MP_MIN_RADIX <= radix <= MP_MAX_RADIX`. */ 327mp_result mp_int_to_string(mp_int z, mp_size radix, char *str, int limit); 328 329/** Reports the minimum number of characters required to represent `z` as a 330 zero-terminated string in the given `radix`. 331 Requires `MP_MIN_RADIX <= radix <= MP_MAX_RADIX`. */ 332mp_result mp_int_string_len(mp_int z, mp_size radix); 333 334/** Reads a string of ASCII digits in the specified `radix` from the zero 335 terminated `str` provided into `z`. For values of `radix > 10`, the letters 336 `A`..`Z` or `a`..`z` are accepted. Letters are interpreted without respect 337 to case. 338 339 Leading whitespace is ignored, and a leading `+` or `-` is interpreted as a 340 sign flag. Processing stops when a NUL or any other character out of range 341 for a digit in the given radix is encountered. 342 343 If the whole string was consumed, `MP_OK` is returned; otherwise 344 `MP_TRUNC`. is returned. 345 346 Requires `MP_MIN_RADIX <= radix <= MP_MAX_RADIX`. */ 347mp_result mp_int_read_string(mp_int z, mp_size radix, const char *str); 348 349/** Reads a string of ASCII digits in the specified `radix` from the zero 350 terminated `str` provided into `z`. For values of `radix > 10`, the letters 351 `A`..`Z` or `a`..`z` are accepted. Letters are interpreted without respect 352 to case. 353 354 Leading whitespace is ignored, and a leading `+` or `-` is interpreted as a 355 sign flag. Processing stops when a NUL or any other character out of range 356 for a digit in the given radix is encountered. 357 358 If the whole string was consumed, `MP_OK` is returned; otherwise 359 `MP_TRUNC`. is returned. If `end` is not NULL, `*end` is set to point to 360 the first unconsumed byte of the input string (the NUL byte if the whole 361 string was consumed). This emulates the behavior of the standard C 362 `strtol()` function. 363 364 Requires `MP_MIN_RADIX <= radix <= MP_MAX_RADIX`. */ 365mp_result mp_int_read_cstring(mp_int z, mp_size radix, const char *str, char **end); 366 367/** Returns the number of significant bits in `z`. */ 368mp_result mp_int_count_bits(mp_int z); 369 370/** Converts `z` to 2's complement binary, writing at most `limit` bytes into 371 the given `buf`. Returns `MP_TRUNC` if the buffer limit was too small to 372 contain the whole value. If this occurs, the contents of buf will be 373 effectively garbage, as the function uses the buffer as scratch space. 374 375 The binary representation of `z` is in base-256 with digits ordered from 376 most significant to least significant (network byte ordering). The 377 high-order bit of the first byte is set for negative values, clear for 378 non-negative values. 379 380 As a result, non-negative values will be padded with a leading zero byte if 381 the high-order byte of the base-256 magnitude is set. This extra byte is 382 accounted for by the `mp_int_binary_len()` function. */ 383mp_result mp_int_to_binary(mp_int z, unsigned char *buf, int limit); 384 385/** Reads a 2's complement binary value from `buf` into `z`, where `len` is the 386 length of the buffer. The contents of `buf` may be overwritten during 387 processing, although they will be restored when the function returns. */ 388mp_result mp_int_read_binary(mp_int z, unsigned char *buf, int len); 389 390/** Returns the number of bytes to represent `z` in 2's complement binary. */ 391mp_result mp_int_binary_len(mp_int z); 392 393/** Converts the magnitude of `z` to unsigned binary, writing at most `limit` 394 bytes into the given `buf`. The sign of `z` is ignored, but `z` is not 395 modified. Returns `MP_TRUNC` if the buffer limit was too small to contain 396 the whole value. If this occurs, the contents of `buf` will be effectively 397 garbage, as the function uses the buffer as scratch space during 398 conversion. 399 400 The binary representation of `z` is in base-256 with digits ordered from 401 most significant to least significant (network byte ordering). */ 402mp_result mp_int_to_unsigned(mp_int z, unsigned char *buf, int limit); 403 404/** Reads an unsigned binary value from `buf` into `z`, where `len` is the 405 length of the buffer. The contents of `buf` are not modified during 406 processing. */ 407mp_result mp_int_read_unsigned(mp_int z, unsigned char *buf, int len); 408 409/** Returns the number of bytes required to represent `z` as an unsigned binary 410 value in base 256. */ 411mp_result mp_int_unsigned_len(mp_int z); 412 413/** Returns a pointer to a brief, human-readable, zero-terminated string 414 describing `res`. The returned string is statically allocated and must not 415 be freed by the caller. */ 416const char *mp_error_string(mp_result res); 417 418#ifdef __cplusplus 419} 420#endif 421#endif /* end IMATH_H_ */ 422