Import sources from 2011-03-30 iso image

59 files changed, 3865 insertions, 0 deletions

diff --git a/sys/src/libmp/386/mkfile b/sys/src/libmp/386/mkfile
new file mode 100755
index 000000000..616e94614
--- /dev/null
+++ b/sys/src/libmp/386/mkfile
@@ -0,0 +1,20 @@
+objtype=386
+</386/mkfile
+
+LIB=/$objtype/lib/libmp.a
+SFILES=\
+	mpvecadd.s\
+	mpvecdigmuladd.s\
+	mpvecdigmulsub.s\
+	mpvecsub.s\
+	mpdigdiv.s\
+
+HFILES=/$objtype/include/u.h /sys/include/mp.h ../port/dat.h
+
+OFILES=${SFILES:%.s=%.$O}
+
+UPDATE=mkfile\
+	$HFILES\
+	$SFILES\
+
+</sys/src/cmd/mksyslib
diff --git a/sys/src/libmp/386/mpdigdiv.s b/sys/src/libmp/386/mpdigdiv.s
new file mode 100755
index 000000000..8ee4d67a3
--- /dev/null
+++ b/sys/src/libmp/386/mpdigdiv.s
@@ -0,0 +1,21 @@
+TEXT	mpdigdiv(SB),$0
+
+	MOVL	dividend+0(FP),BX
+	MOVL	0(BX),AX
+	MOVL	4(BX),DX
+	MOVL	divisor+4(FP),BX
+	MOVL	quotient+8(FP),BP
+	XORL	CX,CX
+	CMPL	DX,BX		/* dividend >= 2^32 * divisor */
+	JHS	_divovfl
+	CMPL	BX,CX		/* divisor == 0 */
+	JE	_divovfl
+	DIVL	BX		/* AX = DX:AX/BX */
+	MOVL	AX,0(BP)
+	RET
+
+	/* return all 1's */
+_divovfl:
+	NOTL	CX
+	MOVL	CX,0(BP)
+	RET
diff --git a/sys/src/libmp/386/mpvecadd.s b/sys/src/libmp/386/mpvecadd.s
new file mode 100755
index 000000000..50a45b140
--- /dev/null
+++ b/sys/src/libmp/386/mpvecadd.s
@@ -0,0 +1,53 @@
+/*
+ *	mpvecadd(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *sum)
+ *
+ *		sum[0:alen] = a[0:alen-1] + b[0:blen-1]
+ *
+ *	prereq: alen >= blen, sum has room for alen+1 digits
+ */
+TEXT	mpvecadd(SB),$0
+
+	MOVL	alen+4(FP),DX
+	MOVL	blen+12(FP),CX
+	MOVL	a+0(FP),SI
+	MOVL	b+8(FP),BX
+	SUBL	CX,DX
+	MOVL	sum+16(FP),DI
+	XORL	BP,BP			/* this also sets carry to 0 */
+
+	/* skip addition if b is zero */
+	TESTL	CX,CX
+	JZ	_add1
+
+	/* sum[0:blen-1],carry = a[0:blen-1] + b[0:blen-1] */
+_addloop1:
+	MOVL	(SI)(BP*4), AX
+	ADCL	(BX)(BP*4), AX
+	MOVL	AX,(DI)(BP*4)
+	INCL	BP
+	LOOP	_addloop1
+
+_add1:
+	/* jump if alen > blen */
+	INCL	DX
+	MOVL	DX,CX
+	LOOP	_addloop2
+
+	/* sum[alen] = carry */
+_addend:
+	JC	_addcarry
+	MOVL	$0,(DI)(BP*4)
+	RET
+_addcarry:
+	MOVL	$1,(DI)(BP*4)
+	RET
+
+	/* sum[blen:alen-1],carry = a[blen:alen-1] + 0 */
+_addloop2:
+	MOVL	(SI)(BP*4),AX
+	ADCL	$0,AX
+	MOVL	AX,(DI)(BP*4)
+	INCL	BP
+	LOOP	_addloop2
+	JMP	_addend
+
diff --git a/sys/src/libmp/386/mpvecdigmuladd.s b/sys/src/libmp/386/mpvecdigmuladd.s
new file mode 100755
index 000000000..5a262ce8b
--- /dev/null
+++ b/sys/src/libmp/386/mpvecdigmuladd.s
@@ -0,0 +1,52 @@
+/*
+ *	mpvecdigmul(mpdigit *b, int n, mpdigit m, mpdigit *p)
+ *
+ *	p += b*m
+ *
+ *	each step look like:
+ *		hi,lo = m*b[i]
+ *		lo += oldhi + carry
+ *		hi += carry
+ *		p[i] += lo
+ *		oldhi = hi
+ *
+ *	the registers are:
+ *		hi = DX		- constrained by hardware
+ *		lo = AX		- constrained by hardware
+ *		b+n = SI	- can't be BP
+ *		p+n = DI	- can't be BP
+ *		i-n = BP
+ *		m = BX
+ *		oldhi = CX
+ *		
+ */
+TEXT	mpvecdigmuladd(SB),$0
+
+	MOVL	b+0(FP),SI
+	MOVL	n+4(FP),CX
+	MOVL	m+8(FP),BX
+	MOVL	p+12(FP),DI
+	MOVL	CX,BP
+	NEGL	BP		/* BP = -n */
+	SHLL	$2,CX
+	ADDL	CX,SI		/* SI = b + n */
+	ADDL	CX,DI		/* DI = p + n */
+	XORL	CX,CX
+_muladdloop:
+	MOVL	(SI)(BP*4),AX	/* lo = b[i] */
+	MULL	BX		/* hi, lo = b[i] * m */
+	ADDL	CX,AX		/* lo += oldhi */
+	JCC	_muladdnocarry1
+	INCL	DX		/* hi += carry */
+_muladdnocarry1:
+	ADDL	AX,(DI)(BP*4)	/* p[i] += lo */
+	JCC	_muladdnocarry2
+	INCL	DX		/* hi += carry */
+_muladdnocarry2:
+	MOVL	DX,CX		/* oldhi = hi */
+	INCL	BP		/* i++ */
+	JNZ	_muladdloop
+	XORL	AX,AX
+	ADDL	CX,(DI)(BP*4)	/* p[n] + oldhi */
+	ADCL	AX,AX		/* return carry out of p[n] */
+	RET
diff --git a/sys/src/libmp/386/mpvecdigmulsub.s b/sys/src/libmp/386/mpvecdigmulsub.s
new file mode 100755
index 000000000..04cbfad30
--- /dev/null
+++ b/sys/src/libmp/386/mpvecdigmulsub.s
@@ -0,0 +1,53 @@
+/*
+ *	mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p)
+ *
+ *	p -= b*m
+ *
+ *	each step look like:
+ *		hi,lo = m*b[i]
+ *		lo += oldhi + carry
+ *		hi += carry
+ *		p[i] += lo
+ *		oldhi = hi
+ *
+ *	the registers are:
+ *		hi = DX		- constrained by hardware
+ *		lo = AX		- constrained by hardware
+ *		b = SI		- can't be BP
+ *		p = DI		- can't be BP
+ *		i = BP
+ *		n = CX		- constrained by LOOP instr
+ *		m = BX
+ *		oldhi = EX
+ *		
+ */
+TEXT	mpvecdigmulsub(SB),$4
+
+	MOVL	b+0(FP),SI
+	MOVL	n+4(FP),CX
+	MOVL	m+8(FP),BX
+	MOVL	p+12(FP),DI
+	XORL	BP,BP
+	MOVL	BP,0(SP)
+_mulsubloop:
+	MOVL	(SI)(BP*4),AX		/* lo = b[i] */
+	MULL	BX			/* hi, lo = b[i] * m */
+	ADDL	0(SP),AX		/* lo += oldhi */
+	JCC	_mulsubnocarry1
+	INCL	DX			/* hi += carry */
+_mulsubnocarry1:
+	SUBL	AX,(DI)(BP*4)
+	JCC	_mulsubnocarry2
+	INCL	DX			/* hi += carry */
+_mulsubnocarry2:
+	MOVL	DX,0(SP)
+	INCL	BP
+	LOOP	_mulsubloop
+	MOVL	0(SP),AX
+	SUBL	AX,(DI)(BP*4)
+	JCC	_mulsubnocarry3
+	MOVL	$-1,AX
+	RET
+_mulsubnocarry3:
+	MOVL	$1,AX
+	RET
diff --git a/sys/src/libmp/386/mpvecsub.s b/sys/src/libmp/386/mpvecsub.s
new file mode 100755
index 000000000..ebe8d29a9
--- /dev/null
+++ b/sys/src/libmp/386/mpvecsub.s
@@ -0,0 +1,44 @@
+/*
+ *	mpvecsub(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *diff)
+ *
+ *		diff[0:alen-1] = a[0:alen-1] - b[0:blen-1]
+ *
+ *	prereq: alen >= blen, diff has room for alen digits
+ */
+TEXT	mpvecsub(SB),$0
+
+	MOVL	a+0(FP),SI
+	MOVL	b+8(FP),BX
+	MOVL	alen+4(FP),DX
+	MOVL	blen+12(FP),CX
+	MOVL	diff+16(FP),DI
+	SUBL	CX,DX
+	XORL	BP,BP			/* this also sets carry to 0 */
+
+	/* skip subraction if b is zero */
+	TESTL	CX,CX
+	JZ	_sub1
+
+	/* diff[0:blen-1],borrow = a[0:blen-1] - b[0:blen-1] */
+_subloop1:
+	MOVL	(SI)(BP*4),AX
+	SBBL	(BX)(BP*4),AX
+	MOVL	AX,(DI)(BP*4)
+	INCL	BP
+	LOOP	_subloop1
+
+_sub1:
+	INCL	DX
+	MOVL	DX,CX
+	LOOP	_subloop2
+	RET
+
+	/* diff[blen:alen-1] = a[blen:alen-1] - 0 */
+_subloop2:
+	MOVL	(SI)(BP*4),AX
+	SBBL	$0,AX
+	MOVL	AX,(DI)(BP*4)
+	INCL	BP
+	LOOP	_subloop2
+	RET
+
diff --git a/sys/src/libmp/alpha/mkfile b/sys/src/libmp/alpha/mkfile
new file mode 100755
index 000000000..d2117d321
--- /dev/null
+++ b/sys/src/libmp/alpha/mkfile
@@ -0,0 +1,13 @@
+objtype=alpha
+</alpha/mkfile
+
+LIB=/$objtype/lib/libmp.a
+OFILES=\
+
+HFILES=/$objtype/include/u.h /sys/include/mp.h ../port/dat.h
+
+UPDATE=\
+	mkfile\
+	$HFILES\
+
+</sys/src/cmd/mksyslib
diff --git a/sys/src/libmp/arm/mkfile b/sys/src/libmp/arm/mkfile
new file mode 100755
index 000000000..e90c8b6ab
--- /dev/null
+++ b/sys/src/libmp/arm/mkfile
@@ -0,0 +1,13 @@
+objtype=arm
+</$objtype/mkfile
+
+LIB=/$objtype/lib/libmp.a
+OFILES=	\
+
+HFILES=/$objtype/include/u.h /sys/include/mp.h ../port/dat.h
+
+UPDATE=\
+	mkfile\
+	$HFILES\
+
+</sys/src/cmd/mksyslib
diff --git a/sys/src/libmp/bigtest.c b/sys/src/libmp/bigtest.c
new file mode 100755
index 000000000..4650916e1
--- /dev/null
+++ b/sys/src/libmp/bigtest.c
@@ -0,0 +1,103 @@
+#include <u.h>
+#include <libc.h>
+#include <mp.h>
+#include <libsec.h>
+
+char *sfactors[] =
+{	"3", "5", "17", "257", "641", "65537", "274177", "2424833", "6700417", "45592577",
+    "6487031809", "67280421310721", "1238926361552897", "59649589127497217",
+    "5704689200685129054721", "4659775785220018543264560743076778192897",
+    "7455602825647884208337395736200454918783366342657",
+    "93461639715357977769163558199606896584051237541638188580280321",
+
+    "741640062627530801524787141901937474059940781097519023905821316144415759504705008092818711693940737",
+
+    "130439874405488189727484768796509903946608530841611892186895295776832416251471863574140227977573104895898783928842923844831149032913798729088601617946094119449010595906710130531906171018354491609619193912488538116080712299672322806217820753127014424577"
+};
+
+long start;
+
+void
+printmp(mpint *b, char *tag)
+{
+	int n;
+	char *p;
+
+	print("%s (%d) ", tag, b->top);
+	p = mptoa(b, 10, nil, 0);
+	write(1, p, strlen(p));
+	free(p);
+	print("\n");
+}
+
+int
+timing(void)
+{
+	long now, span;
+
+	now = time(0);
+	span = now-start;
+	start = now;
+
+	return span;
+}
+
+int expdebug;
+
+
+void
+main(int argc, char **argv)
+{
+	mpint *p, *k, *d, *b, *e, *x, *r;
+	int i;
+
+	start = time(0);
+	fmtinstall('B', mpconv);
+	mpsetminbits(2*Dbits);
+
+	x = mpnew(0);
+	e = mpnew(0);
+	r = mpnew(0);
+	p = mpnew(0);
+
+	// b = 2^32
+	b = mpcopy(mpone);
+	mpleft(b, 32, b);
+
+	// 2^29440
+	p = mpcopy(mpone);
+	mpleft(p, 29440, p);
+	// 2^27392
+	k = mpcopy(mpone);
+	mpleft(k, 27392, k);
+	// k = 2^29440 - 2^27392
+	mpsub(p, k, k);
+	// p = 2^29440 - 2^27392 + 1
+	mpadd(k, mpone, p);
+
+//	if(!probably_prime(p, 18)){
+//		print("not a prime\n");
+//		exits(0);
+//	}
+//	print("probably prime\n");
+
+	mpright(k, 10, k);
+	printmp(k, "k =");
+
+expdebug = 1;
+	mpexp(b, k, p, x);
+	printmp(x, "x =");
+	print("timing %d\n", timing());
+
+	for(i = 0; i < nelem(sfactors); i++){
+		d = strtomp(sfactors[i], nil, 10, nil);
+		// e = k/d
+		mpdiv(k, d, e, r);
+		printmp(r, "r =");
+
+		// x = b^e mod p
+		mpexp(b, e, p, x);
+		printmp(x, "x =");
+		print("timing %d\n", timing());
+	}
+}
diff --git a/sys/src/libmp/mips/mkfile b/sys/src/libmp/mips/mkfile
new file mode 100755
index 000000000..37ea56384
--- /dev/null
+++ b/sys/src/libmp/mips/mkfile
@@ -0,0 +1,20 @@
+objtype=mips
+</mips/mkfile
+
+LIB=/$objtype/lib/libmp.a
+SFILES=\
+	mpvecadd.s\
+	mpvecsub.s\
+	mpvecdigmuladd.s\
+	mpvecdigmulsub.s\
+	mpdigdiv.s\
+
+HFILES=/$objtype/include/u.h /sys/include/mp.h ../port/dat.h
+
+OFILES=${SFILES:%.s=%.$O}
+
+UPDATE=mkfile\
+	$HFILES\
+	$SFILES\
+
+</sys/src/cmd/mksyslib
diff --git a/sys/src/libmp/mips/mpdigdiv.s b/sys/src/libmp/mips/mpdigdiv.s
new file mode 100755
index 000000000..b7fd330e7
--- /dev/null
+++ b/sys/src/libmp/mips/mpdigdiv.s
@@ -0,0 +1,41 @@
+/*
+ *  This only works on R[45]000 chips that allow 64 bit
+ *  integer arithmetic even when uding 32 bit addresses
+ *
+ *	R1 = dividend*
+ *	R2 = dividend[low]
+ *	R3 = dividend[high]
+ *	R4 = 32 bit divisor
+ *	R5 = quotient*
+ */
+TEXT	mpdigdiv(SB),$0
+
+	MOVW	0(R1),R2
+	MOVW	4(R1),R3
+	MOVW	divisor+4(FP),R4
+	MOVW	quotient+8(FP),R5
+
+	/* divisor == 0 */
+	BEQ	R4,_digovfl
+
+	/* dividend >= 2^32 * divisor */
+	SGTU	R4,R3,R7
+	BEQ	R7,_digovfl
+
+_digdiv1:
+	SLLV	$32,R2
+	SLLV	$32,R3
+	SRLV	$32,R2
+	ADDVU	R2,R3
+	SLLV	$32,R4
+	SRLV	$32,R4
+	DIVVU	R4,R3
+	MOVW	LO,R1
+	MOVW	R1,0(R5)
+	RET
+
+_digovfl:
+	MOVW	$-1,R1
+	MOVW	R1,0(R5)
+	RET
+
diff --git a/sys/src/libmp/mips/mpvecadd.s b/sys/src/libmp/mips/mpvecadd.s
new file mode 100755
index 000000000..8d3e8d639
--- /dev/null
+++ b/sys/src/libmp/mips/mpvecadd.s
@@ -0,0 +1,67 @@
+#define	BDNZ	BC	16,0,
+#define	BDNE	BC	0,2,
+
+/*
+ *	mpvecadd(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *sum)
+ *
+ *		sum[0:alen] = a[0:alen-1] + b[0:blen-1]
+ *
+ *	prereq: alen >= blen, sum has room for alen+1 digits
+ *
+ *		R1 == a	(first arg passed in R1)
+ *		R3 == carry
+ *		R4 == alen
+ *		R5 == b
+ *		R6 == blen
+ *		R7 == sum
+ *		R2 == temporary
+ *		R8 == temporary
+ *		R9 == temporary
+ */
+TEXT	mpvecadd(SB),$-4
+
+	MOVW	alen+4(FP), R4
+	MOVW	b+8(FP), R5
+	MOVW	blen+12(FP), R6
+	MOVW	sum+16(FP), R7
+	SUBU	R6, R4		/* calculate counter for second loop (alen > blen) */
+	MOVW	R0, R3
+
+	/* if blen == 0, don't need to add it in */
+	BEQ	R6,_add1
+
+	/* sum[0:blen-1],carry = a[0:blen-1] + b[0:blen-1] */
+_addloop1:
+	MOVW	0(R1), R8
+	ADDU	$4, R1
+	MOVW	0(R5), R9
+	ADDU	$4, R5
+	ADDU	R3, R8
+	SGTU	R3, R8, R3
+	ADDU	R8, R9
+	SGTU	R8, R9, R2
+	ADDU	R2, R3
+	MOVW	R9, 0(R7)
+	ADDU	$4, R7
+	SUBU	$1, R6
+	BNE	R6, _addloop1
+
+_add1:
+	/* if alen == blen, we're done */
+	BEQ	R4, _addend
+
+	/* sum[blen:alen-1],carry = a[blen:alen-1] + 0 + carry */
+_addloop2:
+	MOVW	0(R1), R8
+	ADDU	$4, R1
+	ADDU	R3, R8
+	SGTU	R3, R8, R3
+	MOVW	R8, 0(R7)
+	ADDU	$4, R7
+	SUBU	$1, R4
+	BNE	R4, _addloop2
+
+	/* sum[alen] = carry */
+_addend:
+	MOVW	R3, 0(R7)
+	RET
diff --git a/sys/src/libmp/mips/mpvecdigmuladd.s b/sys/src/libmp/mips/mpvecdigmuladd.s
new file mode 100755
index 000000000..d80073b14
--- /dev/null
+++ b/sys/src/libmp/mips/mpvecdigmuladd.s
@@ -0,0 +1,58 @@
+/*
+ *	mpvecdigmuladd(mpdigit *b, int n, mpdigit m, mpdigit *p)
+ *
+ *	p += b*m
+ *
+ *	each step looks like:
+ *		hi,lo = m*b[i]
+ *		lo += oldhi + carry
+ *		hi += carry
+ *		p[i] += lo
+ *		oldhi = hi
+ *
+ *	the registers are:
+ *		b = R1
+ *		n = R4
+ *		m = R5
+ *		p = R6
+ *		i = R7
+ *		hi = R8		- constrained by hardware
+ *		lo = R9		- constrained by hardware
+ *		oldhi = R10
+ *		tmp = R11
+ *		
+ */
+TEXT	mpvecdigmuladd(SB),$0
+
+	MOVW	n+4(FP),R4
+	MOVW	m+8(FP),R5
+	MOVW	p+12(FP),R6
+
+
+	MOVW	R0, R10		/* oldhi = 0 */
+	BEQ	R6, _muladd1
+_muladdloop:
+	MOVW	0(R1), R9	/* lo = b[i] */
+	ADDU	$4, R1
+	MOVW	0(R6), R11	/* tmp = p[i] */
+	MULU	R9, R5
+	MOVW	HI, R8		/* hi = (b[i] * m)>>32 */
+	MOVW	LO, R9		/* lo = b[i] * m */
+	ADDU	R10, R9		/* lo += oldhi */
+	SGTU	R10, R9, R2
+	ADDU	R2, R8		/* hi += carry */
+	ADDU	R9, R11		/* tmp += lo */
+	SGTU	R9, R11, R2
+	ADDU	R2, R8		/* hi += carry */
+	MOVW	R11, 0(R6)	/* p[i] = tmp */
+	ADDU	$4, R6
+	MOVW	R8, R10		/* oldhi = hi */
+	SUBU	$1, R4
+	BNE	R4, _muladdloop
+
+_muladd1:
+	MOVW	0(R6), R11	/* tmp = p[i] */
+	ADDU	R10, R11	/* tmp += oldhi */
+	MOVW	R11, 0(R6)	/* p[i] = tmp */
+
+	RET
diff --git a/sys/src/libmp/mips/mpvecdigmulsub.s b/sys/src/libmp/mips/mpvecdigmulsub.s
new file mode 100755
index 000000000..c1a1739c3
--- /dev/null
+++ b/sys/src/libmp/mips/mpvecdigmulsub.s
@@ -0,0 +1,61 @@
+/*
+ *	mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p)
+ *
+ *	p -= b*m
+ *
+ *	each step looks like:
+ *		hi,lo = m*b[i]
+ *		lo += oldhi + carry
+ *		hi += carry
+ *		p[i] += lo
+ *		oldhi = hi
+ *
+ *	the registers are:
+ *		b = R1
+ *		n = R4
+ *		m = R5
+ *		p = R6
+ *		i = R7
+ *		hi = R8		- constrained by hardware
+ *		lo = R9		- constrained by hardware
+ *		oldhi = R10
+ *		tmp = R11
+ *		
+ */
+TEXT	mpvecdigmulsub(SB),$0
+
+	MOVW	n+4(FP),R4
+	MOVW	m+8(FP),R5
+	MOVW	p+12(FP),R6
+
+	MOVW	R0, R10		/* oldhi = 0 */
+_mulsubloop:
+	MOVW	0(R1), R9	/* lo = b[i] */
+	ADDU	$4, R1
+	MOVW	0(R6), R11	/* tmp = p[i] */
+	MULU	R9, R5
+	MOVW	HI, R8		/* hi = (b[i] * m)>>32 */
+	MOVW	LO, R9		/* lo = b[i] * m */
+	ADDU	R10, R9		/* lo += oldhi */
+	SGTU	R10, R9, R2
+	ADDU	R2, R8		/* hi += carry */
+	SUBU	R9, R11, R3	/* tmp -= lo */
+	SGTU	R3, R11, R2
+	ADDU	R2, R8		/* hi += carry */
+	MOVW	R3, 0(R6)	/* p[i] = tmp */
+	ADDU	$4, R6
+	MOVW	R8, R10		/* oldhi = hi */
+	SUBU	$1, R4
+	BNE	R4, _mulsubloop
+
+	MOVW	0(R6), R11	/* tmp = p[i] */
+	SUBU	R10, R11, R3	/* tmp -= oldhi */
+	MOVW	R3, 0(R6)	/* p[i] = tmp */
+	SGTU	R3, R11, R1
+	BNE	R1, _mulsub2
+	MOVW	$1, R1		/* return +1 for positive result */
+	RET
+
+_mulsub2:
+	MOVW	$-1, R1		/* return -1 for negative result */
+	RET
diff --git a/sys/src/libmp/mips/mpvecsub.s b/sys/src/libmp/mips/mpvecsub.s
new file mode 100755
index 000000000..ca156fb7c
--- /dev/null
+++ b/sys/src/libmp/mips/mpvecsub.s
@@ -0,0 +1,66 @@
+#define	BDNZ	BC	16,0,
+#define	BDNE	BC	0,2,
+
+/*
+ *	mpvecadd(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *sum)
+ *
+ *		sum[0:alen] = a[0:alen-1] - b[0:blen-1]
+ *
+ *	prereq: alen >= blen, sum has room for alen+1 digits
+ *
+ *		R1 == a	(first arg passed in R1)
+ *		R3 == carry
+ *		R4 == alen
+ *		R5 == b
+ *		R6 == blen
+ *		R7 == sum
+ *		R2 == temporary
+ *		R8 == temporary
+ *		R9 == temporary
+ */
+TEXT	mpvecsub(SB),$-4
+
+	MOVW	alen+4(FP), R4
+	MOVW	b+8(FP), R5
+	MOVW	blen+12(FP), R6
+	MOVW	sum+16(FP), R7
+	SUBU	R6, R4		/* calculate counter for second loop (alen > blen) */
+	MOVW	R0, R3
+
+	/* if blen == 0, don't need to subtract it */
+	BEQ	R6,_sub1
+
+	/* sum[0:blen-1],carry = a[0:blen-1] - b[0:blen-1] */
+_subloop1:
+	MOVW	0(R1), R8
+	ADDU	$4, R1
+	MOVW	0(R5), R9
+	ADDU	$4, R5
+	SUBU	R3, R8, R2
+	SGTU	R2, R8, R3
+	SUBU	R9, R2, R8
+	SGTU	R8, R2, R9
+	ADDU	R9, R3
+	MOVW	R8, 0(R7)
+	ADDU	$4, R7
+	SUBU	$1, R6
+	BNE	R6, _subloop1
+
+_sub1:
+	/* if alen == blen, we're done */
+	BEQ	R4, _subend
+
+	/* sum[blen:alen-1],carry = a[blen:alen-1] + 0 + carry */
+_subloop2:
+	MOVW	0(R1), R8
+	ADDU	$4, R1
+	SUBU	R3, R8, R2
+	SGTU	R2, R8, R3
+	MOVW	R2, 0(R7)
+	ADDU	$4, R7
+	SUBU	$1, R4
+	BNE	R4, _subloop2
+
+	/* sum[alen] = carry */
+_subend:
+	RET
diff --git a/sys/src/libmp/mkfile b/sys/src/libmp/mkfile
new file mode 100755
index 000000000..886e3729f
--- /dev/null
+++ b/sys/src/libmp/mkfile
@@ -0,0 +1,54 @@
+</$objtype/mkfile
+
+DIRS=port $CPUS
+
+default:V:	all
+
+install clean all:V:
+	for(i in port $objtype)@{
+		echo $i
+		cd $i
+		mk $MKFLAGS $target
+	}
+
+nuke:V: clean
+	rm -f /$objtype/lib/libmp.a
+
+update:V:
+	for(i in port $CPUS)@{
+		echo $i
+		cd $i
+		mk $MKFLAGS $target
+	}
+	update /386/lib/libmp.a
+
+installall:V:
+	for(objtype in $CPUS) mk $MKFLAGS install
+
+everything:V:
+	rm -f */*.[012456789kvx]
+	for(objtype in 386)@{
+		echo $objtype
+		mk $MKFLAGS install
+	}
+	rm -f */*.[012456789kvx]
+
+test.$O: test.c /$objtype/include/u.h /sys/include/mp.h port/dat.h
+	$CC -Iport test.c
+
+$O.test: test.$O /$objtype/lib/libmp.a
+	$LD -o $O.test test.$O
+
+bigtest.$O: bigtest.c /$objtype/include/u.h /sys/include/mp.h port/dat.h
+	$CC -Iport bigtest.c
+
+$O.bigtest: bigtest.$O /$objtype/lib/libmp.a
+	$LD -o $O.bigtest bigtest.$O
+
+allout:
+	objtype=386; mk; mk 8.test 8.bigtest
+	objtype=power; mk; mk q.test q.bigtest
+	objtype=mips; mk; mk v.test v.bigtest
+
+cleanout:
+	rm -f [qv8].* *.[qv8]
diff --git a/sys/src/libmp/port/betomp.c b/sys/src/libmp/port/betomp.c
new file mode 100755
index 000000000..9197f3a14
--- /dev/null
+++ b/sys/src/libmp/port/betomp.c
@@ -0,0 +1,42 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+// convert a big-endian byte array (most significant byte first) to an mpint
+mpint*
+betomp(uchar *p, uint n, mpint *b)
+{
+	int m, s;
+	mpdigit x;
+
+	if(b == nil){
+		b = mpnew(0);
+		setmalloctag(b, getcallerpc(&p));
+	}
+
+	// dump leading zeros
+	while(*p == 0 && n > 1){
+		p++;
+		n--;
+	}
+
+	// get the space
+	mpbits(b, n*8);
+	b->top = DIGITS(n*8);
+	m = b->top-1;
+
+	// first digit might not be Dbytes long
+	s = ((n-1)*8)%Dbits;
+	x = 0;
+	for(; n > 0; n--){
+		x |= ((mpdigit)(*p++)) << s;
+		s -= 8;
+		if(s < 0){
+			b->p[m--] = x;
+			s = Dbits-8;
+			x = 0;
+		}
+	}
+
+	return b;
+}
diff --git a/sys/src/libmp/port/crt.c b/sys/src/libmp/port/crt.c
new file mode 100755
index 000000000..ddf26ed57
--- /dev/null
+++ b/sys/src/libmp/port/crt.c
@@ -0,0 +1,122 @@
+#include "os.h"
+#include <mp.h>
+#include <libsec.h>
+
+// chinese remainder theorem
+//
+// handbook of applied cryptography, menezes et al, 1997, pp 610 - 613
+
+struct CRTpre
+{
+	int	n;		// number of moduli
+	mpint	**m;		// pointer to moduli
+	mpint	**c;		// precomputed coefficients
+	mpint	**p;		// precomputed products
+	mpint	*a[1];		// local storage
+};
+
+// setup crt info, returns a newly created structure
+CRTpre*
+crtpre(int n, mpint **m)
+{
+	CRTpre *crt;
+	int i, j;
+	mpint *u;
+
+	crt = malloc(sizeof(CRTpre)+sizeof(mpint)*3*n);
+	if(crt == nil)
+		sysfatal("crtpre: %r");
+	crt->m = crt->a;
+	crt->c = crt->a+n;
+	crt->p = crt->c+n;
+	crt->n = n;
+
+	// make a copy of the moduli
+	for(i = 0; i < n; i++)
+		crt->m[i] = mpcopy(m[i]);
+
+	// precompute the products
+	u = mpcopy(mpone);
+	for(i = 0; i < n; i++){
+		mpmul(u, m[i], u);
+		crt->p[i] = mpcopy(u);
+	}
+
+	// precompute the coefficients
+	for(i = 1; i < n; i++){
+		crt->c[i] = mpcopy(mpone);
+		for(j = 0; j < i; j++){
+			mpinvert(m[j], m[i], u);
+			mpmul(u, crt->c[i], u);
+			mpmod(u, m[i], crt->c[i]);
+		}
+	}
+
+	mpfree(u);
+
+	return crt;		
+}
+
+void
+crtprefree(CRTpre *crt)
+{
+	int i;
+
+	for(i = 0; i < crt->n; i++){
+		if(i != 0)
+			mpfree(crt->c[i]);
+		mpfree(crt->p[i]);
+		mpfree(crt->m[i]);
+	}
+	free(crt);
+}
+
+// convert to residues, returns a newly created structure
+CRTres*
+crtin(CRTpre *crt, mpint *x)
+{
+	int i;
+	CRTres *res;
+
+	res = malloc(sizeof(CRTres)+sizeof(mpint)*crt->n);
+	if(res == nil)
+		sysfatal("crtin: %r");
+	res->n = crt->n;
+	for(i = 0; i < res->n; i++){
+		res->r[i] = mpnew(0);
+		mpmod(x, crt->m[i], res->r[i]);
+	}
+	return res;
+}
+
+// garners algorithm for converting residue form to linear
+void
+crtout(CRTpre *crt, CRTres *res, mpint *x)
+{
+	mpint *u;
+	int i;
+
+	u = mpnew(0);
+	mpassign(res->r[0], x);
+
+	for(i = 1; i < crt->n; i++){
+		mpsub(res->r[i], x, u);
+		mpmul(u, crt->c[i], u);
+		mpmod(u, crt->m[i], u);
+		mpmul(u, crt->p[i-1], u);
+		mpadd(x, u, x);
+	}
+
+	mpfree(u);
+}
+
+// free the residue
+void
+crtresfree(CRTres *res)
+{
+	int i;
+
+	for(i = 0; i < res->n; i++)
+		mpfree(res->r[i]);
+	free(res);
+}
diff --git a/sys/src/libmp/port/crttest.c b/sys/src/libmp/port/crttest.c
new file mode 100755
index 000000000..ef3d18e1f
--- /dev/null
+++ b/sys/src/libmp/port/crttest.c
@@ -0,0 +1,53 @@
+#include "os.h"
+#include <mp.h>
+#include <libsec.h>
+
+void
+testcrt(mpint **p)
+{
+	CRTpre *crt;
+	CRTres *res;
+	mpint *m, *x, *y;
+
+	fmtinstall('B', mpfmt);
+
+	// get a modulus and a test number
+	m = mpnew(1024+160);
+	mpmul(p[0], p[1], m);
+	x = mpnew(1024+160);
+	mpadd(m, mpone, x);
+
+	// do the precomputation for crt conversion
+	crt = crtpre(2, p);
+
+	// convert x to residues
+	res = crtin(crt, x);
+
+	// convert back
+	y = mpnew(1024+160);
+	crtout(crt, res, y);
+	print("x %B\ny %B\n", x, y);
+	mpfree(m);
+	mpfree(x);
+	mpfree(y);
+}
+
+void
+main(void)
+{
+	int i;
+	mpint *p[2];
+	long start;
+
+	start = time(0);
+	for(i = 0; i < 10; i++){
+		p[0] = mpnew(1024);
+		p[1] = mpnew(1024);
+		DSAprimes(p[0], p[1], nil);
+		testcrt(p);
+		mpfree(p[0]);
+		mpfree(p[1]);
+	}
+	print("%ld secs with more\n", time(0)-start);
+	exits(0);
+}
diff --git a/sys/src/libmp/port/dat.h b/sys/src/libmp/port/dat.h
new file mode 100755
index 000000000..7c834ac81
--- /dev/null
+++ b/sys/src/libmp/port/dat.h
@@ -0,0 +1,12 @@
+#define	mpdighi  (mpdigit)(1<<(Dbits-1))
+#define DIGITS(x) ((Dbits - 1 + (x))/Dbits)
+
+// for converting between int's and mpint's
+#define MAXUINT ((uint)-1)
+#define MAXINT (MAXUINT>>1)
+#define MININT (MAXINT+1)
+
+// for converting between vlongs's and mpint's
+#define MAXUVLONG (~0ULL)
+#define MAXVLONG (MAXUVLONG>>1)
+#define MINVLONG (MAXVLONG+1ULL)
diff --git a/sys/src/libmp/port/letomp.c b/sys/src/libmp/port/letomp.c
new file mode 100755
index 000000000..e23fed21e
--- /dev/null
+++ b/sys/src/libmp/port/letomp.c
@@ -0,0 +1,28 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+// convert a little endian byte array (least significant byte first) to an mpint
+mpint*
+letomp(uchar *s, uint n, mpint *b)
+{
+	int i=0, m = 0;
+	mpdigit x=0;
+
+	if(b == nil)
+		b = mpnew(0);
+	mpbits(b, 8*n);
+	for(; n > 0; n--){
+		x |= ((mpdigit)(*s++)) << i;
+		i += 8;
+		if(i == Dbits){
+			b->p[m++] = x;
+			i = 0;
+			x = 0;
+		}
+	}
+	if(i > 0)
+		b->p[m++] = x;
+	b->top = m;
+	return b;
+}
diff --git a/sys/src/libmp/port/mkfile b/sys/src/libmp/port/mkfile
new file mode 100755
index 000000000..05ee9dd33
--- /dev/null
+++ b/sys/src/libmp/port/mkfile
@@ -0,0 +1,52 @@
+</$objtype/mkfile
+
+LIB=/$objtype/lib/libmp.a
+FILES=\
+	mpaux\
+	mpfmt\
+	strtomp\
+	mptobe\
+	mptole\
+	betomp\
+	letomp\
+	mpadd\
+	mpsub\
+	mpcmp\
+	mpfactorial\
+	mpmul\
+	mpleft\
+	mpright\
+	mpvecadd\
+	mpvecsub\
+	mpvecdigmuladd\
+	mpveccmp\
+	mpdigdiv\
+	mpdiv\
+	mpexp\
+	mpmod\
+	mpextendedgcd\
+	mpinvert\
+	mprand\
+	crt\
+	mptoi\
+	mptoui\
+	mptov\
+	mptouv\
+
+ALLOFILES=${FILES:%=%.$O}
+# cull things in the per-machine directories from this list
+OFILES=	`{rc ./reduce $O $objtype $ALLOFILES}
+
+HFILES=\
+	/$objtype/include/u.h\
+	/sys/include/mp.h\
+	dat.h\
+
+CFILES=${FILES:%=%.c}
+
+
+UPDATE=mkfile\
+	$HFILES\
+	$CFILES\
+
+</sys/src/cmd/mksyslib
diff --git a/sys/src/libmp/port/mpadd.c b/sys/src/libmp/port/mpadd.c
new file mode 100755
index 000000000..6022a64ef
--- /dev/null
+++ b/sys/src/libmp/port/mpadd.c
@@ -0,0 +1,54 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+// sum = abs(b1) + abs(b2), i.e., add the magnitudes
+void
+mpmagadd(mpint *b1, mpint *b2, mpint *sum)
+{
+	int m, n;
+	mpint *t;
+
+	// get the sizes right
+	if(b2->top > b1->top){
+		t = b1;
+		b1 = b2;
+		b2 = t;
+	}
+	n = b1->top;
+	m = b2->top;
+	if(n == 0){
+		mpassign(mpzero, sum);
+		return;
+	}
+	if(m == 0){
+		mpassign(b1, sum);
+		return;
+	}
+	mpbits(sum, (n+1)*Dbits);
+	sum->top = n+1;
+
+	mpvecadd(b1->p, n, b2->p, m, sum->p);
+	sum->sign = 1;
+
+	mpnorm(sum);
+}
+
+// sum = b1 + b2
+void
+mpadd(mpint *b1, mpint *b2, mpint *sum)
+{
+	int sign;
+
+	if(b1->sign != b2->sign){
+		if(b1->sign < 0)
+			mpmagsub(b2, b1, sum);
+		else
+			mpmagsub(b1, b2, sum);
+	} else {
+		sign = b1->sign;
+		mpmagadd(b1, b2, sum);
+		if(sum->top != 0)
+			sum->sign = sign;
+	}
+}
diff --git a/sys/src/libmp/port/mpaux.c b/sys/src/libmp/port/mpaux.c
new file mode 100755
index 000000000..66f1524f0
--- /dev/null
+++ b/sys/src/libmp/port/mpaux.c
@@ -0,0 +1,190 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+static mpdigit _mptwodata[1] = { 2 };
+static mpint _mptwo =
+{
+	1,
+	1,
+	1,
+	_mptwodata,
+	MPstatic
+};
+mpint *mptwo = &_mptwo;
+
+static mpdigit _mponedata[1] = { 1 };
+static mpint _mpone =
+{
+	1,
+	1,
+	1,
+	_mponedata,
+	MPstatic
+};
+mpint *mpone = &_mpone;
+
+static mpdigit _mpzerodata[1] = { 0 };
+static mpint _mpzero =
+{
+	1,
+	1,
+	0,
+	_mpzerodata,
+	MPstatic
+};
+mpint *mpzero = &_mpzero;
+
+static int mpmindigits = 33;
+
+// set minimum digit allocation
+void
+mpsetminbits(int n)
+{
+	if(n < 0)
+		sysfatal("mpsetminbits: n < 0");
+	if(n == 0)
+		n = 1;
+	mpmindigits = DIGITS(n);
+}
+
+// allocate an n bit 0'd number 
+mpint*
+mpnew(int n)
+{
+	mpint *b;
+
+	if(n < 0)
+		sysfatal("mpsetminbits: n < 0");
+
+	b = mallocz(sizeof(mpint), 1);
+	setmalloctag(b, getcallerpc(&n));
+	if(b == nil)
+		sysfatal("mpnew: %r");
+	n = DIGITS(n);
+	if(n < mpmindigits)
+		n = mpmindigits;
+	b->p = (mpdigit*)mallocz(n*Dbytes, 1);
+	if(b->p == nil)
+		sysfatal("mpnew: %r");
+	b->size = n;
+	b->sign = 1;
+
+	return b;
+}
+
+// guarantee at least n significant bits
+void
+mpbits(mpint *b, int m)
+{
+	int n;
+
+	n = DIGITS(m);
+	if(b->size >= n){
+		if(b->top >= n)
+			return;
+		memset(&b->p[b->top], 0, Dbytes*(n - b->top));
+		b->top = n;
+		return;
+	}
+	b->p = (mpdigit*)realloc(b->p, n*Dbytes);
+	if(b->p == nil)
+		sysfatal("mpbits: %r");
+	memset(&b->p[b->top], 0, Dbytes*(n - b->top));
+	b->size = n;
+	b->top = n;
+}
+
+void
+mpfree(mpint *b)
+{
+	if(b == nil)
+		return;
+	if(b->flags & MPstatic)
+		sysfatal("freeing mp constant");
+	memset(b->p, 0, b->size*Dbytes);	// information hiding
+	free(b->p);
+	free(b);
+}
+
+void
+mpnorm(mpint *b)
+{
+	int i;
+
+	for(i = b->top-1; i >= 0; i--)
+		if(b->p[i] != 0)
+			break;
+	b->top = i+1;
+	if(b->top == 0)
+		b->sign = 1;
+}
+
+mpint*
+mpcopy(mpint *old)
+{
+	mpint *new;
+
+	new = mpnew(Dbits*old->size);
+	new->top = old->top;
+	new->sign = old->sign;
+	memmove(new->p, old->p, Dbytes*old->top);
+	return new;
+}
+
+void
+mpassign(mpint *old, mpint *new)
+{
+	mpbits(new, Dbits*old->top);
+	new->sign = old->sign;
+	new->top = old->top;
+	memmove(new->p, old->p, Dbytes*old->top);
+}
+
+// number of significant bits in mantissa
+int
+mpsignif(mpint *n)
+{
+	int i, j;
+	mpdigit d;
+
+	if(n->top == 0)
+		return 0;
+	for(i = n->top-1; i >= 0; i--){
+		d = n->p[i];
+		for(j = Dbits-1; j >= 0; j--){
+			if(d & (((mpdigit)1)<<j))
+				return i*Dbits + j + 1;
+		}
+	}
+	return 0;
+}
+
+// k, where n = 2**k * q for odd q
+int
+mplowbits0(mpint *n)
+{
+	int k, bit, digit;
+	mpdigit d;
+
+	if(n->top==0)
+		return 0;
+	k = 0;
+	bit = 0;
+	digit = 0;
+	d = n->p[0];
+	for(;;){
+		if(d & (1<<bit))
+			break;
+		k++;
+		bit++;
+		if(bit==Dbits){
+			if(++digit >= n->top)
+				return 0;
+			d = n->p[digit];
+			bit = 0;
+		}
+	}
+	return k;
+}
+
diff --git a/sys/src/libmp/port/mpcmp.c b/sys/src/libmp/port/mpcmp.c
new file mode 100755
index 000000000..a2e3cf724
--- /dev/null
+++ b/sys/src/libmp/port/mpcmp.c
@@ -0,0 +1,28 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+// return neg, 0, pos as abs(b1)-abs(b2) is neg, 0, pos
+int
+mpmagcmp(mpint *b1, mpint *b2)
+{
+	int i;
+
+	i = b1->top - b2->top;
+	if(i)
+		return i;
+
+	return mpveccmp(b1->p, b1->top, b2->p, b2->top);
+}
+
+// return neg, 0, pos as b1-b2 is neg, 0, pos
+int
+mpcmp(mpint *b1, mpint *b2)
+{
+	if(b1->sign != b2->sign)
+		return b1->sign - b2->sign;
+	if(b1->sign < 0)
+		return mpmagcmp(b2, b1);
+	else
+		return mpmagcmp(b1, b2);
+}
diff --git a/sys/src/libmp/port/mpdigdiv.c b/sys/src/libmp/port/mpdigdiv.c
new file mode 100755
index 000000000..4a73bb3a4
--- /dev/null
+++ b/sys/src/libmp/port/mpdigdiv.c
@@ -0,0 +1,43 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+//
+//	divide two digits by one and return quotient
+//
+void
+mpdigdiv(mpdigit *dividend, mpdigit divisor, mpdigit *quotient)
+{
+	mpdigit hi, lo, q, x, y;
+	int i;
+
+	hi = dividend[1];
+	lo = dividend[0];
+
+	// return highest digit value if the result >= 2**32
+	if(hi >= divisor || divisor == 0){
+		divisor = 0;
+		*quotient = ~divisor;
+		return;
+	}
+
+	// at this point we know that hi < divisor
+	// just shift and subtract till we're done
+	q = 0;
+	x = divisor;
+	for(i = Dbits-1; hi > 0 && i >= 0; i--){
+		x >>= 1;
+		if(x > hi)
+			continue;
+		y = divisor<<i;
+		if(x == hi && y > lo)
+			continue;
+		if(y > lo)
+			hi--;
+		lo -= y;
+		hi -= x;
+		q |= 1<<i;
+	}
+	q += lo/divisor;
+	*quotient = q;
+}
diff --git a/sys/src/libmp/port/mpdiv.c b/sys/src/libmp/port/mpdiv.c
new file mode 100755
index 000000000..92aee03f4
--- /dev/null
+++ b/sys/src/libmp/port/mpdiv.c
@@ -0,0 +1,112 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+// division ala knuth, seminumerical algorithms, pp 237-238
+// the numbers are stored backwards to what knuth expects so j
+// counts down rather than up.
+
+void
+mpdiv(mpint *dividend, mpint *divisor, mpint *quotient, mpint *remainder)
+{
+	int j, s, vn, sign;
+	mpdigit qd, *up, *vp, *qp;
+	mpint *u, *v, *t;
+
+	// divide bv zero
+	if(divisor->top == 0)
+		abort();
+
+	// quick check
+	if(mpmagcmp(dividend, divisor) < 0){
+		if(remainder != nil)
+			mpassign(dividend, remainder);
+		if(quotient != nil)
+			mpassign(mpzero, quotient);
+		return;
+	}
+
+	// D1: shift until divisor, v, has hi bit set (needed to make trial
+	//     divisor accurate)
+	qd = divisor->p[divisor->top-1];
+	for(s = 0; (qd & mpdighi) == 0; s++)
+		qd <<= 1;
+	u = mpnew((dividend->top+2)*Dbits + s);
+	if(s == 0 && divisor != quotient && divisor != remainder) {
+		mpassign(dividend, u);
+		v = divisor;
+	} else {
+		mpleft(dividend, s, u);
+		v = mpnew(divisor->top*Dbits);
+		mpleft(divisor, s, v);
+	}
+	up = u->p+u->top-1;
+	vp = v->p+v->top-1;
+	vn = v->top;
+
+	// D1a: make sure high digit of dividend is less than high digit of divisor
+	if(*up >= *vp){
+		*++up = 0;
+		u->top++;
+	}
+
+	// storage for multiplies
+	t = mpnew(4*Dbits);
+
+	qp = nil;
+	if(quotient != nil){
+		mpbits(quotient, (u->top - v->top)*Dbits);
+		quotient->top = u->top - v->top;
+		qp = quotient->p+quotient->top-1;
+	}
+
+	// D2, D7: loop on length of dividend
+	for(j = u->top; j > vn; j--){
+
+		// D3: calculate trial divisor
+		mpdigdiv(up-1, *vp, &qd);
+
+		// D3a: rule out trial divisors 2 greater than real divisor
+		if(vn > 1) for(;;){
+			memset(t->p, 0, 3*Dbytes);	// mpvecdigmuladd adds to what's there
+			mpvecdigmuladd(vp-1, 2, qd, t->p);
+			if(mpveccmp(t->p, 3, up-2, 3) > 0)
+				qd--;
+			else
+				break;
+		}
+
+		// D4: u -= v*qd << j*Dbits
+		sign = mpvecdigmulsub(v->p, vn, qd, up-vn);
+		if(sign < 0){
+
+			// D6: trial divisor was too high, add back borrowed
+			//     value and decrease divisor
+			mpvecadd(up-vn, vn+1, v->p, vn, up-vn);
+			qd--;
+		}
+
+		// D5: save quotient digit
+		if(qp != nil)
+			*qp-- = qd;
+
+		// push top of u down one
+		u->top--;
+		*up-- = 0;
+	}
+	if(qp != nil){
+		mpnorm(quotient);
+		if(dividend->sign != divisor->sign)
+			quotient->sign = -1;
+	}
+
+	if(remainder != nil){
+		mpright(u, s, remainder);	// u is the remainder shifted
+		remainder->sign = dividend->sign;
+	}
+
+	mpfree(t);
+	mpfree(u);
+	if(v != divisor)
+		mpfree(v);
+}
diff --git a/sys/src/libmp/port/mpeuclid.c b/sys/src/libmp/port/mpeuclid.c
new file mode 100755
index 000000000..80b5983bf
--- /dev/null
+++ b/sys/src/libmp/port/mpeuclid.c
@@ -0,0 +1,85 @@
+#include "os.h"
+#include <mp.h>
+
+// extended euclid
+//
+// For a and b it solves, d = gcd(a,b) and finds x and y s.t.
+// ax + by = d
+//
+// Handbook of Applied Cryptography, Menezes et al, 1997, pg 67
+
+void
+mpeuclid(mpint *a, mpint *b, mpint *d, mpint *x, mpint *y)
+{
+	mpint *tmp, *x0, *x1, *x2, *y0, *y1, *y2, *q, *r;
+
+	if(a->sign<0 || b->sign<0)
+		sysfatal("mpeuclid: negative arg");
+
+	if(mpcmp(a, b) < 0){
+		tmp = a;
+		a = b;
+		b = tmp;
+		tmp = x;
+		x = y;
+		y = tmp;
+	}
+
+	if(b->top == 0){
+		mpassign(a, d);
+		mpassign(mpone, x);
+		mpassign(mpzero, y);
+		return;
+	}
+
+	a = mpcopy(a);
+	b = mpcopy(b);
+	x0 = mpnew(0);
+	x1 = mpcopy(mpzero);
+	x2 = mpcopy(mpone);
+	y0 = mpnew(0);
+	y1 = mpcopy(mpone);
+	y2 = mpcopy(mpzero);
+	q = mpnew(0);
+	r = mpnew(0);
+
+	while(b->top != 0 && b->sign > 0){
+		// q = a/b
+		// r = a mod b
+		mpdiv(a, b, q, r);
+		// x0 = x2 - qx1
+		mpmul(q, x1, x0);
+		mpsub(x2, x0, x0);
+		// y0 = y2 - qy1
+		mpmul(q, y1, y0);
+		mpsub(y2, y0, y0);
+		// rotate values
+		tmp = a;
+		a = b;
+		b = r;
+		r = tmp;
+		tmp = x2;
+		x2 = x1;
+		x1 = x0;
+		x0 = tmp;
+		tmp = y2;
+		y2 = y1;
+		y1 = y0;
+		y0 = tmp;
+	}
+
+	mpassign(a, d);
+	mpassign(x2, x);
+	mpassign(y2, y);
+
+	mpfree(x0);
+	mpfree(x1);
+	mpfree(x2);
+	mpfree(y0);
+	mpfree(y1);
+	mpfree(y2);
+	mpfree(q);
+	mpfree(r);
+	mpfree(a);
+	mpfree(b);
+}
diff --git a/sys/src/libmp/port/mpexp.c b/sys/src/libmp/port/mpexp.c
new file mode 100755
index 000000000..9ec067cb9
--- /dev/null
+++ b/sys/src/libmp/port/mpexp.c
@@ -0,0 +1,94 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+// res = b**e
+//
+// knuth, vol 2, pp 398-400
+
+enum {
+	Freeb=	0x1,
+	Freee=	0x2,
+	Freem=	0x4,
+};
+
+//int expdebug;
+
+void
+mpexp(mpint *b, mpint *e, mpint *m, mpint *res)
+{
+	mpint *t[2];
+	int tofree;
+	mpdigit d, bit;
+	int i, j;
+
+	i = mpcmp(e,mpzero);
+	if(i==0){
+		mpassign(mpone, res);
+		return;
+	}
+	if(i<0)
+		sysfatal("mpexp: negative exponent");
+
+	t[0] = mpcopy(b);
+	t[1] = res;
+
+	tofree = 0;
+	if(res == b){
+		b = mpcopy(b);
+		tofree |= Freeb;
+	}
+	if(res == e){
+		e = mpcopy(e);
+		tofree |= Freee;
+	}
+	if(res == m){
+		m = mpcopy(m);
+		tofree |= Freem;
+	}
+
+	// skip first bit
+	i = e->top-1;
+	d = e->p[i];
+	for(bit = mpdighi; (bit & d) == 0; bit >>= 1)
+		;
+	bit >>= 1;
+
+	j = 0;
+	for(;;){
+		for(; bit != 0; bit >>= 1){
+			mpmul(t[j], t[j], t[j^1]);
+			if(bit & d)
+				mpmul(t[j^1], b, t[j]);
+			else
+				j ^= 1;
+			if(m != nil && t[j]->top > m->top){
+				mpmod(t[j], m, t[j^1]);
+				j ^= 1;
+			}
+		}
+		if(--i < 0)
+			break;
+		bit = mpdighi;
+		d = e->p[i];
+	}
+	if(m != nil){
+		mpmod(t[j], m, t[j^1]);
+		j ^= 1;
+	}
+	if(t[j] == res){
+		mpfree(t[j^1]);
+	} else {
+		mpassign(t[j], res);
+		mpfree(t[j]);
+	}
+
+	if(tofree){
+		if(tofree & Freeb)
+			mpfree(b);
+		if(tofree & Freee)
+			mpfree(e);
+		if(tofree & Freem)
+			mpfree(m);
+	}
+}
diff --git a/sys/src/libmp/port/mpextendedgcd.c b/sys/src/libmp/port/mpextendedgcd.c
new file mode 100755
index 000000000..413a05c2a
--- /dev/null
+++ b/sys/src/libmp/port/mpextendedgcd.c
@@ -0,0 +1,106 @@
+#include "os.h"
+#include <mp.h>
+
+#define iseven(a)	(((a)->p[0] & 1) == 0)
+
+// extended binary gcd
+//
+// For a anv b it solves, v = gcd(a,b) and finds x and y s.t.
+// ax + by = v
+//
+// Handbook of Applied Cryptography, Menezes et al, 1997, pg 608.  
+void
+mpextendedgcd(mpint *a, mpint *b, mpint *v, mpint *x, mpint *y)
+{
+	mpint *u, *A, *B, *C, *D;
+	int g;
+
+	if(a->sign < 0 || b->sign < 0){
+		mpassign(mpzero, v);
+		mpassign(mpzero, y);
+		mpassign(mpzero, x);
+		return;
+	}
+
+	if(a->top == 0){
+		mpassign(b, v);
+		mpassign(mpone, y);
+		mpassign(mpzero, x);
+		return;
+	}
+	if(b->top == 0){
+		mpassign(a, v);
+		mpassign(mpone, x);
+		mpassign(mpzero, y);
+		return;
+	}
+
+	g = 0;
+	a = mpcopy(a);
+	b = mpcopy(b);
+
+	while(iseven(a) && iseven(b)){
+		mpright(a, 1, a);
+		mpright(b, 1, b);
+		g++;
+	}
+
+	u = mpcopy(a);
+	mpassign(b, v);
+	A = mpcopy(mpone);
+	B = mpcopy(mpzero);
+	C = mpcopy(mpzero);
+	D = mpcopy(mpone);
+
+	for(;;) {
+//		print("%B %B %B %B %B %B\n", u, v, A, B, C, D);
+		while(iseven(u)){
+			mpright(u, 1, u);
+			if(!iseven(A) || !iseven(B)) {
+				mpadd(A, b, A);
+				mpsub(B, a, B);
+			}
+			mpright(A, 1, A);
+			mpright(B, 1, B);
+		}
+	
+//		print("%B %B %B %B %B %B\n", u, v, A, B, C, D);
+		while(iseven(v)){
+			mpright(v, 1, v);
+			if(!iseven(C) || !iseven(D)) {
+				mpadd(C, b, C);
+				mpsub(D, a, D);
+			}
+			mpright(C, 1, C);
+			mpright(D, 1, D);
+		}
+	
+//		print("%B %B %B %B %B %B\n", u, v, A, B, C, D);
+		if(mpcmp(u, v) >= 0){
+			mpsub(u, v, u);
+			mpsub(A, C, A);
+			mpsub(B, D, B);
+		} else {
+			mpsub(v, u, v);
+			mpsub(C, A, C);
+			mpsub(D, B, D);
+		}
+
+		if(u->top == 0)
+			break;
+
+	}
+	mpassign(C, x);
+	mpassign(D, y);
+	mpleft(v, g, v);
+
+	mpfree(A);
+	mpfree(B);
+	mpfree(C);
+	mpfree(D);
+	mpfree(u);
+	mpfree(a);
+	mpfree(b);
+
+	return;
+}
diff --git a/sys/src/libmp/port/mpfactorial.c b/sys/src/libmp/port/mpfactorial.c
new file mode 100755
index 000000000..01079c44b
--- /dev/null
+++ b/sys/src/libmp/port/mpfactorial.c
@@ -0,0 +1,75 @@
+#include "os.h"
+#include <mp.h>
+#include <libsec.h>
+#include "dat.h"
+
+mpint*
+mpfactorial(ulong n)
+{
+	int i;
+	ulong k;
+	unsigned cnt;
+	int max, mmax;
+	mpdigit p, pp[2];
+	mpint *r, *s, *stk[31];
+
+	cnt = 0;
+	max = mmax = -1;
+	p = 1;
+	r = mpnew(0);
+	for(k=2; k<=n; k++){
+		pp[0] = 0;
+		pp[1] = 0;
+		mpvecdigmuladd(&p, 1, (mpdigit)k, pp);
+		if(pp[1] == 0)	/* !overflow */
+			p = pp[0];
+		else{
+			cnt++;
+			if((cnt & 1) == 0){
+				s = stk[max];
+				mpbits(r, Dbits*(s->top+1+1));
+				memset(r->p, 0, Dbytes*(s->top+1+1));
+				mpvecmul(s->p, s->top, &p, 1, r->p);
+				r->sign = 1;
+				r->top = s->top+1+1;		/* XXX: norm */
+				mpassign(r, s);
+				for(i=4; (cnt & (i-1)) == 0; i=i<<1){
+					mpmul(stk[max], stk[max-1], r);
+					mpassign(r, stk[max-1]);
+					max--;
+				}
+			}else{
+				max++;
+				if(max > mmax){
+					mmax++;
+					if(max > nelem(stk))
+						abort();
+					stk[max] = mpnew(Dbits);
+				}
+				stk[max]->top = 1;
+				stk[max]->p[0] = p;
+			}
+			p = (mpdigit)k;
+		}
+	}
+	if(max < 0){
+		mpbits(r, Dbits);
+		r->top = 1;
+		r->sign = 1;
+		r->p[0] = p;
+	}else{
+		s = stk[max--];
+		mpbits(r, Dbits*(s->top+1+1));
+		memset(r->p, 0, Dbytes*(s->top+1+1));
+		mpvecmul(s->p, s->top, &p, 1, r->p);
+		r->sign = 1;
+		r->top = s->top+1+1;		/* XXX: norm */
+	}
+
+	while(max >= 0)
+		mpmul(r, stk[max--], r);
+	for(max=mmax; max>=0; max--)
+		mpfree(stk[max]);
+	mpnorm(r);
+	return r;
+}
diff --git a/sys/src/libmp/port/mpfmt.c b/sys/src/libmp/port/mpfmt.c
new file mode 100755
index 000000000..f7c42a7bc
--- /dev/null
+++ b/sys/src/libmp/port/mpfmt.c
@@ -0,0 +1,193 @@
+#include "os.h"
+#include <mp.h>
+#include <libsec.h>
+#include "dat.h"
+
+static int
+to64(mpint *b, char *buf, int len)
+{
+	uchar *p;
+	int n, rv;
+
+	p = nil;
+	n = mptobe(b, nil, 0, &p);
+	if(n < 0)
+		return -1;
+	rv = enc64(buf, len, p, n);
+	free(p);
+	return rv;
+}
+
+static int
+to32(mpint *b, char *buf, int len)
+{
+	uchar *p;
+	int n, rv;
+
+	// leave room for a multiple of 5 buffer size
+	n = b->top*Dbytes + 5;
+	p = malloc(n);
+	if(p == nil)
+		return -1;
+	n = mptobe(b, p, n, nil);
+	if(n < 0)
+		return -1;
+
+	// round up buffer size, enc32 only accepts a multiple of 5
+	if(n%5)
+		n += 5 - (n%5);
+	rv = enc32(buf, len, p, n);
+	free(p);
+	return rv;
+}
+
+static char set16[] = "0123456789ABCDEF";
+
+static int
+to16(mpint *b, char *buf, int len)
+{
+	mpdigit *p, x;
+	int i, j;
+	char *out, *eout;
+
+	if(len < 1)
+		return -1;
+
+	out = buf;
+	eout = buf+len;
+	for(p = &b->p[b->top-1]; p >= b->p; p--){
+		x = *p;
+		for(i = Dbits-4; i >= 0; i -= 4){
+			j = 0xf & (x>>i);
+			if(j != 0 || out != buf){
+				if(out >= eout)
+					return -1;
+				*out++ = set16[j];
+			}
+		}
+	}
+	if(out == buf)
+		*out++ = '0';
+	if(out >= eout)
+		return -1;
+	*out = 0;
+	return 0;
+}
+
+static char*
+modbillion(int rem, ulong r, char *out, char *buf)
+{
+	ulong rr;
+	int i;
+
+	for(i = 0; i < 9; i++){
+		rr = r%10;
+		r /= 10;
+		if(out <= buf)
+			return nil;
+		*--out = '0' + rr;
+		if(rem == 0 && r == 0)
+			break;
+	}
+	return out;
+}
+
+static int
+to10(mpint *b, char *buf, int len)
+{
+	mpint *d, *r, *billion;
+	char *out;
+
+	if(len < 1)
+		return -1;
+
+	d = mpcopy(b);
+	r = mpnew(0);
+	billion = uitomp(1000000000, nil);
+	out = buf+len;
+	*--out = 0;
+	do {
+		mpdiv(d, billion, d, r);
+		out = modbillion(d->top, r->p[0], out, buf);
+		if(out == nil)
+			break;
+	} while(d->top != 0);
+	mpfree(d);
+	mpfree(r);
+	mpfree(billion);
+
+	if(out == nil)
+		return -1;
+	len -= out-buf;
+	if(out != buf)
+		memmove(buf, out, len);
+	return 0;
+}
+
+int
+mpfmt(Fmt *fmt)
+{
+	mpint *b;
+	char *p;
+
+	b = va_arg(fmt->args, mpint*);
+	if(b == nil)
+		return fmtstrcpy(fmt, "*");
+	
+	p = mptoa(b, fmt->prec, nil, 0);
+	fmt->flags &= ~FmtPrec;
+
+	if(p == nil)
+		return fmtstrcpy(fmt, "*");
+	else{
+		fmtstrcpy(fmt, p);
+		free(p);
+		return 0;
+	}
+}
+
+char*
+mptoa(mpint *b, int base, char *buf, int len)
+{
+	char *out;
+	int rv, alloced;
+
+	alloced = 0;
+	if(buf == nil){
+		len = ((b->top+1)*Dbits+2)/3 + 1;
+		buf = malloc(len);
+		if(buf == nil)
+			return nil;
+		alloced = 1;
+	}
+
+	if(len < 2)
+		return nil;
+
+	out = buf;
+	if(b->sign < 0){
+		*out++ = '-';
+		len--;
+	}
+	switch(base){
+	case 64:
+		rv = to64(b, out, len);
+		break;
+	case 32:
+		rv = to32(b, out, len);
+		break;
+	default:
+	case 16:
+		rv = to16(b, out, len);
+		break;
+	case 10:
+		rv = to10(b, out, len);
+		break;
+	}
+	if(rv < 0){
+		if(alloced)
+			free(buf);
+		return nil;
+	}
+	return buf;
+}
diff --git a/sys/src/libmp/port/mpinvert.c b/sys/src/libmp/port/mpinvert.c
new file mode 100755
index 000000000..ee2630702
--- /dev/null
+++ b/sys/src/libmp/port/mpinvert.c
@@ -0,0 +1,21 @@
+#include "os.h"
+#include <mp.h>
+
+#define iseven(a)	(((a)->p[0] & 1) == 0)
+
+// use extended gcd to find the multiplicative inverse
+// res = b**-1 mod m
+void
+mpinvert(mpint *b, mpint *m, mpint *res)
+{
+	mpint *dc1, *dc2;	// don't care
+
+	dc1 = mpnew(0);
+	dc2 = mpnew(0);
+	mpextendedgcd(b, m, dc1, res, dc2);
+	if(mpcmp(dc1, mpone) != 0)
+		abort();
+	mpmod(res, m, res);
+	mpfree(dc1);
+	mpfree(dc2);
+}
diff --git a/sys/src/libmp/port/mpleft.c b/sys/src/libmp/port/mpleft.c
new file mode 100755
index 000000000..cdcdff740
--- /dev/null
+++ b/sys/src/libmp/port/mpleft.c
@@ -0,0 +1,52 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+// res = b << shift
+void
+mpleft(mpint *b, int shift, mpint *res)
+{
+	int d, l, r, i, otop;
+	mpdigit this, last;
+
+	res->sign = b->sign;
+	if(b->top==0){
+		res->top = 0;
+		return;
+	}
+
+	// a negative left shift is a right shift
+	if(shift < 0){
+		mpright(b, -shift, res);
+		return;
+	}
+
+	// b and res may be the same so remember the old top
+	otop = b->top;
+
+	// shift
+	mpbits(res, otop*Dbits + shift);	// overkill
+	res->top = DIGITS(otop*Dbits + shift);
+	d = shift/Dbits;
+	l = shift - d*Dbits;
+	r = Dbits - l;
+
+	if(l == 0){
+		for(i = otop-1; i >= 0; i--)
+			res->p[i+d] = b->p[i];
+	} else {
+		last = 0;
+		for(i = otop-1; i >= 0; i--) {
+			this = b->p[i];
+			res->p[i+d+1] = (last<<l) | (this>>r);
+			last = this;
+		}
+		res->p[d] = last<<l;
+	}
+	for(i = 0; i < d; i++)
+		res->p[i] = 0;
+
+	// normalize
+	while(res->top > 0 && res->p[res->top-1] == 0)
+		res->top--;
+}
diff --git a/sys/src/libmp/port/mpmod.c b/sys/src/libmp/port/mpmod.c
new file mode 100755
index 000000000..91bebfa27
--- /dev/null
+++ b/sys/src/libmp/port/mpmod.c
@@ -0,0 +1,15 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+// remainder = b mod m
+//
+// knuth, vol 2, pp 398-400
+
+void
+mpmod(mpint *b, mpint *m, mpint *remainder)
+{
+	mpdiv(b, m, nil, remainder);
+	if(remainder->sign < 0)
+		mpadd(m, remainder, remainder);
+}
diff --git a/sys/src/libmp/port/mpmul.c b/sys/src/libmp/port/mpmul.c
new file mode 100755
index 000000000..dedd474a7
--- /dev/null
+++ b/sys/src/libmp/port/mpmul.c
@@ -0,0 +1,156 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+//
+//  from knuth's 1969 seminumberical algorithms, pp 233-235 and pp 258-260
+//
+//  mpvecmul is an assembly language routine that performs the inner
+//  loop.
+//
+//  the karatsuba trade off is set empiricly by measuring the algs on
+//  a 400 MHz Pentium II.
+//
+
+// karatsuba like (see knuth pg 258)
+// prereq: p is already zeroed
+static void
+mpkaratsuba(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p)
+{
+	mpdigit *t, *u0, *u1, *v0, *v1, *u0v0, *u1v1, *res, *diffprod;
+	int u0len, u1len, v0len, v1len, reslen;
+	int sign, n;
+
+	// divide each piece in half
+	n = alen/2;
+	if(alen&1)
+		n++;
+	u0len = n;
+	u1len = alen-n;
+	if(blen > n){
+		v0len = n;
+		v1len = blen-n;
+	} else {
+		v0len = blen;
+		v1len = 0;
+	}
+	u0 = a;
+	u1 = a + u0len;
+	v0 = b;
+	v1 = b + v0len;
+
+	// room for the partial products
+	t = mallocz(Dbytes*5*(2*n+1), 1);
+	if(t == nil)
+		sysfatal("mpkaratsuba: %r");
+	u0v0 = t;
+	u1v1 = t + (2*n+1);
+	diffprod = t + 2*(2*n+1);
+	res = t + 3*(2*n+1);
+	reslen = 4*n+1;
+
+	// t[0] = (u1-u0)
+	sign = 1;
+	if(mpveccmp(u1, u1len, u0, u0len) < 0){
+		sign = -1;
+		mpvecsub(u0, u0len, u1, u1len, u0v0);
+	} else
+		mpvecsub(u1, u1len, u0, u1len, u0v0);
+
+	// t[1] = (v0-v1)
+	if(mpveccmp(v0, v0len, v1, v1len) < 0){
+		sign *= -1;
+		mpvecsub(v1, v1len, v0, v1len, u1v1);
+	} else
+		mpvecsub(v0, v0len, v1, v1len, u1v1);
+
+	// t[4:5] = (u1-u0)*(v0-v1)
+	mpvecmul(u0v0, u0len, u1v1, v0len, diffprod);
+
+	// t[0:1] = u1*v1
+	memset(t, 0, 2*(2*n+1)*Dbytes);
+	if(v1len > 0)
+		mpvecmul(u1, u1len, v1, v1len, u1v1);
+
+	// t[2:3] = u0v0
+	mpvecmul(u0, u0len, v0, v0len, u0v0);
+
+	// res = u0*v0<<n + u0*v0
+	mpvecadd(res, reslen, u0v0, u0len+v0len, res);
+	mpvecadd(res+n, reslen-n, u0v0, u0len+v0len, res+n);
+
+	// res += u1*v1<<n + u1*v1<<2*n
+	if(v1len > 0){
+		mpvecadd(res+n, reslen-n, u1v1, u1len+v1len, res+n);
+		mpvecadd(res+2*n, reslen-2*n, u1v1, u1len+v1len, res+2*n);
+	}
+
+	// res += (u1-u0)*(v0-v1)<<n
+	if(sign < 0)
+		mpvecsub(res+n, reslen-n, diffprod, u0len+v0len, res+n);
+	else
+		mpvecadd(res+n, reslen-n, diffprod, u0len+v0len, res+n);
+	memmove(p, res, (alen+blen)*Dbytes);
+
+	free(t);
+}
+
+#define KARATSUBAMIN 32
+
+void
+mpvecmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p)
+{
+	int i;
+	mpdigit d;
+	mpdigit *t;
+
+	// both mpvecdigmuladd and karatsuba are fastest when a is the longer vector
+	if(alen < blen){
+		i = alen;
+		alen = blen;
+		blen = i;
+		t = a;
+		a = b;
+		b = t;
+	}
+	if(blen == 0){
+		memset(p, 0, Dbytes*(alen+blen));
+		return;
+	}
+
+	if(alen >= KARATSUBAMIN && blen > 1){
+		// O(n^1.585)
+		mpkaratsuba(a, alen, b, blen, p);
+	} else {
+		// O(n^2)
+		for(i = 0; i < blen; i++){
+			d = b[i];
+			if(d != 0)
+				mpvecdigmuladd(a, alen, d, &p[i]);
+		}
+	}
+}
+
+void
+mpmul(mpint *b1, mpint *b2, mpint *prod)
+{
+	mpint *oprod;
+
+	oprod = nil;
+	if(prod == b1 || prod == b2){
+		oprod = prod;
+		prod = mpnew(0);
+	}
+
+	prod->top = 0;
+	mpbits(prod, (b1->top+b2->top+1)*Dbits);
+	mpvecmul(b1->p, b1->top, b2->p, b2->top, prod->p);
+	prod->top = b1->top+b2->top+1;
+	prod->sign = b1->sign*b2->sign;
+	mpnorm(prod);
+
+	if(oprod != nil){
+		mpassign(prod, oprod);
+		mpfree(prod);
+	}
+}
diff --git a/sys/src/libmp/port/mprand.c b/sys/src/libmp/port/mprand.c
new file mode 100755
index 000000000..fd288f24e
--- /dev/null
+++ b/sys/src/libmp/port/mprand.c
@@ -0,0 +1,42 @@
+#include "os.h"
+#include <mp.h>
+#include <libsec.h>
+#include "dat.h"
+
+mpint*
+mprand(int bits, void (*gen)(uchar*, int), mpint *b)
+{
+	int n, m;
+	mpdigit mask;
+	uchar *p;
+
+	n = DIGITS(bits);
+	if(b == nil)
+		b = mpnew(bits);
+	else
+		mpbits(b, bits);
+
+	p = malloc(n*Dbytes);
+	if(p == nil)
+		return nil;
+	(*gen)(p, n*Dbytes);
+	betomp(p, n*Dbytes, b);
+	free(p);
+
+	// make sure we don't give too many bits
+	m = bits%Dbits;
+	n--;
+	if(m > 0){
+		mask = 1;
+		mask <<= m;
+		mask--;
+		b->p[n] &= mask;
+	}
+
+	for(; n >= 0; n--)
+		if(b->p[n] != 0)
+			break;
+	b->top = n+1;
+	b->sign = 1;
+	return b;
+}
diff --git a/sys/src/libmp/port/mpright.c b/sys/src/libmp/port/mpright.c
new file mode 100755
index 000000000..03039177b
--- /dev/null
+++ b/sys/src/libmp/port/mpright.c
@@ -0,0 +1,54 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+// res = b >> shift
+void
+mpright(mpint *b, int shift, mpint *res)
+{
+	int d, l, r, i;
+	mpdigit this, last;
+
+	res->sign = b->sign;
+	if(b->top==0){
+		res->top = 0;
+		return;
+	}
+
+	// a negative right shift is a left shift
+	if(shift < 0){
+		mpleft(b, -shift, res);
+		return;
+	}
+
+	if(res != b)
+		mpbits(res, b->top*Dbits - shift);
+	d = shift/Dbits;
+	r = shift - d*Dbits;
+	l = Dbits - r;
+
+	//  shift all the bits out == zero
+	if(d>=b->top){
+		res->top = 0;
+		return;
+	}
+
+	// special case digit shifts
+	if(r == 0){
+		for(i = 0; i < b->top-d; i++)
+			res->p[i] = b->p[i+d];
+	} else {
+		last = b->p[d];
+		for(i = 0; i < b->top-d-1; i++){
+			this = b->p[i+d+1];
+			res->p[i] = (this<<l) | (last>>r);
+			last = this;
+		}
+		res->p[i++] = last>>r;
+	}
+	while(i > 0 && res->p[i-1] == 0)
+		i--;
+	res->top = i;
+	if(i==0)
+		res->sign = 1;
+}
diff --git a/sys/src/libmp/port/mpsub.c b/sys/src/libmp/port/mpsub.c
new file mode 100755
index 000000000..3fe6ca095
--- /dev/null
+++ b/sys/src/libmp/port/mpsub.c
@@ -0,0 +1,52 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+// diff = abs(b1) - abs(b2), i.e., subtract the magnitudes
+void
+mpmagsub(mpint *b1, mpint *b2, mpint *diff)
+{
+	int n, m, sign;
+	mpint *t;
+
+	// get the sizes right
+	if(mpmagcmp(b1, b2) < 0){
+		sign = -1;
+		t = b1;
+		b1 = b2;
+		b2 = t;
+	} else
+		sign = 1;
+	n = b1->top;
+	m = b2->top;
+	if(m == 0){
+		mpassign(b1, diff);
+		diff->sign = sign;
+		return;
+	}
+	mpbits(diff, n*Dbits);
+
+	mpvecsub(b1->p, n, b2->p, m, diff->p);
+	diff->sign = sign;
+	diff->top = n;
+	mpnorm(diff);
+}
+
+// diff = b1 - b2
+void
+mpsub(mpint *b1, mpint *b2, mpint *diff)
+{
+	int sign;
+
+	if(b1->sign != b2->sign){
+		sign = b1->sign;
+		mpmagadd(b1, b2, diff);
+		diff->sign = sign;
+		return;
+	}
+
+	sign = b1->sign;
+	mpmagsub(b1, b2, diff);
+	if(diff->top != 0)
+		diff->sign *= sign;
+}
diff --git a/sys/src/libmp/port/mptobe.c b/sys/src/libmp/port/mptobe.c
new file mode 100755
index 000000000..ed527cc76
--- /dev/null
+++ b/sys/src/libmp/port/mptobe.c
@@ -0,0 +1,58 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+// convert an mpint into a big endian byte array (most significant byte first)
+//   return number of bytes converted
+//   if p == nil, allocate and result array
+int
+mptobe(mpint *b, uchar *p, uint n, uchar **pp)
+{
+	int i, j, suppress;
+	mpdigit x;
+	uchar *e, *s, c;
+
+	if(p == nil){
+		n = (b->top+1)*Dbytes;
+		p = malloc(n);
+		setmalloctag(p, getcallerpc(&b));
+	}
+	if(p == nil)
+		return -1;
+	if(pp != nil)
+		*pp = p;
+	memset(p, 0, n);
+
+	// special case 0
+	if(b->top == 0){
+		if(n < 1)
+			return -1;
+		else
+			return 1;
+	}
+		
+	s = p;
+	e = s+n;
+	suppress = 1;
+	for(i = b->top-1; i >= 0; i--){
+		x = b->p[i];
+		for(j = Dbits-8; j >= 0; j -= 8){
+			c = x>>j;
+			if(c == 0 && suppress)
+				continue;
+			if(p >= e)
+				return -1;
+			*p++ = c;
+			suppress = 0;
+		}
+	}
+
+	// guarantee at least one byte
+	if(s == p){
+		if(p >= e)
+			return -1;
+		*p++ = 0;
+	}
+
+	return p - s;
+}
diff --git a/sys/src/libmp/port/mptoi.c b/sys/src/libmp/port/mptoi.c
new file mode 100755
index 000000000..b3f22b424
--- /dev/null
+++ b/sys/src/libmp/port/mptoi.c
@@ -0,0 +1,46 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+/*
+ *  this code assumes that mpdigit is at least as
+ *  big as an int.
+ */
+
+mpint*
+itomp(int i, mpint *b)
+{
+	if(b == nil)
+		b = mpnew(0);
+	mpassign(mpzero, b);
+	if(i != 0)
+		b->top = 1;
+	if(i < 0){
+		b->sign = -1;
+		*b->p = -i;
+	} else
+		*b->p = i;
+	return b;
+}
+
+int
+mptoi(mpint *b)
+{
+	uint x;
+
+	if(b->top==0)
+		return 0;
+	x = *b->p;
+	if(b->sign > 0){
+		if(b->top > 1 || (x > MAXINT))
+			x = (int)MAXINT;
+		else
+			x = (int)x;
+	} else {
+		if(b->top > 1 || x > MAXINT+1)
+			x = (int)MININT;
+		else
+			x = -(int)x;
+	}
+	return x;
+}
diff --git a/sys/src/libmp/port/mptole.c b/sys/src/libmp/port/mptole.c
new file mode 100755
index 000000000..9421d5f66
--- /dev/null
+++ b/sys/src/libmp/port/mptole.c
@@ -0,0 +1,54 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+// convert an mpint into a little endian byte array (least significant byte first)
+
+//   return number of bytes converted
+//   if p == nil, allocate and result array
+int
+mptole(mpint *b, uchar *p, uint n, uchar **pp)
+{
+	int i, j;
+	mpdigit x;
+	uchar *e, *s;
+
+	if(p == nil){
+		n = (b->top+1)*Dbytes;
+		p = malloc(n);
+	}
+	if(pp != nil)
+		*pp = p;
+	if(p == nil)
+		return -1;
+	memset(p, 0, n);
+
+	// special case 0
+	if(b->top == 0){
+		if(n < 1)
+			return -1;
+		else
+			return 0;
+	}
+		
+	s = p;
+	e = s+n;
+	for(i = 0; i < b->top-1; i++){
+		x = b->p[i];
+		for(j = 0; j < Dbytes; j++){
+			if(p >= e)
+				return -1;
+			*p++ = x;
+			x >>= 8;
+		}
+	}
+	x = b->p[i];
+	while(x > 0){
+		if(p >= e)
+			return -1;
+		*p++ = x;
+		x >>= 8;
+	}
+
+	return p - s;
+}
diff --git a/sys/src/libmp/port/mptoui.c b/sys/src/libmp/port/mptoui.c
new file mode 100755
index 000000000..41c0b0b67
--- /dev/null
+++ b/sys/src/libmp/port/mptoui.c
@@ -0,0 +1,33 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+/*
+ *  this code assumes that mpdigit is at least as
+ *  big as an int.
+ */
+
+mpint*
+uitomp(uint i, mpint *b)
+{
+	if(b == nil)
+		b = mpnew(0);
+	mpassign(mpzero, b);
+	if(i != 0)
+		b->top = 1;
+	*b->p = i;
+	return b;
+}
+
+uint
+mptoui(mpint *b)
+{
+	uint x;
+
+	x = *b->p;
+	if(b->sign < 0)
+		x = 0;
+	else if(b->top > 1 || (sizeof(mpdigit) > sizeof(uint) && x > MAXUINT))
+		x =  MAXUINT;
+	return x;
+}
diff --git a/sys/src/libmp/port/mptouv.c b/sys/src/libmp/port/mptouv.c
new file mode 100755
index 000000000..b2a7632d1
--- /dev/null
+++ b/sys/src/libmp/port/mptouv.c
@@ -0,0 +1,49 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+#define VLDIGITS (sizeof(vlong)/sizeof(mpdigit))
+
+/*
+ *  this code assumes that a vlong is an integral number of
+ *  mpdigits long.
+ */
+mpint*
+uvtomp(uvlong v, mpint *b)
+{
+	int s;
+
+	if(b == nil)
+		b = mpnew(VLDIGITS*sizeof(mpdigit));
+	else
+		mpbits(b, VLDIGITS*sizeof(mpdigit));
+	mpassign(mpzero, b);
+	if(v == 0)
+		return b;
+	for(s = 0; s < VLDIGITS && v != 0; s++){
+		b->p[s] = v;
+		v >>= sizeof(mpdigit)*8;
+	}
+	b->top = s;
+	return b;
+}
+
+uvlong
+mptouv(mpint *b)
+{
+	uvlong v;
+	int s;
+
+	if(b->top == 0)
+		return 0LL;
+
+	mpnorm(b);
+	if(b->top > VLDIGITS)
+		return MAXVLONG;
+
+	v = 0ULL;
+	for(s = 0; s < b->top; s++)
+		v |= (uvlong)b->p[s]<<(s*sizeof(mpdigit)*8);
+
+	return v;
+}
diff --git a/sys/src/libmp/port/mptov.c b/sys/src/libmp/port/mptov.c
new file mode 100755
index 000000000..b09718ef0
--- /dev/null
+++ b/sys/src/libmp/port/mptov.c
@@ -0,0 +1,69 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+#define VLDIGITS (sizeof(vlong)/sizeof(mpdigit))
+
+/*
+ *  this code assumes that a vlong is an integral number of
+ *  mpdigits long.
+ */
+mpint*
+vtomp(vlong v, mpint *b)
+{
+	int s;
+	uvlong uv;
+
+	if(b == nil)
+		b = mpnew(VLDIGITS*sizeof(mpdigit));
+	else
+		mpbits(b, VLDIGITS*sizeof(mpdigit));
+	mpassign(mpzero, b);
+	if(v == 0)
+		return b;
+	if(v < 0){
+		b->sign = -1;
+		uv = -v;
+	} else
+		uv = v;
+	for(s = 0; s < VLDIGITS && uv != 0; s++){
+		b->p[s] = uv;
+		uv >>= sizeof(mpdigit)*8;
+	}
+	b->top = s;
+	return b;
+}
+
+vlong
+mptov(mpint *b)
+{
+	uvlong v;
+	int s;
+
+	if(b->top == 0)
+		return 0LL;
+
+	mpnorm(b);
+	if(b->top > VLDIGITS){
+		if(b->sign > 0)
+			return (vlong)MAXVLONG;
+		else
+			return (vlong)MINVLONG;
+	}
+
+	v = 0ULL;
+	for(s = 0; s < b->top; s++)
+		v |= b->p[s]<<(s*sizeof(mpdigit)*8);
+
+	if(b->sign > 0){
+		if(v > MAXVLONG)
+			v = MAXVLONG;
+	} else {
+		if(v > MINVLONG)
+			v = MINVLONG;
+		else
+			v = -(vlong)v;
+	}
+
+	return (vlong)v;
+}
diff --git a/sys/src/libmp/port/mpvecadd.c b/sys/src/libmp/port/mpvecadd.c
new file mode 100755
index 000000000..98fdcc901
--- /dev/null
+++ b/sys/src/libmp/port/mpvecadd.c
@@ -0,0 +1,35 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+// prereq: alen >= blen, sum has at least blen+1 digits
+void
+mpvecadd(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *sum)
+{
+	int i, carry;
+	mpdigit x, y;
+
+	carry = 0;
+	for(i = 0; i < blen; i++){
+		x = *a++;
+		y = *b++;
+		x += carry;
+		if(x < carry)
+			carry = 1;
+		else
+			carry = 0;
+		x += y;
+		if(x < y)
+			carry++;
+		*sum++ = x;
+	}
+	for(; i < alen; i++){
+		x = *a++ + carry;
+		if(x < carry)
+			carry = 1;
+		else
+			carry = 0;
+		*sum++ = x;
+	}
+	*sum = carry;
+}
diff --git a/sys/src/libmp/port/mpveccmp.c b/sys/src/libmp/port/mpveccmp.c
new file mode 100755
index 000000000..a462b1bd7
--- /dev/null
+++ b/sys/src/libmp/port/mpveccmp.c
@@ -0,0 +1,27 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+int
+mpveccmp(mpdigit *a, int alen, mpdigit *b, int blen)
+{
+	mpdigit x;
+
+	while(alen > blen)
+		if(a[--alen] != 0)
+			return 1;
+	while(blen > alen)
+		if(b[--blen] != 0)
+			return -1;
+	while(alen > 0){
+		--alen;
+		x = a[alen] - b[alen];
+		if(x == 0)
+			continue;
+		if(x > a[alen])
+			return -1;
+		else
+			return 1;
+	}
+	return 0;
+}
diff --git a/sys/src/libmp/port/mpvecdigmuladd.c b/sys/src/libmp/port/mpvecdigmuladd.c
new file mode 100755
index 000000000..6b6c68379
--- /dev/null
+++ b/sys/src/libmp/port/mpvecdigmuladd.c
@@ -0,0 +1,103 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+#define LO(x) ((x) & ((1<<(Dbits/2))-1))
+#define HI(x) ((x) >> (Dbits/2))
+
+static void
+mpdigmul(mpdigit a, mpdigit b, mpdigit *p)
+{
+	mpdigit x, ah, al, bh, bl, p1, p2, p3, p4;
+	int carry;
+
+	// half digits
+	ah = HI(a);
+	al = LO(a);
+	bh = HI(b);
+	bl = LO(b);
+
+	// partial products
+	p1 = ah*bl;
+	p2 = bh*al;
+	p3 = bl*al;
+	p4 = ah*bh;
+
+	// p = ((p1+p2)<<(Dbits/2)) + (p4<<Dbits) + p3
+	carry = 0;
+	x = p1<<(Dbits/2);
+	p3 += x;
+	if(p3 < x)
+		carry++;
+	x = p2<<(Dbits/2);
+	p3 += x;
+	if(p3 < x)
+		carry++;
+	p4 += carry + HI(p1) + HI(p2);	// can't carry out of the high digit
+	p[0] = p3;
+	p[1] = p4;
+}
+
+// prereq: p must have room for n+1 digits
+void
+mpvecdigmuladd(mpdigit *b, int n, mpdigit m, mpdigit *p)
+{
+	int i;
+	mpdigit carry, x, y, part[2];
+
+	carry = 0;
+	part[1] = 0;
+	for(i = 0; i < n; i++){
+		x = part[1] + carry;
+		if(x < carry)
+			carry = 1;
+		else
+			carry = 0;
+		y = *p;
+		mpdigmul(*b++, m, part);
+		x += part[0];
+		if(x < part[0])
+			carry++;
+		x += y;
+		if(x < y)
+			carry++;
+		*p++ = x;
+	}
+	*p = part[1] + carry;
+}
+
+// prereq: p must have room for n+1 digits
+int
+mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p)
+{
+	int i;
+	mpdigit x, y, part[2], borrow;
+
+	borrow = 0;
+	part[1] = 0;
+	for(i = 0; i < n; i++){
+		x = *p;
+		y = x - borrow;
+		if(y > x)
+			borrow = 1;
+		else
+			borrow = 0;
+		x = part[1];
+		mpdigmul(*b++, m, part);
+		x += part[0];
+		if(x < part[0])
+			borrow++;
+		x = y - x;
+		if(x > y)
+			borrow++;
+		*p++ = x;
+	}
+
+	x = *p;
+	y = x - borrow - part[1];
+	*p = y;
+	if(y > x)
+		return -1;
+	else
+		return 1;
+}
diff --git a/sys/src/libmp/port/mpvecsub.c b/sys/src/libmp/port/mpvecsub.c
new file mode 100755
index 000000000..db93b65bb
--- /dev/null
+++ b/sys/src/libmp/port/mpvecsub.c
@@ -0,0 +1,34 @@
+#include "os.h"
+#include <mp.h>
+#include "dat.h"
+
+// prereq: a >= b, alen >= blen, diff has at least alen digits
+void
+mpvecsub(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *diff)
+{
+	int i, borrow;
+	mpdigit x, y;
+
+	borrow = 0;
+	for(i = 0; i < blen; i++){
+		x = *a++;
+		y = *b++;
+		y += borrow;
+		if(y < borrow)
+			borrow = 1;
+		else
+			borrow = 0;
+		if(x < y)
+			borrow++;
+		*diff++ = x - y;
+	}
+	for(; i < alen; i++){
+		x = *a++;
+		y = x - borrow;
+		if(y > x)
+			borrow = 1;
+		else
+			borrow = 0;
+		*diff++ = y;
+	}
+}
diff --git a/sys/src/libmp/port/os.h b/sys/src/libmp/port/os.h
new file mode 100755
index 000000000..dae875d66
--- /dev/null
+++ b/sys/src/libmp/port/os.h
@@ -0,0 +1,3 @@
+#include <u.h>
+#include <libc.h>
+
diff --git a/sys/src/libmp/port/reduce b/sys/src/libmp/port/reduce
new file mode 100755
index 000000000..a857a28cc
--- /dev/null
+++ b/sys/src/libmp/port/reduce
@@ -0,0 +1,16 @@
+O=$1
+shift
+objtype=$1
+shift
+
+ls -p ../$objtype/*.[cs] >[2]/dev/null | sed 's/..$//' > /tmp/reduce.$pid
+#
+#	if empty directory, just return the input files
+#
+if (! ~ $status '|') {
+	echo $*
+	rm /tmp/reduce.$pid
+	exit 0
+}
+echo $* | tr ' ' \012 | grep -v -f /tmp/reduce.$pid | tr \012 ' '
+rm /tmp/reduce.$pid
diff --git a/sys/src/libmp/port/strtomp.c b/sys/src/libmp/port/strtomp.c
new file mode 100755
index 000000000..2ef8c2109
--- /dev/null
+++ b/sys/src/libmp/port/strtomp.c
@@ -0,0 +1,205 @@
+#include "os.h"
+#include <mp.h>
+#include <libsec.h>
+#include "dat.h"
+
+static struct {
+	int	inited;
+
+	uchar	t64[256];
+	uchar	t32[256];
+	uchar	t16[256];
+	uchar	t10[256];
+} tab;
+
+enum {
+	INVAL=	255
+};
+
+static char set64[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/";
+static char set32[] = "23456789abcdefghijkmnpqrstuvwxyz";
+static char set16[] = "0123456789ABCDEF0123456789abcdef";
+static char set10[] = "0123456789";
+
+static void
+init(void)
+{
+	char *p;
+
+	memset(tab.t64, INVAL, sizeof(tab.t64));
+	memset(tab.t32, INVAL, sizeof(tab.t32));
+	memset(tab.t16, INVAL, sizeof(tab.t16));
+	memset(tab.t10, INVAL, sizeof(tab.t10));
+
+	for(p = set64; *p; p++)
+		tab.t64[*p] = p-set64;
+	for(p = set32; *p; p++)
+		tab.t32[*p] = p-set32;
+	for(p = set16; *p; p++)
+		tab.t16[*p] = (p-set16)%16;
+	for(p = set10; *p; p++)
+		tab.t10[*p] = (p-set10);
+
+	tab.inited = 1;
+}
+
+static char*
+from16(char *a, mpint *b)
+{
+	char *p, *next;
+	int i;
+	mpdigit x;
+
+	b->top = 0;
+	for(p = a; *p; p++)
+		if(tab.t16[*(uchar*)p] == INVAL)
+			break;
+	mpbits(b, (p-a)*4);
+	b->top = 0;
+	next = p;
+	while(p > a){
+		x = 0;
+		for(i = 0; i < Dbits; i += 4){
+			if(p <= a)
+				break;
+			x |= tab.t16[*(uchar*)--p]<<i;
+		}
+		b->p[b->top++] = x;
+	}
+	return next;
+}
+
+static ulong mppow10[] = {
+	1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000
+};
+
+static char*
+from10(char *a, mpint *b)
+{
+	ulong x, y;
+	mpint *pow, *r;
+	int i;
+
+	pow = mpnew(0);
+	r = mpnew(0);
+
+	b->top = 0;
+	for(;;){
+		// do a billion at a time in native arithmetic
+		x = 0;
+		for(i = 0; i < 9; i++){
+			y = tab.t10[*(uchar*)a];
+			if(y == INVAL)
+				break;
+			a++;
+			x *= 10;
+			x += y;
+		}
+		if(i == 0)
+			break;
+
+		// accumulate into mpint
+		uitomp(mppow10[i], pow);
+		uitomp(x, r);
+		mpmul(b, pow, b);
+		mpadd(b, r, b);
+		if(i != 9)
+			break;
+	}
+	mpfree(pow);
+	mpfree(r);
+	return a;
+}
+
+static char*
+from64(char *a, mpint *b)
+{
+	char *buf = a;
+	uchar *p;
+	int n, m;
+
+	for(; tab.t64[*(uchar*)a] != INVAL; a++)
+		;
+	n = a-buf;
+	mpbits(b, n*6);
+	p = malloc(n);
+	if(p == nil)
+		return a;
+	m = dec64(p, n, buf, n);
+	betomp(p, m, b);
+	free(p);
+	return a;
+}
+
+static char*
+from32(char *a, mpint *b)
+{
+	char *buf = a;
+	uchar *p;
+	int n, m;
+
+	for(; tab.t64[*(uchar*)a] != INVAL; a++)
+		;
+	n = a-buf;
+	mpbits(b, n*5);
+	p = malloc(n);
+	if(p == nil)
+		return a;
+	m = dec32(p, n, buf, n);
+	betomp(p, m, b);
+	free(p);
+	return a;
+}
+
+mpint*
+strtomp(char *a, char **pp, int base, mpint *b)
+{
+	int sign;
+	char *e;
+
+	if(b == nil)
+		b = mpnew(0);
+
+	if(tab.inited == 0)
+		init();
+
+	while(*a==' ' || *a=='\t')
+		a++;
+
+	sign = 1;
+	for(;; a++){
+		switch(*a){
+		case '-':
+			sign *= -1;
+			continue;
+		}
+		break;
+	}
+
+	switch(base){
+	case 10:
+		e = from10(a, b);
+		break;
+	default:
+	case 16:
+		e = from16(a, b);
+		break;
+	case 32:
+		e = from32(a, b);
+		break;
+	case 64:
+		e = from64(a, b);
+		break;
+	}
+
+	// if no characters parsed, there wasn't a number to convert
+	if(e == a)
+		return nil;
+
+	mpnorm(b);
+	b->sign = sign;
+	if(pp != nil)
+		*pp = e;
+
+	return b;
+}
diff --git a/sys/src/libmp/power/mkfile b/sys/src/libmp/power/mkfile
new file mode 100755
index 000000000..29cc02d39
--- /dev/null
+++ b/sys/src/libmp/power/mkfile
@@ -0,0 +1,19 @@
+objtype=power
+</power/mkfile
+
+LIB=/$objtype/lib/libmp.a
+SFILES=\
+	mpvecadd.s\
+	mpvecsub.s\
+	mpvecdigmuladd.s\
+	mpvecdigmulsub.s\
+
+HFILES=/$objtype/include/u.h /sys/include/mp.h ../port/dat.h
+
+OFILES=${SFILES:%.s=%.$O}
+
+UPDATE=mkfile\
+	$HFILES\
+	$SFILES\
+
+</sys/src/cmd/mksyslib
diff --git a/sys/src/libmp/power/mpvecadd.s b/sys/src/libmp/power/mpvecadd.s
new file mode 100755
index 000000000..abbfe53d7
--- /dev/null
+++ b/sys/src/libmp/power/mpvecadd.s
@@ -0,0 +1,61 @@
+#define	BDNZ	BC	16,0,
+#define	BDNE	BC	0,2,
+
+/*
+ *	mpvecadd(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *sum)
+ *
+ *		sum[0:alen] = a[0:alen-1] + b[0:blen-1]
+ *
+ *	prereq: alen >= blen, sum has room for alen+1 digits
+ *
+ *		R3 == a	(first arg passed in R3)
+ *		R4 == alen
+ *		R5 == b
+ *		R6 == blen
+ *		R7 == sum
+ *		R8 == temporary
+ *		R9 == temporary
+ */
+TEXT	mpvecadd(SB),$-4
+
+	MOVW	alen+4(FP), R4
+	MOVW	b+8(FP), R5
+	MOVW	blen+12(FP), R6
+	MOVW	sum+16(FP), R7
+	SUB	R6, R4		/* calculate counter for second loop (alen > blen) */
+	SUB	$4, R3		/* pre decrement for MOVWU's */
+	SUB	$4, R5		/* pre decrement for MOVWU's */
+	SUB	$4, R7		/* pre decrement for MOVWU's */
+	MOVW	R0, XER		/* zero carry going in */
+
+	/* if blen == 0, don't need to add it in */
+	CMP	R0, R6
+	BEQ	_add1
+
+	/* sum[0:blen-1],carry = a[0:blen-1] + b[0:blen-1] */
+	MOVW	R6, CTR
+_addloop1:
+	MOVWU	4(R3), R8
+	MOVWU	4(R5), R9
+	ADDE	R8, R9
+	MOVWU	R9, 4(R7)
+	BDNZ	_addloop1
+
+_add1:
+	/* if alen == blen, we're done */
+	CMP	R0, R4
+	BEQ	_addend
+
+	/* sum[blen:alen-1],carry = a[blen:alen-1] + 0 + carry */
+	MOVW	R4, CTR
+_addloop2:
+	MOVWU	4(R3), R8
+	ADDE	R0, R8
+	MOVWU	R8, 4(R7)
+	BDNZ	_addloop2
+
+	/* sum[alen] = carry */
+_addend:
+	ADDE	R0, R0, R8
+	MOVW	R8, 4(R7)
+	RETURN
diff --git a/sys/src/libmp/power/mpvecdigmuladd.s b/sys/src/libmp/power/mpvecdigmuladd.s
new file mode 100755
index 000000000..b86f7ec77
--- /dev/null
+++ b/sys/src/libmp/power/mpvecdigmuladd.s
@@ -0,0 +1,56 @@
+#define	BDNZ	BC	16,0,
+#define	BDNE	BC	0,2,
+
+/*
+ *	mpvecdigmuladd(mpdigit *b, int n, mpdigit m, mpdigit *p)
+ *
+ *	p += b*m
+ *
+ *	each step looks like:
+ *		hi,lo = m*b[i]
+ *		lo += oldhi + carry
+ *		hi += carry
+ *		p[i] += lo
+ *		oldhi = hi
+ *
+ *	the registers are:
+ *		b = R3
+ *		n = R4
+ *		m = R5
+ *		p = R6
+ *		i = R7
+ *		hi = R8		- constrained by hardware
+ *		lo = R9		- constrained by hardware
+ *		oldhi = R10
+ *		tmp = R11
+ *		
+ */
+TEXT	mpvecdigmuladd(SB),$0
+
+	MOVW	n+4(FP),R4
+	MOVW	m+8(FP),R5
+	MOVW	p+12(FP),R6
+	SUB	$4, R3		/* pre decrement for MOVWU's */
+	SUB	$4, R6		/* pre decrement for MOVWU's */
+
+	MOVW	R0, R10	
+	MOVW	R0, XER
+	MOVW	R4, CTR
+_muladdloop:
+	MOVWU	4(R3),R9	/* lo = b[i] */
+	MOVW	4(R6),R11	/* tmp = p[i] */
+	MULHWU	R9,R5,R8	/* hi = (b[i] * m)>>32 */
+	MULLW	R9,R5,R9	/* lo = b[i] * m */
+	ADDC	R10,R9		/* lo += oldhi */
+	ADDE	R0,R8		/* hi += carry */
+	ADDC	R9,R11		/* tmp += lo */
+	ADDE	R0,R8		/* hi += carry */
+	MOVWU	R11,4(R6)	/* p[i] = tmp */
+	MOVW	R8,R10		/* oldhi = hi */
+	BDNZ	_muladdloop
+
+	MOVW	4(R6),R11	/* tmp = p[i] */
+	ADDC	R10,R11
+	MOVWU	R11,4(R6)	/* p[i] = tmp */
+
+	RETURN
diff --git a/sys/src/libmp/power/mpvecdigmulsub.s b/sys/src/libmp/power/mpvecdigmulsub.s
new file mode 100755
index 000000000..bfa4b91a6
--- /dev/null
+++ b/sys/src/libmp/power/mpvecdigmulsub.s
@@ -0,0 +1,66 @@
+#define	BDNZ	BC	16,0,
+#define	BDNE	BC	0,2,
+#define	BLT	BC	0xC,0,
+
+/*
+ *	mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p)
+ *
+ *	p -= b*m
+ *
+ *	each step looks like:
+ *		hi,lo = m*b[i]
+ *		lo += oldhi + carry
+ *		hi += carry
+ *		p[i] += lo
+ *		oldhi = hi
+ *
+ *	the registers are:
+ *		b = R3
+ *		n = R4
+ *		m = R5
+ *		p = R6
+ *		i = R7
+ *		hi = R8		- constrained by hardware
+ *		lo = R9		- constrained by hardware
+ *		oldhi = R10
+ *		tmp = R11
+ *		borrow = R12
+ *		
+ */
+TEXT	mpvecdigmulsub(SB),$0
+
+	MOVW	n+4(FP),R10
+	MOVW	R10,CTR
+	MOVW	m+8(FP),R5
+	MOVW	p+12(FP),R6
+	SUB	$4, R3		/* pre decrement for MOVWU's */
+	SUBC	$4, R6		/* pre decrement for MOVWU's and set carry */
+	MOVW	XER,R12
+
+	MOVW	R0, R10	
+
+_mulsubloop:
+	MOVWU	4(R3),R9	/* lo = b[i] */
+	MOVW	4(R6),R11	/* tmp = p[i] */
+	MULHWU	R9,R5,R8	/* hi = (b[i] * m)>>32 */
+	MULLW	R9,R5,R9	/* lo = b[i] * m */
+	ADDC	R10,R9		/* lo += oldhi */
+	ADDE	R0,R8		/* hi += carry */
+	MOVW	R12,XER
+	SUBE	R9,R11		/* tmp -= lo */
+	MOVW	XER,R12
+	MOVWU	R11,4(R6)	/* p[i] = tmp */
+	MOVW	R8,R10		/* oldhi = hi */
+	BDNZ	_mulsubloop
+
+	MOVW	4(R6),R11	/* tmp = p[i] */
+	MOVW	R12,XER
+	SUBE	R10,R11		/* tmp -= lo */
+	MOVWU	R11,4(R6)	/* p[i] = tmp */
+
+	/* return -1 if the result was negative, +1 otherwise */
+	SUBECC	R0,R0,R3
+	BLT	_mulsub2
+	MOVW	$1,R3
+_mulsub2:
+	RETURN
diff --git a/sys/src/libmp/power/mpvecsub.s b/sys/src/libmp/power/mpvecsub.s
new file mode 100755
index 000000000..c2a05bc18
--- /dev/null
+++ b/sys/src/libmp/power/mpvecsub.s
@@ -0,0 +1,57 @@
+#define	BDNZ	BC	16,0,
+#define	BDNE	BC	0,2,
+
+/*
+ *	mpvecsub(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *diff)
+ *
+ *		diff[0:alen-1] = a[0:alen-1] - b[0:blen-1]
+ *
+ *	prereq: alen >= blen, diff has room for alen digits
+ *
+ *		R3 == a
+ *		R4 == alen
+ *		R5 == b
+ *		R6 == blen
+ *		R7 == diff
+ *		R8 == temporary
+ *		R9 == temporary
+ */
+TEXT	mpvecsub(SB),$-4
+
+	MOVW	alen+4(FP),R4
+	MOVW	b+8(FP),R5
+	MOVW	blen+12(FP),R6
+	MOVW	diff+16(FP),R7
+	SUB	R6, R4		/* calculate counter for second loop (alen > blen) */
+	SUB	$4, R3		/* pre decrement for MOVWU's */
+	SUB	$4, R5		/* pre decrement for MOVWU's */
+	SUBC	$4, R7		/* pre decrement for MOVWU's and set carry */
+
+	/* skip subraction if b is zero */
+	CMP	R0,R6
+	BEQ	_sub1
+
+	/* diff[0:blen-1],borrow = a[0:blen-1] - b[0:blen-1] */
+	MOVW	R6, CTR
+_subloop1:
+	MOVWU	4(R3), R8
+	MOVWU	4(R5), R9
+	SUBE	R9, R8, R8
+	MOVWU	R8, 4(R7)
+	BDNZ	_subloop1
+
+_sub1:
+	/* skip subtraction if a is zero */
+	CMP	R0, R4
+	BEQ	_subend
+
+	/* diff[blen:alen-1] = a[blen:alen-1] - 0 + carry */
+	MOVW	R4, CTR
+_subloop2:
+	MOVWU	4(R3), R8
+	SUBE	R0, R8
+	MOVWU	R8, 4(R7)
+	BDNZ	_subloop2
+_subend:
+	RETURN
+
diff --git a/sys/src/libmp/power/placeholder.c b/sys/src/libmp/power/placeholder.c
new file mode 100755
index 000000000..e69de29bb
--- /dev/null
+++ b/sys/src/libmp/power/placeholder.c
diff --git a/sys/src/libmp/test.c b/sys/src/libmp/test.c
new file mode 100755
index 000000000..dbde4154f
--- /dev/null
+++ b/sys/src/libmp/test.c
@@ -0,0 +1,454 @@
+#include <u.h>
+#include <libc.h>
+#include <mp.h>
+#include "dat.h"
+
+int loops = 1;
+
+long randomreg;
+
+void
+srand(long seed)
+{
+	randomreg = seed;
+}
+
+long
+lrand(void)
+{
+	randomreg = randomreg*104381 + 81761;
+	return randomreg;
+}
+
+void
+prng(uchar *p, int n)
+{
+	while(n-- > 0)
+		*p++ = lrand();
+}
+
+void
+testconv(char *str)
+{
+	mpint *b;
+	char *p;
+
+	b = strtomp(str, nil, 16, nil);
+
+	p = mptoa(b, 10, nil, 0);
+	print("%s = ", p);
+	strtomp(p, nil, 10, b);
+	free(p);
+	print("%B\n", b);
+
+	p = mptoa(b, 16, nil, 0);
+	print("%s = ", p);
+	strtomp(p, nil, 16, b);
+	free(p);
+	print("%B\n", b);
+
+	p = mptoa(b, 32, nil, 0);
+	print("%s = ", p);
+	strtomp(p, nil, 32, b);
+	free(p);
+	print("%B\n", b);
+
+	p = mptoa(b, 64, nil, 0);
+	print("%s = ", p);
+	strtomp(p, nil, 64, b);
+	free(p);
+	print("%B\n", b);
+
+	mpfree(b);
+}
+
+void
+testshift(char *str)
+{
+	mpint *b1, *b2;
+	int i;
+
+	b1 = strtomp(str, nil, 16, nil);
+	b2 = mpnew(0);
+	for(i = 0; i < 64; i++){
+		mpleft(b1, i, b2);
+		print("%2.2d %B\n", i, b2);
+	}
+	for(i = 0; i < 64; i++){
+		mpright(b2, i, b1);
+		print("%2.2d %B\n", i, b1);
+	}
+	mpfree(b1);
+	mpfree(b2);
+}
+
+void
+testaddsub(char *str)
+{
+	mpint *b1, *b2;
+	int i;
+
+	b1 = strtomp(str, nil, 16, nil);
+	b2 = mpnew(0);
+	for(i = 0; i < 16; i++){
+		mpadd(b1, b2, b2);
+		print("%2.2d %B\n", i, b2);
+	}
+	for(i = 0; i < 16; i++){
+		mpsub(b2, b1, b2);
+		print("%2.2d %B\n", i, b2);
+	}
+	mpfree(b1);
+	mpfree(b2);
+}
+
+void
+testvecdigmuladd(char *str, mpdigit d)
+{
+	mpint *b, *b2;
+	int i;
+	vlong now;
+
+	b = strtomp(str, nil, 16, nil);
+	b2 = mpnew(0);
+
+	mpbits(b2, (b->top+1)*Dbits);
+	now = nsec();
+	for(i = 0; i < loops; i++){
+		memset(b2->p, 0, b2->top*Dbytes);
+		mpvecdigmuladd(b->p, b->top, d, b2->p);
+	}
+	if(loops > 1)
+		print("%lld ns for a %d*%d vecdigmul\n", (nsec()-now)/loops, b->top*Dbits, Dbits);
+	mpnorm(b2);
+	print("0 + %B * %ux = %B\n", b, d, b2);
+
+	mpfree(b);
+	mpfree(b2);
+}
+
+void
+testvecdigmulsub(char *str, mpdigit d)
+{
+	mpint *b, *b2;
+	int i;
+	vlong now;
+
+	b = strtomp(str, nil, 16, nil);
+	b2 = mpnew(0);
+
+	mpbits(b2, (b->top+1)*Dbits);
+	now = nsec();
+	for(i = 0; i < loops; i++){
+		memset(b2->p, 0, b2->top*Dbytes);
+		mpvecdigmulsub(b->p, b->top, d, b2->p);
+	}
+	if(loops > 1)
+		print("%lld ns for a %d*%d vecdigmul\n", (nsec()-now)/loops, b->top*Dbits, Dbits);
+	mpnorm(b2);
+	print("0 - %B * %ux = %B\n", b, d, b2);
+
+	mpfree(b);
+	mpfree(b2);
+}
+
+void
+testmul(char *str)
+{
+	mpint *b, *b1, *b2;
+	vlong now;
+	int i;
+
+	b = strtomp(str, nil, 16, nil);
+	b1 = mpcopy(b);
+	b2 = mpnew(0);
+
+	now = nsec();
+	for(i = 0; i < loops; i++)
+		mpmul(b, b1, b2);
+	if(loops > 1)
+		print("%lld µs for a %d*%d mult\n", (nsec()-now)/(loops*1000),
+			b->top*Dbits, b1->top*Dbits);
+	print("%B * %B = %B\n", b, b1, b2);
+
+	mpfree(b);
+	mpfree(b1);
+	mpfree(b2);
+}
+
+void
+testmul2(mpint *b, mpint *b1)
+{
+	mpint *b2;
+	vlong now;
+	int i;
+
+	b2 = mpnew(0);
+
+	now = nsec();
+	for(i = 0; i < loops; i++)
+		mpmul(b, b1, b2);
+	if(loops > 1)
+		print("%lld µs for a %d*%d mult\n", (nsec()-now)/(loops*1000), b->top*Dbits, b1->top*Dbits);
+	print("%B * ", b);
+	print("%B = ", b1);
+	print("%B\n", b2);
+
+	mpfree(b2);
+}
+
+void
+testdigdiv(char *str, mpdigit d)
+{
+	mpint *b;
+	mpdigit q;
+	int i;
+	vlong now;
+
+	b = strtomp(str, nil, 16, nil);
+	now = nsec();
+	for(i = 0; i < loops; i++)
+		mpdigdiv(b->p, d, &q);
+	if(loops > 1)
+		print("%lld ns for a %d / %d div\n", (nsec()-now)/loops, 2*Dbits, Dbits);
+	print("%B / %ux = %ux\n", b, d, q);
+	mpfree(b);
+}
+
+void
+testdiv(mpint *x, mpint *y)
+{
+	mpint *b2, *b3;
+	vlong now;
+	int i;
+
+	b2 = mpnew(0);
+	b3 = mpnew(0);
+	now = nsec();
+	for(i = 0; i < loops; i++)
+		mpdiv(x, y, b2, b3);
+	if(loops > 1)
+		print("%lld µs for a %d/%d div\n", (nsec()-now)/(1000*loops),
+			x->top*Dbits, y->top*Dbits);
+	print("%B / %B = %B %B\n", x, y, b2, b3);
+	mpfree(b2);
+	mpfree(b3);
+}
+
+void
+testmod(mpint *x, mpint *y)
+{
+	mpint *r;
+	vlong now;
+	int i;
+
+	r = mpnew(0);
+	now = nsec();
+	for(i = 0; i < loops; i++)
+		mpmod(x, y, r);
+	if(loops > 1)
+		print("%lld µs for a %d/%d mod\n", (nsec()-now)/(1000*loops),
+			x->top*Dbits, y->top*Dbits);
+	print("%B mod %B = %B\n", x, y, r);
+	mpfree(r);
+}
+
+void
+testinvert(mpint *x, mpint *y)
+{
+	mpint *r, *d1, *d2;
+	vlong now;
+	int i;
+
+	r = mpnew(0);
+	d1 = mpnew(0);
+	d2 = mpnew(0);
+	now = nsec();
+	mpextendedgcd(x, y, r, d1, d2);
+	mpdiv(x, r, x, d1);
+	mpdiv(y, r, y, d1);
+	for(i = 0; i < loops; i++)
+		mpinvert(x, y, r);
+	if(loops > 1)
+		print("%lld µs for a %d in %d invert\n", (nsec()-now)/(1000*loops),
+			x->top*Dbits, y->top*Dbits);
+	print("%B**-1 mod %B = %B\n", x, y, r);
+	mpmul(r, x, d1);
+	mpmod(d1, y, d2);
+	print("%B*%B mod %B = %B\n", x, r, y, d2);
+	mpfree(r);
+	mpfree(d1);
+	mpfree(d2);
+}
+
+void
+testsub1(char *a, char *b)
+{
+	mpint *b1, *b2, *b3;
+
+	b1 = strtomp(a, nil, 16, nil);
+	b2 = strtomp(b, nil, 16, nil);
+	b3 = mpnew(0);
+	mpsub(b1, b2, b3);
+	print("%B - %B = %B\n", b1, b2, b3);
+}
+
+void
+testmul1(char *a, char *b)
+{
+	mpint *b1, *b2, *b3;
+
+	b1 = strtomp(a, nil, 16, nil);
+	b2 = strtomp(b, nil, 16, nil);
+	b3 = mpnew(0);
+	mpmul(b1, b2, b3);
+	print("%B * %B = %B\n", b1, b2, b3);
+}
+
+void
+testexp(char *base, char *exp, char *mod)
+{
+	mpint *b, *e, *m, *res;
+	int i;
+	uvlong now;
+
+	b = strtomp(base, nil, 16, nil);
+	e = strtomp(exp, nil, 16, nil);
+	res = mpnew(0);
+	if(mod != nil)
+		m = strtomp(mod, nil, 16, nil);
+	else
+		m = nil;
+	now = nsec();
+	for(i = 0; i < loops; i++)
+		mpexp(b, e, m, res);
+	if(loops > 1)
+		print("%ulldµs for a %d to the %d bit exp\n", (nsec()-now)/(loops*1000),
+			b->top*Dbits, e->top*Dbits);
+	if(m != nil)
+		print("%B ^ %B mod %B == %B\n", b, e, m, res);
+	else
+		print("%B ^ %B == %B\n", b, e, res);
+	mpfree(b);
+	mpfree(e);
+	mpfree(res);
+	if(m != nil)
+		mpfree(m);
+}
+
+void
+testgcd(void)
+{
+	mpint *a, *b, *d, *x, *y, *t1, *t2;
+	int i;
+	uvlong now, then;
+	uvlong etime;
+
+	d = mpnew(0);
+	x = mpnew(0);
+	y = mpnew(0);
+	t1 = mpnew(0);
+	t2 = mpnew(0);
+
+	etime = 0;
+
+	a = strtomp("4EECAB3E04C4E6BC1F49D438731450396BF272B4D7B08F91C38E88ADCD281699889AFF872E2204C80CCAA8E460797103DE539D5DF8335A9B20C0B44886384F134C517287202FCA914D8A5096446B40CD861C641EF9C2730CB057D7B133F4C2B16BBD3D75FDDBD9151AAF0F9144AAA473AC93CF945DBFE0859FB685D5CBD0A8B3", nil, 16, nil);
+	b = strtomp("C41CFBE4D4846F67A3DF7DE9921A49D3B42DC33728427AB159CEC8CBBDB12B5F0C244F1A734AEB9840804EA3C25036AD1B61AFF3ABBC247CD4B384224567A863A6F020E7EE9795554BCD08ABAD7321AF27E1E92E3DB1C6E7E94FAAE590AE9C48F96D93D178E809401ABE8A534A1EC44359733475A36A70C7B425125062B1142D", nil, 16, nil);
+	mpextendedgcd(a, b, d, x, y);
+	print("gcd %B*%B+%B*%B = %B?\n", a, x, b, y, d);
+	mpfree(a);
+	mpfree(b);
+
+	for(i = 0; i < loops; i++){
+		a = mprand(2048, prng, nil);
+		b = mprand(2048, prng, nil);
+		then = nsec();
+		mpextendedgcd(a, b, d, x, y);
+		now = nsec();
+		etime += now-then;
+		mpmul(a, x, t1);
+		mpmul(b, y, t2);
+		mpadd(t1, t2, t2);
+		if(mpcmp(d, t2) != 0)
+			print("%d gcd %B*%B+%B*%B != %B\n", i, a, x, b, y, d);
+//		else
+//			print("%d euclid %B*%B+%B*%B == %B\n", i, a, x, b, y, d);
+		mpfree(a);
+		mpfree(b);
+	}
+
+	mpfree(x);
+	mpfree(y);
+	mpfree(d);
+	mpfree(t1);
+	mpfree(t2);
+
+	if(loops > 1)
+		print("binary %llud\n", etime);
+}
+
+extern int _unnormalizedwarning = 1;
+
+void
+main(int argc, char **argv)
+{
+	mpint *x, *y;
+
+	ARGBEGIN{
+	case 'n':
+		loops = atoi(ARGF());
+		break;
+	}ARGEND;
+
+	fmtinstall('B', mpconv);
+	fmtinstall('Q', mpconv);
+	srand(0);
+	mpsetminbits(2*Dbits);
+	testshift("1111111111111111");
+	testaddsub("fffffffffffffffff");
+	testdigdiv("1234567812345678", 0x76543218);
+	testdigdiv("1ffff", 0xffff);
+	testdigdiv("ffff", 0xffff);
+	testdigdiv("fff", 0xffff);
+	testdigdiv("effffffff", 0xffff);
+	testdigdiv("ffffffff", 0x1);
+	testdigdiv("ffffffff", 0);
+	testdigdiv("200000000", 2);
+	testdigdiv("ffffff00fffffff1", 0xfffffff1);
+	testvecdigmuladd("fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff", 2);
+	testconv("0");
+	testconv("-abc0123456789abcedf");
+	testconv("abc0123456789abcedf");
+	testvecdigmulsub("fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff", 2);
+	testsub1("1FFFFFFFE00000000", "FFFFFFFE00000001");
+	testmul1("ffffffff", "f");
+	testmul("ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff");
+	testmul1("100000000000000000000000000000000000000000000000000000002000000000000000000000000000000000000000000000000000000030000000000000000000000000000000000000000000000000000000400000000000000000000000000000000000000000000000000000004FFFFFFFFFFFFFFFE0000000200000000000000000000000000000003FFFFFFFFFFFFFFFE0000000200000000000000000000000000000002FFFFFFFFFFFFFFFE0000000200000000000000000000000000000001FFFFFFFFFFFFFFFE0000000200000000000000000000000000000000FFFFFFFFFFFFFFFE0000000200000000FFFFFFFE00000001", "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF");
+	testmul1("1000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001", "1000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001");
+	testmul1("1000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001", "1000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001");
+	testmul1("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000", "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");
+	x = mprand(256, prng, nil);
+	y = mprand(128, prng, nil);
+	testdiv(x, y);
+	x = mprand(2048, prng, nil);
+	y = mprand(1024, prng, nil);
+	testdiv(x, y);
+//	x = mprand(4*1024, prng, nil);
+//	y = mprand(4*1024, prng, nil);
+//	testmul2(x, y);
+	testsub1("677132C9", "-A26559B6");
+	testgcd();
+	x = mprand(512, prng, nil);
+	x->sign = -1;
+	y = mprand(256, prng, nil);
+	testdiv(x, y);
+	testmod(x, y);
+	x->sign = 1;
+	testinvert(y, x);
+	testexp("111111111", "222", "1000000000000000000000");
+	testexp("ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff", "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff", "100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");
+#ifdef asdf
+#endif adsf
+	print("done\n");
+	exits(0);
+}