shithub: femtolisp

Info  •  Files  •  Log

don't do recursive makefiling

diff: cannot open b/mp//null: file does not exist: 'b/mp//null'

--- a/.gitignore

+++ b/.gitignore

@@ -1,12 +1,8 @@

-/*.[05678qvt]

-/*.o

-/*.out

-/*.a

+*.[05678qvt]

+*.o

+*.out

+*.a

 /flisp

-/llt/*.[05678qvt]

-/llt/*.o

-/llt/*.out

-/llt/*.a

 /flisp.boot.bak

 /flisp.boot.new

 boot.h

--- a/Makefile

+++ b/Makefile

@@ -3,7 +3,6 @@

 BIN=${DESTDIR}${PREFIX}/bin

 TARG=flisp

-LLT=llt/libllt.a

 CFLAGS?=-O2 -g

 CFLAGS+=-Wall -Wextra -Wno-parentheses -std=c99

 LDFLAGS?=

@@ -16,6 +15,45 @@

 	equalhash.o\

 	table.o\

 	iostream.o\

+	llt/bitvector-ops.o\

+	llt/bitvector.o\

+	llt/dump.o\

+	llt/hashing.o\

+	llt/htable.o\

+	llt/int2str.o\

+	llt/ios.o\

+	llt/lltinit.o\

+	llt/ptrhash.o\

+	llt/random.o\

+	llt/timefuncs.o\

+	llt/utf8.o\

+	mp/mpadd.o\

+	mp/mpaux.o\

+	mp/mpcmp.o\

+	mp/mpdigdiv.o\

+	mp/mpdiv.o\

+	mp/mpfmt.o\

+	mp/mpleft.o\

+	mp/mplogic.o\

+	mp/mpmul.o\

+	mp/mpright.o\

+	mp/mpsub.o\

+	mp/mptobe.o\

+	mp/mptober.o\

+	mp/mptod.o\

+	mp/mptoi.o\

+	mp/mptoui.o\

+	mp/mptouv.o\

+	mp/mptov.o\

+	mp/mpvecadd.o\

+	mp/mpveccmp.o\

+	mp/mpvecdigmuladd.o\

+	mp/mpvecsub.o\

+	mp/mpvectscmp.o\

+	mp/strtomp.o\

+	mp/u16.o\

+	mp/u32.o\

+	mp/u64.o\

 .PHONY: all default test clean

@@ -26,8 +64,8 @@

 test: ${TARG}

 	cd test && ../$(TARG) unittest.lsp

-${TARG}: ${OBJS} ${LLT}

-	${CC} -o $@ ${OBJS} ${LDFLAGS} ${LLT} -lm

+${TARG}: ${OBJS}

+	${CC} -o $@ ${OBJS} ${LDFLAGS} -lm

 .SUFFIXES: .c .o

 .c.o:

@@ -42,9 +80,5 @@

 builtin_fns.h:

 	sed -nE 's/^BUILTIN[_]?(\(".*)/BUILTIN_FN\1/gp' *.c >$@

-${LLT}:

-	${MAKE} -C llt CFLAGS="${CFLAGS} -I../posix" CC="${CC}"

-

 clean:

-	rm -f *.o ${TARG}

-	${MAKE} -C llt clean

+	rm -f ${OBJS} ${TARG}

--- a/llt/Makefile

+++ /dev/null

@@ -1,56 +1,0 @@

-CFLAGS?=-O2 -pipe -g -Wall -falign-functions -Wno-strict-aliasing

-TARG=libllt.a

-

-OBJS=\

-	bitvector-ops.o\

-	bitvector.o\

-	dump.o\

-	hashing.o\

-	htable.o\

-	int2str.o\

-	ios.o\

-	lltinit.o\

-	ptrhash.o\

-	random.o\

-	timefuncs.o\

-	utf8.o\

-    \

-	mpadd.o\

-	mpaux.o\

-	mpcmp.o\

-	mpdigdiv.o\

-	mpdiv.o\

-	mpfmt.o\

-	mpleft.o\

-	mplogic.o\

-	mpmul.o\

-	mpright.o\

-	mpsub.o\

-	mptobe.o\

-	mptod.o\

-	mptoi.o\

-	mptoui.o\

-	mptouv.o\

-	mptov.o\

-	mpvecadd.o\

-	mpvecdigmuladd.o\

-	mpvecsub.o\

-	strtomp.o\

-	u16.o\

-	u32.o\

-	u64.o\

-    mptober.o\

-    mpveccmp.o\

-    mpvectscmp.o\

-

-.PHONY: all default clean

-

-all: default

-

-default: ${TARG}

-

-clean:

-	rm -f *.o ${TARG}

-

-${TARG}: ${OBJS}

-	${AR} crs ${TARG} ${OBJS}

--- a/llt/mkfile

+++ /dev/null

@@ -1,30 +1,0 @@

-</$objtype/mkfile

-LIB=libllt.$O.a

-

-CFLAGS=$CFLAGS -p -I../plan9 -D__plan9__ -D__${objtype}__

-

-OFILES=\

-	bitvector-ops.$O\

-	bitvector.$O\

-	dump.$O\

-	hashing.$O\

-	htable.$O\

-	int2str.$O\

-	ios.$O\

-	lltinit.$O\

-	ptrhash.$O\

-	random.$O\

-	timefuncs.$O\

-	utf8.$O\

-	wcwidth.$O\

-

-HFILES=\

-	../plan9/platform.h\

-    bitvector.h\

-    htable.h\

-    ieee754.h\

-    ios.h\

-    llt.h\

-    utf8.h\

-

-</sys/src/cmd/mklib

--- a/llt/mpadd.c

+++ /dev/null

@@ -1,56 +1,0 @@

-#include "platform.h"

-

-// sum = abs(b1) + abs(b2), i.e., add the magnitudes

-void

-mpmagadd(mpint *b1, mpint *b2, mpint *sum)

-{

-	int m, n;

-	mpint *t;

-

-	sum->flags |= (b1->flags | b2->flags) & MPtimesafe;

-

-	// 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);

-		sum->sign = 1;

-		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){

-		assert(((b1->flags | b2->flags | sum->flags) & MPtimesafe) == 0);

-		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;

-	}

-}

--- a/llt/mpaux.c

+++ /dev/null

@@ -1,201 +1,0 @@

-#include "platform.h"

-

-static mpdigit _mptwodata[1] = { 2 };

-static mpint _mptwo =

-{

-	1, 1, 1,

-	_mptwodata,

-	MPstatic|MPnorm

-};

-mpint *mptwo = &_mptwo;

-

-static mpdigit _mponedata[1] = { 1 };

-static mpint _mpone =

-{

-	1, 1, 1,

-	_mponedata,

-	MPstatic|MPnorm

-};

-mpint *mpone = &_mpone;

-

-static mpdigit _mpzerodata[1] = { 0 };

-static mpint _mpzero =

-{

-	1, 1, 0,

-	_mpzerodata,

-	MPstatic|MPnorm

-};

-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");

-

-	n = DIGITS(n);

-	if(n < mpmindigits)

-		n = mpmindigits;

-	b = calloc(1, sizeof(mpint) + n*Dbytes);

-	if(b == nil)

-		sysfatal("mpnew: %r");

-	b->p = (mpdigit*)&b[1];

-	b->size = n;

-	b->sign = 1;

-	b->flags = MPnorm;

-

-	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;

-	} else {

-		if(b->p == (mpdigit*)&b[1]){

-			b->p = (mpdigit*)malloc(n*Dbytes);

-			if(b->p == nil)

-				sysfatal("mpbits: %r");

-			memmove(b->p, &b[1], Dbytes*b->top);

-			memset(&b[1], 0, Dbytes*b->size);

-		} else {

-			b->p = (mpdigit*)realloc(b->p, n*Dbytes);

-			if(b->p == nil)

-				sysfatal("mpbits: %r");

-		}

-		b->size = n;

-	}

-	memset(&b->p[b->top], 0, Dbytes*(n - b->top));

-	b->top = n;

-	b->flags &= ~MPnorm;

-}

-

-void

-mpfree(mpint *b)

-{

-	if(b == nil)

-		return;

-	if(b->flags & MPstatic)

-		sysfatal("freeing mp constant");

-	memset(b->p, 0, b->size*Dbytes);

-	if(b->p != (mpdigit*)&b[1])

-		free(b->p);

-	free(b);

-}

-

-mpint*

-mpnorm(mpint *b)

-{

-	int i;

-

-	if(b->flags & MPtimesafe){

-		assert(b->sign == 1);

-		b->flags &= ~MPnorm;

-		return b;

-	}

-	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;

-	b->flags |= MPnorm;

-	return b;

-}

-

-mpint*

-mpcopy(mpint *old)

-{

-	mpint *new;

-

-	new = mpnew(Dbits*old->size);

-	new->sign = old->sign;

-	new->top = old->top;

-	new->flags = old->flags & ~(MPstatic|MPfield);

-	memmove(new->p, old->p, Dbytes*old->top);

-	return new;

-}

-

-void

-mpassign(mpint *old, mpint *new)

-{

-	if(new == nil || old == new)

-		return;

-	new->top = 0;

-	mpbits(new, Dbits*old->top);

-	new->sign = old->sign;

-	new->top = old->top;

-	new->flags &= ~MPnorm;

-	new->flags |= old->flags & ~(MPstatic|MPfield);

-	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;

-

-	assert(n->flags & MPnorm);

-	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;

-}

--- a/llt/mpcmp.c

+++ /dev/null

@@ -1,28 +1,0 @@

-#include "platform.h"

-

-// return neg, 0, pos as abs(b1)-abs(b2) is neg, 0, pos

-int

-mpmagcmp(mpint *b1, mpint *b2)

-{

-	int i;

-

-	i = b1->flags | b2->flags;

-	if(i & MPtimesafe)

-		return mpvectscmp(b1->p, b1->top, b2->p, b2->top);

-	if(i & MPnorm){

-		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)

-{

-	int sign;

-

-	sign = (b1->sign - b2->sign) >> 1;	// -1, 0, 1

-	return sign | (sign&1)-1 & mpmagcmp(b1, b2)*b1->sign;

-}

--- a/llt/mpdigdiv.c

+++ /dev/null

@@ -1,54 +1,0 @@

-#include "platform.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;

-	}

-

-	// very common case

-	if(~divisor == 0){

-		lo += hi;

-		if(lo < hi){

-			hi++;

-			lo++;

-		}

-		if(lo+1 == 0)

-			hi++;

-		*quotient = hi;

-		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 |= 1U<<i;

-	}

-	q += lo/divisor;

-	*quotient = q;

-}

--- a/llt/mpdiv.c

+++ /dev/null

@@ -1,140 +1,0 @@

-#include "platform.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, qsign, rsign;

-	mpdigit qd, *up, *vp, *qp;

-	mpint *u, *v, *t;

-

-	assert(quotient != remainder);

-	assert(divisor->flags & MPnorm);

-

-	// divide bv zero

-	if(divisor->top == 0)

-		abort();

-

-	// division by one or small powers of two

-	if(divisor->top == 1 && (divisor->p[0] & divisor->p[0]-1) == 0){

-		vlong r = 0;

-		if(dividend->top > 0)

-			r = (vlong)dividend->sign * (dividend->p[0] & divisor->p[0]-1);

-		if(quotient != nil){

-			sign = divisor->sign;

-			for(s = 0; ((divisor->p[0] >> s) & 1) == 0; s++)

-				;

-			mpright(dividend, s, quotient);

-			if(sign < 0)

-				quotient->sign ^= (-mpmagcmp(quotient, mpzero) >> 31) << 1;

-		}

-		if(remainder != nil){

-			remainder->flags |= dividend->flags & MPtimesafe;

-			vtomp(r, remainder);

-		}

-		return;

-	}

-	assert((dividend->flags & MPtimesafe) == 0);

-

-	// quick check

-	if(mpmagcmp(dividend, divisor) < 0){

-		if(remainder != nil)

-			mpassign(dividend, remainder);

-		if(quotient != nil)

-			mpassign(mpzero, quotient);

-		return;

-	}

-	

-	qsign = divisor->sign * dividend->sign;

-	rsign = dividend->sign;

-

-	// 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){

-		assert((quotient->flags & MPtimesafe) == 0);

-		mpnorm(quotient);

-		if(quotient->top != 0)

-			quotient->sign = qsign;

-	}

-

-	if(remainder != nil){

-		assert((remainder->flags & MPtimesafe) == 0);

-		mpright(u, s, remainder);	// u is the remainder shifted

-		if(remainder->top != 0)

-			remainder->sign = rsign;

-	}

-

-	mpfree(t);

-	mpfree(u);

-	if(v != divisor)

-		mpfree(v);

-}

--- a/llt/mpfmt.c

+++ /dev/null

@@ -1,207 +1,0 @@

-#include "platform.h"

-

-static int

-toencx(mpint *b, char *buf, int len, int (*enc)(char*, int, uchar*, int))

-{

-	uchar *p;

-	int n, rv;

-

-	p = nil;

-	n = mptobe(b, nil, 0, &p);

-	if(n < 0)

-		return -1;

-	rv = (*enc)(buf, len, p, n);

-	free(p);

-	return rv;

-}

-

-static int

-topow2(mpint *b, char *buf, int len, int s)

-{

-	mpdigit *p, x;

-	int i, j, sn;

-	char *out, *eout;

-

-	if(len < 1)

-		return -1;

-

-	sn = 1<<s;

-	out = buf;

-	eout = buf+len;

-	for(p = &b->p[b->top-1]; p >= b->p; p--){

-		x = *p;

-		for(i = Dbits-s; i >= 0; i -= s){

-			j = x >> i & sn - 1;

-			if(j != 0 || out != buf){

-				if(out >= eout)

-					return -1;

-				*out++ = enc16chr(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);

-	d->flags &= ~MPtimesafe;

-	mpnorm(d);

-	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;

-}

-

-static int

-to8(mpint *b, char *buf, int len)

-{

-	mpdigit x, y;

-	char *out;

-	int i, j;

-

-	if(len < 2)

-		return -1;

-

-	out = buf+len;

-	*--out = 0;

-

-	i = j = 0;

-	x = y = 0;

-	while(j < b->top){

-		y = b->p[j++];

-		if(i > 0)

-			x |= y << i;

-		else

-			x = y;

-		i += Dbits;

-		while(i >= 3){

-Digout:			i -= 3;

-			if(out > buf)

-				out--;

-			else if(x != 0)

-				return -1;

-			*out = '0' + (x & 7);

-			x = y >> (Dbits-i);

-		}

-	}

-	if(i > 0)

-		goto Digout;

-

-	while(*out == '0') out++;

-	if(*out == '\0')

-		*--out = '0';

-

-	len -= out-buf;

-	if(out != buf)

-		memmove(buf, out, len);

-	return 0;

-}

-

-char*

-mptoa(mpint *b, int base, char *buf, int len)

-{

-	char *out;

-	int rv, alloced;

-

-	if(base == 0)

-		base = 16;	/* default */

-	alloced = 0;

-	if(buf == nil){

-		/* rv <= log₂(base) */

-		for(rv=1; (base >> rv) > 1; rv++)

-			;

-		len = 10 + (b->top*Dbits / rv);

-		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 = toencx(b, out, len, enc64);

-		break;

-	case 32:

-		rv = toencx(b, out, len, enc32);

-		break;

-	case 16:

-		rv = topow2(b, out, len, 4);

-		break;

-	case 10:

-		rv = to10(b, out, len);

-		break;

-	case 8:

-		rv = to8(b, out, len);

-		break;

-	case 4:

-		rv = topow2(b, out, len, 2);

-		break;

-	case 2:

-		rv = topow2(b, out, len, 1);

-		break;

-	default:

-		abort();

-		return nil;

-	}

-	if(rv < 0){

-		if(alloced)

-			free(buf);

-		return nil;

-	}

-	return buf;

-}

--- a/llt/mpleft.c

+++ /dev/null

@@ -1,49 +1,0 @@

-#include "platform.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 zero or 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;

-

-	res->flags |= b->flags & MPtimesafe;

-	mpnorm(res);

-}

--- a/llt/mplogic.c

+++ /dev/null

@@ -1,210 +1,0 @@

-#include "platform.h"

-

-/*

-	mplogic calculates b1|b2 subject to the

-	following flag bits (fl)

-

-	bit 0: subtract 1 from b1

-	bit 1: invert b1

-	bit 2: subtract 1 from b2

-	bit 3: invert b2

-	bit 4: add 1 to output

-	bit 5: invert output

-	

-	it inverts appropriate bits automatically

-	depending on the signs of the inputs

-*/

-

-static void

-mplogic(mpint *b1, mpint *b2, mpint *sum, int fl)

-{

-	mpint *t;

-	mpdigit *dp1, *dp2, *dpo, d1, d2, d;

-	int c1, c2, co;

-	int i;

-

-	assert(((b1->flags | b2->flags | sum->flags) & MPtimesafe) == 0);

-	if(b1->sign < 0) fl ^= 0x03;

-	if(b2->sign < 0) fl ^= 0x0c;

-	sum->sign = (int)(((fl|fl>>2)^fl>>4)<<30)>>31|1;

-	if(sum->sign < 0) fl ^= 0x30;

-	if(b2->top > b1->top){

-		t = b1;

-		b1 = b2;

-		b2 = t;

-		fl = fl >> 2 & 0x03 | fl << 2 & 0x0c | fl & 0x30;

-	}

-	mpbits(sum, b1->top*Dbits+1);

-	dp1 = b1->p;

-	dp2 = b2->p;

-	dpo = sum->p;

-	c1 = fl & 1;

-	c2 = fl >> 2 & 1;

-	co = fl >> 4 & 1;

-	for(i = 0; i < b1->top; i++){

-		d1 = dp1[i] - c1;

-		if(i < b2->top)

-			d2 = dp2[i] - c2;

-		else

-			d2 = 0;

-		if(d1 != (mpdigit)-1) c1 = 0;

-		if(d2 != (mpdigit)-1) c2 = 0;

-		if((fl & 2) != 0) d1 ^= -1;

-		if((fl & 8) != 0) d2 ^= -1;

-		d = d1 | d2;

-		if((fl & 32) != 0) d ^= -1;

-		d += co;

-		if(d != 0) co = 0;

-		dpo[i] = d;

-	}

-	sum->top = i;

-	if(co)

-		dpo[sum->top++] = co;

-	mpnorm(sum);

-}

-

-void

-mpor(mpint *b1, mpint *b2, mpint *sum)

-{

-	mplogic(b1, b2, sum, 0);

-}

-

-void

-mpand(mpint *b1, mpint *b2, mpint *sum)

-{

-	mplogic(b1, b2, sum, 0x2a);

-}

-

-void

-mpbic(mpint *b1, mpint *b2, mpint *sum)

-{

-	mplogic(b1, b2, sum, 0x22);

-}

-

-void

-mpnot(mpint *b, mpint *r)

-{

-	mpadd(b, mpone, r);

-	if(r->top != 0)

-		r->sign ^= -2;

-}

-

-void

-mpxor(mpint *b1, mpint *b2, mpint *sum)

-{

-	mpint *t;

-	mpdigit *dp1, *dp2, *dpo, d1, d2, d;

-	int c1, c2, co;

-	int i, fl;

-

-	assert(((b1->flags | b2->flags | sum->flags) & MPtimesafe) == 0);

-	if(b2->top > b1->top){

-		t = b1;

-		b1 = b2;

-		b2 = t;

-	}

-	fl = (b1->sign & 10) ^ (b2->sign & 12);

-	sum->sign = (int)(fl << 28) >> 31 | 1;

-	mpbits(sum, b1->top*Dbits+1);

-	dp1 = b1->p;

-	dp2 = b2->p;

-	dpo = sum->p;

-	c1 = fl >> 1 & 1;

-	c2 = fl >> 2 & 1;

-	co = fl >> 3 & 1;

-	for(i = 0; i < b1->top; i++){

-		d1 = dp1[i] - c1;

-		if(i < b2->top)

-			d2 = dp2[i] - c2;

-		else

-			d2 = 0;

-		if(d1 != (mpdigit)-1) c1 = 0;

-		if(d2 != (mpdigit)-1) c2 = 0;

-		d = d1 ^ d2;

-		d += co;

-		if(d != 0) co = 0;

-		dpo[i] = d;

-	}

-	sum->top = i;

-	if(co)

-		dpo[sum->top++] = co;

-	mpnorm(sum);

-}

-

-void

-mptrunc(mpint *b, int n, mpint *r)

-{

-	int d, m, i, c;

-

-	assert(((b->flags | r->flags) & MPtimesafe) == 0);

-	mpbits(r, n);

-	r->top = DIGITS(n);

-	d = n / Dbits;

-	m = n % Dbits;

-	if(b->sign == -1){

-		c = 1;

-		for(i = 0; i < r->top; i++){

-			if(i < b->top)

-				r->p[i] = ~(b->p[i] - c);

-			else

-				r->p[i] = -1;

-			if(r->p[i] != 0)

-				c = 0;

-		}

-		if(m != 0)

-			r->p[d] &= (1<<m) - 1;

-	}else if(b->sign == 1){

-		if(d >= b->top){

-			mpassign(b, r);

-			mpnorm(r);

-			return;

-		}

-		if(b != r)

-			for(i = 0; i < d; i++)

-				r->p[i] = b->p[i];

-		if(m != 0)

-			r->p[d] = b->p[d] & (1<<m)-1;

-	}

-	r->sign = 1;

-	mpnorm(r);

-}

-

-void

-mpxtend(mpint *b, int n, mpint *r)

-{

-	int d, m, c, i;

-

-	d = (n - 1) / Dbits;

-	m = (n - 1) % Dbits;

-	if(d >= b->top){

-		mpassign(b, r);

-		return;

-	}

-	mptrunc(b, n, r);

-	mpbits(r, n);

-	if((r->p[d] & 1<<m) == 0){

-		mpnorm(r);

-		return;

-	}

-	r->p[d] |= -(1<<m);

-	r->sign = -1;

-	c = 1;

-	for(i = 0; i < r->top; i++){

-		r->p[i] = ~(r->p[i] - c);

-		if(r->p[i] != 0)

-			c = 0;

-	}

-	mpnorm(r);

-}

-

-void

-mpasr(mpint *b, int n, mpint *r)

-{

-	if(b->sign > 0 || n <= 0){

-		mpright(b, n, r);

-		return;

-	}

-	mpadd(b, mpone, r);

-	mpright(r, n, r);

-	mpsub(r, mpone, r);

-}

--- a/llt/mpmul.c

+++ /dev/null

@@ -1,174 +1,0 @@

-#include "platform.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 = calloc(1, Dbytes*5*(2*n+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(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

-mpvectsmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p)

-{

-	int i;

-	mpdigit *t;

-

-	if(alen < blen){

-		i = alen;

-		alen = blen;

-		blen = i;

-		t = a;

-		a = b;

-		b = t;

-	}

-	if(blen == 0)

-		return;

-	for(i = 0; i < blen; i++)

-		mpvecdigmuladd(a, alen, b[i], &p[i]);

-}

-

-void

-mpmul(mpint *b1, mpint *b2, mpint *prod)

-{

-	mpint *oprod;

-

-	oprod = prod;

-	if(prod == b1 || prod == b2){

-		prod = mpnew(0);

-		prod->flags = oprod->flags;

-	}

-	prod->flags |= (b1->flags | b2->flags) & MPtimesafe;

-

-	prod->top = 0;

-	mpbits(prod, (b1->top+b2->top+1)*Dbits);

-	if(prod->flags & MPtimesafe)

-		mpvectsmul(b1->p, b1->top, b2->p, b2->top, prod->p);

-	else

-		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 != prod){

-		mpassign(prod, oprod);

-		mpfree(prod);

-	}

-}

--- a/llt/mpright.c

+++ /dev/null

@@ -1,55 +1,0 @@

-#include "platform.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);

-	else if(shift == 0)

-		return;

-

-	d = shift/Dbits;

-	r = shift - d*Dbits;

-	l = Dbits - r;

-

-	//  shift all the bits out == zero

-	if(d>=b->top){

-		res->sign = 1;

-		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;

-	}

-

-	res->top = i;

-	res->flags |= b->flags & MPtimesafe;

-	mpnorm(res);

-}

--- a/llt/mpsub.c

+++ /dev/null

@@ -1,54 +1,0 @@

-#include "platform.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){

-		assert(((b1->flags | b2->flags | diff->flags) & MPtimesafe) == 0);

-		sign = -1;

-		t = b1;

-		b1 = b2;

-		b2 = t;

-	} else {

-		diff->flags |= (b1->flags | b2->flags) & MPtimesafe;

-		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){

-		assert(((b1->flags | b2->flags | diff->flags) & MPtimesafe) == 0);

-		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;

-}

--- a/llt/mptobe.c

+++ /dev/null

@@ -1,29 +1,0 @@

-#include "platform.h"

-

-// convert an mpint into a big endian byte array (most significant byte first; left adjusted)

-//   return number of bytes converted

-//   if p == nil, allocate and result array

-int

-mptobe(mpint *b, uchar *p, uint n, uchar **pp)

-{

-	uint m;

-

-	m = (mpsignif(b)+7)/8;

-	if(m == 0)

-		m++;

-	if(p == nil){

-		n = m;

-		p = malloc(n);

-		if(p == nil)

-			sysfatal("mptobe: %r");

-	} else {

-		if(n < m)

-			return -1;

-		if(n > m)

-			memset(p+m, 0, n-m);

-	}

-	if(pp != nil)

-		*pp = p;

-	mptober(b, p, m);

-	return m;

-}

--- a/llt/mptober.c

+++ /dev/null

@@ -1,32 +1,0 @@

-#include "platform.h"

-

-void

-mptober(mpint *b, uchar *p, int n)

-{

-	int i, j, m;

-	mpdigit x;

-

-	memset(p, 0, n);

-

-	p += n;

-	m = b->top*Dbytes;

-	if(m < n)

-		n = m;

-

-	i = 0;

-	while(n >= Dbytes){

-		n -= Dbytes;

-		x = b->p[i++];

-		for(j = 0; j < Dbytes; j++){

-			*--p = x;

-			x >>= 8;

-		}

-	}

-	if(n > 0){

-		x = b->p[i];

-		for(j = 0; j < n; j++){

-			*--p = x;

-			x >>= 8;

-		}

-	}

-}

--- a/llt/mptod.c

+++ /dev/null

@@ -1,83 +1,0 @@

-#include "platform.h"

-

-extern double D_PINF, D_NINF;

-

-double

-mptod(mpint *a)

-{

-	u64int v;

-	mpdigit w, r;

-	int sf, i, n, m, s;

-	FPdbleword x;

-	

-	if(a->top == 0) return 0.0;

-	sf = mpsignif(a);

-	if(sf > 1024) return a->sign < 0 ? D_NINF : D_PINF;

-	i = a->top - 1;

-	v = a->p[i];

-	n = sf & Dbits - 1;

-	n |= n - 1 & Dbits;

-	r = 0;

-	if(n > 54){

-		s = n - 54;

-		r = v & (1<<s) - 1;

-		v >>= s;

-	}

-	while(n < 54){

-		if(--i < 0)

-			w = 0;

-		else

-			w = a->p[i];

-		m = 54 - n;

-		if(m > Dbits) m = Dbits;

-		s = Dbits - m & Dbits - 1;

-		v = v << m | w >> s;

-		r = w & (1<<s) - 1;

-		n += m;

-	}

-	if((v & 3) == 1){

-		while(--i >= 0)

-			r |= a->p[i];

-		if(r != 0)

-			v++;

-	}else

-		v++;

-	v >>= 1;

-	while((v >> 53) != 0){

-		v >>= 1;

-		if(++sf > 1024)

-			return a->sign < 0 ? D_NINF : D_PINF;

-	}

-	x.lo = v;

-	x.hi = (u32int)(v >> 32) & (1<<20) - 1 | (sf + 1022) << 20 | a->sign & 1<<31;

-	return x.x;

-}

-

-mpint *

-dtomp(double d, mpint *a)

-{

-	FPdbleword x;

-	uvlong v;

-	int e;

-

-	if(a == nil)

-		a = mpnew(0);

-	x.x = d;

-	e = x.hi >> 20 & 2047;

-	assert(e != 2047);

-	if(e < 1022){

-		mpassign(mpzero, a);

-		return a;

-	}

-	v = x.lo | (uvlong)(x.hi & (1<<20) - 1) << 32 | 1ULL<<52;

-	if(e < 1075){

-		v += (1ULL<<(1074 - e)) - (~v >> (1075 - e) & 1);

-		v >>= 1075 - e;

-	}

-	uvtomp(v, a);

-	if(e > 1075)

-		mpleft(a, e - 1075, a);

-	if((int)x.hi < 0)

-		a->sign = -1;

-	return a;

-}

--- a/llt/mptoi.c

+++ /dev/null

@@ -1,41 +1,0 @@

-#include "platform.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);

-	}

-	b->sign = (i >> (sizeof(i)*8 - 1)) | 1;

-	i *= b->sign;

-	*b->p = i;

-	b->top = 1;

-	return mpnorm(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;

-}

--- a/llt/mptoui.c

+++ /dev/null

@@ -1,31 +1,0 @@

-#include "platform.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);

-	}

-	*b->p = i;

-	b->top = 1;

-	b->sign = 1;

-	return mpnorm(b);

-}

-

-uint

-mptoui(mpint *b)

-{

-	uint x;

-

-	x = *b->p;

-	if(b->sign < 0)

-		x = 0;

-	else if(b->top > 1 || (sizeof(mpdigit) > sizeof(uint) && x > MAXUINT))

-		x =  MAXUINT;

-	return x;

-}

--- a/llt/mptouv.c

+++ /dev/null

@@ -1,44 +1,0 @@

-#include "platform.h"

-

-#define VLDIGITS (int)(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*Dbits);

-	}else

-		mpbits(b, VLDIGITS*Dbits);

-	b->sign = 1;

-	for(s = 0; s < VLDIGITS; s++){

-		b->p[s] = v;

-		v >>= sizeof(mpdigit)*8;

-	}

-	b->top = s;

-	return mpnorm(b);

-}

-

-uvlong

-mptouv(mpint *b)

-{

-	uvlong v;

-	int s;

-

-	if(b->top == 0 || b->sign < 0)

-		return 0LL;

-

-	if(b->top > VLDIGITS)

-		return -1LL;

-

-	v = 0ULL;

-	for(s = 0; s < b->top; s++)

-		v |= (uvlong)b->p[s]<<(s*sizeof(mpdigit)*8);

-

-	return v;

-}

--- a/llt/mptov.c

+++ /dev/null

@@ -1,60 +1,0 @@

-#include "platform.h"

-

-#define VLDIGITS (int)(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*Dbits);

-	}else

-		mpbits(b, VLDIGITS*Dbits);

-	b->sign = (v >> (sizeof(v)*8 - 1)) | 1;

-	uv = v * b->sign;

-	for(s = 0; s < VLDIGITS; s++){

-		b->p[s] = uv;

-		uv >>= sizeof(mpdigit)*8;

-	}

-	b->top = s;

-	return mpnorm(b);

-}

-

-vlong

-mptov(mpint *b)

-{

-	uvlong v;

-	int s;

-

-	if(b->top == 0)

-		return 0LL;

-

-	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 |= (uvlong)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;

-}

--- a/llt/mpvecadd.c

+++ /dev/null

@@ -1,34 +1,0 @@

-#include "platform.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;

-	uint 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;

-}

--- a/llt/mpveccmp.c

+++ /dev/null

@@ -1,25 +1,0 @@

-#include "platform.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;

-}

--- a/llt/mpvecdigmuladd.c

+++ /dev/null

@@ -1,101 +1,0 @@

-#include "platform.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;

-}

--- a/llt/mpvecsub.c

+++ /dev/null

@@ -1,32 +1,0 @@

-#include "platform.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 < (mpdigit)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;

-	}

-}

--- a/llt/mpvectscmp.c

+++ /dev/null

@@ -1,32 +1,0 @@

-#include "platform.h"

-

-int

-mpvectscmp(mpdigit *a, int alen, mpdigit *b, int blen)

-{

-	mpdigit x, y, z, v;

-	int m, p;

-

-	if(alen > blen){

-		v = 0;

-		while(alen > blen)

-			v |= a[--alen];

-		m = p = (-v^v|v)>>(Dbits-1);

-	} else if(blen > alen){

-		v = 0;

-		while(blen > alen)

-			v |= b[--blen];

-		m = (-v^v|v)>>(Dbits-1);

-		p = m^1;

-	} else

-		m = p = 0;

-	while(alen-- > 0){

-		x = a[alen];

-		y = b[alen];

-		z = x - y;

-		x = ~x;

-		v = ((-z^z|z)>>(Dbits-1)) & ~m;

-		p = ((~(x&y|x&z|y&z)>>(Dbits-1)) & v) | (p & ~v);

-		m |= v;

-	}

-	return (p-m) | m;

-}

--- a/llt/strtomp.c

+++ /dev/null

@@ -1,174 +1,0 @@

-#include "platform.h"

-

-static char*

-frompow2(char *a, mpint *b, int s)

-{

-	char *p, *next;

-	mpdigit x;

-	int i;

-

-	i = 1<<s;

-	for(p = a; (dec16chr(*p) & 255) < i; p++)

-		;

-

-	mpbits(b, (p-a)*s);

-	b->top = 0;

-	next = p;

-

-	while(p > a){

-		x = 0;

-		for(i = 0; i < Dbits; i += s){

-			if(p <= a)

-				break;

-			x |= dec16chr(*--p)<<i;

-		}

-		b->p[b->top++] = x;

-	}

-	return next;

-}

-

-static char*

-from8(char *a, mpint *b)

-{

-	char *p, *next;

-	mpdigit x, y;

-	int i;

-

-	for(p = a; ((*p - '0') & 255) < 8; p++)

-		;

-

-	mpbits(b, (p-a)*3);

-	b->top = 0;

-	next = p;

-

-	i = 0;

-	x = y = 0;

-	while(p > a){

-		y = *--p - '0';

-		x |= y << i;

-		i += 3;

-		if(i >= Dbits){

-Digout:

-			i -= Dbits;

-			b->p[b->top++] = x;

-			x = y >> (3-i);

-		}

-	}

-	if(i > 0)

-		goto Digout;

-

-	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 = *a - '0';

-			if(y > 9)

-				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;

-}

-

-mpint*

-strtomp(char *a, char **pp, int base, mpint *b)

-{

-	int sign;

-	char *e;

-

-	if(b == nil){

-		b = mpnew(0);

-	}

-

-	while(*a==' ' || *a=='\t')

-		a++;

-

-	sign = 1;

-	for(;; a++){

-		switch(*a){

-		case '-':

-			sign *= -1;

-			continue;

-		}

-		break;

-	}

-

-	if(base == 0){

-		base = 10;

-		if(a[0] == '0'){

-			if(a[1] == 'x' || a[1] == 'X') {

-				a += 2;

-				base = 16;

-			} else if(a[1] == 'b' || a[1] == 'B') {

-				a += 2;

-				base = 2;

-			} else if(a[1] >= '0' && a[1] <= '7') {

-				a++;

-				base = 8;

-			}

-		}

-	}

-

-	switch(base){

-	case 2:

-		e = frompow2(a, b, 1);

-		break;

-	case 4:

-		e = frompow2(a, b, 2);

-		break;

-	case 8:

-		e = from8(a, b);

-		break;

-	case 10:

-		e = from10(a, b);

-		break;

-	case 16:

-		e = frompow2(a, b, 4);

-		break;

-	default:

-		abort();

-		return nil;

-	}

-

-	if(pp != nil)

-		*pp = e;

-

-	// if no characters parsed, there wasn't a number to convert

-	if(e == a)

-		return nil;

-

-	b->sign = sign;

-	return mpnorm(b);

-}

--- a/llt/u16.c

+++ /dev/null

@@ -1,68 +1,0 @@

-#include "platform.h"

-

-#define between(x,min,max)	(((min-1-x) & (x-max-1))>>8)

-

-int

-enc16chr(int o)

-{

-	int c;

-

-	c  = between(o,  0,  9) & ('0'+o);

-	c |= between(o, 10, 15) & ('A'+(o-10));

-	return c;

-}

-

-int

-dec16chr(int c)

-{

-	int o;

-

-	o  = between(c, '0', '9') & (1+(c-'0'));

-	o |= between(c, 'A', 'F') & (1+10+(c-'A'));

-	o |= between(c, 'a', 'f') & (1+10+(c-'a'));

-	return o-1;

-}

-

-int

-dec16(uchar *out, int lim, char *in, int n)

-{

-	int c, w = 0, i = 0;

-	uchar *start = out;

-	uchar *eout = out + lim;

-

-	while(n-- > 0){

-		c = dec16chr(*in++);

-		if(c < 0)

-			continue;

-		w = (w<<4) + c;

-		i++;

-		if(i == 2){

-			if(out + 1 > eout)

-				goto exhausted;

-			*out++ = w;

-			w = 0;

-			i = 0;

-		}

-	}

-exhausted:

-	return out - start;

-}

-

-int

-enc16(char *out, int lim, uchar *in, int n)

-{

-	uint c;

-	char *eout = out + lim;

-	char *start = out;

-

-	while(n-- > 0){

-		c = *in++;

-		if(out + 2 >= eout)

-			goto exhausted;

-		*out++ = enc16chr(c>>4);

-		*out++ = enc16chr(c&15);

-	}

-exhausted:

-	*out = 0;

-	return out - start;

-}

--- a/llt/u32.c

+++ /dev/null

@@ -1,143 +1,0 @@

-#include "platform.h"

-

-#define between(x,min,max)	(((min-1-x) & (x-max-1))>>8)

-

-int

-enc32chr(int o)

-{

-	int c;

-

-	c  = between(o,  0, 25) & ('A'+o);

-	c |= between(o, 26, 31) & ('2'+(o-26));

-	return c;

-}

-

-int

-dec32chr(int c)

-{

-	int o;

-

-	o  = between(c, 'A', 'Z') & (1+(c-'A'));

-	o |= between(c, 'a', 'z') & (1+(c-'a'));

-	o |= between(c, '2', '7') & (1+26+(c-'2'));

-	return o-1;

-}

-

-int

-dec32x(uchar *dest, int ndest, char *src, int nsrc, int (*chr)(int))

-{

-	uchar *start;

-	int i, j, u[8];

-

-	if(ndest+1 < (5*nsrc+7)/8)

-		return -1;

-	start = dest;

-	while(nsrc>=8){

-		for(i=0; i<8; i++){

-			j = chr(src[i]);

-			if(j < 0)

-				j = 0;

-			u[i] = j;

-		}

-		*dest++ = (u[0]<<3) | (0x7 & (u[1]>>2));

-		*dest++ = ((0x3 & u[1])<<6) | (u[2]<<1) | (0x1 & (u[3]>>4));

-		*dest++ = ((0xf & u[3])<<4) | (0xf & (u[4]>>1));

-		*dest++ = ((0x1 & u[4])<<7) | (u[5]<<2) | (0x3 & (u[6]>>3));

-		*dest++ = ((0x7 & u[6])<<5) | u[7];

-		src  += 8;

-		nsrc -= 8;

-	}

-	if(nsrc > 0){

-		if(nsrc == 1 || nsrc == 3 || nsrc == 6)

-			return -1;

-		for(i=0; i<nsrc; i++){

-			j = chr(src[i]);

-			if(j < 0)

-				j = 0;

-			u[i] = j;

-		}

-		*dest++ = (u[0]<<3) | (0x7 & (u[1]>>2));

-		if(nsrc == 2)

-			goto out;

-		*dest++ = ((0x3 & u[1])<<6) | (u[2]<<1) | (0x1 & (u[3]>>4));

-		if(nsrc == 4)

-			goto out;

-		*dest++ = ((0xf & u[3])<<4) | (0xf & (u[4]>>1));

-		if(nsrc == 5)

-			goto out;

-		*dest++ = ((0x1 & u[4])<<7) | (u[5]<<2) | (0x3 & (u[6]>>3));

-	}

-out:

-	return dest-start;

-}

-

-int

-enc32x(char *dest, int ndest, uchar *src, int nsrc, int (*chr)(int))

-{

-	char *start;

-	int j;

-

-	if(ndest <= (8*nsrc+4)/5)

-		return -1;

-	start = dest;

-	while(nsrc>=5){

-		j = (0x1f & (src[0]>>3));

-		*dest++ = chr(j);

-		j = (0x1c & (src[0]<<2)) | (0x03 & (src[1]>>6));

-		*dest++ = chr(j);

-		j = (0x1f & (src[1]>>1));

-		*dest++ = chr(j);

-		j = (0x10 & (src[1]<<4)) | (0x0f & (src[2]>>4));

-		*dest++ = chr(j);

-		j = (0x1e & (src[2]<<1)) | (0x01 & (src[3]>>7));

-		*dest++ = chr(j);

-		j = (0x1f & (src[3]>>2));

-		*dest++ = chr(j);

-		j = (0x18 & (src[3]<<3)) | (0x07 & (src[4]>>5));

-		*dest++ = chr(j);

-		j = (0x1f & (src[4]));

-		*dest++ = chr(j);

-		src  += 5;

-		nsrc -= 5;

-	}

-	if(nsrc){

-		j = (0x1f & (src[0]>>3));

-		*dest++ = chr(j);

-		j = (0x1c & (src[0]<<2));

-		if(nsrc == 1)

-			goto out;

-		j |= (0x03 & (src[1]>>6));

-		*dest++ = chr(j);

-		j = (0x1f & (src[1]>>1));

-		*dest++ = chr(j);

-		j = (0x10 & (src[1]<<4));

-		if(nsrc == 2)

-			goto out;

-		j |= (0x0f & (src[2]>>4));

-		*dest++ = chr(j);

-		j = (0x1e & (src[2]<<1));

-		if(nsrc == 3)

-			goto out;

-		j |= (0x01 & (src[3]>>7));

-		*dest++ = chr(j);

-		j = (0x1f & (src[3]>>2));

-		*dest++ = chr(j);

-		j = (0x18 & (src[3]<<3));

-out:

-		*dest++ = chr(j);

-	}

-	*dest = 0;

-	return dest-start;

-}

-

-int

-enc32(char *dest, int ndest, uchar *src, int nsrc)

-{

-	return enc32x(dest, ndest, src, nsrc, enc32chr);

-}

-

-int

-dec32(uchar *dest, int ndest, char *src, int nsrc)

-{

-	return dec32x(dest, ndest, src, nsrc, dec32chr);

-}

--- a/llt/u64.c

+++ /dev/null

@@ -1,141 +1,0 @@

-#include "platform.h"

-

-#define between(x,min,max)	(((min-1-x) & (x-max-1))>>8)

-

-int

-enc64chr(int o)

-{

-	int c;

-

-	c  = between(o,  0, 25) & ('A'+o);

-	c |= between(o, 26, 51) & ('a'+(o-26));

-	c |= between(o, 52, 61) & ('0'+(o-52));

-	c |= between(o, 62, 62) & ('+');

-	c |= between(o, 63, 63) & ('/');

-	return c;

-}

-

-int

-dec64chr(int c)

-{

-	int o;

-

-	o  = between(c, 'A', 'Z') & (1+(c-'A'));

-	o |= between(c, 'a', 'z') & (1+26+(c-'a'));

-	o |= between(c, '0', '9') & (1+52+(c-'0'));

-	o |= between(c, '+', '+') & (1+62);

-	o |= between(c, '/', '/') & (1+63);

-	return o-1;

-}

-

-int

-dec64x(uchar *out, int lim, char *in, int n, int (*chr)(int))

-{

-	ulong b24;

-	uchar *start = out;

-	uchar *e = out + lim;

-	int i, c;

-

-	b24 = 0;

-	i = 0;

-	while(n-- > 0){

-		c = chr(*in++);

-		if(c < 0)

-			continue;

-		switch(i){

-		case 0:

-			b24 = c<<18;

-			break;

-		case 1:

-			b24 |= c<<12;

-			break;

-		case 2:

-			b24 |= c<<6;

-			break;

-		case 3:

-			if(out + 3 > e)

-				goto exhausted;

-

-			b24 |= c;

-			*out++ = b24>>16;

-			*out++ = b24>>8;

-			*out++ = b24;

-			i = 0;

-			continue;

-		}

-		i++;

-	}

-	switch(i){

-	case 2:

-		if(out + 1 > e)

-			goto exhausted;

-		*out++ = b24>>16;

-		break;

-	case 3:

-		if(out + 2 > e)

-			goto exhausted;

-		*out++ = b24>>16;

-		*out++ = b24>>8;

-		break;

-	}

-exhausted:

-	return out - start;

-}

-

-int

-enc64x(char *out, int lim, uchar *in, int n, int (*chr)(int))

-{

-	int i;

-	ulong b24;

-	char *start = out;

-	char *e = out + lim;

-

-	for(i = n/3; i > 0; i--){

-		b24 = *in++<<16;

-		b24 |= *in++<<8;

-		b24 |= *in++;

-		if(out + 4 >= e)

-			goto exhausted;

-		*out++ = chr(b24>>18);

-		*out++ = chr((b24>>12)&0x3f);

-		*out++ = chr((b24>>6)&0x3f);

-		*out++ = chr(b24&0x3f);

-	}

-

-	switch(n%3){

-	case 2:

-		b24 = *in++<<16;

-		b24 |= *in<<8;

-		if(out + 4 >= e)

-			goto exhausted;

-		*out++ = chr(b24>>18);

-		*out++ = chr((b24>>12)&0x3f);

-		*out++ = chr((b24>>6)&0x3f);

-		*out++ = '=';

-		break;

-	case 1:

-		b24 = *in<<16;

-		if(out + 4 >= e)

-			goto exhausted;

-		*out++ = chr(b24>>18);

-		*out++ = chr((b24>>12)&0x3f);

-		*out++ = '=';

-		*out++ = '=';

-		break;

-	}

-exhausted:

-	*out = 0;

-	return out - start;

-}

-

-int

-enc64(char *out, int lim, uchar *in, int n)

-{

-	return enc64x(out, lim, in, n, enc64chr);

-}

-

-int

-dec64(uchar *out, int lim, char *in, int n)

-{

-	return dec64x(out, lim, in, n, dec64chr);

-}

--- a/mkfile

+++ b/mkfile

@@ -25,16 +25,24 @@

 	iostream.$O\

 	string.$O\

 	table.$O\

+	llt/bitvector-ops.$O\

+	llt/bitvector.$O\

+	llt/dump.$O\

+	llt/hashing.$O\

+	llt/htable.$O\

+	llt/int2str.$O\

+	llt/ios.$O\

+	llt/lltinit.$O\

+	llt/ptrhash.$O\

+	llt/random.$O\

+	llt/timefuncs.$O\

+	llt/utf8.$O\

+	llt/wcwidth.$O\

 default:V: all

 </sys/src/cmd/mkone

-$O.out: llt/libllt.$O.a

-

-llt/libllt.$O.a:

-	cd llt && mk

-

 boot.h: flisp.boot

 	sed 's,\\,\\\\,g;s,",\\",g;s,^,",g;s,$,\\n",g' $prereq >$target

@@ -45,6 +53,9 @@

 flisp.$O: maxstack.inc opcodes.h builtin_fns.h $CHFILES

+%.$O: %.c

+	$CC $CFLAGS -o $target $stem.c

+

 bootstrap:V: $O.out

     ./$O.out gen.lsp && \

 	cp flisp.boot flisp.boot.bak && \

@@ -54,8 +65,6 @@

 nuke:V:

 	rm -f *.[$OS] [$OS].out *.acid $TARG $CLEANFILES

-	cd llt && mk nuke

 clean:V:

 	rm -f *.[$OS] [$OS].out $TARG $CLEANFILES

-	cd llt && mk clean

--- /dev/null

+++ b/mp/mpadd.c

@@ -1,0 +1,56 @@

+#include "platform.h"

+

+// sum = abs(b1) + abs(b2), i.e., add the magnitudes

+void

+mpmagadd(mpint *b1, mpint *b2, mpint *sum)

+{

+	int m, n;

+	mpint *t;

+

+	sum->flags |= (b1->flags | b2->flags) & MPtimesafe;

+

+	// 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);

+		sum->sign = 1;

+		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){

+		assert(((b1->flags | b2->flags | sum->flags) & MPtimesafe) == 0);

+		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;

+	}

+}

--- /dev/null

+++ b/mp/mpaux.c

@@ -1,0 +1,201 @@

+#include "platform.h"

+

+static mpdigit _mptwodata[1] = { 2 };

+static mpint _mptwo =

+{

+	1, 1, 1,

+	_mptwodata,

+	MPstatic|MPnorm

+};

+mpint *mptwo = &_mptwo;

+

+static mpdigit _mponedata[1] = { 1 };

+static mpint _mpone =

+{

+	1, 1, 1,

+	_mponedata,

+	MPstatic|MPnorm

+};

+mpint *mpone = &_mpone;

+

+static mpdigit _mpzerodata[1] = { 0 };

+static mpint _mpzero =

+{

+	1, 1, 0,

+	_mpzerodata,

+	MPstatic|MPnorm

+};

+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");

+

+	n = DIGITS(n);

+	if(n < mpmindigits)

+		n = mpmindigits;

+	b = calloc(1, sizeof(mpint) + n*Dbytes);

+	if(b == nil)

+		sysfatal("mpnew: %r");

+	b->p = (mpdigit*)&b[1];

+	b->size = n;

+	b->sign = 1;

+	b->flags = MPnorm;

+

+	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;

+	} else {

+		if(b->p == (mpdigit*)&b[1]){

+			b->p = (mpdigit*)malloc(n*Dbytes);

+			if(b->p == nil)

+				sysfatal("mpbits: %r");

+			memmove(b->p, &b[1], Dbytes*b->top);

+			memset(&b[1], 0, Dbytes*b->size);

+		} else {

+			b->p = (mpdigit*)realloc(b->p, n*Dbytes);

+			if(b->p == nil)

+				sysfatal("mpbits: %r");

+		}

+		b->size = n;

+	}

+	memset(&b->p[b->top], 0, Dbytes*(n - b->top));

+	b->top = n;

+	b->flags &= ~MPnorm;

+}

+

+void

+mpfree(mpint *b)

+{

+	if(b == nil)

+		return;

+	if(b->flags & MPstatic)

+		sysfatal("freeing mp constant");

+	memset(b->p, 0, b->size*Dbytes);

+	if(b->p != (mpdigit*)&b[1])

+		free(b->p);

+	free(b);

+}

+

+mpint*

+mpnorm(mpint *b)

+{

+	int i;

+

+	if(b->flags & MPtimesafe){

+		assert(b->sign == 1);

+		b->flags &= ~MPnorm;

+		return b;

+	}

+	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;

+	b->flags |= MPnorm;

+	return b;

+}

+

+mpint*

+mpcopy(mpint *old)

+{

+	mpint *new;

+

+	new = mpnew(Dbits*old->size);

+	new->sign = old->sign;

+	new->top = old->top;

+	new->flags = old->flags & ~(MPstatic|MPfield);

+	memmove(new->p, old->p, Dbytes*old->top);

+	return new;

+}

+

+void

+mpassign(mpint *old, mpint *new)

+{

+	if(new == nil || old == new)

+		return;

+	new->top = 0;

+	mpbits(new, Dbits*old->top);

+	new->sign = old->sign;

+	new->top = old->top;

+	new->flags &= ~MPnorm;

+	new->flags |= old->flags & ~(MPstatic|MPfield);

+	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;

+

+	assert(n->flags & MPnorm);

+	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;

+}

--- /dev/null

+++ b/mp/mpcmp.c

@@ -1,0 +1,28 @@

+#include "platform.h"

+

+// return neg, 0, pos as abs(b1)-abs(b2) is neg, 0, pos

+int

+mpmagcmp(mpint *b1, mpint *b2)

+{

+	int i;

+

+	i = b1->flags | b2->flags;

+	if(i & MPtimesafe)

+		return mpvectscmp(b1->p, b1->top, b2->p, b2->top);

+	if(i & MPnorm){

+		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)

+{

+	int sign;

+

+	sign = (b1->sign - b2->sign) >> 1;	// -1, 0, 1

+	return sign | (sign&1)-1 & mpmagcmp(b1, b2)*b1->sign;

+}

--- /dev/null

+++ b/mp/mpdigdiv.c

@@ -1,0 +1,54 @@

+#include "platform.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;

+	}

+

+	// very common case

+	if(~divisor == 0){

+		lo += hi;

+		if(lo < hi){

+			hi++;

+			lo++;

+		}

+		if(lo+1 == 0)

+			hi++;

+		*quotient = hi;

+		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 |= 1U<<i;

+	}

+	q += lo/divisor;

+	*quotient = q;

+}

--- /dev/null

+++ b/mp/mpdiv.c

@@ -1,0 +1,140 @@

+#include "platform.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, qsign, rsign;

+	mpdigit qd, *up, *vp, *qp;

+	mpint *u, *v, *t;

+

+	assert(quotient != remainder);

+	assert(divisor->flags & MPnorm);

+

+	// divide bv zero

+	if(divisor->top == 0)

+		abort();

+

+	// division by one or small powers of two

+	if(divisor->top == 1 && (divisor->p[0] & divisor->p[0]-1) == 0){

+		vlong r = 0;

+		if(dividend->top > 0)

+			r = (vlong)dividend->sign * (dividend->p[0] & divisor->p[0]-1);

+		if(quotient != nil){

+			sign = divisor->sign;

+			for(s = 0; ((divisor->p[0] >> s) & 1) == 0; s++)

+				;

+			mpright(dividend, s, quotient);

+			if(sign < 0)

+				quotient->sign ^= (-mpmagcmp(quotient, mpzero) >> 31) << 1;

+		}

+		if(remainder != nil){

+			remainder->flags |= dividend->flags & MPtimesafe;

+			vtomp(r, remainder);

+		}

+		return;

+	}

+	assert((dividend->flags & MPtimesafe) == 0);

+

+	// quick check

+	if(mpmagcmp(dividend, divisor) < 0){

+		if(remainder != nil)

+			mpassign(dividend, remainder);

+		if(quotient != nil)

+			mpassign(mpzero, quotient);

+		return;

+	}

+	

+	qsign = divisor->sign * dividend->sign;

+	rsign = dividend->sign;

+

+	// 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){

+		assert((quotient->flags & MPtimesafe) == 0);

+		mpnorm(quotient);

+		if(quotient->top != 0)

+			quotient->sign = qsign;

+	}

+

+	if(remainder != nil){

+		assert((remainder->flags & MPtimesafe) == 0);

+		mpright(u, s, remainder);	// u is the remainder shifted

+		if(remainder->top != 0)

+			remainder->sign = rsign;

+	}

+

+	mpfree(t);

+	mpfree(u);

+	if(v != divisor)

+		mpfree(v);

+}

--- /dev/null

+++ b/mp/mpfmt.c

@@ -1,0 +1,207 @@

+#include "platform.h"

+

+static int

+toencx(mpint *b, char *buf, int len, int (*enc)(char*, int, uchar*, int))

+{

+	uchar *p;

+	int n, rv;

+

+	p = nil;

+	n = mptobe(b, nil, 0, &p);

+	if(n < 0)

+		return -1;

+	rv = (*enc)(buf, len, p, n);

+	free(p);

+	return rv;

+}

+

+static int

+topow2(mpint *b, char *buf, int len, int s)

+{

+	mpdigit *p, x;

+	int i, j, sn;

+	char *out, *eout;

+

+	if(len < 1)

+		return -1;

+

+	sn = 1<<s;

+	out = buf;

+	eout = buf+len;

+	for(p = &b->p[b->top-1]; p >= b->p; p--){

+		x = *p;

+		for(i = Dbits-s; i >= 0; i -= s){

+			j = x >> i & sn - 1;

+			if(j != 0 || out != buf){

+				if(out >= eout)

+					return -1;

+				*out++ = enc16chr(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);

+	d->flags &= ~MPtimesafe;

+	mpnorm(d);

+	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;

+}

+

+static int

+to8(mpint *b, char *buf, int len)

+{

+	mpdigit x, y;

+	char *out;

+	int i, j;

+

+	if(len < 2)

+		return -1;

+

+	out = buf+len;

+	*--out = 0;

+

+	i = j = 0;

+	x = y = 0;

+	while(j < b->top){

+		y = b->p[j++];

+		if(i > 0)

+			x |= y << i;

+		else

+			x = y;

+		i += Dbits;

+		while(i >= 3){

+Digout:			i -= 3;

+			if(out > buf)

+				out--;

+			else if(x != 0)

+				return -1;

+			*out = '0' + (x & 7);

+			x = y >> (Dbits-i);

+		}

+	}

+	if(i > 0)

+		goto Digout;

+

+	while(*out == '0') out++;

+	if(*out == '\0')

+		*--out = '0';

+

+	len -= out-buf;

+	if(out != buf)

+		memmove(buf, out, len);

+	return 0;

+}

+

+char*

+mptoa(mpint *b, int base, char *buf, int len)

+{

+	char *out;

+	int rv, alloced;

+

+	if(base == 0)

+		base = 16;	/* default */

+	alloced = 0;

+	if(buf == nil){

+		/* rv <= log₂(base) */

+		for(rv=1; (base >> rv) > 1; rv++)

+			;

+		len = 10 + (b->top*Dbits / rv);

+		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 = toencx(b, out, len, enc64);

+		break;

+	case 32:

+		rv = toencx(b, out, len, enc32);

+		break;

+	case 16:

+		rv = topow2(b, out, len, 4);

+		break;

+	case 10:

+		rv = to10(b, out, len);

+		break;

+	case 8:

+		rv = to8(b, out, len);

+		break;

+	case 4:

+		rv = topow2(b, out, len, 2);

+		break;

+	case 2:

+		rv = topow2(b, out, len, 1);

+		break;

+	default:

+		abort();

+		return nil;

+	}

+	if(rv < 0){

+		if(alloced)

+			free(buf);

+		return nil;

+	}

+	return buf;

+}

--- /dev/null

+++ b/mp/mpleft.c

@@ -1,0 +1,49 @@

+#include "platform.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 zero or 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;

+

+	res->flags |= b->flags & MPtimesafe;

+	mpnorm(res);

+}

--- /dev/null

+++ b/mp/mplogic.c

@@ -1,0 +1,210 @@

+#include "platform.h"

+

+/*

+	mplogic calculates b1|b2 subject to the

+	following flag bits (fl)

+

+	bit 0: subtract 1 from b1

+	bit 1: invert b1

+	bit 2: subtract 1 from b2

+	bit 3: invert b2

+	bit 4: add 1 to output

+	bit 5: invert output

+	

+	it inverts appropriate bits automatically

+	depending on the signs of the inputs

+*/

+

+static void

+mplogic(mpint *b1, mpint *b2, mpint *sum, int fl)

+{

+	mpint *t;

+	mpdigit *dp1, *dp2, *dpo, d1, d2, d;

+	int c1, c2, co;

+	int i;

+

+	assert(((b1->flags | b2->flags | sum->flags) & MPtimesafe) == 0);

+	if(b1->sign < 0) fl ^= 0x03;

+	if(b2->sign < 0) fl ^= 0x0c;

+	sum->sign = (int)(((fl|fl>>2)^fl>>4)<<30)>>31|1;

+	if(sum->sign < 0) fl ^= 0x30;

+	if(b2->top > b1->top){

+		t = b1;

+		b1 = b2;

+		b2 = t;

+		fl = fl >> 2 & 0x03 | fl << 2 & 0x0c | fl & 0x30;

+	}

+	mpbits(sum, b1->top*Dbits+1);

+	dp1 = b1->p;

+	dp2 = b2->p;

+	dpo = sum->p;

+	c1 = fl & 1;

+	c2 = fl >> 2 & 1;

+	co = fl >> 4 & 1;

+	for(i = 0; i < b1->top; i++){

+		d1 = dp1[i] - c1;

+		if(i < b2->top)

+			d2 = dp2[i] - c2;

+		else

+			d2 = 0;

+		if(d1 != (mpdigit)-1) c1 = 0;

+		if(d2 != (mpdigit)-1) c2 = 0;

+		if((fl & 2) != 0) d1 ^= -1;

+		if((fl & 8) != 0) d2 ^= -1;

+		d = d1 | d2;

+		if((fl & 32) != 0) d ^= -1;

+		d += co;

+		if(d != 0) co = 0;

+		dpo[i] = d;

+	}

+	sum->top = i;

+	if(co)

+		dpo[sum->top++] = co;

+	mpnorm(sum);

+}

+

+void

+mpor(mpint *b1, mpint *b2, mpint *sum)

+{

+	mplogic(b1, b2, sum, 0);

+}

+

+void

+mpand(mpint *b1, mpint *b2, mpint *sum)

+{

+	mplogic(b1, b2, sum, 0x2a);

+}

+

+void

+mpbic(mpint *b1, mpint *b2, mpint *sum)

+{

+	mplogic(b1, b2, sum, 0x22);

+}

+

+void

+mpnot(mpint *b, mpint *r)

+{

+	mpadd(b, mpone, r);

+	if(r->top != 0)

+		r->sign ^= -2;

+}

+

+void

+mpxor(mpint *b1, mpint *b2, mpint *sum)

+{

+	mpint *t;

+	mpdigit *dp1, *dp2, *dpo, d1, d2, d;

+	int c1, c2, co;

+	int i, fl;

+

+	assert(((b1->flags | b2->flags | sum->flags) & MPtimesafe) == 0);

+	if(b2->top > b1->top){

+		t = b1;

+		b1 = b2;

+		b2 = t;

+	}

+	fl = (b1->sign & 10) ^ (b2->sign & 12);

+	sum->sign = (int)(fl << 28) >> 31 | 1;

+	mpbits(sum, b1->top*Dbits+1);

+	dp1 = b1->p;

+	dp2 = b2->p;

+	dpo = sum->p;

+	c1 = fl >> 1 & 1;

+	c2 = fl >> 2 & 1;

+	co = fl >> 3 & 1;

+	for(i = 0; i < b1->top; i++){

+		d1 = dp1[i] - c1;

+		if(i < b2->top)

+			d2 = dp2[i] - c2;

+		else

+			d2 = 0;

+		if(d1 != (mpdigit)-1) c1 = 0;

+		if(d2 != (mpdigit)-1) c2 = 0;

+		d = d1 ^ d2;

+		d += co;

+		if(d != 0) co = 0;

+		dpo[i] = d;

+	}

+	sum->top = i;

+	if(co)

+		dpo[sum->top++] = co;

+	mpnorm(sum);

+}

+

+void

+mptrunc(mpint *b, int n, mpint *r)

+{

+	int d, m, i, c;

+

+	assert(((b->flags | r->flags) & MPtimesafe) == 0);

+	mpbits(r, n);

+	r->top = DIGITS(n);

+	d = n / Dbits;

+	m = n % Dbits;

+	if(b->sign == -1){

+		c = 1;

+		for(i = 0; i < r->top; i++){

+			if(i < b->top)

+				r->p[i] = ~(b->p[i] - c);

+			else

+				r->p[i] = -1;

+			if(r->p[i] != 0)

+				c = 0;

+		}

+		if(m != 0)

+			r->p[d] &= (1<<m) - 1;

+	}else if(b->sign == 1){

+		if(d >= b->top){

+			mpassign(b, r);

+			mpnorm(r);

+			return;

+		}

+		if(b != r)

+			for(i = 0; i < d; i++)

+				r->p[i] = b->p[i];

+		if(m != 0)

+			r->p[d] = b->p[d] & (1<<m)-1;

+	}

+	r->sign = 1;

+	mpnorm(r);

+}

+

+void

+mpxtend(mpint *b, int n, mpint *r)

+{

+	int d, m, c, i;

+

+	d = (n - 1) / Dbits;

+	m = (n - 1) % Dbits;

+	if(d >= b->top){

+		mpassign(b, r);

+		return;

+	}

+	mptrunc(b, n, r);

+	mpbits(r, n);

+	if((r->p[d] & 1<<m) == 0){

+		mpnorm(r);

+		return;

+	}

+	r->p[d] |= -(1<<m);

+	r->sign = -1;

+	c = 1;

+	for(i = 0; i < r->top; i++){

+		r->p[i] = ~(r->p[i] - c);

+		if(r->p[i] != 0)

+			c = 0;

+	}

+	mpnorm(r);

+}

+

+void

+mpasr(mpint *b, int n, mpint *r)

+{

+	if(b->sign > 0 || n <= 0){

+		mpright(b, n, r);

+		return;

+	}

+	mpadd(b, mpone, r);

+	mpright(r, n, r);

+	mpsub(r, mpone, r);

+}

--- /dev/null

+++ b/mp/mpmul.c

@@ -1,0 +1,174 @@

+#include "platform.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 = calloc(1, Dbytes*5*(2*n+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(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

+mpvectsmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p)

+{

+	int i;

+	mpdigit *t;

+

+	if(alen < blen){

+		i = alen;

+		alen = blen;

+		blen = i;

+		t = a;

+		a = b;

+		b = t;

+	}

+	if(blen == 0)

+		return;

+	for(i = 0; i < blen; i++)

+		mpvecdigmuladd(a, alen, b[i], &p[i]);

+}

+

+void

+mpmul(mpint *b1, mpint *b2, mpint *prod)

+{

+	mpint *oprod;

+

+	oprod = prod;

+	if(prod == b1 || prod == b2){

+		prod = mpnew(0);

+		prod->flags = oprod->flags;

+	}

+	prod->flags |= (b1->flags | b2->flags) & MPtimesafe;

+

+	prod->top = 0;

+	mpbits(prod, (b1->top+b2->top+1)*Dbits);

+	if(prod->flags & MPtimesafe)

+		mpvectsmul(b1->p, b1->top, b2->p, b2->top, prod->p);

+	else

+		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 != prod){

+		mpassign(prod, oprod);

+		mpfree(prod);

+	}

+}

--- /dev/null

+++ b/mp/mpright.c

@@ -1,0 +1,55 @@

+#include "platform.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);

+	else if(shift == 0)

+		return;

+

+	d = shift/Dbits;

+	r = shift - d*Dbits;

+	l = Dbits - r;

+

+	//  shift all the bits out == zero

+	if(d>=b->top){

+		res->sign = 1;

+		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;

+	}

+

+	res->top = i;

+	res->flags |= b->flags & MPtimesafe;

+	mpnorm(res);

+}

--- /dev/null

+++ b/mp/mpsub.c

@@ -1,0 +1,54 @@

+#include "platform.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){

+		assert(((b1->flags | b2->flags | diff->flags) & MPtimesafe) == 0);

+		sign = -1;

+		t = b1;

+		b1 = b2;

+		b2 = t;

+	} else {

+		diff->flags |= (b1->flags | b2->flags) & MPtimesafe;

+		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){

+		assert(((b1->flags | b2->flags | diff->flags) & MPtimesafe) == 0);

+		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;

+}

--- /dev/null

+++ b/mp/mptobe.c

@@ -1,0 +1,29 @@

+#include "platform.h"

+

+// convert an mpint into a big endian byte array (most significant byte first; left adjusted)

+//   return number of bytes converted

+//   if p == nil, allocate and result array

+int

+mptobe(mpint *b, uchar *p, uint n, uchar **pp)

+{

+	uint m;

+

+	m = (mpsignif(b)+7)/8;

+	if(m == 0)

+		m++;

+	if(p == nil){

+		n = m;

+		p = malloc(n);

+		if(p == nil)

+			sysfatal("mptobe: %r");

+	} else {

+		if(n < m)

+			return -1;

+		if(n > m)

+			memset(p+m, 0, n-m);

+	}

+	if(pp != nil)

+		*pp = p;

+	mptober(b, p, m);

+	return m;

+}

--- /dev/null

+++ b/mp/mptober.c

@@ -1,0 +1,32 @@

+#include "platform.h"

+

+void

+mptober(mpint *b, uchar *p, int n)

+{

+	int i, j, m;

+	mpdigit x;

+

+	memset(p, 0, n);

+

+	p += n;

+	m = b->top*Dbytes;

+	if(m < n)

+		n = m;

+

+	i = 0;

+	while(n >= Dbytes){

+		n -= Dbytes;

+		x = b->p[i++];

+		for(j = 0; j < Dbytes; j++){

+			*--p = x;

+			x >>= 8;

+		}

+	}

+	if(n > 0){

+		x = b->p[i];

+		for(j = 0; j < n; j++){

+			*--p = x;

+			x >>= 8;

+		}

+	}

+}

--- /dev/null

+++ b/mp/mptod.c

@@ -1,0 +1,83 @@

+#include "platform.h"

+

+extern double D_PINF, D_NINF;

+

+double

+mptod(mpint *a)

+{

+	u64int v;

+	mpdigit w, r;

+	int sf, i, n, m, s;

+	FPdbleword x;

+	

+	if(a->top == 0) return 0.0;

+	sf = mpsignif(a);

+	if(sf > 1024) return a->sign < 0 ? D_NINF : D_PINF;

+	i = a->top - 1;

+	v = a->p[i];

+	n = sf & Dbits - 1;

+	n |= n - 1 & Dbits;

+	r = 0;

+	if(n > 54){

+		s = n - 54;

+		r = v & (1<<s) - 1;

+		v >>= s;

+	}

+	while(n < 54){

+		if(--i < 0)

+			w = 0;

+		else

+			w = a->p[i];

+		m = 54 - n;

+		if(m > Dbits) m = Dbits;

+		s = Dbits - m & Dbits - 1;

+		v = v << m | w >> s;

+		r = w & (1<<s) - 1;

+		n += m;

+	}

+	if((v & 3) == 1){

+		while(--i >= 0)

+			r |= a->p[i];

+		if(r != 0)

+			v++;

+	}else

+		v++;

+	v >>= 1;

+	while((v >> 53) != 0){

+		v >>= 1;

+		if(++sf > 1024)

+			return a->sign < 0 ? D_NINF : D_PINF;

+	}

+	x.lo = v;

+	x.hi = (u32int)(v >> 32) & (1<<20) - 1 | (sf + 1022) << 20 | a->sign & 1<<31;

+	return x.x;

+}

+

+mpint *

+dtomp(double d, mpint *a)

+{

+	FPdbleword x;

+	uvlong v;

+	int e;

+

+	if(a == nil)

+		a = mpnew(0);

+	x.x = d;

+	e = x.hi >> 20 & 2047;

+	assert(e != 2047);

+	if(e < 1022){

+		mpassign(mpzero, a);

+		return a;

+	}

+	v = x.lo | (uvlong)(x.hi & (1<<20) - 1) << 32 | 1ULL<<52;

+	if(e < 1075){

+		v += (1ULL<<(1074 - e)) - (~v >> (1075 - e) & 1);

+		v >>= 1075 - e;

+	}

+	uvtomp(v, a);

+	if(e > 1075)

+		mpleft(a, e - 1075, a);

+	if((int)x.hi < 0)

+		a->sign = -1;

+	return a;

+}

--- /dev/null

+++ b/mp/mptoi.c

@@ -1,0 +1,41 @@

+#include "platform.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);

+	}

+	b->sign = (i >> (sizeof(i)*8 - 1)) | 1;

+	i *= b->sign;

+	*b->p = i;

+	b->top = 1;

+	return mpnorm(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;

+}

--- /dev/null

+++ b/mp/mptoui.c

@@ -1,0 +1,31 @@

+#include "platform.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);

+	}

+	*b->p = i;

+	b->top = 1;

+	b->sign = 1;

+	return mpnorm(b);

+}

+

+uint

+mptoui(mpint *b)

+{

+	uint x;

+

+	x = *b->p;

+	if(b->sign < 0)

+		x = 0;

+	else if(b->top > 1 || (sizeof(mpdigit) > sizeof(uint) && x > MAXUINT))

+		x =  MAXUINT;

+	return x;

+}

--- /dev/null

+++ b/mp/mptouv.c

@@ -1,0 +1,44 @@

+#include "platform.h"

+

+#define VLDIGITS (int)(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*Dbits);

+	}else

+		mpbits(b, VLDIGITS*Dbits);

+	b->sign = 1;

+	for(s = 0; s < VLDIGITS; s++){

+		b->p[s] = v;

+		v >>= sizeof(mpdigit)*8;

+	}

+	b->top = s;

+	return mpnorm(b);

+}

+

+uvlong

+mptouv(mpint *b)

+{

+	uvlong v;

+	int s;

+

+	if(b->top == 0 || b->sign < 0)

+		return 0LL;

+

+	if(b->top > VLDIGITS)

+		return -1LL;

+

+	v = 0ULL;

+	for(s = 0; s < b->top; s++)

+		v |= (uvlong)b->p[s]<<(s*sizeof(mpdigit)*8);

+

+	return v;

+}

--- /dev/null

+++ b/mp/mptov.c

@@ -1,0 +1,60 @@

+#include "platform.h"

+

+#define VLDIGITS (int)(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*Dbits);

+	}else

+		mpbits(b, VLDIGITS*Dbits);

+	b->sign = (v >> (sizeof(v)*8 - 1)) | 1;

+	uv = v * b->sign;

+	for(s = 0; s < VLDIGITS; s++){

+		b->p[s] = uv;

+		uv >>= sizeof(mpdigit)*8;

+	}

+	b->top = s;

+	return mpnorm(b);

+}

+

+vlong

+mptov(mpint *b)

+{

+	uvlong v;

+	int s;

+

+	if(b->top == 0)

+		return 0LL;

+

+	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 |= (uvlong)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;

+}

--- /dev/null

+++ b/mp/mpvecadd.c

@@ -1,0 +1,34 @@

+#include "platform.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;

+	uint 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;

+}

--- /dev/null

+++ b/mp/mpveccmp.c

@@ -1,0 +1,25 @@

+#include "platform.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;

+}

--- /dev/null

+++ b/mp/mpvecdigmuladd.c

@@ -1,0 +1,101 @@

+#include "platform.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;

+}

--- /dev/null

+++ b/mp/mpvecsub.c

@@ -1,0 +1,32 @@

+#include "platform.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 < (mpdigit)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;

+	}

+}

--- /dev/null

+++ b/mp/mpvectscmp.c

@@ -1,0 +1,32 @@

+#include "platform.h"

+

+int

+mpvectscmp(mpdigit *a, int alen, mpdigit *b, int blen)

+{

+	mpdigit x, y, z, v;

+	int m, p;

+

+	if(alen > blen){

+		v = 0;

+		while(alen > blen)

+			v |= a[--alen];

+		m = p = (-v^v|v)>>(Dbits-1);

+	} else if(blen > alen){

+		v = 0;

+		while(blen > alen)

+			v |= b[--blen];

+		m = (-v^v|v)>>(Dbits-1);

+		p = m^1;

+	} else

+		m = p = 0;

+	while(alen-- > 0){

+		x = a[alen];

+		y = b[alen];

+		z = x - y;

+		x = ~x;

+		v = ((-z^z|z)>>(Dbits-1)) & ~m;

+		p = ((~(x&y|x&z|y&z)>>(Dbits-1)) & v) | (p & ~v);

+		m |= v;

+	}

+	return (p-m) | m;

+}

--- /dev/null

+++ b/mp/strtomp.c

@@ -1,0 +1,174 @@

+#include "platform.h"

+

+static char*

+frompow2(char *a, mpint *b, int s)

+{

+	char *p, *next;

+	mpdigit x;

+	int i;

+

+	i = 1<<s;

+	for(p = a; (dec16chr(*p) & 255) < i; p++)

+		;

+

+	mpbits(b, (p-a)*s);

+	b->top = 0;

+	next = p;

+

+	while(p > a){

+		x = 0;

+		for(i = 0; i < Dbits; i += s){

+			if(p <= a)

+				break;

+			x |= dec16chr(*--p)<<i;

+		}

+		b->p[b->top++] = x;

+	}

+	return next;

+}

+

+static char*

+from8(char *a, mpint *b)

+{

+	char *p, *next;

+	mpdigit x, y;

+	int i;

+

+	for(p = a; ((*p - '0') & 255) < 8; p++)

+		;

+

+	mpbits(b, (p-a)*3);

+	b->top = 0;

+	next = p;

+

+	i = 0;

+	x = y = 0;

+	while(p > a){

+		y = *--p - '0';

+		x |= y << i;

+		i += 3;

+		if(i >= Dbits){

+Digout:

+			i -= Dbits;

+			b->p[b->top++] = x;

+			x = y >> (3-i);

+		}

+	}

+	if(i > 0)

+		goto Digout;

+

+	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 = *a - '0';

+			if(y > 9)

+				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;

+}

+

+mpint*

+strtomp(char *a, char **pp, int base, mpint *b)

+{

+	int sign;

+	char *e;

+

+	if(b == nil){

+		b = mpnew(0);

+	}

+

+	while(*a==' ' || *a=='\t')

+		a++;

+

+	sign = 1;

+	for(;; a++){

+		switch(*a){

+		case '-':

+			sign *= -1;

+			continue;

+		}

+		break;

+	}

+

+	if(base == 0){

+		base = 10;

+		if(a[0] == '0'){

+			if(a[1] == 'x' || a[1] == 'X') {

+				a += 2;

+				base = 16;

+			} else if(a[1] == 'b' || a[1] == 'B') {

+				a += 2;

+				base = 2;

+			} else if(a[1] >= '0' && a[1] <= '7') {

+				a++;

+				base = 8;

+			}

+		}

+	}

+

+	switch(base){

+	case 2:

+		e = frompow2(a, b, 1);

+		break;

+	case 4:

+		e = frompow2(a, b, 2);

+		break;

+	case 8:

+		e = from8(a, b);

+		break;

+	case 10:

+		e = from10(a, b);

+		break;

+	case 16:

+		e = frompow2(a, b, 4);

+		break;

+	default:

+		abort();

+		return nil;

+	}

+

+	if(pp != nil)

+		*pp = e;

+

+	// if no characters parsed, there wasn't a number to convert

+	if(e == a)

+		return nil;

+

+	b->sign = sign;

+	return mpnorm(b);

+}

--- /dev/null

+++ b/mp/u16.c

@@ -1,0 +1,68 @@

+#include "platform.h"

+

+#define between(x,min,max)	(((min-1-x) & (x-max-1))>>8)

+

+int

+enc16chr(int o)

+{

+	int c;

+

+	c  = between(o,  0,  9) & ('0'+o);

+	c |= between(o, 10, 15) & ('A'+(o-10));

+	return c;

+}

+

+int

+dec16chr(int c)

+{

+	int o;

+

+	o  = between(c, '0', '9') & (1+(c-'0'));

+	o |= between(c, 'A', 'F') & (1+10+(c-'A'));

+	o |= between(c, 'a', 'f') & (1+10+(c-'a'));

+	return o-1;

+}

+

+int

+dec16(uchar *out, int lim, char *in, int n)

+{

+	int c, w = 0, i = 0;

+	uchar *start = out;

+	uchar *eout = out + lim;

+

+	while(n-- > 0){

+		c = dec16chr(*in++);

+		if(c < 0)

+			continue;

+		w = (w<<4) + c;

+		i++;

+		if(i == 2){

+			if(out + 1 > eout)

+				goto exhausted;

+			*out++ = w;

+			w = 0;

+			i = 0;

+		}

+	}

+exhausted:

+	return out - start;

+}

+

+int

+enc16(char *out, int lim, uchar *in, int n)

+{

+	uint c;

+	char *eout = out + lim;

+	char *start = out;

+

+	while(n-- > 0){

+		c = *in++;

+		if(out + 2 >= eout)

+			goto exhausted;

+		*out++ = enc16chr(c>>4);

+		*out++ = enc16chr(c&15);

+	}

+exhausted:

+	*out = 0;

+	return out - start;

+}

--- /dev/null

+++ b/mp/u32.c

@@ -1,0 +1,143 @@

+#include "platform.h"

+

+#define between(x,min,max)	(((min-1-x) & (x-max-1))>>8)

+

+int

+enc32chr(int o)

+{

+	int c;

+

+	c  = between(o,  0, 25) & ('A'+o);

+	c |= between(o, 26, 31) & ('2'+(o-26));

+	return c;

+}

+

+int

+dec32chr(int c)

+{

+	int o;

+

+	o  = between(c, 'A', 'Z') & (1+(c-'A'));

+	o |= between(c, 'a', 'z') & (1+(c-'a'));

+	o |= between(c, '2', '7') & (1+26+(c-'2'));

+	return o-1;

+}

+

+int

+dec32x(uchar *dest, int ndest, char *src, int nsrc, int (*chr)(int))

+{

+	uchar *start;

+	int i, j, u[8];

+

+	if(ndest+1 < (5*nsrc+7)/8)

+		return -1;

+	start = dest;

+	while(nsrc>=8){

+		for(i=0; i<8; i++){

+			j = chr(src[i]);

+			if(j < 0)

+				j = 0;

+			u[i] = j;

+		}

+		*dest++ = (u[0]<<3) | (0x7 & (u[1]>>2));

+		*dest++ = ((0x3 & u[1])<<6) | (u[2]<<1) | (0x1 & (u[3]>>4));

+		*dest++ = ((0xf & u[3])<<4) | (0xf & (u[4]>>1));

+		*dest++ = ((0x1 & u[4])<<7) | (u[5]<<2) | (0x3 & (u[6]>>3));

+		*dest++ = ((0x7 & u[6])<<5) | u[7];

+		src  += 8;

+		nsrc -= 8;

+	}

+	if(nsrc > 0){

+		if(nsrc == 1 || nsrc == 3 || nsrc == 6)

+			return -1;

+		for(i=0; i<nsrc; i++){

+			j = chr(src[i]);

+			if(j < 0)

+				j = 0;

+			u[i] = j;

+		}

+		*dest++ = (u[0]<<3) | (0x7 & (u[1]>>2));

+		if(nsrc == 2)

+			goto out;

+		*dest++ = ((0x3 & u[1])<<6) | (u[2]<<1) | (0x1 & (u[3]>>4));

+		if(nsrc == 4)

+			goto out;

+		*dest++ = ((0xf & u[3])<<4) | (0xf & (u[4]>>1));

+		if(nsrc == 5)

+			goto out;

+		*dest++ = ((0x1 & u[4])<<7) | (u[5]<<2) | (0x3 & (u[6]>>3));

+	}

+out:

+	return dest-start;

+}

+

+int

+enc32x(char *dest, int ndest, uchar *src, int nsrc, int (*chr)(int))

+{

+	char *start;

+	int j;

+

+	if(ndest <= (8*nsrc+4)/5)

+		return -1;

+	start = dest;

+	while(nsrc>=5){

+		j = (0x1f & (src[0]>>3));

+		*dest++ = chr(j);

+		j = (0x1c & (src[0]<<2)) | (0x03 & (src[1]>>6));

+		*dest++ = chr(j);

+		j = (0x1f & (src[1]>>1));

+		*dest++ = chr(j);

+		j = (0x10 & (src[1]<<4)) | (0x0f & (src[2]>>4));

+		*dest++ = chr(j);

+		j = (0x1e & (src[2]<<1)) | (0x01 & (src[3]>>7));

+		*dest++ = chr(j);

+		j = (0x1f & (src[3]>>2));

+		*dest++ = chr(j);

+		j = (0x18 & (src[3]<<3)) | (0x07 & (src[4]>>5));

+		*dest++ = chr(j);

+		j = (0x1f & (src[4]));

+		*dest++ = chr(j);

+		src  += 5;

+		nsrc -= 5;

+	}

+	if(nsrc){

+		j = (0x1f & (src[0]>>3));

+		*dest++ = chr(j);

+		j = (0x1c & (src[0]<<2));

+		if(nsrc == 1)

+			goto out;

+		j |= (0x03 & (src[1]>>6));

+		*dest++ = chr(j);

+		j = (0x1f & (src[1]>>1));

+		*dest++ = chr(j);

+		j = (0x10 & (src[1]<<4));

+		if(nsrc == 2)

+			goto out;

+		j |= (0x0f & (src[2]>>4));

+		*dest++ = chr(j);

+		j = (0x1e & (src[2]<<1));

+		if(nsrc == 3)

+			goto out;

+		j |= (0x01 & (src[3]>>7));

+		*dest++ = chr(j);

+		j = (0x1f & (src[3]>>2));

+		*dest++ = chr(j);

+		j = (0x18 & (src[3]<<3));

+out:

+		*dest++ = chr(j);

+	}

+	*dest = 0;

+	return dest-start;

+}

+

+int

+enc32(char *dest, int ndest, uchar *src, int nsrc)

+{

+	return enc32x(dest, ndest, src, nsrc, enc32chr);

+}

+

+int

+dec32(uchar *dest, int ndest, char *src, int nsrc)

+{

+	return dec32x(dest, ndest, src, nsrc, dec32chr);

+}

--- /dev/null

+++ b/mp/u64.c

@@ -1,0 +1,141 @@

+#include "platform.h"

+

+#define between(x,min,max)	(((min-1-x) & (x-max-1))>>8)

+

+int

+enc64chr(int o)

+{

+	int c;

+

+	c  = between(o,  0, 25) & ('A'+o);

+	c |= between(o, 26, 51) & ('a'+(o-26));

+	c |= between(o, 52, 61) & ('0'+(o-52));

+	c |= between(o, 62, 62) & ('+');

+	c |= between(o, 63, 63) & ('/');

+	return c;

+}

+

+int

+dec64chr(int c)

+{

+	int o;

+

+	o  = between(c, 'A', 'Z') & (1+(c-'A'));

+	o |= between(c, 'a', 'z') & (1+26+(c-'a'));

+	o |= between(c, '0', '9') & (1+52+(c-'0'));

+	o |= between(c, '+', '+') & (1+62);

+	o |= between(c, '/', '/') & (1+63);

+	return o-1;

+}

+

+int

+dec64x(uchar *out, int lim, char *in, int n, int (*chr)(int))

+{

+	ulong b24;

+	uchar *start = out;

+	uchar *e = out + lim;

+	int i, c;

+

+	b24 = 0;

+	i = 0;

+	while(n-- > 0){

+		c = chr(*in++);

+		if(c < 0)

+			continue;

+		switch(i){

+		case 0:

+			b24 = c<<18;

+			break;

+		case 1:

+			b24 |= c<<12;

+			break;

+		case 2:

+			b24 |= c<<6;

+			break;

+		case 3:

+			if(out + 3 > e)

+				goto exhausted;

+

+			b24 |= c;

+			*out++ = b24>>16;

+			*out++ = b24>>8;

+			*out++ = b24;

+			i = 0;

+			continue;

+		}

+		i++;

+	}

+	switch(i){

+	case 2:

+		if(out + 1 > e)

+			goto exhausted;

+		*out++ = b24>>16;

+		break;

+	case 3:

+		if(out + 2 > e)

+			goto exhausted;

+		*out++ = b24>>16;

+		*out++ = b24>>8;

+		break;

+	}

+exhausted:

+	return out - start;

+}

+

+int

+enc64x(char *out, int lim, uchar *in, int n, int (*chr)(int))

+{

+	int i;

+	ulong b24;

+	char *start = out;

+	char *e = out + lim;

+

+	for(i = n/3; i > 0; i--){

+		b24 = *in++<<16;

+		b24 |= *in++<<8;

+		b24 |= *in++;

+		if(out + 4 >= e)

+			goto exhausted;

+		*out++ = chr(b24>>18);

+		*out++ = chr((b24>>12)&0x3f);

+		*out++ = chr((b24>>6)&0x3f);

+		*out++ = chr(b24&0x3f);

+	}

+

+	switch(n%3){

+	case 2:

+		b24 = *in++<<16;

+		b24 |= *in<<8;

+		if(out + 4 >= e)

+			goto exhausted;

+		*out++ = chr(b24>>18);

+		*out++ = chr((b24>>12)&0x3f);

+		*out++ = chr((b24>>6)&0x3f);

+		*out++ = '=';

+		break;

+	case 1:

+		b24 = *in<<16;

+		if(out + 4 >= e)

+			goto exhausted;

+		*out++ = chr(b24>>18);

+		*out++ = chr((b24>>12)&0x3f);

+		*out++ = '=';

+		*out++ = '=';

+		break;

+	}

+exhausted:

+	*out = 0;

+	return out - start;

+}

+

+int

+enc64(char *out, int lim, uchar *in, int n)

+{

+	return enc64x(out, lim, in, n, enc64chr);

+}

+

+int

+dec64(uchar *out, int lim, char *in, int n)

+{

+	return dec64x(out, lim, in, n, dec64chr);

+}