shithub: femtolisp

Info  •  Files  •  Log  •  Branches

initial (unfinished) implementation of mpint number type

--- a/.gitignore

+++ b/.gitignore

@@ -13,3 +13,4 @@

 instructions.lsp

 builtins.lsp

 builtin_fns.h

+*.core

--- a/LICENSE

+++ b/LICENSE

@@ -1,3 +1,5 @@

+9front portions (especially libmp) is under MIT license.

+

 Copyright (c) 2008 Jeff Bezanson

 All rights reserved.

--- a/Makefile

+++ b/Makefile

@@ -5,7 +5,7 @@

 TARG=flisp

 LLT=llt/libllt.a

 CFLAGS?=-O2 -g

-CFLAGS+=-Wall -Wextra -falign-functions -Wno-strict-aliasing -std=c99 -D_DEFAULT_SOURCE

+CFLAGS+=-Wall -Wextra -Wno-shift-op-parentheses -Wno-bitwise-op-parentheses -falign-functions -Wno-strict-aliasing -std=c99 -D_DEFAULT_SOURCE

 LDFLAGS?=

 OBJS=\

--- a/README.md

+++ b/README.md

@@ -19,6 +19,7 @@

  * "boot" image is built into the executable

  * vm opcode definitions and tables are generated from a single file

  * fixed bootstrap (makes it work properly when opcodes change)

+ * bigints (though not completely finished yet)

 ## Characteristics

--- a/cvalues.c

+++ b/cvalues.c

@@ -7,7 +7,7 @@

 #endif

 value_t int8sym, uint8sym, int16sym, uint16sym, int32sym, uint32sym;

-value_t int64sym, uint64sym;

+value_t int64sym, uint64sym, mpintsym;

 value_t longsym, ulongsym, bytesym, wcharsym;

 value_t floatsym, doublesym;

 value_t gftypesym, stringtypesym, wcstringtypesym;

@@ -23,6 +23,7 @@

 static fltype_t *int32type, *uint32type;

 static fltype_t *int64type, *uint64type;

 static fltype_t *longtype, *ulongtype;

+static fltype_t *mpinttype;

 static fltype_t *floattype, *doubletype;

        fltype_t *bytetype, *wchartype;

        fltype_t *stringtype, *wcstringtype;

@@ -127,7 +128,7 @@

     cvalue_t *pcv;

     int str=0;

-    if (valid_numtype(type->numtype))

+    if (valid_numtype(type->numtype) && type->numtype != T_MPINT)

         return cprim(type, sz);

     if (type->eltype == bytetype) {

@@ -308,6 +309,50 @@

 num_ctor(float, float, T_FLOAT)

 num_ctor(double, double, T_DOUBLE)

+static int cvalue_mpint_init(fltype_t *type, value_t arg, void *dest)

+{

+    mpint *n;

+    USED(type);

+    if (isfixnum(arg)) {

+        n = vtomp(numval(arg), nil);

+    }

+    else if (iscprim(arg)) {

+        cprim_t *cp = (cprim_t*)ptr(arg);

+        void *p = cp_data(cp);

+        n = conv_to_mpint(p, cp_numtype(cp));

+    }

+    else {

+        return 1;

+    }

+    *((mpint**)dest) = n;

+    return 0;

+}

+

+/* */ BUILTIN("mpint", mpint)

+{

+    if (nargs==0) { PUSH(fixnum(0)); args = &Stack[SP-1]; }

+    value_t cv = cvalue(mpinttype, sizeof(mpint*));

+    if (cvalue_mpint_init(mpinttype, args[0], cv_data((cvalue_t*)ptr(cv))))

+        type_error("number", args[0]);

+    return cv;

+}

+

+value_t mk_mpint(mpint *n)

+{

+    value_t cv = cvalue(mpinttype, sizeof(mpint*));

+    *(mpint**)cvalue_data(cv) = n;

+    return cv;

+}

+

+static void free_mpint(value_t self)

+{

+    mpint **s = value2c(mpint**, self);

+    if (*s != mpzero && *s != mpone && *s != mptwo)

+        mpfree(*s);

+}

+

+static cvtable_t mpint_vtable = { nil, nil, free_mpint, nil };

+

 value_t size_wrap(size_t sz)

 {

     if (fits_fixnum(sz))

@@ -354,7 +399,7 @@

                 return 0;

             }

         }

-        lerrorf(ArgError, "enum: invalid enum value");

+        lerrorf(ArgError, "invalid enum value");

     }

     if (isfixnum(arg)) {

         n = (int)numval(arg);

@@ -367,7 +412,7 @@

         type_error("number", arg);

     }

     if ((unsigned)n >= vector_size(syms))

-        lerrorf(ArgError, "enum: value out of range");

+        lerrorf(ArgError, "value out of range");

     *(int*)dest = n;

     return 0;

 }

@@ -435,7 +480,7 @@

             arg = cdr_(arg);

         }

         if (i != cnt)

-            lerrorf(ArgError, "array: size mismatch");

+            lerrorf(ArgError, "size mismatch");

         return 0;

     }

     else if (iscvalue(arg)) {

@@ -446,12 +491,12 @@

                 if (cv_len(cv) == sz)

                     memmove(dest, cv_data(cv), sz);

                 else

-                    lerrorf(ArgError, "array: size mismatch");

+                    lerrorf(ArgError, "size mismatch");

                 return 0;

             }

             else {

                 // TODO: initialize array from different type elements

-                lerrorf(ArgError, "array: element type mismatch");

+                lerrorf(ArgError, "element type mismatch");

             }

         }

     }

@@ -553,7 +598,7 @@

         if (hed == arraysym) {

             value_t t = car(cdr_(type));

             if (!iscons(cdr_(cdr_(type))))

-                lerrorf(ArgError, "sizeof: incomplete type");

+                lerrorf(ArgError, "incomplete type");

             value_t n = car_(cdr_(cdr_(type)));

             size_t sz = toulong(n);

             return sz * ctype_sizeof(t, palign);

@@ -570,7 +615,7 @@

         }

     }

-    lerrorf(ArgError, "sizeof: invalid c type");

+    lerrorf(ArgError, "invalid c type");

 }

 extern fltype_t *iostreamtype;

@@ -687,11 +732,11 @@

 {

     argcount(nargs, 1);

     if (iscons(args[0]) || isvector(args[0]))

-        lerrorf(ArgError, "copy: argument must be a leaf atom");

+        lerrorf(ArgError, "argument must be a leaf atom");

     if (!iscvalue(args[0]))

         return args[0];

     if (!cv_isPOD((cvalue_t*)ptr(args[0])))

-        lerrorf(ArgError, "copy: argument must be a plain-old-data type");

+        lerrorf(ArgError, "argument must be a plain-old-data type");

     return cvalue_copy(args[0]);

 }

@@ -708,7 +753,7 @@

     cvinitfunc_t f=type->init;

     if (f == nil)

-        lerrorf(ArgError, "c-value: invalid c type");

+        lerrorf(ArgError, "invalid c type");

     f(type, v, dest);

 }

@@ -912,7 +957,12 @@

     double Faccum=0;

     int32_t inexact = 0;

     uint32_t i;

+    int64_t i64;

     value_t arg;

+    mpint *Maccum = nil, *x;

+    numerictype_t pt;

+    fixnum_t pi;

+    void *a;

     FOR_ARGS(i,0,arg,args) {

         if (isfixnum(arg)) {

@@ -919,17 +969,14 @@

             Saccum += numval(arg);

             continue;

         }

-        else if (iscprim(arg)) {

-            cprim_t *cp = (cprim_t*)ptr(arg);

-            void *a = cp_data(cp);

-            int64_t i64;

-            switch(cp_numtype(cp)) {

+        if (num_to_ptr(arg, &pi, &pt, &a)) {

+            switch(pt) {

             case T_INT8:   Saccum += *(int8_t*)a; break;

-            case T_UINT8:  Saccum += *(uint8_t*)a; break;

+            case T_UINT8:  Uaccum += *(uint8_t*)a; break;

             case T_INT16:  Saccum += *(int16_t*)a; break;

-            case T_UINT16: Saccum += *(uint16_t*)a; break;

+            case T_UINT16: Uaccum += *(uint16_t*)a; break;

             case T_INT32:  Saccum += *(int32_t*)a; break;

-            case T_UINT32: Saccum += *(uint32_t*)a; break;

+            case T_UINT32: Uaccum += *(uint32_t*)a; break;

             case T_INT64:

                 i64 = *(int64_t*)a;

                 if (i64 > 0)

@@ -938,6 +985,11 @@

                     Saccum += i64;

                 break;

             case T_UINT64: Uaccum += *(uint64_t*)a; break;

+            case T_MPINT:

+                if (Maccum == nil)

+                    Maccum = mpnew(0);

+                mpadd(Maccum, *(mpint**)a, Maccum);

+                break;

             case T_FLOAT:  Faccum += *(float*)a; inexact = 1; break;

             case T_DOUBLE: Faccum += *(double*)a; inexact = 1; break;

             default:

@@ -945,15 +997,30 @@

             }

             continue;

         }

-    add_type_error:

+

+add_type_error:

+        mpfree(Maccum);

         type_error("number", arg);

     }

     if (inexact) {

         Faccum += Uaccum;

         Faccum += Saccum;

+        if (Maccum != nil) {

+            Faccum += mptod(Maccum);

+            mpfree(Maccum);

+        }

         return mk_double(Faccum);

     }

-    else if (Saccum < 0) {

+    if (Maccum != nil) {

+        /* FIXME - check if it fits into fixnum first, etc */

+        x = vtomp(Saccum, nil);

+        mpadd(Maccum, x, Maccum);

+        x = uvtomp(Uaccum, x);

+        mpadd(Maccum, x, Maccum);

+        mpfree(x);

+        return mk_mpint(Maccum);

+    }

+    if (Saccum < 0) {

         uint64_t negpart = (uint64_t)(-Saccum);

         if (negpart > Uaccum) {

             Saccum += (int64_t)Uaccum;

@@ -977,20 +1044,23 @@

 static value_t fl_neg(value_t n)

 {

+    uint32_t ui32;

+    int32_t i32;

+    int64_t i64;

+    mpint *mp;

+    numerictype_t pt;

+    fixnum_t pi;

+    void *a;

+

     if (isfixnum(n)) {

         fixnum_t s = fixnum(-numval(n));

         if (__unlikely((value_t)s == n))

             return mk_xlong(-numval(n));

-        else

-            return s;

+        return s;

     }

-    else if (iscprim(n)) {

-        cprim_t *cp = (cprim_t*)ptr(n);

-        void *a = cp_data(cp);

-        uint32_t ui32;

-        int32_t i32;

-        int64_t i64;

-        switch(cp_numtype(cp)) {

+

+    if (num_to_ptr(n, &pi, &pt, &a)) {

+        switch(pt) {

         case T_INT8:   return fixnum(-(int32_t)*(int8_t*)a);

         case T_UINT8:  return fixnum(-(int32_t)*(uint8_t*)a);

         case T_INT16:  return fixnum(-(int32_t)*(int16_t*)a);

@@ -1010,11 +1080,15 @@

                 return mk_uint64((uint64_t)BIT63);

             return mk_int64(-i64);

         case T_UINT64: return mk_int64(-(int64_t)*(uint64_t*)a);

+        case T_MPINT:

+            mp = mpcopy(*(mpint**)a);

+            mpsub(mpzero, mp, mp);

+            return mk_mpint(mp);

         case T_FLOAT:  return mk_float(-*(float*)a);

         case T_DOUBLE: return mk_double(-*(double*)a);

-            break;

         }

     }

+

     type_error("number", n);

 }

@@ -1023,8 +1097,13 @@

     uint64_t Uaccum=1;

     double Faccum=1;

     int32_t inexact = 0;

+    int64_t i64;

     uint32_t i;

     value_t arg;

+    mpint *Maccum=nil, *x;

+    numerictype_t pt;

+    fixnum_t pi;

+    void *a;

     FOR_ARGS(i,0,arg,args) {

         if (isfixnum(arg)) {

@@ -1031,17 +1110,14 @@

             Saccum *= numval(arg);

             continue;

         }

-        else if (iscprim(arg)) {

-            cprim_t *cp = (cprim_t*)ptr(arg);

-            void *a = cp_data(cp);

-            int64_t i64;

-            switch(cp_numtype(cp)) {

+        if (num_to_ptr(arg, &pi, &pt, &a)) {

+            switch(pt) {

             case T_INT8:   Saccum *= *(int8_t*)a; break;

-            case T_UINT8:  Saccum *= *(uint8_t*)a; break;

+            case T_UINT8:  Uaccum *= *(uint8_t*)a; break;

             case T_INT16:  Saccum *= *(int16_t*)a; break;

-            case T_UINT16: Saccum *= *(uint16_t*)a; break;

+            case T_UINT16: Uaccum *= *(uint16_t*)a; break;

             case T_INT32:  Saccum *= *(int32_t*)a; break;

-            case T_UINT32: Saccum *= *(uint32_t*)a; break;

+            case T_UINT32: Uaccum *= *(uint32_t*)a; break;

             case T_INT64:

                 i64 = *(int64_t*)a;

                 if (i64 > 0)

@@ -1050,6 +1126,11 @@

                     Saccum *= i64;

                 break;

             case T_UINT64: Uaccum *= *(uint64_t*)a; break;

+            case T_MPINT:

+                if (Maccum == nil)

+                    Maccum = mpcopy(mpone);

+                mpmul(Maccum, *(mpint**)a, Maccum);

+                break;

             case T_FLOAT:  Faccum *= *(float*)a; inexact = 1; break;

             case T_DOUBLE: Faccum *= *(double*)a; inexact = 1; break;

             default:

@@ -1057,15 +1138,29 @@

             }

             continue;

         }

-    mul_type_error:

+

+mul_type_error:

         type_error("number", arg);

     }

     if (inexact) {

         Faccum *= Uaccum;

         Faccum *= Saccum;

+        if (Maccum != nil) {

+            Faccum *= mptod(Maccum);

+            mpfree(Maccum);

+        }

         return mk_double(Faccum);

     }

-    else if (Saccum < 0) {

+    if (Maccum != nil) {

+        /* FIXME might still fit into a fixnum */

+        x = vtomp(Saccum, nil);

+        mpmul(Maccum, x, Maccum);

+        x = uvtomp(Uaccum, x);

+        mpmul(Maccum, x, Maccum);

+        mpfree(x);

+        return mk_mpint(Maccum);

+    }

+    if (Saccum < 0) {

         Saccum *= (int64_t)Uaccum;

         if (Saccum >= INT32_MIN) {

             if (fits_fixnum(Saccum)) {

@@ -1081,9 +1176,10 @@

     return return_from_uint64(Uaccum);

 }

-static int num_to_ptr(value_t a, fixnum_t *pi, numerictype_t *pt, void **pp)

+int num_to_ptr(value_t a, fixnum_t *pi, numerictype_t *pt, void **pp)

 {

     cprim_t *cp;

+    cvalue_t *cv;

     if (isfixnum(a)) {

         *pi = numval(a);

         *pp = pi;

@@ -1094,6 +1190,12 @@

         *pp = cp_data(cp);

         *pt = cp_numtype(cp);

     }

+    else if (iscvalue(a)) {

+        cv = (cvalue_t*)ptr(a);

+        *pp = cv_data(cv);

+        *pt = cv_class(cv)->numtype;

+        return valid_numtype(*pt);

+    }

     else {

         return 0;

     }

@@ -1465,6 +1567,11 @@

     mk_primtype(wchar, int32_t);

     mk_primtype(float, float);

     mk_primtype(double, double);

+

+    ctor_cv_intern(mpint, T_MPINT, mpint*);

+    mpinttype = get_type(mpintsym);

+    mpinttype->init = cvalue_mpint_init;

+    mpinttype->vtable = &mpint_vtable;

     stringtype = get_type(symbol_value(stringtypesym));

     wcstringtype = get_type(symbol_value(wcstringtypesym));

--- a/equal.c

+++ b/equal.c

@@ -55,6 +55,7 @@

 static value_t bounded_compare(value_t a, value_t b, int bound, int eq)

 {

     value_t d;

+    cvalue_t *cv;

  compare_top:

     if (a == b) return fixnum(0);

@@ -74,6 +75,11 @@

                 return fixnum(1);

             return fixnum(numeric_compare(a, b, eq, 1, 0));

         }

+        if (iscvalue(b)) {

+            cv = ptr(b);

+            if (valid_numtype(cv_class(cv)->numtype))

+                return fixnum(numeric_compare(a, b, eq, 1, 0));

+        }

         return fixnum(-1);

     case TAG_SYM:

         if (eq) return fixnum(1);

@@ -97,6 +103,11 @@

             return fixnum(c);

         break;

     case TAG_CVALUE:

+        cv = ptr(a);

+        if (valid_numtype(cv_class(cv)->numtype)) {

+            if((c = numeric_compare(a, b, eq, 1, 0)) != 2)

+                return fixnum(c);

+        }

         if (iscvalue(b)) {

             if (cv_isPOD((cvalue_t*)ptr(a)) && cv_isPOD((cvalue_t*)ptr(b)))

                 return cvalue_compare(a, b);

@@ -329,7 +340,14 @@

     case TAG_CVALUE:

         cv = (cvalue_t*)ptr(a);

         data = cv_data(cv);

-        return memhash(data, cv_len(cv));

+        if (cv->type == mpinttype) {

+            len = mptobe(*(mpint**)data, nil, 0, (uint8_t**)&data);

+            h = memhash(data, len);

+            free(data);

+        } else {

+            h = memhash(data, cv_len(cv));

+        }

+        return h;

     case TAG_VECTOR:

         if (bound <= 0) {

--- a/flisp.c

+++ b/flisp.c

@@ -711,11 +711,16 @@

 int fl_isnumber(value_t v)

 {

-    if (isfixnum(v)) return 1;

+    if (isfixnum(v))

+        return 1;

     if (iscprim(v)) {

-        cprim_t *c = (cprim_t*)ptr(v);

+        cprim_t *c = ptr(v);

         return c->type != wchartype;

     }

+    if (iscvalue(v)) {

+        cvalue_t *c = ptr(v);

+        return valid_numtype(cv_class(c)->numtype);

+    }

     return 0;

 }

@@ -813,9 +818,9 @@

     value_t s4 = Stack[SP-4];

     value_t s5 = Stack[SP-5];

     if (nargs < nreq)

-        lerrorf(ArgError, "apply: too few arguments");

+        lerrorf(ArgError, "too few arguments");

     if (extr > nelem(args))

-        lerrorf(ArgError, "apply: too many arguments");

+        lerrorf(ArgError, "too many arguments");

     for (i=0; i < extr; i++) args[i] = UNBOUND;

     for (i=nreq; i < nargs; i++) {

         v = Stack[bp+i];

@@ -855,7 +860,7 @@

  no_kw:

     nrestargs = nargs - i;

     if (!va && nrestargs > 0)

-        lerrorf(ArgError, "apply: too many arguments");

+        lerrorf(ArgError, "too many arguments");

     nargs = ntot + nrestargs;

     if (nrestargs)

         memmove(&Stack[bp+ntot], &Stack[bp+i], nrestargs*sizeof(value_t));

@@ -1722,7 +1727,7 @@

                 }

             }

             else if (s < 0) {

-                lerrorf(ArgError, "apply: too few arguments");

+                lerrorf(ArgError, "too few arguments");

             }

             else {

                 PUSH(0);

@@ -1745,10 +1750,10 @@

             i = GET_INT32(ip); ip+=4;

             n = GET_INT32(ip); ip+=4;

             if (nargs < i)

-                lerrorf(ArgError, "apply: too few arguments");

+                lerrorf(ArgError, "too few arguments");

             if ((int32_t)n > 0) {

                 if (nargs > n)

-                    lerrorf(ArgError, "apply: too many arguments");

+                    lerrorf(ArgError, "too many arguments");

             }

             else n = -n;

             if (n > nargs) {

@@ -1834,6 +1839,14 @@

 // builtins -------------------------------------------------------------------

+BUILTIN("gc", gc)

+{

+    USED(args);

+    argcount(nargs, 0);

+    gc(0);

+    return FL_T;

+}

+

 BUILTIN("function", function)

 {

     if (nargs == 1 && issymbol(args[0]))

@@ -1884,7 +1897,7 @@

             }

         }

         if (isgensym(fn->name))

-            lerrorf(ArgError, "function: name should not be a gensym");

+            lerrorf(ArgError, "name should not be a gensym");

     }

     return fv;

 }

@@ -1967,7 +1980,7 @@

 BUILTIN("map", map)

 {

     if (nargs < 2)

-        lerrorf(ArgError, "map: too few arguments");

+        lerrorf(ArgError, "too few arguments");

     if (!iscons(args[1])) return NIL;

     value_t first, last, v;

     int64_t argSP = args-Stack;

--- a/flisp.h

+++ b/flisp.h

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

     T_INT16, T_UINT16,

     T_INT32, T_UINT32,

     T_INT64, T_UINT64,

+    T_MPINT,

     T_FLOAT,

     T_DOUBLE,

 } numerictype_t;

@@ -98,6 +99,8 @@

 // doesn't lead to other values

 #define leafp(a) (((a)&3) != 3)

+int num_to_ptr(value_t a, fixnum_t *pi, numerictype_t *pt, void **pp);

+

 #define isforwarded(v) (((value_t*)ptr(v))[0] == TAG_FWD)

 #define forwardloc(v)  (((value_t*)ptr(v))[1])

 #define forward(v,to) do { (((value_t*)ptr(v))[0] = TAG_FWD); \

@@ -368,6 +371,7 @@

 double conv_to_double(void *data, numerictype_t tag);

 void conv_from_double(void *data, double d, numerictype_t tag);

+mpint *conv_to_mpint(void *data, numerictype_t tag);

 int64_t conv_to_int64(void *data, numerictype_t tag);

 uint64_t conv_to_uint64(void *data, numerictype_t tag);

 int32_t conv_to_int32(void *data, numerictype_t tag);

--- a/llt/Makefile

+++ b/llt/Makefile

@@ -14,6 +14,33 @@

 	random.o\

 	timefuncs.o\

 	utf8.o\

+    \

+	mpadd.o\

+	mpaux.o\

+	mpvecsub.o\

+	mpcmp.o\

+	mpdiv.o\

+	mpfmt.o\

+	mpmul.o\

+	mpsub.o\

+	mpleft.o\

+	mpright.o\

+	mptobe.o\

+	mptod.o\

+	mptoi.o\

+	mptoui.o\

+	mptouv.o\

+	mpdigdiv.o\

+	mpvecdigmuladd.o\

+	mptov.o\

+	mpvecadd.o\

+	strtomp.o\

+	u16.o\

+	u32.o\

+	u64.o\

+    mptober.o\

+    mpveccmp.o\

+    mpvectscmp.o\

 .PHONY: all default clean

--- /dev/null

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

+}

--- a/mkfile

+++ b/mkfile

@@ -6,10 +6,12 @@

 CLEANFILES=boot.h builtin_fns.h

 HFILES=\

-	cvalues.c\

-	equal.c\

 	equalhash.h\

 	flisp.h\

+

+CHFILES=\

+	cvalues.c\

+	equal.c\

 	operators.c\

 	print.c\

 	read.c\

@@ -41,7 +43,7 @@

 flmain.$O: boot.h

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

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

 bootstrap:V: $O.out

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

--- a/operators.c

+++ b/operators.c

@@ -1,5 +1,24 @@

 #include "llt.h"

+mpint *conv_to_mpint(void *data, numerictype_t tag)

+{

+    mpint *i = mpzero;

+    switch (tag) {

+    case T_INT8:   i = itomp(*(int8_t*)data, nil); break;

+    case T_UINT8:  i = uitomp(*(uint8_t*)data, nil); break;

+    case T_INT16:  i = itomp(*(int16_t*)data, nil); break;

+    case T_UINT16: i = uitomp(*(uint16_t*)data, nil); break;

+    case T_INT32:  i = itomp(*(int32_t*)data, nil); break;

+    case T_UINT32: i = uitomp(*(uint32_t*)data, nil); break;

+    case T_INT64:  i = vtomp(*(int64_t*)data, nil); break;

+    case T_UINT64: i = uvtomp(*(int64_t*)data, nil); break;

+    case T_MPINT:  i = mpcopy(*(mpint**)data); break;

+    case T_FLOAT:  i = dtomp(*(float*)data, nil); break;

+    case T_DOUBLE: i = dtomp(*(double*)data, nil); break;

+    }

+    return i;

+}

+

 double conv_to_double(void *data, numerictype_t tag)

 {

     double d=0;

@@ -16,6 +35,7 @@

             d = -d;

         break;

     case T_UINT64: d = (double)*(uint64_t*)data; break;

+    case T_MPINT:  d = mptod(*(mpint**)data); break;

     case T_FLOAT:  d = (double)*(float*)data; break;

     case T_DOUBLE: return *(double*)data;

     }

@@ -37,11 +57,13 @@

             *(int64_t*)dest = INT64_MAX;

         break;

     case T_UINT64: *(uint64_t*)dest = (int64_t)d; break;

+    case T_MPINT:  *(mpint**)dest = dtomp(d, nil); break;

     case T_FLOAT:  *(float*)dest = d; break;

     case T_DOUBLE: *(double*)dest = d; break;

     }

 }

+// FIXME sign with mpint

 #define CONV_TO_INTTYPE(name, ctype)                    \

 ctype conv_to_##name(void *data, numerictype_t tag)     \

 {                                                       \

@@ -54,6 +76,7 @@

     case T_UINT32: return *(uint32_t*)data;             \

     case T_INT64:  return *(int64_t*)data;              \

     case T_UINT64: return *(uint64_t*)data;             \

+    case T_MPINT:  return mptov(*(mpint**)data);        \

     case T_FLOAT:  return *(float*)data;                \

     case T_DOUBLE: return *(double*)data;               \

     }                                                   \

@@ -79,6 +102,7 @@

     case T_UINT32: return *(uint32_t*)data; break;

     case T_INT64:  return *(int64_t*)data; break;

     case T_UINT64: return *(uint64_t*)data; break;

+    case T_MPINT:  return mptouv(*(mpint**)data); break;

     case T_FLOAT:

         if (*(float*)data >= 0)

             return *(float*)data;

@@ -104,6 +128,7 @@

     case T_UINT32: return *(uint32_t*)a < *(uint32_t*)b;

     case T_INT64:  return *(int64_t*)a < *(int64_t*)b;

     case T_UINT64: return *(uint64_t*)a < *(uint64_t*)b;

+    case T_MPINT:  return mpcmp(*(mpint**)a, *(mpint**)b) < 0;

     case T_FLOAT:  return *(float*)a < *(float*)b;

     case T_DOUBLE: return *(double*)a < *(double*)b;

     }

@@ -121,6 +146,7 @@

     case T_UINT32: return *(uint32_t*)a == *(uint32_t*)b;

     case T_INT64:  return *(int64_t*)a == *(int64_t*)b;

     case T_UINT64: return *(uint64_t*)a == *(uint64_t*)b;

+    case T_MPINT:  return mpcmp(*(mpint**)a, *(mpint**)b) == 0;

     case T_FLOAT:  return *(float*)a == *(float*)b;

     case T_DOUBLE: return *(double*)a == *(double*)b;

     }

@@ -127,6 +153,9 @@

     return 0;

 }

+/* FIXME one is allocated for all compare ops */

+static mpint *cmpmpint;

+

 int cmp_lt(void *a, numerictype_t atag, void *b, numerictype_t btag)

 {

     if (atag==btag)

@@ -142,6 +171,9 @@

     if (db < da)

         return 0;

+    if (cmpmpint == nil && (atag == T_MPINT || btag == T_MPINT))

+        cmpmpint = mpnew(0);

+

     if (atag == T_UINT64) {

         if (btag == T_INT64) {

             if (*(int64_t*)b >= 0)

@@ -152,6 +184,9 @@

             if (db != db) return 0;

             return (*(uint64_t*)a < (uint64_t)*(double*)b);

         }

+        else if (btag == T_MPINT) {

+            return mpcmp(uvtomp(*(uint64_t*)a, cmpmpint), *(mpint**)b) < 0;

+        }

     }

     else if (atag == T_INT64) {

         if (btag == T_UINT64) {

@@ -163,6 +198,9 @@

             if (db != db) return 0;

             return (*(int64_t*)a < (int64_t)*(double*)b);

         }

+        else if (btag == T_MPINT) {

+            return mpcmp(vtomp(*(int64_t*)a, cmpmpint), *(mpint**)b) < 0;

+        }

     }

     if (btag == T_UINT64) {

         if (atag == T_DOUBLE) {

@@ -169,6 +207,9 @@

             if (da != da) return 0;

             return (*(uint64_t*)b > (uint64_t)*(double*)a);

         }

+        else if (atag == T_MPINT) {

+            return mpcmp(*(mpint**)a, uvtomp(*(uint64_t*)b, cmpmpint)) < 0;

+        }

     }

     else if (btag == T_INT64) {

         if (atag == T_DOUBLE) {

@@ -175,6 +216,9 @@

             if (da != da) return 0;

             return (*(int64_t*)b > (int64_t)*(double*)a);

         }

+        else if (atag == T_MPINT) {

+            return mpcmp(*(mpint**)a, vtomp(*(int64_t*)b, cmpmpint)) < 0;

+        }

     }

     return 0;

 }

@@ -200,6 +244,9 @@

     if (da != db)

         return 0;

+    if (cmpmpint == nil && (atag == T_MPINT || btag == T_MPINT))

+        cmpmpint = mpnew(0);

+

     if (atag == T_UINT64) {

         // this is safe because if a had been bigger than INT64_MAX,

         // we would already have concluded that it's bigger than b.

@@ -207,6 +254,8 @@

             return ((int64_t)*(uint64_t*)a == *(int64_t*)b);

         else if (btag == T_DOUBLE)

             return (*(uint64_t*)a == (uint64_t)(int64_t)*(double*)b);

+        else if (btag == T_MPINT)

+            return mpcmp(uvtomp(*(uint64_t*)a, cmpmpint), *(mpint**)b) == 0;

     }

     else if (atag == T_INT64) {

         if (btag == T_UINT64)

@@ -213,6 +262,8 @@

             return (*(int64_t*)a == (int64_t)*(uint64_t*)b);

         else if (btag == T_DOUBLE)

             return (*(int64_t*)a == (int64_t)*(double*)b);

+        else if (btag == T_MPINT)

+            return mpcmp(vtomp(*(int64_t*)a, cmpmpint), *(mpint**)b) == 0;

     }

     else if (btag == T_UINT64) {

         if (atag == T_INT64)

@@ -219,6 +270,8 @@

             return ((int64_t)*(uint64_t*)b == *(int64_t*)a);

         else if (atag == T_DOUBLE)

             return (*(uint64_t*)b == (uint64_t)(int64_t)*(double*)a);

+        else if (atag == T_MPINT)

+            return mpcmp(*(mpint**)a, uvtomp(*(uint64_t*)b, cmpmpint)) == 0;

     }

     else if (btag == T_INT64) {

         if (atag == T_UINT64)

@@ -225,6 +278,8 @@

             return (*(int64_t*)b == (int64_t)*(uint64_t*)a);

         else if (atag == T_DOUBLE)

             return (*(int64_t*)b == (int64_t)*(double*)a);

+        else if (atag == T_MPINT)

+            return mpcmp(*(mpint**)a, vtomp(*(int64_t*)b, cmpmpint)) == 0;

     }

     return 1;

 }

--- a/plan9/platform.h

+++ b/plan9/platform.h

@@ -1,6 +1,7 @@

 #include <u.h>

 #include <libc.h>

 #include <ctype.h>

+#include <mp.h>

 #if defined(__amd64__) || \

     defined(__arm64__) || \

--- /dev/null

+++ b/posix/mp.h

@@ -1,0 +1,231 @@

+#ifndef _MPINT_H_

+#define _MPINT_H_

+

+typedef uint32_t mpdigit;

+typedef uint8_t uchar;

+typedef uint32_t ulong;

+typedef uint32_t u32int;

+typedef unsigned int uint;

+typedef int64_t vlong;

+typedef uint64_t uvlong;

+typedef uint64_t u64int;

+

+typedef union FPdbleword FPdbleword;

+union FPdbleword

+{

+	double	x;

+	struct {	/* little endian */

+		uint lo;

+		uint hi;

+	};

+};

+

+#define sysfatal(s) do{ fprintf(stderr, "%s\n", s); abort(); }while(0)

+

+#define mpdighi  (mpdigit)(1U<<(Dbits-1))

+#define DIGITS(x) ((Dbits - 1 + (x))/Dbits)

+

+// for converting between int's and mpint's

+#define MAXUINT ((uint)-1)

+#define MAXINT (MAXUINT>>1)

+#define MININT (MAXINT+1)

+

+// for converting between vlongs's and mpint's

+#define MAXUVLONG (~0ULL)

+#define MAXVLONG (MAXUVLONG>>1)

+#define MINVLONG (MAXVLONG+1ULL)

+

+extern	int	dec64(uchar*, int, char*, int);

+extern	int	enc64(char*, int, uchar*, int);

+extern	int	dec64x(uchar*, int, char*, int, int (*)(int));

+extern	int	enc64x(char*, int, uchar*, int, int (*)(int));

+extern	int	dec32(uchar*, int, char*, int);

+extern	int	enc32(char*, int, uchar*, int);

+extern	int	dec32x(uchar*, int, char*, int, int (*)(int));

+extern	int	enc32x(char*, int, uchar*, int, int (*)(int));

+extern	int	dec16(uchar*, int, char*, int);

+extern	int	enc16(char*, int, uchar*, int);

+extern	int	dec64chr(int);

+extern	int	enc64chr(int);

+extern	int	dec32chr(int);

+extern	int	enc32chr(int);

+extern	int	dec16chr(int);

+extern	int	enc16chr(int);

+

+/*

+ * the code assumes mpdigit to be at least an int

+ * mpdigit must be an atomic type.  mpdigit is defined

+ * in the architecture specific u.h

+ */

+typedef struct mpint mpint;

+

+struct mpint

+{

+	int	sign;	/* +1 or -1 */

+	int	size;	/* allocated digits */

+	int	top;	/* significant digits */

+	mpdigit	*p;

+	char	flags;

+};

+

+enum

+{

+	MPstatic=	0x01,	/* static constant */

+	MPnorm=		0x02,	/* normalization status */

+	MPtimesafe=	0x04,	/* request time invariant computation */

+	MPfield=	0x08,	/* this mpint is a field modulus */

+

+	Dbytes=		sizeof(mpdigit),	/* bytes per digit */

+	Dbits=		Dbytes*8		/* bits per digit */

+};

+

+/* allocation */

+void	mpsetminbits(int n);	/* newly created mpint's get at least n bits */

+mpint*	mpnew(int n);		/* create a new mpint with at least n bits */

+void	mpfree(mpint *b);

+void	mpbits(mpint *b, int n);	/* ensure that b has at least n bits */

+mpint*	mpnorm(mpint *b);		/* dump leading zeros */

+mpint*	mpcopy(mpint *b);

+void	mpassign(mpint *old, mpint *new);

+

+/* random bits */

+mpint*	mprand(int bits, void (*gen)(uchar*, int), mpint *b);

+/* return uniform random [0..n-1] */

+mpint*	mpnrand(mpint *n, void (*gen)(uchar*, int), mpint *b);

+

+/* conversion */

+mpint*	strtomp(char*, char**, int, mpint*);	/* ascii */

+char*	mptoa(mpint*, int, char*, int);

+mpint*	letomp(uchar*, uint, mpint*);	/* byte array, little-endian */

+int	mptole(mpint*, uchar*, uint, uchar**);

+void	mptolel(mpint *b, uchar *p, int n);

+mpint*	betomp(uchar*, uint, mpint*);	/* byte array, big-endian */

+int	mptobe(mpint*, uchar*, uint, uchar**);

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

+uint	mptoui(mpint*);			/* unsigned int */

+mpint*	uitomp(uint, mpint*);

+int	mptoi(mpint*);			/* int */

+mpint*	itomp(int, mpint*);

+uvlong	mptouv(mpint*);			/* unsigned vlong */

+mpint*	uvtomp(uvlong, mpint*);

+vlong	mptov(mpint*);			/* vlong */

+mpint*	vtomp(vlong, mpint*);

+double	mptod(mpint*);			/* double */

+mpint*	dtomp(double, mpint*);

+

+/* divide 2 digits by one */

+void	mpdigdiv(mpdigit *dividend, mpdigit divisor, mpdigit *quotient);

+

+/* in the following, the result mpint may be */

+/* the same as one of the inputs. */

+void	mpadd(mpint *b1, mpint *b2, mpint *sum);	/* sum = b1+b2 */

+void	mpsub(mpint *b1, mpint *b2, mpint *diff);	/* diff = b1-b2 */

+void	mpleft(mpint *b, int shift, mpint *res);	/* res = b<<shift */

+void	mpright(mpint *b, int shift, mpint *res);	/* res = b>>shift */

+void	mpmul(mpint *b1, mpint *b2, mpint *prod);	/* prod = b1*b2 */

+void	mpexp(mpint *b, mpint *e, mpint *m, mpint *res);	/* res = b**e mod m */

+void	mpmod(mpint *b, mpint *m, mpint *remainder);	/* remainder = b mod m */

+

+/* logical operations */

+void	mpand(mpint *b1, mpint *b2, mpint *res);

+void	mpbic(mpint *b1, mpint *b2, mpint *res);

+void	mpor(mpint *b1, mpint *b2, mpint *res);

+void	mpnot(mpint *b, mpint *res);

+void	mpxor(mpint *b1, mpint *b2, mpint *res);

+void	mptrunc(mpint *b, int n, mpint *res);

+void	mpxtend(mpint *b, int n, mpint *res);

+void	mpasr(mpint *b, int shift, mpint *res);

+

+/* modular arithmetic, time invariant when 0≤b1≤m-1 and 0≤b2≤m-1 */

+void	mpmodadd(mpint *b1, mpint *b2, mpint *m, mpint *sum);	/* sum = b1+b2 % m */

+void	mpmodsub(mpint *b1, mpint *b2, mpint *m, mpint *diff);	/* diff = b1-b2 % m */

+void	mpmodmul(mpint *b1, mpint *b2, mpint *m, mpint *prod);	/* prod = b1*b2 % m */

+

+/* quotient = dividend/divisor, remainder = dividend % divisor */

+void	mpdiv(mpint *dividend, mpint *divisor,  mpint *quotient, mpint *remainder);

+

+/* return neg, 0, pos as b1-b2 is neg, 0, pos */

+int	mpcmp(mpint *b1, mpint *b2);

+

+/* res = s != 0 ? b1 : b2 */

+void	mpsel(int s, mpint *b1, mpint *b2, mpint *res);

+

+/* extended gcd return d, x, and y, s.t. d = gcd(a,b) and ax+by = d */

+void	mpextendedgcd(mpint *a, mpint *b, mpint *d, mpint *x, mpint *y);

+

+/* res = b**-1 mod m */

+void	mpinvert(mpint *b, mpint *m, mpint *res);

+

+/* bit counting */

+int	mpsignif(mpint*);	/* number of sigificant bits in mantissa */

+int	mplowbits0(mpint*);	/* k, where n = 2**k * q for odd q */

+

+/* well known constants */

+extern mpint	*mpzero, *mpone, *mptwo;

+

+/* sum[0:alen] = a[0:alen-1] + b[0:blen-1] */

+/* prereq: alen >= blen, sum has room for alen+1 digits */

+void	mpvecadd(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *sum);

+

+/* diff[0:alen-1] = a[0:alen-1] - b[0:blen-1] */

+/* prereq: alen >= blen, diff has room for alen digits */

+void	mpvecsub(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *diff);

+

+/* p[0:n] += m * b[0:n-1] */

+/* prereq: p has room for n+1 digits */

+void	mpvecdigmuladd(mpdigit *b, int n, mpdigit m, mpdigit *p);

+

+/* p[0:n] -= m * b[0:n-1] */

+/* prereq: p has room for n+1 digits */

+int	mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p);

+

+/* p[0:alen+blen-1] = a[0:alen-1] * b[0:blen-1] */

+/* prereq: alen >= blen, p has room for m*n digits */

+void	mpvecmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p);

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

+

+/* sign of a - b or zero if the same */

+int	mpveccmp(mpdigit *a, int alen, mpdigit *b, int blen);

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

+

+/* divide the 2 digit dividend by the one digit divisor and stick in quotient */

+/* we assume that the result is one digit - overflow is all 1's */

+void	mpdigdiv(mpdigit *dividend, mpdigit divisor, mpdigit *quotient);

+

+/* playing with magnitudes */

+int	mpmagcmp(mpint *b1, mpint *b2);

+void	mpmagadd(mpint *b1, mpint *b2, mpint *sum);	/* sum = b1+b2 */

+void	mpmagsub(mpint *b1, mpint *b2, mpint *sum);	/* sum = b1+b2 */

+

+/* chinese remainder theorem */

+typedef struct CRTpre	CRTpre;		/* precomputed values for converting */

+					/*  twixt residues and mpint */

+typedef struct CRTres	CRTres;		/* residue form of an mpint */

+

+struct CRTres

+{

+	int	n;		/* number of residues */

+	mpint	*r[1];		/* residues */

+};

+

+CRTpre*	crtpre(int, mpint**);			/* precompute conversion values */

+CRTres*	crtin(CRTpre*, mpint*);			/* convert mpint to residues */

+void	crtout(CRTpre*, CRTres*, mpint*);	/* convert residues to mpint */

+void	crtprefree(CRTpre*);

+void	crtresfree(CRTres*);

+

+/* fast field arithmetic */

+typedef struct Mfield	Mfield;

+

+struct Mfield

+{

+	mpint	m;

+	int	(*reduce)(Mfield*, mpint*, mpint*);

+};

+

+mpint *mpfield(mpint*);

+

+Mfield *gmfield(mpint*);

+Mfield *cnfield(mpint*);

+

+#endif

--- a/posix/platform.h

+++ b/posix/platform.h

@@ -23,6 +23,7 @@

 #include <unistd.h>

 #include <wctype.h>

 #include <wchar.h>

+#include "mp.h"

 #ifndef __SIZEOF_POINTER__

 #error pointer size unknown

--- a/print.c

+++ b/print.c

@@ -700,6 +700,15 @@

         else

             HPOS += ios_printf(f, "#%s(%"PRIu64")", symbol_name(type), ui64);

     }

+    else if (type == mpintsym) {

+        mpint *i = *(mpint**)data;

+        char *s = mptoa(i, 10, nil, 0);

+        if (weak || print_princ)

+            HPOS += ios_printf(f, "%s", s);

+        else

+            HPOS += ios_printf(f, "#%s(%s)", symbol_name(type), s);

+        free(s);

+    }

     else if (issymbol(type)) {

         // handle other integer prims. we know it's smaller than uint64

         // at this point, so int64 is big enough to capture everything.

--- a/read.c

+++ b/read.c

@@ -6,23 +6,25 @@

 };

 #if defined(__plan9__)

-// FIXME remove all this after adding mp entirely

-#include <mp.h>

-#define ERANGE 34

+static int errno;

 #define VLONG_MAX ~(1LL<<63)

 #define VLONG_MIN (1LL<<63)

 #define UVLONG_MAX (1LL<<63)

-static int errno;

 static mpint *mp_vlong_min, *mp_vlong_max, *mp_uvlong_max;

+#endif

-static vlong strtoll_(char *nptr, char **rptr, int base)

+static int64_t strtoll_mp(char *nptr, char **rptr, int base, mpint **mp)

 {

-    vlong x;

-    mpint *m, *c;

+    int64_t x;

+    mpint *m;

+    *mp = nil;

+    errno = 0;

     x = strtoll(nptr, rptr, base);

+#if defined(__plan9__)

     if((x != VLONG_MAX && x != VLONG_MIN) || *rptr == nptr)

         return x;

+    mpint *c;

     m = strtomp(nptr, rptr, base, nil);

     if(x == VLONG_MAX){

         if(mp_vlong_max == nil) mp_vlong_max = vtomp(VLONG_MAX, nil);

@@ -31,29 +33,45 @@

         if(mp_vlong_min == nil) mp_vlong_min = vtomp(VLONG_MIN, nil);

         c = mp_vlong_min;

     }

-    errno = mpcmp(c, m) == 0 ? 0 : ERANGE;

-    mpfree(m);

+    if (mpcmp(c, m) == 0) {

+        mpfree(m);

+        m = nil;

+    }

+#else

+    m = nil;

+    if (errno == ERANGE && (x == LLONG_MAX || x == LLONG_MIN))

+        m = strtomp(nptr, rptr, base, nil);

+#endif

+    *mp = m;

     return x;

 }

-static uvlong strtoull_(char *nptr, char **rptr, int base)

+static uint64_t strtoull_mp(char *nptr, char **rptr, int base, mpint **mp)

 {

-    uvlong x;

+    uint64_t x;

     mpint *m;

+    *mp = nil;

+    errno = 0;

     x = strtoull(nptr, rptr, base);

+#if defined(__plan9__)

     if(x != UVLONG_MAX || *rptr == nptr)

         return x;

     m = strtomp(nptr, rptr, base, nil);

     if(mp_uvlong_max == nil)

         mp_uvlong_max = uvtomp(UVLONG_MAX, nil);

-    errno = mpcmp(mp_uvlong_max, m) == 0 ? 0 : ERANGE;

-    mpfree(m);

+    if(mpcmp(mp_uvlong_max, m) == 0){

+        mpfree(m);

+        m = nil;

+    }

+#else

+    m = nil;

+    if (errno == ERANGE && x == ULLONG_MAX)

+        m = strtomp(nptr, rptr, base, nil);

+#endif

+    *mp = m;

     return x;

 }

-#define strtoll strtoll_

-#define strtoull strtoull_

-#endif

 #define F value2c(ios_t*,readstate->source)

@@ -73,6 +91,7 @@

     int64_t i64;

     uint64_t ui64;

     double d;

+    mpint *mp = nil;

     if (*tok == '\0')

         return 0;

     if (!((tok[0]=='0' && tok[1]=='x') || (base >= 15)) &&

@@ -110,18 +129,14 @@

             if (pval) *pval = mk_double(D_NINF);

             return 1;

         }

-        errno = 0;

-        i64 = strtoll(tok, &end, base);

-        if (errno)

-            return 0;

-        if (pval) *pval = return_from_int64(i64);

+        i64 = strtoll_mp(tok, &end, base, &mp);

+        if (pval)

+            *pval = mp == nil ? return_from_int64(i64) : mk_mpint(mp);

         return (*end == '\0');

     }

-    errno = 0;

-    ui64 = strtoull(tok, &end, base);

-    if (errno)

-        return 0;

-    if (pval) *pval = return_from_uint64(ui64);

+    ui64 = strtoull_mp(tok, &end, base, &mp);

+    if (pval)

+        *pval = mp == nil ? return_from_uint64(ui64) : mk_mpint(mp);

     return (*end == '\0');

 }

@@ -132,12 +147,7 @@

 static int read_numtok(char *tok, value_t *pval, int base)

 {

-    int result;

-    errno = 0;

-    result = isnumtok_base(tok, pval, base);

-    if (errno == ERANGE)

-        lerrorf(ParseError, "read: overflow in numeric constant %s", tok);

-    return result;

+    return isnumtok_base(tok, pval, base);

 }

 static uint32_t toktype = TOK_NONE;

@@ -182,7 +192,7 @@

 {

     buf[(*pi)++] = c;

     if (*pi >= (int)(sizeof(buf)-1))

-        lerrorf(ParseError, "read: token too long");

+        lerrorf(ParseError, "token too long");

 }

 // return: 1 if escaped (forced to be symbol)

@@ -262,7 +272,7 @@

     else if (c == '#') {

         ch = ios_getc(F); c = (char)ch;

         if (ch == IOS_EOF)

-            lerrorf(ParseError, "read: invalid read macro");

+            lerrorf(ParseError, "invalid read macro");

         if (c == '.') {

             toktype = TOK_SHARPDOT;

         }

@@ -272,7 +282,7 @@

         else if (c == '\\') {

             uint32_t cval;

             if (ios_getutf8(F, &cval) == IOS_EOF)

-                lerrorf(ParseError, "read: end of input in character constant");

+                lerrorf(ParseError, "end of input in character constant");

             if (cval == (uint32_t)'u' || cval == (uint32_t)'U' ||

                 cval == (uint32_t)'x') {

                 read_token('u', 0);

@@ -300,7 +310,7 @@

                 else if (tokval == spacesym)      cval = 0x20;

                 else if (tokval == deletesym)     cval = 0x7F;

                 else

-                    lerrorf(ParseError, "read: unknown character #\\%s", buf);

+                    lerrorf(ParseError, "unknown character #\\%s", buf);

             }

             toktype = TOK_NUM;

             tokval = mk_wchar(cval);

@@ -309,7 +319,7 @@

             toktype = TOK_SHARPOPEN;

         }

         else if (c == '<') {

-            lerrorf(ParseError, "read: unreadable object");

+            lerrorf(ParseError, "unreadable object");

         }

         else if (isdigit(c)) {

             read_token(c, 1);

@@ -319,11 +329,10 @@

             else if (c == '=')

                 toktype = TOK_LABEL;

             else

-                lerrorf(ParseError, "read: invalid label");

-            errno = 0;

+                lerrorf(ParseError, "invalid label");

             x = strtoll(buf, &end, 10);

-            if (*end != '\0' || errno)

-                lerrorf(ParseError, "read: invalid label");

+            if (*end != '\0')

+                lerrorf(ParseError, "invalid label");

             tokval = fixnum(x);

         }

         else if (c == '!') {

@@ -340,7 +349,7 @@

                 ch = ios_getc(F);

             hashpipe_gotc:

                 if (ch == IOS_EOF)

-                    lerrorf(ParseError, "read: eof within comment");

+                    lerrorf(ParseError, "eof within comment");

                 if ((char)ch == '|') {

                     ch = ios_getc(F);

                     if ((char)ch == '#') {

@@ -374,10 +383,9 @@

             if ((char)ch == 'g')

                 ch = ios_getc(F);

             read_token((char)ch, 0);

-            errno = 0;

             x = strtol(buf, &end, 10);

-            if (*end != '\0' || buf[0] == '\0' || errno)

-                lerrorf(ParseError, "read: invalid gensym label");

+            if (*end != '\0' || buf[0] == '\0')

+                lerrorf(ParseError, "invalid gensym label");

             toktype = TOK_GENSYM;

             tokval = fixnum(x);

         }

@@ -391,7 +399,7 @@

                 (isdigit_base(buf[1],base) ||

                  buf[1]=='-')) {

                 if (!read_numtok(&buf[1], &tokval, base))

-                    lerrorf(ParseError, "read: invalid base %d constant", base);

+                    lerrorf(ParseError, "invalid base %d constant", base);

                 return (toktype=TOK_NUM);

             }

@@ -399,7 +407,7 @@

             tokval = symbol(buf);

         }

         else {

-            lerrorf(ParseError, "read: unknown read macro");

+            lerrorf(ParseError, "unknown read macro");

         }

     }

     else if (c == ',') {

@@ -465,7 +473,7 @@

         ptrhash_put(&readstate->backrefs, (void*)label, (void*)v);

     while (peek() != closer) {

         if (ios_eof(F))

-            lerrorf(ParseError, "read: unexpected end of input");

+            lerrorf(ParseError, "unexpected end of input");

         if (i >= vector_size(v)) {

             v = Stack[SP-1] = vector_grow(v);

             if (label != UNBOUND)

@@ -499,7 +507,7 @@

             temp = realloc(buf, sz);

             if (temp == nil) {

                 free(buf);

-                lerrorf(ParseError, "read: out of memory reading string");

+                lerrorf(ParseError, "out of memory reading string");

             }

             buf = temp;

         }

@@ -506,7 +514,7 @@

         c = ios_getc(F);

         if (c == IOS_EOF) {

             free(buf);

-            lerrorf(ParseError, "read: unexpected end of input in string");

+            lerrorf(ParseError, "unexpected end of input in string");

         }

         if (c == '"')

             break;

@@ -514,7 +522,7 @@

             c = ios_getc(F);

             if (c == IOS_EOF) {

                 free(buf);

-                lerrorf(ParseError, "read: end of input in escape sequence");

+                lerrorf(ParseError, "end of input in escape sequence");

             }

             j = 0;

             if (octal_digit(c)) {

@@ -544,7 +552,7 @@

                 if (j) wc = strtol(eseq, nil, 16);

                 if (!j || wc > 0x10ffff) {

                     free(buf);

-                    lerrorf(ParseError, "read: invalid escape sequence");

+                    lerrorf(ParseError, "invalid escape sequence");

                 }

                 if (ndig == 2)

                     buf[i++] = ((char)wc);

@@ -555,7 +563,7 @@

                 char esc = read_escape_control_char((char)c);

                 if (esc == (char)c && !strchr("\\'\"`", esc)) {

                     free(buf);

-                    lerrorf(ParseError, "read: invalid escape sequence: \\%c", (char)c);

+                    lerrorf(ParseError, "invalid escape sequence: \\%c", (char)c);

                 }

                 buf[i++] = esc;

             }

@@ -583,7 +591,7 @@

     t = peek();

     while (t != closer) {

         if (ios_eof(F))

-            lerrorf(ParseError, "read: unexpected end of input");

+            lerrorf(ParseError, "unexpected end of input");

         c = mk_cons(); car_(c) = cdr_(c) = NIL;

         if (iscons(*pc)) {

             cdr_(*pc) = c;

@@ -604,10 +612,10 @@

             cdr_(*pc) = c;

             t = peek();

             if (ios_eof(F))

-                lerrorf(ParseError, "read: unexpected end of input");

+                lerrorf(ParseError, "unexpected end of input");

             if (t != closer) {

                 take();

-                lerrorf(ParseError, "read: expected '%c'", closer==TOK_CLOSEB ? ']' : ')');

+                lerrorf(ParseError, "expected '%c'", closer==TOK_CLOSEB ? ']' : ')');

             }

         }

     }

@@ -628,11 +636,11 @@

     take();

     switch (t) {

     case TOK_CLOSE:

-        lerrorf(ParseError, "read: unexpected ')'");

+        lerrorf(ParseError, "unexpected ')'");

     case TOK_CLOSEB:

-        lerrorf(ParseError, "read: unexpected ']'");

+        lerrorf(ParseError, "unexpected ']'");

     case TOK_DOT:

-        lerrorf(ParseError, "read: unexpected '.'");

+        lerrorf(ParseError, "unexpected '.'");

     case TOK_SYM:

     case TOK_NUM:

         return tokval;

@@ -678,7 +686,7 @@

         c = nextchar();

         if (c != '(') {

             take();

-            lerrorf(ParseError, "read: expected argument list for %s",

+            lerrorf(ParseError, "expected argument list for %s",

                     symbol_name(tokval));

         }

         PUSH(NIL);

@@ -713,7 +721,7 @@

     case TOK_LABEL:

         // create backreference label

         if (ptrhash_has(&readstate->backrefs, (void*)tokval))

-            lerrorf(ParseError, "read: label %"PRIdPTR" redefined", numval(tokval));

+            lerrorf(ParseError, "label %"PRIdPTR" redefined", numval(tokval));

         oldtokval = tokval;

         v = do_read_sexpr(tokval);

         ptrhash_put(&readstate->backrefs, (void*)oldtokval, (void*)v);

@@ -722,7 +730,7 @@

         // look up backreference

         v = (value_t)ptrhash_get(&readstate->backrefs, (void*)tokval);

         if (v == (value_t)HT_NOTFOUND)

-            lerrorf(ParseError, "read: undefined label %"PRIdPTR, numval(tokval));

+            lerrorf(ParseError, "undefined label %"PRIdPTR, numval(tokval));

         return v;

     case TOK_GENSYM:

         pv = (value_t*)ptrhash_bp(&readstate->gensyms, (void*)tokval);

--- a/test/unittest.lsp

+++ b/test/unittest.lsp

@@ -82,6 +82,9 @@

 (assert (> 9223372036854775808 9223372036854775807))

+; mpint

+(assert (> 0x10000000000000000 0x8fffffffffffffff))

+

 ; NaNs

 (assert (equal? +nan.0 +nan.0))

 (assert (not (= +nan.0 +nan.0)))