s4:heimdal: import lorikeet-heimdal-201011102149 (commit 5734d03c20e104c8f45533d07f2a2cbbd3224f29)

15 files changed, 99 insertions, 4892 deletions

diff --git a/source4/heimdal/lib/hcrypto/dh-imath.c b/source4/heimdal/lib/hcrypto/dh-imath.c
deleted file mode 100644
index c2e86fa2fa..0000000000
--- a/source4/heimdal/lib/hcrypto/dh-imath.c
+++ /dev/null
@@ -1,254 +0,0 @@
-/*
- * Copyright (c) 2006 Kungliga Tekniska Högskolan
- * (Royal Institute of Technology, Stockholm, Sweden).
- * 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.
- *
- * 3. Neither the name of the Institute nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE INSTITUTE 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 INSTITUTE 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.
- */
-
-#include <config.h>
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <dh.h>
-
-#ifdef USE_HCRYPTO_IMATH
-
-#include <roken.h>
-
-#include "imath/imath.h"
-
-static void
-BN2mpz(mpz_t *s, const BIGNUM *bn)
-{
-    size_t len;
-    void *p;
-
-    len = BN_num_bytes(bn);
-    p = malloc(len);
-    BN_bn2bin(bn, p);
-    mp_int_read_unsigned(s, p, len);
-    free(p);
-}
-
-
-static BIGNUM *
-mpz2BN(mpz_t *s)
-{
-    size_t size;
-    BIGNUM *bn;
-    void *p;
-
-    size = mp_int_unsigned_len(s);
-    p = malloc(size);
-    if (p == NULL && size != 0)
-	return NULL;
-    mp_int_to_unsigned(s, p, size);
-
-    bn = BN_bin2bn(p, size, NULL);
-    free(p);
-    return bn;
-}
-
-/*
- *
- */
-
-#define DH_NUM_TRIES 10
-
-static int
-dh_generate_key(DH *dh)
-{
-    mpz_t pub, priv_key, g, p;
-    int have_private_key = (dh->priv_key != NULL);
-    int codes, times = 0;
-    mp_result res;
-
-    if (dh->p == NULL || dh->g == NULL)
-	return 0;
-
-    while (times++ < DH_NUM_TRIES) {
-	if (!have_private_key) {
-	    size_t bits = BN_num_bits(dh->p);
-
-	    if (dh->priv_key)
-		BN_free(dh->priv_key);
-
-	    dh->priv_key = BN_new();
-	    if (dh->priv_key == NULL)
-		return 0;
-	    if (!BN_rand(dh->priv_key, bits - 1, 0, 0)) {
-		BN_clear_free(dh->priv_key);
-		dh->priv_key = NULL;
-		return 0;
-	    }
-	}
-	if (dh->pub_key)
-	    BN_free(dh->pub_key);
-
-	mp_int_init(&pub);
-	mp_int_init(&priv_key);
-	mp_int_init(&g);
-	mp_int_init(&p);
-	
-	BN2mpz(&priv_key, dh->priv_key);
-	BN2mpz(&g, dh->g);
-	BN2mpz(&p, dh->p);
-	
-	res = mp_int_exptmod(&g, &priv_key, &p, &pub);
-
-	mp_int_clear(&priv_key);
-	mp_int_clear(&g);
-	mp_int_clear(&p);
-	if (res != MP_OK)
-	    continue;
-
-	dh->pub_key = mpz2BN(&pub);
-	mp_int_clear(&pub);
-	if (dh->pub_key == NULL)
-	    return 0;
-	
-	if (DH_check_pubkey(dh, dh->pub_key, &codes) && codes == 0)
-	    break;
-	if (have_private_key)
-	    return 0;
-    }
-
-    if (times >= DH_NUM_TRIES) {
-	if (!have_private_key && dh->priv_key) {
-	    BN_free(dh->priv_key);
-	    dh->priv_key = NULL;
-	}
-	if (dh->pub_key) {
-	    BN_free(dh->pub_key);
-	    dh->pub_key = NULL;
-	}
-	return 0;
-    }
-
-    return 1;
-}
-
-static int
-dh_compute_key(unsigned char *shared, const BIGNUM * pub, DH *dh)
-{
-    mpz_t s, priv_key, p, peer_pub;
-    size_t size = 0;
-    mp_result res;
-
-    if (dh->pub_key == NULL || dh->g == NULL || dh->priv_key == NULL)
-	return -1;
-
-    mp_int_init(&p);
-    BN2mpz(&p, dh->p);
-
-    mp_int_init(&peer_pub);
-    BN2mpz(&peer_pub, pub);
-
-    /* check if peers pubkey is reasonable */
-    if (MP_SIGN(&peer_pub) == MP_NEG
-	|| mp_int_compare(&peer_pub, &p) >= 0
-	|| mp_int_compare_value(&peer_pub, 1) <= 0)
-    {
-	mp_int_clear(&p);
-	mp_int_clear(&peer_pub);
-	return -1;
-    }
-
-    mp_int_init(&priv_key);
-    BN2mpz(&priv_key, dh->priv_key);
-
-    mp_int_init(&s);
-
-    mp_int_exptmod(&peer_pub, &priv_key, &p, &s);
-
-    mp_int_clear(&p);
-    mp_int_clear(&peer_pub);
-    mp_int_clear(&priv_key);
-
-    size = mp_int_unsigned_len(&s);
-    res = mp_int_to_unsigned(&s, shared, size);
-    mp_int_clear(&s);
-
-    return (res == MP_OK) ? size : -1;
-}
-
-static int
-dh_generate_params(DH *dh, int a, int b, BN_GENCB *callback)
-{
-    /* groups should already be known, we don't care about this */
-    return 0;
-}
-
-static int
-dh_init(DH *dh)
-{
-    return 1;
-}
-
-static int
-dh_finish(DH *dh)
-{
-    return 1;
-}
-
-
-/*
- *
- */
-
-const DH_METHOD _hc_dh_imath_method = {
-    "hcrypto imath DH",
-    dh_generate_key,
-    dh_compute_key,
-    NULL,
-    dh_init,
-    dh_finish,
-    0,
-    NULL,
-    dh_generate_params
-};
-#endif /* USE_HCRYPTO_DH_IMATH */
-
-/**
- * DH implementation using libimath.
- *
- * @return the DH_METHOD for the DH implementation using libimath.
- *
- * @ingroup hcrypto_dh
- */
-
-const DH_METHOD *
-DH_imath_method(void)
-{
-#ifdef USE_HCRYPTO_DH_IMATH
-    return &_hc_dh_imath_method;
-#else
-    return NULL;
-#endif
-}
diff --git a/source4/heimdal/lib/hcrypto/dh-ltm.c b/source4/heimdal/lib/hcrypto/dh-ltm.c
index 6166100b08..f66cd5aff2 100644
--- a/source4/heimdal/lib/hcrypto/dh-ltm.c
+++ b/source4/heimdal/lib/hcrypto/dh-ltm.c
@@ -119,14 +119,12 @@ ltm_dh_generate_key(DH *dh)
 	
 	res = mp_exptmod(&g, &priv_key, &p, &pub);
 
-	mp_zero(&priv_key);
-	mp_zero(&g);
-	mp_zero(&p);
+	mp_clear_multi(&priv_key, &g, &p, NULL);
 	if (res != 0)
 	    continue;
 
 	dh->pub_key = mpz2BN(&pub);
-	mp_zero(&pub);
+	mp_clear(&pub);
 	if (dh->pub_key == NULL)
 	    return 0;
 	
@@ -155,16 +153,13 @@ static int
 ltm_dh_compute_key(unsigned char *shared, const BIGNUM * pub, DH *dh)
 {
     mp_int s, priv_key, p, peer_pub;
-    size_t size = 0;
     int ret;
 
     if (dh->pub_key == NULL || dh->g == NULL || dh->priv_key == NULL)
 	return -1;
 
-    mp_init(&p);
+    mp_init_multi(&s, &priv_key, &p, &peer_pub, NULL);
     BN2mpz(&p, dh->p);
-
-    mp_init(&peer_pub);
     BN2mpz(&peer_pub, pub);
 
     /* check if peers pubkey is reasonable */
@@ -172,30 +167,26 @@ ltm_dh_compute_key(unsigned char *shared, const BIGNUM * pub, DH *dh)
 	|| mp_cmp(&peer_pub, &p) >= 0
 	|| mp_cmp_d(&peer_pub, 1) <= 0)
     {
-	mp_zero(&p);
-	mp_zero(&peer_pub);
-	return -1;
+	ret = -1;
+	goto out;
     }
 
-    mp_init(&priv_key);
     BN2mpz(&priv_key, dh->priv_key);
 
-    mp_init(&s);
-
     ret = mp_exptmod(&peer_pub, &priv_key, &p, &s);
 
-    mp_zero(&p);
-    mp_zero(&peer_pub);
-    mp_zero(&priv_key);
-
-    if (ret != 0)
-	return -1;
+    if (ret != 0) {
+	ret = -1;
+	goto out;
+    }
 
-    size = mp_unsigned_bin_size(&s);
+    ret = mp_unsigned_bin_size(&s);
     mp_to_unsigned_bin(&s, shared);
-    mp_zero(&s);
 
-    return size;
+ out:
+    mp_clear_multi(&s, &priv_key, &p, &peer_pub, NULL);
+
+    return ret;
 }
 
 static int
@@ -235,9 +226,9 @@ const DH_METHOD _hc_dh_ltm_method = {
 };
 
 /**
- * DH implementation using libimath.
+ * DH implementation using libtommath.
  *
- * @return the DH_METHOD for the DH implementation using libimath.
+ * @return the DH_METHOD for the DH implementation using libtommath.
  *
  * @ingroup hcrypto_dh
  */
diff --git a/source4/heimdal/lib/hcrypto/dh.c b/source4/heimdal/lib/hcrypto/dh.c
index 3ad37f87a7..43e1d6ac1b 100644
--- a/source4/heimdal/lib/hcrypto/dh.c
+++ b/source4/heimdal/lib/hcrypto/dh.c
@@ -304,7 +304,7 @@ DH_check_pubkey(const DH *dh, const BIGNUM *pub_key, int *codes)
     if (!BN_set_word(bn, 2))
 	goto out;
 
-    if (BN_cmp(bn, pub_key) == 0) {
+    if (BN_cmp(bn, dh->g) == 0) {
 	unsigned i, n = BN_num_bits(pub_key);
 	unsigned bits = 0;
 
@@ -312,7 +312,7 @@ DH_check_pubkey(const DH *dh, const BIGNUM *pub_key, int *codes)
 	    if (BN_is_bit_set(pub_key, i))
 		bits++;
 
-	if (bits > 1) {
+	if (bits < 2) {
 	    *codes |= DH_CHECK_PUBKEY_TOO_SMALL;
 	    goto out;
 	}
@@ -445,8 +445,8 @@ static const DH_METHOD dh_null_method = {
     dh_null_generate_params
 };
 
-extern const DH_METHOD _hc_dh_imath_method;
-static const DH_METHOD *dh_default_method = &_hc_dh_imath_method;
+extern const DH_METHOD _hc_dh_ltm_method;
+static const DH_METHOD *dh_default_method = &_hc_dh_ltm_method;
 
 /**
  * Return the dummy DH implementation.
diff --git a/source4/heimdal/lib/hcrypto/dh.h b/source4/heimdal/lib/hcrypto/dh.h
index 3a24f9dfdf..637f218bcf 100644
--- a/source4/heimdal/lib/hcrypto/dh.h
+++ b/source4/heimdal/lib/hcrypto/dh.h
@@ -40,7 +40,6 @@
 
 /* symbol renaming */
 #define DH_null_method hc_DH_null_method
-#define DH_imath_method hc_DH_imath_method
 #define DH_tfm_method hc_DH_tfm_method
 #define DH_ltm_method hc_DH_ltm_method
 #define DH_new hc_DH_new
@@ -119,7 +118,6 @@ struct DH {
 const DH_METHOD *DH_null_method(void);
 const DH_METHOD *DH_tfm_method(void);
 const DH_METHOD *DH_ltm_method(void);
-const DH_METHOD *DH_imath_method(void);
 
 DH *	DH_new(void);
 DH *	DH_new_method(ENGINE *);
diff --git a/source4/heimdal/lib/hcrypto/engine.c b/source4/heimdal/lib/hcrypto/engine.c
index 6e3e5e3939..15853420f6 100644
--- a/source4/heimdal/lib/hcrypto/engine.c
+++ b/source4/heimdal/lib/hcrypto/engine.c
@@ -284,26 +284,6 @@ ENGINE_load_builtin_engines(void)
 	ENGINE_finish(engine);
 #endif
 
-#ifdef USE_HCRYPTO_IMATH
-    /*
-     * imath
-     */
-
-    engine = ENGINE_new();
-    if (engine == NULL)
-	return;
-
-    ENGINE_set_id(engine, "imath");
-    ENGINE_set_name(engine,
-		    "Heimdal crypto imath engine version " PACKAGE_VERSION);
-    ENGINE_set_RSA(engine, RSA_imath_method());
-    ENGINE_set_DH(engine, DH_imath_method());
-
-    ret = add_engine(engine);
-    if (ret != 1)
-	ENGINE_finish(engine);
-#endif
-
 #ifdef HAVE_GMP
     /*
      * gmp
diff --git a/source4/heimdal/lib/hcrypto/imath/LICENSE b/source4/heimdal/lib/hcrypto/imath/LICENSE
deleted file mode 100644
index 5b0104fa1b..0000000000
--- a/source4/heimdal/lib/hcrypto/imath/LICENSE
+++ /dev/null
@@ -1,21 +0,0 @@
-IMath is Copyright © 2002-2008 Michael J. Fromberger
-You may use it subject to the following Licensing Terms:
-
-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.
diff --git a/source4/heimdal/lib/hcrypto/imath/imath.c b/source4/heimdal/lib/hcrypto/imath/imath.c
deleted file mode 100644
index 0079bafd02..0000000000
--- a/source4/heimdal/lib/hcrypto/imath/imath.c
+++ /dev/null
@@ -1,3353 +0,0 @@
-/*
-  Name:     imath.c
-  Purpose:  Arbitrary precision integer arithmetic routines.
-  Author:   M. J. Fromberger <http://spinning-yarns.org/michael/>
-  Info:     $Id: imath.c 826 2009-02-11 16:21:04Z sting $
-
-  Copyright (C) 2002-2007 Michael J. Fromberger, All Rights Reserved.
-
-  Permission is hereby granted, free of charge, to any person
-  obtaining a copy of this software and associated documentation files
-  (the "Software"), to deal in the Software without restriction,
-  including without limitation the rights to use, copy, modify, merge,
-  publish, distribute, sublicense, and/or sell copies of the Software,
-  and to permit persons to whom the Software is furnished to do so,
-  subject to the following conditions:
-
-  The above copyright notice and this permission notice shall be
-  included in all copies or substantial portions of the Software.
-
-  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
-  NONINFRINGEMENT.  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
-  BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
-  ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
-  CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-  SOFTWARE.
- */
-
-#include "imath.h"
-
-#if DEBUG
-#include <stdio.h>
-#endif
-
-#include <stdlib.h>
-#include <string.h>
-#include <ctype.h>
-
-#include <assert.h>
-
-#if DEBUG
-#define STATIC /* public */
-#else
-#define STATIC static
-#endif
-
-/* {{{ 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_result MP_MINERR = -6;
-
-const mp_sign   MP_NEG  = 1;    /* value is strictly negative */
-const mp_sign   MP_ZPOS = 0;    /* value is non-negative      */
-
-STATIC const char *s_unknown_err = "unknown result code";
-STATIC const char *s_error_msg[] = {
-  "error code 0",
-  "boolean true",
-  "out of memory",
-  "argument out of range",
-  "result undefined",
-  "output truncated",
-  "invalid argument",
-  NULL
-};
-
-/* }}} */
-
-/* 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).
-
-   The use of this table eliminates a dependency upon linkage against
-   the standard math libraries.
- */
-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,                                         /* 36          */
-};
-
-/* }}} */
-/* {{{ 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)
-
-#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)
-
-/* Select min/max.  Do not provide expressions for which multiple
-   evaluation would be problematic, e.g. x++ */
-#define MIN(A, B) ((B)<(A)?(B):(A))
-#define MAX(A, B) ((B)>(A)?(B):(A))
-
-/* Exchange lvalues A and B of type T, e.g.
-   SWAP(int, x, y) where x and y are variables of type int. */
-#define SWAP(T, A, B) do{T t_=(A);A=(B);B=t_;}while(0)
-
-/* Used to set up and access simple temp stacks within functions. */
-#define TEMP(K) (temp + (K))
-#define SETUP(E, C) \
-do{if((res = (E)) != MP_OK) goto CLEANUP; ++(C);}while(0)
-
-/* Compare value to zero. */
-#define CMPZ(Z) \
-(((Z)->used==1&&(Z)->digits[0]==0)?0:((Z)->sign==MP_NEG)?-1:1)
-
-/* Multiply X by Y into Z, ignoring signs.  Requires that Z have
-   enough storage preallocated to hold the result. */
-#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)
-
-/* Square X into Z.  Requires that Z have enough storage to hold the
-   result. */
-#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 configuration settings */
-
-/* Default number of digits allocated to a new mp_int */
-#if IMATH_TEST
-mp_size default_precision = MP_DEFAULT_PREC;
-#else
-STATIC const mp_size default_precision = MP_DEFAULT_PREC;
-#endif
-
-/* Minimum number of digits to invoke recursive multiply */
-#if IMATH_TEST
-mp_size multiply_threshold = MP_MULT_THRESH;
-#else
-STATIC const mp_size multiply_threshold = MP_MULT_THRESH;
-#endif
-
-/* }}} */
-
-/* Allocate a buffer of (at least) num digits, or return
-   NULL if that couldn't be done.  */
-STATIC mp_digit *s_alloc(mp_size num);
-
-/* Release a buffer of digits allocated by s_alloc(). */
-STATIC void s_free(void *ptr);
-
-/* Insure that z has at least min digits allocated, resizing if
-   necessary.  Returns true if successful, false if out of memory. */
-STATIC int  s_pad(mp_int z, mp_size min);
-
-/* Fill in a "fake" mp_int on the stack with a given value */
-STATIC void      s_fake(mp_int z, mp_small value, mp_digit vbuf[]);
-
-/* 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(mp_small v, mp_digit t[]);
-
-/* Compare magnitudes of a and b, returns <0, 0, >0 */
-STATIC int       s_ucmp(mp_int a, mp_int b);
-
-/* Compare magnitudes of a and v, returns <0, 0, >0 */
-STATIC int       s_vcmp(mp_int a, mp_small v);
-
-/* Unsigned magnitude addition; assumes dc is big enough.
-   Carry out is returned (no memory allocated). */
-STATIC mp_digit  s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc, 
-		        mp_size size_a, mp_size size_b);
-
-/* Unsigned magnitude subtraction.  Assumes dc is big enough. */
-STATIC void      s_usub(mp_digit *da, mp_digit *db, mp_digit *dc,
-		        mp_size size_a, mp_size size_b);
-
-/* Unsigned recursive multiplication.  Assumes dc is big enough. */
-STATIC int       s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc,
-			mp_size size_a, mp_size size_b);
-
-/* Unsigned magnitude multiplication.  Assumes dc is big enough. */
-STATIC void      s_umul(mp_digit *da, mp_digit *db, mp_digit *dc,
-			mp_size size_a, mp_size size_b);
-
-/* Unsigned recursive squaring.  Assumes dc is big enough. */
-STATIC int       s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a);
-
-/* Unsigned magnitude squaring.  Assumes dc is big enough. */
-STATIC void      s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a);
-
-/* Single digit addition.  Assumes a is big enough. */
-STATIC void      s_dadd(mp_int a, mp_digit b);
-
-/* Single digit multiplication.  Assumes a is big enough. */
-STATIC void      s_dmul(mp_int a, mp_digit b);
-
-/* Single digit multiplication on buffers; assumes dc is big enough. */
-STATIC void      s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc,
-			 mp_size size_a);
-
-/* Single digit division.  Replaces a with the quotient, 
-   returns the remainder.  */
-STATIC mp_digit  s_ddiv(mp_int a, mp_digit b);
-
-/* Quick division by a power of 2, replaces z (no allocation) */
-STATIC void      s_qdiv(mp_int z, mp_size p2);
-
-/* Quick remainder by a power of 2, replaces z (no allocation) */
-STATIC void      s_qmod(mp_int z, mp_size p2);
-
-/* Quick multiplication by a power of 2, replaces z. 
-   Allocates if necessary; returns false in case this fails. */
-STATIC int       s_qmul(mp_int z, mp_size p2);
-
-/* Quick subtraction from a power of 2, replaces z.
-   Allocates if necessary; returns false in case this fails. */
-STATIC int       s_qsub(mp_int z, mp_size p2);
-
-/* Return maximum k such that 2^k divides z. */
-STATIC int       s_dp2k(mp_int z);
-
-/* Return k >= 0 such that z = 2^k, or -1 if there is no such k. */
-STATIC int       s_isp2(mp_int z);
-
-/* Set z to 2^k.  May allocate; returns false in case this fails. */
-STATIC int       s_2expt(mp_int z, mp_small k);
-
-/* Normalize a and b for division, returns normalization constant */
-STATIC int       s_norm(mp_int a, mp_int b);
-
-/* Compute constant mu for Barrett reduction, given modulus m, result
-   replaces z, m is untouched. */
-STATIC mp_result s_brmu(mp_int z, mp_int m);
-
-/* Reduce a modulo m, using Barrett's algorithm. */
-STATIC int       s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2);
-
-/* Modular exponentiation, using Barrett reduction */
-STATIC mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c);
-
-/* Unsigned magnitude division.  Assumes |a| > |b|.  Allocates
-   temporaries; overwrites a with quotient, b with remainder. */
-STATIC mp_result s_udiv(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 DEBUG
-/* 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
-
-/* {{{ mp_int_init(z) */
-
-mp_result mp_int_init(mp_int z)
-{
-  if(z == NULL)
-    return MP_BADARG;
-
-  z->single = 0;
-  z->digits = &(z->single);
-  z->alloc  = 1;
-  z->used   = 1;
-  z->sign   = MP_ZPOS;
-
-  return MP_OK;
-}
-
-/* }}} */
-
-/* {{{ mp_int_alloc() */
-
-mp_int    mp_int_alloc(void)
-{
-  mp_int out = malloc(sizeof(mpz_t));
-  
-  if(out != NULL) 
-    mp_int_init(out);
-
-  return out;
-}
-
-/* }}} */
-
-/* {{{ mp_int_init_size(z, prec) */
-
-mp_result mp_int_init_size(mp_int z, mp_size prec)
-{
-  CHECK(z != NULL);
-
-  if(prec == 0)
-    prec = default_precision;
-  else if(prec == 1) 
-    return mp_int_init(z);
-  else 
-    prec = (mp_size) ROUND_PREC(prec);
-  
-  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;
-
-  CHECK(z != NULL && old != NULL);
-
-  uold = MP_USED(old);
-  if(uold == 1) {
-    mp_int_init(z);
-  }
-  else {
-    mp_size 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, mp_small value)
-{
-  mpz_t     vtmp;
-  mp_digit  vbuf[MP_VALUE_DIGITS(value)];
-
-  s_fake(&vtmp, value, vbuf);
-  return mp_int_init_copy(z, &vtmp);
-}
-
-/* }}} */
-
-/* {{{ mp_int_set_value(z, value) */
-
-mp_result  mp_int_set_value(mp_int z, mp_small value)
-{
-  mpz_t    vtmp;
-  mp_digit vbuf[MP_VALUE_DIGITS(value)];
-
-  s_fake(&vtmp, value, vbuf);
-  return mp_int_copy(&vtmp, z);
-}
-
-/* }}} */
-
-/* {{{ mp_int_clear(z) */
-
-void      mp_int_clear(mp_int z)
-{
-  if(z == NULL)
-    return;
-
-  if(MP_DIGITS(z) != NULL) {
-    if((void *) MP_DIGITS(z) != (void *) z)
-      s_free(MP_DIGITS(z));
-
-    MP_DIGITS(z) = NULL;
-  }
-}
-
-/* }}} */
-
-/* {{{ mp_int_free(z) */
-
-void      mp_int_free(mp_int z)
-{
-  NRCHECK(z != NULL);
-
-  mp_int_clear(z);
-  free(z); /* note: NOT s_free() */
-}
-
-/* }}} */
-
-/* {{{ 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, max;
-
-  CHECK(a != NULL && b != NULL && c != NULL);
-
-  ua = MP_USED(a); ub = MP_USED(b);
-  max = MAX(ua, ub);
-
-  if(MP_SIGN(a) == MP_SIGN(b)) {
-    /* Same sign -- add magnitudes, preserve sign of addends */
-    mp_digit carry;
-    mp_size uc;
-
-    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 comparision, sign ignored */
-
-    /* Set x to max(a, b), y to min(a, b) to simplify later code.
-       A special case yields zero for equal magnitudes.
-    */
-    if(cmp == 0) {
-      mp_int_zero(c);
-      return MP_OK;
-    }
-    else if(cmp < 0) {
-      x = b; y = a;
-    }
-    else {
-      x = a; y = b;
-    }
-
-    if(!s_pad(c, MP_USED(x)))
-      return MP_MEMORY;
-
-    /* Subtract smaller from larger */
-    s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y));
-    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, mp_small value, mp_int c)
-{
-  mpz_t     vtmp;
-  mp_digit  vbuf[MP_VALUE_DIGITS(value)];
-
-  s_fake(&vtmp, value, vbuf);
-
-  return mp_int_add(a, &vtmp, c);
-}
-
-/* }}} */
-
-/* {{{ mp_int_sub(a, b, c) */
-
-mp_result mp_int_sub(mp_int a, mp_int b, mp_int c)
-{
-  mp_size  ua, ub, max;
-
-  CHECK(a != NULL && b != NULL && c != NULL);
-
-  ua = MP_USED(a); ub = MP_USED(b);
-  max = MAX(ua, ub);
-
-  if(MP_SIGN(a) != MP_SIGN(b)) {
-    /* Different signs -- add magnitudes and keep sign of a */
-    mp_digit carry;
-    mp_size uc;
-
-    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, mp_small value, mp_int c)
-{
-  mpz_t     vtmp;
-  mp_digit  vbuf[MP_VALUE_DIGITS(value)];
-
-  s_fake(&vtmp, value, vbuf);
-
-  return mp_int_sub(a, &vtmp, c);
-}
-
-/* }}} */
-
-/* {{{ mp_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 identical to any of the inputs, we'll write
-     the results directly; otherwise, allocate a temporary space. */
-  ua = MP_USED(a); ub = MP_USED(b);
-  osize = MAX(ua, ub);
-  osize = 4 * ((osize + 1) / 2);
-
-  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)) {
-    if((void *) MP_DIGITS(c) != (void *) 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, mp_small value, mp_int c)
-{
-  mpz_t     vtmp;
-  mp_digit  vbuf[MP_VALUE_DIGITS(value)];
-
-  s_fake(&vtmp, value, vbuf);
-
-  return mp_int_mul(a, &vtmp, c);
-}
-
-/* }}} */
-
-/* {{{ mp_int_mul_pow2(a, p2, c) */
-
-mp_result mp_int_mul_pow2(mp_int a, mp_small 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) 4 * ((MP_USED(a) + 1) / 2);
-  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)) {
-    if((void *) MP_DIGITS(c) != (void *) 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) { 
-      if((res = mp_int_copy(a, q)) != MP_OK)
-	goto CLEANUP;
-      else
-	qout = q;
-    } 
-    else {
-      qout = TEMP(last);
-      SETUP(mp_int_init_copy(TEMP(last), a), last);
-    }
-
-    if(r && a != r) {
-      if((res = mp_int_copy(b, r)) != MP_OK)
-	goto CLEANUP;
-      else
-	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) {
-    mp_int_init(&tmp);
-    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, mp_small value, mp_int q, mp_small *r)
-{
-  mpz_t     vtmp, rtmp;
-  mp_digit  vbuf[MP_VALUE_DIGITS(value)];
-  mp_result res;
-
-  mp_int_init(&rtmp);
-  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, mp_small 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, mp_small 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(mp_small a, mp_small 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, mp_small 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, mp_small value, mp_int m, mp_int c)
-{
-  mpz_t    vtmp;
-  mp_digit vbuf[MP_VALUE_DIGITS(value)];
-
-  s_fake(&vtmp, value, vbuf);
-
-  return mp_int_exptmod(a, &vtmp, m, c);
-}
-
-/* }}} */
-
-/* {{{ mp_int_exptmod_bvalue(v, b, m, c) */
-
-mp_result mp_int_exptmod_bvalue(mp_small value, mp_int b,
-				mp_int m, mp_int c)
-{
-  mpz_t    vtmp;
-  mp_digit vbuf[MP_VALUE_DIGITS(value)];
-
-  s_fake(&vtmp, value, vbuf);
-
-  return mp_int_exptmod(&vtmp, b, m, c);
-}
-
-/* }}} */
-
-/* {{{ mp_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)
-    mp_int_init(TEMP(last));
-
-  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);
-
-  mp_int_init(&t);
-  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) 
-    mp_int_init(TEMP(last));
-  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_lcm(a, b, c) */
-
-mp_result mp_int_lcm(mp_int a, mp_int b, mp_int c)
-{
-  mpz_t     lcm;
-  mp_result res;
-
-  CHECK(a != NULL && b != NULL && c != NULL);
-
-  /* Since a * b = gcd(a, b) * lcm(a, b), we can compute 
-     lcm(a, b) = (a / gcd(a, b)) * b.  
-
-     This formulation insures everything works even if the input
-     variables share space.
-   */
-  if((res = mp_int_init(&lcm)) != MP_OK)
-    return res;
-  if((res = mp_int_gcd(a, b, &lcm)) != MP_OK)
-    goto CLEANUP;
-  if((res = mp_int_div(a, &lcm, &lcm, NULL)) != MP_OK)
-    goto CLEANUP;
-  if((res = mp_int_mul(&lcm, b, &lcm)) != MP_OK)
-    goto CLEANUP;
-
-  res = mp_int_copy(&lcm, c);
-
-  CLEANUP:
-    mp_int_clear(&lcm);
-
-  return res;
-}
-
-/* }}} */
-
-/* {{{ mp_int_divisible_value(a, v) */
-
-int       mp_int_divisible_value(mp_int a, mp_small v)
-{
-  mp_small rem = 0;
-
-  if(mp_int_div_value(a, v, NULL, &rem) != MP_OK)
-    return 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_root(a, b, c) */
-
-/* Implementation of Newton's root finding method, based loosely on a
-   patch contributed by Hal Finkel <half@halssoftware.com>
-   modified by M. J. Fromberger.
- */
-mp_result mp_int_root(mp_int a, mp_small b, mp_int c)
-{
-  mp_result  res = MP_OK;
-  mpz_t      temp[5];
-  int        last = 0;
-  int        flips = 0;
-
-  CHECK(a != NULL && c != NULL && b > 0);
-
-  if(b == 1) {
-    return mp_int_copy(a, c);
-  }
-  if(MP_SIGN(a) == MP_NEG) {
-    if(b % 2 == 0)
-      return MP_UNDEF; /* root does not exist for negative a with even b */
-    else
-      flips = 1;
-  }
-
-  SETUP(mp_int_init_copy(TEMP(last), a), last);
-  SETUP(mp_int_init_copy(TEMP(last), a), last);
-  SETUP(mp_int_init(TEMP(last)), last);
-  SETUP(mp_int_init(TEMP(last)), last);
-  SETUP(mp_int_init(TEMP(last)), last);
-
-  (void) mp_int_abs(TEMP(0), TEMP(0));
-  (void) mp_int_abs(TEMP(1), TEMP(1));
-
-  for(;;) {
-    if((res = mp_int_expt(TEMP(1), b, TEMP(2))) != MP_OK)
-      goto CLEANUP;
-
-    if(mp_int_compare_unsigned(TEMP(2), TEMP(0)) <= 0)
-      break;
-
-    if((res = mp_int_sub(TEMP(2), TEMP(0), TEMP(2))) != MP_OK)
-      goto CLEANUP;
-    if((res = mp_int_expt(TEMP(1), b - 1, TEMP(3))) != MP_OK)
-      goto CLEANUP;
-    if((res = mp_int_mul_value(TEMP(3), b, TEMP(3))) != MP_OK)
-      goto CLEANUP;
-    if((res = mp_int_div(TEMP(2), TEMP(3), TEMP(4), NULL)) != MP_OK)
-      goto CLEANUP;
-    if((res = mp_int_sub(TEMP(1), TEMP(4), TEMP(4))) != MP_OK)
-      goto CLEANUP;
-
-    if(mp_int_compare_unsigned(TEMP(1), TEMP(4)) == 0) {
-      if((res = mp_int_sub_value(TEMP(4), 1, TEMP(4))) != MP_OK)
-	goto CLEANUP;
-    }
-    if((res = mp_int_copy(TEMP(4), TEMP(1))) != MP_OK)
-      goto CLEANUP;
-  }
-  
-  if((res = mp_int_copy(TEMP(1), c)) != MP_OK)
-    goto CLEANUP;
-
-  /* If the original value of a was negative, flip the output sign. */
-  if(flips)
-    (void) mp_int_neg(c, c); /* cannot fail */
-
- 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, mp_small *out)
-{
-  mp_usmall 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, MP_SMALL_MAX) > 0) ||
-     mp_int_compare_value(z, MP_SMALL_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) ? -(mp_small)uv : (mp_small)uv;
-
-  return MP_OK;
-}
-
-/* }}} */
-
-/* {{{ mp_int_to_uint(z, *out) */
-
-mp_result mp_int_to_uint(mp_int z, mp_usmall *out)
-{
-  mp_usmall 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, UINT_MAX) <= 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 = 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, 1);
-  } 
-  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, 1);
-    }
-    t = str - 1;
-
-    /* Put digits back in correct output order */
-    while(h < t) {
-      char tc = *h;
-      *h++ = *t;
-      *t-- = tc;
-    }
-
-    mp_int_clear(&tmp);
-  }
-
-  *str = '\0';
-  if(cmp == 0)
-    return MP_OK;
-  else
-    return MP_TRUNC;
-}
-
-/* }}} */
-
-/* {{{ mp_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((int)*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;
-
-  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 = malloc(num * sizeof(mp_digit));
-
-  assert(out != NULL); /* for debugging */
-#if DEBUG > 1
-  {
-    mp_digit v = (mp_digit) 0xdeadbeef;
-    int      ix;
-
-    for(ix = 0; ix < num; ++ix)
-      out[ix] = v;
-  }
-#endif
-
-  return out;
-}
-
-/* }}} */
-
-/* {{{ s_realloc(old, osize, nsize) */
-
-STATIC mp_digit *s_realloc(mp_digit *old, mp_size osize, mp_size nsize)
-{
-#if DEBUG > 1
-  mp_digit *new = s_alloc(nsize);
-  int       ix;
-
-  for(ix = 0; ix < nsize; ++ix)
-    new[ix] = (mp_digit) 0xdeadbeef;
-
-  memcpy(new, old, osize * sizeof(mp_digit));
-#else
-  mp_digit *new = realloc(old, nsize * sizeof(mp_digit));
-
-  assert(new != NULL); /* for debugging */
-#endif
-  return new;
-}
-
-/* }}} */
-
-/* {{{ s_free(ptr) */
-
-STATIC void s_free(void *ptr)
-{
-  free(ptr);
-}
-
-/* }}} */
-
-/* {{{ 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;
-
-    if((void *)z->digits == (void *)z) {
-      if((tmp = s_alloc(nsize)) == NULL)
-        return 0;
-
-      COPY(MP_DIGITS(z), tmp, MP_USED(z));
-    }
-    else if((tmp = s_realloc(MP_DIGITS(z), MP_ALLOC(z), nsize)) == NULL)
-      return 0;
-    
-    MP_DIGITS(z) = tmp;
-    MP_ALLOC(z) = nsize;
-  }
-
-  return 1;
-}
-
-/* }}} */
-
-/* {{{ s_fake(z, value, vbuf) */
-
-STATIC void      s_fake(mp_int z, mp_small 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(mp_small v, mp_digit t[])
-{
-  mp_usmall    uv = (mp_usmall) ((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, mp_small 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, buf_size);
-    ZERO(t2, buf_size);
-    (void) s_kmul(da, db, t1, bot_size, bot_size);     /* t1 = a0 * b0 */
-    (void) s_kmul(a_top, b_top, t2, at_size, bt_size); /* t2 = a1 * b1 */
-
-    /* Subtract out t1 and t2 to get the inner product */
-    s_usub(t3, t1, t3, buf_size + 2, buf_size);
-    s_usub(t3, t2, t3, buf_size + 2, buf_size);
-
-    /* Assemble the output value */
-    COPY(t1, dc, buf_size);
-    carry = s_uadd(t3, dc + bot_size, dc + bot_size,
-		   buf_size + 1, buf_size); 
-    assert(carry == 0);
-    
-    carry = s_uadd(t2, dc + 2*bot_size, dc + 2*bot_size,
-		   buf_size, buf_size); 
-    assert(carry == 0);
-    
-    s_free(t1); /* note t2 and t3 are just internal pointers to t1 */
-  } 
-  else {
-    s_umul(da, db, dc, size_a, size_b);
-  }
-
-  return 1;
-}
-
-/* }}} */
-
-/* {{{ 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, carry;
-    mp_size    at_size = size_a - bot_size;
-    mp_size    buf_size = 2 * bot_size;
-
-    if((t1 = s_alloc(4 * buf_size)) == NULL) return 0;
-    t2 = t1 + buf_size;
-    t3 = t2 + buf_size;
-    ZERO(t1, 4 * buf_size);
-
-    (void) s_ksqr(da, t1, bot_size);    /* t1 = a0 ^ 2 */
-    (void) s_ksqr(a_top, t2, at_size);  /* t2 = a1 ^ 2 */
-
-    (void) s_kmul(da, a_top, t3, bot_size, at_size);  /* t3 = a0 * a1 */
-
-    /* Quick multiply t3 by 2, shifting left (can't overflow) */
-    {
-      int     i, top = bot_size + at_size;
-      mp_word w, save = 0;
-
-      for(i = 0; i < top; ++i) {
-	w = t3[i];
-	w = (w << 1) | save;
-	t3[i] = LOWER_HALF(w);
-	save = UPPER_HALF(w);
-      }
-      t3[i] = LOWER_HALF(save);
-    }
-
-    /* Assemble the output value */
-    COPY(t1, dc, 2 * bot_size);
-    carry = s_uadd(t3, dc + bot_size, dc + bot_size,
-		   buf_size + 1, buf_size);
-    assert(carry == 0);
-
-    carry = s_uadd(t2, dc + 2*bot_size, dc + 2*bot_size,
-		   buf_size, buf_size);
-    assert(carry == 0);
-
-    s_free(t1); /* note that t2 and t2 are internal pointers only */
-
-  } 
-  else {
-    s_usqr(da, dc, size_a);
-  }
-
-  return 1;
-}
-
-/* }}} */
-
-/* {{{ 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) */
-
-/* Compute z = 2^p2 - |z|; requires that 2^p2 >= |z|
-   The sign of the result is always zero/positive.
- */
-STATIC int       s_qsub(mp_int z, mp_size p2)
-{
-  mp_digit hi = (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, mp_small k)
-{
-  mp_size  ndig, rest;
-  mp_digit *dz;
-
-  ndig = (k + MP_DIGIT_BIT) / MP_DIGIT_BIT;
-  rest = k % MP_DIGIT_BIT;
-
-  if(!s_pad(z, ndig))
-    return 0;
-
-  dz = MP_DIGITS(z);
-  ZERO(dz, ndig);
-  *(dz + ndig - 1) = (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) (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 = 0;
-  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), 4 * umu), last);
-    ZERO(MP_DIGITS(TEMP(last - 1)), MP_ALLOC(TEMP(last - 1)));
-  }
-
-  (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 && qpos > 0) 
-	q.digits[qpos++] = 0;
-      
-      CLAMP(&r);
-    }
-    else {
-      mp_word  pfx = r.digits[r.used - 1];
-      mp_word  qdigit;
-      
-      if(r.used > 1 && pfx <= btop) {
-	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 = MP_DIGIT_MAX;
-      }
-      
-      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) */
-
-STATIC int       s_outlen(mp_int z, mp_size r)
-{
-  mp_result  bits;
-  double     raw;
-
-  assert(r >= MP_MIN_RADIX && r <= MP_MAX_RADIX);
-
-  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(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(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 DEBUG
-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 */
diff --git a/source4/heimdal/lib/hcrypto/imath/imath.h b/source4/heimdal/lib/hcrypto/imath/imath.h
deleted file mode 100644
index cb877959e9..0000000000
--- a/source4/heimdal/lib/hcrypto/imath/imath.h
+++ /dev/null
@@ -1,231 +0,0 @@
-/*
-  Name:     imath.h
-  Purpose:  Arbitrary precision integer arithmetic routines.
-  Author:   M. J. Fromberger <http://spinning-yarns.org/michael/>
-  Info:     $Id: imath.h 635 2008-01-08 18:19:40Z sting $
-
-  Copyright (C) 2002-2007 Michael J. Fromberger, All Rights Reserved.
-
-  Permission is hereby granted, free of charge, to any person
-  obtaining a copy of this software and associated documentation files
-  (the "Software"), to deal in the Software without restriction,
-  including without limitation the rights to use, copy, modify, merge,
-  publish, distribute, sublicense, and/or sell copies of the Software,
-  and to permit persons to whom the Software is furnished to do so,
-  subject to the following conditions:
-
-  The above copyright notice and this permission notice shall be
-  included in all copies or substantial portions of the Software.
-
-  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
-  NONINFRINGEMENT.  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
-  BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
-  ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
-  CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-  SOFTWARE.
- */
-
-#ifndef IMATH_H_
-#define IMATH_H_
-
-#include <limits.h>
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-typedef unsigned char      mp_sign;
-typedef unsigned int       mp_size;
-typedef int                mp_result;
-typedef long               mp_small;  /* must be a signed type */
-typedef unsigned long      mp_usmall; /* must be an unsigned type */
-#ifdef USE_LONG_LONG
-typedef unsigned int       mp_digit;
-typedef unsigned long long mp_word;
-#else
-typedef unsigned short     mp_digit;
-typedef unsigned int       mp_word;
-#endif
-
-typedef struct mpz {
-  mp_digit    single;
-  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;
-extern const mp_result MP_MINERR;
-
-#define MP_DIGIT_BIT    (sizeof(mp_digit) * CHAR_BIT)
-#define MP_WORD_BIT     (sizeof(mp_word) * CHAR_BIT)
-#define MP_SMALL_MIN    LONG_MIN
-#define MP_SMALL_MAX    LONG_MAX
-#define MP_USMALL_MIN   ULONG_MIN
-#define MP_USMALL_MAX   ULONG_MAX
-
-#ifdef USE_LONG_LONG
-#  ifndef ULONG_LONG_MAX
-#    ifdef ULLONG_MAX
-#      define ULONG_LONG_MAX   ULLONG_MAX
-#    else
-#      error "Maximum value of unsigned long long not defined!"
-#    endif
-#  endif
-#  define MP_DIGIT_MAX   (ULONG_MAX * 1ULL)
-#  define MP_WORD_MAX    ULONG_LONG_MAX
-#else
-#  define MP_DIGIT_MAX    (USHRT_MAX * 1UL)
-#  define MP_WORD_MAX     (UINT_MAX * 1UL)
-#endif
-
-#define MP_MIN_RADIX    2
-#define MP_MAX_RADIX    36
-
-/* Values with fewer than this many significant digits use the
-   standard multiplication algorithm; otherwise, a recursive algorithm
-   is used.  Choose a value to suit your platform.  
- */
-#define MP_MULT_THRESH  22
-
-#define MP_DEFAULT_PREC 8   /* default memory allocation, in digits */
-
-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_result mp_int_init(mp_int z);
-mp_int    mp_int_alloc(void);
-mp_result mp_int_init_size(mp_int z, mp_size prec);
-mp_result mp_int_init_copy(mp_int z, mp_int old);
-mp_result mp_int_init_value(mp_int z, mp_small value);
-mp_result mp_int_set_value(mp_int z, mp_small 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, mp_small 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, mp_small 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, mp_small value, mp_int c);
-mp_result mp_int_mul_pow2(mp_int a, mp_small p2, mp_int c);
-mp_result mp_int_sqr(mp_int a, mp_int c);            /* 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, mp_small value, /* q = a / value */
-			   mp_int q, mp_small *r);   /* r = a % value */
-mp_result mp_int_div_pow2(mp_int a, mp_small 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, mp_small b, mp_int c);         /* c = a^b */
-mp_result mp_int_expt_value(mp_small a, mp_small 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, mp_small value); /* a <=> v  */
-
-/* Returns true if v|a, false otherwise (including errors) */
-int       mp_int_divisible_value(mp_int a, mp_small 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, mp_small value, 
-				mp_int m, mp_int c);   /* c = a^v (mod m) */
-mp_result mp_int_exptmod_bvalue(mp_small 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_lcm(mp_int a, mp_int b, mp_int c);    /* c = lcm(a, b)   */
-
-mp_result mp_int_root(mp_int a, mp_small b, mp_int c); /* c = floor(a^{1/b}) */
-#define   mp_int_sqrt(a, c) mp_int_root(a, 2, c)       /* c = floor(sqrt(a)) */
-
-/* Convert to a small int, if representable; else MP_RANGE */
-mp_result mp_int_to_int(mp_int z, mp_small *out);
-mp_result mp_int_to_uint(mp_int z, mp_usmall *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 DEBUG
-void      s_print(char *tag, mp_int z);
-void      s_print_buf(char *tag, mp_digit *buf, mp_size num);
-#endif
-
-#ifdef __cplusplus
-}
-#endif
-#endif /* end IMATH_H_ */
diff --git a/source4/heimdal/lib/hcrypto/imath/iprime.c b/source4/heimdal/lib/hcrypto/imath/iprime.c
deleted file mode 100644
index 2bc9e7a6d1..0000000000
--- a/source4/heimdal/lib/hcrypto/imath/iprime.c
+++ /dev/null
@@ -1,189 +0,0 @@
-/*
-  Name:     iprime.c
-  Purpose:  Pseudoprimality testing routines
-  Author:   M. J. Fromberger <http://spinning-yarns.org/michael/>
-  Info:     $Id: iprime.c 635 2008-01-08 18:19:40Z sting $
-
-  Copyright (C) 2002-2008 Michael J. Fromberger, All Rights Reserved.
-
-  Permission is hereby granted, free of charge, to any person
-  obtaining a copy of this software and associated documentation files
-  (the "Software"), to deal in the Software without restriction,
-  including without limitation the rights to use, copy, modify, merge,
-  publish, distribute, sublicense, and/or sell copies of the Software,
-  and to permit persons to whom the Software is furnished to do so,
-  subject to the following conditions:
-
-  The above copyright notice and this permission notice shall be
-  included in all copies or substantial portions of the Software.
-
-  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
-  NONINFRINGEMENT.  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
-  BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
-  ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
-  CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-  SOFTWARE.
- */
-
-#include "iprime.h"
-#include <stdlib.h>
-
-static const int s_ptab[] = {
-    3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
-    47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101,
-    103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
-    157, 163, 167, 173, 179, 181, 191, 193, 197, 199,
-    211, 223, 227, 229, 233, 239, 241, 251, 257, 263,
-    269, 271, 277, 281, 283, 293, 307, 311, 313, 317,
-    331, 337, 347, 349, 353, 359, 367, 373, 379, 383,
-    389, 397, 401, 409, 419, 421, 431, 433, 439, 443,
-    449, 457, 461, 463, 467, 479, 487, 491, 499, 503,
-    509, 521, 523, 541, 547, 557, 563, 569, 571, 577,
-    587, 593, 599, 601, 607, 613, 617, 619, 631, 641,
-    643, 647, 653, 659, 661, 673, 677, 683, 691, 701,
-    709, 719, 727, 733, 739, 743, 751, 757, 761, 769,
-    773, 787, 797, 809, 811, 821, 823, 827, 829, 839,
-    853, 857, 859, 863, 877, 881, 883, 887, 907, 911,
-    919, 929, 937, 941, 947, 953, 967, 971, 977, 983,
-    991, 997
-#ifdef IMATH_LARGE_PRIME_TABLE
-    , 1009, 1013, 1019, 1021, 1031, 1033,
-    1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091,
-    1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151,
-    1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213,
-    1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277,
-    1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307,
-    1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399,
-    1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451,
-    1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493,
-    1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559,
-    1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609,
-    1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667,
-    1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733,
-    1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789,
-    1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871,
-    1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931,
-    1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997,
-    1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053,
-    2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111,
-    2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161,
-    2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243,
-    2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297,
-    2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357,
-    2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411,
-    2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473,
-    2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551,
-    2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633,
-    2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687,
-    2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729,
-    2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791,
-    2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851,
-    2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917,
-    2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999,
-    3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061,
-    3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137,
-    3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209,
-    3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271,
-    3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331,
-    3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391,
-    3407, 3413, 3433, 3449, 3457, 3461, 3463, 3467,
-    3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533,
-    3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583,
-    3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643,
-    3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709,
-    3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779,
-    3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851,
-    3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917,
-    3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989,
-    4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049,
-    4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111,
-    4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177,
-    4201, 4211, 4217, 4219, 4229, 4231, 4241, 4243,
-    4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297,
-    4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391,
-    4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457,
-    4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519,
-    4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597,
-    4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657,
-    4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729,
-    4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799,
-    4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889,
-    4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951,
-    4957, 4967, 4969, 4973, 4987, 4993, 4999
-#endif
-};
-static const int s_ptab_size = sizeof(s_ptab)/sizeof(s_ptab[0]);
-
-/* {{{ mp_int_is_prime(z) */
-
-/* Test whether z is likely to be prime:
-   MP_TRUE  means it is probably prime
-   MP_FALSE means it is definitely composite
- */
-mp_result mp_int_is_prime(mp_int z)
-{
-  int       i;
-  mp_small  rem;
-  mp_result res;
-
-  /* First check for divisibility by small primes; this eliminates a
-     large number of composite candidates quickly
-   */
-  for(i = 0; i < s_ptab_size; ++i) {
-    if((res = mp_int_div_value(z, s_ptab[i], NULL, &rem)) != MP_OK)
-      return res;
-
-    if(rem == 0)
-      return MP_FALSE;
-  }
-
-  /* Now try Fermat's test for several prime witnesses (since we now
-     know from the above that z is not a multiple of any of them)
-   */
-  {
-    mpz_t  tmp;
-
-    if((res = mp_int_init(&tmp)) != MP_OK) return res;
-
-    for(i = 0; i < 10 && i < s_ptab_size; ++i) {
-      if((res = mp_int_exptmod_bvalue(s_ptab[i], z, z, &tmp)) != MP_OK)
-	return res;
-
-      if(mp_int_compare_value(&tmp, s_ptab[i]) != 0) {
-	mp_int_clear(&tmp);
-	return MP_FALSE;
-      }
-    }
-
-    mp_int_clear(&tmp);
-  }
-
-  return MP_TRUE;
-}
-
-/* }}} */
-
-/* {{{ mp_int_find_prime(z) */
-
-/* Find the first apparent prime in ascending order from z */
-mp_result mp_int_find_prime(mp_int z)
-{
-  mp_result  res;
-
-  if(mp_int_is_even(z) && ((res = mp_int_add_value(z, 1, z)) != MP_OK))
-    return res;
-
-  while((res = mp_int_is_prime(z)) == MP_FALSE) {
-    if((res = mp_int_add_value(z, 2, z)) != MP_OK)
-      break;
-
-  }
-
-  return res;
-}
-
-/* }}} */
-
-/* Here there be dragons */
diff --git a/source4/heimdal/lib/hcrypto/imath/iprime.h b/source4/heimdal/lib/hcrypto/imath/iprime.h
deleted file mode 100644
index 6110dccb55..0000000000
--- a/source4/heimdal/lib/hcrypto/imath/iprime.h
+++ /dev/null
@@ -1,51 +0,0 @@
-/*
-  Name:     iprime.h
-  Purpose:  Pseudoprimality testing routines
-  Author:   M. J. Fromberger <http://spinning-yarns.org/michael/>
-  Info:     $Id: iprime.h 635 2008-01-08 18:19:40Z sting $
-
-  Copyright (C) 2002-2008 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.
- */
-
-#ifndef IPRIME_H_
-#define IPRIME_H_
-
-#include "imath.h"
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-/* Test whether z is likely to be prime
-   MP_YES means it is probably prime
-   MP_NO  means it is definitely composite
- */
-mp_result mp_int_is_prime(mp_int z);
-
-/* Find the first apparent prime in ascending order from z */
-mp_result mp_int_find_prime(mp_int z);
-
-#ifdef __cplusplus
-}
-#endif
-#endif /* IPRIME_H_ */
diff --git a/source4/heimdal/lib/hcrypto/rsa-imath.c b/source4/heimdal/lib/hcrypto/rsa-imath.c
deleted file mode 100644
index 23d5352700..0000000000
--- a/source4/heimdal/lib/hcrypto/rsa-imath.c
+++ /dev/null
@@ -1,688 +0,0 @@
-/*
- * Copyright (c) 2006 - 2007 Kungliga Tekniska Högskolan
- * (Royal Institute of Technology, Stockholm, Sweden).
- * 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.
- *
- * 3. Neither the name of the Institute nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE INSTITUTE 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 INSTITUTE 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.
- */
-
-#include <config.h>
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <krb5-types.h>
-#include <assert.h>
-
-#include <rsa.h>
-
-#include <roken.h>
-
-#ifdef USE_HCRYPTO_IMATH
-
-#include "imath/imath.h"
-#include "imath/iprime.h"
-
-static void
-BN2mpz(mpz_t *s, const BIGNUM *bn)
-{
-    size_t len;
-    void *p;
-
-    mp_int_init(s);
-
-    len = BN_num_bytes(bn);
-    p = malloc(len);
-    BN_bn2bin(bn, p);
-    mp_int_read_unsigned(s, p, len);
-    free(p);
-}
-
-static BIGNUM *
-mpz2BN(mpz_t *s)
-{
-    size_t size;
-    BIGNUM *bn;
-    void *p;
-
-    size = mp_int_unsigned_len(s);
-    p = malloc(size);
-    if (p == NULL && size != 0)
-	return NULL;
-    mp_int_to_unsigned(s, p, size);
-
-    bn = BN_bin2bn(p, size, NULL);
-    free(p);
-    return bn;
-}
-
-static int random_num(mp_int, size_t);
-
-static void
-setup_blind(mp_int n, mp_int b, mp_int bi)
-{
-    mp_int_init(b);
-    mp_int_init(bi);
-    random_num(b, mp_int_count_bits(n));
-    mp_int_mod(b, n, b);
-    mp_int_invmod(b, n, bi);
-}
-
-static void
-blind(mp_int in, mp_int b, mp_int e, mp_int n)
-{
-    mpz_t t1;
-    mp_int_init(&t1);
-    /* in' = (in * b^e) mod n */
-    mp_int_exptmod(b, e, n, &t1);
-    mp_int_mul(&t1, in, in);
-    mp_int_mod(in, n, in);
-    mp_int_clear(&t1);
-}
-
-static void
-unblind(mp_int out, mp_int bi, mp_int n)
-{
-    /* out' = (out * 1/b) mod n */
-    mp_int_mul(out, bi, out);
-    mp_int_mod(out, n, out);
-}
-
-static mp_result
-rsa_private_calculate(mp_int in, mp_int p,  mp_int q,
-		      mp_int dmp1, mp_int dmq1, mp_int iqmp,
-		      mp_int out)
-{
-    mpz_t vp, vq, u;
-    mp_int_init(&vp); mp_int_init(&vq); mp_int_init(&u);
-
-    /* vq = c ^ (d mod (q - 1)) mod q */
-    /* vp = c ^ (d mod (p - 1)) mod p */
-    mp_int_mod(in, p, &u);
-    mp_int_exptmod(&u, dmp1, p, &vp);
-    mp_int_mod(in, q, &u);
-    mp_int_exptmod(&u, dmq1, q, &vq);
-
-    /* C2 = 1/q mod p  (iqmp) */
-    /* u = (vp - vq)C2 mod p. */
-    mp_int_sub(&vp, &vq, &u);
-    if (mp_int_compare_zero(&u) < 0)
-	mp_int_add(&u, p, &u);
-    mp_int_mul(&u, iqmp, &u);
-    mp_int_mod(&u, p, &u);
-
-    /* c ^ d mod n = vq + u q */
-    mp_int_mul(&u, q, &u);
-    mp_int_add(&u, &vq, out);
-
-    mp_int_clear(&vp);
-    mp_int_clear(&vq);
-    mp_int_clear(&u);
-
-    return MP_OK;
-}
-
-/*
- *
- */
-
-static int
-imath_rsa_public_encrypt(int flen, const unsigned char* from,
-			unsigned char* to, RSA* rsa, int padding)
-{
-    unsigned char *p, *p0;
-    mp_result res;
-    size_t size, padlen;
-    mpz_t enc, dec, n, e;
-
-    if (padding != RSA_PKCS1_PADDING)
-	return -1;
-
-    size = RSA_size(rsa);
-
-    if (size < RSA_PKCS1_PADDING_SIZE || size - RSA_PKCS1_PADDING_SIZE < flen)
-	return -2;
-
-    BN2mpz(&n, rsa->n);
-    BN2mpz(&e, rsa->e);
-
-    p = p0 = malloc(size - 1);
-    if (p0 == NULL) {
-	mp_int_clear(&e);
-	mp_int_clear(&n);
-	return -3;
-    }
-
-    padlen = size - flen - 3;
-
-    *p++ = 2;
-    if (RAND_bytes(p, padlen) != 1) {
-	mp_int_clear(&e);
-	mp_int_clear(&n);
-	free(p0);
-	return -4;
-    }
-    while(padlen) {
-	if (*p == 0)
-	    *p = 1;
-	padlen--;
-	p++;
-    }
-    *p++ = 0;
-    memcpy(p, from, flen);
-    p += flen;
-    assert((p - p0) == size - 1);
-
-    mp_int_init(&enc);
-    mp_int_init(&dec);
-    mp_int_read_unsigned(&dec, p0, size - 1);
-    free(p0);
-
-    res = mp_int_exptmod(&dec, &e, &n, &enc);
-
-    mp_int_clear(&dec);
-    mp_int_clear(&e);
-    mp_int_clear(&n);
-
-    if (res != MP_OK)
-	return -4;
-
-    {
-	size_t ssize;
-	ssize = mp_int_unsigned_len(&enc);
-	assert(size >= ssize);
-	mp_int_to_unsigned(&enc, to, ssize);
-	size = ssize;
-    }
-    mp_int_clear(&enc);
-
-    return size;
-}
-
-static int
-imath_rsa_public_decrypt(int flen, const unsigned char* from,
-			 unsigned char* to, RSA* rsa, int padding)
-{
-    unsigned char *p;
-    mp_result res;
-    size_t size;
-    mpz_t s, us, n, e;
-
-    if (padding != RSA_PKCS1_PADDING)
-	return -1;
-
-    if (flen > RSA_size(rsa))
-	return -2;
-
-    BN2mpz(&n, rsa->n);
-    BN2mpz(&e, rsa->e);
-
-#if 0
-    /* Check that the exponent is larger then 3 */
-    if (mp_int_compare_value(&e, 3) <= 0) {
-	mp_int_clear(&n);
-	mp_int_clear(&e);
-	return -3;
-    }
-#endif
-
-    mp_int_init(&s);
-    mp_int_init(&us);
-    mp_int_read_unsigned(&s, rk_UNCONST(from), flen);
-
-    if (mp_int_compare(&s, &n) >= 0) {
-	mp_int_clear(&n);
-	mp_int_clear(&e);
-	return -4;
-    }
-
-    res = mp_int_exptmod(&s, &e, &n, &us);
-
-    mp_int_clear(&s);
-    mp_int_clear(&n);
-    mp_int_clear(&e);
-
-    if (res != MP_OK)
-	return -5;
-    p = to;
-
-
-    size = mp_int_unsigned_len(&us);
-    assert(size <= RSA_size(rsa));
-    mp_int_to_unsigned(&us, p, size);
-
-    mp_int_clear(&us);
-
-    /* head zero was skipped by mp_int_to_unsigned */
-    if (*p == 0)
-	return -6;
-    if (*p != 1)
-	return -7;
-    size--; p++;
-    while (size && *p == 0xff) {
-	size--; p++;
-    }
-    if (size == 0 || *p != 0)
-	return -8;
-    size--; p++;
-
-    memmove(to, p, size);
-
-    return size;
-}
-
-static int
-imath_rsa_private_encrypt(int flen, const unsigned char* from,
-			  unsigned char* to, RSA* rsa, int padding)
-{
-    unsigned char *p, *p0;
-    mp_result res;
-    int size;
-    mpz_t in, out, n, e, b, bi;
-    int blinding = (rsa->flags & RSA_FLAG_NO_BLINDING) == 0;
-    int do_unblind = 0;
-
-    if (padding != RSA_PKCS1_PADDING)
-	return -1;
-
-    size = RSA_size(rsa);
-
-    if (size < RSA_PKCS1_PADDING_SIZE || size - RSA_PKCS1_PADDING_SIZE < flen)
-	return -2;
-
-    p0 = p = malloc(size);
-    *p++ = 0;
-    *p++ = 1;
-    memset(p, 0xff, size - flen - 3);
-    p += size - flen - 3;
-    *p++ = 0;
-    memcpy(p, from, flen);
-    p += flen;
-    assert((p - p0) == size);
-
-    BN2mpz(&n, rsa->n);
-    BN2mpz(&e, rsa->e);
-
-    mp_int_init(&in);
-    mp_int_init(&out);
-    mp_int_read_unsigned(&in, p0, size);
-    free(p0);
-
-    if(mp_int_compare_zero(&in) < 0 ||
-       mp_int_compare(&in, &n) >= 0) {
-	size = -3;
-	goto out;
-    }
-
-    if (blinding) {
-	setup_blind(&n, &b, &bi);
-	blind(&in, &b, &e, &n);
-	do_unblind = 1;
-    }
-
-    if (rsa->p && rsa->q && rsa->dmp1 && rsa->dmq1 && rsa->iqmp) {
-	mpz_t p, q, dmp1, dmq1, iqmp;
-
-	BN2mpz(&p, rsa->p);
-	BN2mpz(&q, rsa->q);
-	BN2mpz(&dmp1, rsa->dmp1);
-	BN2mpz(&dmq1, rsa->dmq1);
-	BN2mpz(&iqmp, rsa->iqmp);
-
-	res = rsa_private_calculate(&in, &p, &q, &dmp1, &dmq1, &iqmp, &out);
-
-	mp_int_clear(&p);
-	mp_int_clear(&q);
-	mp_int_clear(&dmp1);
-	mp_int_clear(&dmq1);
-	mp_int_clear(&iqmp);
-
-	if (res != MP_OK) {
-	    size = -4;
-	    goto out;
-	}
-    } else {
-	mpz_t d;
-
-	BN2mpz(&d, rsa->d);
-	res = mp_int_exptmod(&in, &d, &n, &out);
-	mp_int_clear(&d);
-	if (res != MP_OK) {
-	    size = -5;
-	    goto out;
-	}
-    }
-
-    if (do_unblind)
-	unblind(&out, &bi, &n);
-
-    if (size > 0) {
-	size_t ssize;
-	ssize = mp_int_unsigned_len(&out);
-	assert(size >= ssize);
-	mp_int_to_unsigned(&out, to, size);
-	size = ssize;
-    }
-
- out:
-    if (do_unblind) {
-	mp_int_clear(&b);
-	mp_int_clear(&bi);
-    }
-
-    mp_int_clear(&e);
-    mp_int_clear(&n);
-    mp_int_clear(&in);
-    mp_int_clear(&out);
-
-    return size;
-}
-
-static int
-imath_rsa_private_decrypt(int flen, const unsigned char* from,
-			  unsigned char* to, RSA* rsa, int padding)
-{
-    unsigned char *ptr;
-    mp_result res;
-    int size;
-    mpz_t in, out, n, e, b, bi;
-    int blinding = (rsa->flags & RSA_FLAG_NO_BLINDING) == 0;
-    int do_unblind = 0;
-
-    if (padding != RSA_PKCS1_PADDING)
-	return -1;
-
-    size = RSA_size(rsa);
-    if (flen > size)
-	return -2;
-
-    mp_int_init(&in);
-    mp_int_init(&out);
-
-    BN2mpz(&n, rsa->n);
-    BN2mpz(&e, rsa->e);
-
-    res = mp_int_read_unsigned(&in, rk_UNCONST(from), flen);
-    if (res != MP_OK) {
-	size = -1;
-	goto out;
-    }
-
-    if(mp_int_compare_zero(&in) < 0 ||
-       mp_int_compare(&in, &n) >= 0) {
-	size = -2;
-	goto out;
-    }
-
-    if (blinding) {
-	setup_blind(&n, &b, &bi);
-	blind(&in, &b, &e, &n);
-	do_unblind = 1;
-    }
-
-    if (rsa->p && rsa->q && rsa->dmp1 && rsa->dmq1 && rsa->iqmp) {
-	mpz_t p, q, dmp1, dmq1, iqmp;
-
-	BN2mpz(&p, rsa->p);
-	BN2mpz(&q, rsa->q);
-	BN2mpz(&dmp1, rsa->dmp1);
-	BN2mpz(&dmq1, rsa->dmq1);
-	BN2mpz(&iqmp, rsa->iqmp);
-
-	res = rsa_private_calculate(&in, &p, &q, &dmp1, &dmq1, &iqmp, &out);
-
-	mp_int_clear(&p);
-	mp_int_clear(&q);
-	mp_int_clear(&dmp1);
-	mp_int_clear(&dmq1);
-	mp_int_clear(&iqmp);
-
-	if (res != MP_OK) {
-	    size = -3;
-	    goto out;
-	}
-
-    } else {
-	mpz_t d;
-
-	if(mp_int_compare_zero(&in) < 0 ||
-	   mp_int_compare(&in, &n) >= 0)
-	    return MP_RANGE;
-
-	BN2mpz(&d, rsa->d);
-	res = mp_int_exptmod(&in, &d, &n, &out);
-	mp_int_clear(&d);
-	if (res != MP_OK) {
-	    size = -4;
-	    goto out;
-	}
-    }
-
-    if (do_unblind)
-	unblind(&out, &bi, &n);
-
-    ptr = to;
-    {
-	size_t ssize;
-	ssize = mp_int_unsigned_len(&out);
-	assert(size >= ssize);
-	mp_int_to_unsigned(&out, ptr, ssize);
-	size = ssize;
-    }
-
-    /* head zero was skipped by mp_int_to_unsigned */
-    if (*ptr != 2) {
-	size = -5;
-	goto out;
-    }
-    size--; ptr++;
-    while (size && *ptr != 0) {
-	size--; ptr++;
-    }
-    if (size == 0)
-	return -6;
-    size--; ptr++;
-
-    memmove(to, ptr, size);
-
- out:
-    if (do_unblind) {
-	mp_int_clear(&b);
-	mp_int_clear(&bi);
-    }
-
-    mp_int_clear(&e);
-    mp_int_clear(&n);
-    mp_int_clear(&in);
-    mp_int_clear(&out);
-
-    return size;
-}
-
-static int
-random_num(mp_int num, size_t len)
-{
-    unsigned char *p;
-    mp_result res;
-
-    len = (len + 7) / 8;
-    p = malloc(len);
-    if (p == NULL)
-	return 1;
-    if (RAND_bytes(p, len) != 1) {
-	free(p);
-	return 1;
-    }
-    res = mp_int_read_unsigned(num, p, len);
-    free(p);
-    if (res != MP_OK)
-	return 1;
-    return 0;
-}
-
-#define CHECK(f, v) if ((f) != (v)) { goto out; }
-
-static int
-imath_rsa_generate_key(RSA *rsa, int bits, BIGNUM *e, BN_GENCB *cb)
-{
-    mpz_t el, p, q, n, d, dmp1, dmq1, iqmp, t1, t2, t3;
-    int counter, ret;
-
-    if (bits < 789)
-	return -1;
-
-    ret = -1;
-
-    mp_int_init(&el);
-    mp_int_init(&p);
-    mp_int_init(&q);
-    mp_int_init(&n);
-    mp_int_init(&d);
-    mp_int_init(&dmp1);
-    mp_int_init(&dmq1);
-    mp_int_init(&iqmp);
-    mp_int_init(&t1);
-    mp_int_init(&t2);
-    mp_int_init(&t3);
-
-    BN2mpz(&el, e);
-
-    /* generate p and q so that p != q and bits(pq) ~ bits */
-    counter = 0;
-    do {
-	BN_GENCB_call(cb, 2, counter++);
-	CHECK(random_num(&p, bits / 2 + 1), 0);
-	CHECK(mp_int_find_prime(&p), MP_TRUE);
-
-	CHECK(mp_int_sub_value(&p, 1, &t1), MP_OK);
-	CHECK(mp_int_gcd(&t1, &el, &t2), MP_OK);
-    } while(mp_int_compare_value(&t2, 1) != 0);
-
-    BN_GENCB_call(cb, 3, 0);
-
-    counter = 0;
-    do {
-	BN_GENCB_call(cb, 2, counter++);
-	CHECK(random_num(&q, bits / 2 + 1), 0);
-	CHECK(mp_int_find_prime(&q), MP_TRUE);
-
-	if (mp_int_compare(&p, &q) == 0) /* don't let p and q be the same */
-	    continue;
-
-	CHECK(mp_int_sub_value(&q, 1, &t1), MP_OK);
-	CHECK(mp_int_gcd(&t1, &el, &t2), MP_OK);
-    } while(mp_int_compare_value(&t2, 1) != 0);
-
-    /* make p > q */
-    if (mp_int_compare(&p, &q) < 0)
-	mp_int_swap(&p, &q);
-
-    BN_GENCB_call(cb, 3, 1);
-
-    /* calculate n,  		n = p * q */
-    CHECK(mp_int_mul(&p, &q, &n), MP_OK);
-
-    /* calculate d, 		d = 1/e mod (p - 1)(q - 1) */
-    CHECK(mp_int_sub_value(&p, 1, &t1), MP_OK);
-    CHECK(mp_int_sub_value(&q, 1, &t2), MP_OK);
-    CHECK(mp_int_mul(&t1, &t2, &t3), MP_OK);
-    CHECK(mp_int_invmod(&el, &t3, &d), MP_OK);
-
-    /* calculate dmp1		dmp1 = d mod (p-1) */
-    CHECK(mp_int_mod(&d, &t1, &dmp1), MP_OK);
-    /* calculate dmq1		dmq1 = d mod (q-1) */
-    CHECK(mp_int_mod(&d, &t2, &dmq1), MP_OK);
-    /* calculate iqmp 		iqmp = 1/q mod p */
-    CHECK(mp_int_invmod(&q, &p, &iqmp), MP_OK);
-
-    /* fill in RSA key */
-
-    rsa->e = mpz2BN(&el);
-    rsa->p = mpz2BN(&p);
-    rsa->q = mpz2BN(&q);
-    rsa->n = mpz2BN(&n);
-    rsa->d = mpz2BN(&d);
-    rsa->dmp1 = mpz2BN(&dmp1);
-    rsa->dmq1 = mpz2BN(&dmq1);
-    rsa->iqmp = mpz2BN(&iqmp);
-
-    ret = 1;
-out:
-    mp_int_clear(&el);
-    mp_int_clear(&p);
-    mp_int_clear(&q);
-    mp_int_clear(&n);
-    mp_int_clear(&d);
-    mp_int_clear(&dmp1);
-    mp_int_clear(&dmq1);
-    mp_int_clear(&iqmp);
-    mp_int_clear(&t1);
-    mp_int_clear(&t2);
-    mp_int_clear(&t3);
-
-    return ret;
-}
-
-static int
-imath_rsa_init(RSA *rsa)
-{
-    return 1;
-}
-
-static int
-imath_rsa_finish(RSA *rsa)
-{
-    return 1;
-}
-
-const RSA_METHOD hc_rsa_imath_method = {
-    "hcrypto imath RSA",
-    imath_rsa_public_encrypt,
-    imath_rsa_public_decrypt,
-    imath_rsa_private_encrypt,
-    imath_rsa_private_decrypt,
-    NULL,
-    NULL,
-    imath_rsa_init,
-    imath_rsa_finish,
-    0,
-    NULL,
-    NULL,
-    NULL,
-    imath_rsa_generate_key
-};
-#endif
-
-const RSA_METHOD *
-RSA_imath_method(void)
-{
-#ifdef USE_HCRYPTO_IMATH
-    return &hc_rsa_imath_method;
-#else
-    return NULL;
-#endif
-}
diff --git a/source4/heimdal/lib/hcrypto/rsa.c b/source4/heimdal/lib/hcrypto/rsa.c
index a6e09fe283..235686999e 100644
--- a/source4/heimdal/lib/hcrypto/rsa.c
+++ b/source4/heimdal/lib/hcrypto/rsa.c
@@ -56,15 +56,17 @@
  * Speed for RSA in seconds
  *   no key blinding
  *   1000 iteration, 
- *   same rsa key
+ *   same rsa keys (1024 and 2048)
  *   operation performed each eteration sign, verify, encrypt, decrypt on a random bit pattern
  *
- * gmp: 	 0.733615
- * tfm: 	 2.450173
- * ltm:		 3.79 (default in hcrypto)
- * openssl:	 4.04
- * cdsa:	15.89
- * imath: 	40.62
+ * name		1024	2048	4098
+ * =================================
+ * gmp: 	 0.73	  6.60	 44.80
+ * tfm: 	 2.45	    --	    --
+ * ltm:		 3.79	 20.74	105.41	(default in hcrypto)
+ * openssl:	 4.04	 11.90	 82.59
+ * cdsa:	15.89	102.89	721.40
+ * imath: 	40.62	    --	    --
  *
  * See the library functions here: @ref hcrypto_rsa
  */
@@ -516,7 +518,6 @@ RSA_null_method(void)
 }
 
 extern const RSA_METHOD hc_rsa_gmp_method;
-extern const RSA_METHOD hc_rsa_imath_method;
 extern const RSA_METHOD hc_rsa_tfm_method;
 extern const RSA_METHOD hc_rsa_ltm_method;
 static const RSA_METHOD *default_rsa_method = &hc_rsa_ltm_method;
diff --git a/source4/heimdal/lib/hcrypto/rsa.h b/source4/heimdal/lib/hcrypto/rsa.h
index 3fd805fcf0..798852d3f1 100644
--- a/source4/heimdal/lib/hcrypto/rsa.h
+++ b/source4/heimdal/lib/hcrypto/rsa.h
@@ -40,7 +40,6 @@
 
 /* symbol renaming */
 #define RSA_null_method hc_RSA_null_method
-#define RSA_imath_method hc_RSA_imath_method
 #define RSA_ltm_method hc_RSA_ltm_method
 #define RSA_gmp_method hc_RSA_gmp_method
 #define RSA_tfm_method hc_RSA_tfm_method
@@ -136,7 +135,6 @@ struct RSA {
  */
 
 const RSA_METHOD *RSA_null_method(void);
-const RSA_METHOD *RSA_imath_method(void);
 const RSA_METHOD *RSA_gmp_method(void);
 const RSA_METHOD *RSA_tfm_method(void);
 const RSA_METHOD *RSA_ltm_method(void);
diff --git a/source4/heimdal/lib/hcrypto/sha.h b/source4/heimdal/lib/hcrypto/sha.h
index a1f5a99318..fbc1810b4a 100644
--- a/source4/heimdal/lib/hcrypto/sha.h
+++ b/source4/heimdal/lib/hcrypto/sha.h
@@ -43,6 +43,12 @@
 #define SHA256_Init hc_SHA256_Init
 #define SHA256_Update hc_SHA256_Update
 #define SHA256_Final hc_SHA256_Final
+#define SHA384_Init hc_SHA384_Init
+#define SHA384_Update hc_SHA384_Update
+#define SHA384_Final hc_SHA384_Final
+#define SHA512_Init hc_SHA512_Init
+#define SHA512_Update hc_SHA512_Update
+#define SHA512_Final hc_SHA512_Final
 
 /*
  * SHA-1
diff --git a/source4/heimdal/lib/hcrypto/sha512.c b/source4/heimdal/lib/hcrypto/sha512.c
index 18447b6c67..fb38cadb6f 100644
--- a/source4/heimdal/lib/hcrypto/sha512.c
+++ b/source4/heimdal/lib/hcrypto/sha512.c
@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2006 Kungliga Tekniska Högskolan
+ * Copyright (c) 2006, 2010 Kungliga Tekniska Högskolan
  * (Royal Institute of Technology, Stockholm, Sweden).
  * All rights reserved.
  *
@@ -56,26 +56,46 @@
 #define H m->counter[7]
 
 static const uint64_t constant_512[80] = {
-    0x428a2f98d728ae22, 0x7137449123ef65cd, 0xb5c0fbcfec4d3b2f, 0xe9b5dba58189dbbc,
-    0x3956c25bf348b538, 0x59f111f1b605d019, 0x923f82a4af194f9b, 0xab1c5ed5da6d8118,
-    0xd807aa98a3030242, 0x12835b0145706fbe, 0x243185be4ee4b28c, 0x550c7dc3d5ffb4e2,
-    0x72be5d74f27b896f, 0x80deb1fe3b1696b1, 0x9bdc06a725c71235, 0xc19bf174cf692694,
-    0xe49b69c19ef14ad2, 0xefbe4786384f25e3, 0x0fc19dc68b8cd5b5, 0x240ca1cc77ac9c65,
-    0x2de92c6f592b0275, 0x4a7484aa6ea6e483, 0x5cb0a9dcbd41fbd4, 0x76f988da831153b5,
-    0x983e5152ee66dfab, 0xa831c66d2db43210, 0xb00327c898fb213f, 0xbf597fc7beef0ee4,
-    0xc6e00bf33da88fc2, 0xd5a79147930aa725, 0x06ca6351e003826f, 0x142929670a0e6e70,
-    0x27b70a8546d22ffc, 0x2e1b21385c26c926, 0x4d2c6dfc5ac42aed, 0x53380d139d95b3df,
-    0x650a73548baf63de, 0x766a0abb3c77b2a8, 0x81c2c92e47edaee6, 0x92722c851482353b,
-    0xa2bfe8a14cf10364, 0xa81a664bbc423001, 0xc24b8b70d0f89791, 0xc76c51a30654be30,
-    0xd192e819d6ef5218, 0xd69906245565a910, 0xf40e35855771202a, 0x106aa07032bbd1b8,
-    0x19a4c116b8d2d0c8, 0x1e376c085141ab53, 0x2748774cdf8eeb99, 0x34b0bcb5e19b48a8,
-    0x391c0cb3c5c95a63, 0x4ed8aa4ae3418acb, 0x5b9cca4f7763e373, 0x682e6ff3d6b2b8a3,
-    0x748f82ee5defb2fc, 0x78a5636f43172f60, 0x84c87814a1f0ab72, 0x8cc702081a6439ec,
-    0x90befffa23631e28, 0xa4506cebde82bde9, 0xbef9a3f7b2c67915, 0xc67178f2e372532b,
-    0xca273eceea26619c, 0xd186b8c721c0c207, 0xeada7dd6cde0eb1e, 0xf57d4f7fee6ed178,
-    0x06f067aa72176fba, 0x0a637dc5a2c898a6, 0x113f9804bef90dae, 0x1b710b35131c471b,
-    0x28db77f523047d84, 0x32caab7b40c72493, 0x3c9ebe0a15c9bebc, 0x431d67c49c100d4c,
-    0x4cc5d4becb3e42b6, 0x597f299cfc657e2a, 0x5fcb6fab3ad6faec, 0x6c44198c4a475817
+    0x428a2f98d728ae22ULL, 0x7137449123ef65cdULL,
+    0xb5c0fbcfec4d3b2fULL, 0xe9b5dba58189dbbcULL,
+    0x3956c25bf348b538ULL, 0x59f111f1b605d019ULL,
+    0x923f82a4af194f9bULL, 0xab1c5ed5da6d8118ULL,
+    0xd807aa98a3030242ULL, 0x12835b0145706fbeULL,
+    0x243185be4ee4b28cULL, 0x550c7dc3d5ffb4e2ULL,
+    0x72be5d74f27b896fULL, 0x80deb1fe3b1696b1ULL,
+    0x9bdc06a725c71235ULL, 0xc19bf174cf692694ULL,
+    0xe49b69c19ef14ad2ULL, 0xefbe4786384f25e3ULL,
+    0x0fc19dc68b8cd5b5ULL, 0x240ca1cc77ac9c65ULL,
+    0x2de92c6f592b0275ULL, 0x4a7484aa6ea6e483ULL,
+    0x5cb0a9dcbd41fbd4ULL, 0x76f988da831153b5ULL,
+    0x983e5152ee66dfabULL, 0xa831c66d2db43210ULL,
+    0xb00327c898fb213fULL, 0xbf597fc7beef0ee4ULL,
+    0xc6e00bf33da88fc2ULL, 0xd5a79147930aa725ULL,
+    0x06ca6351e003826fULL, 0x142929670a0e6e70ULL,
+    0x27b70a8546d22ffcULL, 0x2e1b21385c26c926ULL,
+    0x4d2c6dfc5ac42aedULL, 0x53380d139d95b3dfULL,
+    0x650a73548baf63deULL, 0x766a0abb3c77b2a8ULL,
+    0x81c2c92e47edaee6ULL, 0x92722c851482353bULL,
+    0xa2bfe8a14cf10364ULL, 0xa81a664bbc423001ULL,
+    0xc24b8b70d0f89791ULL, 0xc76c51a30654be30ULL,
+    0xd192e819d6ef5218ULL, 0xd69906245565a910ULL,
+    0xf40e35855771202aULL, 0x106aa07032bbd1b8ULL,
+    0x19a4c116b8d2d0c8ULL, 0x1e376c085141ab53ULL,
+    0x2748774cdf8eeb99ULL, 0x34b0bcb5e19b48a8ULL,
+    0x391c0cb3c5c95a63ULL, 0x4ed8aa4ae3418acbULL,
+    0x5b9cca4f7763e373ULL, 0x682e6ff3d6b2b8a3ULL,
+    0x748f82ee5defb2fcULL, 0x78a5636f43172f60ULL,
+    0x84c87814a1f0ab72ULL, 0x8cc702081a6439ecULL,
+    0x90befffa23631e28ULL, 0xa4506cebde82bde9ULL,
+    0xbef9a3f7b2c67915ULL, 0xc67178f2e372532bULL,
+    0xca273eceea26619cULL, 0xd186b8c721c0c207ULL,
+    0xeada7dd6cde0eb1eULL, 0xf57d4f7fee6ed178ULL,
+    0x06f067aa72176fbaULL, 0x0a637dc5a2c898a6ULL,
+    0x113f9804bef90daeULL, 0x1b710b35131c471bULL,
+    0x28db77f523047d84ULL, 0x32caab7b40c72493ULL,
+    0x3c9ebe0a15c9bebcULL, 0x431d67c49c100d4cULL,
+    0x4cc5d4becb3e42b6ULL, 0x597f299cfc657e2aULL,
+    0x5fcb6fab3ad6faecULL, 0x6c44198c4a475817ULL
 };
 
 void
@@ -83,14 +103,14 @@ SHA512_Init (SHA512_CTX *m)
 {
     m->sz[0] = 0;
     m->sz[1] = 0;
-    A = 0x6a09e667f3bcc908;
-    B = 0xbb67ae8584caa73b;
-    C = 0x3c6ef372fe94f82b;
-    D = 0xa54ff53a5f1d36f1;
-    E = 0x510e527fade682d1;
-    F = 0x9b05688c2b3e6c1f;
-    G = 0x1f83d9abfb41bd6b;
-    H = 0x5be0cd19137e2179;
+    A = 0x6a09e667f3bcc908ULL;
+    B = 0xbb67ae8584caa73bULL;
+    C = 0x3c6ef372fe94f82bULL;
+    D = 0xa54ff53a5f1d36f1ULL;
+    E = 0x510e527fade682d1ULL;
+    F = 0x9b05688c2b3e6c1fULL;
+    G = 0x1f83d9abfb41bd6bULL;
+    H = 0x5be0cd19137e2179ULL;
 }
 
 static void
@@ -152,17 +172,17 @@ swap_uint64_t (uint64_t t)
     uint64_t temp;
 
     temp   = cshift64(t, 32);
-    temp = ((temp & 0xff00ff00ff00ff00) >> 8) |
-           ((temp & 0x00ff00ff00ff00ff) << 8);
-    return ((temp & 0xffff0000ffff0000) >> 16) |
-           ((temp & 0x0000ffff0000ffff) << 16);
+    temp = ((temp & 0xff00ff00ff00ff00ULL) >> 8) |
+           ((temp & 0x00ff00ff00ff00ffULL) << 8);
+    return ((temp & 0xffff0000ffff0000ULL) >> 16) |
+           ((temp & 0x0000ffff0000ffffULL) << 16);
 }
-#endif
 
 struct x64{
-    uint64_t a:64;
-    uint64_t b:64;
+    uint64_t a;
+    uint64_t b;
 };
+#endif
 
 void
 SHA512_Update (SHA512_CTX *m, const void *v, size_t len)
@@ -248,14 +268,14 @@ SHA384_Init(SHA384_CTX *m)
 {
     m->sz[0] = 0;
     m->sz[1] = 0;
-    A = 0xcbbb9d5dc1059ed8;
-    B = 0x629a292a367cd507;
-    C = 0x9159015a3070dd17;
-    D = 0x152fecd8f70e5939;
-    E = 0x67332667ffc00b31;
-    F = 0x8eb44a8768581511;
-    G = 0xdb0c2e0d64f98fa7;
-    H = 0x47b5481dbefa4fa4;
+    A = 0xcbbb9d5dc1059ed8ULL;
+    B = 0x629a292a367cd507ULL;
+    C = 0x9159015a3070dd17ULL;
+    D = 0x152fecd8f70e5939ULL;
+    E = 0x67332667ffc00b31ULL;
+    F = 0x8eb44a8768581511ULL;
+    G = 0xdb0c2e0d64f98fa7ULL;
+    H = 0x47b5481dbefa4fa4ULL;
 }
 
 void