ref: 92980154dee84747ca798bafe6e05d87519ab665
parent: 1b0eb3b4437a74231451a95225233fada1039876
author: Sigrid Solveig Haflínudóttir <[email protected]>
date: Sat Apr 1 11:46:00 EDT 2023
don't do recursive makefiling
--- 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}/binTARG=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);
+}