r20640: Commit part 2/2

Update Heimdal to match current lorikeet-heimdal.  This includes
integrated PAC hooks, so Samba doesn't have to handle this any more.

This also brings in the PKINIT code, hence so many new files.

Andrew Bartlett
(This used to be commit 351f7040f7bb73b9a60b22b564686f7c2f98a729)

15 files changed, 6480 insertions, 2 deletions

diff --git a/source4/heimdal/lib/des/bn.c b/source4/heimdal/lib/des/bn.c
new file mode 100644
index 0000000000..c4230b6abc
--- /dev/null
+++ b/source4/heimdal/lib/des/bn.c
@@ -0,0 +1,445 @@
+/*
+ * 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. 
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+RCSID("$Id: bn.c,v 1.9 2006/10/14 09:21:09 lha Exp $");
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <limits.h>
+
+#include <krb5-types.h>
+#include <rfc2459_asn1.h> /* XXX */
+#include <der.h>
+
+#include <bn.h>
+#include <rand.h>
+#include <hex.h>
+
+BIGNUM *
+BN_new(void)
+{
+    heim_integer *hi;
+    hi = calloc(1, sizeof(*hi));
+    return (BIGNUM *)hi;
+}
+
+void
+BN_free(BIGNUM *bn)
+{
+    BN_clear(bn);
+    free(bn);
+}
+
+void
+BN_clear(BIGNUM *bn)
+{
+    heim_integer *hi = (heim_integer *)bn;
+    if (hi->data) {
+	memset(hi->data, 0, hi->length);
+	free(hi->data);
+    }
+    memset(hi, 0, sizeof(*hi));
+}
+
+void
+BN_clear_free(BIGNUM *bn)
+{
+    BN_free(bn);
+}
+
+BIGNUM *
+BN_dup(const BIGNUM *bn)
+{
+    BIGNUM *b = BN_new();
+    if (der_copy_heim_integer((const heim_integer *)bn, (heim_integer *)b)) {
+	BN_free(b);
+	return NULL;
+    }
+    return b;
+}
+
+/*
+ * If the caller really want to know the number of bits used, subtract
+ * one from the length, multiply by 8, and then lookup in the table
+ * how many bits the hightest byte uses.
+ */
+int
+BN_num_bits(const BIGNUM *bn)
+{
+    static unsigned char num2bits[256] = {
+	0,1,2,2,3,3,3,3,4,4,4,4,4,4,4,4,  5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
+	6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,  6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
+	7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,  7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
+	7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,  7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
+	8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,  8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,
+	8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,  8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,
+	8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,  8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,
+	8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,  8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,
+    };
+    const heim_integer *i = (const void *)bn;
+    if (i->length == 0)
+	return 0;
+    return (i->length - 1) * 8 + num2bits[((unsigned char *)i->data)[0]];
+}
+
+int
+BN_num_bytes(const BIGNUM *bn)
+{
+    return ((const heim_integer *)bn)->length;
+}
+
+/*
+ * Ignore negative flag.
+ */
+
+BIGNUM *
+BN_bin2bn(const void *s, int len, BIGNUM *bn)
+{
+    heim_integer *hi = (void *)bn;
+
+    if (len < 0)
+	return NULL;
+
+    if (hi == NULL) {
+	hi = (heim_integer *)BN_new();
+	if (hi == NULL)
+	    return NULL;
+    }
+    if (hi->data)
+	BN_clear((BIGNUM *)hi);
+    hi->negative = 0;
+    hi->data = malloc(len);
+    if (hi->data == NULL && len != 0) {
+	if (bn == NULL)
+	    BN_free((BIGNUM *)hi);
+	return NULL;
+    }
+    hi->length = len;
+    memcpy(hi->data, s, len);
+    return (BIGNUM *)hi;
+}
+
+int
+BN_bn2bin(const BIGNUM *bn, void *to)
+{
+    const heim_integer *hi = (const void *)bn;
+    memcpy(to, hi->data, hi->length);
+    return hi->length;
+}
+
+int
+BN_hex2bn(BIGNUM **bnp, const char *in)
+{
+    int negative;
+    ssize_t ret;
+    size_t len;
+    void *data;
+
+    len = strlen(in);
+    data = malloc(len);
+    if (data == NULL)
+	return 0;
+
+    if (*in == '-') {
+	negative = 1;
+	in++;
+    } else
+	negative = 0;
+
+    ret = hex_decode(in, data, len);
+    if (ret < 0) {
+	free(data);
+	return 0;
+    }
+
+    *bnp = BN_bin2bn(data, ret, NULL);
+    free(data);
+    if (*bnp == NULL)
+	return 0;
+    BN_set_negative(*bnp, negative);
+    return 1;
+}
+
+char *
+BN_bn2hex(const BIGNUM *bn)
+{
+    ssize_t ret;
+    size_t len;
+    void *data;
+    char *str;
+
+    len = BN_num_bytes(bn);
+    data = malloc(len);
+    if (data == NULL)
+	return 0;
+
+    len = BN_bn2bin(bn, data);
+
+    ret = hex_encode(data, len, &str);
+    free(data);
+    if (ret < 0)
+	return 0;
+
+    return str;
+}
+
+int
+BN_cmp(const BIGNUM *bn1, const BIGNUM *bn2)
+{
+    return der_heim_integer_cmp((const heim_integer *)bn1, 
+				(const heim_integer *)bn2);
+}
+
+void
+BN_set_negative(BIGNUM *bn, int flag)
+{
+    ((heim_integer *)bn)->negative = (flag ? 1 : 0);
+}
+
+int
+BN_is_negative(BIGNUM *bn)
+{
+    return ((heim_integer *)bn)->negative ? 1 : 0;
+}
+
+static const unsigned char is_set[8] = { 1, 2, 4, 8, 16, 32, 64, 128 };
+
+int
+BN_is_bit_set(const BIGNUM *bn, int bit)
+{
+    heim_integer *hi = (heim_integer *)bn;
+    unsigned char *p = hi->data;
+
+    if ((bit / 8) > hi->length || hi->length == 0)
+	return 0;
+    
+    return p[hi->length - 1 - (bit / 8)] & is_set[bit % 8];
+}
+
+int
+BN_set_bit(BIGNUM *bn, int bit)
+{
+    heim_integer *hi = (heim_integer *)bn;
+    unsigned char *p;
+
+    if ((bit / 8) > hi->length || hi->length == 0) {
+	size_t len = (bit + 7) / 8;
+	void *d = realloc(hi->data, len);
+	if (d == NULL)
+	    return 0;
+	hi->data = d;
+	p = hi->data;
+	memset(&p[hi->length], 0, len);
+	hi->length = len;
+    } else
+	p = hi->data;
+
+    p[hi->length - 1 - (bit / 8)] |= is_set[bit % 8];
+    return 1;
+}
+
+int
+BN_clear_bit(BIGNUM *bn, int bit)
+{
+    heim_integer *hi = (heim_integer *)bn;
+    unsigned char *p = hi->data;
+
+    if ((bit / 8) > hi->length || hi->length == 0)
+	return 0;
+
+    p[hi->length - 1 - (bit / 8)] &= (unsigned char)(~(is_set[bit % 8]));
+
+    return 1;
+}
+
+int
+BN_set_word(BIGNUM *bn, unsigned long num)
+{
+    unsigned char p[sizeof(num)];
+    unsigned long num2;
+    int i, len;
+
+    for (num2 = num, i = 0; num2 > 0; i++)
+	num2 = num2 >> 8;
+
+    len = i - 1;
+    for (; i > 0; i--) {
+	p[i - 1] = (num & 0xff);
+	num = num >> 8;
+    }
+
+    bn = BN_bin2bn(p, len + 1, bn);
+    return bn != NULL;
+}
+
+unsigned long
+BN_get_word(const BIGNUM *bn)
+{
+    heim_integer *hi = (heim_integer *)bn;
+    unsigned long num = 0;
+    int i;
+
+    if (hi->negative || hi->length > sizeof(num))
+	return ULONG_MAX;
+
+    for (i = 0; i < hi->length; i++)
+	num = ((unsigned char *)hi->data)[i] | (num << 8);
+    return num;
+}
+
+int
+BN_rand(BIGNUM *bn, int bits, int top, int bottom)
+{
+    size_t len = (bits + 7) / 8;
+    heim_integer *i = (heim_integer *)bn;
+
+    BN_clear(bn);
+
+    i->negative = 0;
+    i->data = malloc(len);
+    if (i->data == NULL && len != 0)
+	return 0;
+    i->length = len;
+
+    if (RAND_bytes(i->data, i->length) != 1) {
+	free(i->data);
+	i->data = NULL;
+	return 0;
+    }
+
+    {
+	size_t j = len * 8;
+	while(j > bits) {
+	    BN_clear_bit(bn, j - 1);
+	    j--;
+	}
+    }
+
+    if (top == -1) {
+	;
+    } else if (top == 0 && bits > 0) {
+	BN_set_bit(bn, bits - 1);
+    } else if (top == 1 && bits > 1) {
+	BN_set_bit(bn, bits - 1);
+	BN_set_bit(bn, bits - 2);
+    } else {
+	BN_clear(bn);
+	return 0;
+    }
+
+    if (bottom && bits > 0)
+	BN_set_bit(bn, 0);
+
+    return 1;
+}
+
+/*
+ *
+ */
+
+int
+BN_uadd(BIGNUM *res, const BIGNUM *a, const BIGNUM *b)
+{
+    const heim_integer *ai = (const heim_integer *)a;
+    const heim_integer *bi = (const heim_integer *)b;
+    const unsigned char *ap, *bp;
+    unsigned char *cp;
+    heim_integer ci;
+    int carry = 0;
+    ssize_t len;
+
+    if (ai->negative && bi->negative)
+	return 0;
+    if (ai->length < bi->length) {
+	const heim_integer *si = bi;
+	bi = ai; ai = si;
+    }
+    
+    ci.negative = 0;
+    ci.length = ai->length + 1;
+    ci.data = malloc(ci.length);
+    if (ci.data == NULL)
+	return 0;
+
+    ap = &((const unsigned char *)ai->data)[ai->length - 1];
+    bp = &((const unsigned char *)bi->data)[bi->length - 1];
+    cp = &((unsigned char *)ci.data)[ci.length - 1];
+
+    for (len = bi->length; len > 0; len--) {
+	carry = *ap + *bp + carry;
+	*cp = carry & 0xff;
+	carry = (carry & ~0xff) ? 1 : 0;
+	ap--; bp--; cp--;
+    }
+    for (len = ai->length - bi->length; len > 0; len--) {
+	carry = *ap + carry;
+	*cp = carry & 0xff;
+	carry = (carry & ~0xff) ? 1 : 0;
+	ap--; cp--;
+    }
+    if (!carry)
+	memmove(cp, cp + 1, --ci.length);
+    else
+	*cp = carry;
+    
+    BN_clear(res);
+    *((heim_integer *)res) = ci;
+
+    return 1;
+}
+
+
+/*
+ * Callback when doing slow generation of numbers, like primes.
+ */
+
+void
+BN_GENCB_set(BN_GENCB *gencb, int (*cb_2)(int, int, BN_GENCB *), void *ctx)
+{
+    gencb->ver = 2;
+    gencb->cb.cb_2 = cb_2;
+    gencb->arg = ctx;
+}
+
+int
+BN_GENCB_call(BN_GENCB *cb, int a, int b)
+{
+    if (cb == NULL || cb->cb.cb_2 == NULL)
+	return 1;
+    return cb->cb.cb_2(a, b, cb);
+}
diff --git a/source4/heimdal/lib/des/dh-imath.c b/source4/heimdal/lib/des/dh-imath.c
new file mode 100644
index 0000000000..ebf02c72be
--- /dev/null
+++ b/source4/heimdal/lib/des/dh-imath.c
@@ -0,0 +1,243 @@
+/*
+ * 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. 
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <dh.h>
+
+#include <roken.h>
+
+#include "imath/imath.h"
+
+RCSID("$Id: dh-imath.c,v 1.6 2006/10/20 06:56:57 lha Exp $");
+
+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
+};
+
+const DH_METHOD *
+DH_imath_method(void)
+{
+    return &hc_dh_imath_method;
+}
diff --git a/source4/heimdal/lib/des/dh.c b/source4/heimdal/lib/des/dh.c
new file mode 100644
index 0000000000..66d611f6d4
--- /dev/null
+++ b/source4/heimdal/lib/des/dh.c
@@ -0,0 +1,294 @@
+/*
+ * 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. 
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+RCSID("$Id: dh.c,v 1.10 2006/10/19 17:31:51 lha Exp $");
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <dh.h>
+
+#include <roken.h>
+
+/*
+ *
+ */
+
+DH *
+DH_new(void)
+{
+    return DH_new_method(NULL);
+}
+
+DH *
+DH_new_method(ENGINE *engine)
+{
+    DH *dh;
+
+    dh = calloc(1, sizeof(*dh));
+    if (dh == NULL)
+	return NULL;
+
+    dh->references = 1;
+
+    if (engine) {
+	ENGINE_up_ref(engine);
+	dh->engine = engine;
+    } else {
+	dh->engine = ENGINE_get_default_DH();
+    }
+
+    if (dh->engine) {
+	dh->meth = ENGINE_get_DH(dh->engine);
+	if (dh->meth == NULL) {
+	    ENGINE_finish(engine);
+	    free(dh);
+	    return 0;
+	}
+    }
+
+    if (dh->meth == NULL)
+	dh->meth = DH_get_default_method();
+
+    (*dh->meth->init)(dh);
+
+    return dh;
+}
+
+void
+DH_free(DH *dh)
+{
+    if (dh->references <= 0)
+	abort();
+
+    if (--dh->references > 0)
+	return;
+
+    (*dh->meth->finish)(dh);
+
+    if (dh->engine)
+	ENGINE_finish(dh->engine);
+
+#define free_if(f) if (f) { BN_free(f); }
+    free_if(dh->p);
+    free_if(dh->g);
+    free_if(dh->pub_key);
+    free_if(dh->priv_key);
+    free_if(dh->q);
+    free_if(dh->j);
+    free_if(dh->counter);
+#undef free_if
+
+    memset(dh, 0, sizeof(*dh));
+    free(dh);
+}    
+
+int
+DH_up_ref(DH *dh)
+{
+    return ++dh->references;
+}
+
+int
+DH_size(const DH *dh)
+{
+    return BN_num_bytes(dh->p);
+}
+
+int
+DH_set_ex_data(DH *dh, int idx, void *data)
+{
+    dh->ex_data.sk = data;
+    return 1;
+}
+
+void *
+DH_get_ex_data(DH *dh, int idx)
+{
+    return dh->ex_data.sk;
+}
+
+int
+DH_generate_parameters_ex(DH *dh, int prime_len, int generator, BN_GENCB *cb)
+{
+    if (dh->meth->generate_params)
+	return dh->meth->generate_params(dh, prime_len, generator, cb);
+    return 0;
+}
+
+/*
+ * Check that
+ *
+ * 	pub_key > 1    and    pub_key < p - 1
+ *
+ * to avoid small subgroups attack.
+ */
+
+int
+DH_check_pubkey(const DH *dh, const BIGNUM *pub_key, int *codes)
+{
+    BIGNUM *bn = NULL, *sum = NULL;
+    int ret = 0;
+
+    *codes = 0;
+
+    bn = BN_new();
+    if (bn == NULL)
+	goto out;
+
+    if (!BN_set_word(bn, 1))
+	goto out;
+
+    if (BN_cmp(bn, pub_key) >= 0)
+	*codes |= DH_CHECK_PUBKEY_TOO_SMALL;
+
+    sum = BN_new();
+    if (sum == NULL)
+	goto out;
+
+    BN_uadd(sum, pub_key, bn);
+
+    if (BN_cmp(sum, dh->p) >= 0)
+	*codes |= DH_CHECK_PUBKEY_TOO_LARGE;
+
+    ret = 1;
+out:
+    if (bn)
+	BN_free(bn);
+    if (sum)
+	BN_free(sum);
+
+    return ret;
+}
+
+int
+DH_generate_key(DH *dh)
+{
+    return dh->meth->generate_key(dh);
+}
+
+int
+DH_compute_key(unsigned char *shared_key,
+	       const BIGNUM *peer_pub_key, DH *dh)
+{
+    int codes;
+
+    if (!DH_check_pubkey(dh, peer_pub_key, &codes) || codes != 0)
+	return -1;
+
+    return dh->meth->compute_key(shared_key, peer_pub_key, dh);
+}
+
+int
+DH_set_method(DH *dh, const DH_METHOD *method)
+{
+    (*dh->meth->finish)(dh);
+    if (dh->engine) {
+	ENGINE_finish(dh->engine);
+	dh->engine = NULL;
+    }
+    dh->meth = method;
+    (*dh->meth->init)(dh);
+    return 1;
+}
+
+/*
+ *
+ */
+
+static int
+dh_null_generate_key(DH *dh)
+{
+    return 0;
+}
+
+static int
+dh_null_compute_key(unsigned char *shared,const BIGNUM *pub, DH *dh)
+{
+    return 0;
+}
+
+static int
+dh_null_init(DH *dh)
+{
+    return 1;
+}
+
+static int
+dh_null_finish(DH *dh)
+{
+    return 1;
+}
+
+static int
+dh_null_generate_params(DH *dh, int prime_num, int len, BN_GENCB *cb)
+{
+    return 0;
+}
+
+static const DH_METHOD dh_null_method = {
+    "hcrypto null DH",
+    dh_null_generate_key,
+    dh_null_compute_key,
+    NULL,
+    dh_null_init,
+    dh_null_finish,
+    0,
+    NULL,
+    dh_null_generate_params
+};
+
+extern const DH_METHOD hc_dh_imath_method;
+static const DH_METHOD *dh_default_method = &hc_dh_imath_method;
+
+const DH_METHOD *
+DH_null_method(void)
+{
+    return &dh_null_method;
+}
+
+void
+DH_set_default_method(const DH_METHOD *meth)
+{
+    dh_default_method = meth;
+}
+
+const DH_METHOD *
+DH_get_default_method(void)
+{
+    return dh_default_method;
+}
+
diff --git a/source4/heimdal/lib/des/dsa.c b/source4/heimdal/lib/des/dsa.c
new file mode 100644
index 0000000000..411597b1c6
--- /dev/null
+++ b/source4/heimdal/lib/des/dsa.c
@@ -0,0 +1,125 @@
+/*
+ * 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. 
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+RCSID("$Id: dsa.c,v 1.2 2006/05/07 11:31:58 lha Exp $");
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <dsa.h>
+
+#include <roken.h>
+
+/*
+ *
+ */
+
+DSA *
+DSA_new(void)
+{
+    DSA *dsa = calloc(1, sizeof(*dsa));
+    dsa->meth = rk_UNCONST(DSA_get_default_method());
+    dsa->references = 1;
+    return dsa;
+}
+
+void
+DSA_free(DSA *dsa)
+{
+    if (dsa->references <= 0)
+	abort();
+
+    if (--dsa->references > 0)
+	return;
+
+    (*dsa->meth->finish)(dsa);
+
+#define free_if(f) if (f) { BN_free(f); }
+    free_if(dsa->p);
+    free_if(dsa->q);
+    free_if(dsa->g);
+    free_if(dsa->pub_key);
+    free_if(dsa->priv_key);
+    free_if(dsa->kinv);
+    free_if(dsa->r);
+#undef free_if
+
+    memset(dsa, 0, sizeof(*dsa));
+    free(dsa);
+
+}
+
+int
+DSA_up_ref(DSA *dsa)
+{
+    return ++dsa->references;
+}
+
+/*
+ *
+ */
+
+static const DSA_METHOD dsa_null_method = {
+    "hcrypto null DSA"
+};
+
+const DSA_METHOD *
+DSA_null_method(void)
+{
+    return &dsa_null_method;
+}
+
+
+const DSA_METHOD *dsa_default_mech = &dsa_null_method;
+
+void
+DSA_set_default_method(const DSA_METHOD *mech)
+{
+    dsa_default_mech = mech;
+}
+
+const DSA_METHOD *
+DSA_get_default_method(void)
+{
+    return dsa_default_mech;
+}
+
+int
+DSA_verify(int type, const unsigned char * digest, int digest_len,
+	   const unsigned char *sig, int sig_len, DSA *dsa)
+{
+    return -1;
+}
diff --git a/source4/heimdal/lib/des/engine.c b/source4/heimdal/lib/des/engine.c
new file mode 100644
index 0000000000..b72339c362
--- /dev/null
+++ b/source4/heimdal/lib/des/engine.c
@@ -0,0 +1,345 @@
+/*
+ * 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. 
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+RCSID("$Id: engine.c,v 1.11 2006/10/19 14:23:00 lha Exp $");
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <engine.h>
+
+#ifdef HAVE_DLFCN_H
+#include <dlfcn.h>
+#ifndef RTLD_NOW
+#define RTLD_NOW 0
+#endif
+#endif
+
+struct hc_engine {
+    int references;
+    char *name;
+    char *id;
+    void (*destroy)(ENGINE *);
+    const RSA_METHOD *rsa;
+    const DH_METHOD *dh;
+    const RAND_METHOD *rand;
+};
+
+int
+ENGINE_finish(ENGINE *engine)
+{
+    if (engine->references-- <= 0)
+	abort();
+    if (engine->references > 0)
+	return 1;
+
+    if (engine->name)
+	free(engine->name);
+    if (engine->id)
+	free(engine->id);
+    if(engine->destroy)
+	(*engine->destroy)(engine);
+
+    memset(engine, 0, sizeof(engine));
+    engine->references = -1;
+
+
+    free(engine);
+    return 1;
+}
+
+int
+ENGINE_up_ref(ENGINE *engine)
+{
+    if (engine->references < 0)
+	abort();
+    engine->references++;
+    return 1;
+}
+
+int
+ENGINE_set_id(ENGINE *engine, const char *id)
+{
+    engine->id = strdup(id);
+    return (engine->id == NULL) ? 0 : 1;
+}
+
+int
+ENGINE_set_name(ENGINE *engine, const char *name)
+{
+    engine->name = strdup(name);
+    return (engine->name == NULL) ? 0 : 1;
+}
+
+int
+ENGINE_set_RSA(ENGINE *engine, const RSA_METHOD *method)
+{
+    engine->rsa = method;
+    return 1;
+}
+
+int
+ENGINE_set_DH(ENGINE *engine, const DH_METHOD *method)
+{
+    engine->dh = method;
+    return 1;
+}
+
+int
+ENGINE_set_destroy_function(ENGINE *e, void (*destroy)(ENGINE *))
+{
+    e->destroy = destroy;
+    return 1;
+}
+
+const char *
+ENGINE_get_id(const ENGINE *engine)
+{
+    return engine->id;
+}
+
+const char *
+ENGINE_get_name(const ENGINE *engine)
+{
+    return engine->name;
+}
+
+const RSA_METHOD *
+ENGINE_get_RSA(const ENGINE *engine)
+{
+    return engine->rsa;
+}
+
+const DH_METHOD *
+ENGINE_get_DH(const ENGINE *engine)
+{
+    return engine->dh;
+}
+
+const RAND_METHOD *
+ENGINE_get_RAND(const ENGINE *engine)
+{
+    return engine->rand;
+}
+
+/*
+ *
+ */
+
+#define SG_default_engine(type)			\
+static ENGINE *type##_engine;			\
+int						\
+ENGINE_set_default_##type(ENGINE *engine)	\
+{						\
+    if (type##_engine)				\
+	ENGINE_finish(type##_engine);		\
+    type##_engine = engine;			\
+    if (type##_engine)				\
+	ENGINE_up_ref(type##_engine);		\
+    return 1;					\
+}						\
+ENGINE *					\
+ENGINE_get_default_##type(void)			\
+{						\
+    if (type##_engine)				\
+	ENGINE_up_ref(type##_engine);		\
+    return type##_engine;			\
+}
+
+SG_default_engine(RSA)
+SG_default_engine(DH)
+
+#undef SG_default_engine
+
+/*
+ *
+ */
+
+static ENGINE **engines;
+static unsigned int num_engines;
+
+static int
+add_engine(ENGINE *engine)
+{
+    ENGINE **d, *dup;
+
+    dup = ENGINE_by_id(engine->id);
+    if (dup) {
+	ENGINE_finish(dup);
+	return 0;
+    }
+
+    d = realloc(engines, (num_engines + 1) * sizeof(*engines));
+    if (d == NULL)
+	return 1;
+    engines = d;
+    engines[num_engines++] = engine;
+    
+    return 1;
+}
+
+void
+ENGINE_load_builtin_engines(void)
+{
+    ENGINE *engine;
+    int ret;
+
+    engine = calloc(1, sizeof(*engine));
+    if (engine == NULL)
+	return;
+
+    ENGINE_set_id(engine, "builtin");
+    ENGINE_set_name(engine, 
+		    "Heimdal crypto builtin 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);
+}
+
+ENGINE *
+ENGINE_by_dso(const char *path, const char *id)
+{
+#ifdef HAVE_DLOPEN
+    ENGINE *engine;
+    void *handle;
+    int ret;
+
+    engine = calloc(1, sizeof(*engine));
+    if (engine == NULL)
+	return NULL;
+
+    handle = dlopen(path, RTLD_NOW);
+    if (handle == NULL) {
+	/* printf("error: %s\n", dlerror()); */
+	free(engine);
+	return NULL;
+    }
+
+    {
+	unsigned long version;
+	openssl_v_check v_check;
+
+	v_check = (openssl_v_check)dlsym(handle, "v_check");
+	if (v_check == NULL) {
+	    dlclose(handle);
+	    free(engine);
+	    return NULL;
+	}
+
+	version = (*v_check)(OPENSSL_DYNAMIC_VERSION);
+	if (version == 0) {
+	    dlclose(handle);
+	    free(engine);
+	    return NULL;
+	}	    
+    }
+
+    {
+	openssl_bind_engine bind_engine;
+
+	bind_engine = (openssl_bind_engine)dlsym(handle, "bind_engine");
+	if (bind_engine == NULL) {
+	    dlclose(handle);
+	    free(engine);
+	    return NULL;
+	}
+
+	ret = (*bind_engine)(engine, id, NULL); /* XXX fix third arg */
+	if (ret != 1) {
+	    dlclose(handle);
+	    free(engine);
+	    return NULL;
+	}	    
+    }
+
+    ENGINE_up_ref(engine);
+
+    ret = add_engine(engine);
+    if (ret != 1) {
+	dlclose(handle);
+	ENGINE_finish(engine);
+	return NULL;
+    }
+
+    return engine;
+#else
+    return NULL;
+#endif
+}
+
+ENGINE *
+ENGINE_by_id(const char *id)
+{
+    int i;
+
+    for (i = 0; i < num_engines; i++) {
+	if (strcmp(id, engines[i]->id) == 0) {
+	    ENGINE_up_ref(engines[i]);
+	    return engines[i];
+	}
+    }
+    return NULL;
+}
+
+void
+ENGINE_add_conf_module(void)
+{
+    ENGINE *engine;
+
+    /*
+     * XXX Parse configuration file instead
+     */
+
+    engine = ENGINE_by_dso("/usr/heimdal/lib/hc-modules/hc-gmp.so", NULL);
+    if (engine == NULL)
+	return;
+    {
+	const RSA_METHOD *method = ENGINE_get_RSA(engine);
+	if (method)
+	    RSA_set_default_method(method);
+    }
+    {
+	const DH_METHOD *method = ENGINE_get_DH(engine);
+	if (method)
+	    DH_set_default_method(method);
+    }
+
+}
diff --git a/source4/heimdal/lib/des/imath/LICENSE b/source4/heimdal/lib/des/imath/LICENSE
new file mode 100644
index 0000000000..cecfb11404
--- /dev/null
+++ b/source4/heimdal/lib/des/imath/LICENSE
@@ -0,0 +1,21 @@
+IMath is Copyright 2002-2006 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/des/imath/imath.c b/source4/heimdal/lib/des/imath/imath.c
new file mode 100755
index 0000000000..0a124fa13f
--- /dev/null
+++ b/source4/heimdal/lib/des/imath/imath.c
@@ -0,0 +1,3246 @@
+/*
+  Name:     imath.c
+  Purpose:  Arbitrary precision integer arithmetic routines.
+  Author:   M. J. Fromberger <http://www.dartmouth.edu/~sting/>
+  Info:     $Id: imath.c,v 1.6 2007/01/08 10:17:31 lha Exp $
+
+  Copyright (C) 2002 Michael J. Fromberger, All Rights Reserved.
+
+  Permission is hereby granted, free of charge, to any person
+  obtaining a copy of this software and associated documentation files
+  (the "Software"), to deal in the Software without restriction,
+  including without limitation the rights to use, copy, modify, merge,
+  publish, distribute, sublicense, and/or sell copies of the Software,
+  and to permit persons to whom the Software is furnished to do so,
+  subject to the following conditions:
+
+  The above copyright notice and this permission notice shall be
+  included in all copies or substantial portions of the Software.
+
+  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+  NONINFRINGEMENT.  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+  BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+  ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+  CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+  SOFTWARE.
+ */
+
+#include "imath.h"
+
+#if DEBUG
+#include <stdio.h>
+#endif
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <ctype.h>
+
+#include <assert.h>
+
+/* {{{ Constants */
+
+const mp_result MP_OK     = 0;  /* no error, all is well  */
+const mp_result MP_FALSE  = 0;  /* boolean false          */
+const mp_result MP_TRUE   = -1; /* boolean true           */
+const mp_result MP_MEMORY = -2; /* out of memory          */
+const mp_result MP_RANGE  = -3; /* argument out of range  */
+const mp_result MP_UNDEF  = -4; /* result undefined       */
+const mp_result MP_TRUNC  = -5; /* output truncated       */
+const mp_result MP_BADARG = -6; /* invalid null argument  */
+
+const mp_sign   MP_NEG  = 1;    /* value is strictly negative */
+const mp_sign   MP_ZPOS = 0;    /* value is non-negative      */
+
+static const char *s_unknown_err = "unknown result code";
+static const char *s_error_msg[] = {
+  "error code 0",
+  "boolean true",
+  "out of memory",
+  "argument out of range",
+  "result undefined",
+  "output truncated",
+  "invalid null argument",
+  NULL
+};
+
+/* }}} */
+
+/* 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, 0.191958720, 0.190551412, 0.189200360, 	/* 36 37 38 39 */
+   0.187901825, 0.186652411, 0.185449023, 0.184288833, 	/* 40 41 42 43 */
+   0.183169251, 0.182087900, 0.181042597, 0.180031327, 	/* 44 45 46 47 */
+   0.179052232, 0.178103594, 0.177183820, 0.176291434, 	/* 48 49 50 51 */
+   0.175425064, 0.174583430, 0.173765343, 0.172969690, 	/* 52 53 54 55 */
+   0.172195434, 0.171441601, 0.170707280, 0.169991616, 	/* 56 57 58 59 */
+   0.169293808, 0.168613099, 0.167948779, 0.167300179, 	/* 60 61 62 63 */
+   0.166666667
+};
+
+/* }}} */
+/* {{{ Various macros */
+
+/* Return the number of digits needed to represent a static value */
+#define MP_VALUE_DIGITS(V) \
+((sizeof(V)+(sizeof(mp_digit)-1))/sizeof(mp_digit))
+
+/* Round precision P to nearest word boundary */
+#define ROUND_PREC(P) ((mp_size)(2*(((P)+1)/2)))
+
+/* Set array P of S digits to zero */
+#define ZERO(P, S) \
+do{mp_size i__=(S)*sizeof(mp_digit);mp_digit *p__=(P);memset(p__,0,i__);}while(0)
+
+/* Copy S digits from array P to array Q */
+#define COPY(P, Q, S) \
+do{mp_size i__=(S)*sizeof(mp_digit);mp_digit *p__=(P),*q__=(Q);\
+memcpy(q__,p__,i__);}while(0)
+
+/* Reverse N elements of type T in array A */
+#define REV(T, A, N) \
+do{T *u_=(A),*v_=u_+(N)-1;while(u_<v_){T xch=*u_;*u_++=*v_;*v_--=xch;}}while(0)
+
+#if TRACEABLE_CLAMP
+#define CLAMP(Z) s_clamp(Z)
+#else
+#define CLAMP(Z) \
+do{mp_int z_=(Z);mp_size uz_=MP_USED(z_);mp_digit *dz_=MP_DIGITS(z_)+uz_-1;\
+while(uz_ > 1 && (*dz_-- == 0)) --uz_;MP_USED(z_)=uz_;}while(0)
+#endif
+
+#define MIN(A, B) ((B)<(A)?(B):(A))
+#define MAX(A, B) ((B)>(A)?(B):(A))
+#define SWAP(T, A, B) do{T t_=(A);A=(B);B=t_;}while(0)
+
+#define TEMP(K) (temp + (K))
+#define SETUP(E, C) \
+do{if((res = (E)) != MP_OK) goto CLEANUP; ++(C);}while(0)
+
+#define CMPZ(Z) \
+(((Z)->used==1&&(Z)->digits[0]==0)?0:((Z)->sign==MP_NEG)?-1:1)
+
+#define UMUL(X, Y, Z) \
+do{mp_size ua_=MP_USED(X),ub_=MP_USED(Y);mp_size o_=ua_+ub_;\
+ZERO(MP_DIGITS(Z),o_);\
+(void) s_kmul(MP_DIGITS(X),MP_DIGITS(Y),MP_DIGITS(Z),ua_,ub_);\
+MP_USED(Z)=o_;CLAMP(Z);}while(0)
+
+#define USQR(X, Z) \
+do{mp_size ua_=MP_USED(X),o_=ua_+ua_;ZERO(MP_DIGITS(Z),o_);\
+(void) s_ksqr(MP_DIGITS(X),MP_DIGITS(Z),ua_);MP_USED(Z)=o_;CLAMP(Z);}while(0)
+
+#define UPPER_HALF(W)           ((mp_word)((W) >> MP_DIGIT_BIT))
+#define LOWER_HALF(W)           ((mp_digit)(W))
+#define HIGH_BIT_SET(W)         ((W) >> (MP_WORD_BIT - 1))
+#define ADD_WILL_OVERFLOW(W, V) ((MP_WORD_MAX - (V)) < (W))
+
+/* }}} */
+/* {{{ Default 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);
+#if TRACEABLE_FREE
+static void s_free(void *ptr);
+#else
+#define s_free(P) free(P)
+#endif
+
+/* Insure that z has at least min digits allocated, resizing if
+   necessary.  Returns true if successful, false if out of memory. */
+int       s_pad(mp_int z, mp_size min);
+
+/* Normalize by removing leading zeroes (except when z = 0) */
+#if TRACEABLE_CLAMP
+static void      s_clamp(mp_int z);
+#endif
+
+/* Fill in a "fake" mp_int on the stack with a given value */
+static void      s_fake(mp_int z, int value, mp_digit vbuf[]);
+
+/* Compare two runs of digits of given length, returns <0, 0, >0 */
+static int       s_cdig(mp_digit *da, mp_digit *db, mp_size len);
+
+/* Pack the unsigned digits of v into array t */
+static int       s_vpack(int v, mp_digit t[]);
+
+/* Compare magnitudes of a and b, returns <0, 0, >0 */
+static int       s_ucmp(mp_int a, mp_int b);
+
+/* Compare magnitudes of a and v, returns <0, 0, >0 */
+static int       s_vcmp(mp_int a, int v);
+
+/* Unsigned magnitude addition; assumes dc is big enough.
+   Carry out is returned (no memory allocated). */
+static mp_digit  s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc, 
+		        mp_size size_a, mp_size size_b);
+
+/* Unsigned magnitude subtraction.  Assumes dc is big enough. */
+static void      s_usub(mp_digit *da, mp_digit *db, mp_digit *dc,
+		        mp_size size_a, mp_size size_b);
+
+/* Unsigned recursive multiplication.  Assumes dc is big enough. */
+static int       s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc,
+			mp_size size_a, mp_size size_b);
+
+/* Unsigned magnitude multiplication.  Assumes dc is big enough. */
+static void      s_umul(mp_digit *da, mp_digit *db, mp_digit *dc,
+			mp_size size_a, mp_size size_b);
+
+/* Unsigned recursive squaring.  Assumes dc is big enough. */
+static int       s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a);
+
+/* Unsigned magnitude squaring.  Assumes dc is big enough. */
+static void      s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a);
+
+/* Single digit addition.  Assumes a is big enough. */
+static void      s_dadd(mp_int a, mp_digit b);
+
+/* Single digit multiplication.  Assumes a is big enough. */
+static void      s_dmul(mp_int a, mp_digit b);
+
+/* Single digit multiplication on buffers; assumes dc is big enough. */
+static void      s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc,
+			 mp_size size_a);
+
+/* Single digit division.  Replaces a with the quotient, 
+   returns the remainder.  */
+static mp_digit  s_ddiv(mp_int a, mp_digit b);
+
+/* Quick division by a power of 2, replaces z (no allocation) */
+static void      s_qdiv(mp_int z, mp_size p2);
+
+/* Quick remainder by a power of 2, replaces z (no allocation) */
+static void      s_qmod(mp_int z, mp_size p2);
+
+/* Quick multiplication by a power of 2, replaces z. 
+   Allocates if necessary; returns false in case this fails. */
+static int       s_qmul(mp_int z, mp_size p2);
+
+/* Quick subtraction from a power of 2, replaces z.
+   Allocates if necessary; returns false in case this fails. */
+static int       s_qsub(mp_int z, mp_size p2);
+
+/* Return maximum k such that 2^k divides z. */
+static int       s_dp2k(mp_int z);
+
+/* Return k >= 0 such that z = 2^k, or -1 if there is no such k. */
+static int       s_isp2(mp_int z);
+
+/* Set z to 2^k.  May allocate; returns false in case this fails. */
+static int       s_2expt(mp_int z, int k);
+
+/* Normalize a and b for division, returns normalization constant */
+static int       s_norm(mp_int a, mp_int b);
+
+/* Compute constant mu for Barrett reduction, given modulus m, result
+   replaces z, m is untouched. */
+static mp_result s_brmu(mp_int z, mp_int m);
+
+/* Reduce a modulo m, using Barrett's algorithm. */
+static int       s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2);
+
+/* Modular exponentiation, using Barrett reduction */
+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, int 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, int 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);
+}
+
+/* }}} */
+
+/* {{{ mp_int_copy(a, c) */
+
+mp_result mp_int_copy(mp_int a, mp_int c)
+{
+  CHECK(a != NULL && c != NULL);
+
+  if(a != c) {
+    mp_size   ua = MP_USED(a);
+    mp_digit *da, *dc;
+
+    if(!s_pad(c, ua))
+      return MP_MEMORY;
+
+    da = MP_DIGITS(a); dc = MP_DIGITS(c);
+    COPY(da, dc, ua);
+
+    MP_USED(c) = ua;
+    MP_SIGN(c) = MP_SIGN(a);
+  }
+
+  return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_swap(a, c) */
+
+void      mp_int_swap(mp_int a, mp_int c)
+{
+  if(a != c) {
+    mpz_t tmp = *a;
+
+    *a = *c;
+    *c = tmp;
+  }
+}
+
+/* }}} */
+
+/* {{{ mp_int_zero(z) */
+
+void      mp_int_zero(mp_int z)
+{
+  NRCHECK(z != NULL);
+
+  z->digits[0] = 0;
+  MP_USED(z) = 1;
+  MP_SIGN(z) = MP_ZPOS;
+}
+
+/* }}} */
+
+/* {{{ mp_int_abs(a, c) */
+
+mp_result mp_int_abs(mp_int a, mp_int c)
+{
+  mp_result res;
+
+  CHECK(a != NULL && c != NULL);
+
+  if((res = mp_int_copy(a, c)) != MP_OK)
+    return res;
+
+  MP_SIGN(c) = MP_ZPOS;
+  return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_neg(a, c) */
+
+mp_result mp_int_neg(mp_int a, mp_int c)
+{
+  mp_result res;
+
+  CHECK(a != NULL && c != NULL);
+
+  if((res = mp_int_copy(a, c)) != MP_OK)
+    return res;
+
+  if(CMPZ(c) != 0)
+    MP_SIGN(c) = 1 - MP_SIGN(a);
+
+  return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_add(a, b, c) */
+
+mp_result mp_int_add(mp_int a, mp_int b, mp_int c)
+{ 
+  mp_size  ua, ub, uc, max;
+
+  CHECK(a != NULL && b != NULL && c != NULL);
+
+  ua = MP_USED(a); ub = MP_USED(b); uc = MP_USED(c);
+  max = MAX(ua, ub);
+
+  if(MP_SIGN(a) == MP_SIGN(b)) {
+    /* Same sign -- add magnitudes, preserve sign of addends */
+    mp_digit carry;
+
+    if(!s_pad(c, max))
+      return MP_MEMORY;
+
+    carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub);
+    uc = max;
+
+    if(carry) {
+      if(!s_pad(c, max + 1))
+	return MP_MEMORY;
+
+      c->digits[max] = carry;
+      ++uc;
+    }
+
+    MP_USED(c) = uc;
+    MP_SIGN(c) = MP_SIGN(a);
+
+  } 
+  else {
+    /* Different signs -- subtract magnitudes, preserve sign of greater */
+    mp_int  x, y;
+    int     cmp = s_ucmp(a, b); /* magnitude comparision, sign ignored */
+
+    /* Set x to max(a, b), y to min(a, b) to simplify later code */
+    if(cmp >= 0) {
+      x = a; y = b;
+    } 
+    else {
+      x = b; y = a; 
+    }
+
+    if(!s_pad(c, MP_USED(x)))
+      return MP_MEMORY;
+
+    /* Subtract smaller from larger */
+    s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y));
+    MP_USED(c) = MP_USED(x);
+    CLAMP(c);
+    
+    /* Give result the sign of the larger */
+    MP_SIGN(c) = MP_SIGN(x);
+  }
+
+  return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_add_value(a, value, c) */
+
+mp_result mp_int_add_value(mp_int a, int value, mp_int c)
+{
+  mpz_t     vtmp;
+  mp_digit  vbuf[MP_VALUE_DIGITS(value)];
+
+  s_fake(&vtmp, value, vbuf);
+
+  return mp_int_add(a, &vtmp, c);
+}
+
+/* }}} */
+
+/* {{{ mp_int_sub(a, b, c) */
+
+mp_result mp_int_sub(mp_int a, mp_int b, mp_int c)
+{
+  mp_size  ua, ub, uc, max;
+
+  CHECK(a != NULL && b != NULL && c != NULL);
+
+  ua = MP_USED(a); ub = MP_USED(b); uc = MP_USED(c);
+  max = MAX(ua, ub);
+
+  if(MP_SIGN(a) != MP_SIGN(b)) {
+    /* Different signs -- add magnitudes and keep sign of a */
+    mp_digit carry;
+
+    if(!s_pad(c, max))
+      return MP_MEMORY;
+
+    carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub);
+    uc = max;
+
+    if(carry) {
+      if(!s_pad(c, max + 1))
+	return MP_MEMORY;
+
+      c->digits[max] = carry;
+      ++uc;
+    }
+
+    MP_USED(c) = uc;
+    MP_SIGN(c) = MP_SIGN(a);
+
+  } 
+  else {
+    /* Same signs -- subtract magnitudes */
+    mp_int  x, y;
+    mp_sign osign;
+    int     cmp = s_ucmp(a, b);
+
+    if(!s_pad(c, max))
+      return MP_MEMORY;
+
+    if(cmp >= 0) {
+      x = a; y = b; osign = MP_ZPOS;
+    } 
+    else {
+      x = b; y = a; osign = MP_NEG;
+    }
+
+    if(MP_SIGN(a) == MP_NEG && cmp != 0)
+      osign = 1 - osign;
+
+    s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y));
+    MP_USED(c) = MP_USED(x);
+    CLAMP(c);
+
+    MP_SIGN(c) = osign;
+  }
+
+  return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_sub_value(a, value, c) */
+
+mp_result mp_int_sub_value(mp_int a, int value, mp_int c)
+{
+  mpz_t     vtmp;
+  mp_digit  vbuf[MP_VALUE_DIGITS(value)];
+
+  s_fake(&vtmp, value, vbuf);
+
+  return mp_int_sub(a, &vtmp, c);
+}
+
+/* }}} */
+
+/* {{{ mp_int_mul(a, b, c) */
+
+mp_result mp_int_mul(mp_int a, mp_int b, mp_int c)
+{ 
+  mp_digit *out;
+  mp_size   osize, ua, ub, p = 0;
+  mp_sign   osign;
+
+  CHECK(a != NULL && b != NULL && c != NULL);
+
+  /* If either input is zero, we can shortcut multiplication */
+  if(mp_int_compare_zero(a) == 0 || mp_int_compare_zero(b) == 0) {
+    mp_int_zero(c);
+    return MP_OK;
+  }
+  
+  /* Output is positive if inputs have same sign, otherwise negative */
+  osign = (MP_SIGN(a) == MP_SIGN(b)) ? MP_ZPOS : MP_NEG;
+
+  /* If the output is not equal to any of the inputs, we'll write the
+     results there directly; otherwise, allocate a temporary space. */
+  ua = MP_USED(a); ub = MP_USED(b);
+  osize = ua + ub;
+
+  if(c == a || c == b) {
+    p = ROUND_PREC(osize);
+    p = MAX(p, default_precision);
+
+    if((out = s_alloc(p)) == NULL)
+      return MP_MEMORY;
+  } 
+  else {
+    if(!s_pad(c, osize))
+      return MP_MEMORY;
+    
+    out = MP_DIGITS(c);
+  }
+  ZERO(out, osize);
+
+  if(!s_kmul(MP_DIGITS(a), MP_DIGITS(b), out, ua, ub))
+    return MP_MEMORY;
+
+  /* If we allocated a new buffer, get rid of whatever memory c was
+     already using, and fix up its fields to reflect that.
+   */
+  if(out != MP_DIGITS(c)) {
+    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, int value, mp_int c)
+{
+  mpz_t     vtmp;
+  mp_digit  vbuf[MP_VALUE_DIGITS(value)];
+
+  s_fake(&vtmp, value, vbuf);
+
+  return mp_int_mul(a, &vtmp, c);
+}
+
+/* }}} */
+
+/* {{{ mp_int_mul_pow2(a, p2, c) */
+
+mp_result mp_int_mul_pow2(mp_int a, int p2, mp_int c)
+{
+  mp_result res;
+  CHECK(a != NULL && c != NULL && p2 >= 0);
+
+  if((res = mp_int_copy(a, c)) != MP_OK)
+    return res;
+
+  if(s_qmul(c, (mp_size) p2))
+    return MP_OK;
+  else
+    return MP_MEMORY;
+}
+
+/* }}} */
+
+/* {{{ mp_int_sqr(a, c) */
+
+mp_result mp_int_sqr(mp_int a, mp_int c)
+{ 
+  mp_digit *out;
+  mp_size   osize, p = 0;
+
+  CHECK(a != NULL && c != NULL);
+
+  /* Get a temporary buffer big enough to hold the result */
+  osize = (mp_size) 2 * MP_USED(a);
+  if(a == c) {
+    p = ROUND_PREC(osize);
+    p = MAX(p, default_precision);
+
+    if((out = s_alloc(p)) == NULL)
+      return MP_MEMORY;
+  } 
+  else {
+    if(!s_pad(c, osize)) 
+      return MP_MEMORY;
+
+    out = MP_DIGITS(c);
+  }
+  ZERO(out, osize);
+
+  s_ksqr(MP_DIGITS(a), out, MP_USED(a));
+
+  /* Get rid of whatever memory c was already using, and fix up its
+     fields to reflect the new digit array it's using
+   */
+  if(out != MP_DIGITS(c)) {
+    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 && (res = mp_int_copy(a, q)) == MP_OK) {
+      qout = q;
+    } 
+    else {
+      qout = TEMP(last);
+      SETUP(mp_int_init_copy(TEMP(last), a), last);
+    }
+
+    if(r && a != r && (res = mp_int_copy(b, r)) == MP_OK) {
+      rout = r;
+    } 
+    else {
+      rout = TEMP(last);
+      SETUP(mp_int_init_copy(TEMP(last), b), last);
+    }
+
+    if((res = s_udiv(qout, rout)) != MP_OK) goto CLEANUP;
+  } 
+  else {
+    if(q && (res = mp_int_copy(a, q)) != MP_OK) goto CLEANUP;
+    if(r && (res = mp_int_copy(a, r)) != MP_OK) goto CLEANUP;
+
+    if(q) s_qdiv(q, (mp_size) lg); qout = q;
+    if(r) s_qmod(r, (mp_size) lg); rout = r;
+  }
+
+  /* Recompute signs for output */
+  if(rout) { 
+    MP_SIGN(rout) = sa;
+    if(CMPZ(rout) == 0)
+      MP_SIGN(rout) = MP_ZPOS;
+  }
+  if(qout) {
+    MP_SIGN(qout) = (sa == sb) ? MP_ZPOS : MP_NEG;
+    if(CMPZ(qout) == 0)
+      MP_SIGN(qout) = MP_ZPOS;
+  }
+
+  if(q && (res = mp_int_copy(qout, q)) != MP_OK) goto CLEANUP;
+  if(r && (res = mp_int_copy(rout, r)) != MP_OK) goto CLEANUP;
+
+ CLEANUP:
+  while(--last >= 0)
+    mp_int_clear(TEMP(last));
+
+  return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_mod(a, m, c) */
+
+mp_result mp_int_mod(mp_int a, mp_int m, mp_int c)
+{
+  mp_result res;
+  mpz_t     tmp;
+  mp_int    out;
+
+  if(m == c) {
+    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, int value, mp_int q, int *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, int p2, mp_int q, mp_int r)
+{
+  mp_result res = MP_OK;
+
+  CHECK(a != NULL && p2 >= 0 && q != r);
+
+  if(q != NULL && (res = mp_int_copy(a, q)) == MP_OK)
+    s_qdiv(q, (mp_size) p2);
+  
+  if(res == MP_OK && r != NULL && (res = mp_int_copy(a, r)) == MP_OK)
+    s_qmod(r, (mp_size) p2);
+
+  return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_expt(a, b, c) */
+
+mp_result mp_int_expt(mp_int a, int b, mp_int c)
+{
+  mpz_t     t;
+  mp_result res;
+  unsigned int v = abs(b);
+  
+  CHECK(b >= 0 && c != NULL);
+
+  if((res = mp_int_init_copy(&t, a)) != MP_OK)
+    return res;
+
+  (void) mp_int_set_value(c, 1);
+  while(v != 0) {
+    if(v & 1) {
+      if((res = mp_int_mul(c, &t, c)) != MP_OK)
+	goto CLEANUP;
+    }
+
+    v >>= 1;
+    if(v == 0) break;
+
+    if((res = mp_int_sqr(&t, &t)) != MP_OK)
+      goto CLEANUP;
+  }
+  
+ CLEANUP:
+  mp_int_clear(&t);
+  return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_expt_value(a, b, c) */
+
+mp_result mp_int_expt_value(int a, int b, mp_int c)
+{
+  mpz_t     t;
+  mp_result res;
+  unsigned int v = abs(b);
+  
+  CHECK(b >= 0 && c != NULL);
+
+  if((res = mp_int_init_value(&t, a)) != MP_OK)
+    return res;
+
+  (void) mp_int_set_value(c, 1);
+  while(v != 0) {
+    if(v & 1) {
+      if((res = mp_int_mul(c, &t, c)) != MP_OK)
+	goto CLEANUP;
+    }
+
+    v >>= 1;
+    if(v == 0) break;
+
+    if((res = mp_int_sqr(&t, &t)) != MP_OK)
+      goto CLEANUP;
+  }
+  
+ CLEANUP:
+  mp_int_clear(&t);
+  return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_compare(a, b) */
+
+int       mp_int_compare(mp_int a, mp_int b)
+{ 
+  mp_sign sa;
+
+  CHECK(a != NULL && b != NULL);
+
+  sa = MP_SIGN(a);
+  if(sa == MP_SIGN(b)) {
+    int cmp = s_ucmp(a, b);
+
+    /* If they're both zero or positive, the normal comparison
+       applies; if both negative, the sense is reversed. */
+    if(sa == MP_ZPOS) 
+      return cmp;
+    else
+      return -cmp;
+
+  } 
+  else {
+    if(sa == MP_ZPOS)
+      return 1;
+    else
+      return -1;
+  }
+}
+
+/* }}} */
+
+/* {{{ mp_int_compare_unsigned(a, b) */
+
+int       mp_int_compare_unsigned(mp_int a, mp_int b)
+{ 
+  NRCHECK(a != NULL && b != NULL);
+
+  return s_ucmp(a, b);
+}
+
+/* }}} */
+
+/* {{{ mp_int_compare_zero(z) */
+
+int       mp_int_compare_zero(mp_int z)
+{ 
+  NRCHECK(z != NULL);
+
+  if(MP_USED(z) == 1 && z->digits[0] == 0)
+    return 0;
+  else if(MP_SIGN(z) == MP_ZPOS)
+    return 1;
+  else 
+    return -1;
+}
+
+/* }}} */
+
+/* {{{ mp_int_compare_value(z, value) */
+
+int       mp_int_compare_value(mp_int z, int value)
+{
+  mp_sign vsign = (value < 0) ? MP_NEG : MP_ZPOS;
+  int     cmp;
+
+  CHECK(z != NULL);
+
+  if(vsign == MP_SIGN(z)) {
+    cmp = s_vcmp(z, value);
+
+    if(vsign == MP_ZPOS)
+      return cmp;
+    else
+      return -cmp;
+  } 
+  else {
+    if(value < 0)
+      return 1;
+    else
+      return -1;
+  }
+}
+
+/* }}} */
+
+/* {{{ mp_int_exptmod(a, b, m, c) */
+
+mp_result mp_int_exptmod(mp_int a, mp_int b, mp_int m, mp_int c)
+{ 
+  mp_result res;
+  mp_size   um;
+  mpz_t     temp[3];
+  mp_int    s;
+  int       last = 0;
+
+  CHECK(a != NULL && b != NULL && c != NULL && m != NULL);
+
+  /* Zero moduli and negative exponents are not considered. */
+  if(CMPZ(m) == 0)
+    return MP_UNDEF;
+  if(CMPZ(b) < 0)
+    return MP_RANGE;
+
+  um = MP_USED(m);
+  SETUP(mp_int_init_size(TEMP(0), 2 * um), last);
+  SETUP(mp_int_init_size(TEMP(1), 2 * um), last);
+
+  if(c == b || c == m) {
+    SETUP(mp_int_init_size(TEMP(2), 2 * um), last);
+    s = TEMP(2);
+  } 
+  else {
+    s = c;
+  }
+  
+  if((res = mp_int_mod(a, m, TEMP(0))) != MP_OK) goto CLEANUP;
+
+  if((res = s_brmu(TEMP(1), m)) != MP_OK) goto CLEANUP;
+
+  if((res = s_embar(TEMP(0), b, m, TEMP(1), s)) != MP_OK)
+    goto CLEANUP;
+
+  res = mp_int_copy(s, c);
+
+ CLEANUP:
+  while(--last >= 0)
+    mp_int_clear(TEMP(last));
+
+  return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_exptmod_evalue(a, value, m, c) */
+
+mp_result mp_int_exptmod_evalue(mp_int a, int value, mp_int m, mp_int c)
+{
+  mpz_t    vtmp;
+  mp_digit vbuf[MP_VALUE_DIGITS(value)];
+
+  s_fake(&vtmp, value, vbuf);
+
+  return mp_int_exptmod(a, &vtmp, m, c);
+}
+
+/* }}} */
+
+/* {{{ mp_int_exptmod_bvalue(v, b, m, c) */
+
+mp_result mp_int_exptmod_bvalue(int value, mp_int b,
+				mp_int m, mp_int c)
+{
+  mpz_t    vtmp;
+  mp_digit vbuf[MP_VALUE_DIGITS(value)];
+
+  s_fake(&vtmp, value, vbuf);
+
+  return mp_int_exptmod(&vtmp, b, m, c);
+}
+
+/* }}} */
+
+/* {{{ mp_int_exptmod_known(a, b, m, mu, c) */
+
+mp_result mp_int_exptmod_known(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c)
+{
+  mp_result res;
+  mp_size   um;
+  mpz_t     temp[2];
+  mp_int    s;
+  int       last = 0;
+
+  CHECK(a && b && m && c);
+
+  /* Zero moduli and negative exponents are not considered. */
+  if(CMPZ(m) == 0)
+    return MP_UNDEF;
+  if(CMPZ(b) < 0)
+    return MP_RANGE;
+
+  um = MP_USED(m);
+  SETUP(mp_int_init_size(TEMP(0), 2 * um), last);
+
+  if(c == b || c == m) {
+    SETUP(mp_int_init_size(TEMP(1), 2 * um), last);
+    s = TEMP(1);
+  } 
+  else {
+    s = c;
+  }
+  
+  if((res = mp_int_mod(a, m, TEMP(0))) != MP_OK) goto CLEANUP;
+
+  if((res = s_embar(TEMP(0), b, m, mu, s)) != MP_OK)
+    goto CLEANUP;
+
+  res = mp_int_copy(s, c);
+
+ CLEANUP:
+  while(--last >= 0)
+    mp_int_clear(TEMP(last));
+
+  return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_redux_const(m, c) */
+
+mp_result mp_int_redux_const(mp_int m, mp_int c)
+{
+  CHECK(m != NULL && c != NULL && m != c);
+
+  return s_brmu(c, m);
+}
+
+/* }}} */
+
+/* {{{ mp_int_invmod(a, m, c) */
+
+mp_result mp_int_invmod(mp_int a, mp_int m, mp_int c)
+{
+  mp_result res;
+  mp_sign   sa;
+  int       last = 0;
+  mpz_t     temp[2];
+
+  CHECK(a != NULL && m != NULL && c != NULL);
+
+  if(CMPZ(a) == 0 || CMPZ(m) <= 0)
+    return MP_RANGE;
+
+  sa = MP_SIGN(a); /* need this for the result later */
+
+  for(last = 0; last < 2; ++last)
+    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_divisible_value(a, v) */
+
+int       mp_int_divisible_value(mp_int a, int v)
+{
+  int       rem = 0;
+
+  if(mp_int_div_value(a, v, NULL, &rem) != MP_OK)
+    return 0;
+
+  return rem == 0;
+}
+
+/* }}} */
+
+/* {{{ mp_int_is_pow2(z) */
+
+int       mp_int_is_pow2(mp_int z)
+{
+  CHECK(z != NULL);
+
+  return s_isp2(z);
+}
+
+/* }}} */
+
+/* {{{ mp_int_sqrt(a, c) */
+
+mp_result mp_int_sqrt(mp_int a, mp_int c)
+{
+  mp_result  res = MP_OK;
+  mpz_t      temp[2];
+  int        last = 0;
+
+  CHECK(a != NULL && c != NULL);
+
+  /* The square root of a negative value does not exist in the integers. */
+  if(MP_SIGN(a) == MP_NEG)
+    return MP_UNDEF;
+
+  SETUP(mp_int_init_copy(TEMP(last), a), last);
+  SETUP(mp_int_init(TEMP(last)), last);
+
+  for(;;) {
+    if((res = mp_int_sqr(TEMP(0), TEMP(1))) != MP_OK)
+      goto CLEANUP;
+
+    if(mp_int_compare_unsigned(a, TEMP(1)) == 0) break;
+
+    if((res = mp_int_copy(a, TEMP(1))) != MP_OK) 
+      goto CLEANUP;
+    if((res = mp_int_div(TEMP(1), TEMP(0), TEMP(1), NULL)) != MP_OK) 
+      goto CLEANUP;
+    if((res = mp_int_add(TEMP(0), TEMP(1), TEMP(1))) != MP_OK) 
+      goto CLEANUP;
+    if((res = mp_int_div_pow2(TEMP(1), 1, TEMP(1), NULL)) != MP_OK)
+      goto CLEANUP;
+
+    if(mp_int_compare_unsigned(TEMP(0), TEMP(1)) == 0) break;
+    if((res = mp_int_sub_value(TEMP(0), 1, TEMP(0))) != MP_OK) goto CLEANUP;
+    if(mp_int_compare_unsigned(TEMP(0), TEMP(1)) == 0) break;
+
+    if((res = mp_int_copy(TEMP(1), TEMP(0))) != MP_OK) goto CLEANUP;
+  }
+  
+  res = mp_int_copy(TEMP(0), c);
+
+ CLEANUP:
+  while(--last >= 0)
+    mp_int_clear(TEMP(last));
+  
+  return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_to_int(z, out) */
+
+mp_result mp_int_to_int(mp_int z, int *out)
+{
+  unsigned int uv = 0;
+  mp_size   uz;
+  mp_digit *dz;
+  mp_sign   sz;
+
+  CHECK(z != NULL);
+
+  /* Make sure the value is representable as an int */
+  sz = MP_SIGN(z);
+  if((sz == MP_ZPOS && mp_int_compare_value(z, INT_MAX) > 0) ||
+     mp_int_compare_value(z, INT_MIN) < 0)
+    return MP_RANGE;
+     
+  uz = MP_USED(z);
+  dz = MP_DIGITS(z) + uz - 1;
+  
+  while(uz > 0) {
+    uv <<= MP_DIGIT_BIT/2;
+    uv = (uv << (MP_DIGIT_BIT/2)) | *dz--;
+    --uz;
+  }
+
+  if(out)
+    *out = (sz == MP_NEG) ? -(int)uv : (int)uv;
+
+  return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_to_string(z, radix, str, limit) */
+
+mp_result mp_int_to_string(mp_int z, mp_size radix, 
+			   char *str, int limit)
+{
+  mp_result res;
+  int       cmp = 0;
+
+  CHECK(z != NULL && str != NULL && limit >= 2);
+
+  if(radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
+    return MP_RANGE;
+
+  if(CMPZ(z) == 0) {
+    *str++ = s_val2ch(0, 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 = mp_int_unsigned_len(z);
+
+  if(res <= 0)
+    return res;
+
+  bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT;
+
+  /* If the highest-order bit falls exactly on a byte boundary, we
+     need to pad with an extra byte so that the sign will be read
+     correctly when reading it back in. */
+  if(bytes * CHAR_BIT == res)
+    ++bytes;
+
+  return bytes;
+}
+
+/* }}} */
+
+/* {{{ mp_int_to_unsigned(z, buf, limit) */
+
+mp_result mp_int_to_unsigned(mp_int z, unsigned char *buf, int limit)
+{
+  static const int NO_PADDING = 0;
+
+  CHECK(z != NULL && buf != NULL);
+
+  return s_tobin(z, buf, &limit, NO_PADDING);
+}
+
+/* }}} */
+
+/* {{{ mp_int_read_unsigned(z, buf, len) */
+
+mp_result mp_int_read_unsigned(mp_int z, unsigned char *buf, int len)
+{
+  mp_size need, i;
+  unsigned char *tmp;
+  mp_digit *dz;
+
+  CHECK(z != NULL && buf != NULL && len > 0);
+
+  /* Figure out how many digits are needed to represent this value */
+  need = ((len * CHAR_BIT) + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT;
+  if(!s_pad(z, need))
+    return MP_MEMORY;
+
+  mp_int_zero(z);
+
+  dz = MP_DIGITS(z);
+  for(tmp = buf, i = len; i > 0; --i, ++tmp) {
+    (void) s_qmul(z, CHAR_BIT);
+    *dz |= *tmp;
+  }
+
+  return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_unsigned_len(z) */
+
+mp_result mp_int_unsigned_len(mp_int z)
+{
+  mp_result  res = mp_int_count_bits(z);
+  int        bytes;
+
+  if(res <= 0)
+    return res;
+
+  bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT;
+
+  return bytes;
+}
+
+/* }}} */
+
+/* {{{ mp_error_string(res) */
+
+const char *mp_error_string(mp_result res)
+{
+  int ix;
+  if(res > 0)
+    return s_unknown_err;
+
+  res = -res;
+  for(ix = 0; ix < res && s_error_msg[ix] != NULL; ++ix)
+    ;
+
+  if(s_error_msg[ix] != NULL)
+    return s_error_msg[ix];
+  else
+    return s_unknown_err;
+}
+
+/* }}} */
+
+/*------------------------------------------------------------------------*/
+/* Private functions for internal use.  These make assumptions.           */
+
+/* {{{ s_alloc(num) */
+
+static mp_digit *s_alloc(mp_size num)
+{
+  mp_digit *out = malloc(num * sizeof(mp_digit));
+
+  assert(out != NULL); /* for debugging */
+
+  return out;
+}
+
+/* }}} */
+
+/* {{{ s_realloc(old, num) */
+
+static mp_digit *s_realloc(mp_digit *old, mp_size num)
+{
+  mp_digit *new = realloc(old, num * sizeof(mp_digit));
+
+  assert(new != NULL); /* for debugging */
+
+  return new;
+}
+
+/* }}} */
+
+/* {{{ s_free(ptr) */
+
+#if TRACEABLE_FREE
+static void s_free(void *ptr)
+{
+  free(ptr);
+}
+#endif
+
+/* }}} */
+
+/* {{{ s_pad(z, min) */
+
+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), nsize)) == NULL)
+      return 0;
+    
+    MP_DIGITS(z) = tmp;
+    MP_ALLOC(z) = nsize;
+  }
+
+  return 1;
+}
+
+/* }}} */
+
+/* {{{ s_clamp(z) */
+
+#if TRACEABLE_CLAMP
+static void     s_clamp(mp_int z)
+{
+  mp_size   uz = MP_USED(z);
+  mp_digit *zd = MP_DIGITS(z) + uz - 1;
+
+  while(uz > 1 && (*zd-- == 0))
+    --uz;
+
+  MP_USED(z) = uz;
+}
+#endif
+
+/* }}} */
+
+/* {{{ s_fake(z, value, vbuf) */
+
+static void      s_fake(mp_int z, int value, mp_digit vbuf[])
+{
+  mp_size uv = (mp_size) s_vpack(value, vbuf);
+
+  z->used = uv;
+  z->alloc = MP_VALUE_DIGITS(value);
+  z->sign = (value < 0) ? MP_NEG : MP_ZPOS;
+  z->digits = vbuf;
+}
+
+/* }}} */
+
+/* {{{ s_cdig(da, db, len) */
+
+static int      s_cdig(mp_digit *da, mp_digit *db, mp_size len)
+{
+  mp_digit *dat = da + len - 1, *dbt = db + len - 1;
+
+  for(/* */; len != 0; --len, --dat, --dbt) {
+    if(*dat > *dbt)
+      return 1;
+    else if(*dat < *dbt)
+      return -1;
+  }
+
+  return 0;
+}
+
+/* }}} */
+
+/* {{{ s_vpack(v, t[]) */
+
+static int       s_vpack(int v, mp_digit t[])
+{
+  unsigned int uv = (unsigned int)((v < 0) ? -v : v);
+  int          ndig = 0;
+  
+  if(uv == 0)
+    t[ndig++] = 0;
+  else {
+    while(uv != 0) {
+      t[ndig++] = (mp_digit) uv;
+      uv >>= MP_DIGIT_BIT/2;
+      uv >>= MP_DIGIT_BIT/2;
+    }
+  }
+
+  return ndig;
+}
+
+/* }}} */
+
+/* {{{ s_ucmp(a, b) */
+
+static int      s_ucmp(mp_int a, mp_int b)
+{
+  mp_size  ua = MP_USED(a), ub = MP_USED(b);
+  
+  if(ua > ub)
+    return 1;
+  else if(ub > ua) 
+    return -1;
+  else 
+    return s_cdig(MP_DIGITS(a), MP_DIGITS(b), ua);
+}
+
+/* }}} */
+
+/* {{{ s_vcmp(a, v) */
+
+static int      s_vcmp(mp_int a, int v)
+{
+  mp_digit     vdig[MP_VALUE_DIGITS(v)];
+  int          ndig = 0;
+  mp_size      ua = MP_USED(a);
+
+  ndig = s_vpack(v, vdig);
+
+  if(ua > ndig)
+    return 1;
+  else if(ua < ndig)
+    return -1;
+  else
+    return s_cdig(MP_DIGITS(a), vdig, ndig);
+}
+
+/* }}} */
+
+/* {{{ s_uadd(da, db, dc, size_a, size_b) */
+
+static mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc, 
+		       mp_size size_a, mp_size size_b)
+{
+  mp_size pos;
+  mp_word w = 0;
+
+  /* Insure that da is the longer of the two to simplify later code */
+  if(size_b > size_a) {
+    SWAP(mp_digit *, da, db);
+    SWAP(mp_size, size_a, size_b);
+  }
+
+  /* Add corresponding digits until the shorter number runs out */
+  for(pos = 0; pos < size_b; ++pos, ++da, ++db, ++dc) {
+    w = w + (mp_word) *da + (mp_word) *db;
+    *dc = LOWER_HALF(w);
+    w = UPPER_HALF(w);
+  }
+
+  /* Propagate carries as far as necessary */
+  for(/* */; pos < size_a; ++pos, ++da, ++dc) {
+    w = w + *da;
+
+    *dc = LOWER_HALF(w);
+    w = UPPER_HALF(w);
+  }
+
+  /* Return carry out */
+  return (mp_digit)w;
+}
+
+/* }}} */
+
+/* {{{ s_usub(da, db, dc, size_a, size_b) */
+
+static void     s_usub(mp_digit *da, mp_digit *db, mp_digit *dc,
+		       mp_size size_a, mp_size size_b)
+{
+  mp_size pos;
+  mp_word w = 0;
+
+  /* We assume that |a| >= |b| so this should definitely hold */
+  assert(size_a >= size_b);
+
+  /* Subtract corresponding digits and propagate borrow */
+  for(pos = 0; pos < size_b; ++pos, ++da, ++db, ++dc) {
+    w = ((mp_word)MP_DIGIT_MAX + 1 +  /* MP_RADIX */
+	 (mp_word)*da) - w - (mp_word)*db;
+
+    *dc = LOWER_HALF(w);
+    w = (UPPER_HALF(w) == 0);
+  }
+
+  /* Finish the subtraction for remaining upper digits of da */
+  for(/* */; pos < size_a; ++pos, ++da, ++dc) {
+    w = ((mp_word)MP_DIGIT_MAX + 1 +  /* MP_RADIX */
+	 (mp_word)*da) - w; 
+
+    *dc = LOWER_HALF(w);
+    w = (UPPER_HALF(w) == 0);
+  }
+
+  /* If there is a borrow out at the end, it violates the precondition */
+  assert(w == 0);
+}
+
+/* }}} */
+
+/* {{{ s_kmul(da, db, dc, size_a, size_b) */
+
+static int       s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc,
+			mp_size size_a, mp_size size_b)
+{
+  mp_size  bot_size;
+
+  /* Make sure b is the smaller of the two input values */
+  if(size_b > size_a) {
+    SWAP(mp_digit *, da, db);
+    SWAP(mp_size, size_a, size_b);
+  }
+
+  /* Insure that the bottom is the larger half in an odd-length split;
+     the code below relies on this being true.
+   */
+  bot_size = (size_a + 1) / 2;
+
+  /* If the values are big enough to bother with recursion, use the
+     Karatsuba algorithm to compute the product; otherwise use the
+     normal multiplication algorithm
+   */
+  if(multiply_threshold && 
+     size_a >= multiply_threshold && 
+     size_b > bot_size) {
+
+    mp_digit *t1, *t2, *t3, carry;
+
+    mp_digit *a_top = da + bot_size; 
+    mp_digit *b_top = db + bot_size;
+
+    mp_size  at_size = size_a - bot_size;
+    mp_size  bt_size = size_b - bot_size;
+    mp_size  buf_size = 2 * bot_size;
+
+    /* Do a single allocation for all three temporary buffers needed;
+       each buffer must be big enough to hold the product of two
+       bottom halves, and one buffer needs space for the completed 
+       product; twice the space is plenty.
+     */
+    if((t1 = s_alloc(4 * buf_size)) == NULL) return 0;
+    t2 = t1 + buf_size; 
+    t3 = t2 + buf_size;
+    ZERO(t1, 4 * buf_size);
+
+    /* t1 and t2 are initially used as temporaries to compute the inner product
+       (a1 + a0)(b1 + b0) = a1b1 + a1b0 + a0b1 + a0b0
+     */
+    carry = s_uadd(da, a_top, t1, bot_size, at_size); /* t1 = a1 + a0 */
+    t1[bot_size] = carry;
+
+    carry = s_uadd(db, b_top, t2, bot_size, bt_size); /* t2 = b1 + b0 */
+    t2[bot_size] = carry;
+
+    (void) s_kmul(t1, t2, t3, bot_size + 1, bot_size + 1);   /* t3 = t1 * t2 */
+
+    /* Now we'll get t1 = a0b0 and t2 = a1b1, and subtract them out so that
+       we're left with only the pieces we want:  t3 = a1b0 + a0b1
+     */
+    ZERO(t1, bot_size + 1);
+    ZERO(t2, bot_size + 1);
+    (void) s_kmul(da, db, t1, bot_size, bot_size);     /* t1 = a0 * b0 */
+    (void) s_kmul(a_top, b_top, t2, at_size, bt_size); /* t2 = a1 * b1 */
+
+    /* Subtract out t1 and t2 to get the inner product */
+    s_usub(t3, t1, t3, buf_size + 2, buf_size);
+    s_usub(t3, t2, t3, buf_size + 2, buf_size);
+
+    /* Assemble the output value */
+    COPY(t1, dc, buf_size);
+    (void) s_uadd(t3, dc + bot_size, dc + bot_size,
+		  buf_size + 1, buf_size + 1); 
+    
+    (void) s_uadd(t2, dc + 2*bot_size, dc + 2*bot_size,
+		  buf_size, buf_size); 
+    
+    s_free(t1); /* note t2 and t3 are just internal pointers to t1 */
+  } 
+  else {
+    s_umul(da, db, dc, size_a, size_b);
+  }
+
+  return 1;
+}
+
+/* }}} */
+
+/* {{{ s_umul(da, db, dc, size_a, size_b) */
+
+static void     s_umul(mp_digit *da, mp_digit *db, mp_digit *dc,
+		       mp_size size_a, mp_size size_b)
+{
+  mp_size   a, b;
+  mp_word   w;
+
+  for(a = 0; a < size_a; ++a, ++dc, ++da) {
+    mp_digit *dct = dc;
+    mp_digit *dbt = db;
+
+    if(*da == 0)
+      continue;
+
+    w = 0;
+    for(b = 0; b < size_b; ++b, ++dbt, ++dct) {
+      w = (mp_word)*da * (mp_word)*dbt + w + (mp_word)*dct;
+
+      *dct = LOWER_HALF(w);
+      w = UPPER_HALF(w);
+    }
+
+    *dct = (mp_digit)w;
+  }
+}
+
+/* }}} */
+
+/* {{{ s_ksqr(da, dc, size_a) */
+
+static int       s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a)
+{
+  if(multiply_threshold && size_a > multiply_threshold) {
+    mp_size    bot_size = (size_a + 1) / 2;
+    mp_digit  *a_top = da + bot_size;
+    mp_digit  *t1, *t2, *t3;
+    mp_size    at_size = size_a - bot_size;
+    mp_size    buf_size = 2 * bot_size;
+
+    if((t1 = s_alloc(4 * buf_size)) == NULL) return 0;
+    t2 = t1 + buf_size;
+    t3 = t2 + buf_size;
+    ZERO(t1, 4 * buf_size);
+
+    (void) s_ksqr(da, t1, bot_size);    /* t1 = a0 ^ 2 */
+    (void) s_ksqr(a_top, t2, at_size);  /* t2 = a1 ^ 2 */
+
+    (void) s_kmul(da, a_top, t3, bot_size, at_size);  /* t3 = a0 * a1 */
+
+    /* Quick multiply t3 by 2, shifting left (can't overflow) */
+    {
+      int     i, top = bot_size + at_size;
+      mp_word w, save = 0;
+
+      for(i = 0; i < top; ++i) {
+	w = t3[i];
+	w = (w << 1) | save;
+	t3[i] = LOWER_HALF(w);
+	save = UPPER_HALF(w);
+      }
+      t3[i] = LOWER_HALF(save);
+    }
+
+    /* Assemble the output value */
+    COPY(t1, dc, 2 * bot_size);
+    (void) s_uadd(t3, dc + bot_size, dc + bot_size,
+		  buf_size + 1, buf_size + 1);
+    
+    (void) s_uadd(t2, dc + 2*bot_size, dc + 2*bot_size,
+		  buf_size, buf_size);
+
+    free(t1); /* note that t2 and t2 are internal pointers only */
+
+  } 
+  else {
+    s_usqr(da, dc, size_a);
+  }
+
+  return 1;
+}
+
+/* }}} */
+
+/* {{{ s_usqr(da, dc, size_a) */
+
+static void      s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a)
+{
+  mp_size  i, j;
+  mp_word  w;
+
+  for(i = 0; i < size_a; ++i, dc += 2, ++da) {
+    mp_digit  *dct = dc, *dat = da;
+
+    if(*da == 0)
+      continue;
+
+    /* Take care of the first digit, no rollover */
+    w = (mp_word)*dat * (mp_word)*dat + (mp_word)*dct;
+    *dct = LOWER_HALF(w);
+    w = UPPER_HALF(w);
+    ++dat; ++dct;
+
+    for(j = i + 1; j < size_a; ++j, ++dat, ++dct) {
+      mp_word  t = (mp_word)*da * (mp_word)*dat;
+      mp_word  u = w + (mp_word)*dct, ov = 0;
+
+      /* Check if doubling t will overflow a word */
+      if(HIGH_BIT_SET(t))
+	ov = 1;
+
+      w = t + t;
+
+      /* Check if adding u to w will overflow a word */
+      if(ADD_WILL_OVERFLOW(w, u))
+	ov = 1;
+
+      w += u;
+
+      *dct = LOWER_HALF(w);
+      w = UPPER_HALF(w);
+      if(ov) {
+	w += MP_DIGIT_MAX; /* MP_RADIX */
+	++w;
+      }
+    }
+
+    w = w + *dct;
+    *dct = (mp_digit)w; 
+    while((w = UPPER_HALF(w)) != 0) {
+      ++dct; w = w + *dct;
+      *dct = LOWER_HALF(w);
+    }
+
+    assert(w == 0);
+  }
+}
+
+/* }}} */
+
+/* {{{ s_dadd(a, b) */
+
+static void      s_dadd(mp_int a, mp_digit b)
+{
+  mp_word   w = 0;
+  mp_digit *da = MP_DIGITS(a);
+  mp_size   ua = MP_USED(a);
+
+  w = (mp_word)*da + b;
+  *da++ = LOWER_HALF(w);
+  w = UPPER_HALF(w);
+
+  for(ua -= 1; ua > 0; --ua, ++da) {
+    w = (mp_word)*da + w;
+
+    *da = LOWER_HALF(w);
+    w = UPPER_HALF(w);
+  }
+
+  if(w) {
+    *da = (mp_digit)w;
+    MP_USED(a) += 1;
+  }
+}
+
+/* }}} */
+
+/* {{{ s_dmul(a, b) */
+
+static void      s_dmul(mp_int a, mp_digit b)
+{
+  mp_word   w = 0;
+  mp_digit *da = MP_DIGITS(a);
+  mp_size   ua = MP_USED(a);
+
+  while(ua > 0) {
+    w = (mp_word)*da * b + w;
+    *da++ = LOWER_HALF(w);
+    w = UPPER_HALF(w);
+    --ua;
+  }
+
+  if(w) {
+    *da = (mp_digit)w;
+    MP_USED(a) += 1;
+  }
+}
+
+/* }}} */
+
+/* {{{ s_dbmul(da, b, dc, size_a) */
+
+static void      s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a)
+{
+  mp_word  w = 0;
+
+  while(size_a > 0) {
+    w = (mp_word)*da++ * (mp_word)b + w;
+
+    *dc++ = LOWER_HALF(w);
+    w = UPPER_HALF(w);
+    --size_a;
+  }
+
+  if(w)
+    *dc = LOWER_HALF(w);
+}
+
+/* }}} */
+
+/* {{{ s_ddiv(da, d, dc, size_a) */
+
+static mp_digit  s_ddiv(mp_int a, mp_digit b)
+{
+  mp_word   w = 0, qdigit;
+  mp_size   ua = MP_USED(a);
+  mp_digit *da = MP_DIGITS(a) + ua - 1;
+  
+  for(/* */; ua > 0; --ua, --da) {
+    w = (w << MP_DIGIT_BIT) | *da;
+
+    if(w >= b) {
+      qdigit = w / b;
+      w = w % b;
+    } 
+    else {
+      qdigit = 0;
+    }
+      
+    *da = (mp_digit)qdigit;
+  }
+
+  CLAMP(a);
+  return (mp_digit)w;
+}
+
+/* }}} */
+
+/* {{{ s_qdiv(z, p2) */
+
+static void     s_qdiv(mp_int z, mp_size p2)
+{
+  mp_size ndig = p2 / MP_DIGIT_BIT, nbits = p2 % MP_DIGIT_BIT;
+  mp_size uz = MP_USED(z);
+
+  if(ndig) {
+    mp_size  mark;
+    mp_digit *to, *from;
+
+    if(ndig >= uz) {
+      mp_int_zero(z);
+      return;
+    }
+
+    to = MP_DIGITS(z); from = to + ndig;
+
+    for(mark = ndig; mark < uz; ++mark) 
+      *to++ = *from++;
+
+    MP_USED(z) = uz - ndig;
+  }
+
+  if(nbits) {
+    mp_digit d = 0, *dz, save;
+    mp_size  up = MP_DIGIT_BIT - nbits;
+
+    uz = MP_USED(z);
+    dz = MP_DIGITS(z) + uz - 1;
+
+    for(/* */; uz > 0; --uz, --dz) {
+      save = *dz;
+
+      *dz = (*dz >> nbits) | (d << up);
+      d = save;
+    }
+
+    CLAMP(z);
+  }
+
+  if(MP_USED(z) == 1 && z->digits[0] == 0)
+    MP_SIGN(z) = MP_ZPOS;
+}
+
+/* }}} */
+
+/* {{{ s_qmod(z, p2) */
+
+static void     s_qmod(mp_int z, mp_size p2)
+{
+  mp_size   start = p2 / MP_DIGIT_BIT + 1, rest = p2 % MP_DIGIT_BIT;
+  mp_size   uz = MP_USED(z);
+  mp_digit  mask = (1 << rest) - 1;
+
+  if(start <= uz) {
+    MP_USED(z) = start;
+    z->digits[start - 1] &= mask;
+    CLAMP(z);
+  }
+}
+
+/* }}} */
+
+/* {{{ s_qmul(z, p2) */
+
+static int      s_qmul(mp_int z, mp_size p2)
+{
+  mp_size   uz, need, rest, extra, i;
+  mp_digit *from, *to, d;
+
+  if(p2 == 0)
+    return 1;
+
+  uz = MP_USED(z); 
+  need = p2 / MP_DIGIT_BIT; rest = p2 % MP_DIGIT_BIT;
+
+  /* Figure out if we need an extra digit at the top end; this occurs
+     if the topmost `rest' bits of the high-order digit of z are not
+     zero, meaning they will be shifted off the end if not preserved */
+  extra = 0;
+  if(rest != 0) {
+    mp_digit *dz = MP_DIGITS(z) + uz - 1;
+
+    if((*dz >> (MP_DIGIT_BIT - rest)) != 0)
+      extra = 1;
+  }
+
+  if(!s_pad(z, uz + need + extra))
+    return 0;
+
+  /* If we need to shift by whole digits, do that in one pass, then
+     to back and shift by partial digits.
+   */
+  if(need > 0) {
+    from = MP_DIGITS(z) + uz - 1;
+    to = from + need;
+
+    for(i = 0; i < uz; ++i)
+      *to-- = *from--;
+
+    ZERO(MP_DIGITS(z), need);
+    uz += need;
+  }
+
+  if(rest) {
+    d = 0;
+    for(i = need, from = MP_DIGITS(z) + need; i < uz; ++i, ++from) {
+      mp_digit save = *from;
+      
+      *from = (*from << rest) | (d >> (MP_DIGIT_BIT - rest));
+      d = save;
+    }
+
+    d >>= (MP_DIGIT_BIT - rest);
+    if(d != 0) {
+      *from = d;
+      uz += extra;
+    }
+  }
+
+  MP_USED(z) = uz;
+  CLAMP(z);
+
+  return 1;
+}
+
+/* }}} */
+
+/* {{{ s_qsub(z, p2) */
+
+/* Subtract |z| from 2^p2, assuming 2^p2 > |z|, and set z to be positive */
+static int       s_qsub(mp_int z, mp_size p2)
+{
+  mp_digit hi = (1 << (p2 % MP_DIGIT_BIT)), *zp;
+  mp_size  tdig = (p2 / MP_DIGIT_BIT), pos;
+  mp_word  w = 0;
+
+  if(!s_pad(z, tdig + 1))
+    return 0;
+
+  for(pos = 0, zp = MP_DIGITS(z); pos < tdig; ++pos, ++zp) {
+    w = ((mp_word) MP_DIGIT_MAX + 1) - w - (mp_word)*zp;
+
+    *zp = LOWER_HALF(w);
+    w = UPPER_HALF(w) ? 0 : 1;
+  }
+
+  w = ((mp_word) MP_DIGIT_MAX + 1 + hi) - w - (mp_word)*zp;
+  *zp = LOWER_HALF(w);
+
+  assert(UPPER_HALF(w) != 0); /* no borrow out should be possible */
+  
+  MP_SIGN(z) = MP_ZPOS;
+  CLAMP(z);
+
+  return 1;
+}
+
+/* }}} */
+
+/* {{{ s_dp2k(z) */
+
+static int      s_dp2k(mp_int z)
+{
+  int       k = 0;
+  mp_digit *dp = MP_DIGITS(z), d;
+
+  if(MP_USED(z) == 1 && *dp == 0)
+    return 1;
+
+  while(*dp == 0) {
+    k += MP_DIGIT_BIT;
+    ++dp;
+  }
+  
+  d = *dp;
+  while((d & 1) == 0) {
+    d >>= 1;
+    ++k;
+  }
+
+  return k;
+}
+
+/* }}} */
+
+/* {{{ s_isp2(z) */
+
+static int       s_isp2(mp_int z)
+{
+  mp_size uz = MP_USED(z), k = 0;
+  mp_digit *dz = MP_DIGITS(z), d;
+
+  while(uz > 1) {
+    if(*dz++ != 0)
+      return -1;
+    k += MP_DIGIT_BIT;
+    --uz;
+  }
+
+  d = *dz;
+  while(d > 1) {
+    if(d & 1)
+      return -1;
+    ++k; d >>= 1;
+  }
+
+  return (int) k;
+}
+
+/* }}} */
+
+/* {{{ s_2expt(z, k) */
+
+static int       s_2expt(mp_int z, int k)
+{
+  mp_size  ndig, rest;
+  mp_digit *dz;
+
+  ndig = (k + MP_DIGIT_BIT) / MP_DIGIT_BIT;
+  rest = k % MP_DIGIT_BIT;
+
+  if(!s_pad(z, ndig))
+    return 0;
+
+  dz = MP_DIGITS(z);
+  ZERO(dz, ndig);
+  *(dz + ndig - 1) = (1 << rest);
+  MP_USED(z) = ndig;
+
+  return 1;
+}
+
+/* }}} */
+
+/* {{{ s_norm(a, b) */
+
+static int      s_norm(mp_int a, mp_int b)
+{
+  mp_digit d = b->digits[MP_USED(b) - 1];
+  int      k = 0;
+
+  while(d < (mp_digit) (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. */
+mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c)
+{
+  mp_digit  *db, *dbt, umu, d;
+  mpz_t     temp[3]; 
+  mp_result res;
+  int       last = 0;
+
+  umu = MP_USED(mu); db = MP_DIGITS(b); dbt = db + MP_USED(b) - 1;
+
+  while(last < 3) 
+    SETUP(mp_int_init_size(TEMP(last), 4 * umu), last);
+
+  (void) mp_int_set_value(c, 1);
+
+  /* Take care of low-order digits */
+  while(db < dbt) {
+    int      i;
+
+    for(d = *db, i = MP_DIGIT_BIT; i > 0; --i, d >>= 1) {
+      if(d & 1) {
+	/* The use of a second temporary avoids allocation */
+	UMUL(c, a, TEMP(0));
+	if(!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
+	  res = MP_MEMORY; goto CLEANUP;
+	}
+	mp_int_copy(TEMP(0), c);
+      }
+
+
+      USQR(a, TEMP(0));
+      assert(MP_SIGN(TEMP(0)) == MP_ZPOS);
+      if(!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
+	res = MP_MEMORY; goto CLEANUP;
+      }
+      assert(MP_SIGN(TEMP(0)) == MP_ZPOS);
+      mp_int_copy(TEMP(0), a);
+
+
+    }
+
+    ++db;
+  }
+
+  /* Take care of highest-order digit */
+  d = *dbt;
+  for(;;) {
+    if(d & 1) {
+      UMUL(c, a, TEMP(0));
+      if(!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
+	res = MP_MEMORY; goto CLEANUP;
+      }
+      mp_int_copy(TEMP(0), c);
+    }
+    
+    d >>= 1;
+    if(!d) break;
+
+    USQR(a, TEMP(0));
+    if(!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
+      res = MP_MEMORY; goto CLEANUP;
+    }
+    (void) mp_int_copy(TEMP(0), a);
+  }
+
+ CLEANUP:
+  while(--last >= 0)
+    mp_int_clear(TEMP(last));
+  
+  return res;
+}
+
+/* }}} */
+
+/* {{{ s_udiv(a, b) */
+
+/* Precondition:  a >= b and b > 0
+   Postcondition: a' = a / b, b' = a % b
+ */
+static mp_result s_udiv(mp_int a, mp_int b)
+{
+  mpz_t     q, r, t;
+  mp_size   ua, ub, qpos = 0;
+  mp_digit *da, btop;
+  mp_result res = MP_OK;
+  int       k, skip = 0;
+
+  /* Force signs to positive */
+  MP_SIGN(a) = MP_ZPOS;
+  MP_SIGN(b) = MP_ZPOS;
+
+  /* Normalize, per Knuth */
+  k = s_norm(a, b);
+
+  ua = MP_USED(a); ub = MP_USED(b); btop = b->digits[ub - 1];
+  if((res = mp_int_init_size(&q, ua)) != MP_OK) return res;
+  if((res = mp_int_init_size(&t, ua + 1)) != MP_OK) goto CLEANUP;
+
+  da = MP_DIGITS(a);
+  r.digits = da + ua - 1;  /* The contents of r are shared with a */
+  r.used   = 1;
+  r.sign   = MP_ZPOS;
+  r.alloc  = MP_ALLOC(a);
+  ZERO(t.digits, t.alloc);
+
+  /* Solve for quotient digits, store in q.digits in reverse order */
+  while(r.digits >= da) {
+    if (qpos > q.alloc) {
+      char buf[1024];
+      printf("qpos = %d q.alloc = %d da = %d ua = %d\n",
+	     (int)qpos, (int)q.alloc, (int)da, (int)ua);
+      mp_int_to_string(a, 10, buf, sizeof(buf));
+      printf("a = %s\n", buf);
+      mp_int_to_string(b, 10, buf, sizeof(buf));
+      printf("b = %s\n", buf);
+      assert(qpos <= q.alloc);
+    }	
+
+    if(s_ucmp(b, &r) > 0) {
+      r.digits -= 1;
+      r.used += 1;
+      
+      if(++skip > 1)
+	q.digits[qpos++] = 0;
+      
+      CLAMP(&r);
+    }
+    else {
+      mp_word  pfx = r.digits[r.used - 1];
+      mp_word  qdigit;
+      
+      if(r.used > 1 && (pfx < btop || r.digits[r.used - 2] == 0)) {
+	pfx <<= MP_DIGIT_BIT / 2;
+	pfx <<= MP_DIGIT_BIT / 2;
+	pfx |= r.digits[r.used - 2];
+      }
+
+      qdigit = pfx / btop;
+      if(qdigit > MP_DIGIT_MAX) 
+	qdigit = 1;
+      
+      s_dbmul(MP_DIGITS(b), (mp_digit) qdigit, t.digits, ub);
+      t.used = ub + 1; CLAMP(&t);
+      while(s_ucmp(&t, &r) > 0) {
+	--qdigit;
+	(void) mp_int_sub(&t, b, &t); /* cannot fail */
+      }
+      
+      s_usub(r.digits, t.digits, r.digits, r.used, t.used);
+      CLAMP(&r);
+      
+      q.digits[qpos++] = (mp_digit) qdigit;
+      ZERO(t.digits, t.used);
+      skip = 0;
+    }
+  }
+  
+  /* Put quotient digits in the correct order, and discard extra zeroes */
+  q.used = qpos;
+  REV(mp_digit, q.digits, qpos);
+  CLAMP(&q);
+
+  /* Denormalize the remainder */
+  CLAMP(a);
+  if(k != 0)
+    s_qdiv(a, k);
+  
+  mp_int_copy(a, b);  /* ok:  0 <= r < b */
+  mp_int_copy(&q, a); /* ok:  q <= a     */
+  
+  mp_int_clear(&t);
+ CLEANUP:
+  mp_int_clear(&q);
+  return res;
+}
+
+/* }}} */
+
+/* {{{ s_outlen(z, r) */
+
+/* Precondition:  2 <= r < 64 */
+static int       s_outlen(mp_int z, mp_size r)
+{
+  mp_result  bits;
+  double     raw;
+
+  bits = mp_int_count_bits(z);
+  raw = (double)bits * s_log2[r];
+
+  return (int)(raw + 0.999999);
+}
+
+/* }}} */
+
+/* {{{ s_inlen(len, r) */
+
+static mp_size   s_inlen(int len, mp_size r)
+{
+  double  raw = (double)len / s_log2[r];
+  mp_size bits = (mp_size)(raw + 0.5);
+
+  return (mp_size)((bits + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT);
+}
+
+/* }}} */
+
+/* {{{ s_ch2val(c, r) */
+
+static int       s_ch2val(char c, int r)
+{
+  int out;
+
+  if(isdigit((unsigned char) c))
+    out = c - '0';
+  else if(r > 10 && isalpha((unsigned char) c))
+    out = toupper(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/des/imath/imath.h b/source4/heimdal/lib/des/imath/imath.h
new file mode 100755
index 0000000000..93cc35654d
--- /dev/null
+++ b/source4/heimdal/lib/des/imath/imath.h
@@ -0,0 +1,220 @@
+/*
+  Name:     imath.h
+  Purpose:  Arbitrary precision integer arithmetic routines.
+  Author:   M. J. Fromberger <http://www.dartmouth.edu/~sting/>
+  Info:     $Id: imath.h,v 1.3 2006/10/21 16:32:15 lha Exp $
+
+  Copyright (C) 2002 Michael J. Fromberger, All Rights Reserved.
+
+  Permission is hereby granted, free of charge, to any person
+  obtaining a copy of this software and associated documentation files
+  (the "Software"), to deal in the Software without restriction,
+  including without limitation the rights to use, copy, modify, merge,
+  publish, distribute, sublicense, and/or sell copies of the Software,
+  and to permit persons to whom the Software is furnished to do so,
+  subject to the following conditions:
+
+  The above copyright notice and this permission notice shall be
+  included in all copies or substantial portions of the Software.
+
+  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+  NONINFRINGEMENT.  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+  BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+  ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+  CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+  SOFTWARE.
+ */
+
+#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;
+#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;
+
+#define MP_DIGIT_BIT    (sizeof(mp_digit) * CHAR_BIT)
+#define MP_WORD_BIT     (sizeof(mp_word) * CHAR_BIT)
+
+#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  32
+
+#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, int value);
+mp_result mp_int_set_value(mp_int z, int value);
+void      mp_int_clear(mp_int z);
+void      mp_int_free(mp_int z);
+
+mp_result mp_int_copy(mp_int a, mp_int c);           /* c = a     */
+void      mp_int_swap(mp_int a, mp_int c);           /* swap a, c */
+void      mp_int_zero(mp_int z);                     /* z = 0     */
+mp_result mp_int_abs(mp_int a, mp_int c);            /* c = |a|   */
+mp_result mp_int_neg(mp_int a, mp_int c);            /* c = -a    */
+mp_result mp_int_add(mp_int a, mp_int b, mp_int c);  /* c = a + b */
+mp_result mp_int_add_value(mp_int a, int value, mp_int c);
+mp_result mp_int_sub(mp_int a, mp_int b, mp_int c);  /* c = a - b */
+mp_result mp_int_sub_value(mp_int a, int value, mp_int c);
+mp_result mp_int_mul(mp_int a, mp_int b, mp_int c);  /* c = a * b */
+mp_result mp_int_mul_value(mp_int a, int value, mp_int c);
+mp_result mp_int_mul_pow2(mp_int a, int p2, mp_int c);
+mp_result mp_int_sqr(mp_int a, mp_int c);            /* c = a * a */
+mp_result mp_int_div(mp_int a, mp_int b,             /* q = a / b */
+		     mp_int q, mp_int r);            /* r = a % b */
+mp_result mp_int_div_value(mp_int a, int value,      /* q = a / value */
+			   mp_int q, int *r);        /* r = a % value */
+mp_result mp_int_div_pow2(mp_int a, int p2,          /* q = a / 2^p2  */
+			  mp_int q, mp_int r);       /* r = q % 2^p2  */
+mp_result mp_int_mod(mp_int a, mp_int m, mp_int c);  /* c = a % m */
+#define   mp_int_mod_value(A, V, R) mp_int_div_value((A), (V), 0, (R))
+mp_result mp_int_expt(mp_int a, int b, mp_int c);    /* c = a^b   */
+mp_result mp_int_expt_value(int a, int b, mp_int c); /* c = a^b   */
+
+int       mp_int_compare(mp_int a, mp_int b);          /* a <=> b     */
+int       mp_int_compare_unsigned(mp_int a, mp_int b); /* |a| <=> |b| */
+int       mp_int_compare_zero(mp_int z);               /* a <=> 0     */
+int       mp_int_compare_value(mp_int z, int value);   /* a <=> v     */
+
+/* Returns true if v|a, false otherwise (including errors) */
+int       mp_int_divisible_value(mp_int a, int v);
+
+/* Returns k >= 0 such that z = 2^k, if one exists; otherwise < 0 */
+int       mp_int_is_pow2(mp_int z);
+
+mp_result mp_int_exptmod(mp_int a, mp_int b, mp_int m,
+			 mp_int c);                    /* c = a^b (mod m) */
+mp_result mp_int_exptmod_evalue(mp_int a, int value, 
+				mp_int m, mp_int c);   /* c = a^v (mod m) */
+mp_result mp_int_exptmod_bvalue(int value, mp_int b,
+				mp_int m, mp_int c);   /* c = v^b (mod m) */
+mp_result mp_int_exptmod_known(mp_int a, mp_int b,
+			       mp_int m, mp_int mu,
+			       mp_int c);              /* c = a^b (mod m) */
+mp_result mp_int_redux_const(mp_int m, mp_int c); 
+
+mp_result mp_int_invmod(mp_int a, mp_int m, mp_int c); /* c = 1/a (mod m) */
+
+mp_result mp_int_gcd(mp_int a, mp_int b, mp_int c);    /* c = gcd(a, b)   */
+
+mp_result mp_int_egcd(mp_int a, mp_int b, mp_int c,    /* c = gcd(a, b)   */
+		      mp_int x, mp_int y);             /* c = ax + by     */
+
+mp_result mp_int_sqrt(mp_int a, mp_int c);          /* c = floor(sqrt(q)) */
+
+/* Convert to an int, if representable (returns MP_RANGE if not). */
+mp_result mp_int_to_int(mp_int z, int *out);
+
+/* Convert to nul-terminated string with the specified radix, writing at
+   most limit characters including the nul terminator  */
+mp_result mp_int_to_string(mp_int z, mp_size radix, 
+			   char *str, int limit);
+
+/* Return the number of characters required to represent 
+   z in the given radix.  May over-estimate. */
+mp_result mp_int_string_len(mp_int z, mp_size radix);
+
+/* Read zero-terminated string into z */
+mp_result mp_int_read_string(mp_int z, mp_size radix, const char *str);
+mp_result mp_int_read_cstring(mp_int z, mp_size radix, const char *str, 
+			      char **end);
+
+/* Return the number of significant bits in z */
+mp_result mp_int_count_bits(mp_int z);
+
+/* Convert z to two's complement binary, writing at most limit bytes */
+mp_result mp_int_to_binary(mp_int z, unsigned char *buf, int limit);
+
+/* Read a two's complement binary value into z from the given buffer */
+mp_result mp_int_read_binary(mp_int z, unsigned char *buf, int len);
+
+/* Return the number of bytes required to represent z in binary. */
+mp_result mp_int_binary_len(mp_int z);
+
+/* Convert z to unsigned binary, writing at most limit bytes */
+mp_result mp_int_to_unsigned(mp_int z, unsigned char *buf, int limit);
+
+/* Read an unsigned binary value into z from the given buffer */
+mp_result mp_int_read_unsigned(mp_int z, unsigned char *buf, int len);
+
+/* Return the number of bytes required to represent z as unsigned output */
+mp_result mp_int_unsigned_len(mp_int z);
+
+/* Return a statically allocated string describing error code res */
+const char *mp_error_string(mp_result res);
+
+#if 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/des/imath/iprime.c b/source4/heimdal/lib/des/imath/iprime.c
new file mode 100755
index 0000000000..582ade0f54
--- /dev/null
+++ b/source4/heimdal/lib/des/imath/iprime.c
@@ -0,0 +1,186 @@
+/*
+  Name:     iprime.c
+  Purpose:  Pseudoprimality testing routines
+  Author:   M. J. Fromberger <http://www.dartmouth.edu/~sting/>
+  Info:     $Id: iprime.c,v 1.5 2007/01/05 21:01:48 lha Exp $
+
+  Copyright (C) 2002 Michael J. Fromberger, All Rights Reserved.
+
+  Permission is hereby granted, free of charge, to any person
+  obtaining a copy of this software and associated documentation files
+  (the "Software"), to deal in the Software without restriction,
+  including without limitation the rights to use, copy, modify, merge,
+  publish, distribute, sublicense, and/or sell copies of the Software,
+  and to permit persons to whom the Software is furnished to do so,
+  subject to the following conditions:
+
+  The above copyright notice and this permission notice shall be
+  included in all copies or substantial portions of the Software.
+
+  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+  NONINFRINGEMENT.  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+  BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+  ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+  CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+  SOFTWARE.
+ */
+
+#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, 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
+};
+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, 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/des/imath/iprime.h b/source4/heimdal/lib/des/imath/iprime.h
new file mode 100755
index 0000000000..cd54a73127
--- /dev/null
+++ b/source4/heimdal/lib/des/imath/iprime.h
@@ -0,0 +1,51 @@
+/*
+  Name:     iprime.h
+  Purpose:  Pseudoprimality testing routines
+  Author:   M. J. Fromberger <http://www.dartmouth.edu/~sting/>
+  Info:     $Id: iprime.h,v 1.3 2006/10/21 16:32:30 lha Exp $
+
+  Copyright (C) 2002 Michael J. Fromberger, All Rights Reserved.
+
+  Permission is hereby granted, free of charge, to any person
+  obtaining a copy of this software and associated documentation files
+  (the "Software"), to deal in the Software without restriction,
+  including without limitation the rights to use, copy, modify, merge,
+  publish, distribute, sublicense, and/or sell copies of the Software,
+  and to permit persons to whom the Software is furnished to do so,
+  subject to the following conditions:
+
+  The above copyright notice and this permission notice shall be
+  included in all copies or substantial portions of the Software.
+
+  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+  NONINFRINGEMENT.  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+  BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+  ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+  CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+  SOFTWARE.
+ */
+
+#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/des/pkcs12.c b/source4/heimdal/lib/des/pkcs12.c
new file mode 100644
index 0000000000..cc92285754
--- /dev/null
+++ b/source4/heimdal/lib/des/pkcs12.c
@@ -0,0 +1,145 @@
+/*
+ * 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. 
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+RCSID("$Id: pkcs12.c,v 1.1 2006/01/13 08:26:49 lha Exp $");
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+
+#include <pkcs12.h>
+#include <bn.h>
+
+#include <roken.h>
+
+int
+PKCS12_key_gen(const void *key, size_t keylen, 
+	       const void *salt, size_t saltlen,
+	       int id, int iteration, size_t outkeysize, 
+	       void *out, const EVP_MD *md)
+{
+    unsigned char *v, *I, hash[EVP_MAX_MD_SIZE];
+    unsigned int size, size_I = 0;
+    unsigned char idc = id;
+    EVP_MD_CTX ctx;
+    unsigned char *outp = out;
+    int i, vlen;
+
+    EVP_MD_CTX_init(&ctx);
+
+    vlen = EVP_MD_block_size(md);
+    v = malloc(vlen + 1);
+    if (v == NULL)
+	return 0;
+
+    I = calloc(1, vlen * 2);
+    if (I == NULL) {
+	free(v);
+	return 0;
+    }
+
+    if (salt && saltlen > 0) {
+	for (i = 0; i < vlen; i++)
+	    I[i] = ((unsigned char*)salt)[i % saltlen];
+	size_I += vlen;
+    }
+    if (key && keylen > 0) {
+	for (i = 0; i < vlen / 2; i++) {
+	    I[(i * 2) + size_I] = 0;
+	    I[(i * 2) + size_I + 1] = ((unsigned char*)key)[i % (keylen + 1)];
+	}
+	size_I += vlen;
+    }
+
+    while (1) {
+	BIGNUM *bnB, *bnOne;
+
+	if (!EVP_DigestInit_ex(&ctx, md, NULL))
+	    return 0;
+	for (i = 0; i < vlen; i++)
+	    EVP_DigestUpdate(&ctx, &idc, 1);
+	EVP_DigestUpdate(&ctx, I, size_I);
+	EVP_DigestFinal_ex(&ctx, hash, &size);
+
+	for (i = 1; i < iteration; i++)
+	    EVP_Digest(hash, size, hash, &size, md, NULL);
+
+	memcpy(outp, hash, min(outkeysize, size));
+	if (outkeysize < size)
+	    break;
+	outkeysize -= size;
+	outp += size;
+
+	for (i = 0; i < vlen; i++) 
+	    v[i] = hash[i % size];
+
+	bnB = BN_bin2bn(v, vlen, NULL);
+	bnOne = BN_new();
+	BN_set_word(bnOne, 1);
+
+	BN_uadd(bnB, bnB, bnOne);
+
+	for (i = 0; i < vlen * 2; i += vlen) {
+	    BIGNUM *bnI;
+	    int j;
+
+	    bnI = BN_bin2bn(I + i, vlen, NULL);
+
+	    BN_uadd(bnI, bnI, bnB);
+
+	    j = BN_num_bytes(bnI);
+	    if (j > vlen) {
+		assert(j == vlen + 1);
+		BN_bn2bin(bnI, v);
+		memcpy(I + i, v + 1, vlen);
+	    } else {
+		memset(I + i, 0, vlen - j);
+		BN_bn2bin(bnI, I + i + vlen - j);
+	    }
+	    BN_free(bnI);
+	}	    
+	BN_free(bnB);
+	BN_free(bnOne);
+	size_I = vlen * 2;
+    }
+
+    EVP_MD_CTX_cleanup(&ctx);
+    free(I);
+    free(v);
+
+    return 1;
+}
diff --git a/source4/heimdal/lib/des/resource.h b/source4/heimdal/lib/des/resource.h
new file mode 100644
index 0000000000..02c6a7c6d9
--- /dev/null
+++ b/source4/heimdal/lib/des/resource.h
@@ -0,0 +1,18 @@
+//{{NO_DEPENDENCIES}}
+// Microsoft Developer Studio generated include file.
+// Used by passwd_dialog.rc
+//
+#define IDD_PASSWD_DIALOG               101
+#define IDC_EDIT1                       1000
+#define IDC_PASSWD_EDIT                 1001
+
+// Next default values for new objects
+// 
+#ifdef APSTUDIO_INVOKED
+#ifndef APSTUDIO_READONLY_SYMBOLS
+#define _APS_NEXT_RESOURCE_VALUE        102
+#define _APS_NEXT_COMMAND_VALUE         40001
+#define _APS_NEXT_CONTROL_VALUE         1002
+#define _APS_NEXT_SYMED_VALUE           101
+#endif
+#endif
diff --git a/source4/heimdal/lib/des/rsa-imath.c b/source4/heimdal/lib/des/rsa-imath.c
new file mode 100644
index 0000000000..298affadfe
--- /dev/null
+++ b/source4/heimdal/lib/des/rsa-imath.c
@@ -0,0 +1,661 @@
+/*
+ * 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. 
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+RCSID("$Id: rsa-imath.c,v 1.23 2007/01/06 13:45:25 lha Exp $");
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <krb5-types.h>
+#include <assert.h>
+
+#include <rsa.h>
+
+#include <roken.h>
+
+#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;
+    assert(padlen >= 8);
+
+    *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);
+    {
+	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;
+    size_t size;
+    mpz_t in, out, n, e, b, bi;
+    int blinding = (rsa->flags & RSA_FLAG_NO_BLINDING) == 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 = 0;
+	goto out;
+    }
+
+    if (blinding) {
+	setup_blind(&n, &b, &bi);
+	blind(&in, &b, &e, &n);
+    }
+
+    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);
+    } 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 = 0;
+	    goto out;
+	}
+    }
+
+    if (blinding) {
+	unblind(&out, &bi, &n);
+	mp_int_clear(&b);
+	mp_int_clear(&bi);
+    }
+
+    {
+	size_t ssize;
+	ssize = mp_int_unsigned_len(&out);
+	assert(size >= ssize);
+	mp_int_to_unsigned(&out, to, size);
+	size = ssize;
+    }
+    
+out:
+    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;
+    size_t size;
+    mpz_t in, out, n, e, b, bi;
+    int blinding = (rsa->flags & RSA_FLAG_NO_BLINDING) == 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 = 0;
+	goto out;
+    }
+
+    if (blinding) {
+	setup_blind(&n, &b, &bi);
+	blind(&in, &b, &e, &n);
+    }
+
+    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);
+    } 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 = 0;
+	    goto out;
+	}
+    }
+
+    if (blinding) {
+	unblind(&out, &bi, &n);
+	mp_int_clear(&b);
+	mp_int_clear(&bi);
+    }
+
+    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)
+	return -3;
+    size--; ptr++;
+    while (size && *ptr != 0) {
+	size--; ptr++;
+    }
+    if (size == 0)
+	return -4;
+    size--; ptr++;
+
+    memmove(to, ptr, size);
+
+out:
+    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
+};
+
+const RSA_METHOD *
+RSA_imath_method(void)
+{
+    return &hc_rsa_imath_method;
+}
diff --git a/source4/heimdal/lib/des/rsa.c b/source4/heimdal/lib/des/rsa.c
new file mode 100644
index 0000000000..241afb2e46
--- /dev/null
+++ b/source4/heimdal/lib/des/rsa.c
@@ -0,0 +1,471 @@
+/*
+ * 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. 
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+RCSID("$Id: rsa.c,v 1.19 2007/01/09 10:04:20 lha Exp $");
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <krb5-types.h>
+#include <rfc2459_asn1.h>
+
+#include <rsa.h>
+
+#include <roken.h>
+
+RSA *
+RSA_new(void)
+{
+    return RSA_new_method(NULL);
+}
+
+RSA *
+RSA_new_method(ENGINE *engine)
+{
+    RSA *rsa;
+
+    rsa = calloc(1, sizeof(*rsa));
+    if (rsa == NULL)
+	return NULL;
+
+    rsa->references = 1;
+
+    if (engine) {
+	ENGINE_up_ref(engine);
+	rsa->engine = engine;
+    } else {
+	rsa->engine = ENGINE_get_default_RSA();
+    }
+
+    if (rsa->engine) {
+	rsa->meth = ENGINE_get_RSA(rsa->engine);
+	if (rsa->meth == NULL) {
+	    ENGINE_finish(engine);
+	    free(rsa);
+	    return 0;
+	}
+    }
+
+    if (rsa->meth == NULL)
+	rsa->meth = rk_UNCONST(RSA_get_default_method());
+
+    (*rsa->meth->init)(rsa);
+
+    return rsa;
+}
+
+
+void
+RSA_free(RSA *rsa)
+{
+    if (rsa->references <= 0)
+	abort();
+
+    if (--rsa->references > 0)
+	return;
+
+    (*rsa->meth->finish)(rsa);
+
+    if (rsa->engine)
+	ENGINE_finish(rsa->engine);
+
+#define free_if(f) if (f) { BN_free(f); }
+    free_if(rsa->n);
+    free_if(rsa->e);
+    free_if(rsa->d);
+    free_if(rsa->p);
+    free_if(rsa->q);
+    free_if(rsa->dmp1);
+    free_if(rsa->dmq1);
+#undef free_if
+
+    memset(rsa, 0, sizeof(*rsa));
+    free(rsa);
+}
+
+int
+RSA_up_ref(RSA *rsa)
+{
+    return ++rsa->references;
+}
+
+const RSA_METHOD *
+RSA_get_method(const RSA *rsa)
+{
+    return rsa->meth;
+}
+
+int
+RSA_set_method(RSA *rsa, const RSA_METHOD *method)
+{
+    (*rsa->meth->finish)(rsa);
+
+    if (rsa->engine) {
+	ENGINE_finish(rsa->engine);
+	rsa->engine = NULL;
+    }
+
+    rsa->meth = method;
+    (*rsa->meth->init)(rsa);
+    return 1;
+}
+
+int
+RSA_set_app_data(RSA *rsa, void *arg)
+{
+    rsa->ex_data.sk = arg;
+    return 1;
+}
+
+void *
+RSA_get_app_data(RSA *rsa)
+{
+    return rsa->ex_data.sk;
+}
+
+int
+RSA_check_key(const RSA *key)
+{
+    static const unsigned char inbuf[] = "hello, world!";
+    RSA *rsa = rk_UNCONST(key);
+    void *buffer;
+    int ret;
+
+    /* 
+     * XXX I have no clue how to implement this w/o a bignum library.
+     * Well, when we have a RSA key pair, we can try to encrypt/sign
+     * and then decrypt/verify.
+     */
+
+    if ((rsa->d == NULL || rsa->n == NULL) &&
+	(rsa->p == NULL || rsa->q || rsa->dmp1 == NULL || rsa->dmq1 == NULL || rsa->iqmp == NULL))
+	return 0;
+
+    buffer = malloc(RSA_size(rsa));
+    if (buffer == NULL)
+	return 0;
+    
+    ret = RSA_private_encrypt(sizeof(inbuf), inbuf, buffer, 
+			     rsa, RSA_PKCS1_PADDING);
+    if (ret == -1) {
+	free(buffer);
+	return 0;
+    }
+
+    ret = RSA_public_decrypt(ret, buffer, buffer,
+			      rsa, RSA_PKCS1_PADDING);
+    if (ret == -1) {
+	free(buffer);
+	return 0;
+    }
+
+    if (ret == sizeof(inbuf) && memcmp(buffer, inbuf, sizeof(inbuf)) == 0) {
+	free(buffer);
+	return 1;
+    }
+    free(buffer);
+    return 0; 
+}
+
+int
+RSA_size(const RSA *rsa)
+{
+    return BN_num_bytes(rsa->n);
+}
+
+#define RSAFUNC(name, body) \
+int \
+name(int flen,const unsigned char* f, unsigned char* t, RSA* r, int p){\
+    return body; \
+}
+
+RSAFUNC(RSA_public_encrypt, (r)->meth->rsa_pub_enc(flen, f, t, r, p))
+RSAFUNC(RSA_public_decrypt, (r)->meth->rsa_pub_dec(flen, f, t, r, p))
+RSAFUNC(RSA_private_encrypt, (r)->meth->rsa_priv_enc(flen, f, t, r, p))
+RSAFUNC(RSA_private_decrypt, (r)->meth->rsa_priv_dec(flen, f, t, r, p))
+
+/* XXX */
+int
+RSA_sign(int type, const unsigned char *from, unsigned int flen,
+	 unsigned char *to, unsigned int *tlen, RSA *rsa)
+{
+    return -1;
+}
+
+int
+RSA_verify(int type, const unsigned char *from, unsigned int flen,
+	   unsigned char *to, unsigned int tlen, RSA *rsa)
+{
+    return -1;
+}
+
+/*
+ * A NULL RSA_METHOD that returns failure for all operations. This is
+ * used as the default RSA method is we don't have any native
+ * support.
+ */
+
+static RSAFUNC(null_rsa_public_encrypt, -1)
+static RSAFUNC(null_rsa_public_decrypt, -1)
+static RSAFUNC(null_rsa_private_encrypt, -1)
+static RSAFUNC(null_rsa_private_decrypt, -1)
+
+/*
+ *
+ */
+
+int
+RSA_generate_key_ex(RSA *r, int bits, BIGNUM *e, BN_GENCB *cb)
+{
+    if (r->meth->rsa_keygen)
+	return (*r->meth->rsa_keygen)(r, bits, e, cb);
+    return 0;
+}
+
+
+/*
+ *
+ */
+
+static int 
+null_rsa_init(RSA *rsa)
+{
+    return 1;
+}
+
+static int
+null_rsa_finish(RSA *rsa)
+{
+    return 1;
+}
+
+static const RSA_METHOD rsa_null_method = {
+    "hcrypto null RSA",
+    null_rsa_public_encrypt,
+    null_rsa_public_decrypt,
+    null_rsa_private_encrypt,
+    null_rsa_private_decrypt,
+    NULL,
+    NULL,
+    null_rsa_init,
+    null_rsa_finish,
+    0,
+    NULL,
+    NULL,
+    NULL
+};
+
+const RSA_METHOD *
+RSA_null_method(void)
+{
+    return &rsa_null_method;
+}
+
+extern const RSA_METHOD hc_rsa_imath_method;
+static const RSA_METHOD *default_rsa_method = &hc_rsa_imath_method;
+
+const RSA_METHOD *
+RSA_get_default_method(void)
+{
+    return default_rsa_method;
+}
+
+void
+RSA_set_default_method(const RSA_METHOD *meth)
+{
+    default_rsa_method = meth;
+}
+
+/*
+ *
+ */
+
+static BIGNUM *
+heim_int2BN(const heim_integer *i)
+{
+    BIGNUM *bn;
+
+    bn = BN_bin2bn(i->data, i->length, NULL);
+    if (bn)
+	BN_set_negative(bn, i->negative);
+    return bn;
+}
+
+static int
+bn2heim_int(BIGNUM *bn, heim_integer *integer)
+{
+    integer->length = BN_num_bytes(bn);
+    integer->data = malloc(integer->length);
+    if (integer->data == NULL) {
+	integer->length = 0;
+	return ENOMEM;
+    }
+    BN_bn2bin(bn, integer->data);
+    integer->negative = BN_is_negative(bn);
+    return 0;
+}
+
+
+RSA *
+d2i_RSAPrivateKey(RSA *rsa, const unsigned char **pp, size_t len)
+{
+    RSAPrivateKey data;
+    RSA *k = rsa;
+    size_t size;
+    int ret;
+
+    ret = decode_RSAPrivateKey(*pp, len, &data, &size);
+    if (ret)
+	return NULL;
+    
+    *pp += size;
+
+    if (k == NULL) {
+	k = RSA_new();
+	if (k == NULL) {
+	    free_RSAPrivateKey(&data);
+	    return NULL;
+	}
+    }
+
+    k->n = heim_int2BN(&data.modulus);
+    k->e = heim_int2BN(&data.publicExponent);
+    k->d = heim_int2BN(&data.privateExponent);
+    k->p = heim_int2BN(&data.prime1);
+    k->q = heim_int2BN(&data.prime2);
+    k->dmp1 = heim_int2BN(&data.exponent1);
+    k->dmq1 = heim_int2BN(&data.exponent2);
+    k->iqmp = heim_int2BN(&data.coefficient);
+    free_RSAPrivateKey(&data);
+
+    if (k->n == NULL || k->e == NULL || k->d == NULL || k->p == NULL ||
+	k->q == NULL || k->dmp1 == NULL || k->dmq1 == NULL || k->iqmp == NULL) 
+    {
+	RSA_free(k);
+	return NULL;
+    }
+	
+    return k;
+}
+
+int
+i2d_RSAPrivateKey(RSA *rsa, unsigned char **pp)
+{
+    RSAPrivateKey data;
+    size_t size;
+    int ret;
+
+    if (rsa->n == NULL || rsa->e == NULL || rsa->d == NULL || rsa->p == NULL ||
+	rsa->q == NULL || rsa->dmp1 == NULL || rsa->dmq1 == NULL ||
+	rsa->iqmp == NULL)
+	return -1;
+
+    memset(&data, 0, sizeof(data));
+
+    ret  = bn2heim_int(rsa->n, &data.modulus);
+    ret |= bn2heim_int(rsa->e, &data.publicExponent);
+    ret |= bn2heim_int(rsa->d, &data.privateExponent);
+    ret |= bn2heim_int(rsa->p, &data.prime1);
+    ret |= bn2heim_int(rsa->q, &data.prime2);
+    ret |= bn2heim_int(rsa->dmp1, &data.exponent1);
+    ret |= bn2heim_int(rsa->dmq1, &data.exponent2);
+    ret |= bn2heim_int(rsa->iqmp, &data.coefficient);
+    if (ret) {
+	free_RSAPrivateKey(&data);
+	return -1;
+    }
+
+    if (pp == NULL) {
+	size = length_RSAPrivateKey(&data);
+	free_RSAPrivateKey(&data);
+    } else {
+	void *p;
+	size_t len;
+
+	ASN1_MALLOC_ENCODE(RSAPrivateKey, p, len, &data, &size, ret);
+	free_RSAPrivateKey(&data);
+	if (ret)
+	    return -1;
+	if (len != size)
+	    abort();
+
+	memcpy(*pp, p, size);
+	free(p);
+
+	*pp += size;
+
+    }
+    return size;
+}
+
+int
+i2d_RSAPublicKey(RSA *rsa, unsigned char **pp)
+{
+    RSAPublicKey data;
+    size_t size;
+    int ret;
+
+    memset(&data, 0, sizeof(data));
+
+    if (bn2heim_int(rsa->n, &data.modulus) ||
+	bn2heim_int(rsa->e, &data.publicExponent))
+    {
+	free_RSAPublicKey(&data);
+	return -1;
+    }
+
+    if (pp == NULL) {
+	size = length_RSAPublicKey(&data);
+	free_RSAPublicKey(&data);
+    } else {
+	void *p;
+	size_t len;
+
+	ASN1_MALLOC_ENCODE(RSAPublicKey, p, len, &data, &size, ret);
+	free_RSAPublicKey(&data);
+	if (ret)
+	    return -1;
+	if (len != size)
+	    abort();
+    
+	memcpy(*pp, p, size);
+	free(p);
+
+	*pp += size;
+    }
+    
+    return size;
+}
diff --git a/source4/heimdal/lib/des/rsa.h b/source4/heimdal/lib/des/rsa.h
index 137dd9894b..0aceb9f9da 100644
--- a/source4/heimdal/lib/des/rsa.h
+++ b/source4/heimdal/lib/des/rsa.h
@@ -32,7 +32,7 @@
  */
 
 /*
- * $Id: rsa.h,v 1.5 2006/05/07 11:34:02 lha Exp $
+ * $Id: rsa.h,v 1.9 2007/01/05 20:26:23 lha Exp $
  */
 
 #ifndef _HEIM_RSA_H
@@ -59,7 +59,9 @@
 #define RSA_private_decrypt hc_RSA_private_decrypt
 #define RSA_sign hc_RSA_sign
 #define RSA_verify hc_RSA_verify
+#define RSA_generate_key_ex hc_RSA_generate_key_ex
 #define d2i_RSAPrivateKey hc_d2i_RSAPrivateKey
+#define i2d_RSAPrivateKey hc_i2d_RSAPrivateKey
 #define i2d_RSAPublicKey hc_i2d_RSAPublicKey
 
 /*
@@ -119,9 +121,10 @@ struct RSA {
     void *mt_blinding;
 };
 
-#define RSA_FLAG_SIGN_VER		0x40
+#define RSA_FLAG_NO_BLINDING		0x0080
 
 #define RSA_PKCS1_PADDING		1
+#define RSA_PKCS1_OAEP_PADDING		4
 #define RSA_PKCS1_PADDING_SIZE		11
 
 /*
@@ -162,7 +165,11 @@ int RSA_sign(int, const unsigned char *, unsigned int,
 int RSA_verify(int, const unsigned char *, unsigned int,
 	       unsigned char *, unsigned int, RSA *);
 
+int	RSA_generate_key_ex(RSA *, int, BIGNUM *, BN_GENCB *);
+
 RSA *	d2i_RSAPrivateKey(RSA *, const unsigned char **, size_t);
+int	i2d_RSAPrivateKey(RSA *, unsigned char **);
+
 int	i2d_RSAPublicKey(RSA *, unsigned char **);
 
 #endif /* _HEIM_RSA_H */