1 files changed, 85 insertions, 0 deletions
diff --git a/src/libmp/port/mpeuclid.c b/src/libmp/port/mpeuclid.c
new file mode 100644
index 00000000..80b5983b
--- /dev/null
+++ b/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);
+}