shithub: femtolisp

Download patch

ref: 8837293e1085d71e10ef9216c457458450419801
parent: bdb4b15192e2452ec6019675517599b970755c7c
author: Sigrid Solveig Haflínudóttir <sigrid@ftrv.se>
date: Fri Mar 31 19:35:47 EDT 2023

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