/* Name: imath.c Purpose: Arbitrary precision integer arithmetic routines. Author: M. J. Fromberger Copyright (C) 2002-2007 Michael J. Fromberger, All Rights Reserved. Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ #include "imath.h" #include <assert.h> #include <ctype.h> #include <stdlib.h> #include <string.h> const mp_result MP_OK = …; /* no error, all is well */ const mp_result MP_FALSE = …; /* boolean false */ const mp_result MP_TRUE = …; /* boolean true */ const mp_result MP_MEMORY = …; /* out of memory */ const mp_result MP_RANGE = …; /* argument out of range */ const mp_result MP_UNDEF = …; /* result undefined */ const mp_result MP_TRUNC = …; /* output truncated */ const mp_result MP_BADARG = …; /* invalid null argument */ const mp_result MP_MINERR = …; const mp_sign MP_NEG = …; /* value is strictly negative */ const mp_sign MP_ZPOS = …; /* value is non-negative */ static const char *s_unknown_err = …; static const char *s_error_msg[] = …; /* The ith entry of this table gives the value of log_i(2). An integer value n requires ceil(log_i(n)) digits to be represented in base i. Since it is easy to compute lg(n), by counting bits, we can compute log_i(n) = lg(n) * log_i(2). The use of this table eliminates a dependency upon linkage against the standard math libraries. If MP_MAX_RADIX is increased, this table should be expanded too. */ static const double s_log2[] = …; /* Return the number of digits needed to represent a static value */ #define MP_VALUE_DIGITS(V) … /* Round precision P to nearest word boundary */ static inline mp_size s_round_prec(mp_size P) { … } /* Set array P of S digits to zero */ static inline void ZERO(mp_digit *P, mp_size S) { … } /* Copy S digits from array P to array Q */ static inline void COPY(mp_digit *P, mp_digit *Q, mp_size S) { … } /* Reverse N elements of unsigned char in A. */ static inline void REV(unsigned char *A, int N) { … } /* Strip leading zeroes from z_ in-place. */ static inline void CLAMP(mp_int z_) { … } /* Select min/max. */ static inline int MIN(int A, int B) { … } static inline mp_size MAX(mp_size A, mp_size B) { … } /* Exchange lvalues A and B of type T, e.g. SWAP(int, x, y) where x and y are variables of type int. */ #define SWAP(T, A, B) … /* Declare a block of N temporary mpz_t values. These values are initialized to zero. You must add CLEANUP_TEMP() at the end of the function. Use TEMP(i) to access a pointer to the ith value. */ #define DECLARE_TEMP(N) … /* Clear all allocated temp values. */ #define CLEANUP_TEMP() … /* A pointer to the kth temp value. */ #define TEMP(K) … /* Evaluate E, an expression of type mp_result expected to return MP_OK. If the value is not MP_OK, the error is cached and control resumes at the cleanup handler, which returns it. */ #define REQUIRE(E) … /* Compare value to zero. */ static inline int CMPZ(mp_int Z) { … } static inline mp_word UPPER_HALF(mp_word W) { … } static inline mp_digit LOWER_HALF(mp_word W) { … } /* Report whether the highest-order bit of W is 1. */ static inline bool HIGH_BIT_SET(mp_word W) { … } /* Report whether adding W + V will carry out. */ static inline bool ADD_WILL_OVERFLOW(mp_word W, mp_word V) { … } /* Default number of digits allocated to a new mp_int */ static mp_size default_precision = …; void mp_int_default_precision(mp_size size) { … } /* Minimum number of digits to invoke recursive multiply */ static mp_size multiply_threshold = …; void mp_int_multiply_threshold(mp_size thresh) { … } /* Allocate a buffer of (at least) num digits, or return NULL if that couldn't be done. */ static mp_digit *s_alloc(mp_size num); /* Release a buffer of digits allocated by s_alloc(). */ static void s_free(void *ptr); /* Insure that z has at least min digits allocated, resizing if necessary. Returns true if successful, false if out of memory. */ static bool s_pad(mp_int z, mp_size min); /* Ensure Z has at least N digits allocated. */ static inline mp_result GROW(mp_int Z, mp_size N) { … } /* Fill in a "fake" mp_int on the stack with a given value */ static void s_fake(mp_int z, mp_small value, mp_digit vbuf[]); static void s_ufake(mp_int z, mp_usmall value, mp_digit vbuf[]); /* Compare two runs of digits of given length, returns <0, 0, >0 */ static int s_cdig(mp_digit *da, mp_digit *db, mp_size len); /* Pack the unsigned digits of v into array t */ static int s_uvpack(mp_usmall v, mp_digit t[]); /* Compare magnitudes of a and b, returns <0, 0, >0 */ static int s_ucmp(mp_int a, mp_int b); /* Compare magnitudes of a and v, returns <0, 0, >0 */ static int s_vcmp(mp_int a, mp_small v); static int s_uvcmp(mp_int a, mp_usmall uv); /* Unsigned magnitude addition; assumes dc is big enough. Carry out is returned (no memory allocated). */ static mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, mp_size size_b); /* Unsigned magnitude subtraction. Assumes dc is big enough. */ static void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, mp_size size_b); /* Unsigned recursive multiplication. Assumes dc is big enough. */ static int s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, mp_size size_b); /* Unsigned magnitude multiplication. Assumes dc is big enough. */ static void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, mp_size size_b); /* Unsigned recursive squaring. Assumes dc is big enough. */ static int s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a); /* Unsigned magnitude squaring. Assumes dc is big enough. */ static void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a); /* Single digit addition. Assumes a is big enough. */ static void s_dadd(mp_int a, mp_digit b); /* Single digit multiplication. Assumes a is big enough. */ static void s_dmul(mp_int a, mp_digit b); /* Single digit multiplication on buffers; assumes dc is big enough. */ static void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a); /* Single digit division. Replaces a with the quotient, returns the remainder. */ static mp_digit s_ddiv(mp_int a, mp_digit b); /* Quick division by a power of 2, replaces z (no allocation) */ static void s_qdiv(mp_int z, mp_size p2); /* Quick remainder by a power of 2, replaces z (no allocation) */ static void s_qmod(mp_int z, mp_size p2); /* Quick multiplication by a power of 2, replaces z. Allocates if necessary; returns false in case this fails. */ static int s_qmul(mp_int z, mp_size p2); /* Quick subtraction from a power of 2, replaces z. Allocates if necessary; returns false in case this fails. */ static int s_qsub(mp_int z, mp_size p2); /* Return maximum k such that 2^k divides z. */ static int s_dp2k(mp_int z); /* Return k >= 0 such that z = 2^k, or -1 if there is no such k. */ static int s_isp2(mp_int z); /* Set z to 2^k. May allocate; returns false in case this fails. */ static int s_2expt(mp_int z, mp_small k); /* Normalize a and b for division, returns normalization constant */ static int s_norm(mp_int a, mp_int b); /* Compute constant mu for Barrett reduction, given modulus m, result replaces z, m is untouched. */ static mp_result s_brmu(mp_int z, mp_int m); /* Reduce a modulo m, using Barrett's algorithm. */ static int s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2); /* Modular exponentiation, using Barrett reduction */ static mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c); /* Unsigned magnitude division. Assumes |a| > |b|. Allocates temporaries; overwrites a with quotient, b with remainder. */ static mp_result s_udiv_knuth(mp_int a, mp_int b); /* Compute the number of digits in radix r required to represent the given value. Does not account for sign flags, terminators, etc. */ static int s_outlen(mp_int z, mp_size r); /* Guess how many digits of precision will be needed to represent a radix r value of the specified number of digits. Returns a value guaranteed to be no smaller than the actual number required. */ static mp_size s_inlen(int len, mp_size r); /* Convert a character to a digit value in radix r, or -1 if out of range */ static int s_ch2val(char c, int r); /* Convert a digit value to a character */ static char s_val2ch(int v, int caps); /* Take 2's complement of a buffer in place */ static void s_2comp(unsigned char *buf, int len); /* Convert a value to binary, ignoring sign. On input, *limpos is the bound on how many bytes should be written to buf; on output, *limpos is set to the number of bytes actually written. */ static mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad); /* Multiply X by Y into Z, ignoring signs. Requires that Z have enough storage preallocated to hold the result. */ static inline void UMUL(mp_int X, mp_int Y, mp_int Z) { … } /* Square X into Z. Requires that Z have enough storage to hold the result. */ static inline void USQR(mp_int X, mp_int Z) { … } mp_result mp_int_init(mp_int z) { … } mp_int mp_int_alloc(void) { … } mp_result mp_int_init_size(mp_int z, mp_size prec) { … } mp_result mp_int_init_copy(mp_int z, mp_int old) { … } mp_result mp_int_init_value(mp_int z, mp_small value) { … } mp_result mp_int_init_uvalue(mp_int z, mp_usmall uvalue) { … } mp_result mp_int_set_value(mp_int z, mp_small value) { … } mp_result mp_int_set_uvalue(mp_int z, mp_usmall uvalue) { … } void mp_int_clear(mp_int z) { … } void mp_int_free(mp_int z) { … } mp_result mp_int_copy(mp_int a, mp_int c) { … } void mp_int_swap(mp_int a, mp_int c) { … } void mp_int_zero(mp_int z) { … } mp_result mp_int_abs(mp_int a, mp_int c) { … } mp_result mp_int_neg(mp_int a, mp_int c) { … } mp_result mp_int_add(mp_int a, mp_int b, mp_int c) { … } mp_result mp_int_add_value(mp_int a, mp_small value, mp_int c) { … } mp_result mp_int_sub(mp_int a, mp_int b, mp_int c) { … } mp_result mp_int_sub_value(mp_int a, mp_small value, mp_int c) { … } mp_result mp_int_mul(mp_int a, mp_int b, mp_int c) { … } mp_result mp_int_mul_value(mp_int a, mp_small value, mp_int c) { … } mp_result mp_int_mul_pow2(mp_int a, mp_small p2, mp_int c) { … } mp_result mp_int_sqr(mp_int a, mp_int c) { … } mp_result mp_int_div(mp_int a, mp_int b, mp_int q, mp_int r) { … } mp_result mp_int_mod(mp_int a, mp_int m, mp_int c) { … } mp_result mp_int_div_value(mp_int a, mp_small value, mp_int q, mp_small *r) { … } mp_result mp_int_div_pow2(mp_int a, mp_small p2, mp_int q, mp_int r) { … } mp_result mp_int_expt(mp_int a, mp_small b, mp_int c) { … } mp_result mp_int_expt_value(mp_small a, mp_small b, mp_int c) { … } mp_result mp_int_expt_full(mp_int a, mp_int b, mp_int c) { … } int mp_int_compare(mp_int a, mp_int b) { … } int mp_int_compare_unsigned(mp_int a, mp_int b) { … } int mp_int_compare_zero(mp_int z) { … } int mp_int_compare_value(mp_int z, mp_small value) { … } int mp_int_compare_uvalue(mp_int z, mp_usmall uv) { … } mp_result mp_int_exptmod(mp_int a, mp_int b, mp_int m, mp_int c) { … } mp_result mp_int_exptmod_evalue(mp_int a, mp_small value, mp_int m, mp_int c) { … } mp_result mp_int_exptmod_bvalue(mp_small value, mp_int b, mp_int m, mp_int c) { … } mp_result mp_int_exptmod_known(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c) { … } mp_result mp_int_redux_const(mp_int m, mp_int c) { … } mp_result mp_int_invmod(mp_int a, mp_int m, mp_int c) { … } /* Binary GCD algorithm due to Josef Stein, 1961 */ mp_result mp_int_gcd(mp_int a, mp_int b, mp_int c) { … } /* This is the binary GCD algorithm again, but this time we keep track of the elementary matrix operations as we go, so we can get values x and y satisfying c = ax + by. */ mp_result mp_int_egcd(mp_int a, mp_int b, mp_int c, mp_int x, mp_int y) { … } mp_result mp_int_lcm(mp_int a, mp_int b, mp_int c) { … } bool mp_int_divisible_value(mp_int a, mp_small v) { … } int mp_int_is_pow2(mp_int z) { … } /* Implementation of Newton's root finding method, based loosely on a patch contributed by Hal Finkel <[email protected]> modified by M. J. Fromberger. */ mp_result mp_int_root(mp_int a, mp_small b, mp_int c) { … } mp_result mp_int_to_int(mp_int z, mp_small *out) { … } mp_result mp_int_to_uint(mp_int z, mp_usmall *out) { … } mp_result mp_int_to_string(mp_int z, mp_size radix, char *str, int limit) { … } mp_result mp_int_string_len(mp_int z, mp_size radix) { … } /* Read zero-terminated string into z */ mp_result mp_int_read_string(mp_int z, mp_size radix, const char *str) { … } mp_result mp_int_read_cstring(mp_int z, mp_size radix, const char *str, char **end) { … } mp_result mp_int_count_bits(mp_int z) { … } mp_result mp_int_to_binary(mp_int z, unsigned char *buf, int limit) { … } mp_result mp_int_read_binary(mp_int z, unsigned char *buf, int len) { … } mp_result mp_int_binary_len(mp_int z) { … } mp_result mp_int_to_unsigned(mp_int z, unsigned char *buf, int limit) { … } mp_result mp_int_read_unsigned(mp_int z, unsigned char *buf, int len) { … } mp_result mp_int_unsigned_len(mp_int z) { … } const char *mp_error_string(mp_result res) { … } /*------------------------------------------------------------------------*/ /* Private functions for internal use. These make assumptions. */ #if DEBUG static const mp_digit fill = (mp_digit)0xdeadbeefabad1dea; #endif static mp_digit *s_alloc(mp_size num) { … } static mp_digit *s_realloc(mp_digit *old, mp_size osize, mp_size nsize) { … } static void s_free(void *ptr) { … } static bool s_pad(mp_int z, mp_size min) { … } /* Note: This will not work correctly when value == MP_SMALL_MIN */ static void s_fake(mp_int z, mp_small value, mp_digit vbuf[]) { … } static void s_ufake(mp_int z, mp_usmall value, mp_digit vbuf[]) { … } static int s_cdig(mp_digit *da, mp_digit *db, mp_size len) { … } static int s_uvpack(mp_usmall uv, mp_digit t[]) { … } static int s_ucmp(mp_int a, mp_int b) { … } static int s_vcmp(mp_int a, mp_small v) { … } static int s_uvcmp(mp_int a, mp_usmall uv) { … } static mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, mp_size size_b) { … } static void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, mp_size size_b) { … } static int s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, mp_size size_b) { … } static void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, mp_size size_b) { … } static int s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a) { … } static void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a) { … } static void s_dadd(mp_int a, mp_digit b) { … } static void s_dmul(mp_int a, mp_digit b) { … } static void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a) { … } static mp_digit s_ddiv(mp_int a, mp_digit b) { … } static void s_qdiv(mp_int z, mp_size p2) { … } static void s_qmod(mp_int z, mp_size p2) { … } static int s_qmul(mp_int z, mp_size p2) { … } /* Compute z = 2^p2 - |z|; requires that 2^p2 >= |z| The sign of the result is always zero/positive. */ static int s_qsub(mp_int z, mp_size p2) { … } static int s_dp2k(mp_int z) { … } static int s_isp2(mp_int z) { … } static int s_2expt(mp_int z, mp_small k) { … } static int s_norm(mp_int a, mp_int b) { … } static mp_result s_brmu(mp_int z, mp_int m) { … } static int s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2) { … } /* Perform modular exponentiation using Barrett's method, where mu is the reduction constant for m. Assumes a < m, b > 0. */ static mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c) { … } /* Division of nonnegative integers This function implements division algorithm for unsigned multi-precision integers. The algorithm is based on Algorithm D from Knuth's "The Art of Computer Programming", 3rd ed. 1998, pg 272-273. We diverge from Knuth's algorithm in that we do not perform the subtraction from the remainder until we have determined that we have the correct quotient digit. This makes our algorithm less efficient that Knuth because we might have to perform multiple multiplication and comparison steps before the subtraction. The advantage is that it is easy to implement and ensure correctness without worrying about underflow from the subtraction. inputs: u a n+m digit integer in base b (b is 2^MP_DIGIT_BIT) v a n digit integer in base b (b is 2^MP_DIGIT_BIT) n >= 1 m >= 0 outputs: u / v stored in u u % v stored in v */ static mp_result s_udiv_knuth(mp_int u, mp_int v) { … } static int s_outlen(mp_int z, mp_size r) { … } static mp_size s_inlen(int len, mp_size r) { … } static int s_ch2val(char c, int r) { … } static char s_val2ch(int v, int caps) { … } static void s_2comp(unsigned char *buf, int len) { … } static mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad) { … } /* Here there be dragons */
Codebase Browser
llvm