diff options
Diffstat (limited to 'contrib/pgcrypto/imath.c')
| -rw-r--r-- | contrib/pgcrypto/imath.c | 3588 |
1 files changed, 0 insertions, 3588 deletions
diff --git a/contrib/pgcrypto/imath.c b/contrib/pgcrypto/imath.c deleted file mode 100644 index 0bfa080fa5..0000000000 --- a/contrib/pgcrypto/imath.c +++ /dev/null @@ -1,3588 +0,0 @@ -/*------------------------------------------------------------------------- - * - * imath.c - * - * Last synchronized from https://github.com/creachadair/imath/tree/v1.29, - * using the following procedure: - * - * 1. Download imath.c and imath.h of the last synchronized version. Remove - * "#ifdef __cplusplus" blocks, which upset pgindent. Run pgindent on the - * two files. Filter the two files through "unexpand -t4 --first-only". - * Diff the result against the PostgreSQL versions. As of the last - * synchronization, changes were as follows: - * - * - replace malloc(), realloc() and free() with px_ versions - * - redirect assert() to Assert() - * - #undef MIN, #undef MAX before defining them - * - remove includes covered by c.h - * - rename DEBUG to IMATH_DEBUG - * - replace stdint.h usage with c.h equivalents - * - suppress MSVC warning 4146 - * - add required PG_USED_FOR_ASSERTS_ONLY - * - * 2. Download a newer imath.c and imath.h. Transform them like in step 1. - * Apply to these files the diff you saved in step 1. Look for new lines - * requiring the same kind of change, such as new malloc() calls. - * - * 3. Configure PostgreSQL using --without-openssl. Run "make -C - * contrib/pgcrypto check". - * - * 4. Update this header comment. - * - * Portions Copyright (c) 1996-2021, PostgreSQL Global Development Group - * - * IDENTIFICATION - * contrib/pgcrypto/imath.c - * - * Upstream copyright terms follow. - *------------------------------------------------------------------------- - */ - -/* - 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 "postgres.h" - -#include "imath.h" -#include "px.h" - -#undef assert -#define assert(TEST) Assert(TEST) - -const mp_result MP_OK = 0; /* no error, all is well */ -const mp_result MP_FALSE = 0; /* boolean false */ -const mp_result MP_TRUE = -1; /* boolean true */ -const mp_result MP_MEMORY = -2; /* out of memory */ -const mp_result MP_RANGE = -3; /* argument out of range */ -const mp_result MP_UNDEF = -4; /* result undefined */ -const mp_result MP_TRUNC = -5; /* output truncated */ -const mp_result MP_BADARG = -6; /* invalid null argument */ -const mp_result MP_MINERR = -6; - -const mp_sign MP_NEG = 1; /* value is strictly negative */ -const mp_sign MP_ZPOS = 0; /* value is non-negative */ - -static const char *s_unknown_err = "unknown result code"; -static const char *s_error_msg[] = {"error code 0", "boolean true", - "out of memory", "argument out of range", - "result undefined", "output truncated", -"invalid argument", NULL}; - -/* 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[] = { - 0.000000000, 0.000000000, 1.000000000, 0.630929754, /* (D)(D) 2 3 */ - 0.500000000, 0.430676558, 0.386852807, 0.356207187, /* 4 5 6 7 */ - 0.333333333, 0.315464877, 0.301029996, 0.289064826, /* 8 9 10 11 */ - 0.278942946, 0.270238154, 0.262649535, 0.255958025, /* 12 13 14 15 */ - 0.250000000, 0.244650542, 0.239812467, 0.235408913, /* 16 17 18 19 */ - 0.231378213, 0.227670249, 0.224243824, 0.221064729, /* 20 21 22 23 */ - 0.218104292, 0.215338279, 0.212746054, 0.210309918, /* 24 25 26 27 */ - 0.208014598, 0.205846832, 0.203795047, 0.201849087, /* 28 29 30 31 */ - 0.200000000, 0.198239863, 0.196561632, 0.194959022, /* 32 33 34 35 */ - 0.193426404, /* 36 */ -}; - -/* Return the number of digits needed to represent a static value */ -#define MP_VALUE_DIGITS(V) \ - ((sizeof(V) + (sizeof(mp_digit) - 1)) / sizeof(mp_digit)) - -/* Round precision P to nearest word boundary */ -static inline mp_size -s_round_prec(mp_size P) -{ - return 2 * ((P + 1) / 2); -} - -/* Set array P of S digits to zero */ -static inline void -ZERO(mp_digit *P, mp_size S) -{ - mp_size i__ = S * sizeof(mp_digit); - mp_digit *p__ = P; - - memset(p__, 0, i__); -} - -/* Copy S digits from array P to array Q */ -static inline void -COPY(mp_digit *P, mp_digit *Q, mp_size S) -{ - mp_size i__ = S * sizeof(mp_digit); - mp_digit *p__ = P; - mp_digit *q__ = Q; - - memcpy(q__, p__, i__); -} - -/* Reverse N elements of unsigned char in A. */ -static inline void -REV(unsigned char *A, int N) -{ - unsigned char *u_ = A; - unsigned char *v_ = u_ + N - 1; - - while (u_ < v_) - { - unsigned char xch = *u_; - - *u_++ = *v_; - *v_-- = xch; - } -} - -/* Strip leading zeroes from z_ in-place. */ -static inline void -CLAMP(mp_int z_) -{ - mp_size uz_ = MP_USED(z_); - mp_digit *dz_ = MP_DIGITS(z_) + uz_ - 1; - - while (uz_ > 1 && (*dz_-- == 0)) - --uz_; - z_->used = uz_; -} - -/* Select min/max. */ -#undef MIN -#undef MAX -static inline int -MIN(int A, int B) -{ - return (B < A ? B : A); -} -static inline mp_size -MAX(mp_size A, mp_size B) -{ - return (B > A ? B : A); -} - -/* 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) \ - do { \ - T t_ = (A); \ - A = (B); \ - B = t_; \ - } while (0) - -/* 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) \ - struct { \ - mpz_t value[(N)]; \ - int len; \ - mp_result err; \ - } temp_ = { \ - .len = (N), \ - .err = MP_OK, \ - }; \ - do { \ - for (int i = 0; i < temp_.len; i++) { \ - mp_int_init(TEMP(i)); \ - } \ - } while (0) - -/* Clear all allocated temp values. */ -#define CLEANUP_TEMP() \ - CLEANUP: \ - do { \ - for (int i = 0; i < temp_.len; i++) { \ - mp_int_clear(TEMP(i)); \ - } \ - if (temp_.err != MP_OK) { \ - return temp_.err; \ - } \ - } while (0) - -/* A pointer to the kth temp value. */ -#define TEMP(K) (temp_.value + (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) \ - do { \ - temp_.err = (E); \ - if (temp_.err != MP_OK) goto CLEANUP; \ - } while (0) - -/* Compare value to zero. */ -static inline int -CMPZ(mp_int Z) -{ - if (Z->used == 1 && Z->digits[0] == 0) - return 0; - return (Z->sign == MP_NEG) ? -1 : 1; -} - -static inline mp_word -UPPER_HALF(mp_word W) -{ - return (W >> MP_DIGIT_BIT); -} -static inline mp_digit -LOWER_HALF(mp_word W) -{ - return (mp_digit) (W); -} - -/* Report whether the highest-order bit of W is 1. */ -static inline bool -HIGH_BIT_SET(mp_word W) -{ - return (W >> (MP_WORD_BIT - 1)) != 0; -} - -/* Report whether adding W + V will carry out. */ -static inline bool -ADD_WILL_OVERFLOW(mp_word W, mp_word V) -{ - return ((MP_WORD_MAX - V) < W); -} - -/* Default number of digits allocated to a new mp_int */ -static mp_size default_precision = 8; - -void -mp_int_default_precision(mp_size size) -{ - assert(size > 0); - default_precision = size; -} - -/* Minimum number of digits to invoke recursive multiply */ -static mp_size multiply_threshold = 32; - -void -mp_int_multiply_threshold(mp_size thresh) -{ - assert(thresh >= sizeof(mp_word)); - multiply_threshold = 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) -{ - return s_pad(Z, N) ? MP_OK : MP_MEMORY; -} - -/* 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) -{ - mp_size ua_ = MP_USED(X); - mp_size ub_ = MP_USED(Y); - mp_size o_ = ua_ + ub_; - - ZERO(MP_DIGITS(Z), o_); - (void) s_kmul(MP_DIGITS(X), MP_DIGITS(Y), MP_DIGITS(Z), ua_, ub_); - Z->used = o_; - CLAMP(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_size ua_ = MP_USED(X); - mp_size o_ = ua_ + ua_; - - ZERO(MP_DIGITS(Z), o_); - (void) s_ksqr(MP_DIGITS(X), MP_DIGITS(Z), ua_); - Z->used = o_; - CLAMP(Z); -} - -mp_result -mp_int_init(mp_int z) -{ - if (z == NULL) - return MP_BADARG; - - z->single = 0; - z->digits = &(z->single); - z->alloc = 1; - z->used = 1; - z->sign = MP_ZPOS; - - return MP_OK; -} - -mp_int -mp_int_alloc(void) -{ - mp_int out = palloc(sizeof(mpz_t)); - - if (out != NULL) - mp_int_init(out); - - return out; -} - -mp_result -mp_int_init_size(mp_int z, mp_size prec) -{ - assert(z != NULL); - - if (prec == 0) - { - prec = default_precision; - } - else if (prec == 1) - { - return mp_int_init(z); - } - else - { - prec = s_round_prec(prec); - } - - z->digits = s_alloc(prec); - if (MP_DIGITS(z) == NULL) - return MP_MEMORY; - - z->digits[0] = 0; - z->used = 1; - z->alloc = prec; - z->sign = MP_ZPOS; - - return MP_OK; -} - -mp_result -mp_int_init_copy(mp_int z, mp_int old) -{ - assert(z != NULL && old != NULL); - - mp_size uold = MP_USED(old); - - if (uold == 1) - { - mp_int_init(z); - } - else - { - mp_size target = MAX(uold, default_precision); - mp_result res = mp_int_init_size(z, target); - - if (res != MP_OK) - return res; - } - - z->used = uold; - z->sign = old->sign; - COPY(MP_DIGITS(old), MP_DIGITS(z), uold); - - return MP_OK; -} - -mp_result -mp_int_init_value(mp_int z, mp_small value) -{ - mpz_t vtmp; - mp_digit vbuf[MP_VALUE_DIGITS(value)]; - - s_fake(&vtmp, value, vbuf); - return mp_int_init_copy(z, &vtmp); -} - -mp_result -mp_int_init_uvalue(mp_int z, mp_usmall uvalue) -{ - mpz_t vtmp; - mp_digit vbuf[MP_VALUE_DIGITS(uvalue)]; - - s_ufake(&vtmp, uvalue, vbuf); - return mp_int_init_copy(z, &vtmp); -} - -mp_result -mp_int_set_value(mp_int z, mp_small value) -{ - mpz_t vtmp; - mp_digit vbuf[MP_VALUE_DIGITS(value)]; - - s_fake(&vtmp, value, vbuf); - return mp_int_copy(&vtmp, z); -} - -mp_result -mp_int_set_uvalue(mp_int z, mp_usmall uvalue) -{ - mpz_t vtmp; - mp_digit vbuf[MP_VALUE_DIGITS(uvalue)]; - - s_ufake(&vtmp, uvalue, vbuf); - return mp_int_copy(&vtmp, z); -} - -void -mp_int_clear(mp_int z) -{ - if (z == NULL) - return; - - if (MP_DIGITS(z) != NULL) - { - if (MP_DIGITS(z) != &(z->single)) - s_free(MP_DIGITS(z)); - - z->digits = NULL; - } -} - -void -mp_int_free(mp_int z) -{ - assert(z != NULL); - - mp_int_clear(z); - pfree(z); /* note: NOT s_free() */ -} - -mp_result -mp_int_copy(mp_int a, mp_int c) -{ - assert(a != NULL && c != NULL); - - if (a != c) - { - mp_size ua = MP_USED(a); - mp_digit *da, - *dc; - - if (!s_pad(c, ua)) - return MP_MEMORY; - - da = MP_DIGITS(a); - dc = MP_DIGITS(c); - COPY(da, dc, ua); - - c->used = ua; - c->sign = a->sign; - } - - return MP_OK; -} - -void -mp_int_swap(mp_int a, mp_int c) -{ - if (a != c) - { - mpz_t tmp = *a; - - *a = *c; - *c = tmp; - - if (MP_DIGITS(a) == &(c->single)) - a->digits = &(a->single); - if (MP_DIGITS(c) == &(a->single)) - c->digits = &(c->single); - } -} - -void -mp_int_zero(mp_int z) -{ - assert(z != NULL); - - z->digits[0] = 0; - z->used = 1; - z->sign = MP_ZPOS; -} - -mp_result -mp_int_abs(mp_int a, mp_int c) -{ - assert(a != NULL && c != NULL); - - mp_result res; - - if ((res = mp_int_copy(a, c)) != MP_OK) - return res; - - c->sign = MP_ZPOS; - return MP_OK; -} - -mp_result -mp_int_neg(mp_int a, mp_int c) -{ - assert(a != NULL && c != NULL); - - mp_result res; - - if ((res = mp_int_copy(a, c)) != MP_OK) - return res; - - if (CMPZ(c) != 0) - c->sign = 1 - MP_SIGN(a); - - return MP_OK; -} - -mp_result -mp_int_add(mp_int a, mp_int b, mp_int c) -{ - assert(a != NULL && b != NULL && c != NULL); - - mp_size ua = MP_USED(a); - mp_size ub = MP_USED(b); - mp_size max = MAX(ua, ub); - - if (MP_SIGN(a) == MP_SIGN(b)) - { - /* Same sign -- add magnitudes, preserve sign of addends */ - if (!s_pad(c, max)) - return MP_MEMORY; - - mp_digit carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub); - mp_size uc = max; - - if (carry) - { - if (!s_pad(c, max + 1)) - return MP_MEMORY; - - c->digits[max] = carry; - ++uc; - } - - c->used = uc; - c->sign = a->sign; - - } - else - { - /* Different signs -- subtract magnitudes, preserve sign of greater */ - int cmp = s_ucmp(a, b); /* magnitude comparision, sign ignored */ - - /* - * Set x to max(a, b), y to min(a, b) to simplify later code. A - * special case yields zero for equal magnitudes. - */ - mp_int x, - y; - - if (cmp == 0) - { - mp_int_zero(c); - return MP_OK; - } - else if (cmp < 0) - { - x = b; - y = a; - } - else - { - x = a; - y = b; - } - - if (!s_pad(c, MP_USED(x))) - return MP_MEMORY; - - /* Subtract smaller from larger */ - s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y)); - c->used = x->used; - CLAMP(c); - - /* Give result the sign of the larger */ - c->sign = x->sign; - } - - return MP_OK; -} - -mp_result -mp_int_add_value(mp_int a, mp_small value, mp_int c) -{ - mpz_t vtmp; - mp_digit vbuf[MP_VALUE_DIGITS(value)]; - - s_fake(&vtmp, value, vbuf); - - return mp_int_add(a, &vtmp, c); -} - -mp_result -mp_int_sub(mp_int a, mp_int b, mp_int c) -{ - assert(a != NULL && b != NULL && c != NULL); - - mp_size ua = MP_USED(a); - mp_size ub = MP_USED(b); - mp_size max = MAX(ua, ub); - - if (MP_SIGN(a) != MP_SIGN(b)) - { - /* Different signs -- add magnitudes and keep sign of a */ - if (!s_pad(c, max)) - return MP_MEMORY; - - mp_digit carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub); - mp_size uc = max; - - if (carry) - { - if (!s_pad(c, max + 1)) - return MP_MEMORY; - - c->digits[max] = carry; - ++uc; - } - - c->used = uc; - c->sign = a->sign; - - } - else - { - /* Same signs -- subtract magnitudes */ - if (!s_pad(c, max)) - return MP_MEMORY; - mp_int x, - y; - mp_sign osign; - - int cmp = s_ucmp(a, b); - - if (cmp >= 0) - { - x = a; - y = b; - osign = MP_ZPOS; - } - else - { - x = b; - y = a; - osign = MP_NEG; - } - - if (MP_SIGN(a) == MP_NEG && cmp != 0) - osign = 1 - osign; - - s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y)); - c->used = x->used; - CLAMP(c); - - c->sign = osign; - } - - return MP_OK; -} - -mp_result -mp_int_sub_value(mp_int a, mp_small value, mp_int c) -{ - mpz_t vtmp; - mp_digit vbuf[MP_VALUE_DIGITS(value)]; - - s_fake(&vtmp, value, vbuf); - - return mp_int_sub(a, &vtmp, c); -} - -mp_result -mp_int_mul(mp_int a, mp_int b, mp_int c) -{ - assert(a != NULL && b != NULL && c != NULL); - - /* If either input is zero, we can shortcut multiplication */ - if (mp_int_compare_zero(a) == 0 || mp_int_compare_zero(b) == 0) - { - mp_int_zero(c); - return MP_OK; - } - - /* Output is positive if inputs have same sign, otherwise negative */ - mp_sign osign = (MP_SIGN(a) == MP_SIGN(b)) ? MP_ZPOS : MP_NEG; - - /* - * If the output is not identical to any of the inputs, we'll write the - * results directly; otherwise, allocate a temporary space. - */ - mp_size ua = MP_USED(a); - mp_size ub = MP_USED(b); - mp_size osize = MAX(ua, ub); - - osize = 4 * ((osize + 1) / 2); - - mp_digit *out; - mp_size p = 0; - - if (c == a || c == b) - { - p = MAX(s_round_prec(osize), default_precision); - - if ((out = s_alloc(p)) == NULL) - return MP_MEMORY; - } - else - { - if (!s_pad(c, osize)) - return MP_MEMORY; - - out = MP_DIGITS(c); - } - ZERO(out, osize); - - if (!s_kmul(MP_DIGITS(a), MP_DIGITS(b), out, ua, ub)) - return MP_MEMORY; - - /* - * If we allocated a new buffer, get rid of whatever memory c was already - * using, and fix up its fields to reflect that. - */ - if (out != MP_DIGITS(c)) - { - if ((void *) MP_DIGITS(c) != (void *) c) - s_free(MP_DIGITS(c)); - c->digits = out; - c->alloc = p; - } - - c->used = osize; /* might not be true, but we'll fix it ... */ - CLAMP(c); /* ... right here */ - c->sign = osign; - - return MP_OK; -} - -mp_result -mp_int_mul_value(mp_int a, mp_small value, mp_int c) -{ - mpz_t vtmp; - mp_digit vbuf[MP_VALUE_DIGITS(value)]; - - s_fake(&vtmp, value, vbuf); - - return mp_int_mul(a, &vtmp, c); -} - -mp_result -mp_int_mul_pow2(mp_int a, mp_small p2, mp_int c) -{ - assert(a != NULL && c != NULL && p2 >= 0); - - mp_result res = mp_int_copy(a, c); - - if (res != MP_OK) - return res; - - if (s_qmul(c, (mp_size) p2)) - { - return MP_OK; - } - else - { - return MP_MEMORY; - } -} - -mp_result -mp_int_sqr(mp_int a, mp_int c) -{ - assert(a != NULL && c != NULL); - - /* Get a temporary buffer big enough to hold the result */ - mp_size osize = (mp_size) 4 * ((MP_USED(a) + 1) / 2); - mp_size p = 0; - mp_digit *out; - - if (a == c) - { - p = s_round_prec(osize); - p = MAX(p, default_precision); - - if ((out = s_alloc(p)) == NULL) - return MP_MEMORY; - } - else - { - if (!s_pad(c, osize)) - return MP_MEMORY; - - out = MP_DIGITS(c); - } - ZERO(out, osize); - - s_ksqr(MP_DIGITS(a), out, MP_USED(a)); - - /* - * Get rid of whatever memory c was already using, and fix up its fields - * to reflect the new digit array it's using - */ - if (out != MP_DIGITS(c)) - { - if ((void *) MP_DIGITS(c) != (void *) c) - s_free(MP_DIGITS(c)); - c->digits = out; - c->alloc = p; - } - - c->used = osize; /* might not be true, but we'll fix it ... */ - CLAMP(c); /* ... right here */ - c->sign = MP_ZPOS; - - return MP_OK; -} - -mp_result -mp_int_div(mp_int a, mp_int b, mp_int q, mp_int r) -{ - assert(a != NULL && b != NULL && q != r); - - int cmp; - mp_result res = MP_OK; - mp_int qout, - rout; - mp_sign sa = MP_SIGN(a); - mp_sign sb = MP_SIGN(b); - - if (CMPZ(b) == 0) - { - return MP_UNDEF; - } - else if ((cmp = s_ucmp(a, b)) < 0) - { - /* - * If |a| < |b|, no division is required: q = 0, r = a - */ - if (r && (res = mp_int_copy(a, r)) != MP_OK) - return res; - - if (q) - mp_int_zero(q); - - return MP_OK; - } - else if (cmp == 0) - { - /* - * If |a| = |b|, no division is required: q = 1 or -1, r = 0 - */ - if (r) - mp_int_zero(r); - - if (q) - { - mp_int_zero(q); - q->digits[0] = 1; - - if (sa != sb) - q->sign = MP_NEG; - } - - return MP_OK; - } - - /* - * When |a| > |b|, real division is required. We need someplace to store - * quotient and remainder, but q and r are allowed to be NULL or to - * overlap with the inputs. - */ - DECLARE_TEMP(2); - int lg; - - if ((lg = s_isp2(b)) < 0) - { - if (q && b != q) - { - REQUIRE(mp_int_copy(a, q)); - qout = q; - } - else - { - REQUIRE(mp_int_copy(a, TEMP(0))); - qout = TEMP(0); - } - - if (r && a != r) - { - REQUIRE(mp_int_copy(b, r)); - rout = r; - } - else - { - REQUIRE(mp_int_copy(b, TEMP(1))); - rout = TEMP(1); - } - - REQUIRE(s_udiv_knuth(qout, rout)); - } - else - { - if (q) - REQUIRE(mp_int_copy(a, q)); - if (r) - REQUIRE(mp_int_copy(a, r)); - - if (q) - s_qdiv(q, (mp_size) lg); - qout = q; - if (r) - s_qmod(r, (mp_size) lg); - rout = r; - } - - /* Recompute signs for output */ - if (rout) - { - rout->sign = sa; - if (CMPZ(rout) == 0) - rout->sign = MP_ZPOS; - } - if (qout) - { - qout->sign = (sa == sb) ? MP_ZPOS : MP_NEG; - if (CMPZ(qout) == 0) - qout->sign = MP_ZPOS; - } - - if (q) - REQUIRE(mp_int_copy(qout, q)); - if (r) - REQUIRE(mp_int_copy(rout, r)); - CLEANUP_TEMP(); - return res; -} - -mp_result -mp_int_mod(mp_int a, mp_int m, mp_int c) -{ - DECLARE_TEMP(1); - mp_int out = (m == c) ? TEMP(0) : c; - - REQUIRE(mp_int_div(a, m, NULL, out)); - if (CMPZ(out) < 0) - { - REQUIRE(mp_int_add(out, m, c)); - } - else - { - REQUIRE(mp_int_copy(out, c)); - } - CLEANUP_TEMP(); - return MP_OK; -} - -mp_result -mp_int_div_value(mp_int a, mp_small value, mp_int q, mp_small *r) -{ - mpz_t vtmp; - mp_digit vbuf[MP_VALUE_DIGITS(value)]; - - s_fake(&vtmp, value, vbuf); - - DECLARE_TEMP(1); - REQUIRE(mp_int_div(a, &vtmp, q, TEMP(0))); - - if (r) - (void) mp_int_to_int(TEMP(0), r); /* can't fail */ - - CLEANUP_TEMP(); - return MP_OK; -} - -mp_result -mp_int_div_pow2(mp_int a, mp_small p2, mp_int q, mp_int r) -{ - assert(a != NULL && p2 >= 0 && q != r); - - mp_result res = MP_OK; - - if (q != NULL && (res = mp_int_copy(a, q)) == MP_OK) - { - s_qdiv(q, (mp_size) p2); - } - - if (res == MP_OK && r != NULL && (res = mp_int_copy(a, r)) == MP_OK) - { - s_qmod(r, (mp_size) p2); - } - - return res; -} - -mp_result -mp_int_expt(mp_int a, mp_small b, mp_int c) -{ - assert(c != NULL); - if (b < 0) - return MP_RANGE; - - DECLARE_TEMP(1); - REQUIRE(mp_int_copy(a, TEMP(0))); - - (void) mp_int_set_value(c, 1); - unsigned int v = labs(b); - - while (v != 0) - { - if (v & 1) - { - REQUIRE(mp_int_mul(c, TEMP(0), c)); - } - - v >>= 1; - if (v == 0) - break; - - REQUIRE(mp_int_sqr(TEMP(0), TEMP(0))); - } - - CLEANUP_TEMP(); - return MP_OK; -} - -mp_result -mp_int_expt_value(mp_small a, mp_small b, mp_int c) -{ - assert(c != NULL); - if (b < 0) - return MP_RANGE; - - DECLARE_TEMP(1); - REQUIRE(mp_int_set_value(TEMP(0), a)); - - (void) mp_int_set_value(c, 1); - unsigned int v = labs(b); - - while (v != 0) - { - if (v & 1) - { - REQUIRE(mp_int_mul(c, TEMP(0), c)); - } - - v >>= 1; - if (v == 0) - break; - - REQUIRE(mp_int_sqr(TEMP(0), TEMP(0))); - } - - CLEANUP_TEMP(); - return MP_OK; -} - -mp_result -mp_int_expt_full(mp_int a, mp_int b, mp_int c) -{ - assert(a != NULL && b != NULL && c != NULL); - if (MP_SIGN(b) == MP_NEG) - return MP_RANGE; - - DECLARE_TEMP(1); - REQUIRE(mp_int_copy(a, TEMP(0))); - - (void) mp_int_set_value(c, 1); - for (unsigned ix = 0; ix < MP_USED(b); ++ix) - { - mp_digit d = b->digits[ix]; - - for (unsigned jx = 0; jx < MP_DIGIT_BIT; ++jx) - { - if (d & 1) - { - REQUIRE(mp_int_mul(c, TEMP(0), c)); - } - - d >>= 1; - if (d == 0 && ix + 1 == MP_USED(b)) - break; - REQUIRE(mp_int_sqr(TEMP(0), TEMP(0))); - } - } - - CLEANUP_TEMP(); - return MP_OK; -} - -int -mp_int_compare(mp_int a, mp_int b) -{ - assert(a != NULL && b != NULL); - - mp_sign sa = MP_SIGN(a); - - if (sa == MP_SIGN(b)) - { - int cmp = s_ucmp(a, b); - - /* - * If they're both zero or positive, the normal comparison applies; if - * both negative, the sense is reversed. - */ - if (sa == MP_ZPOS) - { - return cmp; - } - else - { - return -cmp; - } - } - else if (sa == MP_ZPOS) - { - return 1; - } - else - { - return -1; - } -} - -int -mp_int_compare_unsigned(mp_int a, mp_int b) -{ - assert(a != NULL && b != NULL); - - return s_ucmp(a, b); -} - -int -mp_int_compare_zero(mp_int z) -{ - assert(z != NULL); - - if (MP_USED(z) == 1 && z->digits[0] == 0) - { - return 0; - } - else if (MP_SIGN(z) == MP_ZPOS) - { - return 1; - } - else - { - return -1; - } -} - -int -mp_int_compare_value(mp_int z, mp_small value) -{ - assert(z != NULL); - - mp_sign vsign = (value < 0) ? MP_NEG : MP_ZPOS; - - if (vsign == MP_SIGN(z)) - { - int cmp = s_vcmp(z, value); - - return (vsign == MP_ZPOS) ? cmp : -cmp; - } - else - { - return (value < 0) ? 1 : -1; - } -} - -int -mp_int_compare_uvalue(mp_int z, mp_usmall uv) -{ - assert(z != NULL); - - if (MP_SIGN(z) == MP_NEG) - { - return -1; - } - else - { - return s_uvcmp(z, uv); - } -} - -mp_result -mp_int_exptmod(mp_int a, mp_int b, mp_int m, mp_int c) -{ - assert(a != NULL && b != NULL && c != NULL && m != NULL); - - /* Zero moduli and negative exponents are not considered. */ - if (CMPZ(m) == 0) - return MP_UNDEF; - if (CMPZ(b) < 0) - return MP_RANGE; - - mp_size um = MP_USED(m); - - DECLARE_TEMP(3); - REQUIRE(GROW(TEMP(0), 2 * um)); - REQUIRE(GROW(TEMP(1), 2 * um)); - - mp_int s; - - if (c == b || c == m) - { - REQUIRE(GROW(TEMP(2), 2 * um)); - s = TEMP(2); - } - else - { - s = c; - } - - REQUIRE(mp_int_mod(a, m, TEMP(0))); - REQUIRE(s_brmu(TEMP(1), m)); - REQUIRE(s_embar(TEMP(0), b, m, TEMP(1), s)); - REQUIRE(mp_int_copy(s, c)); - - CLEANUP_TEMP(); - return MP_OK; -} - -mp_result -mp_int_exptmod_evalue(mp_int a, mp_small value, mp_int m, mp_int c) -{ - mpz_t vtmp; - mp_digit vbuf[MP_VALUE_DIGITS(value)]; - - s_fake(&vtmp, value, vbuf); - - return mp_int_exptmod(a, &vtmp, m, c); -} - -mp_result -mp_int_exptmod_bvalue(mp_small value, mp_int b, mp_int m, mp_int c) -{ - mpz_t vtmp; - mp_digit vbuf[MP_VALUE_DIGITS(value)]; - - s_fake(&vtmp, value, vbuf); - - return mp_int_exptmod(&vtmp, b, m, c); -} - -mp_result -mp_int_exptmod_known(mp_int a, mp_int b, mp_int m, mp_int mu, - mp_int c) -{ - assert(a && b && m && c); - - /* Zero moduli and negative exponents are not considered. */ - if (CMPZ(m) == 0) - return MP_UNDEF; - if (CMPZ(b) < 0) - return MP_RANGE; - - DECLARE_TEMP(2); - mp_size um = MP_USED(m); - - REQUIRE(GROW(TEMP(0), 2 * um)); - - mp_int s; - - if (c == b || c == m) - { - REQUIRE(GROW(TEMP(1), 2 * um)); - s = TEMP(1); - } - else - { - s = c; - } - - REQUIRE(mp_int_mod(a, m, TEMP(0))); - REQUIRE(s_embar(TEMP(0), b, m, mu, s)); - REQUIRE(mp_int_copy(s, c)); - - CLEANUP_TEMP(); - return MP_OK; -} - -mp_result -mp_int_redux_const(mp_int m, mp_int c) -{ - assert(m != NULL && c != NULL && m != c); - - return s_brmu(c, m); -} - -mp_result -mp_int_invmod(mp_int a, mp_int m, mp_int c) -{ - assert(a != NULL && m != NULL && c != NULL); - - if (CMPZ(a) == 0 || CMPZ(m) <= 0) - return MP_RANGE; - - DECLARE_TEMP(2); - - REQUIRE(mp_int_egcd(a, m, TEMP(0), TEMP(1), NULL)); - - if (mp_int_compare_value(TEMP(0), 1) != 0) - { - REQUIRE(MP_UNDEF); - } - - /* It is first necessary to constrain the value to the proper range */ - REQUIRE(mp_int_mod(TEMP(1), m, TEMP(1))); - - /* - * Now, if 'a' was originally negative, the value we have is actually the - * magnitude of the negative representative; to get the positive value we - * have to subtract from the modulus. Otherwise, the value is okay as it - * stands. - */ - if (MP_SIGN(a) == MP_NEG) - { - REQUIRE(mp_int_sub(m, TEMP(1), c)); - } - else - { - REQUIRE(mp_int_copy(TEMP(1), c)); - } - - CLEANUP_TEMP(); - return MP_OK; -} - -/* Binary GCD algorithm due to Josef Stein, 1961 */ -mp_result -mp_int_gcd(mp_int a, mp_int b, mp_int c) -{ - assert(a != NULL && b != NULL && c != NULL); - - int ca = CMPZ(a); - int cb = CMPZ(b); - - if (ca == 0 && cb == 0) - { - return MP_UNDEF; - } - else if (ca == 0) - { - return mp_int_abs(b, c); - } - else if (cb == 0) - { - return mp_int_abs(a, c); - } - - DECLARE_TEMP(3); - REQUIRE(mp_int_copy(a, TEMP(0))); - REQUIRE(mp_int_copy(b, TEMP(1))); - - TEMP(0)->sign = MP_ZPOS; - TEMP(1)->sign = MP_ZPOS; - - int k = 0; - - { /* Divide out common factors of 2 from u and v */ - int div2_u = s_dp2k(TEMP(0)); - int div2_v = s_dp2k(TEMP(1)); - - k = MIN(div2_u, div2_v); - s_qdiv(TEMP(0), (mp_size) k); - s_qdiv(TEMP(1), (mp_size) k); - } - - if (mp_int_is_odd(TEMP(0))) - { - REQUIRE(mp_int_neg(TEMP(1), TEMP(2))); - } - else - { - REQUIRE(mp_int_copy(TEMP(0), TEMP(2))); - } - - for (;;) - { - s_qdiv(TEMP(2), s_dp2k(TEMP(2))); - - if (CMPZ(TEMP(2)) > 0) - { - REQUIRE(mp_int_copy(TEMP(2), TEMP(0))); - } - else - { - REQUIRE(mp_int_neg(TEMP(2), TEMP(1))); - } - - REQUIRE(mp_int_sub(TEMP(0), TEMP(1), TEMP(2))); - - if (CMPZ(TEMP(2)) == 0) - break; - } - - REQUIRE(mp_int_abs(TEMP(0), c)); - if (!s_qmul(c, (mp_size) k)) - REQUIRE(MP_MEMORY); - - CLEANUP_TEMP(); - return MP_OK; -} - -/* 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) -{ - assert(a != NULL && b != NULL && c != NULL && (x != NULL || y != NULL)); - - mp_result res = MP_OK; - int ca = CMPZ(a); - int cb = CMPZ(b); - - if (ca == 0 && cb == 0) - { - return MP_UNDEF; - } - else if (ca == 0) - { - if ((res = mp_int_abs(b, c)) != MP_OK) - return res; - mp_int_zero(x); - (void) mp_int_set_value(y, 1); - return MP_OK; - } - else if (cb == 0) - { - if ((res = mp_int_abs(a, c)) != MP_OK) - return res; - (void) mp_int_set_value(x, 1); - mp_int_zero(y); - return MP_OK; - } - - /* - * Initialize temporaries: A:0, B:1, C:2, D:3, u:4, v:5, ou:6, ov:7 - */ - DECLARE_TEMP(8); - REQUIRE(mp_int_set_value(TEMP(0), 1)); - REQUIRE(mp_int_set_value(TEMP(3), 1)); - REQUIRE(mp_int_copy(a, TEMP(4))); - REQUIRE(mp_int_copy(b, TEMP(5))); - - /* We will work with absolute values here */ - TEMP(4)->sign = MP_ZPOS; - TEMP(5)->sign = MP_ZPOS; - - int k = 0; - - { /* Divide out common factors of 2 from u and v */ - int div2_u = s_dp2k(TEMP(4)), - div2_v = s_dp2k(TEMP(5)); - - k = MIN(div2_u, div2_v); - s_qdiv(TEMP(4), k); - s_qdiv(TEMP(5), k); - } - - REQUIRE(mp_int_copy(TEMP(4), TEMP(6))); - REQUIRE(mp_int_copy(TEMP(5), TEMP(7))); - - for (;;) - { - while (mp_int_is_even(TEMP(4))) - { - s_qdiv(TEMP(4), 1); - - if (mp_int_is_odd(TEMP(0)) || mp_int_is_odd(TEMP(1))) - { - REQUIRE(mp_int_add(TEMP(0), TEMP(7), TEMP(0))); - REQUIRE(mp_int_sub(TEMP(1), TEMP(6), TEMP(1))); - } - - s_qdiv(TEMP(0), 1); - s_qdiv(TEMP(1), 1); - } - - while (mp_int_is_even(TEMP(5))) - { - s_qdiv(TEMP(5), 1); - - if (mp_int_is_odd(TEMP(2)) || mp_int_is_odd(TEMP(3))) - { - REQUIRE(mp_int_add(TEMP(2), TEMP(7), TEMP(2))); - REQUIRE(mp_int_sub(TEMP(3), TEMP(6), TEMP(3))); - } - - s_qdiv(TEMP(2), 1); - s_qdiv(TEMP(3), 1); - } - - if (mp_int_compare(TEMP(4), TEMP(5)) >= 0) - { - REQUIRE(mp_int_sub(TEMP(4), TEMP(5), TEMP(4))); - REQUIRE(mp_int_sub(TEMP(0), TEMP(2), TEMP(0))); - REQUIRE(mp_int_sub(TEMP(1), TEMP(3), TEMP(1))); - } - else - { - REQUIRE(mp_int_sub(TEMP(5), TEMP(4), TEMP(5))); - REQUIRE(mp_int_sub(TEMP(2), TEMP(0), TEMP(2))); - REQUIRE(mp_int_sub(TEMP(3), TEMP(1), TEMP(3))); - } - - if (CMPZ(TEMP(4)) == 0) - { - if (x) - REQUIRE(mp_int_copy(TEMP(2), x)); - if (y) - REQUIRE(mp_int_copy(TEMP(3), y)); - if (c) - { - if (!s_qmul(TEMP(5), k)) - { - REQUIRE(MP_MEMORY); - } - REQUIRE(mp_int_copy(TEMP(5), c)); - } - - break; - } - } - - CLEANUP_TEMP(); - return MP_OK; -} - -mp_result -mp_int_lcm(mp_int a, mp_int b, mp_int c) -{ - assert(a != NULL && b != NULL && c != NULL); - - /* - * Since a * b = gcd(a, b) * lcm(a, b), we can compute lcm(a, b) = (a / - * gcd(a, b)) * b. - * - * This formulation insures everything works even if the input variables - * share space. - */ - DECLARE_TEMP(1); - REQUIRE(mp_int_gcd(a, b, TEMP(0))); - REQUIRE(mp_int_div(a, TEMP(0), TEMP(0), NULL)); - REQUIRE(mp_int_mul(TEMP(0), b, TEMP(0))); - REQUIRE(mp_int_copy(TEMP(0), c)); - - CLEANUP_TEMP(); - return MP_OK; -} - -bool -mp_int_divisible_value(mp_int a, mp_small v) -{ - mp_small rem = 0; - - if (mp_int_div_value(a, v, NULL, &rem) != MP_OK) - { - return false; - } - return rem == 0; -} - -int -mp_int_is_pow2(mp_int z) -{ - assert(z != NULL); - - return s_isp2(z); -} - -/* Implementation of Newton's root finding method, based loosely on a patch - contributed by Hal Finkel <half@halssoftware.com> - modified by M. J. Fromberger. - */ -mp_result -mp_int_root(mp_int a, mp_small b, mp_int c) -{ - assert(a != NULL && c != NULL && b > 0); - - if (b == 1) - { - return mp_int_copy(a, c); - } - bool flips = false; - - if (MP_SIGN(a) == MP_NEG) - { - if (b % 2 == 0) - { - return MP_UNDEF; /* root does not exist for negative a with - * even b */ - } - else - { - flips = true; - } - } - - DECLARE_TEMP(5); - REQUIRE(mp_int_copy(a, TEMP(0))); - REQUIRE(mp_int_copy(a, TEMP(1))); - TEMP(0)->sign = MP_ZPOS; - TEMP(1)->sign = MP_ZPOS; - - for (;;) - { - REQUIRE(mp_int_expt(TEMP(1), b, TEMP(2))); - - if (mp_int_compare_unsigned(TEMP(2), TEMP(0)) <= 0) - break; - - REQUIRE(mp_int_sub(TEMP(2), TEMP(0), TEMP(2))); - REQUIRE(mp_int_expt(TEMP(1), b - 1, TEMP(3))); - REQUIRE(mp_int_mul_value(TEMP(3), b, TEMP(3))); - REQUIRE(mp_int_div(TEMP(2), TEMP(3), TEMP(4), NULL)); - REQUIRE(mp_int_sub(TEMP(1), TEMP(4), TEMP(4))); - - if (mp_int_compare_unsigned(TEMP(1), TEMP(4)) == 0) - { - REQUIRE(mp_int_sub_value(TEMP(4), 1, TEMP(4))); - } - REQUIRE(mp_int_copy(TEMP(4), TEMP(1))); - } - - REQUIRE(mp_int_copy(TEMP(1), c)); - - /* If the original value of a was negative, flip the output sign. */ - if (flips) - (void) mp_int_neg(c, c); /* cannot fail */ - - CLEANUP_TEMP(); - return MP_OK; -} - -mp_result -mp_int_to_int(mp_int z, mp_small *out) -{ - assert(z != NULL); - - /* Make sure the value is representable as a small integer */ - mp_sign sz = MP_SIGN(z); - - if ((sz == MP_ZPOS && mp_int_compare_value(z, MP_SMALL_MAX) > 0) || - mp_int_compare_value(z, MP_SMALL_MIN) < 0) - { - return MP_RANGE; - } - - mp_usmall uz = MP_USED(z); - mp_digit *dz = MP_DIGITS(z) + uz - 1; - mp_small uv = 0; - - while (uz > 0) - { - uv <<= MP_DIGIT_BIT / 2; - uv = (uv << (MP_DIGIT_BIT / 2)) | *dz--; - --uz; - } - - if (out) - *out = (mp_small) ((sz == MP_NEG) ? -uv : uv); - - return MP_OK; -} - -mp_result -mp_int_to_uint(mp_int z, mp_usmall *out) -{ - assert(z != NULL); - - /* Make sure the value is representable as an unsigned small integer */ - mp_size sz = MP_SIGN(z); - - if (sz == MP_NEG || mp_int_compare_uvalue(z, MP_USMALL_MAX) > 0) - { - return MP_RANGE; - } - - mp_size uz = MP_USED(z); - mp_digit *dz = MP_DIGITS(z) + uz - 1; - mp_usmall uv = 0; - - while (uz > 0) - { - uv <<= MP_DIGIT_BIT / 2; - uv = (uv << (MP_DIGIT_BIT / 2)) | *dz--; - --uz; - } - - if (out) - *out = uv; - - return MP_OK; -} - -mp_result -mp_int_to_string(mp_int z, mp_size radix, char *str, int limit) -{ - assert(z != NULL && str != NULL && limit >= 2); - assert(radix >= MP_MIN_RADIX && radix <= MP_MAX_RADIX); - - int cmp = 0; - - if (CMPZ(z) == 0) - { - *str++ = s_val2ch(0, 1); - } - else - { - mp_result res; - mpz_t tmp; - char *h, - *t; - - if ((res = mp_int_init_copy(&tmp, z)) != MP_OK) - return res; - - if (MP_SIGN(z) == MP_NEG) - { - *str++ = '-'; - --limit; - } - h = str; - - /* Generate digits in reverse order until finished or limit reached */ - for ( /* */ ; limit > 0; --limit) - { - mp_digit d; - - if ((cmp = CMPZ(&tmp)) == 0) - break; - - d = s_ddiv(&tmp, (mp_digit) radix); - *str++ = s_val2ch(d, 1); - } - t = str - 1; - - /* Put digits back in correct output order */ - while (h < t) - { - char tc = *h; - - *h++ = *t; - *t-- = tc; - } - - mp_int_clear(&tmp); - } - - *str = '\0'; - if (cmp == 0) - { - return MP_OK; - } - else - { - return MP_TRUNC; - } -} - -mp_result -mp_int_string_len(mp_int z, mp_size radix) -{ - assert(z != NULL); - assert(radix >= MP_MIN_RADIX && radix <= MP_MAX_RADIX); - - int len = s_outlen(z, radix) + 1; /* for terminator */ - - /* Allow for sign marker on negatives */ - if (MP_SIGN(z) == MP_NEG) - len += 1; - - return len; -} - -/* Read zero-terminated string into z */ -mp_result -mp_int_read_string(mp_int z, mp_size radix, const char *str) -{ - return mp_int_read_cstring(z, radix, str, NULL); -} - -mp_result -mp_int_read_cstring(mp_int z, mp_size radix, const char *str, - char **end) -{ - assert(z != NULL && str != NULL); - assert(radix >= MP_MIN_RADIX && radix <= MP_MAX_RADIX); - - /* Skip leading whitespace */ - while (isspace((unsigned char) *str)) - ++str; - - /* Handle leading sign tag (+/-, positive default) */ - switch (*str) - { - case '-': - z->sign = MP_NEG; - ++str; - break; - case '+': - ++str; /* fallthrough */ - default: - z->sign = MP_ZPOS; - break; - } - - /* Skip leading zeroes */ - int ch; - - while ((ch = s_ch2val(*str, radix)) == 0) - ++str; - - /* Make sure there is enough space for the value */ - if (!s_pad(z, s_inlen(strlen(str), radix))) - return MP_MEMORY; - - z->used = 1; - z->digits[0] = 0; - - while (*str != '\0' && ((ch = s_ch2val(*str, radix)) >= 0)) - { - s_dmul(z, (mp_digit) radix); - s_dadd(z, (mp_digit) ch); - ++str; - } - - CLAMP(z); - - /* Override sign for zero, even if negative specified. */ - if (CMPZ(z) == 0) - z->sign = MP_ZPOS; - - if (end != NULL) - *end = unconstify(char *, str); - - /* - * Return a truncation error if the string has unprocessed characters - * remaining, so the caller can tell if the whole string was done - */ - if (*str != '\0') - { - return MP_TRUNC; - } - else - { - return MP_OK; - } -} - -mp_result -mp_int_count_bits(mp_int z) -{ - assert(z != NULL); - - mp_size uz = MP_USED(z); - - if (uz == 1 && z->digits[0] == 0) - return 1; - - --uz; - mp_size nbits = uz * MP_DIGIT_BIT; - mp_digit d = z->digits[uz]; - - while (d != 0) - { - d >>= 1; - ++nbits; - } - - return nbits; -} - -mp_result -mp_int_to_binary(mp_int z, unsigned char *buf, int limit) -{ - static const int PAD_FOR_2C = 1; - - assert(z != NULL && buf != NULL); - - int limpos = limit; - mp_result res = s_tobin(z, buf, &limpos, PAD_FOR_2C); - - if (MP_SIGN(z) == MP_NEG) - s_2comp(buf, limpos); - - return res; -} - -mp_result -mp_int_read_binary(mp_int z, unsigned char *buf, int len) -{ - assert(z != NULL && buf != NULL && len > 0); - - /* Figure out how many digits are needed to represent this value */ - mp_size need = ((len * CHAR_BIT) + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT; - - if (!s_pad(z, need)) - return MP_MEMORY; - - mp_int_zero(z); - - /* - * If the high-order bit is set, take the 2's complement before reading - * the value (it will be restored afterward) - */ - if (buf[0] >> (CHAR_BIT - 1)) - { - z->sign = MP_NEG; - s_2comp(buf, len); - } - - mp_digit *dz = MP_DIGITS(z); - unsigned char *tmp = buf; - - for (int i = len; i > 0; --i, ++tmp) - { - s_qmul(z, (mp_size) CHAR_BIT); - *dz |= *tmp; - } - - /* Restore 2's complement if we took it before */ - if (MP_SIGN(z) == MP_NEG) - s_2comp(buf, len); - - return MP_OK; -} - -mp_result -mp_int_binary_len(mp_int z) -{ - mp_result res = mp_int_count_bits(z); - - if (res <= 0) - return res; - - int bytes = mp_int_unsigned_len(z); - - /* - * If the highest-order bit falls exactly on a byte boundary, we need to - * pad with an extra byte so that the sign will be read correctly when - * reading it back in. - */ - if (bytes * CHAR_BIT == res) - ++bytes; - - return bytes; -} - -mp_result -mp_int_to_unsigned(mp_int z, unsigned char *buf, int limit) -{ - static const int NO_PADDING = 0; - - assert(z != NULL && buf != NULL); - - return s_tobin(z, buf, &limit, NO_PADDING); -} - -mp_result -mp_int_read_unsigned(mp_int z, unsigned char *buf, int len) -{ - assert(z != NULL && buf != NULL && len > 0); - - /* Figure out how many digits are needed to represent this value */ - mp_size need = ((len * CHAR_BIT) + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT; - - if (!s_pad(z, need)) - return MP_MEMORY; - - mp_int_zero(z); - - unsigned char *tmp = buf; - - for (int i = len; i > 0; --i, ++tmp) - { - (void) s_qmul(z, CHAR_BIT); - *MP_DIGITS(z) |= *tmp; - } - - return MP_OK; -} - -mp_result -mp_int_unsigned_len(mp_int z) -{ - mp_result res = mp_int_count_bits(z); - - if (res <= 0) - return res; - - int bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT; - - return bytes; -} - -const char * -mp_error_string(mp_result res) -{ - if (res > 0) - return s_unknown_err; - - res = -res; - int ix; - - for (ix = 0; ix < res && s_error_msg[ix] != NULL; ++ix) - ; - - if (s_error_msg[ix] != NULL) - { - return s_error_msg[ix]; - } - else - { - return s_unknown_err; - } -} - -/*------------------------------------------------------------------------*/ -/* Private functions for internal use. These make assumptions. */ - -#if IMATH_DEBUG -static const mp_digit fill = (mp_digit) 0xdeadbeefabad1dea; -#endif - -static mp_digit * -s_alloc(mp_size num) -{ - mp_digit *out = palloc(num * sizeof(mp_digit)); - - assert(out != NULL); - -#if IMATH_DEBUG - for (mp_size ix = 0; ix < num; ++ix) - out[ix] = fill; -#endif - return out; -} - -static mp_digit * -s_realloc(mp_digit *old, mp_size osize, mp_size nsize) -{ -#if IMATH_DEBUG - mp_digit *new = s_alloc(nsize); - - assert(new != NULL); - - for (mp_size ix = 0; ix < nsize; ++ix) - new[ix] = fill; - memcpy(new, old, osize * sizeof(mp_digit)); -#else - mp_digit *new = repalloc(old, nsize * sizeof(mp_digit)); - - assert(new != NULL); -#endif - - return new; -} - -static void -s_free(void *ptr) -{ - pfree(ptr); -} - -static bool -s_pad(mp_int z, mp_size min) -{ - if (MP_ALLOC(z) < min) - { - mp_size nsize = s_round_prec(min); - mp_digit *tmp; - - if (z->digits == &(z->single)) - { - if ((tmp = s_alloc(nsize)) == NULL) - return false; - tmp[0] = z->single; - } - else if ((tmp = s_realloc(MP_DIGITS(z), MP_ALLOC(z), nsize)) == NULL) - { - return false; - } - - z->digits = tmp; - z->alloc = nsize; - } - - return true; -} - -/* Note: This will not work correctly when value == MP_SMALL_MIN */ -static void -s_fake(mp_int z, mp_small value, mp_digit vbuf[]) -{ - mp_usmall uv = (mp_usmall) (value < 0) ? -value : value; - - s_ufake(z, uv, vbuf); - if (value < 0) - z->sign = MP_NEG; -} - -static void -s_ufake(mp_int z, mp_usmall value, mp_digit vbuf[]) -{ - mp_size ndig = (mp_size) s_uvpack(value, vbuf); - - z->used = ndig; - z->alloc = MP_VALUE_DIGITS(value); - z->sign = MP_ZPOS; - z->digits = vbuf; -} - -static int -s_cdig(mp_digit *da, mp_digit *db, mp_size len) -{ - mp_digit *dat = da + len - 1, - *dbt = db + len - 1; - - for ( /* */ ; len != 0; --len, --dat, --dbt) - { - if (*dat > *dbt) - { - return 1; - } - else if (*dat < *dbt) - { - return -1; - } - } - - return 0; -} - -static int -s_uvpack(mp_usmall uv, mp_digit t[]) -{ - int ndig = 0; - - if (uv == 0) - t[ndig++] = 0; - else - { - while (uv != 0) - { - t[ndig++] = (mp_digit) uv; - uv >>= MP_DIGIT_BIT / 2; - uv >>= MP_DIGIT_BIT / 2; - } - } - - return ndig; -} - -static int -s_ucmp(mp_int a, mp_int b) -{ - mp_size ua = MP_USED(a), - ub = MP_USED(b); - - if (ua > ub) - { - return 1; - } - else if (ub > ua) - { - return -1; - } - else - { - return s_cdig(MP_DIGITS(a), MP_DIGITS(b), ua); - } -} - -static int -s_vcmp(mp_int a, mp_small v) -{ -#ifdef _MSC_VER -#pragma warning(push) -#pragma warning(disable: 4146) -#endif - mp_usmall uv = (v < 0) ? -(mp_usmall) v : (mp_usmall) v; -#ifdef _MSC_VER -#pragma warning(pop) -#endif - - return s_uvcmp(a, uv); -} - -static int -s_uvcmp(mp_int a, mp_usmall uv) -{ - mpz_t vtmp; - mp_digit vdig[MP_VALUE_DIGITS(uv)]; - - s_ufake(&vtmp, uv, vdig); - return s_ucmp(a, &vtmp); -} - -static mp_digit -s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, - mp_size size_b) -{ - mp_size pos; - mp_word w = 0; - - /* Insure that da is the longer of the two to simplify later code */ - if (size_b > size_a) - { - SWAP(mp_digit *, da, db); - SWAP(mp_size, size_a, size_b); - } - - /* Add corresponding digits until the shorter number runs out */ - for (pos = 0; pos < size_b; ++pos, ++da, ++db, ++dc) - { - w = w + (mp_word) *da + (mp_word) *db; - *dc = LOWER_HALF(w); - w = UPPER_HALF(w); - } - - /* Propagate carries as far as necessary */ - for ( /* */ ; pos < size_a; ++pos, ++da, ++dc) - { - w = w + *da; - - *dc = LOWER_HALF(w); - w = UPPER_HALF(w); - } - - /* Return carry out */ - return (mp_digit) w; -} - -static void -s_usub(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, - mp_size size_b) -{ - mp_size pos; - mp_word w = 0; - - /* We assume that |a| >= |b| so this should definitely hold */ - assert(size_a >= size_b); - - /* Subtract corresponding digits and propagate borrow */ - for (pos = 0; pos < size_b; ++pos, ++da, ++db, ++dc) - { - w = ((mp_word) MP_DIGIT_MAX + 1 + /* MP_RADIX */ - (mp_word) *da) - - w - (mp_word) *db; - - *dc = LOWER_HALF(w); - w = (UPPER_HALF(w) == 0); - } - - /* Finish the subtraction for remaining upper digits of da */ - for ( /* */ ; pos < size_a; ++pos, ++da, ++dc) - { - w = ((mp_word) MP_DIGIT_MAX + 1 + /* MP_RADIX */ - (mp_word) *da) - - w; - - *dc = LOWER_HALF(w); - w = (UPPER_HALF(w) == 0); - } - - /* If there is a borrow out at the end, it violates the precondition */ - assert(w == 0); -} - -static int -s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, - mp_size size_b) -{ - mp_size bot_size; - - /* Make sure b is the smaller of the two input values */ - if (size_b > size_a) - { - SWAP(mp_digit *, da, db); - SWAP(mp_size, size_a, size_b); - } - - /* - * Insure that the bottom is the larger half in an odd-length split; the - * code below relies on this being true. - */ - bot_size = (size_a + 1) / 2; - - /* - * If the values are big enough to bother with recursion, use the - * Karatsuba algorithm to compute the product; otherwise use the normal - * multiplication algorithm - */ - if (multiply_threshold && size_a >= multiply_threshold && size_b > bot_size) - { - mp_digit *t1, - *t2, - *t3, - carry; - - mp_digit *a_top = da + bot_size; - mp_digit *b_top = db + bot_size; - - mp_size at_size = size_a - bot_size; - mp_size bt_size = size_b - bot_size; - mp_size buf_size = 2 * bot_size; - - /* - * Do a single allocation for all three temporary buffers needed; each - * buffer must be big enough to hold the product of two bottom halves, - * and one buffer needs space for the completed product; twice the - * space is plenty. - */ - if ((t1 = s_alloc(4 * buf_size)) == NULL) - return 0; - t2 = t1 + buf_size; - t3 = t2 + buf_size; - ZERO(t1, 4 * buf_size); - - /* - * t1 and t2 are initially used as temporaries to compute the inner - * product (a1 + a0)(b1 + b0) = a1b1 + a1b0 + a0b1 + a0b0 - */ - carry = s_uadd(da, a_top, t1, bot_size, at_size); /* t1 = a1 + a0 */ - t1[bot_size] = carry; - - carry = s_uadd(db, b_top, t2, bot_size, bt_size); /* t2 = b1 + b0 */ - t2[bot_size] = carry; - - (void) s_kmul(t1, t2, t3, bot_size + 1, bot_size + 1); /* t3 = t1 * t2 */ - - /* - * Now we'll get t1 = a0b0 and t2 = a1b1, and subtract them out so - * that we're left with only the pieces we want: t3 = a1b0 + a0b1 - */ - ZERO(t1, buf_size); - ZERO(t2, buf_size); - (void) s_kmul(da, db, t1, bot_size, bot_size); /* t1 = a0 * b0 */ - (void) s_kmul(a_top, b_top, t2, at_size, bt_size); /* t2 = a1 * b1 */ - - /* Subtract out t1 and t2 to get the inner product */ - s_usub(t3, t1, t3, buf_size + 2, buf_size); - s_usub(t3, t2, t3, buf_size + 2, buf_size); - - /* Assemble the output value */ - COPY(t1, dc, buf_size); - carry = s_uadd(t3, dc + bot_size, dc + bot_size, buf_size + 1, buf_size); - assert(carry == 0); - - carry = - s_uadd(t2, dc + 2 * bot_size, dc + 2 * bot_size, buf_size, buf_size); - assert(carry == 0); - - s_free(t1); /* note t2 and t3 are just internal pointers - * to t1 */ - } - else - { - s_umul(da, db, dc, size_a, size_b); - } - - return 1; -} - -static void -s_umul(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, - mp_size size_b) -{ - mp_size a, - b; - mp_word w; - - for (a = 0; a < size_a; ++a, ++dc, ++da) - { - mp_digit *dct = dc; - mp_digit *dbt = db; - - if (*da == 0) - continue; - - w = 0; - for (b = 0; b < size_b; ++b, ++dbt, ++dct) - { - w = (mp_word) *da * (mp_word) *dbt + w + (mp_word) *dct; - - *dct = LOWER_HALF(w); - w = UPPER_HALF(w); - } - - *dct = (mp_digit) w; - } -} - -static int -s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a) -{ - if (multiply_threshold && size_a > multiply_threshold) - { - mp_size bot_size = (size_a + 1) / 2; - mp_digit *a_top = da + bot_size; - mp_digit *t1, - *t2, - *t3, - carry PG_USED_FOR_ASSERTS_ONLY; - mp_size at_size = size_a - bot_size; - mp_size buf_size = 2 * bot_size; - - if ((t1 = s_alloc(4 * buf_size)) == NULL) - return 0; - t2 = t1 + buf_size; - t3 = t2 + buf_size; - ZERO(t1, 4 * buf_size); - - (void) s_ksqr(da, t1, bot_size); /* t1 = a0 ^ 2 */ - (void) s_ksqr(a_top, t2, at_size); /* t2 = a1 ^ 2 */ - - (void) s_kmul(da, a_top, t3, bot_size, at_size); /* t3 = a0 * a1 */ - - /* Quick multiply t3 by 2, shifting left (can't overflow) */ - { - int i, - top = bot_size + at_size; - mp_word w, - save = 0; - - for (i = 0; i < top; ++i) - { - w = t3[i]; - w = (w << 1) | save; - t3[i] = LOWER_HALF(w); - save = UPPER_HALF(w); - } - t3[i] = LOWER_HALF(save); - } - - /* Assemble the output value */ - COPY(t1, dc, 2 * bot_size); - carry = s_uadd(t3, dc + bot_size, dc + bot_size, buf_size + 1, buf_size); - assert(carry == 0); - - carry = - s_uadd(t2, dc + 2 * bot_size, dc + 2 * bot_size, buf_size, buf_size); - assert(carry == 0); - - s_free(t1); /* note that t2 and t2 are internal pointers - * only */ - - } - else - { - s_usqr(da, dc, size_a); - } - - return 1; -} - -static void -s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a) -{ - mp_size i, - j; - mp_word w; - - for (i = 0; i < size_a; ++i, dc += 2, ++da) - { - mp_digit *dct = dc, - *dat = da; - - if (*da == 0) - continue; - - /* Take care of the first digit, no rollover */ - w = (mp_word) *dat * (mp_word) *dat + (mp_word) *dct; - *dct = LOWER_HALF(w); - w = UPPER_HALF(w); - ++dat; - ++dct; - - for (j = i + 1; j < size_a; ++j, ++dat, ++dct) - { - mp_word t = (mp_word) *da * (mp_word) *dat; - mp_word u = w + (mp_word) *dct, - ov = 0; - - /* Check if doubling t will overflow a word */ - if (HIGH_BIT_SET(t)) - ov = 1; - - w = t + t; - - /* Check if adding u to w will overflow a word */ - if (ADD_WILL_OVERFLOW(w, u)) - ov = 1; - - w += u; - - *dct = LOWER_HALF(w); - w = UPPER_HALF(w); - if (ov) - { - w += MP_DIGIT_MAX; /* MP_RADIX */ - ++w; - } - } - - w = w + *dct; - *dct = (mp_digit) w; - while ((w = UPPER_HALF(w)) != 0) - { - ++dct; - w = w + *dct; - *dct = LOWER_HALF(w); - } - - assert(w == 0); - } -} - -static void -s_dadd(mp_int a, mp_digit b) -{ - mp_word w = 0; - mp_digit *da = MP_DIGITS(a); - mp_size ua = MP_USED(a); - - w = (mp_word) *da + b; - *da++ = LOWER_HALF(w); - w = UPPER_HALF(w); - - for (ua -= 1; ua > 0; --ua, ++da) - { - w = (mp_word) *da + w; - - *da = LOWER_HALF(w); - w = UPPER_HALF(w); - } - - if (w) - { - *da = (mp_digit) w; - a->used += 1; - } -} - -static void -s_dmul(mp_int a, mp_digit b) -{ - mp_word w = 0; - mp_digit *da = MP_DIGITS(a); - mp_size ua = MP_USED(a); - - while (ua > 0) - { - w = (mp_word) *da * b + w; - *da++ = LOWER_HALF(w); - w = UPPER_HALF(w); - --ua; - } - - if (w) - { - *da = (mp_digit) w; - a->used += 1; - } -} - -static void -s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a) -{ - mp_word w = 0; - - while (size_a > 0) - { - w = (mp_word) *da++ * (mp_word) b + w; - - *dc++ = LOWER_HALF(w); - w = UPPER_HALF(w); - --size_a; - } - - if (w) - *dc = LOWER_HALF(w); -} - -static mp_digit -s_ddiv(mp_int a, mp_digit b) -{ - mp_word w = 0, - qdigit; - mp_size ua = MP_USED(a); - mp_digit *da = MP_DIGITS(a) + ua - 1; - - for ( /* */ ; ua > 0; --ua, --da) - { - w = (w << MP_DIGIT_BIT) | *da; - - if (w >= b) - { - qdigit = w / b; - w = w % b; - } - else - { - qdigit = 0; - } - - *da = (mp_digit) qdigit; - } - - CLAMP(a); - return (mp_digit) w; -} - -static void -s_qdiv(mp_int z, mp_size p2) -{ - mp_size ndig = p2 / MP_DIGIT_BIT, - nbits = p2 % MP_DIGIT_BIT; - mp_size uz = MP_USED(z); - - if (ndig) - { - mp_size mark; - mp_digit *to, - *from; - - if (ndig >= uz) - { - mp_int_zero(z); - return; - } - - to = MP_DIGITS(z); - from = to + ndig; - - for (mark = ndig; mark < uz; ++mark) - { - *to++ = *from++; - } - - z->used = uz - ndig; - } - - if (nbits) - { - mp_digit d = 0, - *dz, - save; - mp_size up = MP_DIGIT_BIT - nbits; - - uz = MP_USED(z); - dz = MP_DIGITS(z) + uz - 1; - - for ( /* */ ; uz > 0; --uz, --dz) - { - save = *dz; - - *dz = (*dz >> nbits) | (d << up); - d = save; - } - - CLAMP(z); - } - - if (MP_USED(z) == 1 && z->digits[0] == 0) - z->sign = MP_ZPOS; -} - -static void -s_qmod(mp_int z, mp_size p2) -{ - mp_size start = p2 / MP_DIGIT_BIT + 1, - rest = p2 % MP_DIGIT_BIT; - mp_size uz = MP_USED(z); - mp_digit mask = (1u << rest) - 1; - - if (start <= uz) - { - z->used = start; - z->digits[start - 1] &= mask; - CLAMP(z); - } -} - -static int -s_qmul(mp_int z, mp_size p2) -{ - mp_size uz, - need, - rest, - extra, - i; - mp_digit *from, - *to, - d; - - if (p2 == 0) - return 1; - - uz = MP_USED(z); - need = p2 / MP_DIGIT_BIT; - rest = p2 % MP_DIGIT_BIT; - - /* - * Figure out if we need an extra digit at the top end; this occurs if the - * topmost `rest' bits of the high-order digit of z are not zero, meaning - * they will be shifted off the end if not preserved - */ - extra = 0; - if (rest != 0) - { - mp_digit *dz = MP_DIGITS(z) + uz - 1; - - if ((*dz >> (MP_DIGIT_BIT - rest)) != 0) - extra = 1; - } - - if (!s_pad(z, uz + need + extra)) - return 0; - - /* - * If we need to shift by whole digits, do that in one pass, then to back - * and shift by partial digits. - */ - if (need > 0) - { - from = MP_DIGITS(z) + uz - 1; - to = from + need; - - for (i = 0; i < uz; ++i) - *to-- = *from--; - - ZERO(MP_DIGITS(z), need); - uz += need; - } - - if (rest) - { - d = 0; - for (i = need, from = MP_DIGITS(z) + need; i < uz; ++i, ++from) - { - mp_digit save = *from; - - *from = (*from << rest) | (d >> (MP_DIGIT_BIT - rest)); - d = save; - } - - d >>= (MP_DIGIT_BIT - rest); - if (d != 0) - { - *from = d; - uz += extra; - } - } - - z->used = uz; - CLAMP(z); - - return 1; -} - -/* 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) -{ - mp_digit hi = (1u << (p2 % MP_DIGIT_BIT)), - *zp; - mp_size tdig = (p2 / MP_DIGIT_BIT), - pos; - mp_word w = 0; - - if (!s_pad(z, tdig + 1)) - return 0; - - for (pos = 0, zp = MP_DIGITS(z); pos < tdig; ++pos, ++zp) - { - w = ((mp_word) MP_DIGIT_MAX + 1) - w - (mp_word) *zp; - - *zp = LOWER_HALF(w); - w = UPPER_HALF(w) ? 0 : 1; - } - - w = ((mp_word) MP_DIGIT_MAX + 1 + hi) - w - (mp_word) *zp; - *zp = LOWER_HALF(w); - - assert(UPPER_HALF(w) != 0); /* no borrow out should be possible */ - - z->sign = MP_ZPOS; - CLAMP(z); - - return 1; -} - -static int -s_dp2k(mp_int z) -{ - int k = 0; - mp_digit *dp = MP_DIGITS(z), - d; - - if (MP_USED(z) == 1 && *dp == 0) - return 1; - - while (*dp == 0) - { - k += MP_DIGIT_BIT; - ++dp; - } - - d = *dp; - while ((d & 1) == 0) - { - d >>= 1; - ++k; - } - - return k; -} - -static int -s_isp2(mp_int z) -{ - mp_size uz = MP_USED(z), - k = 0; - mp_digit *dz = MP_DIGITS(z), - d; - - while (uz > 1) - { - if (*dz++ != 0) - return -1; - k += MP_DIGIT_BIT; - --uz; - } - - d = *dz; - while (d > 1) - { - if (d & 1) - return -1; - ++k; - d >>= 1; - } - - return (int) k; -} - -static int -s_2expt(mp_int z, mp_small k) -{ - mp_size ndig, - rest; - mp_digit *dz; - - ndig = (k + MP_DIGIT_BIT) / MP_DIGIT_BIT; - rest = k % MP_DIGIT_BIT; - - if (!s_pad(z, ndig)) - return 0; - - dz = MP_DIGITS(z); - ZERO(dz, ndig); - *(dz + ndig - 1) = (1u << rest); - z->used = ndig; - - return 1; -} - -static int -s_norm(mp_int a, mp_int b) -{ - mp_digit d = b->digits[MP_USED(b) - 1]; - int k = 0; - - while (d < (1u << (mp_digit) (MP_DIGIT_BIT - 1))) - { /* d < (MP_RADIX / 2) */ - d <<= 1; - ++k; - } - - /* These multiplications can't fail */ - if (k != 0) - { - (void) s_qmul(a, (mp_size) k); - (void) s_qmul(b, (mp_size) k); - } - - return k; -} - -static mp_result -s_brmu(mp_int z, mp_int m) -{ - mp_size um = MP_USED(m) * 2; - - if (!s_pad(z, um)) - return MP_MEMORY; - - s_2expt(z, MP_DIGIT_BIT * um); - return mp_int_div(z, m, z, NULL); -} - -static int -s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2) -{ - mp_size um = MP_USED(m), - umb_p1, - umb_m1; - - umb_p1 = (um + 1) * MP_DIGIT_BIT; - umb_m1 = (um - 1) * MP_DIGIT_BIT; - - if (mp_int_copy(x, q1) != MP_OK) - return 0; - - /* Compute q2 = floor((floor(x / b^(k-1)) * mu) / b^(k+1)) */ - s_qdiv(q1, umb_m1); - UMUL(q1, mu, q2); - s_qdiv(q2, umb_p1); - - /* Set x = x mod b^(k+1) */ - s_qmod(x, umb_p1); - - /* - * Now, q is a guess for the quotient a / m. Compute x - q * m mod - * b^(k+1), replacing x. This may be off by a factor of 2m, but no more - * than that. - */ - UMUL(q2, m, q1); - s_qmod(q1, umb_p1); - (void) mp_int_sub(x, q1, x); /* can't fail */ - - /* - * The result may be < 0; if it is, add b^(k+1) to pin it in the proper - * range. - */ - if ((CMPZ(x) < 0) && !s_qsub(x, umb_p1)) - return 0; - - /* - * If x > m, we need to back it off until it is in range. This will be - * required at most twice. - */ - if (mp_int_compare(x, m) >= 0) - { - (void) mp_int_sub(x, m, x); - if (mp_int_compare(x, m) >= 0) - { - (void) mp_int_sub(x, m, x); - } - } - - /* At this point, x has been properly reduced. */ - return 1; -} - -/* 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) -{ - mp_digit umu = MP_USED(mu); - mp_digit *db = MP_DIGITS(b); - mp_digit *dbt = db + MP_USED(b) - 1; - - DECLARE_TEMP(3); - REQUIRE(GROW(TEMP(0), 4 * umu)); - REQUIRE(GROW(TEMP(1), 4 * umu)); - REQUIRE(GROW(TEMP(2), 4 * umu)); - ZERO(TEMP(0)->digits, TEMP(0)->alloc); - ZERO(TEMP(1)->digits, TEMP(1)->alloc); - ZERO(TEMP(2)->digits, TEMP(2)->alloc); - - (void) mp_int_set_value(c, 1); - - /* Take care of low-order digits */ - while (db < dbt) - { - mp_digit d = *db; - - for (int i = MP_DIGIT_BIT; i > 0; --i, d >>= 1) - { - if (d & 1) - { - /* The use of a second temporary avoids allocation */ - UMUL(c, a, TEMP(0)); - if (!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) - { - REQUIRE(MP_MEMORY); - } - mp_int_copy(TEMP(0), c); - } - - USQR(a, TEMP(0)); - assert(MP_SIGN(TEMP(0)) == MP_ZPOS); - if (!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) - { - REQUIRE(MP_MEMORY); - } - assert(MP_SIGN(TEMP(0)) == MP_ZPOS); - mp_int_copy(TEMP(0), a); - } - - ++db; - } - - /* Take care of highest-order digit */ - mp_digit d = *dbt; - - for (;;) - { - if (d & 1) - { - UMUL(c, a, TEMP(0)); - if (!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) - { - REQUIRE(MP_MEMORY); - } - mp_int_copy(TEMP(0), c); - } - - d >>= 1; - if (!d) - break; - - USQR(a, TEMP(0)); - if (!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) - { - REQUIRE(MP_MEMORY); - } - (void) mp_int_copy(TEMP(0), a); - } - - CLEANUP_TEMP(); - return MP_OK; -} - -/* 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) -{ - /* Force signs to positive */ - u->sign = MP_ZPOS; - v->sign = MP_ZPOS; - - /* Use simple division algorithm when v is only one digit long */ - if (MP_USED(v) == 1) - { - mp_digit d, - rem; - - d = v->digits[0]; - rem = s_ddiv(u, d); - mp_int_set_value(v, rem); - return MP_OK; - } - - /* - * Algorithm D - * - * The n and m variables are defined as used by Knuth. u is an n digit - * number with digits u_{n-1}..u_0. v is an n+m digit number with digits - * from v_{m+n-1}..v_0. We require that n > 1 and m >= 0 - */ - mp_size n = MP_USED(v); - mp_size m = MP_USED(u) - n; - - assert(n > 1); - /* assert(m >= 0) follows because m is unsigned. */ - - /* - * D1: Normalize. The normalization step provides the necessary condition - * for Theorem B, which states that the quotient estimate for q_j, call it - * qhat - * - * qhat = u_{j+n}u_{j+n-1} / v_{n-1} - * - * is bounded by - * - * qhat - 2 <= q_j <= qhat. - * - * That is, qhat is always greater than the actual quotient digit q, and - * it is never more than two larger than the actual quotient digit. - */ - int k = s_norm(u, v); - - /* - * Extend size of u by one if needed. - * - * The algorithm begins with a value of u that has one more digit of - * input. The normalization step sets u_{m+n}..u_0 = 2^k * u_{m+n-1}..u_0. - * If the multiplication did not increase the number of digits of u, we - * need to add a leading zero here. - */ - if (k == 0 || MP_USED(u) != m + n + 1) - { - if (!s_pad(u, m + n + 1)) - return MP_MEMORY; - u->digits[m + n] = 0; - u->used = m + n + 1; - } - - /* - * Add a leading 0 to v. - * - * The multiplication in step D4 multiplies qhat * 0v_{n-1}..v_0. We need - * to add the leading zero to v here to ensure that the multiplication - * will produce the full n+1 digit result. - */ - if (!s_pad(v, n + 1)) - return MP_MEMORY; - v->digits[n] = 0; - - /* - * Initialize temporary variables q and t. q allocates space for m+1 - * digits to store the quotient digits t allocates space for n+1 digits to - * hold the result of q_j*v - */ - DECLARE_TEMP(2); - REQUIRE(GROW(TEMP(0), m + 1)); - REQUIRE(GROW(TEMP(1), n + 1)); - - /* D2: Initialize j */ - int j = m; - mpz_t r; - - r.digits = MP_DIGITS(u) + j; /* The contents of r are shared with u */ - r.used = n + 1; - r.sign = MP_ZPOS; - r.alloc = MP_ALLOC(u); - ZERO(TEMP(1)->digits, TEMP(1)->alloc); - - /* Calculate the m+1 digits of the quotient result */ - for (; j >= 0; j--) - { - /* D3: Calculate q' */ - /* r->digits is aligned to position j of the number u */ - mp_word pfx, - qhat; - - pfx = r.digits[n]; - pfx <<= MP_DIGIT_BIT / 2; - pfx <<= MP_DIGIT_BIT / 2; - pfx |= r.digits[n - 1]; /* pfx = u_{j+n}{j+n-1} */ - - qhat = pfx / v->digits[n - 1]; - - /* - * Check to see if qhat > b, and decrease qhat if so. Theorem B - * guarantess that qhat is at most 2 larger than the actual value, so - * it is possible that qhat is greater than the maximum value that - * will fit in a digit - */ - if (qhat > MP_DIGIT_MAX) - qhat = MP_DIGIT_MAX; - - /* - * D4,D5,D6: Multiply qhat * v and test for a correct value of q - * - * We proceed a bit different than the way described by Knuth. This - * way is simpler but less efficent. Instead of doing the multiply and - * subtract then checking for underflow, we first do the multiply of - * qhat * v and see if it is larger than the current remainder r. If - * it is larger, we decrease qhat by one and try again. We may need to - * decrease qhat one more time before we get a value that is smaller - * than r. - * - * This way is less efficent than Knuth becuase we do more multiplies, - * but we do not need to worry about underflow this way. - */ - /* t = qhat * v */ - s_dbmul(MP_DIGITS(v), (mp_digit) qhat, TEMP(1)->digits, n + 1); - TEMP(1)->used = n + 1; - CLAMP(TEMP(1)); - - /* Clamp r for the comparison. Comparisons do not like leading zeros. */ - CLAMP(&r); - if (s_ucmp(TEMP(1), &r) > 0) - { /* would the remainder be negative? */ - qhat -= 1; /* try a smaller q */ - s_dbmul(MP_DIGITS(v), (mp_digit) qhat, TEMP(1)->digits, n + 1); - TEMP(1)->used = n + 1; - CLAMP(TEMP(1)); - if (s_ucmp(TEMP(1), &r) > 0) - { /* would the remainder be negative? */ - assert(qhat > 0); - qhat -= 1; /* try a smaller q */ - s_dbmul(MP_DIGITS(v), (mp_digit) qhat, TEMP(1)->digits, n + 1); - TEMP(1)->used = n + 1; - CLAMP(TEMP(1)); - } - assert(s_ucmp(TEMP(1), &r) <= 0 && "The mathematics failed us."); - } - - /* - * Unclamp r. The D algorithm expects r = u_{j+n}..u_j to always be - * n+1 digits long. - */ - r.used = n + 1; - - /* - * D4: Multiply and subtract - * - * Note: The multiply was completed above so we only need to subtract - * here. - */ - s_usub(r.digits, TEMP(1)->digits, r.digits, r.used, TEMP(1)->used); - - /* - * D5: Test remainder - * - * Note: Not needed because we always check that qhat is the correct - * value before performing the subtract. Value cast to mp_digit to - * prevent warning, qhat has been clamped to MP_DIGIT_MAX - */ - TEMP(0)->digits[j] = (mp_digit) qhat; - - /* - * D6: Add back Note: Not needed because we always check that qhat is - * the correct value before performing the subtract. - */ - - /* D7: Loop on j */ - r.digits--; - ZERO(TEMP(1)->digits, TEMP(1)->alloc); - } - - /* Get rid of leading zeros in q */ - TEMP(0)->used = m + 1; - CLAMP(TEMP(0)); - - /* Denormalize the remainder */ - CLAMP(u); /* use u here because the r.digits pointer is - * off-by-one */ - if (k != 0) - s_qdiv(u, k); - - mp_int_copy(u, v); /* ok: 0 <= r < v */ - mp_int_copy(TEMP(0), u); /* ok: q <= u */ - - CLEANUP_TEMP(); - return MP_OK; -} - -static int -s_outlen(mp_int z, mp_size r) -{ - assert(r >= MP_MIN_RADIX && r <= MP_MAX_RADIX); - - mp_result bits = mp_int_count_bits(z); - double raw = (double) bits * s_log2[r]; - - return (int) (raw + 0.999999); -} - -static mp_size -s_inlen(int len, mp_size r) -{ - double raw = (double) len / s_log2[r]; - mp_size bits = (mp_size) (raw + 0.5); - - return (mp_size) ((bits + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT) + 1; -} - -static int -s_ch2val(char c, int r) -{ - int out; - - /* - * In some locales, isalpha() accepts characters outside the range A-Z, - * producing out<0 or out>=36. The "out >= r" check will always catch - * out>=36. Though nothing explicitly catches out<0, our caller reacts - * the same way to every negative return value. - */ - if (isdigit((unsigned char) c)) - out = c - '0'; - else if (r > 10 && isalpha((unsigned char) c)) - out = toupper((unsigned char) c) - 'A' + 10; - else - return -1; - - return (out >= r) ? -1 : out; -} - -static char -s_val2ch(int v, int caps) -{ - assert(v >= 0); - - if (v < 10) - { - return v + '0'; - } - else - { - char out = (v - 10) + 'a'; - - if (caps) - { - return toupper((unsigned char) out); - } - else - { - return out; - } - } -} - -static void -s_2comp(unsigned char *buf, int len) -{ - unsigned short s = 1; - - for (int i = len - 1; i >= 0; --i) - { - unsigned char c = ~buf[i]; - - s = c + s; - c = s & UCHAR_MAX; - s >>= CHAR_BIT; - - buf[i] = c; - } - - /* last carry out is ignored */ -} - -static mp_result -s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad) -{ - int pos = 0, - limit = *limpos; - mp_size uz = MP_USED(z); - mp_digit *dz = MP_DIGITS(z); - - while (uz > 0 && pos < limit) - { - mp_digit d = *dz++; - int i; - - for (i = sizeof(mp_digit); i > 0 && pos < limit; --i) - { - buf[pos++] = (unsigned char) d; - d >>= CHAR_BIT; - - /* Don't write leading zeroes */ - if (d == 0 && uz == 1) - i = 0; /* exit loop without signaling truncation */ - } - - /* Detect truncation (loop exited with pos >= limit) */ - if (i > 0) - break; - - --uz; - } - - if (pad != 0 && (buf[pos - 1] >> (CHAR_BIT - 1))) - { - if (pos < limit) - { - buf[pos++] = 0; - } - else - { - uz = 1; - } - } - - /* Digits are in reverse order, fix that */ - REV(buf, pos); - - /* Return the number of bytes actually written */ - *limpos = pos; - - return (uz == 0) ? MP_OK : MP_TRUNC; -} - -/* Here there be dragons */ |