shithub: femtolisp

Info  •  Files  •  Log

div0 and bitwise logic for bigints

--- a/Makefile

+++ b/Makefile

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

 TARG=flisp

 LLT=llt/libllt.a

 CFLAGS?=-O2 -g

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

+CFLAGS+=-Wall -Wextra -std=c99 -D_DEFAULT_SOURCE

 LDFLAGS?=

 OBJS=\

--- a/README.md

+++ b/README.md

@@ -19,7 +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)

+ * bigints

 ## Characteristics

--- a/cvalues.c

+++ b/cvalues.c

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

     numerictype_t ta, tb;

     void *aptr, *bptr;

     int64_t a64, b64;

+    mpint *x;

     if (!num_to_ptr(a, &ai, &ta, &aptr))

         type_error("number", a);

@@ -1285,6 +1286,22 @@

     if (!num_to_ptr(b, &bi, &tb, &bptr))

         type_error("number", b);

+    if (ta == T_MPINT) {

+        if (tb == T_MPINT) {

+            if (mpsignif(*(mpint**)bptr) == 0)

+                goto div_error;

+            x = mpnew(0);

+            mpdiv(*(mpint**)aptr, *(mpint**)bptr, x, nil);

+            return mk_mpint(x);

+        } else {

+            b64 = conv_to_int64(bptr, tb);

+            if (b64 == 0)

+                goto div_error;

+            x = tb == T_UINT64 ? uvtomp(b64, nil) : vtomp(b64, nil);

+            mpdiv(*(mpint**)aptr, x, x, nil);

+            return mk_mpint(x);

+        }

+    }

     if (ta == T_UINT64) {

         if (tb == T_UINT64) {

             if (*(uint64_t*)bptr == 0) goto div_error;

@@ -1322,6 +1339,7 @@

     lltint_t ai, bi;

     numerictype_t ta, tb, itmp;

     void *aptr=nil, *bptr=nil, *ptmp;

+    mpint *bmp, *resmp;

     int64_t b64;

     if (!num_to_ptr(a, &ai, &ta, &aptr) || ta >= T_FLOAT)

@@ -1334,7 +1352,17 @@

         ptmp = aptr; aptr = bptr; bptr = ptmp;

     }

     // now a's type is larger than or same as b's

-    b64 = conv_to_int64(bptr, tb);

+    if (ta == T_MPINT) {

+        if (tb == T_MPINT) {

+            bmp = *(mpint**)bptr;

+            resmp = mpnew(0);

+        } else {

+            bmp = conv_to_mpint(bptr, tb);

+            resmp = bmp;

+        }

+    }

+    else

+        b64 = conv_to_int64(bptr, tb);

     switch (opcode) {

     case 0:

     switch (ta) {

@@ -1346,6 +1374,7 @@

     case T_UINT32: return mk_uint32(*(uint32_t*)aptr & (uint32_t)b64);

     case T_INT64:  return mk_int64( *(int64_t*)aptr  & (int64_t )b64);

     case T_UINT64: return mk_uint64(*(uint64_t*)aptr & (uint64_t)b64);

+    case T_MPINT:  mpand(*(mpint**)aptr, bmp, resmp); return mk_mpint(resmp);

     case T_FLOAT:

     case T_DOUBLE: assert(0);

     }

@@ -1360,6 +1389,7 @@

     case T_UINT32: return mk_uint32(*(uint32_t*)aptr | (uint32_t)b64);

     case T_INT64:  return mk_int64( *(int64_t*)aptr  | (int64_t )b64);

     case T_UINT64: return mk_uint64(*(uint64_t*)aptr | (uint64_t)b64);

+    case T_MPINT:  mpor(*(mpint**)aptr, bmp, resmp); return mk_mpint(resmp);

     case T_FLOAT:

     case T_DOUBLE: assert(0);

     }

@@ -1374,6 +1404,7 @@

     case T_UINT32: return mk_uint32(*(uint32_t*)aptr ^ (uint32_t)b64);

     case T_INT64:  return mk_int64( *(int64_t*)aptr  ^ (int64_t )b64);

     case T_UINT64: return mk_uint64(*(uint64_t*)aptr ^ (uint64_t)b64);

+    case T_MPINT:  mpxor(*(mpint**)aptr, bmp, resmp); return mk_mpint(resmp);

     case T_FLOAT:

     case T_DOUBLE: assert(0);

     }

--- a/llt/Makefile

+++ b/llt/Makefile

@@ -17,23 +17,24 @@

     \

 	mpadd.o\

 	mpaux.o\

-	mpvecsub.o\

 	mpcmp.o\

+	mpdigdiv.o\

 	mpdiv.o\

 	mpfmt.o\

-	mpmul.o\

-	mpsub.o\

 	mpleft.o\

+	mplogic.o\

+	mpmul.o\

 	mpright.o\

+	mpsub.o\

 	mptobe.o\

 	mptod.o\

 	mptoi.o\

 	mptoui.o\

 	mptouv.o\

-	mpdigdiv.o\

-	mpvecdigmuladd.o\

 	mptov.o\

 	mpvecadd.o\

+	mpvecdigmuladd.o\

+	mpvecsub.o\

 	strtomp.o\

 	u16.o\

 	u32.o\

--- /dev/null

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

+}