Initial revision

4 files changed, 3349 insertions, 0 deletions

diff --git a/lisp/mp/mp.c b/lisp/mp/mp.c
new file mode 100644
index 0000000..78b7a0e
--- /dev/null
+++ b/lisp/mp/mp.c
@@ -0,0 +1,822 @@
+/*
+ * Copyright (c) 2002 by The XFree86 Project, Inc.
+ *
+ * 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 XFREE86 PROJECT 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.
+ *
+ * Except as contained in this notice, the name of the XFree86 Project shall
+ * not be used in advertising or otherwise to promote the sale, use or other
+ * dealings in this Software without prior written authorization from the
+ * XFree86 Project.
+ *
+ * Author: Paulo César Pereira de Andrade
+ */
+
+/* $XFree86: xc/programs/xedit/lisp/mp/mp.c,v 1.3 2002/11/20 07:44:43 paulo Exp $ */
+
+#include "mp.h"
+
+/*
+ * TODO:
+ *	o Optimize squaring
+ *	o Write better division code and move from mpi.c to here
+ *	o Make multiplication code don't required memory to be zeroed
+ *		+ The first step is easy, just multiply the low word,
+ *		then the high word, that may overlap with the result
+ *		of the first multiply (in case of carry), and then
+ *		just make sure carry is properly propagated in the
+ *		subsequent multiplications.
+ *		+ Some code needs also to be rewritten because some
+ *		intermediate addition code in mp_mul, mp_karatsuba_mul,
+ *		and mp_toom_mul is assuming the memory is zeroed.
+ */
+
+/*
+ * Prototypes
+ */
+	/* out of memory handler */
+static void mp_outmem(void);
+
+	/* memory allocation fallback functions */
+static void *_mp_malloc(size_t);
+static void *_mp_calloc(size_t, size_t);
+static void *_mp_realloc(void*, size_t);
+static void _mp_free(void*);
+
+/*
+ * Initialization
+ */
+static mp_malloc_fun __mp_malloc = _mp_malloc;
+static mp_calloc_fun __mp_calloc = _mp_calloc;
+static mp_realloc_fun __mp_realloc = _mp_realloc;
+static mp_free_fun __mp_free = _mp_free;
+
+/*
+ * Implementation
+ */
+static void
+mp_outmem(void)
+{
+    fprintf(stderr, "out of memory in MP library.\n");
+    exit(1);
+}
+
+static void *
+_mp_malloc(size_t size)
+{
+    return (malloc(size));
+}
+
+void *
+mp_malloc(size_t size)
+{
+    void *pointer = (*__mp_malloc)(size);
+
+    if (pointer == NULL)
+	mp_outmem();
+
+    return (pointer);
+}
+
+mp_malloc_fun
+mp_set_malloc(mp_malloc_fun fun)
+{
+    mp_malloc_fun old = __mp_malloc;
+
+    __mp_malloc = fun;
+
+    return (old);
+}
+
+static void *
+_mp_calloc(size_t nmemb, size_t size)
+{
+    return (calloc(nmemb, size));
+}
+
+void *
+mp_calloc(size_t nmemb, size_t size)
+{
+    void *pointer = (*__mp_calloc)(nmemb, size);
+
+    if (pointer == NULL)
+	mp_outmem();
+
+    return (pointer);
+}
+
+mp_calloc_fun
+mp_set_calloc(mp_calloc_fun fun)
+{
+    mp_calloc_fun old = __mp_calloc;
+
+    __mp_calloc = fun;
+
+    return (old);
+}
+
+static void *
+_mp_realloc(void *old, size_t size)
+{
+    return (realloc(old, size));
+}
+
+void *
+mp_realloc(void *old, size_t size)
+{
+    void *pointer = (*__mp_realloc)(old, size);
+
+    if (pointer == NULL)
+	mp_outmem();
+
+    return (pointer);
+}
+
+mp_realloc_fun
+mp_set_realloc(mp_realloc_fun fun)
+{
+    mp_realloc_fun old = __mp_realloc;
+
+    __mp_realloc = fun;
+
+    return (old);
+}
+
+static void
+_mp_free(void *pointer)
+{
+    free(pointer);
+}
+
+void
+mp_free(void *pointer)
+{
+    (*__mp_free)(pointer);
+}
+
+mp_free_fun
+mp_set_free(mp_free_fun fun)
+{
+    mp_free_fun old = __mp_free;
+
+    __mp_free = fun;
+
+    return (old);
+}
+
+long
+mp_add(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2)
+{
+    BNI value;			/* intermediate result */
+    BNS carry;			/* carry flag */
+    long size;			/* result size */
+
+    if (len1 < len2)
+	MP_SWAP(op1, op2, len1, len2);
+
+    /* unroll start of loop */
+    value = op1[0] + op2[0];
+    rop[0] = value;
+    carry = value >> BNSBITS;
+
+    /* add op1 and op2 */
+    for (size = 1; size < len2; size++) {
+	value = op1[size] + op2[size] + carry;
+	rop[size] = value;
+	carry = value >> BNSBITS;
+    }
+    if (rop != op1) {
+	for (; size < len1; size++) {
+	    value = op1[size] + carry;
+	    rop[size] = value;
+	    carry = value >> BNSBITS;
+	}
+    }
+    else {
+	/* if rop == op1, than just adjust carry */
+	for (; carry && size < len1; size++) {
+	    value = op1[size] + carry;
+	    rop[size] = value;
+	    carry = value >> BNSBITS;
+	}
+	size = len1;
+    }
+    if (carry)
+	rop[size++] = carry;
+
+    return (size);
+}
+
+long
+mp_sub(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2)
+{
+    long svalue;		/* intermediate result */
+    BNS carry;			/* carry flag */
+    long size;			/* result size */
+
+    /* special case */
+    if (op1 == op2) {
+	rop[0] = 0;
+
+	return (1);
+    }
+
+    /* unroll start of loop */
+    svalue = op1[0] - op2[0];
+    rop[0] = svalue;
+    carry = svalue < 0;
+
+    /* subtracts op2 from op1 */
+    for (size = 1; size < len2; size++) {
+	svalue = (long)(op1[size]) - op2[size] - carry;
+	rop[size] = svalue;
+	carry = svalue < 0;
+    }
+    if (rop != op1) {
+	for (; size < len1; size++) {
+	    svalue = op1[size] - carry;
+	    rop[size] = svalue;
+	    carry = svalue < 0;
+	}
+    }
+    else {
+	/* if rop == op1, than just adjust carry */
+	for (; carry && size < len1; size++) {
+	    svalue = op1[size] - carry;
+	    rop[size] = svalue;
+	    carry = svalue < 0;
+	}
+	size = len1;
+    }
+
+    /* calculate result size */
+    while (size > 1 && rop[size - 1] == 0)
+	--size;
+
+    return (size);
+}
+
+long
+mp_lshift(BNS *rop, BNS *op, BNI len, long shift)
+{
+    long i, size;
+    BNI words, bits;		/* how many word and bit shifts */
+
+    words = shift / BNSBITS;
+    bits = shift % BNSBITS;
+    size = len + words;
+
+    if (bits) {
+	BNS hi, lo;
+	BNI carry;
+	int adj;
+
+	for (i = 1, carry = CARRY >> 1; carry; i++, carry >>= 1)
+	    if (op[len - 1] & carry)
+		break;
+	adj = (bits + (BNSBITS - i)) / BNSBITS;
+	size += adj;
+
+	lo = hi = op[0];
+	rop[words] = lo << bits;
+	for (i = 1; i < len; i++) {
+	    hi = op[i];
+	    rop[words + i] = hi << bits | (lo >> (BNSBITS - bits));
+	    lo = hi;
+	}
+	if (adj)
+	    rop[size - 1] = hi >> (BNSBITS - bits);
+    }
+    else
+	memmove(rop + size - len, op, sizeof(BNS) * len);
+
+    if (words)
+	memset(rop, '\0', sizeof(BNS) * words);
+
+    return (size);
+}
+
+long
+mp_rshift(BNS *rop, BNS *op, BNI len, long shift)
+{
+    int adj = 0;
+    long i, size;
+    BNI words, bits;		/* how many word and bit shifts */
+
+    words = shift / BNSBITS;
+    bits = shift % BNSBITS;
+    size = len - words;
+
+    if (bits) {
+	BNS hi, lo;
+	BNI carry;
+
+	for (i = 0, carry = CARRY >> 1; carry; i++, carry >>= 1)
+	    if (op[len - 1] & carry)
+		break;
+	adj = (bits + i) / BNSBITS;
+	if (size - adj == 0) {
+	    rop[0] = 0;
+
+	    return (1);
+	}
+
+	hi = lo = op[words + size - 1];
+	rop[size - 1] = hi >> bits;
+	for (i = size - 2; i >= 0; i--) {
+	    lo = op[words + i];
+	    rop[i] = (lo >> bits) | (hi << (BNSBITS - bits));
+	    hi = lo;
+	}
+	if (adj)
+	    rop[0] |= lo << (BNSBITS - bits);
+    }
+    else
+	memmove(rop, op + len - size, size * sizeof(BNS));
+
+    return (size - adj);
+}
+
+	/* rop must be a pointer to len1 + len2 elements
+	 * rop cannot be either op1 or op2
+	 * rop must be all zeros */
+long
+mp_base_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2)
+{
+    long i, j;			/* counters */
+    BNI value;			/* intermediate result */
+    BNS carry;			/* carry value */
+    long size = len1 + len2;
+
+    /* simple optimization: first pass does not need to deference rop[i+j] */
+    if (op1[0]) {
+	value = (BNI)(op1[0]) * op2[0];
+	rop[0] = value;
+	carry = (BNS)(value >> BNSBITS);
+	for (j = 1; j < len2; j++) {
+	    value = (BNI)(op1[0]) * op2[j] + carry;
+	    rop[j] = value;
+	    carry = (BNS)(value >> BNSBITS);
+	}
+	rop[j] = carry;
+    }
+
+    /* do the multiplication */
+    for (i = 1; i < len1; i++) {
+	if (op1[i]) {
+	    /* unrool loop initialization */
+	    value = (BNI)(op1[i]) * op2[0] + rop[i];
+	    rop[i] = value;
+	    carry = (BNS)(value >> BNSBITS);
+	    /* multiply */
+	    for (j = 1; j < len2; j++) {
+		value = (BNI)(op1[i]) * op2[j] + rop[i + j] + carry;
+		rop[i + j] = value;
+		carry = (BNS)(value >> BNSBITS);
+	    }
+	    rop[i + j] = carry;
+	}
+    }
+
+    if (size > 1 && rop[size - 1] == 0)
+	--size;
+
+    return (size);
+}
+
+	/* Karatsuba method
+	 * t + ((a0 + a1) (b0 + b1) - t - u) x + ux²
+	 * where t = a0b0 and u = a1b1
+	 *
+	 * Karatsuba method reduces the number of multiplications. Example:
+	 *	Square a 40 length number
+	 *	instead of a plain 40*40 = 1600 multiplies/adds, it does:
+	 *	20*20+20*20+20*20 = 1200
+	 *	but since it is recursive, every 20*20=400 is reduced to
+	 *	10*10+10*10+10*10=300
+	 *	and so on.
+	 * The multiplication by x and x² is a just a shift, as it is a
+	 * power of two, and is implemented below by just writting at the
+	 * correct offset */
+long
+mp_karatsuba_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2)
+{
+    BNI x;				/* shift count */
+    BNI la0, la1, lb0, lb1;		/* length of a0, a1, b0, and b1 */
+    BNS *t;				/* temporary memory for t product */
+    BNS *u;				/* temporary memory for u product */
+    BNS *r;				/* pointer to rop */
+    long xlen, tlen, ulen;
+
+    /* calculate value of x, that is 2^(BNSBITS*x) */
+    if (len1 >= len2)
+	x = (len1 + 1) >> 1;
+    else
+	x = (len2 + 1) >> 1;
+
+    /* calculate length of operands */
+    la0 = x;
+    la1 = len1 - x;
+    lb0 = x;
+    lb1 = len2 - x;
+
+    /* allocate buffer for t and (a0 + a1) */
+    tlen = la0 + lb0;
+    t = mp_malloc(sizeof(BNS) * tlen);
+
+    /* allocate buffer for u and (b0 + b1) */
+    if (la1 + lb1 < lb0 + lb1 + 1)
+	ulen = lb0 + lb1 + 1;
+    else
+	ulen = la1 + lb1;
+    u = mp_malloc(sizeof(BNS) * ulen);
+
+    /* calculate a0 + a1, store result in t */
+    tlen = mp_add(t, op1, op1 + x, la0, la1);
+
+    /* calculate b0 + b1, store result in u */
+    ulen = mp_add(u, op2, op2 + x, lb0, lb1);
+
+    /* store (a0 + a1) * (b0 + b1) in rop */
+
+    r = rop + x;		/* multiplied by 2^(BNSBITS*x) */
+    xlen = mp_mul(r, t, u, tlen, ulen);
+
+    /* must zero t and u memory, this is required for mp_mul */
+
+    /* calculate t = a0 * b0 */
+    tlen = la0 + lb0;
+    memset(t, '\0', sizeof(BNS) * tlen);
+    tlen = mp_mul(t, op1, op2, la0, lb0);
+
+    /* calculate u = a1 * b1 */
+    ulen = la1 + lb1;
+    memset(u, '\0', sizeof(BNS) * ulen);
+    ulen = mp_mul(u, op1 + x, op2 + x, la1, lb1);
+
+    /* subtract t from partial result */
+    xlen = mp_sub(r, r, t, xlen, tlen);
+
+    /* subtract u form partial result */
+    xlen = mp_sub(r, r, u, xlen, ulen);
+
+    /* add ux^2 to partial result */
+
+    r = rop + (x << 1);		/* multiplied by x^2 = 2^(BNSBITS*x*2) */
+    xlen = len1 + len2;
+    xlen = mp_add(r, r, u, xlen, ulen);
+
+    /* now add t to final result */
+    xlen = mp_add(rop, rop, t, xlen, tlen);
+
+    mp_free(t);
+    mp_free(u);
+
+    if (xlen > 1 && rop[xlen - 1] == 0)
+	--xlen;
+
+    return (xlen);
+}
+
+	/* Toom method	(partially based on GMP documentation)
+	 * Evaluation at k = [ 0 1/2 1 2 oo ]
+	 * U(x) = (U2k + U1)k + U0
+	 * V(x) = (V2k + V1)k + V0
+	 * W(x) = U(x)V(x)
+	 *
+	 * Sample:
+	 *	123 * 456
+	 *
+	 *	EVALUATION:
+	 * U(0) = (1*0+2)*0+3	=> 3
+	 * U(1) = 1+(2+3*2)*2	=> 17
+	 * U(2) = 1+2+3		=> 6
+	 * U(3) = (1*2+2)*2+3	=> 11
+	 * U(4)	= 1+(2+3*0)*0	=> 1
+	 *
+	 * V(0) = (4*0+5)*0+6	=> 6
+	 * V(1) = 4+(5+6*2)*2	=> 38
+	 * V(2) = 4+5+6		=> 15
+	 * V(3) = (4*2+5)*2+6	=> 32
+	 * V(4)	= 4+(5+6*0)*0	=> 4
+	 *
+	 *	U = [ 3   17  6  11 1 ]
+	 *	V = [ 6   38 15  32 4 ]
+	 *	W = [ 18 646 90 352 4 ]
+	 *
+	 * After that, we have:
+	 *	a = 18					(w0 already known)
+	 *	b = 16w0 + 8w1 + 4w2 + 2w3 +   w4
+	 *	c =   w0 +  w1 +  w2 +  w3 +   w4
+	 *	d =   w0 + 2w1 + 4w2 + 8w3 + 16w4
+	 *	e = 4					(w4 already known)
+	 *
+	 *	INTERPOLATION:
+	 *	b = b -16a - e		(354)
+	 *	c = c - a - e		(68)
+	 *	d = d - a - 16e		(270)
+	 *
+	 *	w = (b + d) - 8c = (10w1+8w2+10w3) - (8w1+8w2+8w3) = 2w1+2w3
+	 *	w = 2c - w		(56)
+	 *	b = b/2 = 4w1+w+w3
+	 *	b = b-c = 4w1+w+w3 - w1+w2+w3 = 3w1+w2
+	 *	c = w/2			(w2 = 28)
+	 *	b = b-c = 3w1+c - c = 3w1
+	 *	b = b/3			(w1 = 27)
+	 *	d = d/2
+	 *	d = d-b-w = b+w+4w3 - b-w = 4w3
+	 *	d = d/4			(w3 = 13)
+	 *
+	 *	RESULT:
+	 *	w4*10^4 + w3*10³ + w2*10² + w1*10 + w0
+	 *	40000   + 13000   + 2800   + 270    + 18
+	 *	10 is the base where the calculation was done
+	 *
+	 *	This sample uses small numbers, so it does not show the
+	 * advantage of the method. But for example (in base 10), when squaring
+	 *	123456789012345678901234567890
+	 *	The normal method would do 30*30=900 multiplications
+	 *	Karatsuba method would do 15*15*3=675 multiplications
+	 *	Toom method would do 10*10*5=500 multiplications
+	 * Toom method has a larger overhead if compared with Karatsuba method,
+	 * due to evaluation and interpolation, so it should be used for larger
+	 * numbers, so that the computation time of evaluation/interpolation
+	 * would be smaller than the time spent using other methods.
+	 *
+	 *	Note that Karatsuba method can be seen as a special case of
+	 * Toom method, i.e:
+	 *	U1U0 * V1V0
+	 *	with k = [ 0 1 oo ]
+	 *	U = [ U0 U1+U0 U1 ]
+	 *	V = [ V0 V1+V0 V1 ]
+	 *	W = [ U0*V0 (U1+U0)*(V1+V0) (U1+V1) ]
+	 *
+	 *	w0 = U0*V0
+	 *	w = (U1+U0)*(V1+V0)
+	 *	w2 = (U1*V1)
+	 *
+	 *	w1 = w - w0 - w2
+	 * w2x² + w1x + w0
+	 *
+	 *	See Knuth's Seminumerical Algorithms for a sample implemention
+	 * using 4 stacks and k = [ 0 1 2 3 ... ], based on the size of the
+	 * input.
+	 */
+long
+mp_toom_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2)
+{
+    long size, xsize, i;
+    BNI value;				/* used in division */
+    BNS carry;
+    BNI x;				/* shift count */
+    BNI l1, l2;
+    BNI al, bl, cl, dl, el, Ul[3], Vl[3];
+    BNS *a, *b, *c, *d, *e, *U[3], *V[3];
+
+    /* x is the base i.e. 2^(BNSBITS*x) */
+    x = (len1 + len2 + 4) / 6;
+    l1 = len1 - (x << 1);	/* length of remaining piece of op1 */
+    l2 = len2 - (x << 1);	/* length of remaining piece of op2 */
+
+    /* allocate memory for storing U and V */
+    U[0] = mp_malloc(sizeof(BNS) * (x + 2));
+    V[0] = mp_malloc(sizeof(BNS) * (x + 2));
+    U[1] = mp_malloc(sizeof(BNS) * (x + 1));
+    V[1] = mp_malloc(sizeof(BNS) * (x + 1));
+    U[2] = mp_malloc(sizeof(BNS) * (x + 2));
+    V[2] = mp_malloc(sizeof(BNS) * (x + 2));
+
+	/* EVALUATE U AND V */
+
+    /* Numbers are in the format U2x²+U1x+U0 and V2x²+V1x+V0 */
+
+	/* U[0] = U2+U1*2+U0*4 */
+
+    /* store U1*2 in U[1], this value is used twice */
+    Ul[1] = mp_lshift(U[1], op1 + x, x, 1);
+
+    /* store U0*4 in U[0] */
+    Ul[0] = mp_lshift(U[0], op1, x, 2);
+    /* add U1*2 to U[0] */
+    Ul[0] = mp_add(U[0], U[0], U[1], Ul[0], Ul[1]);
+    /* add U2 to U[0] */
+    Ul[0] = mp_add(U[0], U[0], op1 + x + x, Ul[0], l1);
+
+	/* U[2] = U2*4+U1*2+U0 */
+
+    /* store U2*4 in U[2] */
+    Ul[2] = mp_lshift(U[2], op1 + x + x, l1, 2);
+    /* add U1*2 to U[2] */
+    Ul[2] = mp_add(U[2], U[2], U[1], Ul[2], Ul[1]);
+    /* add U0 to U[2] */
+    Ul[2] = mp_add(U[2], U[2], op1, Ul[2], x);
+
+	/* U[1] = U2+U1+U0 */
+
+    Ul[1] = mp_add(U[1], op1, op1 + x, x, x);
+    Ul[1] = mp_add(U[1], U[1], op1 + x + x, Ul[1], l1);
+
+
+	/* Evaluate V[x], same code as U[x] */
+    Vl[1] = mp_lshift(V[1], op2 + x, x, 1);
+    Vl[0] = mp_lshift(V[0], op2, x, 2);
+    Vl[0] = mp_add(V[0], V[0], V[1], Vl[0], Vl[1]);
+    Vl[0] = mp_add(V[0], V[0], op2 + x + x, Vl[0], l2);
+    Vl[2] = mp_lshift(V[2], op2 + x + x, l2, 2);
+    Vl[2] = mp_add(V[2], V[2], V[1], Vl[2], Vl[1]);
+    Vl[2] = mp_add(V[2], V[2], op2, Vl[2], x);
+    Vl[1] = mp_add(V[1], op2, op2 + x, x, x);
+    Vl[1] = mp_add(V[1], V[1], op2 + x + x, Vl[1], l2);
+
+
+	/* MULTIPLY U[] AND V[] */
+
+	/* calculate (U2+U1*2+U0*4) * (V2+V1*2+V0*4) */
+    b = mp_calloc(1, sizeof(BNS) * (Ul[0] * Vl[0]));
+    bl = mp_mul(b, U[0], V[0], Ul[0], Vl[0]);
+    mp_free(U[0]);
+    mp_free(V[0]);
+
+	/* calculate (U2+U1+U0) * (V2+V1+V0) */
+    c = mp_calloc(1, sizeof(BNS) * (Ul[1] * Vl[1]));
+    cl = mp_mul(c, U[1], V[1], Ul[1], Vl[1]);
+    mp_free(U[1]);
+    mp_free(V[1]);
+
+	/* calculate (U2*4+U1*2+U0) * (V2*4+V1*2+V0) */
+    d = mp_calloc(1, sizeof(BNS) * (Ul[2] * Vl[2]));
+    dl = mp_mul(d, U[2], V[2], Ul[2], Vl[2]);
+    mp_free(U[2]);
+    mp_free(V[2]);
+
+	/* calculate U0 * V0 */
+    a = mp_calloc(1, sizeof(BNS) * (x + x));
+    al = mp_mul(a, op1, op2, x, x);
+
+	/* calculate U2 * V2 */
+    e = mp_calloc(1, sizeof(BNS) * (l1 + l2));
+    el = mp_mul(e, op1 + x + x, op2 + x + x, l1, l2);
+
+
+	/* INTERPOLATE COEFFICIENTS */
+
+    /* b = b - 16a - e */
+    size = mp_lshift(rop, a, al, 4);
+    bl = mp_sub(b, b, rop, bl, size);
+    bl = mp_sub(b, b, e, bl, el);
+
+    /* c = c - a - e*/
+    cl = mp_sub(c, c, a, cl, al);
+    cl = mp_sub(c, c, e, cl, el);
+
+    /* d = d - a - 16e */
+    dl = mp_sub(d, d, a, dl, al);
+    size = mp_lshift(rop, e, el, 4);
+    dl = mp_sub(d, d, rop, dl, size);
+
+    /* w = (b + d) - 8c */
+    size = mp_add(rop, b, d, bl, dl);
+    xsize = mp_lshift(rop + size, c, cl, 3);	/* rop has enough storage */
+    size = mp_sub(rop, rop, rop + size, size, xsize);
+
+    /* w = 2c - w*/
+    xsize = mp_lshift(rop + size, c, cl, 1);
+    size = mp_sub(rop, rop + size, rop, xsize, size);
+
+    /* b = b/2 */
+    bl = mp_rshift(b, b, bl, 1);
+
+    /* b = b - c */
+    bl = mp_sub(b, b, c, bl, cl);
+
+    /* c = w / 2 */
+    cl = mp_rshift(c, rop, size, 1);
+
+    /* b = b - c */
+    bl = mp_sub(b, b, c, bl, cl);
+
+    /* b = b/3 */
+	/* maybe the most expensive calculation */
+    i = bl - 1;
+    value = b[i];
+    b[i] = value / 3;
+    for (--i; i >= 0; i--) {
+	carry = value % 3;
+	value = ((BNI)carry << BNSBITS) + b[i];
+	b[i] = (BNS)(value / 3);
+    }
+
+    /* d = d/2 */
+    dl = mp_rshift(d, d, dl, 1);
+
+    /* d = d - b - w */
+    dl = mp_sub(d, d, b, dl, bl);
+    dl = mp_sub(d, d, rop, dl, size);
+
+    /* d = d/4 */
+    dl = mp_rshift(d, d, dl, 2);
+
+
+	/* STORE RESULT IN ROP */
+    /* first clear memory used as temporary variable w and 8c */
+    memset(rop, '\0', sizeof(BNS) * (len1 + len2));
+
+    i = x * 4;
+    xsize = (len1 + len2) - i;
+    size = mp_add(rop + i, rop + i, e, xsize, el) + i;
+    i = x * 3;
+    xsize = size - i;
+    size = mp_add(rop + i, rop + i, d, xsize, dl) + i;
+    i = x * 2;
+    xsize = size - i;
+    size = mp_add(rop + i, rop + i, c, xsize, cl) + i;
+    i = x;
+    xsize = size - i;
+    size = mp_add(rop + i, rop + i, b, xsize, bl) + i;
+    size = mp_add(rop, rop, a, size, al);
+
+    mp_free(e);
+    mp_free(d);
+    mp_free(c);
+    mp_free(b);
+    mp_free(a);
+
+    if (size > 1 && rop[size - 1] == 0)
+	--size;
+
+    return (size);
+}
+
+long
+mp_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2)
+{
+    if (len1 < len2)
+	MP_SWAP(op1, op2, len1, len2);
+
+    if (len1 < KARATSUBA || len2 < KARATSUBA)
+	return (mp_base_mul(rop, op1, op2, len1, len2));
+    else if (len1 < TOOM && len2 < TOOM && len2 > ((len1 + 1) >> 1))
+	return (mp_karatsuba_mul(rop, op1, op2, len1, len2));
+    else if (len1 >= TOOM && len2 >= TOOM && (len2 + 2) / 3 == (len1 + 2) / 3)
+	return (mp_toom_mul(rop, op1, op2, len1, len2));
+    else {
+	long xsize, psize, isize;
+	BNS *ptr;
+
+	/* adjust index pointer and estimated size of result */
+	isize = 0;
+	xsize = len1 + len2;
+	mp_mul(rop, op1, op2, len2, len2);
+	/* adjust pointers */
+	len1 -= len2;
+	op1 += len2;
+
+	/* allocate buffer for intermediate multiplications */
+	if (len1 > len2)
+	    ptr = mp_calloc(1, sizeof(BNS) * (len2 + len2));
+	else
+	    ptr = mp_calloc(1, sizeof(BNS) * (len1 + len2));
+
+	/* loop multiplying len2 size operands at a time */
+	while (len1 >= len2) {
+	    isize += len2;
+	    psize = mp_mul(ptr, op1, op2, len2, len2);
+	    mp_add(rop + isize, rop + isize, ptr, xsize - isize, psize);
+	    len1 -= len2;
+	    op1 += len2;
+
+	    /* multiplication routines require zeroed memory */
+	    memset(ptr, '\0', sizeof(BNS) * (MIN(len1, len2) + len2));
+	}
+
+	/* len1 was not a multiple of len2 */
+	if (len1) {
+	    isize += len2;
+	    psize = mp_mul(ptr, op2, op1, len2, len1);
+	    mp_add(rop + isize, rop + isize, ptr, xsize, psize);
+	}
+
+	/* adjust result size */
+	if (rop[xsize - 1] == 0)
+	    --xsize;
+
+	mp_free(ptr);
+
+	return (xsize);
+    }
+}
diff --git a/lisp/mp/mp.h b/lisp/mp/mp.h
new file mode 100644
index 0000000..88f1b24
--- /dev/null
+++ b/lisp/mp/mp.h
@@ -0,0 +1,435 @@
+/*
+ * Copyright (c) 2002 by The XFree86 Project, Inc.
+ *
+ * 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 XFREE86 PROJECT 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.
+ *
+ * Except as contained in this notice, the name of the XFree86 Project shall
+ * not be used in advertising or otherwise to promote the sale, use or other
+ * dealings in this Software without prior written authorization from the
+ * XFree86 Project.
+ *
+ * Author: Paulo César Pereira de Andrade
+ */
+
+/* $XFree86: xc/programs/xedit/lisp/mp/mp.h,v 1.6 2003/01/12 03:55:51 tsi Exp $ */
+
+#include <stdio.h>
+#include <math.h>
+#ifdef sun
+#include <ieeefp.h>
+#endif
+#include <float.h>
+#include <stdlib.h>
+#include <limits.h>
+#include <ctype.h>
+#include <string.h>
+
+#ifndef __mp_h_
+#define __mp_h_
+
+#ifdef __GNUC__
+#define INLINE __inline__
+#else
+#define INLINE /**/
+#endif
+
+/* this normally is better for multiplication and also
+ * simplify addition loops putting the larger value first */
+#define MP_SWAP(op1, op2, len1, len2) {	\
+    BNS *top = op1;			\
+    BNI tlen = len1;			\
+					\
+    op1 = op2;				\
+    len1 = len2;			\
+    op2 = top;				\
+    len2 = tlen;			\
+}
+
+/*
+ * At least this length to use Karatsuba multiplication method
+ */
+#define KARATSUBA		32
+
+/*
+ * At least this length to use Toom multiplication method
+ */
+#define TOOM			128
+
+#if ULONG_MAX > 4294967295UL
+  /* sizeof(long) == 8 and sizeof(int) == 4 */
+# define BNI		unsigned long
+# define BNS		unsigned int
+# define MINSLONG	0x8000000000000000UL
+# define CARRY		0x100000000
+# define LMASK		0xffffffff00000000UL
+# define SMASK		0x00000000ffffffffUL
+# define BNIBITS	64
+# define BNSBITS	32
+# ifndef LONG64
+#  define LONG64
+# endif
+#else
+  /* sizeof(long) == 4 and sizeof(short) == 2 */
+# define BNI		unsigned long
+# define BNS		unsigned short
+# define MINSLONG	0x80000000UL
+# define CARRY		0x10000
+# define LMASK		0xffff0000UL
+# define SMASK		0x0000ffffUL
+# define BNIBITS	32
+# define BNSBITS	16
+#endif
+
+#ifdef MAX
+#undef MAX
+#endif
+#define MAX(a, b)	((a) > (b) ? (a) : (b))
+
+#ifdef MIN
+#undef MIN
+#endif
+#define MIN(a, b)	((a) < (b) ? (a) : (b))
+
+/*
+ * Types
+ */
+typedef struct _mpi {
+    unsigned int size : 31;
+    unsigned int sign : 1;
+    BNI alloc;
+    BNS *digs;		/* LSF format */
+} mpi;
+
+typedef struct _mpr {
+    mpi num;
+    mpi den;
+} mpr;
+
+typedef void *(*mp_malloc_fun)(size_t);
+typedef void *(*mp_calloc_fun)(size_t, size_t);
+typedef void *(*mp_realloc_fun)(void*, size_t);
+typedef void (*mp_free_fun)(void*);
+
+/*
+ * Prototypes
+ */
+/* GENERIC FUNCTIONS */
+	/* memory allocation wrappers */
+void *mp_malloc(size_t size);
+void *mp_calloc(size_t nmemb, size_t size);
+void *mp_realloc(void *pointer, size_t size);
+void mp_free(void *pointer);
+mp_malloc_fun mp_set_malloc(mp_malloc_fun);
+mp_calloc_fun mp_set_calloc(mp_calloc_fun);
+mp_realloc_fun mp_set_realloc(mp_realloc_fun);
+mp_free_fun mp_set_free(mp_free_fun);
+
+	/* adds op1 and op2, stores result in rop
+	 * rop must pointer to at least len1 + len2 + 1 elements
+	 * rop can be either op1 or op2 */
+long mp_add(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2);
+
+	/* subtracts op2 from op1, stores result in rop
+	 * rop must pointer to at least len1 + len2 elements
+	 * op1 must be >= op2
+	 * rop can be either op1 or op2 */
+long mp_sub(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2);
+
+	/* shift op to the left shift bits
+	 * rop must have enough storage for result
+	 * rop can be op */
+long mp_lshift(BNS *rop, BNS *op, BNI len, long shift);
+
+	/* shift op to the right shift bits
+	 * shift must be positive
+	 * rop can be op */
+long mp_rshift(BNS *rop, BNS *op, BNI len, long shift);
+
+	/* use simple generic multiplication method
+	 * rop cannot be the same as op1 or op2
+	 * rop must be zeroed
+	 * op1 can be op2 */
+long mp_base_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2);
+
+	/* use Karatsuba method
+	 * MIN(len1, len2) must be larger than (MAX(len1, len2) + 1) >> 1
+	 * MAX(len1, len2) should be at least 2
+	 * rop cannot be the same as op1 or op2
+	 * rop must be zeroed
+	 * op1 can be op2 */
+long mp_karatsuba_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2);
+
+	/* use Toom method
+	 * len1 / 3 should be equal to len2 / 3
+	 * len1 / 3 should be at least 1
+	 * rop cannot be the same as op1 or op2
+	 * rop must be zeroed
+	 * op1 can be op2 */
+long mp_toom_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2);
+
+	/* chooses the available multiplication methods based on it's input
+	 * rop must be a pointer to len1 + len2 elements
+	 * rop cannot be the same as op1 or op2
+	 * rop must be zeroed
+	 * op1 can be op2 */
+long mp_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2);
+
+/* INTEGER FUNCTIONS */
+	/* initialize op and set it to 0 */
+void mpi_init(mpi *op);
+
+	/* clear memory associated to op */
+void mpi_clear(mpi *op);
+
+	/* set rop to the value of op */
+void mpi_set(mpi *rop, mpi *op);
+
+	/* set rop to the value of si */
+void mpi_seti(mpi *rop, long si);
+
+	/* set rop to the floor(fabs(d)) */
+void mpi_setd(mpi *rop, double d);
+
+	/* initialize rop to number representation in str in the given base.
+	 * leading zeros are skipped.
+	 * if sign present, it is processed.
+	 * base must be in the range 2 to 36. */
+void mpi_setstr(mpi *rop, char *str, int base);
+
+	/* adds two mp integers */
+void mpi_add(mpi *rop, mpi *op1, mpi *op2);
+
+	/* adds op1 and op2 */
+void mpi_addi(mpi *rop, mpi *op1, long op2);
+
+	/* subtracts two mp integers */
+void mpi_sub(mpi *rop, mpi *op1, mpi *op2);
+
+	/* subtracts op2 from op1 */
+void mpi_subi(mpi *rop, mpi *op1, long op2);
+
+	/* multiply two mp integers */
+void mpi_mul(mpi *rop, mpi *op1, mpi *op2);
+
+	/* multiply op1 by op2 */
+void mpi_muli(mpi *rop, mpi *op1, long op2);
+
+	/* divides num by den and sets rop to result */
+void mpi_div(mpi *rop, mpi *num, mpi *den);
+
+	/* divides num by den and sets rop to the remainder */
+void mpi_rem(mpi *rop, mpi *num, mpi *den);
+
+	/* divides num by den, sets quotient to qrop and remainder to rrop
+	 * qrop is truncated towards zero.
+	 * qrop and rrop are optional
+	 * qrop and rrop cannot be the same variable */
+void mpi_divqr(mpi *qrop, mpi *rrop, mpi *num, mpi *den);
+
+	/* divides num by then and stores result in rop */
+void mpi_divi(mpi *rop, mpi *num, long den);
+
+	/* divides num by den and returns remainder */
+long mpi_remi(mpi *num, long den);
+
+	/* divides num by den
+	 * stores quotient in qrop and returns remainder */
+long mpi_divqri(mpi *qrop, mpi *num, long den);
+
+	/* sets rop to num modulo den */
+void mpi_mod(mpi *rop, mpi *num, mpi *den);
+
+	/* returns num modulo den */
+long mpi_modi(mpi *num, long den);
+
+	/* sets rop to the greatest common divisor of num and den
+	 * result is always positive */
+void mpi_gcd(mpi *rop, mpi *num, mpi *den);
+
+	/* sets rop to the least common multiple of num and den
+	 * result is always positive */
+void mpi_lcm(mpi *rop, mpi *num, mpi *den);
+
+	/* sets rop to op raised to exp */
+void mpi_pow(mpi *rop, mpi *op, unsigned long exp);
+
+	/* sets rop to the integer part of the nth root of op.
+	 * returns 1 if result is exact, 0 otherwise */
+int mpi_root(mpi *rop, mpi *op, unsigned long nth);
+
+	/* sets rop to the integer part of the square root of op.
+	 * returns 1 if result is exact, 0 otherwise */
+int mpi_sqrt(mpi *rop, mpi *op);
+
+	/* bit shift, left if shift positive, right if negative
+	 * a fast way to multiply and divide by powers of two */
+void mpi_ash(mpi *rop, mpi *op, long shift);
+
+	/* sets rop to op1 logand op2 */
+void mpi_and(mpi *rop, mpi *op1, mpi *op2);
+
+	/* sets rop to op1 logior op2 */
+void mpi_ior(mpi *rop, mpi *op1, mpi *op2);
+
+	/* sets rop to op1 logxor op2 */
+void mpi_xor(mpi *rop, mpi *op1, mpi *op2);
+
+	/* sets rop to one's complement of op */
+void mpi_com(mpi *rop, mpi *op);
+
+	/* sets rop to -op */
+void mpi_neg(mpi *rop, mpi *op);
+
+	/* sets rop to the absolute value of op */
+void mpi_abs(mpi *rop, mpi *op);
+
+	/* compares op1 and op2
+	 * returns >0 if op1 > op2, 0 if op1 = op2, and  <0 if op1 < op2 */
+int mpi_cmp(mpi *op1, mpi *op2);
+
+	/* mpi_cmp with a long integer operand */
+int mpi_cmpi(mpi *op1, long op2);
+
+	/* compares absolute value of op1 and op2
+	 * returns >0 if abs(op1) > abs(op2), 0 if abs(op1) = abs(op2),
+	 * and  <0 if abs(op1) < abs(op2) */
+int mpi_cmpabs(mpi *op1, mpi *op2);
+
+	/* mpi_cmpabs with a long integer operand */
+int mpi_cmpabsi(mpi *op1, long op2);
+
+	/* returns 1 if op1 > 0, 0 if op1 = 0, and  -1 if op1 < 0 */
+int mpi_sgn(mpi *op);
+
+	/* fastly swaps contents of op1 and op2 */
+void mpi_swap(mpi *op1, mpi *op2);
+
+	/* returns 1 if op fits in a signed long int, 0 otherwise */
+int mpi_fiti(mpi *op);
+
+	/* converts mp integer to long int
+	 * to know if the value will fit, call mpi_fiti */
+long mpi_geti(mpi *op);
+
+	/* convert mp integer to double */
+double mpi_getd(mpi *op);
+
+	/* returns exact number of characters to represent mp integer
+	 * in given base, excluding sign and ending null character.
+	 * base must be in the range 2 to 36 */
+unsigned long mpi_getsize(mpi *op, int base);
+
+	/* returns pointer to string with representation of mp integer
+	 * if str is not NULL, it must have enough space to store integer
+	 * representation, if NULL a newly allocated string is returned.
+	 * base must be in the range 2 to 36 */
+char *mpi_getstr(char *str, mpi *op, int base);
+
+/* RATIO FUNCTIONS */
+#define mpr_num(op)	(&((op)->num))
+#define mpr_den(op)	(&((op)->den))
+
+	/* initialize op and set it to 0/1 */
+void mpr_init(mpr *op);
+
+	/* clear memory associated to op */
+void mpr_clear(mpr *op);
+
+	/* set rop to the value of op */
+void mpr_set(mpr *rop, mpr *op);
+
+	/* set rop to num/den */
+void mpr_seti(mpr *rop, long num, long den);
+
+	/* set rop to the value of d */
+void mpr_setd(mpr *rop, double d);
+
+	/* initialize rop to number representation in str in the given base.
+	 * leading zeros are skipped.
+	 * if sign present, it is processed.
+	 * base must be in the range 2 to 36. */
+void mpr_setstr(mpr *rop, char *str, int base);
+
+	/* remove common factors of op */
+void mpr_canonicalize(mpr *op);
+
+	/* adds two mp rationals */
+void mpr_add(mpr *rop, mpr *op1, mpr *op2);
+
+	/* adds op1 and op2 */
+void mpr_addi(mpr *rop, mpr *op1, long op2);
+
+	/* subtracts two mp rationals */
+void mpr_sub(mpr *rop, mpr *op1, mpr *op2);
+
+	/* subtracts op2 from op1 */
+void mpr_subi(mpr *rop, mpr *op1, long op2);
+
+	/* multiply two mp rationals */
+void mpr_mul(mpr *rop, mpr *op1, mpr *op2);
+
+	/* multiply op1 by op2 */
+void mpr_muli(mpr *rop, mpr *op1, long op2);
+
+	/* divide two mp rationals */
+void mpr_div(mpr *rop, mpr *op1, mpr *op2);
+
+	/* divides op1 by op2 */
+void mpr_divi(mpr *rop, mpr *op1, long op2);
+
+	/* sets rop to 1/op */
+void mpr_inv(mpr *rop, mpr *op);
+
+	/* sets rop to -op */
+void mpr_neg(mpr *rop, mpr *op);
+
+	/* sets rop to the absolute value of op */
+void mpr_abs(mpr *rop, mpr *op);
+
+	/* compares op1 and op2
+	 * returns >0 if op1 > op2, 0 if op1 = op2, and  <0 if op1 < op2 */
+int mpr_cmp(mpr *op1, mpr *op2);
+
+	/* mpr_cmp with a long integer operand */
+int mpr_cmpi(mpr *op1, long op2);
+
+	/* compares absolute value of op1 and op2
+	 * returns >0 if abs(op1) > abs(op2), 0 if abs(op1) = abs(op2),
+	 * and  <0 if abs(op1) < abs(op2) */
+int mpr_cmpabs(mpr *op1, mpr *op2);
+
+	/* mpr_cmpabs with a long integer operand */
+int mpr_cmpabsi(mpr *op1, long op2);
+
+	/* fastly swaps contents of op1 and op2 */
+void mpr_swap(mpr *op1, mpr *op2);
+
+	/* returns 1 if op fits in a signed long int, 0 otherwise */
+int mpr_fiti(mpr *op);
+
+	/* convert mp rational to double */
+double mpr_getd(mpr *op);
+
+	/* returns pointer to string with representation of mp rational
+	 * if str is not NULL, it must have enough space to store rational
+	 * representation, if NULL a newly allocated string is returned.
+	 * base must be in the range 2 to 36 */
+char *mpr_getstr(char *str, mpr *op, int base);
+
+#endif /* __mp_h_ */
diff --git a/lisp/mp/mpi.c b/lisp/mp/mpi.c
new file mode 100644
index 0000000..506dc7e
--- /dev/null
+++ b/lisp/mp/mpi.c
@@ -0,0 +1,1656 @@
+/*
+ * Copyright (c) 2002 by The XFree86 Project, Inc.
+ *
+ * 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 XFREE86 PROJECT 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.
+ *
+ * Except as contained in this notice, the name of the XFree86 Project shall
+ * not be used in advertising or otherwise to promote the sale, use or other
+ * dealings in this Software without prior written authorization from the
+ * XFree86 Project.
+ *
+ * Author: Paulo César Pereira de Andrade
+ */
+
+/* $XFree86: xc/programs/xedit/lisp/mp/mpi.c,v 1.12 2002/11/20 07:44:43 paulo Exp $ */
+
+#include "mp.h"
+
+/*
+ * Prototypes
+ */
+	/* do the hard work of mpi_add and mpi_sub */
+static void mpi_addsub(mpi *rop, mpi *op1, mpi *op2, int sub);
+
+	/* logical functions implementation */
+static INLINE BNS mpi_logic(BNS op1, BNS op2, BNS op);
+static void mpi_log(mpi *rop, mpi *op1, mpi *op2,  BNS op);
+
+	/* internal mpi_seti, whithout memory allocation */
+static void _mpi_seti(mpi *rop, long si);
+
+/*
+ * Initialization
+ */
+static BNS onedig[1] = { 1 };
+static mpi mpone = { 1, 1, 0, (BNS*)&onedig };
+
+/*
+ * Implementation
+ */
+void
+mpi_init(mpi *op)
+{
+    op->sign = 0;
+    op->size = op->alloc = 1;
+    op->digs = mp_malloc(sizeof(BNS));
+    op->digs[0] = 0;
+}
+
+void
+mpi_clear(mpi *op)
+{
+    op->sign = 0;
+    op->size = op->alloc = 0;
+    mp_free(op->digs);
+}
+
+void
+mpi_set(mpi *rop, mpi *op)
+{
+    if (rop != op) {
+	if (rop->alloc < op->size) {
+	    rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op->size);
+	    rop->alloc = op->size;
+	}
+	rop->size = op->size;
+	memcpy(rop->digs, op->digs, sizeof(BNS) * op->size);
+	rop->sign = op->sign;
+    }
+}
+
+void
+mpi_seti(mpi *rop, long si)
+{
+    unsigned long ui;
+    int sign = si < 0;
+    int size;
+
+    if (si == MINSLONG) {
+	ui = MINSLONG;
+	size = 2;
+    }
+    else {
+	if (sign)
+	    ui = -si;
+	else
+	    ui = si;
+	if (ui < CARRY)
+	    size = 1;
+	else
+	    size = 2;
+    }
+
+    if (rop->alloc < size) {
+	rop->digs = mp_realloc(rop->digs, sizeof(BNS) * size);
+	rop->alloc = size;
+    }
+    rop->size = size;
+
+    /* store data in small mp integer */
+    rop->digs[0] = (BNS)ui;
+    if (size > 1)
+	rop->digs[1] = (BNS)(ui >> BNSBITS);
+    rop->size = size;
+
+    /* adjust result sign */
+    rop->sign = sign;
+}
+
+static void
+_mpi_seti(mpi *rop, long si)
+{
+    unsigned long ui;
+    int sign = si < 0;
+    int size;
+
+    if (si == MINSLONG) {
+	ui = MINSLONG;
+	size = 2;
+    }
+    else {
+	if (sign)
+	    ui = -si;
+	else
+	    ui = si;
+	if (ui < CARRY)
+	    size = 1;
+	else
+	    size = 2;
+    }
+
+    rop->digs[0] = (BNS)ui;
+    if (size > 1)
+	rop->digs[1] = (BNS)(ui >> BNSBITS);
+    rop->size = size;
+
+    rop->sign = sign;
+}
+
+void
+mpi_setd(mpi *rop, double d)
+{
+    long i;
+    double mantissa;
+    int shift, exponent;
+    BNI size;
+
+    if (isnan(d))
+	d = 0.0;
+    else if (!finite(d))
+	d = copysign(1.0, d) * DBL_MAX;
+
+    /* check if number is larger than 1 */
+    if (fabs(d) < 1.0) {
+	rop->digs[0] = 0;
+	rop->size = 1;
+	rop->sign = d < 0.0;
+
+	return;
+    }
+
+    mantissa = frexp(d, &exponent);
+    if (mantissa < 0)
+	mantissa = -mantissa;
+
+    size = (exponent + (BNSBITS - 1)) / BNSBITS;
+    shift = BNSBITS - (exponent & (BNSBITS - 1));
+
+    /* adjust amount of memory */
+    if (rop->alloc < size) {
+	rop->digs = mp_realloc(rop->digs, sizeof(BNS) * size);
+	rop->alloc = size;
+    }
+    rop->size = size;
+
+    /* adjust the exponent */
+    if (shift < BNSBITS)
+	mantissa = ldexp(mantissa, -shift);
+
+    /* convert double */
+    for (i = size - 1; i >= 0 && mantissa != 0.0; i--) {
+	mantissa = ldexp(mantissa, BNSBITS);
+	rop->digs[i] = (BNS)mantissa;
+	mantissa -= rop->digs[i];
+    }
+    for (; i >= 0; i--)
+	rop->digs[i] = 0;
+
+    /* normalize */
+    if (size > 1 && rop->digs[size - 1] == 0)
+	--rop->size;
+
+    rop->sign = d < 0.0;
+}
+
+/* how many BNS in the given base, log(base) / log(CARRY) */
+#ifdef LONG64
+static double str_bases[37] = {
+    0.0000000000000000, 0.0000000000000000, 0.0312500000000000,
+    0.0495300781475362, 0.0625000000000000, 0.0725602529652301,
+    0.0807800781475362, 0.0877298413143002, 0.0937500000000000,
+    0.0990601562950723, 0.1038102529652301, 0.1081072380824156,
+    0.1120300781475362, 0.1156387411919092, 0.1189798413143002,
+    0.1220903311127662, 0.1250000000000000, 0.1277332137890731,
+    0.1303101562950723, 0.1327477347951120, 0.1350602529652300,
+    0.1372599194618363, 0.1393572380824156, 0.1413613111267817,
+    0.1432800781475362, 0.1451205059304602, 0.1468887411919092,
+    0.1485902344426084, 0.1502298413143002, 0.1518119060977367,
+    0.1533403311127662, 0.1548186346995899, 0.1562500000000000,
+    0.1576373162299517, 0.1589832137890731, 0.1602900942795302,
+    0.1615601562950723,
+};
+#else
+static double str_bases[37] = {
+    0.0000000000000000, 0.0000000000000000, 0.0625000000000000,
+    0.0990601562950723, 0.1250000000000000, 0.1451205059304602,
+    0.1615601562950723, 0.1754596826286003, 0.1875000000000000,
+    0.1981203125901446, 0.2076205059304602, 0.2162144761648311,
+    0.2240601562950723, 0.2312774823838183, 0.2379596826286003,
+    0.2441806622255325, 0.2500000000000000, 0.2554664275781462,
+    0.2606203125901445, 0.2654954695902241, 0.2701205059304602,
+    0.2745198389236725, 0.2787144761648311, 0.2827226222535633,
+    0.2865601562950723, 0.2902410118609203, 0.2937774823838183,
+    0.2971804688852168, 0.3004596826286003, 0.3036238121954733,
+    0.3066806622255324, 0.3096372693991797, 0.3125000000000000,
+    0.3152746324599034, 0.3179664275781462, 0.3205801885590604,
+    0.3231203125901446,
+};
+#endif
+
+void
+mpi_setstr(mpi *rop, char *str, int base)
+{
+    long i;			/* counter */
+    int sign;			/* result sign */
+    BNI carry;			/* carry value */
+    BNI value;			/* temporary value */
+    BNI size;			/* size of result */
+    char *ptr;			/* end of valid input */
+
+    /* initialization */
+    sign = 0;
+    carry = 0;
+
+    /* skip leading spaces */
+    while (isspace(*str))
+	++str;
+
+    /* check if sign supplied */
+    if (*str == '-') {
+	sign = 1;
+	++str;
+    }
+    else if (*str == '+')
+	++str;
+
+    /* skip leading zeros */
+    while (*str == '0')
+	++str;
+
+    ptr = str;
+    while (*ptr) {
+	if (*ptr >= '0' && *ptr <= '9') {
+	    if (*ptr - '0' >= base)
+		break;
+	}
+	else if (*ptr >= 'A' && *ptr <= 'Z') {
+	    if (*ptr - 'A' + 10 >= base)
+		break;
+	}
+	else if (*ptr >= 'a' && *ptr <= 'z') {
+	    if (*ptr - 'a' + 10 >= base)
+		break;
+	}
+	else
+	    break;
+	++ptr;
+    }
+
+    /* resulting size */
+    size = (ptr - str) * str_bases[base] + 1;
+
+    /* make sure rop has enough storage */
+    if (rop->alloc < size) {
+	rop->digs = mp_realloc(rop->digs, size * sizeof(BNS));
+	rop->alloc = size;
+    }
+    rop->size = size;
+
+    /* initialize rop to zero */
+    memset(rop->digs, '\0', size * sizeof(BNS));
+
+    /* set result sign */
+    rop->sign = sign;
+
+    /* convert string */
+    for (; str < ptr; str++) {
+	value = *str;
+	if (islower(value))
+	    value = toupper(value);
+	value = value > '9' ? value - 'A' + 10 : value - '0';
+	value += rop->digs[0] * base;
+	carry = value >> BNSBITS;
+	rop->digs[0] = (BNS)value;
+	for (i = 1; i < size; i++) {
+	    value = (BNI)rop->digs[i] * base + carry;
+	    carry = value >> BNSBITS;
+	    rop->digs[i] = (BNS)value;
+	}
+    }
+
+    /* normalize */
+    if (rop->size > 1 && rop->digs[rop->size - 1] == 0)
+	--rop->size;
+}
+
+void
+mpi_add(mpi *rop, mpi *op1, mpi *op2)
+{
+    mpi_addsub(rop, op1, op2, 0);
+}
+
+void
+mpi_addi(mpi *rop, mpi *op1, long op2)
+{
+    BNS digs[2];
+    mpi op;
+
+    op.digs = (BNS*)digs;
+    _mpi_seti(&op, op2);
+
+    mpi_addsub(rop, op1, &op, 0);
+}
+
+void
+mpi_sub(mpi *rop, mpi *op1, mpi *op2)
+{
+    mpi_addsub(rop, op1, op2, 1);
+}
+
+void
+mpi_subi(mpi *rop, mpi *op1, long op2)
+{
+    BNS digs[2];
+    mpi op;
+
+    op.digs = (BNS*)digs;
+    _mpi_seti(&op, op2);
+
+    mpi_addsub(rop, op1, &op, 1);
+}
+
+static void
+mpi_addsub(mpi *rop, mpi *op1, mpi *op2, int sub)
+{
+    long xlen;				/* maximum result size */
+
+    if (sub ^ (op1->sign == op2->sign)) {
+	/* plus one for possible carry */
+	xlen = MAX(op1->size, op2->size) + 1;
+	if (rop->alloc < xlen) {
+	    rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xlen);
+	    rop->alloc = xlen;
+	}
+	rop->size = mp_add(rop->digs, op1->digs, op2->digs,
+			   op1->size, op2->size);
+	rop->sign = op1->sign;
+    }
+    else {
+	long cmp;			/* check for larger operator */
+
+	cmp = mpi_cmpabs(op1, op2);
+	if (cmp == 0) {
+	    rop->digs[0] = 0;
+	    rop->size = 1;
+	    rop->sign = 0;
+	}
+	else {
+	    xlen = MAX(op1->size, op2->size);
+	    if (rop->alloc < xlen) {
+		rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xlen);
+		rop->alloc = xlen;
+	    }
+	    if (cmp > 0) {
+		rop->size = mp_sub(rop->digs, op1->digs, op2->digs,
+				   op1->size, op2->size);
+		rop->sign = op1->sign;
+	    }
+	    else {
+		rop->size = mp_sub(rop->digs, op2->digs, op1->digs,
+				   op2->size, op1->size);
+		rop->sign = sub ^ op2->sign;
+	    }
+	}
+    }
+}
+
+void
+mpi_mul(mpi *rop, mpi *op1, mpi *op2)
+{
+    int sign;				/* sign flag */
+    BNS *digs;				/* result data */
+    long xsize;				/* result size */
+
+    /* get result sign */
+    sign = op1->sign ^ op2->sign;
+
+    /* check for special cases */
+    if (op1->size == 1) {
+	if (*op1->digs == 0) {
+	    /* multiply by 0 */
+	    mpi_seti(rop, 0);
+	    return;
+	}
+	else if (*op1->digs == 1) {
+	    /* multiply by +-1 */
+	    if (rop->alloc < op2->size) {
+		rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op2->size);
+		rop->alloc = op2->size;
+	    }
+	    rop->size = op2->size;
+	    memmove(rop->digs, op2->digs, sizeof(BNS) * op2->size);
+	    rop->sign = op2->size > 1 || *op2->digs ? sign : 0;
+
+	    return;
+	}
+    }
+    else if (op2->size == 1) {
+	if (*op2->digs == 0) {
+	    /* multiply by 0 */
+	    mpi_seti(rop, 0);
+	    return;
+	}
+	else if (*op2->digs == 1) {
+	    /* multiply by +-1 */
+	    if (rop->alloc < op1->size) {
+		rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op1->size);
+		rop->alloc = op1->size;
+	    }
+	    rop->size = op1->size;
+	    memmove(rop->digs, op1->digs, sizeof(BNS) * op1->size);
+	    rop->sign = op1->size > 1 || *op1->digs ? sign : 0;
+
+	    return;
+	}
+    }
+
+    /* allocate result data and set it to zero */
+    xsize = op1->size + op2->size;
+    if (rop->digs == op1->digs || rop->digs == op2->digs)
+	/* rop is also an operand */
+	digs = mp_calloc(1, sizeof(BNS) * xsize);
+    else {
+	if (rop->alloc < xsize) {
+	    rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xsize);
+	    rop->alloc = xsize;
+	}
+	digs = rop->digs;
+	memset(digs, '\0', sizeof(BNS) * xsize);
+    }
+
+    /* multiply operands */
+    xsize = mp_mul(digs, op1->digs, op2->digs, op1->size, op2->size);
+
+    /* store result in rop */
+    if (digs != rop->digs) {
+	/* if rop was an operand, free old data */
+	mp_free(rop->digs);
+	rop->digs = digs;
+    }
+    rop->size = xsize;
+
+    /* set result sign */
+    rop->sign = sign;
+}
+
+void
+mpi_muli(mpi *rop, mpi *op1, long op2)
+{
+    BNS digs[2];
+    mpi op;
+
+    op.digs = (BNS*)digs;
+    _mpi_seti(&op, op2);
+
+    mpi_mul(rop, op1, &op);
+}
+
+void
+mpi_div(mpi *rop, mpi *num, mpi *den)
+{
+    mpi_divqr(rop, NULL, num, den);
+}
+
+void
+mpi_rem(mpi *rop, mpi *num, mpi *den)
+{
+    mpi_divqr(NULL, rop, num, den);
+}
+
+/*
+ * Could/should be changed to not allocate qdigs if qrop is NULL
+ * Performance wouldn't suffer too much with a test on every loop iteration.
+ */
+void
+mpi_divqr(mpi *qrop, mpi *rrop, mpi *num, mpi *den)
+{
+    long i, j;			/* counters */
+    int qsign;			/* sign of quotient */
+    int rsign;			/* sign of remainder */
+    BNI qsize;			/* size of quotient */
+    BNI rsize;			/* size of remainder */
+    BNS qest;			/* estimative of quotient value */
+    BNS *qdigs, *rdigs;		/* work copy or result */
+    BNS *ndigs, *ddigs;		/* work copy or divisor and dividend */
+    BNI value;			/* temporary result */
+    long svalue;		/* signed temporary result (2's complement) */
+    BNS carry, scarry, denorm;	/* carry and normalization */
+    BNI dpos, npos;		/* offsets in data */
+
+    /* get signs */
+    rsign = num->sign;
+    qsign = rsign ^ den->sign;
+
+    /* check for special case */
+    if (num->size < den->size) {
+	/* quotient is zero and remainder is numerator */
+	if (rrop && rrop->digs != num->digs) {
+	    if (rrop->alloc < num->size) {
+		rrop->digs = mp_realloc(rrop->digs, sizeof(BNS) * num->size);
+		rrop->alloc = num->size;
+	    }
+	    rrop->size = num->size;
+	    memcpy(rrop->digs, num->digs, sizeof(BNS) * num->size);
+	    rrop->sign = rsign;
+	}
+	if (qrop)
+	    mpi_seti(qrop, 0);
+
+	return;
+    }
+
+    /* estimate result sizes */
+    rsize = den->size;
+    qsize = num->size - den->size + 1;
+
+    /* offsets */
+    npos = num->size - 1;
+    dpos = den->size - 1;
+
+    /* allocate space for quotient and remainder */
+    if (qrop == NULL || qrop->digs == num->digs || qrop->digs == den->digs)
+	qdigs = mp_calloc(1, sizeof(BNS) * qsize);
+    else {
+	if (qrop->alloc < qsize) {
+	    qrop->digs = mp_realloc(qrop->digs, sizeof(BNS) * qsize);
+	    qrop->alloc = qsize;
+	}
+	memset(qrop->digs, '\0', sizeof(BNS) * qsize);
+	qdigs = qrop->digs;
+    }
+    if (rrop) {
+	if (rrop->digs == num->digs || rrop->digs == den->digs)
+	    rdigs = mp_calloc(1, sizeof(BNS) * rsize);
+	else {
+	    if (rrop->alloc < rsize) {
+		rrop->digs = mp_realloc(rrop->digs, sizeof(BNS) * rsize);
+		rrop->alloc = rsize;
+	    }
+	    memset(rrop->digs, '\0', sizeof(BNS) * rsize);
+	    rdigs = rrop->digs;
+	}
+    }
+    else
+	rdigs = NULL;	/* fix gcc warning */
+
+    /* special case, only one word in divisor */
+    if (dpos == 0) {
+	for (carry = 0, i = npos; i >= 0; i--) {
+	    value = ((BNI)carry << BNSBITS) + num->digs[i];
+	    qdigs[i] = (BNS)(value / den->digs[0]);
+	    carry = (BNS)(value % den->digs[0]);
+	}
+	if (rrop)
+	    rdigs[0] = carry;
+
+	goto mpi_divqr_done;
+    }
+
+    /* make work copy of numerator */
+    ndigs = mp_malloc(sizeof(BNS) * (num->size + 1));
+    /* allocate one extra word an update offset */
+    memcpy(ndigs, num->digs, sizeof(BNS) * num->size);
+    ndigs[num->size] = 0;
+    ++npos;
+
+    /* normalize */
+    denorm = (BNS)((BNI)CARRY / ((BNI)(den->digs[dpos]) + 1));
+
+    if (denorm > 1) {
+	/* i <= num->size because ndigs has an extra word */
+	for (carry = 0, i = 0; i <= num->size; i++) {
+	    value = ndigs[i] * (BNI)denorm + carry;
+	    ndigs[i] = (BNS)value;
+	    carry = (BNS)(value >> BNSBITS);
+	}
+	/* make work copy of denominator */
+	ddigs = mp_malloc(sizeof(BNS) * den->size);
+	memcpy(ddigs, den->digs, sizeof(BNS) * den->size);
+	for (carry = 0, i = 0; i < den->size; i++) {
+	    value = ddigs[i] * (BNI)denorm + carry;
+	    ddigs[i] = (BNS)value;
+	    carry = (BNS)(value >> BNSBITS);
+	}
+    }
+    else
+	/* only allocate copy of denominator if going to change it */
+	ddigs = den->digs;
+
+    /* divide mp integers */
+    for (j = qsize - 1; j >= 0; j--, npos--) {
+	/* estimate quotient */
+	if (ndigs[npos] == ddigs[dpos])
+	    qest = (BNS)SMASK;
+	else
+	    qest = (BNS)((((BNI)(ndigs[npos]) << 16) + ndigs[npos - 1]) /
+			 ddigs[dpos]);
+
+	while ((value = ((BNI)(ndigs[npos]) << 16) + ndigs[npos - 1] -
+	        qest * (BNI)(ddigs[dpos])) < CARRY &&
+		ddigs[dpos - 1] * (BNI)qest >
+		(value << BNSBITS) + ndigs[npos - 2])
+	       --qest;
+
+	/* multiply and subtract */
+	carry = scarry = 0;
+	for (i = 0; i < den->size; i++) {
+	    value = qest * (BNI)ddigs[i] + carry;
+	    carry = (BNS)(value >> BNSBITS);
+	    svalue = (long)ndigs[npos - dpos + i - 1] - (long)(value & SMASK) -
+		     (long)scarry;
+	    ndigs[npos - dpos + i - 1] = (BNS)svalue;
+	    scarry = svalue < 0;
+	}
+
+	svalue = (long)ndigs[npos] - (long)(carry & SMASK) - (long)scarry;
+	ndigs[npos] = (BNS)svalue;
+
+	if (svalue & LMASK) {
+	    /* quotient too big */
+	    --qest;
+	    carry = 0;
+	    for (i = 0; i < den->size; i++) {
+		value = ndigs[npos - dpos + i - 1] + (BNI)carry + (BNI)ddigs[i];
+		ndigs[npos - dpos + i - 1] = (BNS)value;
+		carry = (BNS)(value >> BNSBITS);
+	    }
+	    ndigs[npos] += carry;
+	}
+
+	qdigs[j] = qest;
+    }
+
+    /* calculate remainder */
+    if (rrop) {
+	for (carry = 0, j = dpos; j >= 0; j--) {
+	    value = ((BNI)carry << BNSBITS) + ndigs[j];
+	    rdigs[j] = (BNS)(value / denorm);
+	    carry = (BNS)(value % denorm);
+	}
+    }
+
+    mp_free(ndigs);
+    if (ddigs != den->digs)
+	mp_free(ddigs);
+
+mpi_divqr_done:
+    if (rrop) {
+	if (rrop->digs != rdigs)
+	    mp_free(rrop->digs);
+	/* normalize remainder */
+	for (i = rsize - 1; i >= 0; i--)
+	    if (rdigs[i] != 0)
+		break;
+	if (i != rsize - 1) {
+	    if (i < 0) {
+		rsign = 0;
+		rsize = 1;
+	    }
+	    else
+		rsize = i + 1;
+	}
+	rrop->digs = rdigs;
+	rrop->sign = rsign;
+	rrop->size = rsize;
+    }
+
+    /* normalize quotient */
+    if (qrop) {
+	if (qrop->digs != qdigs)
+	    mp_free(qrop->digs);
+	for (i = qsize - 1; i >= 0; i--)
+	    if (qdigs[i] != 0)
+		break;
+	if (i != qsize - 1) {
+	    if (i < 0) {
+		qsign = 0;
+		qsize = 1;
+	    }
+	    else
+		qsize = i + 1;
+	}
+	qrop->digs = qdigs;
+	qrop->sign = qsign;
+	qrop->size = qsize;
+    }
+    else
+	mp_free(qdigs);
+}
+
+long
+mpi_divqri(mpi *qrop, mpi *num, long den)
+{
+    BNS ddigs[2];
+    mpi dop, rrop;
+    long remainder;
+
+    dop.digs = (BNS*)ddigs;
+    _mpi_seti(&dop, den);
+
+    memset(&rrop, '\0', sizeof(mpi));
+    mpi_init(&rrop);
+    mpi_divqr(qrop, &rrop, num, &dop);
+    remainder = rrop.digs[0];
+    if (rrop.size > 1)
+	remainder |= (BNI)(rrop.digs[1]) << BNSBITS;
+    if (rrop.sign)
+	remainder = -remainder;
+    mpi_clear(&rrop);
+
+    return (remainder);
+}
+
+void
+mpi_divi(mpi *rop, mpi *num, long den)
+{
+    BNS ddigs[2];
+    mpi dop;
+
+    dop.digs = (BNS*)ddigs;
+    _mpi_seti(&dop, den);
+
+    mpi_divqr(rop, NULL, num, &dop);
+}
+
+long
+mpi_remi(mpi *num, long den)
+{
+    return (mpi_divqri(NULL, num, den));
+}
+
+void
+mpi_mod(mpi *rop, mpi *num, mpi *den)
+{
+    mpi_rem(rop, num, den);
+    if (num->sign ^ den->sign)
+	mpi_add(rop, rop, den);
+}
+
+long
+mpi_modi(mpi *num, long den)
+{
+    long remainder;
+
+    remainder = mpi_remi(num, den);
+    if (num->sign ^ (den < 0))
+	remainder += den;
+
+    return (remainder);
+}
+
+void
+mpi_gcd(mpi *rop, mpi *num, mpi *den)
+{
+    long cmp;
+    mpi rem;
+
+    /* check if result already given */
+    cmp = mpi_cmpabs(num, den);
+
+    /* check if num is equal to den or if num is zero */
+    if (cmp == 0 || (num->size == 1 && num->digs[0] == 0)) {
+	mpi_set(rop, den);
+	rop->sign = 0;
+	return;
+    }
+    /* check if den is not zero */
+    if (den->size == 1 && den->digs[0] == 0) {
+	mpi_set(rop, num);
+	rop->sign = 0;
+	return;
+    }
+
+    /* don't call mpi_init, relies on realloc(0, size) == malloc(size) */
+    memset(&rem, '\0', sizeof(mpi));
+
+    /* if num larger than den */
+    if (cmp > 0) {
+	mpi_rem(&rem, num, den);
+	if (rem.size == 1 && rem.digs[0] == 0) {
+	    /* exact division */
+	    mpi_set(rop, den);
+	    rop->sign = 0;
+	    mpi_clear(&rem);
+	    return;
+	}
+	mpi_set(rop, den);
+    }
+    else {
+	mpi_rem(&rem, den, num);
+	if (rem.size == 1 && rem.digs[0] == 0) {
+	    /* exact division */
+	    mpi_set(rop, num);
+	    rop->sign = 0;
+	    mpi_clear(&rem);
+	    return;
+	}
+	mpi_set(rop, num);
+    }
+
+    /* loop using positive values */
+    rop->sign = rem.sign = 0;
+
+    /* cannot optimize this inverting rem/rop assignment earlier
+     * because rop mais be an operand */
+    mpi_swap(rop, &rem);
+
+    /* Euclides algorithm */
+    for (;;) {
+	mpi_rem(&rem, &rem, rop);
+	if (rem.size == 1 && rem.digs[0] == 0)
+	    break;
+	mpi_swap(rop, &rem);
+    }
+    mpi_clear(&rem);
+}
+
+void
+mpi_lcm(mpi *rop, mpi *num, mpi *den)
+{
+    mpi gcd;
+
+    /* check for zero operand */
+    if ((num->size == 1 && num->digs[0] == 0) ||
+	(den->size == 1 && den->digs[0] == 0)) {
+	rop->digs[0] = 0;
+	rop->sign = 0;
+	return;
+    }
+
+    /* don't call mpi_init, relies on realloc(0, size) == malloc(size) */
+    memset(&gcd, '\0', sizeof(mpi));
+
+    mpi_gcd(&gcd, num, den);
+    mpi_div(&gcd, den, &gcd);
+    mpi_mul(rop, &gcd, num);
+    rop->sign = 0;
+
+    mpi_clear(&gcd);
+}
+
+void
+mpi_pow(mpi *rop, mpi *op, unsigned long exp)
+{
+    mpi zop, top;
+
+    if (exp == 2) {
+	mpi_mul(rop, op, op);
+	return;
+    }
+    /* check for op**0 */
+    else if (exp == 0) {
+	rop->digs[0] = 1;
+	rop->size = 1;
+	rop->sign = 0;
+	return;
+    }
+    else if (exp == 1) {
+	mpi_set(rop, op);
+	return;
+    }
+    else if (op->size == 1) {
+	if (op->digs[0] == 0) {
+	    mpi_seti(rop, 0);
+	    return;
+	}
+	else if (op->digs[0] == 1) {
+	    mpi_seti(rop, op->sign && (exp & 1) ? -1 : 1);
+	    return;
+	}
+    }
+
+    memset(&zop, '\0', sizeof(mpi));
+    memset(&top, '\0', sizeof(mpi));
+    mpi_set(&zop, op);
+    mpi_set(&top, op);
+    for (--exp; exp; exp >>= 1) {
+	if (exp & 1)
+	    mpi_mul(&zop, &top, &zop);
+	mpi_mul(&top, &top, &top);
+    }
+
+    mpi_clear(&top);
+    rop->sign = zop.sign;
+    mp_free(rop->digs);
+    rop->digs = zop.digs;
+    rop->size = zop.size;
+}
+
+/* Find integer root of given number using the iteration
+ *	x{n+1} = ((K-1) * x{n} + N / x{n}^(K-1)) / K
+ */
+int
+mpi_root(mpi *rop, mpi *op, unsigned long nth)
+{
+    long bits, cmp;
+    int exact;
+    int sign;
+    mpi *r, t, temp, quot, old, rem;
+
+    sign = op->sign;
+
+    /* divide by zero op**1/0 */
+    if (nth == 0) {
+	int one = 1, zero = 0;
+	one = one / zero;
+    }
+    /* result is complex */
+    if (sign && !(nth & 1)) {
+	int one = 1, zero = 0;
+	one = one / zero;
+    }
+
+    /* special case op**1/1 = op */
+    if (nth == 1) {
+	mpi_set(rop, op);
+	return (1);
+    }
+
+    bits = mpi_getsize(op, 2) - 2;
+
+    if (bits < 0 || bits / nth == 0) {
+	/* integral root is surely less than 2 */
+	exact = op->size == 1 && (op->digs[0] == 1 || op->digs[0] == 0);
+	mpi_seti(rop, sign ? -1 : op->digs[0] == 0 ? 0 : 1);
+
+	return (exact == 1);
+    }
+
+    /* initialize */
+    if (rop != op)
+	r = rop;
+    else {
+	r = &t;
+	memset(r, '\0', sizeof(mpi));
+    }
+    memset(&temp, '\0', sizeof(mpi));
+    memset(&quot, '\0', sizeof(mpi));
+    memset(&old, '\0', sizeof(mpi));
+    memset(&rem, '\0', sizeof(mpi));
+
+    if (sign)
+	r->sign = 0;
+
+    /* root aproximation */
+    mpi_ash(r, op, -(bits - (bits / nth)));
+
+    for (;;) {
+	mpi_pow(&temp, r, nth - 1);
+	mpi_divqr(&quot, &rem, op, &temp);
+	cmp = mpi_cmpabs(r, &quot);
+	if (cmp == 0) {
+	    exact = mpi_cmpi(&rem, 0) == 0;
+	    break;
+	}
+	else if (cmp < 0) {
+	    if (mpi_cmpabs(r, &old) == 0) {
+		exact = 0;
+		break;
+	    }
+	    mpi_set(&old, r);
+	}
+	mpi_muli(&temp, r, nth - 1);
+	mpi_add(&quot, &quot, &temp);
+	mpi_divi(r, &quot, nth);
+    }
+
+    mpi_clear(&temp);
+    mpi_clear(&quot);
+    mpi_clear(&old);
+    mpi_clear(&rem);
+    if (r != rop) {
+	mpi_set(rop, r);
+	mpi_clear(r);
+    }
+    rop->sign = sign;
+
+    return (exact);
+}
+
+/*
+ * Find square root using the iteration:
+ *	x{n+1} = (x{n}+N/x{n})/2
+ */
+int
+mpi_sqrt(mpi *rop, mpi *op)
+{
+    long bits, cmp;
+    int exact;
+    mpi *r, t, quot, rem, old;
+
+    /* result is complex */
+    if (op->sign) {
+	int one = 1, zero = 0;
+	one = one / zero;
+    }
+
+    bits = mpi_getsize(op, 2) - 2;
+
+    if (bits < 2) {
+	/* integral root is surely less than 2 */
+	exact = op->size == 1 && (op->digs[0] == 1 || op->digs[0] == 0);
+	mpi_seti(rop, op->digs[0] == 0 ? 0 : 1);
+
+	return (exact == 1);
+    }
+
+    /* initialize */
+    if (rop != op)
+	r = rop;
+    else {
+	r = &t;
+	memset(r, '\0', sizeof(mpi));
+    }
+    memset(&quot, '\0', sizeof(mpi));
+    memset(&rem, '\0', sizeof(mpi));
+    memset(&old, '\0', sizeof(mpi));
+
+    /* root aproximation */
+    mpi_ash(r, op, -(bits - (bits / 2)));
+
+    for (;;) {
+	if (mpi_cmpabs(r, &old) == 0) {
+	    exact = 0;
+	    break;
+	}
+	mpi_divqr(&quot, &rem, op, r);
+	cmp = mpi_cmpabs(&quot, r);
+	if (cmp == 0) {
+	    exact = mpi_cmpi(&rem, 0) == 0;
+	    break;
+	}
+	else if (cmp > 0 && rem.size == 1 && rem.digs[0] == 0)
+	    mpi_subi(&quot, &quot, 1);
+	mpi_set(&old, r);
+	mpi_add(r, r, &quot);
+	mpi_ash(r, r, -1);
+    }
+    mpi_clear(&quot);
+    mpi_clear(&rem);
+    mpi_clear(&old);
+    if (r != rop) {
+	mpi_set(rop, r);
+	mpi_clear(r);
+    }
+
+    return (exact);
+}
+
+void
+mpi_ash(mpi *rop, mpi *op, long shift)
+{
+    long i;			/* counter */
+    long xsize;			/* maximum result size */
+    BNS *digs;
+
+    /* check for 0 shift, multiply/divide by 1 */
+    if (shift == 0) {
+	if (rop != op) {
+	    if (rop->alloc < op->size) {
+		rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op->size);
+		rop->alloc = op->size;
+	    }
+	    rop->size = op->size;
+	    memcpy(rop->digs, op->digs, sizeof(BNS) * op->size);
+	}
+
+	return;
+    }
+    else if (op->size == 1 && op->digs[0] == 0) {
+	rop->sign = 0;
+	rop->size = 1;
+	rop->digs[0] = 0;
+
+	return;
+    }
+
+    /* check shift and initialize */
+    if (shift > 0)
+	xsize = op->size + (shift / BNSBITS) + 1;
+    else {
+	xsize = op->size - ((-shift) / BNSBITS);
+	if (xsize <= 0) {
+	    rop->size = 1;
+	    rop->sign = op->sign;
+	    rop->digs[0] = op->sign ? 1 : 0;
+
+	    return;
+	}
+    }
+
+    /* allocate/adjust memory for result */
+    if (rop == op)
+	digs = mp_malloc(sizeof(BNS) * xsize);
+    else {
+	if (rop->alloc < xsize) {
+	    rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xsize);
+	    rop->alloc = xsize;
+	}
+	digs = rop->digs;
+    }
+
+    /* left shift, multiply by power of two */
+    if (shift > 0)
+	rop->size = mp_lshift(digs, op->digs, op->size, shift);
+    /* right shift, divide by power of two */
+    else {
+	long carry = 0;
+
+	if (op->sign) {
+	    BNI words, bits;
+
+	    words = -shift / BNSBITS;
+	    bits = -shift % BNSBITS;
+	    for (i = 0; i < words; i++)
+		carry |= op->digs[xsize + i];
+	    if (!carry) {
+		for (i = 0; i < bits; i++)
+		    if (op->digs[op->size - xsize] & (1 << i)) {
+			carry = 1;
+			break;
+		    }
+	    }
+	}
+	rop->size = mp_rshift(digs, op->digs, op->size, -shift);
+
+	if (carry)
+	    /* emulates two's complement subtracting 1 from the result */
+	    rop->size = mp_add(digs, digs, mpone.digs, rop->size, 1);
+    }
+
+    if (rop->digs != digs) {
+	mp_free(rop->digs);
+	rop->alloc = rop->size;
+	rop->digs = digs;
+    }
+    rop->sign = op->sign;
+}
+
+static INLINE BNS
+mpi_logic(BNS op1, BNS op2, BNS op)
+{
+    switch (op) {
+	case '&':
+	    return (op1 & op2);
+	case '|':
+	    return (op1 | op2);
+	case '^':
+	    return (op1 ^ op2);
+    }
+
+    return (SMASK);
+}
+
+static void
+mpi_log(mpi *rop, mpi *op1, mpi *op2, BNS op)
+{
+    long i;			/* counter */
+    long c, c1, c2;		/* carry */
+    BNS *digs, *digs1, *digs2;	/* pointers to mp data */
+    BNI size, size1, size2;
+    BNS sign, sign1, sign2;
+    BNS n, n1, n2;		/* logical operands */
+    BNI sum;
+
+    /* initialize */
+    size1 = op1->size;
+    size2 = op2->size;
+
+    sign1 = op1->sign ? SMASK : 0;
+    sign2 = op2->sign ? SMASK : 0;
+
+    sign = mpi_logic(sign1, sign2, op);
+
+    size = MAX(size1, size2);
+    if (sign)
+	++size;
+    if (rop->alloc < size) {
+	rop->digs = mp_realloc(rop->digs, sizeof(BNS) * size);
+	rop->alloc = size;
+    }
+
+    digs = rop->digs;
+    digs1 = op1->digs;
+    digs2 = op2->digs;
+
+    c = c1 = c2 = 1;
+
+    /* apply logical operation */
+    for (i = 0; i < size; i++) {
+	if (i >= size1)
+	    n1 = sign1;
+	else if (sign1) {
+	    sum = (BNI)(BNS)(~digs1[i]) + c1;
+	    c1 = (long)(sum >> BNSBITS);
+	    n1 = (BNS)sum;
+	}
+	else
+	    n1 = digs1[i];
+
+	if (i >= size2)
+	    n2 = sign2;
+	else if (sign2) {
+	    sum = (BNI)(BNS)(~digs2[i]) + c2;
+	    c2 = (long)(sum >> BNSBITS);
+	    n2 = (BNS)sum;
+	}
+	else
+	    n2 = digs2[i];
+
+	n = mpi_logic(n1, n2, op);
+	if (sign) {
+	    sum = (BNI)(BNS)(~n) + c;
+	    c = (long)(sum >> BNSBITS);
+	    digs[i] = (BNS)sum;
+	}
+	else
+	    digs[i] = n;
+    }
+
+    /* normalize */
+    for (i = size - 1; i >= 0; i--)
+	if (digs[i] != 0)
+	    break;
+    if (i != size - 1) {
+	if (i < 0) {
+	    sign = 0;
+	    size = 1;
+	}
+	else
+	    size = i + 1;
+    }
+
+    rop->sign = sign != 0;
+    rop->size = size;
+}
+
+void
+mpi_and(mpi *rop, mpi *op1, mpi *op2)
+{
+    mpi_log(rop, op1, op2, '&');
+}
+
+void
+mpi_ior(mpi *rop, mpi *op1, mpi *op2)
+{
+    mpi_log(rop, op1, op2, '|');
+}
+
+void
+mpi_xor(mpi *rop, mpi *op1, mpi *op2)
+{
+    mpi_log(rop, op1, op2, '^');
+}
+
+void
+mpi_com(mpi *rop, mpi *op)
+{
+    static BNS digs[1] = { 1 };
+    static mpi one = { 0, 1, 1, (BNS*)&digs };
+
+    mpi_log(rop, rop, &one, '^');
+}
+
+void
+mpi_neg(mpi *rop, mpi *op)
+{
+    int sign = op->sign ^ 1;
+
+    if (rop->digs != op->digs) {
+	if (rop->alloc < op->size) {
+	    rop->digs = mp_realloc(rop->digs, sizeof(BNS) * rop->size);
+	    rop->alloc = op->size;
+	}
+	rop->size = op->size;
+	memcpy(rop->digs, op->digs, sizeof(BNS) * rop->size);
+    }
+
+    rop->sign = sign;
+}
+
+void
+mpi_abs(mpi *rop, mpi *op)
+{
+    if (rop->digs != op->digs) {
+	if (rop->alloc < op->size) {
+	    rop->digs = mp_realloc(rop->digs, sizeof(BNS) * rop->size);
+	    rop->alloc = op->size;
+	}
+	rop->size = op->size;
+	memcpy(rop->digs, op->digs, sizeof(BNS) * rop->size);
+    }
+
+    rop->sign = 0;
+}
+
+int
+mpi_cmp(mpi *op1, mpi *op2)
+{
+    if (op1->sign ^ op2->sign)
+	return (op1->sign ? -1 : 1);
+
+    if (op1->size == op2->size) {
+	long i, cmp = 0;
+
+	for (i = op1->size - 1; i >= 0; i--)
+	    if ((cmp = (long)op1->digs[i] - (long)op2->digs[i]) != 0)
+		break;
+
+	return (cmp == 0 ? 0 : (cmp < 0) ^ op1->sign ? -1 : 1);
+    }
+
+    return ((op1->size < op2->size) ^ op1->sign ? -1 : 1);
+}
+
+int
+mpi_cmpi(mpi *op1, long op2)
+{
+    long cmp;
+
+    if (op1->size > 2)
+	return (op1->sign ? -1 : 1);
+
+    cmp = op1->digs[0];
+    if (op1->size == 2) {
+	cmp |= (long)op1->digs[1] << BNSBITS;
+	if (cmp == MINSLONG)
+	    return (op2 == cmp && op1->sign ? 0 : op1->sign ? -1 : 1);
+    }
+    if (op1->sign)
+	cmp = -cmp;
+
+    return (cmp - op2);
+}
+
+int
+mpi_cmpabs(mpi *op1, mpi *op2)
+{
+    if (op1->size == op2->size) {
+	long i, cmp = 0;
+
+	for (i = op1->size - 1; i >= 0; i--)
+	    if ((cmp = (long)op1->digs[i] - (long)op2->digs[i]) != 0)
+		break;
+
+	return (cmp);
+    }
+
+    return ((op1->size < op2->size) ? -1 : 1);
+}
+
+int
+mpi_cmpabsi(mpi *op1, long op2)
+{
+    unsigned long cmp;
+
+    if (op1->size > 2)
+	return (1);
+
+    cmp = op1->digs[0];
+    if (op1->size == 2)
+	cmp |= (unsigned long)op1->digs[1] << BNSBITS;
+
+    return (cmp > op2 ? 1 : cmp == op2 ? 0 : -1);
+}
+
+int
+mpi_sgn(mpi *op)
+{
+    return (op->sign ? -1 : op->size > 1 || op->digs[0] ? 1 : 0);
+}
+
+void
+mpi_swap(mpi *op1, mpi *op2)
+{
+    if (op1 != op2) {
+	mpi swap;
+
+	memcpy(&swap, op1, sizeof(mpi));
+	memcpy(op1, op2, sizeof(mpi));
+	memcpy(op2, &swap, sizeof(mpi));
+    }
+}
+
+int
+mpi_fiti(mpi *op)
+{
+    if (op->size == 1)
+	return (1);
+    else if (op->size == 2) {
+	unsigned long value = ((BNI)(op->digs[1]) << BNSBITS) | op->digs[0];
+
+	if (value & MINSLONG)
+	    return (op->sign && value == MINSLONG) ? 1 : 0;
+
+	return (1);
+    }
+
+    return (0);
+}
+
+long
+mpi_geti(mpi *op)
+{
+    long value;
+
+    value = op->digs[0];
+    if (op->size > 1)
+	value |= (BNI)(op->digs[1]) << BNSBITS;
+
+    return (op->sign && value != MINSLONG ? -value : value);
+}
+
+double
+mpi_getd(mpi *op)
+{
+    long i, len;
+    double d = 0.0;
+    int exponent;
+
+#define FLOATDIGS	sizeof(double) / sizeof(BNS)
+
+    switch (op->size) {
+	case 2:
+	    d = (BNI)(op->digs[1]) << BNSBITS;
+	case 1:
+	    d += op->digs[0];
+	    return (op->sign ? -d : d);
+	default:
+	    break;
+    }
+
+    for (i = 0, len = op->size; len > 0 && i < FLOATDIGS; i++)
+	d = ldexp(d, BNSBITS) + op->digs[--len];
+    d = frexp(d, &exponent);
+    if (len > 0)
+	exponent += len * BNSBITS;
+
+    if (d == 0.0)
+	return (d);
+
+    d = ldexp(d, exponent);
+
+    return (op->sign ? -d : d);
+}
+
+/* how many digits in a given base, floor(log(CARRY) / log(base)) */
+#ifdef LONG64
+static char dig_bases[37] = {
+     0,  0, 32, 20, 16, 13, 12, 11, 10, 10,  9,  9,  8,  8,  8,  8,
+     8,  7,  7,  7,  7,  7,  7,  7,  6,  6,  6,  6,  6,  6,  6,  6,
+     6,  6,  6,  6,  6,
+};
+#else
+static char dig_bases[37] = {
+     0,  0, 16, 10,  8,  6,  6,  5,  5,  5,  4,  4,  4,  4,  4,  4,
+     4,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,
+     3,  3,  3,  3,  3,
+};
+#endif
+
+/* how many digits per bit in a given base, log(2) / log(base) */
+static double bit_bases[37] = {
+    0.0000000000000000, 0.0000000000000000, 1.0000000000000000,
+    0.6309297535714575, 0.5000000000000000, 0.4306765580733931,
+    0.3868528072345416, 0.3562071871080222, 0.3333333333333334,
+    0.3154648767857287, 0.3010299956639811, 0.2890648263178878,
+    0.2789429456511298, 0.2702381544273197, 0.2626495350371936,
+    0.2559580248098155, 0.2500000000000000, 0.2446505421182260,
+    0.2398124665681315, 0.2354089133666382, 0.2313782131597592,
+    0.2276702486969530, 0.2242438242175754, 0.2210647294575037,
+    0.2181042919855316, 0.2153382790366965, 0.2127460535533632,
+    0.2103099178571525, 0.2080145976765095, 0.2058468324604344,
+    0.2037950470905062, 0.2018490865820999, 0.2000000000000000,
+    0.1982398631705605, 0.1965616322328226, 0.1949590218937863,
+    0.1934264036172708,
+};
+
+/* normalization base for string conversion, pow(base, dig_bases[base]) & ~CARRY */
+#ifdef LONG64
+static BNS big_bases[37] = {
+    0x00000001, 0x00000001, 0x00000000, 0xCFD41B91, 0x00000000, 0x48C27395,
+    0x81BF1000, 0x75DB9C97, 0x40000000, 0xCFD41B91, 0x3B9ACA00, 0x8C8B6D2B,
+    0x19A10000, 0x309F1021, 0x57F6C100, 0x98C29B81, 0x00000000, 0x18754571,
+    0x247DBC80, 0x3547667B, 0x4C4B4000, 0x6B5A6E1D, 0x94ACE180, 0xCAF18367,
+    0x0B640000, 0x0E8D4A51, 0x1269AE40, 0x17179149, 0x1CB91000, 0x23744899,
+    0x2B73A840, 0x34E63B41, 0x40000000, 0x4CFA3CC1, 0x5C13D840, 0x6D91B519,
+    0x81BF1000,
+};
+#else
+static BNS big_bases[37] = {
+    0x0001, 0x0001, 0x0000, 0xE6A9, 0x0000, 0x3D09, 0xB640, 0x41A7, 0x8000,
+    0xE6A9, 0x2710, 0x3931, 0x5100, 0x6F91, 0x9610, 0xC5C1, 0x0000, 0x1331,
+    0x16C8, 0x1ACB, 0x1F40, 0x242D, 0x2998, 0x2F87, 0x3600, 0x3D09, 0x44A8,
+    0x4CE3, 0x55C0, 0x5F45, 0x6978, 0x745F, 0x8000, 0x8C61, 0x9988, 0xA77B,
+    0xb640,
+};
+#endif
+
+unsigned long
+mpi_getsize(mpi *op, int base)
+{
+    unsigned long value, bits;
+
+    value = op->digs[op->size - 1];
+
+    /* count leading bits */
+    if (value) {
+	unsigned long count, carry;
+
+	for (count = 0, carry = CARRY >> 1; carry; count++, carry >>= 1)
+	    if (value & carry)
+		break;
+
+	bits = BNSBITS - count;
+    }
+    else
+	bits = 0;
+
+    return ((bits + (op->size - 1) * BNSBITS) * bit_bases[base] + 1);
+}
+
+char *
+mpi_getstr(char *str, mpi *op, int base)
+{
+    long i;			/* counter */
+    BNS *digs, *xdigs;		/* copy of op data */
+    BNI size;			/* size of op */
+    BNI digits;			/* digits per word in given base */
+    BNI bigbase;		/* big base of given base */
+    BNI strsize;		/* size of resulting string */
+    char *cp;			/* pointer in str for conversion */
+
+    /* initialize */
+    size = op->size;
+    strsize = mpi_getsize(op, base) + op->sign + 1;
+
+    if (str == NULL)
+	str = mp_malloc(strsize);
+
+    /* check for zero */
+    if (size == 1 && op->digs[0] == 0) {
+	str[0] = '0';
+	str[1] = '\0';
+
+	return (str);
+    }
+
+    digits = dig_bases[base];
+    bigbase = big_bases[base];
+
+    cp = str + strsize;
+    *--cp = '\0';
+
+    /* make copy of op data and adjust digs */
+    xdigs = mp_malloc(size * sizeof(BNS));
+    memcpy(xdigs, op->digs, size * sizeof(unsigned short));
+    digs = xdigs + size - 1;
+
+    /* convert to string */
+    for (;;) {
+	long count = -1;
+	BNI value;
+	BNS quotient, remainder = 0;
+
+	/* if power of two base */
+	if ((base & (base - 1)) == 0) {
+	    for (i = 0; i < size; i++) {
+		quotient = remainder;
+		remainder = digs[-i];
+		digs[-i] = quotient;
+		if (count < 0 && quotient)
+		    count = i;
+	    }
+	}
+	else {
+	    for (i = 0; i < size; i++) {
+		value = digs[-i] + ((BNI)remainder << BNSBITS);
+		quotient = (BNS)(value / bigbase);
+		remainder = (BNS)(value % bigbase);
+		digs[-i] = quotient;
+		if (count < 0 && quotient)
+		    count = i;
+	    }
+	}
+	quotient = remainder;
+	for (i = 0; i < digits; i++) {
+	    if (quotient == 0 && count < 0)
+		break;
+	    remainder = quotient % base;
+	    quotient /= base;
+	    *--cp = remainder < 10 ? remainder + '0' : remainder - 10 + 'A';
+	}
+	if (count < 0)
+	    break;
+	digs -= count;
+	size -= count;
+    }
+
+    /* adjust sign */
+    if (op->sign)
+	*--cp = '-';
+
+    /* remove any extra characters */
+    if (cp > str)
+	strcpy(str, cp);
+
+    mp_free(xdigs);
+
+    return (str);
+}
diff --git a/lisp/mp/mpr.c b/lisp/mp/mpr.c
new file mode 100644
index 0000000..8b26fe0
--- /dev/null
+++ b/lisp/mp/mpr.c
@@ -0,0 +1,436 @@
+/*
+ * Copyright (c) 2002 by The XFree86 Project, Inc.
+ *
+ * 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 XFREE86 PROJECT 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.
+ *
+ * Except as contained in this notice, the name of the XFree86 Project shall
+ * not be used in advertising or otherwise to promote the sale, use or other
+ * dealings in this Software without prior written authorization from the
+ * XFree86 Project.
+ *
+ * Author: Paulo César Pereira de Andrade
+ */
+
+/* $XFree86: xc/programs/xedit/lisp/mp/mpr.c,v 1.2 2002/11/08 08:01:00 paulo Exp $ */
+
+#include "mp.h"
+
+/*
+ * TODO:
+ *  Implement a fast gcd and divexact for integers, so that this code
+ * could be changed to do intermediary calculations faster using smaller
+ * numbers.
+ */
+
+/*
+ * Prototypes
+ */
+	/* do the hard work of mpr_add and mpr_sub */
+static void mpr_addsub(mpr *rop, mpr *op1, mpr *op2, int sub);
+
+	/* do the hard work of mpr_cmp and mpr_cmpabs */
+static int mpr_docmp(mpr *op1, mpr *op2, int sign);
+
+/*
+ * Implementation
+ */
+void
+mpr_init(mpr *op)
+{
+    op->num.digs = mp_malloc(sizeof(BNS));
+    op->num.sign = 0;
+    op->num.size = op->num.alloc = 1;
+    op->num.digs[0] = 0;
+
+    op->den.digs = mp_malloc(sizeof(BNS));
+    op->den.sign = 0;
+    op->den.size = op->den.alloc = 1;
+    op->den.digs[0] = 1;
+}
+
+void
+mpr_clear(mpr *op)
+{
+    op->num.sign = 0;
+    op->num.size = op->num.alloc = 0;
+    mp_free(op->num.digs);
+
+    op->den.sign = 0;
+    op->den.size = op->den.alloc = 0;
+    mp_free(op->den.digs);
+}
+
+void
+mpr_set(mpr *rop, mpr *op)
+{
+    if (rop != op) {
+	mpi_set(mpr_num(rop), mpr_num(op));
+	mpi_set(mpr_den(rop), mpr_den(op));
+    }
+}
+
+void
+mpr_seti(mpr *rop, long num, long den)
+{
+    mpi_seti(mpr_num(rop), num);
+    mpi_seti(mpr_den(rop), den);
+}
+
+void
+mpr_setd(mpr *rop, double d)
+{
+    double val, num;
+    int e, sign;
+
+    /* initialize */
+    if (d < 0) {
+	sign = 1;
+	val = -d;
+    }
+    else {
+	sign = 0;
+	val = d;
+    }
+
+    val = frexp(val, &e);
+    while (modf(val, &num) != 0.0 && val <= DBL_MAX / 2.0) {
+	--e;
+	val *= 2.0;
+    }
+    if (e >= 0) {
+	mpi_setd(mpr_num(rop), d);
+	mpi_seti(mpr_den(rop), 1);
+    }
+    else {
+	mpi_setd(mpr_num(rop), sign ? -num : num);
+	mpi_setd(mpr_den(rop), ldexp(1.0, -e));
+    }
+}
+
+void
+mpr_setstr(mpr *rop, char *str, int base)
+{
+    char *slash = strchr(str, '/');
+
+    mpi_setstr(mpr_num(rop), str, base);
+    if (slash != NULL)
+	mpi_setstr(mpr_den(rop), slash + 1, base);
+    else
+	mpi_seti(mpr_den(rop), 1);
+}
+
+void
+mpr_canonicalize(mpr *op)
+{
+    mpi gcd;
+
+    memset(&gcd, '\0', sizeof(mpi));
+
+    mpi_gcd(&gcd, mpr_num(op), mpr_den(op));
+    if (mpi_cmpabsi(&gcd, 1)) {
+	mpi_div(mpr_num(op), mpr_num(op), &gcd);
+	mpi_div(mpr_den(op), mpr_den(op), &gcd);
+    }
+
+    if (op->den.sign) {
+	op->num.sign = !op->num.sign;
+	op->den.sign = 0;
+    }
+
+    mpi_clear(&gcd);
+}
+
+void
+mpr_add(mpr *rop, mpr *op1, mpr *op2)
+{
+    mpr_addsub(rop, op1, op2, 0);
+}
+
+void
+mpr_addi(mpr *rop, mpr *op1, long op2)
+{
+    mpi prod;
+
+    memset(&prod, '\0', sizeof(mpi));
+
+    mpi_muli(&prod, mpr_den(op1), op2);
+    mpi_add(mpr_num(rop), mpr_num(op1), &prod);
+    mpi_clear(&prod);
+}
+
+void
+mpr_sub(mpr *rop, mpr *op1, mpr *op2)
+{
+    mpr_addsub(rop, op1, op2, 1);
+}
+
+void
+mpr_subi(mpr *rop, mpr *op1, long op2)
+{
+    mpi prod;
+
+    memset(&prod, '\0', sizeof(mpi));
+
+    mpi_muli(&prod, mpr_den(op1), op2);
+    mpi_sub(mpr_num(rop), mpr_num(op1), &prod);
+    mpi_clear(&prod);
+}
+
+static void
+mpr_addsub(mpr *rop, mpr *op1, mpr *op2, int sub)
+{
+    mpi prod1, prod2;
+
+    memset(&prod1, '\0', sizeof(mpi));
+    memset(&prod2, '\0', sizeof(mpi));
+
+    mpi_mul(&prod1, mpr_num(op1), mpr_den(op2));
+    mpi_mul(&prod2, mpr_num(op2), mpr_den(op1));
+
+    if (sub)
+	mpi_sub(mpr_num(rop), &prod1, &prod2);
+    else
+	mpi_add(mpr_num(rop), &prod1, &prod2);
+
+    mpi_clear(&prod1);
+    mpi_clear(&prod2);
+
+    mpi_mul(mpr_den(rop), mpr_den(op1), mpr_den(op2));
+}
+
+void
+mpr_mul(mpr *rop, mpr *op1, mpr *op2)
+{
+    /* check if temporary storage is required */
+    if (op1 == op2 && rop == op1) {
+	mpi prod;
+
+	memset(&prod, '\0', sizeof(mpi));
+
+	mpi_mul(&prod, mpr_num(op1), mpr_num(op2));
+	mpi_mul(mpr_den(rop), mpr_den(op1), mpr_den(op2));
+	mpi_set(mpr_num(rop), &prod);
+
+	mpi_clear(&prod);
+    }
+    else {
+	mpi_mul(mpr_num(rop), mpr_num(op1), mpr_num(op2));
+	mpi_mul(mpr_den(rop), mpr_den(op1), mpr_den(op2));
+    }
+}
+
+void
+mpr_muli(mpr *rop, mpr *op1, long op2)
+{
+    mpi_muli(mpr_num(rop), mpr_num(op1), op2);
+}
+
+void
+mpr_div(mpr *rop, mpr *op1, mpr *op2)
+{
+    /* check if temporary storage is required */
+    if (op1 == op2 && rop == op1) {
+	mpi prod;
+
+	memset(&prod, '\0', sizeof(mpi));
+
+	mpi_mul(&prod, mpr_num(op1), mpr_den(op2));
+	mpi_mul(mpr_den(rop), mpr_num(op2), mpr_den(op1));
+	mpi_set(mpr_num(rop), &prod);
+
+	mpi_clear(&prod);
+    }
+    else {
+	mpi_mul(mpr_num(rop), mpr_num(op1), mpr_den(op2));
+	mpi_mul(mpr_den(rop), mpr_num(op2), mpr_den(op1));
+    }
+}
+
+void
+mpr_divi(mpr *rop, mpr *op1, long op2)
+{
+    mpi_muli(mpr_den(rop), mpr_den(op1), op2);
+}
+
+void
+mpr_inv(mpr *rop, mpr *op)
+{
+    if (rop == op)
+	mpi_swap(mpr_num(op), mpr_den(op));
+    else {
+	mpi_set(mpr_num(rop), mpr_den(op));
+	mpi_set(mpr_den(rop), mpr_num(op));
+    }
+}
+
+void
+mpr_neg(mpr *rop, mpr *op)
+{
+    mpi_neg(mpr_num(rop), mpr_num(op));
+    mpi_set(mpr_den(rop), mpr_den(op));
+}
+
+void
+mpr_abs(mpr *rop, mpr *op)
+{
+    if (mpr_num(op)->sign)
+	mpi_neg(mpr_num(rop), mpr_num(op));
+    else
+	mpi_set(mpr_num(rop), mpr_num(op));
+
+    /* op may not be canonicalized */
+    if (mpr_den(op)->sign)
+	mpi_neg(mpr_den(rop), mpr_den(op));
+    else
+	mpi_set(mpr_den(rop), mpr_den(op));
+}
+
+int
+mpr_cmp(mpr *op1, mpr *op2)
+{
+    return (mpr_docmp(op1, op2, 1));
+}
+
+int
+mpr_cmpi(mpr *op1, long op2)
+{
+    int cmp;
+    mpr rat;
+
+    mpr_init(&rat);
+    mpi_seti(mpr_num(&rat), op2);
+    cmp = mpr_docmp(op1, &rat, 1);
+    mpr_clear(&rat);
+
+    return (cmp);
+}
+
+int
+mpr_cmpabs(mpr *op1, mpr *op2)
+{
+    return (mpr_docmp(op1, op2, 0));
+}
+
+int
+mpr_cmpabsi(mpr *op1, long op2)
+{
+    int cmp;
+    mpr rat;
+
+    mpr_init(&rat);
+    mpi_seti(mpr_num(&rat), op2);
+    cmp = mpr_docmp(op1, &rat, 0);
+    mpr_clear(&rat);
+
+    return (cmp);
+}
+
+static int
+mpr_docmp(mpr *op1, mpr *op2, int sign)
+{
+    int cmp, neg;
+    mpi prod1, prod2;
+
+    neg = 0;
+    if (sign) {
+	/* if op1 is negative */
+	if (mpr_num(op1)->sign ^ mpr_den(op1)->sign) {
+	    /* if op2 is positive */
+	    if (!(mpr_num(op2)->sign ^ mpr_den(op2)->sign))
+		return (-1);
+	    else
+		neg = 1;
+	}
+	/* if op2 is negative */
+	else if (mpr_num(op2)->sign ^ mpr_den(op2)->sign)
+	    return (1);
+	/* else same sign */
+    }
+
+    /* if denominators are equal, compare numerators */
+    if (mpi_cmpabs(mpr_den(op1), mpr_den(op2)) == 0) {
+	cmp = mpi_cmpabs(mpr_num(op1), mpr_num(op2));
+	if (cmp == 0)
+	    return (0);
+	if (sign && neg)
+	    return (cmp < 0 ? 1 : -1);
+	return (cmp);
+    }
+
+    memset(&prod1, '\0', sizeof(mpi));
+    memset(&prod2, '\0', sizeof(mpi));
+
+    /* "divide" op1 by op2
+     * if result is smaller than 1, op1 is smaller than op2 */
+    mpi_mul(&prod1, mpr_num(op1), mpr_den(op2));
+    mpi_mul(&prod2, mpr_num(op2), mpr_den(op1));
+
+    cmp = mpi_cmpabs(&prod1, &prod2);
+
+    mpi_clear(&prod1);
+    mpi_clear(&prod2);
+
+    if (sign && neg)
+	return (cmp < 0 ? 1 : -1);
+    return (cmp);
+}
+
+void
+mpr_swap(mpr *op1, mpr *op2)
+{
+    if (op1 != op2) {
+	mpr swap;
+
+	memcpy(&swap, op1, sizeof(mpr));
+	memcpy(op1, op2, sizeof(mpr));
+	memcpy(op2, &swap, sizeof(mpr));
+    }
+}
+
+int
+mpr_fiti(mpr *op)
+{
+    return (mpi_fiti(mpr_num(op)) && mpi_fiti(mpr_den(op)));
+}
+
+double
+mpr_getd(mpr *op)
+{
+    return (mpi_getd(mpr_num(op)) / mpi_getd(mpr_den(op)));
+}
+
+char *
+mpr_getstr(char *str, mpr *op, int base)
+{
+    int len;
+
+    if (str == NULL) {
+	len = mpi_getsize(mpr_num(op), base) + mpr_num(op)->sign + 1 +
+	      mpi_getsize(mpr_den(op), base) + mpr_den(op)->sign + 1;
+
+	str = mp_malloc(len);
+    }
+
+    (void)mpi_getstr(str, mpr_num(op), base);
+    len = strlen(str);
+    str[len] = '/';
+    (void)mpi_getstr(str + len + 1, mpr_den(op), base);
+
+    return (str);
+}