From b3f61791f1e9095ce8ae9c6d6415b4ee94e2f7eb Mon Sep 17 00:00:00 2001 From: rsc Date: Sun, 21 Mar 2004 14:06:38 +0000 Subject: Add libmp. --- src/libmp/port/betomp.c | 40 ++++++++ src/libmp/port/crt.c | 122 ++++++++++++++++++++++++ src/libmp/port/crttest.c | 54 +++++++++++ src/libmp/port/dat.h | 12 +++ src/libmp/port/letomp.c | 28 ++++++ src/libmp/port/mkfile | 47 +++++++++ src/libmp/port/mpadd.c | 54 +++++++++++ src/libmp/port/mpaux.c | 189 ++++++++++++++++++++++++++++++++++++ src/libmp/port/mpcmp.c | 28 ++++++ src/libmp/port/mpdigdiv.c | 43 +++++++++ src/libmp/port/mpdiv.c | 112 ++++++++++++++++++++++ src/libmp/port/mpeuclid.c | 85 +++++++++++++++++ src/libmp/port/mpexp.c | 94 ++++++++++++++++++ src/libmp/port/mpextendedgcd.c | 106 +++++++++++++++++++++ src/libmp/port/mpfactorial.c | 75 +++++++++++++++ src/libmp/port/mpfmt.c | 193 +++++++++++++++++++++++++++++++++++++ src/libmp/port/mpinvert.c | 21 ++++ src/libmp/port/mpleft.c | 52 ++++++++++ src/libmp/port/mpmod.c | 15 +++ src/libmp/port/mpmul.c | 156 ++++++++++++++++++++++++++++++ src/libmp/port/mprand.c | 42 ++++++++ src/libmp/port/mpright.c | 54 +++++++++++ src/libmp/port/mpsub.c | 52 ++++++++++ src/libmp/port/mptobe.c | 57 +++++++++++ src/libmp/port/mptoi.c | 46 +++++++++ src/libmp/port/mptole.c | 54 +++++++++++ src/libmp/port/mptoui.c | 35 +++++++ src/libmp/port/mptouv.c | 49 ++++++++++ src/libmp/port/mptov.c | 69 ++++++++++++++ src/libmp/port/mpvecadd.c | 35 +++++++ src/libmp/port/mpveccmp.c | 27 ++++++ src/libmp/port/mpvecdigmuladd.c | 103 ++++++++++++++++++++ src/libmp/port/mpvecsub.c | 34 +++++++ src/libmp/port/os.h | 3 + src/libmp/port/reduce | 16 ++++ src/libmp/port/strtomp.c | 205 ++++++++++++++++++++++++++++++++++++++++ 36 files changed, 2407 insertions(+) create mode 100644 src/libmp/port/betomp.c create mode 100644 src/libmp/port/crt.c create mode 100644 src/libmp/port/crttest.c create mode 100644 src/libmp/port/dat.h create mode 100644 src/libmp/port/letomp.c create mode 100644 src/libmp/port/mkfile create mode 100644 src/libmp/port/mpadd.c create mode 100644 src/libmp/port/mpaux.c create mode 100644 src/libmp/port/mpcmp.c create mode 100644 src/libmp/port/mpdigdiv.c create mode 100644 src/libmp/port/mpdiv.c create mode 100644 src/libmp/port/mpeuclid.c create mode 100644 src/libmp/port/mpexp.c create mode 100644 src/libmp/port/mpextendedgcd.c create mode 100644 src/libmp/port/mpfactorial.c create mode 100644 src/libmp/port/mpfmt.c create mode 100644 src/libmp/port/mpinvert.c create mode 100644 src/libmp/port/mpleft.c create mode 100644 src/libmp/port/mpmod.c create mode 100644 src/libmp/port/mpmul.c create mode 100644 src/libmp/port/mprand.c create mode 100644 src/libmp/port/mpright.c create mode 100644 src/libmp/port/mpsub.c create mode 100644 src/libmp/port/mptobe.c create mode 100644 src/libmp/port/mptoi.c create mode 100644 src/libmp/port/mptole.c create mode 100644 src/libmp/port/mptoui.c create mode 100644 src/libmp/port/mptouv.c create mode 100644 src/libmp/port/mptov.c create mode 100644 src/libmp/port/mpvecadd.c create mode 100644 src/libmp/port/mpveccmp.c create mode 100644 src/libmp/port/mpvecdigmuladd.c create mode 100644 src/libmp/port/mpvecsub.c create mode 100644 src/libmp/port/os.h create mode 100644 src/libmp/port/reduce create mode 100644 src/libmp/port/strtomp.c (limited to 'src/libmp/port') diff --git a/src/libmp/port/betomp.c b/src/libmp/port/betomp.c new file mode 100644 index 00000000..95935fcd --- /dev/null +++ b/src/libmp/port/betomp.c @@ -0,0 +1,40 @@ +#include "os.h" +#include +#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 00000000..ddf26ed5 --- /dev/null +++ b/src/libmp/port/crt.c @@ -0,0 +1,122 @@ +#include "os.h" +#include +#include + +// 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 00000000..b58e1736 --- /dev/null +++ b/src/libmp/port/crttest.c @@ -0,0 +1,54 @@ +#include "os.h" +#include +#include + +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 00000000..7c834ac8 --- /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 00000000..e23fed21 --- /dev/null +++ b/src/libmp/port/letomp.c @@ -0,0 +1,28 @@ +#include "os.h" +#include +#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 00000000..43670918 --- /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 00000000..6022a64e --- /dev/null +++ b/src/libmp/port/mpadd.c @@ -0,0 +1,54 @@ +#include "os.h" +#include +#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 00000000..c395d837 --- /dev/null +++ b/src/libmp/port/mpaux.c @@ -0,0 +1,189 @@ +#include "os.h" +#include +#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)<top==0) + return 0; + k = 0; + bit = 0; + digit = 0; + d = n->p[0]; + for(;;){ + if(d & (1<= 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 00000000..a2e3cf72 --- /dev/null +++ b/src/libmp/port/mpcmp.c @@ -0,0 +1,28 @@ +#include "os.h" +#include +#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 00000000..4a73bb3a --- /dev/null +++ b/src/libmp/port/mpdigdiv.c @@ -0,0 +1,43 @@ +#include "os.h" +#include +#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< lo) + continue; + if(y > lo) + hi--; + lo -= y; + hi -= x; + q |= 1< +#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 00000000..80b5983b --- /dev/null +++ b/src/libmp/port/mpeuclid.c @@ -0,0 +1,85 @@ +#include "os.h" +#include + +// 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 00000000..9ec067cb --- /dev/null +++ b/src/libmp/port/mpexp.c @@ -0,0 +1,94 @@ +#include "os.h" +#include +#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 00000000..413a05c2 --- /dev/null +++ b/src/libmp/port/mpextendedgcd.c @@ -0,0 +1,106 @@ +#include "os.h" +#include + +#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 00000000..01079c44 --- /dev/null +++ b/src/libmp/port/mpfactorial.c @@ -0,0 +1,75 @@ +#include "os.h" +#include +#include +#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 00000000..f7c42a7b --- /dev/null +++ b/src/libmp/port/mpfmt.c @@ -0,0 +1,193 @@ +#include "os.h" +#include +#include +#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 00000000..ee263070 --- /dev/null +++ b/src/libmp/port/mpinvert.c @@ -0,0 +1,21 @@ +#include "os.h" +#include + +#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 00000000..cdcdff74 --- /dev/null +++ b/src/libmp/port/mpleft.c @@ -0,0 +1,52 @@ +#include "os.h" +#include +#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<>r); + last = this; + } + res->p[d] = last<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 00000000..91bebfa2 --- /dev/null +++ b/src/libmp/port/mpmod.c @@ -0,0 +1,15 @@ +#include "os.h" +#include +#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 00000000..dedd474a --- /dev/null +++ b/src/libmp/port/mpmul.c @@ -0,0 +1,156 @@ +#include "os.h" +#include +#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< 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)<= 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 00000000..fd288f24 --- /dev/null +++ b/src/libmp/port/mprand.c @@ -0,0 +1,42 @@ +#include "os.h" +#include +#include +#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 00000000..03039177 --- /dev/null +++ b/src/libmp/port/mpright.c @@ -0,0 +1,54 @@ +#include "os.h" +#include +#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<>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 00000000..3fe6ca09 --- /dev/null +++ b/src/libmp/port/mpsub.c @@ -0,0 +1,52 @@ +#include "os.h" +#include +#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 00000000..e08ae56b --- /dev/null +++ b/src/libmp/port/mptobe.c @@ -0,0 +1,57 @@ +#include "os.h" +#include +#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 00000000..b3f22b42 --- /dev/null +++ b/src/libmp/port/mptoi.c @@ -0,0 +1,46 @@ +#include "os.h" +#include +#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 00000000..9421d5f6 --- /dev/null +++ b/src/libmp/port/mptole.c @@ -0,0 +1,54 @@ +#include "os.h" +#include +#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 00000000..9d80c1df --- /dev/null +++ b/src/libmp/port/mptoui.c @@ -0,0 +1,35 @@ +#include "os.h" +#include +#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 00000000..2548303c --- /dev/null +++ b/src/libmp/port/mptouv.c @@ -0,0 +1,49 @@ +#include "os.h" +#include +#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 00000000..b09718ef --- /dev/null +++ b/src/libmp/port/mptov.c @@ -0,0 +1,69 @@ +#include "os.h" +#include +#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 00000000..98fdcc90 --- /dev/null +++ b/src/libmp/port/mpvecadd.c @@ -0,0 +1,35 @@ +#include "os.h" +#include +#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 00000000..a462b1bd --- /dev/null +++ b/src/libmp/port/mpveccmp.c @@ -0,0 +1,27 @@ +#include "os.h" +#include +#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 00000000..6b6c6837 --- /dev/null +++ b/src/libmp/port/mpvecdigmuladd.c @@ -0,0 +1,103 @@ +#include "os.h" +#include +#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< 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 00000000..db93b65b --- /dev/null +++ b/src/libmp/port/mpvecsub.c @@ -0,0 +1,34 @@ +#include "os.h" +#include +#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 00000000..dae875d6 --- /dev/null +++ b/src/libmp/port/os.h @@ -0,0 +1,3 @@ +#include +#include + diff --git a/src/libmp/port/reduce b/src/libmp/port/reduce new file mode 100644 index 00000000..a857a28c --- /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 00000000..8e07f009 --- /dev/null +++ b/src/libmp/port/strtomp.c @@ -0,0 +1,205 @@ +#include "os.h" +#include +#include +#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]<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; +} -- cgit v1.2.3