polly/lib/External/isl/imath/imath.h

   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
  35 extern "C" {
  36 #endif
  37
  38 typedef unsigned char  mp_sign;
  39 typedef unsigned int   mp_size;
  40 typedef int            mp_result;
  41 typedef long           mp_small;  /* must be a signed type */
  42 typedef 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
  47 typedef uint16_t        mp_digit;
  48 typedef uint32_t        mp_word;
  49 #  define MP_DIGIT_MAX  (UINT16_MAX * 1UL)
  50 #  define MP_WORD_MAX   (UINT32_MAX * 1UL)
  51 #else
  52 typedef uint32_t        mp_digit;
  53 typedef uint64_t        mp_word;
  54 #  define MP_DIGIT_MAX  (UINT32_MAX * UINT64_C(1))
  55 #  define MP_WORD_MAX   (UINT64_MAX)
  56 #endif
  57
  58 typedef 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
  66 static inline mp_digit* MP_DIGITS(mp_int Z) { return Z->digits; }
  67 static inline mp_size   MP_ALLOC(mp_int Z)  { return Z->alloc; }
  68 static inline mp_size   MP_USED(mp_int Z)   { return Z->used; }
  69 static inline mp_sign   MP_SIGN(mp_int Z)   { return Z->sign; }
  70
  71 extern const mp_result MP_OK;
  72 extern const mp_result MP_FALSE;
  73 extern const mp_result MP_TRUE;
  74 extern const mp_result MP_MEMORY;
  75 extern const mp_result MP_RANGE;
  76 extern const mp_result MP_UNDEF;
  77 extern const mp_result MP_TRUNC;
  78 extern const mp_result MP_BADARG;
  79 extern 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`. */
  94 void 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)`. */
  99 void mp_int_multiply_threshold(mp_size ndigits);
 100
 101 /** A sign indicating a (strictly) negative value. */
 102 extern const mp_sign MP_NEG;
 103
 104 /** A sign indicating a zero or positive value. */
 105 extern const mp_sign MP_ZPOS;
 106
 107 /** Reports whether `z` is odd, having remainder 1 when divided by 2. */
 108 static 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. */
 111 static 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`. */
 115 mp_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. */
 119 mp_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. */
 124 mp_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. */
 128 mp_result mp_int_init_copy(mp_int z, mp_int old);
 129
 130 /** Initializes `z` to the specified signed `value` at default precision. */
 131 mp_result mp_int_init_value(mp_int z, mp_small value);
 132
 133 /** Initializes `z` to the specified unsigned `value` at default precision. */
 134 mp_result mp_int_init_uvalue(mp_int z, mp_usmall uvalue);
 135
 136 /** Sets `z` to the value of the specified signed `value`. */
 137 mp_result mp_int_set_value(mp_int z, mp_small value);
 138
 139 /** Sets `z` to the value of the specified unsigned `value`. */
 140 mp_result mp_int_set_uvalue(mp_int z, mp_usmall uvalue);
 141
 142 /** Releases the storage used by `z`. */
 143 void 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()`. */
 147 void 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. */
 151 mp_result mp_int_copy(mp_int a, mp_int c);
 152
 153 /** Swaps the values and storage between `a` and `c`. */
 154 void mp_int_swap(mp_int a, mp_int c);
 155
 156 /** Sets `z` to zero. The allocated storage of `z` is not changed. */
 157 void mp_int_zero(mp_int z);
 158
 159 /** Sets `c` to the absolute value of `a`. */
 160 mp_result mp_int_abs(mp_int a, mp_int c);
 161
 162 /** Sets `c` to the additive inverse (negation) of `a`. */
 163 mp_result mp_int_neg(mp_int a, mp_int c);
 164
 165 /** Sets `c` to the sum of `a` and `b`. */
 166 mp_result mp_int_add(mp_int a, mp_int b, mp_int c);
 167
 168 /** Sets `c` to the sum of `a` and `value`. */
 169 mp_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`. */
 172 mp_result mp_int_sub(mp_int a, mp_int b, mp_int c);
 173
 174 /** Sets `c` to the difference of `a` less `value`. */
 175 mp_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`. */
 178 mp_result mp_int_mul(mp_int a, mp_int b, mp_int c);
 179
 180 /** Sets `c` to the product of `a` and `value`. */
 181 mp_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`. */
 184 mp_result mp_int_mul_pow2(mp_int a, mp_small p2, mp_int c);
 185
 186 /** Sets `c` to the square of `a`. */
 187 mp_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. */
 195 mp_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. */
 200 mp_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. */
 206 mp_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`. */
 210 mp_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`. */
 214 mp_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`. */
 218 mp_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`. */
 222 mp_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`. */
 226 static inline
 227 mp_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`. */
 232 int 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. */
 236 int mp_int_compare_unsigned(mp_int a, mp_int b);
 237
 238 /** Returns the comparator of `z` and zero. */
 239 int mp_int_compare_zero(mp_int z);
 240
 241 /** Returns the comparator of `z` and the signed value `v`. */
 242 int mp_int_compare_value(mp_int z, mp_small v);
 243
 244 /** Returns the comparator of `z` and the unsigned value `uv`. */
 245 int mp_int_compare_uvalue(mp_int z, mp_usmall uv);
 246
 247 /** Reports whether `a` is divisible by `v`. */
 248 bool 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. */
 252 int 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`. */
 256 mp_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`. */
 260 mp_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`. */
 264 mp_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`. */
 271 mp_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. */
 275 mp_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`. */
 282 mp_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. */
 288 mp_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. */
 295 mp_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. */
 301 mp_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. */
 306 mp_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()`. */
 310 static inline
 311 mp_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`. */
 315 mp_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`. */
 319 mp_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`. */
 327 mp_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`. */
 332 mp_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`. */
 347 mp_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`. */
 365 mp_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`. */
 368 mp_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. */
 383 mp_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. */
 388 mp_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. */
 391 mp_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). */
 402 mp_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. */
 407 mp_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. */
 411 mp_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. */
 416 const char *mp_error_string(mp_result res);
 417
 418 #ifdef __cplusplus
 419 }
 420 #endif
 421 #endif /* end IMATH_H_ */