Add libmp.

1 files changed, 112 insertions, 0 deletions

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