Add libmp.
diff --git a/src/libmp/mkfile b/src/libmp/mkfile
new file mode 100644
index 0000000..5e3998a
--- /dev/null
+++ b/src/libmp/mkfile
@@ -0,0 +1,8 @@
+PLAN9=../..
+<$PLAN9/src/mkhdr
+
+DIRS=\
+ port\
+# $systype-$objtype\
+
+<$PLAN9/src/mkdirs
diff --git a/src/libmp/port/betomp.c b/src/libmp/port/betomp.c
new file mode 100644
index 0000000..95935fc
--- /dev/null
+++ b/src/libmp/port/betomp.c
@@ -0,0 +1,40 @@
+#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);
+
+ // 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/src/libmp/port/crt.c b/src/libmp/port/crt.c
new file mode 100644
index 0000000..ddf26ed
--- /dev/null
+++ b/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/src/libmp/port/crttest.c b/src/libmp/port/crttest.c
new file mode 100644
index 0000000..b58e173
--- /dev/null
+++ b/src/libmp/port/crttest.c
@@ -0,0 +1,54 @@
+#include "os.h"
+#include <mp.h>
+#include <libsec.h>
+
+void
+testcrt(mpint **p)
+{
+ CRTpre *crt;
+ CRTres *res;
+ mpint *m, *x, *y;
+ int i;
+
+ fmtinstall('B', mpconv);
+
+ // 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("%d secs with more\n", time(0)-start);
+ exits(0);
+}
diff --git a/src/libmp/port/dat.h b/src/libmp/port/dat.h
new file mode 100644
index 0000000..7c834ac
--- /dev/null
+++ b/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/src/libmp/port/letomp.c b/src/libmp/port/letomp.c
new file mode 100644
index 0000000..e23fed2
--- /dev/null
+++ b/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/src/libmp/port/mkfile b/src/libmp/port/mkfile
new file mode 100644
index 0000000..4367091
--- /dev/null
+++ b/src/libmp/port/mkfile
@@ -0,0 +1,47 @@
+PLAN9=../../..
+<$PLAN9/src/mkhdr
+
+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\
+
+OFILES=${FILES:%=%.$O}
+
+# cull things in the per-machine directories from this list
+# OFILES= `{sh ./reduce $O $objtype $ALLOFILES}
+
+HFILES=\
+ $PLAN9/include/lib9.h\
+ $PLAN9/include/mp.h\
+ dat.h\
+
+<$PLAN9/src/mksyslib
diff --git a/src/libmp/port/mpadd.c b/src/libmp/port/mpadd.c
new file mode 100644
index 0000000..6022a64
--- /dev/null
+++ b/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/src/libmp/port/mpaux.c b/src/libmp/port/mpaux.c
new file mode 100644
index 0000000..c395d83
--- /dev/null
+++ b/src/libmp/port/mpaux.c
@@ -0,0 +1,189 @@
+#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);
+ 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/src/libmp/port/mpcmp.c b/src/libmp/port/mpcmp.c
new file mode 100644
index 0000000..a2e3cf7
--- /dev/null
+++ b/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/src/libmp/port/mpdigdiv.c b/src/libmp/port/mpdigdiv.c
new file mode 100644
index 0000000..4a73bb3
--- /dev/null
+++ b/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/src/libmp/port/mpdiv.c b/src/libmp/port/mpdiv.c
new file mode 100644
index 0000000..92aee03
--- /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);
+}
diff --git a/src/libmp/port/mpeuclid.c b/src/libmp/port/mpeuclid.c
new file mode 100644
index 0000000..80b5983
--- /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);
+}
diff --git a/src/libmp/port/mpexp.c b/src/libmp/port/mpexp.c
new file mode 100644
index 0000000..9ec067c
--- /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);
+ }
+}
diff --git a/src/libmp/port/mpextendedgcd.c b/src/libmp/port/mpextendedgcd.c
new file mode 100644
index 0000000..413a05c
--- /dev/null
+++ b/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/src/libmp/port/mpfactorial.c b/src/libmp/port/mpfactorial.c
new file mode 100644
index 0000000..01079c4
--- /dev/null
+++ b/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/src/libmp/port/mpfmt.c b/src/libmp/port/mpfmt.c
new file mode 100644
index 0000000..f7c42a7
--- /dev/null
+++ b/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/src/libmp/port/mpinvert.c b/src/libmp/port/mpinvert.c
new file mode 100644
index 0000000..ee26307
--- /dev/null
+++ b/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/src/libmp/port/mpleft.c b/src/libmp/port/mpleft.c
new file mode 100644
index 0000000..cdcdff7
--- /dev/null
+++ b/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/src/libmp/port/mpmod.c b/src/libmp/port/mpmod.c
new file mode 100644
index 0000000..91bebfa
--- /dev/null
+++ b/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/src/libmp/port/mpmul.c b/src/libmp/port/mpmul.c
new file mode 100644
index 0000000..dedd474
--- /dev/null
+++ b/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/src/libmp/port/mprand.c b/src/libmp/port/mprand.c
new file mode 100644
index 0000000..fd288f2
--- /dev/null
+++ b/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/src/libmp/port/mpright.c b/src/libmp/port/mpright.c
new file mode 100644
index 0000000..0303917
--- /dev/null
+++ b/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/src/libmp/port/mpsub.c b/src/libmp/port/mpsub.c
new file mode 100644
index 0000000..3fe6ca0
--- /dev/null
+++ b/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/src/libmp/port/mptobe.c b/src/libmp/port/mptobe.c
new file mode 100644
index 0000000..e08ae56
--- /dev/null
+++ b/src/libmp/port/mptobe.c
@@ -0,0 +1,57 @@
+#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);
+ }
+ 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/src/libmp/port/mptoi.c b/src/libmp/port/mptoi.c
new file mode 100644
index 0000000..b3f22b4
--- /dev/null
+++ b/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/src/libmp/port/mptole.c b/src/libmp/port/mptole.c
new file mode 100644
index 0000000..9421d5f
--- /dev/null
+++ b/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/src/libmp/port/mptoui.c b/src/libmp/port/mptoui.c
new file mode 100644
index 0000000..9d80c1d
--- /dev/null
+++ b/src/libmp/port/mptoui.c
@@ -0,0 +1,35 @@
+#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 || x > MAXUINT)
+ x = MAXUINT;
+ }
+ return x;
+}
diff --git a/src/libmp/port/mptouv.c b/src/libmp/port/mptouv.c
new file mode 100644
index 0000000..2548303
--- /dev/null
+++ b/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 |= b->p[s]<<(s*sizeof(mpdigit)*8);
+
+ return v;
+}
diff --git a/src/libmp/port/mptov.c b/src/libmp/port/mptov.c
new file mode 100644
index 0000000..b09718e
--- /dev/null
+++ b/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/src/libmp/port/mpvecadd.c b/src/libmp/port/mpvecadd.c
new file mode 100644
index 0000000..98fdcc9
--- /dev/null
+++ b/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/src/libmp/port/mpveccmp.c b/src/libmp/port/mpveccmp.c
new file mode 100644
index 0000000..a462b1b
--- /dev/null
+++ b/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/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;
+}
diff --git a/src/libmp/port/mpvecsub.c b/src/libmp/port/mpvecsub.c
new file mode 100644
index 0000000..db93b65
--- /dev/null
+++ b/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/src/libmp/port/os.h b/src/libmp/port/os.h
new file mode 100644
index 0000000..dae875d
--- /dev/null
+++ b/src/libmp/port/os.h
@@ -0,0 +1,3 @@
+#include <u.h>
+#include <libc.h>
+
diff --git a/src/libmp/port/reduce b/src/libmp/port/reduce
new file mode 100644
index 0000000..a857a28
--- /dev/null
+++ b/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/src/libmp/port/strtomp.c b/src/libmp/port/strtomp.c
new file mode 100644
index 0000000..8e07f00
--- /dev/null
+++ b/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[(uchar)*p] = p-set64;
+ for(p = set32; *p; p++)
+ tab.t32[(uchar)*p] = p-set32;
+ for(p = set16; *p; p++)
+ tab.t16[(uchar)*p] = (p-set16)%16;
+ for(p = set10; *p; p++)
+ tab.t10[(uchar)*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;
+}
