diff options
Diffstat (limited to '')
| -rw-r--r-- | contrib/pgcrypto/imath.c | 3588 |
1 files changed, 3588 insertions, 0 deletions
diff --git a/contrib/pgcrypto/imath.c b/contrib/pgcrypto/imath.c new file mode 100644 index 0000000..da4cded --- /dev/null +++ b/contrib/pgcrypto/imath.c @@ -0,0 +1,3588 @@ +/*------------------------------------------------------------------------- + * + * 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-2020, 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 = px_alloc(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); + px_free(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 = px_alloc(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 = px_realloc(old, nsize * sizeof(mp_digit)); + + assert(new != NULL); +#endif + + return new; +} + +static void +s_free(void *ptr) +{ + px_free(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 */ |