Add libmp.
diff --git a/src/libmp/port/mpvecdigmuladd.c b/src/libmp/port/mpvecdigmuladd.c
new file mode 100644
index 0000000..6b6c683
--- /dev/null
+++ b/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;
+}
