Add libmp.

1 files changed, 94 insertions, 0 deletions

diff --git a/src/libmp/port/mpexp.c b/src/libmp/port/mpexp.c
new file mode 100644
index 00000000..9ec067cb
--- /dev/null
+++ b/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);
+	}
+}