http://git-wip-us.apache.org/repos/asf/incubator-hawq/blob/801100ed/contrib/pgcrypto/imath.c
----------------------------------------------------------------------
diff --git a/contrib/pgcrypto/imath.c b/contrib/pgcrypto/imath.c
new file mode 100644
index 0000000..3286cd9
--- /dev/null
+++ b/contrib/pgcrypto/imath.c
@@ -0,0 +1,3681 @@
+/* imath version 1.3 */
+/*
+  Name:                imath.c
+  Purpose:     Arbitrary precision integer arithmetic routines.
+  Author:      M. J. Fromberger <http://spinning-yarns.org/michael/sw/>
+  Info:                Id: imath.c 21 2006-04-02 18:58:36Z sting
+
+  Copyright (C) 2002 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.
+ */
+/* contrib/pgcrypto/imath.c */
+
+#include "postgres.h"
+#include "px.h"
+#include "imath.h"
+
+#undef assert
+#define assert(TEST) Assert(TEST)
+#define TRACEABLE_CLAMP 0
+#define TRACEABLE_FREE 0
+
+/* {{{ Constants */
+
+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_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 null argument",
+       NULL
+};
+
+/* }}} */
+
+/* Optional library flags */
+#define MP_CAP_DIGITS  1               /* flag bit to capitalize letter digits 
*/
+
+/* Argument checking macros
+   Use CHECK() where a return value is required; NRCHECK() elsewhere */
+#define CHECK(TEST)   assert(TEST)
+#define NRCHECK(TEST) assert(TEST)
+
+/* {{{ Logarithm table for computing output sizes */
+
+/* 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).
+ */
+static const double s_log2[] = {
+       0.000000000, 0.000000000, 1.000000000, 0.630929754, /* 0  1  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, 0.191958720, 0.190551412, 0.189200360, /* 36 37 38 39 */
+       0.187901825, 0.186652411, 0.185449023, 0.184288833, /* 40 41 42 43 */
+       0.183169251, 0.182087900, 0.181042597, 0.180031327, /* 44 45 46 47 */
+       0.179052232, 0.178103594, 0.177183820, 0.176291434, /* 48 49 50 51 */
+       0.175425064, 0.174583430, 0.173765343, 0.172969690, /* 52 53 54 55 */
+       0.172195434, 0.171441601, 0.170707280, 0.169991616, /* 56 57 58 59 */
+       0.169293808, 0.168613099, 0.167948779, 0.167300179, /* 60 61 62 63 */
+       0.166666667
+};
+
+/* }}} */
+/* {{{ Various macros */
+
+/* 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 */
+#define ROUND_PREC(P) ((mp_size)(2*(((P)+1)/2)))
+
+/* Set array P of S digits to zero */
+#define ZERO(P, S) \
+do{mp_size i__=(S)*sizeof(mp_digit);mp_digit 
*p__=(P);memset(p__,0,i__);}while(0)
+
+/* Copy S digits from array P to array Q */
+#define COPY(P, Q, S) \
+do{mp_size i__=(S)*sizeof(mp_digit);mp_digit *p__=(P),*q__=(Q);\
+memcpy(q__,p__,i__);}while(0)
+
+/* Reverse N elements of type T in array A */
+#define REV(T, A, N) \
+do{T *u_=(A),*v_=u_+(N)-1;while(u_<v_){T xch=*u_;*u_++=*v_;*v_--=xch;}}while(0)
+
+#if TRACEABLE_CLAMP
+#define CLAMP(Z) s_clamp(Z)
+#else
+#define CLAMP(Z) \
+do{mp_int z_=(Z);mp_size uz_=MP_USED(z_);mp_digit *dz_=MP_DIGITS(z_)+uz_-1;\
+while(uz_ > 1 && (*dz_-- == 0)) --uz_;MP_USED(z_)=uz_;}while(0)
+#endif
+
+#undef MIN
+#undef MAX
+#define MIN(A, B) ((B)<(A)?(B):(A))
+#define MAX(A, B) ((B)>(A)?(B):(A))
+#define SWAP(T, A, B) do{T t_=(A);A=(B);B=t_;}while(0)
+
+#define TEMP(K) (temp + (K))
+#define SETUP(E, C) \
+do{if((res = (E)) != MP_OK) goto CLEANUP; ++(C);}while(0)
+
+#define CMPZ(Z) \
+(((Z)->used==1&&(Z)->digits[0]==0)?0:((Z)->sign==MP_NEG)?-1:1)
+
+#define UMUL(X, Y, Z) \
+do{mp_size ua_=MP_USED(X),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_);\
+MP_USED(Z)=o_;CLAMP(Z);}while(0)
+
+#define USQR(X, Z) \
+do{mp_size ua_=MP_USED(X),o_=ua_+ua_;ZERO(MP_DIGITS(Z),o_);\
+(void) s_ksqr(MP_DIGITS(X),MP_DIGITS(Z),ua_);MP_USED(Z)=o_;CLAMP(Z);}while(0)
+
+#define UPPER_HALF(W)                  ((mp_word)((W) >> MP_DIGIT_BIT))
+#define LOWER_HALF(W)                  ((mp_digit)(W))
+#define HIGH_BIT_SET(W)                        ((W) >> (MP_WORD_BIT - 1))
+#define ADD_WILL_OVERFLOW(W, V) ((MP_WORD_MAX - (V)) < (W))
+
+/* }}} */
+
+/* Default number of digits allocated to a new mp_int */
+static mp_size default_precision = 64;
+
+/* Minimum number of digits to invoke recursive multiply */
+static mp_size multiply_threshold = 32;
+
+/* Default library configuration flags */
+static mp_word mp_flags = MP_CAP_DIGITS;
+
+/* 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);
+
+#if TRACEABLE_FREE
+static void s_free(void *ptr);
+#else
+#define s_free(P) px_free(P)
+#endif
+
+/* Insure that z has at least min digits allocated, resizing if
+   necessary.  Returns true if successful, false if out of memory. */
+static int     s_pad(mp_int z, mp_size min);
+
+/* Normalize by removing leading zeroes (except when z = 0) */
+#if TRACEABLE_CLAMP
+static void s_clamp(mp_int z);
+#endif
+
+/* Fill in a "fake" mp_int on the stack with a given value */
+static void s_fake(mp_int z, int 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_vpack(int 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, int v);
+
+/* 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, int 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(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);
+
+#if 0
+/* Dump a representation of the mp_int to standard output */
+void           s_print(char *tag, mp_int z);
+void           s_print_buf(char *tag, mp_digit *buf, mp_size num);
+#endif
+
+/* {{{ get_default_precision() */
+
+mp_size
+mp_get_default_precision(void)
+{
+       return default_precision;
+}
+
+/* }}} */
+
+/* {{{ mp_set_default_precision(s) */
+
+void
+mp_set_default_precision(mp_size s)
+{
+       NRCHECK(s > 0);
+
+       default_precision = (mp_size) ROUND_PREC(s);
+}
+
+/* }}} */
+
+/* {{{ mp_get_multiply_threshold() */
+
+mp_size
+mp_get_multiply_threshold(void)
+{
+       return multiply_threshold;
+}
+
+/* }}} */
+
+/* {{{ mp_set_multiply_threshold(s) */
+
+void
+mp_set_multiply_threshold(mp_size s)
+{
+       multiply_threshold = s;
+}
+
+/* }}} */
+
+/* {{{ mp_int_init(z) */
+
+mp_result
+mp_int_init(mp_int z)
+{
+       return mp_int_init_size(z, default_precision);
+}
+
+/* }}} */
+
+/* {{{ mp_int_alloc() */
+
+mp_int
+mp_int_alloc(void)
+{
+       mp_int          out = px_alloc(sizeof(mpz_t));
+
+       assert(out != NULL);
+       out->digits = NULL;
+       out->used = 0;
+       out->alloc = 0;
+       out->sign = 0;
+
+       return out;
+}
+
+/* }}} */
+
+/* {{{ mp_int_init_size(z, prec) */
+
+mp_result
+mp_int_init_size(mp_int z, mp_size prec)
+{
+       CHECK(z != NULL);
+
+       prec = (mp_size) ROUND_PREC(prec);
+       prec = MAX(prec, default_precision);
+
+       if ((MP_DIGITS(z) = s_alloc(prec)) == NULL)
+               return MP_MEMORY;
+
+       z->digits[0] = 0;
+       MP_USED(z) = 1;
+       MP_ALLOC(z) = prec;
+       MP_SIGN(z) = MP_ZPOS;
+
+       return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_init_copy(z, old) */
+
+mp_result
+mp_int_init_copy(mp_int z, mp_int old)
+{
+       mp_result       res;
+       mp_size         uold,
+                               target;
+
+       CHECK(z != NULL && old != NULL);
+
+       uold = MP_USED(old);
+       target = MAX(uold, default_precision);
+
+       if ((res = mp_int_init_size(z, target)) != MP_OK)
+               return res;
+
+       MP_USED(z) = uold;
+       MP_SIGN(z) = MP_SIGN(old);
+       COPY(MP_DIGITS(old), MP_DIGITS(z), uold);
+
+       return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_init_value(z, value) */
+
+mp_result
+mp_int_init_value(mp_int z, int value)
+{
+       mp_result       res;
+
+       CHECK(z != NULL);
+
+       if ((res = mp_int_init(z)) != MP_OK)
+               return res;
+
+       return mp_int_set_value(z, value);
+}
+
+/* }}} */
+
+/* {{{ mp_int_set_value(z, value) */
+
+mp_result
+mp_int_set_value(mp_int z, int value)
+{
+       mp_size         ndig;
+
+       CHECK(z != NULL);
+
+       /* How many digits to copy */
+       ndig = (mp_size) MP_VALUE_DIGITS(value);
+
+       if (!s_pad(z, ndig))
+               return MP_MEMORY;
+
+       MP_USED(z) = (mp_size) s_vpack(value, MP_DIGITS(z));
+       MP_SIGN(z) = (value < 0) ? MP_NEG : MP_ZPOS;
+
+       return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_clear(z) */
+
+void
+mp_int_clear(mp_int z)
+{
+       if (z == NULL)
+               return;
+
+       if (MP_DIGITS(z) != NULL)
+       {
+               s_free(MP_DIGITS(z));
+               MP_DIGITS(z) = NULL;
+       }
+}
+
+/* }}} */
+
+/* {{{ mp_int_free(z) */
+
+void
+mp_int_free(mp_int z)
+{
+       NRCHECK(z != NULL);
+
+       if (z->digits != NULL)
+               mp_int_clear(z);
+
+       px_free(z);
+}
+
+/* }}} */
+
+/* {{{ mp_int_copy(a, c) */
+
+mp_result
+mp_int_copy(mp_int a, mp_int c)
+{
+       CHECK(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);
+
+               MP_USED(c) = ua;
+               MP_SIGN(c) = MP_SIGN(a);
+       }
+
+       return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_swap(a, c) */
+
+void
+mp_int_swap(mp_int a, mp_int c)
+{
+       if (a != c)
+       {
+               mpz_t           tmp = *a;
+
+               *a = *c;
+               *c = tmp;
+       }
+}
+
+/* }}} */
+
+/* {{{ mp_int_zero(z) */
+
+void
+mp_int_zero(mp_int z)
+{
+       NRCHECK(z != NULL);
+
+       z->digits[0] = 0;
+       MP_USED(z) = 1;
+       MP_SIGN(z) = MP_ZPOS;
+}
+
+/* }}} */
+
+/* {{{ mp_int_abs(a, c) */
+
+mp_result
+mp_int_abs(mp_int a, mp_int c)
+{
+       mp_result       res;
+
+       CHECK(a != NULL && c != NULL);
+
+       if ((res = mp_int_copy(a, c)) != MP_OK)
+               return res;
+
+       MP_SIGN(c) = MP_ZPOS;
+       return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_neg(a, c) */
+
+mp_result
+mp_int_neg(mp_int a, mp_int c)
+{
+       mp_result       res;
+
+       CHECK(a != NULL && c != NULL);
+
+       if ((res = mp_int_copy(a, c)) != MP_OK)
+               return res;
+
+       if (CMPZ(c) != 0)
+               MP_SIGN(c) = 1 - MP_SIGN(a);
+
+       return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_add(a, b, c) */
+
+mp_result
+mp_int_add(mp_int a, mp_int b, mp_int c)
+{
+       mp_size         ua,
+                               ub,
+                               uc,
+                               max;
+
+       CHECK(a != NULL && b != NULL && c != NULL);
+
+       ua = MP_USED(a);
+       ub = MP_USED(b);
+       uc = MP_USED(c);
+       max = MAX(ua, ub);
+
+       if (MP_SIGN(a) == MP_SIGN(b))
+       {
+               /* Same sign -- add magnitudes, preserve sign of addends */
+               mp_digit        carry;
+
+               if (!s_pad(c, max))
+                       return MP_MEMORY;
+
+               carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, 
ub);
+               uc = max;
+
+               if (carry)
+               {
+                       if (!s_pad(c, max + 1))
+                               return MP_MEMORY;
+
+                       c->digits[max] = carry;
+                       ++uc;
+               }
+
+               MP_USED(c) = uc;
+               MP_SIGN(c) = MP_SIGN(a);
+
+       }
+       else
+       {
+               /* Different signs -- subtract magnitudes, preserve sign of 
greater */
+               mp_int          x,
+                                       y;
+               int                     cmp = s_ucmp(a, b); /* magnitude 
comparison, sign ignored */
+
+               /* Set x to max(a, b), y to min(a, b) to simplify later code */
+               if (cmp >= 0)
+               {
+                       x = a;
+                       y = b;
+               }
+               else
+               {
+                       x = b;
+                       y = a;
+               }
+
+               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));
+               MP_USED(c) = MP_USED(x);
+               CLAMP(c);
+
+               /* Give result the sign of the larger */
+               MP_SIGN(c) = MP_SIGN(x);
+       }
+
+       return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_add_value(a, value, c) */
+
+mp_result
+mp_int_add_value(mp_int a, int 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_int_sub(a, b, c) */
+
+mp_result
+mp_int_sub(mp_int a, mp_int b, mp_int c)
+{
+       mp_size         ua,
+                               ub,
+                               uc,
+                               max;
+
+       CHECK(a != NULL && b != NULL && c != NULL);
+
+       ua = MP_USED(a);
+       ub = MP_USED(b);
+       uc = MP_USED(c);
+       max = MAX(ua, ub);
+
+       if (MP_SIGN(a) != MP_SIGN(b))
+       {
+               /* Different signs -- add magnitudes and keep sign of a */
+               mp_digit        carry;
+
+               if (!s_pad(c, max))
+                       return MP_MEMORY;
+
+               carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, 
ub);
+               uc = max;
+
+               if (carry)
+               {
+                       if (!s_pad(c, max + 1))
+                               return MP_MEMORY;
+
+                       c->digits[max] = carry;
+                       ++uc;
+               }
+
+               MP_USED(c) = uc;
+               MP_SIGN(c) = MP_SIGN(a);
+
+       }
+       else
+       {
+               /* Same signs -- subtract magnitudes */
+               mp_int          x,
+                                       y;
+               mp_sign         osign;
+               int                     cmp = s_ucmp(a, b);
+
+               if (!s_pad(c, max))
+                       return MP_MEMORY;
+
+               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));
+               MP_USED(c) = MP_USED(x);
+               CLAMP(c);
+
+               MP_SIGN(c) = osign;
+       }
+
+       return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_sub_value(a, value, c) */
+
+mp_result
+mp_int_sub_value(mp_int a, int 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_int_mul(a, b, c) */
+
+mp_result
+mp_int_mul(mp_int a, mp_int b, mp_int c)
+{
+       mp_digit   *out;
+       mp_size         osize,
+                               ua,
+                               ub,
+                               p = 0;
+       mp_sign         osign;
+
+       CHECK(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 */
+       osign = (MP_SIGN(a) == MP_SIGN(b)) ? MP_ZPOS : MP_NEG;
+
+       /*
+        * If the output is not equal to any of the inputs, we'll write the
+        * results there directly; otherwise, allocate a temporary space.
+        */
+       ua = MP_USED(a);
+       ub = MP_USED(b);
+       osize = ua + ub;
+
+       if (c == a || c == b)
+       {
+               p = 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);
+
+       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))
+       {
+               s_free(MP_DIGITS(c));
+               MP_DIGITS(c) = out;
+               MP_ALLOC(c) = p;
+       }
+
+       MP_USED(c) = osize;                     /* might not be true, but we'll 
fix it ... */
+       CLAMP(c);                                       /* ... right here */
+       MP_SIGN(c) = osign;
+
+       return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_mul_value(a, value, c) */
+
+mp_result
+mp_int_mul_value(mp_int a, int 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_int_mul_pow2(a, p2, c) */
+
+mp_result
+mp_int_mul_pow2(mp_int a, int p2, mp_int c)
+{
+       mp_result       res;
+
+       CHECK(a != NULL && c != NULL && p2 >= 0);
+
+       if ((res = mp_int_copy(a, c)) != MP_OK)
+               return res;
+
+       if (s_qmul(c, (mp_size) p2))
+               return MP_OK;
+       else
+               return MP_MEMORY;
+}
+
+/* }}} */
+
+/* {{{ mp_int_sqr(a, c) */
+
+mp_result
+mp_int_sqr(mp_int a, mp_int c)
+{
+       mp_digit   *out;
+       mp_size         osize,
+                               p = 0;
+
+       CHECK(a != NULL && c != NULL);
+
+       /* Get a temporary buffer big enough to hold the result */
+       osize = (mp_size) 2 *MP_USED(a);
+
+       if (a == c)
+       {
+               p = 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))
+       {
+               s_free(MP_DIGITS(c));
+               MP_DIGITS(c) = out;
+               MP_ALLOC(c) = p;
+       }
+
+       MP_USED(c) = osize;                     /* might not be true, but we'll 
fix it ... */
+       CLAMP(c);                                       /* ... right here */
+       MP_SIGN(c) = MP_ZPOS;
+
+       return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_div(a, b, q, r) */
+
+mp_result
+mp_int_div(mp_int a, mp_int b, mp_int q, mp_int r)
+{
+       int                     cmp,
+                               last = 0,
+                               lg;
+       mp_result       res = MP_OK;
+       mpz_t           temp[2];
+       mp_int          qout,
+                               rout;
+       mp_sign         sa = MP_SIGN(a),
+                               sb = MP_SIGN(b);
+
+       CHECK(a != NULL && b != NULL && q != r);
+
+       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)
+                               MP_SIGN(q) = 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.
+        */
+       if ((lg = s_isp2(b)) < 0)
+       {
+               if (q && b != q && (res = mp_int_copy(a, q)) == MP_OK)
+               {
+                       qout = q;
+               }
+               else
+               {
+                       qout = TEMP(last);
+                       SETUP(mp_int_init_copy(TEMP(last), a), last);
+               }
+
+               if (r && a != r && (res = mp_int_copy(b, r)) == MP_OK)
+               {
+                       rout = r;
+               }
+               else
+               {
+                       rout = TEMP(last);
+                       SETUP(mp_int_init_copy(TEMP(last), b), last);
+               }
+
+               if ((res = s_udiv(qout, rout)) != MP_OK)
+                       goto CLEANUP;
+       }
+       else
+       {
+               if (q && (res = mp_int_copy(a, q)) != MP_OK)
+                       goto CLEANUP;
+               if (r && (res = mp_int_copy(a, r)) != MP_OK)
+                       goto CLEANUP;
+
+               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)
+       {
+               MP_SIGN(rout) = sa;
+               if (CMPZ(rout) == 0)
+                       MP_SIGN(rout) = MP_ZPOS;
+       }
+       if (qout)
+       {
+               MP_SIGN(qout) = (sa == sb) ? MP_ZPOS : MP_NEG;
+               if (CMPZ(qout) == 0)
+                       MP_SIGN(qout) = MP_ZPOS;
+       }
+
+       if (q && (res = mp_int_copy(qout, q)) != MP_OK)
+               goto CLEANUP;
+       if (r && (res = mp_int_copy(rout, r)) != MP_OK)
+               goto CLEANUP;
+
+CLEANUP:
+       while (--last >= 0)
+               mp_int_clear(TEMP(last));
+
+       return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_mod(a, m, c) */
+
+mp_result
+mp_int_mod(mp_int a, mp_int m, mp_int c)
+{
+       mp_result       res;
+       mpz_t           tmp;
+       mp_int          out;
+
+       if (m == c)
+       {
+               if ((res = mp_int_init(&tmp)) != MP_OK)
+                       return res;
+
+               out = &tmp;
+       }
+       else
+       {
+               out = c;
+       }
+
+       if ((res = mp_int_div(a, m, NULL, out)) != MP_OK)
+               goto CLEANUP;
+
+       if (CMPZ(out) < 0)
+               res = mp_int_add(out, m, c);
+       else
+               res = mp_int_copy(out, c);
+
+CLEANUP:
+       if (out != c)
+               mp_int_clear(&tmp);
+
+       return res;
+}
+
+/* }}} */
+
+
+/* {{{ mp_int_div_value(a, value, q, r) */
+
+mp_result
+mp_int_div_value(mp_int a, int value, mp_int q, int *r)
+{
+       mpz_t           vtmp,
+                               rtmp;
+       mp_digit        vbuf[MP_VALUE_DIGITS(value)];
+       mp_result       res;
+
+       if ((res = mp_int_init(&rtmp)) != MP_OK)
+               return res;
+       s_fake(&vtmp, value, vbuf);
+
+       if ((res = mp_int_div(a, &vtmp, q, &rtmp)) != MP_OK)
+               goto CLEANUP;
+
+       if (r)
+               (void) mp_int_to_int(&rtmp, r); /* can't fail */
+
+CLEANUP:
+       mp_int_clear(&rtmp);
+       return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_div_pow2(a, p2, q, r) */
+
+mp_result
+mp_int_div_pow2(mp_int a, int p2, mp_int q, mp_int r)
+{
+       mp_result       res = MP_OK;
+
+       CHECK(a != NULL && p2 >= 0 && q != r);
+
+       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_int_expt(a, b, c) */
+
+mp_result
+mp_int_expt(mp_int a, int b, mp_int c)
+{
+       mpz_t           t;
+       mp_result       res;
+       unsigned int v = abs(b);
+
+       CHECK(b >= 0 && c != NULL);
+
+       if ((res = mp_int_init_copy(&t, a)) != MP_OK)
+               return res;
+
+       (void) mp_int_set_value(c, 1);
+       while (v != 0)
+       {
+               if (v & 1)
+               {
+                       if ((res = mp_int_mul(c, &t, c)) != MP_OK)
+                               goto CLEANUP;
+               }
+
+               v >>= 1;
+               if (v == 0)
+                       break;
+
+               if ((res = mp_int_sqr(&t, &t)) != MP_OK)
+                       goto CLEANUP;
+       }
+
+CLEANUP:
+       mp_int_clear(&t);
+       return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_expt_value(a, b, c) */
+
+mp_result
+mp_int_expt_value(int a, int b, mp_int c)
+{
+       mpz_t           t;
+       mp_result       res;
+       unsigned int v = abs(b);
+
+       CHECK(b >= 0 && c != NULL);
+
+       if ((res = mp_int_init_value(&t, a)) != MP_OK)
+               return res;
+
+       (void) mp_int_set_value(c, 1);
+       while (v != 0)
+       {
+               if (v & 1)
+               {
+                       if ((res = mp_int_mul(c, &t, c)) != MP_OK)
+                               goto CLEANUP;
+               }
+
+               v >>= 1;
+               if (v == 0)
+                       break;
+
+               if ((res = mp_int_sqr(&t, &t)) != MP_OK)
+                       goto CLEANUP;
+       }
+
+CLEANUP:
+       mp_int_clear(&t);
+       return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_compare(a, b) */
+
+int
+mp_int_compare(mp_int a, mp_int b)
+{
+       mp_sign         sa;
+
+       CHECK(a != NULL && b != NULL);
+
+       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;
+       }
+}
+
+/* }}} */
+
+/* {{{ mp_int_compare_unsigned(a, b) */
+
+int
+mp_int_compare_unsigned(mp_int a, mp_int b)
+{
+       NRCHECK(a != NULL && b != NULL);
+
+       return s_ucmp(a, b);
+}
+
+/* }}} */
+
+/* {{{ mp_int_compare_zero(z) */
+
+int
+mp_int_compare_zero(mp_int z)
+{
+       NRCHECK(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;
+}
+
+/* }}} */
+
+/* {{{ mp_int_compare_value(z, value) */
+
+int
+mp_int_compare_value(mp_int z, int value)
+{
+       mp_sign         vsign = (value < 0) ? MP_NEG : MP_ZPOS;
+       int                     cmp;
+
+       CHECK(z != NULL);
+
+       if (vsign == MP_SIGN(z))
+       {
+               cmp = s_vcmp(z, value);
+
+               if (vsign == MP_ZPOS)
+                       return cmp;
+               else
+                       return -cmp;
+       }
+       else
+       {
+               if (value < 0)
+                       return 1;
+               else
+                       return -1;
+       }
+}
+
+/* }}} */
+
+/* {{{ mp_int_exptmod(a, b, m, c) */
+
+mp_result
+mp_int_exptmod(mp_int a, mp_int b, mp_int m, mp_int c)
+{
+       mp_result       res;
+       mp_size         um;
+       mpz_t           temp[3];
+       mp_int          s;
+       int                     last = 0;
+
+       CHECK(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;
+
+       um = MP_USED(m);
+       SETUP(mp_int_init_size(TEMP(0), 2 * um), last);
+       SETUP(mp_int_init_size(TEMP(1), 2 * um), last);
+
+       if (c == b || c == m)
+       {
+               SETUP(mp_int_init_size(TEMP(2), 2 * um), last);
+               s = TEMP(2);
+       }
+       else
+       {
+               s = c;
+       }
+
+       if ((res = mp_int_mod(a, m, TEMP(0))) != MP_OK)
+               goto CLEANUP;
+
+       if ((res = s_brmu(TEMP(1), m)) != MP_OK)
+               goto CLEANUP;
+
+       if ((res = s_embar(TEMP(0), b, m, TEMP(1), s)) != MP_OK)
+               goto CLEANUP;
+
+       res = mp_int_copy(s, c);
+
+CLEANUP:
+       while (--last >= 0)
+               mp_int_clear(TEMP(last));
+
+       return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_exptmod_evalue(a, value, m, c) */
+
+mp_result
+mp_int_exptmod_evalue(mp_int a, int 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_int_exptmod_bvalue(v, b, m, c) */
+
+mp_result
+mp_int_exptmod_bvalue(int 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_int_exptmod_known(a, b, m, mu, c) */
+
+mp_result
+mp_int_exptmod_known(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c)
+{
+       mp_result       res;
+       mp_size         um;
+       mpz_t           temp[2];
+       mp_int          s;
+       int                     last = 0;
+
+       CHECK(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;
+
+       um = MP_USED(m);
+       SETUP(mp_int_init_size(TEMP(0), 2 * um), last);
+
+       if (c == b || c == m)
+       {
+               SETUP(mp_int_init_size(TEMP(1), 2 * um), last);
+               s = TEMP(1);
+       }
+       else
+       {
+               s = c;
+       }
+
+       if ((res = mp_int_mod(a, m, TEMP(0))) != MP_OK)
+               goto CLEANUP;
+
+       if ((res = s_embar(TEMP(0), b, m, mu, s)) != MP_OK)
+               goto CLEANUP;
+
+       res = mp_int_copy(s, c);
+
+CLEANUP:
+       while (--last >= 0)
+               mp_int_clear(TEMP(last));
+
+       return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_redux_const(m, c) */
+
+mp_result
+mp_int_redux_const(mp_int m, mp_int c)
+{
+       CHECK(m != NULL && c != NULL && m != c);
+
+       return s_brmu(c, m);
+}
+
+/* }}} */
+
+/* {{{ mp_int_invmod(a, m, c) */
+
+mp_result
+mp_int_invmod(mp_int a, mp_int m, mp_int c)
+{
+       mp_result       res;
+       mp_sign         sa;
+       int                     last = 0;
+       mpz_t           temp[2];
+
+       CHECK(a != NULL && m != NULL && c != NULL);
+
+       if (CMPZ(a) == 0 || CMPZ(m) <= 0)
+               return MP_RANGE;
+
+       sa = MP_SIGN(a);                        /* need this for the result 
later */
+
+       for (last = 0; last < 2; ++last)
+               if ((res = mp_int_init(TEMP(last))) != MP_OK)
+                       goto CLEANUP;
+
+       if ((res = mp_int_egcd(a, m, TEMP(0), TEMP(1), NULL)) != MP_OK)
+               goto CLEANUP;
+
+       if (mp_int_compare_value(TEMP(0), 1) != 0)
+       {
+               res = MP_UNDEF;
+               goto CLEANUP;
+       }
+
+       /* It is first necessary to constrain the value to the proper range */
+       if ((res = mp_int_mod(TEMP(1), m, TEMP(1))) != MP_OK)
+               goto CLEANUP;
+
+       /*
+        * 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 (sa == MP_NEG)
+               res = mp_int_sub(m, TEMP(1), c);
+       else
+               res = mp_int_copy(TEMP(1), c);
+
+CLEANUP:
+       while (--last >= 0)
+               mp_int_clear(TEMP(last));
+
+       return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_gcd(a, b, c) */
+
+/* Binary GCD algorithm due to Josef Stein, 1961 */
+mp_result
+mp_int_gcd(mp_int a, mp_int b, mp_int c)
+{
+       int                     ca,
+                               cb,
+                               k = 0;
+       mpz_t           u,
+                               v,
+                               t;
+       mp_result       res;
+
+       CHECK(a != NULL && b != NULL && c != NULL);
+
+       ca = CMPZ(a);
+       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);
+
+       if ((res = mp_int_init(&t)) != MP_OK)
+               return res;
+       if ((res = mp_int_init_copy(&u, a)) != MP_OK)
+               goto U;
+       if ((res = mp_int_init_copy(&v, b)) != MP_OK)
+               goto V;
+
+       MP_SIGN(&u) = MP_ZPOS;
+       MP_SIGN(&v) = MP_ZPOS;
+
+       {                                                       /* Divide out 
common factors of 2 from u and v */
+               int                     div2_u = s_dp2k(&u),
+                                       div2_v = s_dp2k(&v);
+
+               k = MIN(div2_u, div2_v);
+               s_qdiv(&u, (mp_size) k);
+               s_qdiv(&v, (mp_size) k);
+       }
+
+       if (mp_int_is_odd(&u))
+       {
+               if ((res = mp_int_neg(&v, &t)) != MP_OK)
+                       goto CLEANUP;
+       }
+       else
+       {
+               if ((res = mp_int_copy(&u, &t)) != MP_OK)
+                       goto CLEANUP;
+       }
+
+       for (;;)
+       {
+               s_qdiv(&t, s_dp2k(&t));
+
+               if (CMPZ(&t) > 0)
+               {
+                       if ((res = mp_int_copy(&t, &u)) != MP_OK)
+                               goto CLEANUP;
+               }
+               else
+               {
+                       if ((res = mp_int_neg(&t, &v)) != MP_OK)
+                               goto CLEANUP;
+               }
+
+               if ((res = mp_int_sub(&u, &v, &t)) != MP_OK)
+                       goto CLEANUP;
+
+               if (CMPZ(&t) == 0)
+                       break;
+       }
+
+       if ((res = mp_int_abs(&u, c)) != MP_OK)
+               goto CLEANUP;
+       if (!s_qmul(c, (mp_size) k))
+               res = MP_MEMORY;
+
+CLEANUP:
+       mp_int_clear(&v);
+V: mp_int_clear(&u);
+U: mp_int_clear(&t);
+
+       return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_egcd(a, b, c, x, y) */
+
+/* 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)
+{
+       int                     k,
+                               last = 0,
+                               ca,
+                               cb;
+       mpz_t           temp[8];
+       mp_result       res;
+
+       CHECK(a != NULL && b != NULL && c != NULL &&
+                 (x != NULL || y != NULL));
+
+       ca = CMPZ(a);
+       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
+        */
+       for (last = 0; last < 4; ++last)
+       {
+               if ((res = mp_int_init(TEMP(last))) != MP_OK)
+                       goto CLEANUP;
+       }
+       TEMP(0)->digits[0] = 1;
+       TEMP(3)->digits[0] = 1;
+
+       SETUP(mp_int_init_copy(TEMP(4), a), last);
+       SETUP(mp_int_init_copy(TEMP(5), b), last);
+
+       /* We will work with absolute values here */
+       MP_SIGN(TEMP(4)) = MP_ZPOS;
+       MP_SIGN(TEMP(5)) = MP_ZPOS;
+
+       {                                                       /* 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);
+       }
+
+       SETUP(mp_int_init_copy(TEMP(6), TEMP(4)), last);
+       SETUP(mp_int_init_copy(TEMP(7), TEMP(5)), last);
+
+       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)))
+                       {
+                               if ((res = mp_int_add(TEMP(0), TEMP(7), 
TEMP(0))) != MP_OK)
+                                       goto CLEANUP;
+                               if ((res = mp_int_sub(TEMP(1), TEMP(6), 
TEMP(1))) != MP_OK)
+                                       goto CLEANUP;
+                       }
+
+                       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)))
+                       {
+                               if ((res = mp_int_add(TEMP(2), TEMP(7), 
TEMP(2))) != MP_OK)
+                                       goto CLEANUP;
+                               if ((res = mp_int_sub(TEMP(3), TEMP(6), 
TEMP(3))) != MP_OK)
+                                       goto CLEANUP;
+                       }
+
+                       s_qdiv(TEMP(2), 1);
+                       s_qdiv(TEMP(3), 1);
+               }
+
+               if (mp_int_compare(TEMP(4), TEMP(5)) >= 0)
+               {
+                       if ((res = mp_int_sub(TEMP(4), TEMP(5), TEMP(4))) != 
MP_OK)
+                               goto CLEANUP;
+                       if ((res = mp_int_sub(TEMP(0), TEMP(2), TEMP(0))) != 
MP_OK)
+                               goto CLEANUP;
+                       if ((res = mp_int_sub(TEMP(1), TEMP(3), TEMP(1))) != 
MP_OK)
+                               goto CLEANUP;
+               }
+               else
+               {
+                       if ((res = mp_int_sub(TEMP(5), TEMP(4), TEMP(5))) != 
MP_OK)
+                               goto CLEANUP;
+                       if ((res = mp_int_sub(TEMP(2), TEMP(0), TEMP(2))) != 
MP_OK)
+                               goto CLEANUP;
+                       if ((res = mp_int_sub(TEMP(3), TEMP(1), TEMP(3))) != 
MP_OK)
+                               goto CLEANUP;
+               }
+
+               if (CMPZ(TEMP(4)) == 0)
+               {
+                       if (x && (res = mp_int_copy(TEMP(2), x)) != MP_OK)
+                               goto CLEANUP;
+                       if (y && (res = mp_int_copy(TEMP(3), y)) != MP_OK)
+                               goto CLEANUP;
+                       if (c)
+                       {
+                               if (!s_qmul(TEMP(5), k))
+                               {
+                                       res = MP_MEMORY;
+                                       goto CLEANUP;
+                               }
+
+                               res = mp_int_copy(TEMP(5), c);
+                       }
+
+                       break;
+               }
+       }
+
+CLEANUP:
+       while (--last >= 0)
+               mp_int_clear(TEMP(last));
+
+       return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_divisible_value(a, v) */
+
+int
+mp_int_divisible_value(mp_int a, int v)
+{
+       int                     rem = 0;
+
+       if (mp_int_div_value(a, v, NULL, &rem) != MP_OK)
+               return 0;
+
+       return rem == 0;
+}
+
+/* }}} */
+
+/* {{{ mp_int_is_pow2(z) */
+
+int
+mp_int_is_pow2(mp_int z)
+{
+       CHECK(z != NULL);
+
+       return s_isp2(z);
+}
+
+/* }}} */
+
+/* {{{ mp_int_sqrt(a, c) */
+
+mp_result
+mp_int_sqrt(mp_int a, mp_int c)
+{
+       mp_result       res = MP_OK;
+       mpz_t           temp[2];
+       int                     last = 0;
+
+       CHECK(a != NULL && c != NULL);
+
+       /* The square root of a negative value does not exist in the integers. 
*/
+       if (MP_SIGN(a) == MP_NEG)
+               return MP_UNDEF;
+
+       SETUP(mp_int_init_copy(TEMP(last), a), last);
+       SETUP(mp_int_init(TEMP(last)), last);
+
+       for (;;)
+       {
+               if ((res = mp_int_sqr(TEMP(0), TEMP(1))) != MP_OK)
+                       goto CLEANUP;
+
+               if (mp_int_compare_unsigned(a, TEMP(1)) == 0)
+                       break;
+
+               if ((res = mp_int_copy(a, TEMP(1))) != MP_OK)
+                       goto CLEANUP;
+               if ((res = mp_int_div(TEMP(1), TEMP(0), TEMP(1), NULL)) != 
MP_OK)
+                       goto CLEANUP;
+               if ((res = mp_int_add(TEMP(0), TEMP(1), TEMP(1))) != MP_OK)
+                       goto CLEANUP;
+               if ((res = mp_int_div_pow2(TEMP(1), 1, TEMP(1), NULL)) != MP_OK)
+                       goto CLEANUP;
+
+               if (mp_int_compare_unsigned(TEMP(0), TEMP(1)) == 0)
+                       break;
+               if ((res = mp_int_sub_value(TEMP(0), 1, TEMP(0))) != MP_OK)
+                       goto CLEANUP;
+               if (mp_int_compare_unsigned(TEMP(0), TEMP(1)) == 0)
+                       break;
+
+               if ((res = mp_int_copy(TEMP(1), TEMP(0))) != MP_OK)
+                       goto CLEANUP;
+       }
+
+       res = mp_int_copy(TEMP(0), c);
+
+CLEANUP:
+       while (--last >= 0)
+               mp_int_clear(TEMP(last));
+
+       return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_to_int(z, out) */
+
+mp_result
+mp_int_to_int(mp_int z, int *out)
+{
+       unsigned int uv = 0;
+       mp_size         uz;
+       mp_digit   *dz;
+       mp_sign         sz;
+
+       CHECK(z != NULL);
+
+       /* Make sure the value is representable as an int */
+       sz = MP_SIGN(z);
+       if ((sz == MP_ZPOS && mp_int_compare_value(z, INT_MAX) > 0) ||
+               mp_int_compare_value(z, INT_MIN) < 0)
+               return MP_RANGE;
+
+       uz = MP_USED(z);
+       dz = MP_DIGITS(z) + uz - 1;
+
+       while (uz > 0)
+       {
+               uv <<= MP_DIGIT_BIT / 2;
+               uv = (uv << (MP_DIGIT_BIT / 2)) | *dz--;
+               --uz;
+       }
+
+       if (out)
+               *out = (sz == MP_NEG) ? -(int) uv : (int) uv;
+
+       return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_to_string(z, radix, str, limit) */
+
+mp_result
+mp_int_to_string(mp_int z, mp_size radix,
+                                char *str, int limit)
+{
+       mp_result       res;
+       int                     cmp = 0;
+
+       CHECK(z != NULL && str != NULL && limit >= 2);
+
+       if (radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
+               return MP_RANGE;
+
+       if (CMPZ(z) == 0)
+       {
+               *str++ = s_val2ch(0, mp_flags & MP_CAP_DIGITS);
+       }
+       else
+       {
+               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, mp_flags & MP_CAP_DIGITS);
+               }
+               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_int_string_len(z, radix) */
+
+mp_result
+mp_int_string_len(mp_int z, mp_size radix)
+{
+       int                     len;
+
+       CHECK(z != NULL);
+
+       if (radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
+               return MP_RANGE;
+
+       len = s_outlen(z, radix) + 1;           /* for terminator */
+
+       /* Allow for sign marker on negatives */
+       if (MP_SIGN(z) == MP_NEG)
+               len += 1;
+
+       return len;
+}
+
+/* }}} */
+
+/* {{{ mp_int_read_string(z, radix, *str) */
+
+/* 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_int_read_cstring(z, radix, *str, **end) */
+
+mp_result
+mp_int_read_cstring(mp_int z, mp_size radix, const char *str, char **end)
+{
+       int                     ch;
+
+       CHECK(z != NULL && str != NULL);
+
+       if (radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
+               return MP_RANGE;
+
+       /* Skip leading whitespace */
+       while (isspace((unsigned char) *str))
+               ++str;
+
+       /* Handle leading sign tag (+/-, positive default) */
+       switch (*str)
+       {
+               case '-':
+                       MP_SIGN(z) = MP_NEG;
+                       ++str;
+                       break;
+               case '+':
+                       ++str;                          /* fallthrough */
+               default:
+                       MP_SIGN(z) = MP_ZPOS;
+                       break;
+       }
+
+       /* Skip leading zeroes */
+       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;
+
+       MP_USED(z) = 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)
+               MP_SIGN(z) = MP_ZPOS;
+
+       if (end != NULL)
+               *end = (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_int_count_bits(z) */
+
+mp_result
+mp_int_count_bits(mp_int z)
+{
+       mp_size         nbits = 0,
+                               uz;
+       mp_digit        d;
+
+       CHECK(z != NULL);
+
+       uz = MP_USED(z);
+       if (uz == 1 && z->digits[0] == 0)
+               return 1;
+
+       --uz;
+       nbits = uz * MP_DIGIT_BIT;
+       d = z->digits[uz];
+
+       while (d != 0)
+       {
+               d >>= 1;
+               ++nbits;
+       }
+
+       return nbits;
+}
+
+/* }}} */
+
+/* {{{ mp_int_to_binary(z, buf, limit) */
+
+mp_result
+mp_int_to_binary(mp_int z, unsigned char *buf, int limit)
+{
+       static const int PAD_FOR_2C = 1;
+
+       mp_result       res;
+       int                     limpos = limit;
+
+       CHECK(z != NULL && buf != NULL);
+
+       res = s_tobin(z, buf, &limpos, PAD_FOR_2C);
+
+       if (MP_SIGN(z) == MP_NEG)
+               s_2comp(buf, limpos);
+
+       return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_read_binary(z, buf, len) */
+
+mp_result
+mp_int_read_binary(mp_int z, unsigned char *buf, int len)
+{
+       mp_size         need,
+                               i;
+       unsigned char *tmp;
+       mp_digit   *dz;
+
+       CHECK(z != NULL && buf != NULL && len > 0);
+
+       /* Figure out how many digits are needed to represent this value */
+       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))
+       {
+               MP_SIGN(z) = MP_NEG;
+               s_2comp(buf, len);
+       }
+
+       dz = MP_DIGITS(z);
+       for (tmp = buf, 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_int_binary_len(z) */
+
+mp_result
+mp_int_binary_len(mp_int z)
+{
+       mp_result       res = mp_int_count_bits(z);
+       int                     bytes = mp_int_unsigned_len(z);
+
+       if (res <= 0)
+               return res;
+
+       bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT;
+
+       /*
+        * 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_int_to_unsigned(z, buf, limit) */
+
+mp_result
+mp_int_to_unsigned(mp_int z, unsigned char *buf, int limit)
+{
+       static const int NO_PADDING = 0;
+
+       CHECK(z != NULL && buf != NULL);
+
+       return s_tobin(z, buf, &limit, NO_PADDING);
+}
+
+/* }}} */
+
+/* {{{ mp_int_read_unsigned(z, buf, len) */
+
+mp_result
+mp_int_read_unsigned(mp_int z, unsigned char *buf, int len)
+{
+       mp_size         need,
+                               i;
+       unsigned char *tmp;
+       mp_digit   *dz;
+
+       CHECK(z != NULL && buf != NULL && len > 0);
+
+       /* Figure out how many digits are needed to represent this value */
+       need = ((len * CHAR_BIT) + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT;
+       if (!s_pad(z, need))
+               return MP_MEMORY;
+
+       mp_int_zero(z);
+
+       dz = MP_DIGITS(z);
+       for (tmp = buf, i = len; i > 0; --i, ++tmp)
+       {
+               (void) s_qmul(z, CHAR_BIT);
+               *dz |= *tmp;
+       }
+
+       return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_unsigned_len(z) */
+
+mp_result
+mp_int_unsigned_len(mp_int z)
+{
+       mp_result       res = mp_int_count_bits(z);
+       int                     bytes;
+
+       if (res <= 0)
+               return res;
+
+       bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT;
+
+       return bytes;
+}
+
+/* }}} */
+
+/* {{{ mp_error_string(res) */
+
+const char *
+mp_error_string(mp_result res)
+{
+       int                     ix;
+
+       if (res > 0)
+               return s_unknown_err;
+
+       res = -res;
+       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.                 
  */
+
+/* {{{ s_alloc(num) */
+
+static mp_digit *
+s_alloc(mp_size num)
+{
+       mp_digit   *out = px_alloc(num * sizeof(mp_digit));
+
+       assert(out != NULL);            /* for debugging */
+
+       return out;
+}
+
+/* }}} */
+
+/* {{{ s_realloc(old, num) */
+
+static mp_digit *
+s_realloc(mp_digit *old, mp_size num)
+{
+       mp_digit   *new = px_realloc(old, num * sizeof(mp_digit));
+
+       assert(new != NULL);            /* for debugging */
+
+       return new;
+}
+
+/* }}} */
+
+/* {{{ s_free(ptr) */
+
+#if TRACEABLE_FREE
+static void
+s_free(void *ptr)
+{
+       px_free(ptr);
+}
+#endif
+
+/* }}} */
+
+/* {{{ s_pad(z, min) */
+
+static int
+s_pad(mp_int z, mp_size min)
+{
+       if (MP_ALLOC(z) < min)
+       {
+               mp_size         nsize = ROUND_PREC(min);
+               mp_digit   *tmp = s_realloc(MP_DIGITS(z), nsize);
+
+               if (tmp == NULL)
+                       return 0;
+
+               MP_DIGITS(z) = tmp;
+               MP_ALLOC(z) = nsize;
+       }
+
+       return 1;
+}
+
+/* }}} */
+
+/* {{{ s_clamp(z) */
+
+#if TRACEABLE_CLAMP
+static void
+s_clamp(mp_int z)
+{
+       mp_size         uz = MP_USED(z);
+       mp_digit   *zd = MP_DIGITS(z) + uz - 1;
+
+       while (uz > 1 && (*zd-- == 0))
+               --uz;
+
+       MP_USED(z) = uz;
+}
+#endif
+
+/* }}} */
+
+/* {{{ s_fake(z, value, vbuf) */
+
+static void
+s_fake(mp_int z, int value, mp_digit vbuf[])
+{
+       mp_size         uv = (mp_size) s_vpack(value, vbuf);
+
+       z->used = uv;
+       z->alloc = MP_VALUE_DIGITS(value);
+       z->sign = (value < 0) ? MP_NEG : MP_ZPOS;
+       z->digits = vbuf;
+}
+
+/* }}} */
+
+/* {{{ s_cdig(da, db, len) */
+
+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;
+}
+
+/* }}} */
+
+/* {{{ s_vpack(v, t[]) */
+
+static int
+s_vpack(int v, mp_digit t[])
+{
+       unsigned int uv = (unsigned int) ((v < 0) ? -v : v);
+       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;
+}
+
+/* }}} */
+
+/* {{{ s_ucmp(a, b) */
+
+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);
+}
+
+/* }}} */
+
+/* {{{ s_vcmp(a, v) */
+
+static int
+s_vcmp(mp_int a, int v)
+{
+       mp_digit        vdig[MP_VALUE_DIGITS(v)];
+       int                     ndig = 0;
+       mp_size         ua = MP_USED(a);
+
+       ndig = s_vpack(v, vdig);
+
+       if (ua > ndig)
+               return 1;
+       else if (ua < ndig)
+               return -1;
+       else
+               return s_cdig(MP_DIGITS(a), vdig, ndig);
+}
+
+/* }}} */
+
+/* {{{ s_uadd(da, db, dc, size_a, size_b) */
+
+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;
+}
+
+/* }}} */
+
+/* {{{ s_usub(da, db, dc, size_a, size_b) */
+
+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);
+}
+
+/* }}} */
+
+/* {{{ s_kmul(da, db, dc, size_a, size_b) */
+
+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, bot_size + 1);
+               ZERO(t2, bot_size + 1);
+               (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);
+               (void) s_uadd(t3, dc + bot_size, dc + bot_size,
+                                         buf_size + 1, buf_size + 1);
+
+               (void) s_uadd(t2, dc + 2 * bot_size, dc + 2 * bot_size,
+                                         buf_size, buf_size);
+
+               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;
+}
+
+/* }}} */
+
+/* {{{ s_umul(da, db, dc, size_a, size_b) */
+
+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;
+       }
+}
+
+/* }}} */
+
+/* {{{ s_ksqr(da, dc, size_a) */
+
+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;
+               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);
+               (void) s_uadd(t3, dc + bot_size, dc + bot_size,
+                                         buf_size + 1, buf_size + 1);
+
+               (void) s_uadd(t2, dc + 2 * bot_size, dc + 2 * bot_size,
+                                         buf_size, buf_size);
+
+               px_free(t1);                    /* note that t2 and t2 are 
internal pointers
+                                                                * only */
+
+       }
+       else
+       {
+               s_usqr(da, dc, size_a);
+       }
+
+       return 1;
+}
+
+/* }}} */
+
+/* {{{ s_usqr(da, dc, size_a) */
+
+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);
+       }
+}
+
+/* }}} */
+
+/* {{{ s_dadd(a, b) */
+
+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;
+               MP_USED(a) += 1;
+       }
+}
+
+/* }}} */
+
+/* {{{ s_dmul(a, b) */
+
+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;
+               MP_USED(a) += 1;
+       }
+}
+
+/* }}} */
+
+/* {{{ s_dbmul(da, b, dc, size_a) */
+
+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);
+}
+
+/* }}} */
+
+/* {{{ s_ddiv(da, d, dc, size_a) */
+
+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;
+}
+
+/* }}} */
+
+/* {{{ s_qdiv(z, p2) */
+
+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++;
+
+               MP_USED(z) = 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)
+               MP_SIGN(z) = MP_ZPOS;
+}
+
+/* }}} */
+
+/* {{{ s_qmod(z, p2) */
+
+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 = (1 << rest) - 1;
+
+       if (start <= uz)
+       {
+               MP_USED(z) = start;
+               z->digits[start - 1] &= mask;
+               CLAMP(z);
+       }
+}
+
+/* }}} */
+
+/* {{{ s_qmul(z, p2) */
+
+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;
+               }
+       }
+
+       MP_USED(z) = uz;
+       CLAMP(z);
+
+       return 1;
+}
+
+/* }}} */
+
+/* {{{ s_qsub(z, p2) */
+
+/* Subtract |z| from 2^p2, assuming 2^p2 > |z|, and set z to be positive */
+static int
+s_qsub(mp_int z, mp_size p2)
+{
+       mp_digit        hi = (1 << (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 */
+
+       MP_SIGN(z) = MP_ZPOS;
+       CLAMP(z);
+
+       return 1;
+}
+
+/* }}} */
+
+/* {{{ s_dp2k(z) */
+
+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;
+}
+
+/* }}} */
+
+/* {{{ s_isp2(z) */
+
+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;
+}
+
+/* }}} */
+
+/* {{{ s_2expt(z, k) */
+
+static int
+s_2expt(mp_int z, int 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) = (1 << rest);
+       MP_USED(z) = ndig;
+
+       return 1;
+}
+
+/* }}} */
+
+/* {{{ s_norm(a, b) */
+
+static int
+s_norm(mp_int a, mp_int b)
+{
+       mp_digit        d = b->digits[MP_USED(b) - 1];
+       int                     k = 0;
+
+       while (d < (mp_digit) ((mp_digit) 1 << (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;
+}
+
+/* }}} */
+
+/* {{{ s_brmu(z, m) */
+
+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);
+}
+
+/* }}} */
+
+/* {{{ s_reduce(x, m, mu, q1, q2) */
+
+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;
+}
+
+/* }}} */
+
+/* {{{ s_embar(a, b, m, mu, c) */
+
+/* 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   *db,
+                          *dbt,
+                               umu,
+                               d;
+       mpz_t           temp[3];
+       mp_result       res;
+       int                     last = 0;
+
+       umu = MP_USED(mu);
+       db = MP_DIGITS(b);
+       dbt = db + MP_USED(b) - 1;
+
+       while (last < 3)
+               SETUP(mp_int_init_size(TEMP(last), 2 * umu), last);
+
+       (void) mp_int_set_value(c, 1);
+
+       /* Take care of low-order digits */
+       while (db < dbt)
+       {
+               int                     i;
+
+               for (d = *db, 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)))
+                               {
+                                       res = MP_MEMORY;
+                                       goto CLEANUP;
+                               }
+                               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)))
+                       {
+                               res = MP_MEMORY;
+                               goto CLEANUP;
+                       }
+                       assert(MP_SIGN(TEMP(0)) == MP_ZPOS);
+                       mp_int_copy(TEMP(0), a);
+
+
+               }
+
+               ++db;
+       }
+
+       /* Take care of highest-order digit */
+       d = *dbt;
+       for (;;)
+       {
+               if (d & 1)
+               {
+                       UMUL(c, a, TEMP(0));
+                       if (!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2)))
+                       {
+                               res = MP_MEMORY;
+                               goto CLEANUP;
+                       }
+                       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)))
+               {
+                       res = MP_MEMORY;
+                       goto CLEANUP;
+               }
+               (void) mp_int_copy(TEMP(0), a);
+       }
+
+CLEANUP:
+       while (--last >= 0)
+               mp_int_clear(TEMP(last));
+
+       return res;
+}
+
+/* }}} */
+
+/* {{{ s_udiv(a, b) */
+
+/* Precondition:  a >= b and b > 0
+   Postcondition: a' = a / b, b' = a % b
+ */
+static mp_result
+s_udiv(mp_int a, mp_int b)
+{
+       mpz_t           q,
+                               r,
+                               t;
+       mp_size         ua,
+                               ub,
+                               qpos = 0;
+       mp_digit   *da,
+                               btop;
+       mp_result       res = MP_OK;
+       int                     k,
+                               skip = 0;
+
+       /* Force signs to positive */
+       MP_SIGN(a) = MP_ZPOS;
+       MP_SIGN(b) = MP_ZPOS;
+
+       /* Normalize, per Knuth */
+       k = s_norm(a, b);
+
+       ua = MP_USED(a);
+       ub = MP_USED(b);
+       btop = b->digits[ub - 1];
+       if ((res = mp_int_init_size(&q, ua)) != MP_OK)
+               return res;
+       if ((res = mp_int_init_size(&t, ua + 1)) != MP_OK)
+               goto CLEANUP;
+
+       da = MP_DIGITS(a);
+       r.digits = da + ua - 1;         /* The contents of r are shared with a 
*/
+       r.used = 1;
+       r.sign = MP_ZPOS;
+       r.alloc = MP_ALLOC(a);
+       ZERO(t.digits, t.alloc);
+
+       /* Solve for quotient digits, store in q.digits in reverse order */
+       while (r.digits >= da)
+       {
+               assert(qpos <= q.alloc);
+
+               if (s_ucmp(b, &r) > 0)
+               {
+                       r.digits -= 1;
+                       r.used += 1;
+
+                       if (++skip > 1)
+                               q.digits[qpos++] = 0;
+
+                       CLAMP(&r);
+               }
+               else
+               {
+                       mp_word         pfx = r.digits[r.used - 1];
+                       mp_word         qdigit;
+
+                       if (r.used > 1 && (pfx < btop || r.digits[r.used - 2] 
== 0))
+                       {
+                               pfx <<= MP_DIGIT_BIT / 2;
+                               pfx <<= MP_DIGIT_BIT / 2;
+                               pfx |= r.digits[r.used - 2];
+                       }
+
+                       qdigit = pfx / btop;
+                       if (qdigit > MP_DIGIT_MAX)
+                               qdigit = 1;
+
+                       s_dbmul(MP_DIGITS(b), (mp_digit) qdigit, t.digits, ub);
+                       t.used = ub + 1;
+                       CLAMP(&t);
+                       while (s_ucmp(&t, &r) > 0)
+                       {
+                               --qdigit;
+                               (void) mp_int_sub(&t, b, &t);   /* cannot fail 
*/
+                       }
+
+                       s_usub(r.digits, t.digits, r.digits, r.used, t.used);
+                       CLAMP(&r);
+
+                       q.digits[qpos++] = (mp_digit) qdigit;
+                       ZERO(t.digits, t.used);
+                       skip = 0;
+               }
+       }
+
+       /* Put quotient digits in the correct order, and discard extra zeroes */
+       q.used = qpos;
+       REV(mp_digit, q.digits, qpos);
+       CLAMP(&q);
+
+       /* Denormalize the remainder */
+       CLAMP(a);
+       if (k != 0)
+               s_qdiv(a, k);
+
+       mp_int_copy(a, b);                      /* ok:  0 <= r < b */
+       mp_int_copy(&q, a);                     /* ok:  q <= a     */
+
+       mp_int_clear(&t);
+CLEANUP:
+       mp_int_clear(&q);
+       return res;
+}
+
+/* }}} */
+
+/* {{{ s_outlen(z, r) */
+
+/* Precondition:  2 <= r < 64 */
+static int
+s_outlen(mp_int z, mp_size r)
+{
+       mp_result       bits;
+       double          raw;
+
+       bits = mp_int_count_bits(z);
+       raw = (double) bits *s_log2[r];
+
+       return (int) (raw + 0.999999);
+}
+
+/* }}} */
+
+/* {{{ s_inlen(len, r) */
+
+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);
+}
+
+/* }}} */
+
+/* {{{ s_ch2val(c, r) */
+
+static int
+s_ch2val(char c, int r)
+{
+       int                     out;
+
+       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;
+}
+
+/* }}} */
+
+/* {{{ s_val2ch(v, caps) */
+
+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;
+       }
+}
+
+/* }}} */
+
+/* {{{ s_2comp(buf, len) */
+
+static void
+s_2comp(unsigned char *buf, int len)
+{
+       int                     i;
+       unsigned short s = 1;
+
+       for (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 */
+}
+
+/* }}} */
+
+/* {{{ s_tobin(z, buf, *limpos) */
+
+static mp_result
+s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad)
+{
+       mp_size         uz;
+       mp_digit   *dz;
+       int                     pos = 0,
+                               limit = *limpos;
+
+       uz = MP_USED(z);
+       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(unsigned char, buf, pos);
+
+       /* Return the number of bytes actually written */
+       *limpos = pos;
+
+       return (uz == 0) ? MP_OK : MP_TRUNC;
+}
+
+/* }}} */
+
+/* {{{ s_print(tag, z) */
+
+#if 0
+void
+s_print(char *tag, mp_int z)
+{
+       int                     i;
+
+       fprintf(stderr, "%s: %c ", tag,
+                       (MP_SIGN(z) == MP_NEG) ? '-' : '+');
+
+       for (i = MP_USED(z) - 1; i >= 0; --i)
+               fprintf(stderr, "%0*X", (int) (MP_DIGIT_BIT / 4), z->digits[i]);
+
+       fputc('\n', stderr);
+
+}
+
+void
+s_print_buf(char *tag, mp_digit *buf, mp_size num)
+{
+       int                     i;
+
+       fprintf(stderr, "%s: ", tag);
+
+       for (i = num - 1; i >= 0; --i)
+               fprintf(stderr, "%0*X", (int) (MP_DIGIT_BIT / 4), buf[i]);
+
+       fputc('\n', stderr);
+}
+#endif
+
+/* }}} */
+
+/* HERE THERE BE DRAGONS */

http://git-wip-us.apache.org/repos/asf/incubator-hawq/blob/801100ed/contrib/pgcrypto/imath.h
----------------------------------------------------------------------
diff --git a/contrib/pgcrypto/imath.h b/contrib/pgcrypto/imath.h
new file mode 100644
index 0000000..f2b02d0
--- /dev/null
+++ b/contrib/pgcrypto/imath.h
@@ -0,0 +1,217 @@
+/*
+  Name:                imath.h
+  Purpose:     Arbitrary precision integer arithmetic routines.
+  Author:      M. J. Fromberger <http://spinning-yarns.org/michael/sw/>
+  Info:                Id: imath.h 21 2006-04-02 18:58:36Z sting
+
+  Copyright (C) 2002 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.
+ */
+/* contrib/pgcrypto/imath.h */
+
+#ifndef IMATH_H_
+#define IMATH_H_
+
+/* use always 32bit digits - should some arch use 16bit digits? */
+#define USE_LONG_LONG
+
+#include <limits.h>
+
+typedef unsigned char mp_sign;
+typedef unsigned int mp_size;
+typedef int mp_result;
+
+#ifdef USE_LONG_LONG
+typedef uint32 mp_digit;
+typedef uint64 mp_word;
+
+#define MP_DIGIT_MAX      0xFFFFFFFFULL
+#define MP_WORD_MAX               0xFFFFFFFFFFFFFFFFULL
+#else
+typedef uint16 mp_digit;
+typedef uint32 mp_word;
+
+#define MP_DIGIT_MAX      0xFFFFUL
+#define MP_WORD_MAX               0xFFFFFFFFUL
+#endif
+
+typedef struct mpz
+{
+       mp_digit   *digits;
+       mp_size         alloc;
+       mp_size         used;
+       mp_sign         sign;
+} mpz_t    ,
+                  *mp_int;
+
+#define MP_DIGITS(Z) ((Z)->digits)
+#define MP_ALLOC(Z)  ((Z)->alloc)
+#define MP_USED(Z)      ((Z)->used)
+#define MP_SIGN(Z)      ((Z)->sign)
+
+extern const mp_result MP_OK;
+extern const mp_result MP_FALSE;
+extern const mp_result MP_TRUE;
+extern const mp_result MP_MEMORY;
+extern const mp_result MP_RANGE;
+extern const mp_result MP_UNDEF;
+extern const mp_result MP_TRUNC;
+extern const mp_result MP_BADARG;
+
+#define MP_DIGIT_BIT   (sizeof(mp_digit) * CHAR_BIT)
+#define MP_WORD_BIT            (sizeof(mp_word) * CHAR_BIT)
+
+#define MP_MIN_RADIX   2
+#define MP_MAX_RADIX   36
+
+extern const mp_sign MP_NEG;
+extern const mp_sign MP_ZPOS;
+
+#define mp_int_is_odd(Z)  ((Z)->digits[0] & 1)
+#define mp_int_is_even(Z) !((Z)->digits[0] & 1)
+
+mp_size                mp_get_default_precision(void);
+void           mp_set_default_precision(mp_size s);
+mp_size                mp_get_multiply_threshold(void);
+void           mp_set_multiply_threshold(mp_size s);
+
+mp_result      mp_int_init(mp_int z);
+mp_int         mp_int_alloc(void);
+mp_result      mp_int_init_size(mp_int z, mp_size prec);
+mp_result      mp_int_init_copy(mp_int z, mp_int old);
+mp_result      mp_int_init_value(mp_int z, int value);
+mp_result      mp_int_set_value(mp_int z, int value);
+void           mp_int_clear(mp_int z);
+void           mp_int_free(mp_int z);
+
+mp_result      mp_int_copy(mp_int a, mp_int c);        /* c = a         */
+void           mp_int_swap(mp_int a, mp_int c);        /* swap a, c */
+void           mp_int_zero(mp_int z);          /* z = 0         */
+mp_result      mp_int_abs(mp_int a, mp_int c);         /* c = |a|       */
+mp_result      mp_int_neg(mp_int a, mp_int c);         /* c = -a        */
+mp_result      mp_int_add(mp_int a, mp_int b, mp_int c);       /* c = a + b */
+mp_result      mp_int_add_value(mp_int a, int value, mp_int c);
+mp_result      mp_int_sub(mp_int a, mp_int b, mp_int c);       /* c = a - b */
+mp_result      mp_int_sub_value(mp_int a, int value, mp_int c);
+mp_result      mp_int_mul(mp_int a, mp_int b, mp_int c);       /* c = a * b */
+mp_result      mp_int_mul_value(mp_int a, int value, mp_int c);
+mp_result      mp_int_mul_pow2(mp_int a, int p2, mp_int c);
+mp_result      mp_int_sqr(mp_int a, mp_int c);         /* c = a * a */
+
+mp_result
+mp_int_div(mp_int a, mp_int b, /* q = a / b */
+                  mp_int q, mp_int r); /* r = a % b */
+mp_result
+mp_int_div_value(mp_int a, int value,  /* q = a / value */
+                                mp_int q, int *r);             /* r = a % 
value */
+mp_result
+mp_int_div_pow2(mp_int a, int p2,              /* q = a / 2^p2  */
+                               mp_int q, mp_int r);    /* r = q % 2^p2  */
+mp_result      mp_int_mod(mp_int a, mp_int m, mp_int c);       /* c = a % m */
+
+#define   mp_int_mod_value(A, V, R) mp_int_div_value((A), (V), 0, (R))
+mp_result      mp_int_expt(mp_int a, int b, mp_int c);         /* c = a^b      
 */
+mp_result      mp_int_expt_value(int a, int b, mp_int c);      /* c = a^b      
 */
+
+int                    mp_int_compare(mp_int a, mp_int b); /* a <=> b     */
+int                    mp_int_compare_unsigned(mp_int a, mp_int b);            
/* |a| <=> |b| */
+int                    mp_int_compare_zero(mp_int z);          /* a <=> 0      
   */
+int                    mp_int_compare_value(mp_int z, int value);      /* a 
<=> v         */
+
+/* Returns true if v|a, false otherwise (including errors) */
+int                    mp_int_divisible_value(mp_int a, int v);
+
+/* Returns k >= 0 such that z = 2^k, if one exists; otherwise < 0 */
+int                    mp_int_is_pow2(mp_int z);
+
+mp_result
+mp_int_exptmod(mp_int a, mp_int b, mp_int m,
+                          mp_int c);           /* c = a^b (mod m) */
+mp_result
+mp_int_exptmod_evalue(mp_int a, int value,
+                                         mp_int m, mp_int c);          /* c = 
a^v (mod m) */
+mp_result
+mp_int_exptmod_bvalue(int value, mp_int b,
+                                         mp_int m, mp_int c);          /* c = 
v^b (mod m) */
+mp_result
+mp_int_exptmod_known(mp_int a, mp_int b,
+                                        mp_int m, mp_int mu,
+                                        mp_int c); /* c = a^b (mod m) */
+mp_result      mp_int_redux_const(mp_int m, mp_int c);
+
+mp_result      mp_int_invmod(mp_int a, mp_int m, mp_int c);            /* c = 
1/a (mod m) */
+
+mp_result      mp_int_gcd(mp_int a, mp_int b, mp_int c);       /* c = gcd(a, 
b)   */
+
+mp_result
+mp_int_egcd(mp_int a, mp_int b, mp_int c,              /* c = gcd(a, b)   */
+                       mp_int x, mp_int y);    /* c = ax + by     */
+
+mp_result      mp_int_sqrt(mp_int a, mp_int c);        /* c = floor(sqrt(q)) */
+
+/* Convert to an int, if representable (returns MP_RANGE if not). */
+mp_result      mp_int_to_int(mp_int z, int *out);
+
+/* Convert to nul-terminated string with the specified radix, writing at
+   most limit characters including the nul terminator  */
+mp_result mp_int_to_string(mp_int z, mp_size radix,
+                                char *str, int limit);
+
+/* Return the number of characters required to represent
+   z in the given radix.  May over-estimate. */
+mp_result      mp_int_string_len(mp_int z, mp_size radix);
+
+/* Read zero-terminated string into z */
+mp_result      mp_int_read_string(mp_int z, mp_size radix, const char *str);
+mp_result mp_int_read_cstring(mp_int z, mp_size radix, const char *str,
+                                       char **end);
+
+/* Return the number of significant bits in z */
+mp_result      mp_int_count_bits(mp_int z);
+
+/* Convert z to two's complement binary, writing at most limit bytes */
+mp_result      mp_int_to_binary(mp_int z, unsigned char *buf, int limit);
+
+/* Read a two's complement binary value into z from the given buffer */
+mp_result      mp_int_read_binary(mp_int z, unsigned char *buf, int len);
+
+/* Return the number of bytes required to represent z in binary. */
+mp_result      mp_int_binary_len(mp_int z);
+
+/* Convert z to unsigned binary, writing at most limit bytes */
+mp_result      mp_int_to_unsigned(mp_int z, unsigned char *buf, int limit);
+
+/* Read an unsigned binary value into z from the given buffer */
+mp_result      mp_int_read_unsigned(mp_int z, unsigned char *buf, int len);
+
+/* Return the number of bytes required to represent z as unsigned output */
+mp_result      mp_int_unsigned_len(mp_int z);
+
+/* Return a statically allocated string describing error code res */
+const char *mp_error_string(mp_result res);
+
+#if 0
+void           s_print(char *tag, mp_int z);
+void           s_print_buf(char *tag, mp_digit *buf, mp_size num);
+#endif
+
+#endif   /* end IMATH_H_ */

http://git-wip-us.apache.org/repos/asf/incubator-hawq/blob/801100ed/contrib/pgcrypto/internal-sha2.c
----------------------------------------------------------------------
diff --git a/contrib/pgcrypto/internal-sha2.c b/contrib/pgcrypto/internal-sha2.c
new file mode 100644
index 0000000..f86b478
--- /dev/null
+++ b/contrib/pgcrypto/internal-sha2.c
@@ -0,0 +1,316 @@
+/*
+ * internal.c
+ *             Wrapper for builtin functions
+ *
+ * Copyright (c) 2001 Marko Kreen
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *       notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *       notice, this list of conditions and the following disclaimer in the
+ *       documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.     IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * contrib/pgcrypto/internal-sha2.c
+ */
+
+#include "postgres.h"
+
+#include <time.h>
+
+#include "px.h"
+#include "sha2.h"
+
+void           init_sha224(PX_MD *h);
+void           init_sha256(PX_MD *h);
+void           init_sha384(PX_MD *h);
+void           init_sha512(PX_MD *h);
+
+/* SHA224 */
+
+static unsigned
+int_sha224_len(PX_MD *h)
+{
+       return SHA224_DIGEST_LENGTH;
+}
+
+static unsigned
+int_sha224_block_len(PX_MD *h)
+{
+       return SHA224_BLOCK_LENGTH;
+}
+
+static void
+int_sha224_update(PX_MD *h, const uint8 *data, unsigned dlen)
+{
+       SHA224_CTX *ctx = (SHA224_CTX *) h->p.ptr;
+
+       SHA224_Update(ctx, data, dlen);
+}
+
+static void
+int_sha224_reset(PX_MD *h)
+{
+       SHA224_CTX *ctx = (SHA224_CTX *) h->p.ptr;
+
+       SHA224_Init(ctx);
+}
+
+static void
+int_sha224_finish(PX_MD *h, uint8 *dst)
+{
+       SHA224_CTX *ctx = (SHA224_CTX *) h->p.ptr;
+
+       SHA224_Final(dst, ctx);
+}
+
+static void
+int_sha224_free(PX_MD *h)
+{
+       SHA224_CTX *ctx = (SHA224_CTX *) h->p.ptr;
+
+       memset(ctx, 0, sizeof(*ctx));
+       px_free(ctx);
+       px_free(h);
+}
+
+/* SHA256 */
+
+static unsigned
+int_sha256_len(PX_MD *h)
+{
+       return SHA256_DIGEST_LENGTH;
+}
+
+static unsigned
+int_sha256_block_len(PX_MD *h)
+{
+       return SHA256_BLOCK_LENGTH;
+}
+
+static void
+int_sha256_update(PX_MD *h, const uint8 *data, unsigned dlen)
+{
+       SHA256_CTX *ctx = (SHA256_CTX *) h->p.ptr;
+
+       SHA256_Update(ctx, data, dlen);
+}
+
+static void
+int_sha256_reset(PX_MD *h)
+{
+       SHA256_CTX *ctx = (SHA256_CTX *) h->p.ptr;
+
+       SHA256_Init(ctx);
+}
+
+static void
+int_sha256_finish(PX_MD *h, uint8 *dst)
+{
+       SHA256_CTX *ctx = (SHA256_CTX *) h->p.ptr;
+
+       SHA256_Final(dst, ctx);
+}
+
+static void
+int_sha256_free(PX_MD *h)
+{
+       SHA256_CTX *ctx = (SHA256_CTX *) h->p.ptr;
+
+       memset(ctx, 0, sizeof(*ctx));
+       px_free(ctx);
+       px_free(h);
+}
+
+/* SHA384 */
+
+static unsigned
+int_sha384_len(PX_MD *h)
+{
+       return SHA384_DIGEST_LENGTH;
+}
+
+static unsigned
+int_sha384_block_len(PX_MD *h)
+{
+       return SHA384_BLOCK_LENGTH;
+}
+
+static void
+int_sha384_update(PX_MD *h, const uint8 *data, unsigned dlen)
+{
+       SHA384_CTX *ctx = (SHA384_CTX *) h->p.ptr;
+
+       SHA384_Update(ctx, data, dlen);
+}
+
+static void
+int_sha384_reset(PX_MD *h)
+{
+       SHA384_CTX *ctx = (SHA384_CTX *) h->p.ptr;
+
+       SHA384_Init(ctx);
+}
+
+static void
+int_sha384_finish(PX_MD *h, uint8 *dst)
+{
+       SHA384_CTX *ctx = (SHA384_CTX *) h->p.ptr;
+
+       SHA384_Final(dst, ctx);
+}
+
+static void
+int_sha384_free(PX_MD *h)
+{
+       SHA384_CTX *ctx = (SHA384_CTX *) h->p.ptr;
+
+       memset(ctx, 0, sizeof(*ctx));
+       px_free(ctx);
+       px_free(h);
+}
+
+/* SHA512 */
+
+static unsigned
+int_sha512_len(PX_MD *h)
+{
+       return SHA512_DIGEST_LENGTH;
+}
+
+static unsigned
+int_sha512_block_len(PX_MD *h)
+{
+       return SHA512_BLOCK_LENGTH;
+}
+
+static void
+int_sha512_update(PX_MD *h, const uint8 *data, unsigned dlen)
+{
+       SHA512_CTX *ctx = (SHA512_CTX *) h->p.ptr;
+
+       SHA512_Update(ctx, data, dlen);
+}
+
+static void
+int_sha512_reset(PX_MD *h)
+{
+       SHA512_CTX *ctx = (SHA512_CTX *) h->p.ptr;
+
+       SHA512_Init(ctx);
+}
+
+static void
+int_sha512_finish(PX_MD *h, uint8 *dst)
+{
+       SHA512_CTX *ctx = (SHA512_CTX *) h->p.ptr;
+
+       SHA512_Final(dst, ctx);
+}
+
+static void
+int_sha512_free(PX_MD *h)
+{
+       SHA512_CTX *ctx = (SHA512_CTX *) h->p.ptr;
+
+       memset(ctx, 0, sizeof(*ctx));
+       px_free(ctx);
+       px_free(h);
+}
+
+/* init functions */
+
+void
+init_sha224(PX_MD *md)
+{
+       SHA224_CTX *ctx;
+
+       ctx = px_alloc(sizeof(*ctx));
+       memset(ctx, 0, sizeof(*ctx));
+
+       md->p.ptr = ctx;
+
+       md->result_size = int_sha224_len;
+       md->block_size = int_sha224_block_len;
+       md->reset = int_sha224_reset;
+       md->update = int_sha224_update;
+       md->finish = int_sha224_finish;
+       md->free = int_sha224_free;
+
+       md->reset(md);
+}
+
+void
+init_sha256(PX_MD *md)
+{
+       SHA256_CTX *ctx;
+
+       ctx = px_alloc(sizeof(*ctx));
+       memset(ctx, 0, sizeof(*ctx));
+
+       md->p.ptr = ctx;
+
+       md->result_size = int_sha256_len;
+       md->block_size = int_sha256_block_len;
+       md->reset = int_sha256_reset;
+       md->update = int_sha256_update;
+       md->finish = int_sha256_finish;
+       md->free = int_sha256_free;
+
+       md->reset(md);
+}
+
+void
+init_sha384(PX_MD *md)
+{
+       SHA384_CTX *ctx;
+
+       ctx = px_alloc(sizeof(*ctx));
+       memset(ctx, 0, sizeof(*ctx));
+
+       md->p.ptr = ctx;
+
+       md->result_size = int_sha384_len;
+       md->block_size = int_sha384_block_len;
+       md->reset = int_sha384_reset;
+       md->update = int_sha384_update;
+       md->finish = int_sha384_finish;
+       md->free = int_sha384_free;
+
+       md->reset(md);
+}
+
+void
+init_sha512(PX_MD *md)
+{
+       SHA512_CTX *ctx;
+
+       ctx = px_alloc(sizeof(*ctx));
+       memset(ctx, 0, sizeof(*ctx));
+
+       md->p.ptr = ctx;
+
+       md->result_size = int_sha512_len;
+       md->block_size = int_sha512_block_len;
+       md->reset = int_sha512_reset;
+       md->update = int_sha512_update;
+       md->finish = int_sha512_finish;
+       md->free = int_sha512_free;
+
+       md->reset(md);
+}

Reply via email to