ref: 467826415fe2cc2846097dfd37ec879b9d62cf62
parent: 076b3b75c0f116421c17e7c2afd52b3079521ef4
author: Sigrid Solveig Haflínudóttir <[email protected]>
date: Sat Nov 16 20:52:20 EST 2024
mp: import a test and fix the issues it found
--- a/3rd/mp/mpaux.c
+++ b/3rd/mp/mpaux.c
@@ -33,8 +33,10 @@
void
mpsetminbits(int n)
{- if(n < 0)
- sysfatal("mpsetminbits: n < 0");+ if(n < 0){+ fprintf(stderr, "mpsetminbits: n < 0\n");
+ exit(2);
+ }
if(n == 0)
n = 1;
mpmindigits = DIGITS(n);
@@ -46,15 +48,19 @@
{mpint *b;
- if(n < 0)
- sysfatal("mpnew: n < 0");+ if(n < 0){+ fprintf(stderr, "mpnew: n < 0\n");
+ exit(2);
+ }
n = DIGITS(n);
if(n < mpmindigits)
n = mpmindigits;
b = calloc(1, sizeof(mpint) + n*Dbytes);
- if(b == nil)
- sysfatal("mpnew: %r");+ if(b == nil){+ fprintf(stderr, "mpnew: no memory\n");
+ exit(2);
+ }
b->p = (mpdigit*)&b[1];
b->size = n;
b->sign = 1;
@@ -76,14 +82,18 @@
} else { if(b->p == (mpdigit*)&b[1]){b->p = (mpdigit*)MEM_ALLOC(n*Dbytes);
- if(b->p == nil)
- sysfatal("mpbits: %r");+ if(b->p == nil){+ fprintf(stderr, "mpbits: no memory\n");
+ exit(2);
+ }
memmove(b->p, &b[1], Dbytes*b->top);
memset(&b[1], 0, Dbytes*b->size);
} else {b->p = (mpdigit*)MEM_REALLOC(b->p, n*Dbytes);
- if(b->p == nil)
- sysfatal("mpbits: %r");+ if(b->p == nil){+ fprintf(stderr, "mpbits: no memory\n");
+ exit(2);
+ }
}
b->size = n;
}
@@ -97,8 +107,10 @@
{if(b == nil)
return;
- if(b->flags & MPstatic)
- sysfatal("freeing mp constant");+ if(b->flags & MPstatic){+ fprintf(stderr, "freeing mp constant\n");
+ exit(2);
+ }
memset(b->p, 0, b->size*Dbytes);
if(b->p != (mpdigit*)&b[1])
MEM_FREE(b->p);
--- a/3rd/mp/mpdiv.c
+++ b/3rd/mp/mpdiv.c
@@ -14,9 +14,11 @@
assert(quotient != remainder);
assert(divisor->flags & MPnorm);
- // divide bv zero
- if(divisor->top == 0)
- abort();
+ // divide by zero
+ if(divisor->top == 0){+ fprintf(stderr, "mpdiv: by zero\n");
+ exit(2);
+ }
// division by one or small powers of two
if(divisor->top == 1 && (divisor->p[0] & divisor->p[0]-1) == 0){--- /dev/null
+++ b/3rd/mp/mpexp.c
@@ -1,0 +1,97 @@
+#include "platform.h"
+#include "mp.h"
+
+// res = b**e
+//
+// knuth, vol 2, pp 398-400
+
+enum {+ Freeb= 0x1,
+ Freee= 0x2,
+ Freem= 0x4,
+};
+
+//int expdebug;
+
+void
+mpexp(mpint *b, mpint *e, mpint *m, mpint *res)
+{+ mpint *t[2];
+ int tofree;
+ mpdigit d, bit;
+ int i, j;
+
+ assert(m == nil || m->flags & MPnorm);
+ assert((e->flags & MPtimesafe) == 0);
+ res->flags |= b->flags & MPtimesafe;
+
+ i = mpcmp(e,mpzero);
+ if(i==0){+ mpassign(mpone, res);
+ return;
+ }
+ if(i < 0){+ fprintf(stderr, "mpexp: negative exponent\n");
+ exit(2);
+ }
+
+ t[0] = mpcopy(b);
+ t[1] = res;
+
+ tofree = 0;
+ if(res == b){+ b = mpcopy(b);
+ tofree |= Freeb;
+ }
+ if(res == e){+ e = mpcopy(e);
+ tofree |= Freee;
+ }
+ if(res == m){+ m = mpcopy(m);
+ tofree |= Freem;
+ }
+
+ // skip first bit
+ i = e->top-1;
+ d = e->p[i];
+ for(bit = mpdighi; (bit & d) == 0; bit >>= 1)
+ ;
+ bit >>= 1;
+
+ j = 0;
+ for(;;){+ for(; bit != 0; bit >>= 1){+ if(m != nil)
+ mpmodmul(t[j], t[j], m, t[j^1]);
+ else
+ mpmul(t[j], t[j], t[j^1]);
+ if(bit & d) {+ if(m != nil)
+ mpmodmul(t[j^1], b, m, t[j]);
+ else
+ mpmul(t[j^1], b, t[j]);
+ } else
+ j ^= 1;
+ }
+ if(--i < 0)
+ break;
+ bit = mpdighi;
+ d = e->p[i];
+ }
+ if(t[j] == res){+ mpfree(t[j^1]);
+ } else {+ mpassign(t[j], res);
+ mpfree(t[j]);
+ }
+
+ if(tofree){+ if(tofree & Freeb)
+ mpfree(b);
+ if(tofree & Freee)
+ mpfree(e);
+ if(tofree & Freem)
+ mpfree(m);
+ }
+}
--- /dev/null
+++ b/3rd/mp/mpextendedgcd.c
@@ -1,0 +1,115 @@
+#include "platform.h"
+#include "mp.h"
+
+#define iseven(a) (((a)->p[0] & 1) == 0)
+
+// extended binary gcd
+//
+// For a and b it solves, v = gcd(a,b) and finds x and y s.t.
+// ax + by = v
+//
+// Handbook of Applied Cryptography, Menezes et al, 1997, pg 608.
+void
+mpextendedgcd(mpint *a, mpint *b, mpint *v, mpint *x, mpint *y)
+{+ mpint *u, *A, *B, *C, *D;
+ int g;
+
+ if(v == nil){+ v = mpnew(0);
+ mpextendedgcd(a, b, v, x, y);
+ mpfree(v);
+ return;
+ }
+ assert(x == nil || (x->flags & MPtimesafe) == 0);
+ assert(y == nil || (y->flags & MPtimesafe) == 0);
+ assert((a->flags&b->flags) & MPnorm);
+ assert(((a->flags|b->flags|v->flags) & MPtimesafe) == 0);
+
+ if(a->sign < 0 || b->sign < 0){+ mpassign(mpzero, v);
+ mpassign(mpzero, y);
+ mpassign(mpzero, x);
+ return;
+ }
+
+ if(a->top == 0){+ mpassign(b, v);
+ mpassign(mpone, y);
+ mpassign(mpzero, x);
+ return;
+ }
+ if(b->top == 0){+ mpassign(a, v);
+ mpassign(mpone, x);
+ mpassign(mpzero, y);
+ return;
+ }
+
+ g = 0;
+ a = mpcopy(a);
+ b = mpcopy(b);
+
+ while(iseven(a) && iseven(b)){+ mpright(a, 1, a);
+ mpright(b, 1, b);
+ g++;
+ }
+
+ u = mpcopy(a);
+ mpassign(b, v);
+ A = mpcopy(mpone);
+ B = mpcopy(mpzero);
+ C = mpcopy(mpzero);
+ D = mpcopy(mpone);
+
+ for(;;) {+// print("%B %B %B %B %B %B\n", u, v, A, B, C, D);+ while(iseven(u)){+ mpright(u, 1, u);
+ if(!iseven(A) || !iseven(B)) {+ mpadd(A, b, A);
+ mpsub(B, a, B);
+ }
+ mpright(A, 1, A);
+ mpright(B, 1, B);
+ }
+
+// print("%B %B %B %B %B %B\n", u, v, A, B, C, D);+ while(iseven(v)){+ mpright(v, 1, v);
+ if(!iseven(C) || !iseven(D)) {+ mpadd(C, b, C);
+ mpsub(D, a, D);
+ }
+ mpright(C, 1, C);
+ mpright(D, 1, D);
+ }
+
+// print("%B %B %B %B %B %B\n", u, v, A, B, C, D);+ if(mpcmp(u, v) >= 0){+ mpsub(u, v, u);
+ mpsub(A, C, A);
+ mpsub(B, D, B);
+ } else {+ mpsub(v, u, v);
+ mpsub(C, A, C);
+ mpsub(D, B, D);
+ }
+
+ if(u->top == 0)
+ break;
+
+ }
+ mpassign(C, x);
+ mpassign(D, y);
+ mpleft(v, g, v);
+
+ mpfree(A);
+ mpfree(B);
+ mpfree(C);
+ mpfree(D);
+ mpfree(u);
+ mpfree(a);
+ mpfree(b);
+}
--- a/3rd/mp/mpfmt.c
+++ b/3rd/mp/mpfmt.c
@@ -13,7 +13,7 @@
sn = 1<<s;
out = buf;
eout = buf+len;
- for(p = &b->p[b->top-1]; p >= b->p; p--){+ 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;
@@ -89,7 +89,8 @@
{mpdigit x, y;
char *out;
- uint32_t i, j;
+ uint32_t j;
+ int i;
if(len < 2)
return -1;
@@ -107,7 +108,8 @@
x = y;
i += Dbits;
while(i >= 3){-Digout: i -= 3;
+Digout:
+ i -= 3;
if(out > buf)
out--;
else if(x != 0)
@@ -174,8 +176,8 @@
rv = topow2(b, out, len, 1);
break;
default:
- abort();
- return nil;
+ fprintf(stderr, "mptoa: invalid base: %d\n", base);
+ exit(2);
}
if(rv < 0){if(alloced)
--- /dev/null
+++ b/3rd/mp/mpinvert.c
@@ -1,0 +1,19 @@
+#include "platform.h"
+#include "mp.h"
+
+// use extended gcd to find the multiplicative inverse
+// res = b**-1 mod m
+void
+mpinvert(mpint *b, mpint *m, mpint *res)
+{+ mpint *v;
+
+ v = mpnew(0);
+ mpextendedgcd(b, m, v, res, nil);
+ if(mpcmp(v, mpone) != 0){+ fprintf(stderr, "mpinvert: impossible\n");
+ exit(2);
+ }
+ mpfree(v);
+ mpmod(res, m, res);
+}
--- /dev/null
+++ b/3rd/mp/mpmod.c
@@ -1,0 +1,19 @@
+#include "platform.h"
+#include "mp.h"
+
+void
+mpmod(mpint *x, mpint *n, mpint *r)
+{+ int sign;
+ mpint *ns;
+
+ sign = x->sign;
+ ns = sign < 0 && n == r ? mpcopy(n) : n;
+ if((n->flags & MPfield) == 0
+ || ((Mfield*)n)->reduce((Mfield*)n, x, r) != 0)
+ mpdiv(x, n, nil, r);
+ if(sign < 0){+ mpmagsub(ns, r, r);
+ if(ns != n) mpfree(ns);
+ }
+}
--- /dev/null
+++ b/3rd/mp/mpmodop.c
@@ -1,0 +1,95 @@
+#include "platform.h"
+#include "mp.h"
+
+/* operands need to have m->top+1 digits of space and satisfy 0 ≤ a ≤ m-1 */
+static mpint*
+modarg(mpint *a, mpint *m)
+{+ if(a->size <= m->top || a->sign < 0 || mpmagcmp(a, m) >= 0){+ a = mpcopy(a);
+ mpmod(a, m, a);
+ mpbits(a, Dbits*(m->top+1));
+ a->top = m->top;
+ } else if(a->top < m->top){+ memset(&a->p[a->top], 0, (m->top - a->top)*Dbytes);
+ }
+ return a;
+}
+
+void
+mpmodadd(mpint *b1, mpint *b2, mpint *m, mpint *sum)
+{+ mpint *a, *b;
+ mpdigit d;
+ uint32_t i, j;
+
+ a = modarg(b1, m);
+ b = modarg(b2, m);
+
+ sum->flags |= (a->flags | b->flags) & MPtimesafe;
+ mpbits(sum, Dbits*2*(m->top+1));
+
+ mpvecadd(a->p, m->top, b->p, m->top, sum->p);
+ mpvecsub(sum->p, m->top+1, m->p, m->top, sum->p+m->top+1);
+
+ d = sum->p[2*m->top+1];
+ for(i = 0, j = m->top+1; i < m->top; i++, j++)
+ sum->p[i] = (sum->p[i] & d) | (sum->p[j] & ~d);
+
+ sum->top = m->top;
+ sum->sign = 1;
+ mpnorm(sum);
+
+ if(a != b1)
+ mpfree(a);
+ if(b != b2)
+ mpfree(b);
+}
+
+void
+mpmodsub(mpint *b1, mpint *b2, mpint *m, mpint *diff)
+{+ mpint *a, *b;
+ mpdigit d;
+ uint32_t i, j;
+
+ a = modarg(b1, m);
+ b = modarg(b2, m);
+
+ diff->flags |= (a->flags | b->flags) & MPtimesafe;
+ mpbits(diff, Dbits*2*(m->top+1));
+
+ a->p[m->top] = 0;
+ mpvecsub(a->p, m->top+1, b->p, m->top, diff->p);
+ mpvecadd(diff->p, m->top, m->p, m->top, diff->p+m->top+1);
+
+ d = ~diff->p[m->top];
+ for(i = 0, j = m->top+1; i < m->top; i++, j++)
+ diff->p[i] = (diff->p[i] & d) | (diff->p[j] & ~d);
+
+ diff->top = m->top;
+ diff->sign = 1;
+ mpnorm(diff);
+
+ if(a != b1)
+ mpfree(a);
+ if(b != b2)
+ mpfree(b);
+}
+
+void
+mpmodmul(mpint *b1, mpint *b2, mpint *m, mpint *prod)
+{+ mpint *a, *b;
+
+ a = modarg(b1, m);
+ b = modarg(b2, m);
+
+ mpmul(a, b, prod);
+ mpmod(prod, m, prod);
+
+ if(a != b1)
+ mpfree(a);
+ if(b != b2)
+ mpfree(b);
+}
--- a/3rd/mp/mpmul.c
+++ b/3rd/mp/mpmul.c
@@ -39,8 +39,10 @@
// room for the partial products
t = calloc(1, Dbytes*5*(2*n+1));
- if(t == nil)
- sysfatal("mpkaratsuba: %r");+ if(t == nil){+ fprintf(stderr, "mpkaratsuba: no memory\n");
+ exit(2);
+ }
u0v0 = t;
u1v1 = t + (2*n+1);
diffprod = t + 2*(2*n+1);
--- /dev/null
+++ b/3rd/mp/mprand.c
@@ -1,0 +1,23 @@
+#include "platform.h"
+#include "mp.h"
+
+mpint *
+mprand(int bits, void (*gen)(uint8_t*, int), mpint *b)
+{+ mpdigit mask;
+
+ if(b == nil)
+ b = mpnew(bits);
+ else
+ mpbits(b, bits);
+
+ b->sign = 1;
+ b->top = DIGITS(bits);
+ (*gen)((uint8_t*)b->p, b->top*Dbytes);
+
+ mask = ((mpdigit)1 << (bits%Dbits))-1;
+ if(mask != 0)
+ b->p[b->top-1] &= mask;
+
+ return mpnorm(b);
+}
--- a/3rd/mp/mptobe.c
+++ b/3rd/mp/mptobe.c
@@ -14,8 +14,10 @@
if(p == nil){n = m;
p = MEM_ALLOC(n);
- if(p == nil)
- sysfatal("mptobe: %r");+ if(p == nil){+ fprintf(stderr, "mptobe: no memory\n");
+ exit(2);
+ }
} else {if(n < m)
return -1;
--- a/3rd/mp/strtomp.c
+++ b/3rd/mp/strtomp.c
@@ -108,9 +108,8 @@
int sign;
const char *e;
- if(b == nil){+ if(b == nil)
b = mpnew(0);
- }
while(*a==' ' || *a=='\t')
a++;
@@ -131,10 +130,10 @@
if(a[1] == 'x' || a[1] == 'X') {a += 2;
base = 16;
- } else if(a[1] == 'b' || a[1] == 'B') {+ }else if(a[1] == 'b' || a[1] == 'B') {a += 2;
base = 2;
- } else if(a[1] >= '0' && a[1] <= '7') {+ }else if(a[1] >= '0' && a[1] <= '7') {a++;
base = 8;
}
--- /dev/null
+++ b/3rd/mp/test.c
@@ -1,0 +1,500 @@
+#include "platform.h"
+#include "mp.h"
+#include "ieee754.h"
+
+double D_PNAN, D_NNAN, D_PINF, D_NINF;
+float F_PNAN, F_NNAN, F_PINF, F_NINF;
+
+static int loops = 1;
+static char str[16][8192];
+static int istr = -1;
+
+static char *
+MPB(int base, mpint *m)
+{+ char *s, *b;
+ bool minus;
+ istr = (istr+1) % nelem(str);
+ b = str[istr];
+ s = mptoa(m, base, b+2, sizeof(str[istr])-2);
+ if(base == 10)
+ return s;
+ minus = s[0] == '-';
+ if(minus)
+ s++;
+ if(base == 2){+ *--s = 'b';
+ *--s = '0';
+ }else if(base == 8){+ *--s = '0';
+ }else if(base == 16){+ *--s = 'x';
+ *--s = '0';
+ }
+ if(minus)
+ *--s = '-';
+ return s;
+}
+
+static char *
+MP(mpint *m)
+{+ char *b;
+ istr = (istr+1) % nelem(str);
+ b = str[istr];
+ return mptoa(m, 10, b, sizeof(str[istr]));
+}
+
+static int64_t
+nsec(void)
+{+ return 0;
+}
+
+static void
+prng(uint8_t *p, int n)
+{+ while(n-- > 0)
+ *p++ = rand();
+}
+
+static void
+testconv(char *str)
+{+ int i, base[] = {2,8,10,16/*,32,64*/};+ mpint *b;
+ char *p;
+
+ printf("testconv \"%s\":\n", str);+ b = strtomp(str, nil, 16, nil);
+
+ for(i=0; i < nelem(base); i++){+ p = mptoa(b, base[i], nil, 0);
+ if(p == nil){+ fprintf(stderr, "mptoa(base=%d) -> nil\n", base[i]);
+ exit(2);
+ }
+ printf("base%d: %s = ", base[i], p);+ if(strtomp(p, nil, base[i], b) == nil){+ fprintf(stderr, "A strtomp(%s) -> nil\n", p);
+ exit(2);
+ }
+ free(p);
+ printf("%s\n", MPB(base[i], b));+
+ switch(base[i]){+ case 2:
+ case 8:
+ case 10:
+ case 16:
+ asprintf(&p, "%s", MPB(base[i], b));
+ printf("# %s = ", p);+ if(strtomp(p, nil, 0, b) == nil){+ fprintf(stderr, "B strtomp(%s) -> nil\n", p);
+ exit(2);
+ }
+ free(p);
+ printf("%s\n", MPB(base[i], b));+ break;
+ }
+
+ }
+
+ mpfree(b);
+}
+
+static void
+testshift(char *str)
+{+ mpint *b1, *b2;
+ int i;
+
+ b1 = strtomp(str, nil, 16, nil);
+ b2 = mpnew(0);
+ for(i = 0; i < 64; i++){+ mpleft(b1, i, b2);
+ printf("%2.2d %s\n", i, MP(b2));+ }
+ for(i = 0; i < 64; i++){+ mpright(b2, i, b1);
+ printf("%2.2d %s\n", i, MP(b1));+ }
+ mpfree(b1);
+ mpfree(b2);
+}
+
+static void
+testaddsub(char *str)
+{+ mpint *b1, *b2;
+ int i;
+
+ b1 = strtomp(str, nil, 16, nil);
+ b2 = mpnew(0);
+ for(i = 0; i < 16; i++){+ mpadd(b1, b2, b2);
+ printf("%2.2d %s\n", i, MP(b2));+ }
+ for(i = 0; i < 16; i++){+ mpsub(b2, b1, b2);
+ printf("%2.2d %s\n", i, MP(b2));+ }
+ mpfree(b1);
+ mpfree(b2);
+}
+
+static void
+testvecdigmuladd(char *str, mpdigit d)
+{+ mpint *b, *b2;
+ int i;
+ int64_t now;
+
+ b = strtomp(str, nil, 16, nil);
+ b2 = mpnew(0);
+
+ mpbits(b2, (b->top+1)*Dbits);
+ now = nsec();
+ for(i = 0; i < loops; i++){+ memset(b2->p, 0, b2->top*Dbytes);
+ mpvecdigmuladd(b->p, b->top, d, b2->p);
+ }
+ if(loops > 1)
+ printf("%"PRId64" ns for a %d*%d vecdigmul\n", (nsec()-now)/loops, b->top*Dbits, Dbits);+ mpnorm(b2);
+ printf("0 + %s * %x = %s\n", MP(b), d, MP(b2));+
+ mpfree(b);
+ mpfree(b2);
+}
+
+static void
+testvecdigmulsub(char *str, mpdigit d)
+{+ mpint *b, *b2;
+ int i;
+ int64_t now;
+
+ b = strtomp(str, nil, 16, nil);
+ b2 = mpnew(0);
+
+ mpbits(b2, (b->top+1)*Dbits);
+ now = nsec();
+ for(i = 0; i < loops; i++){+ memset(b2->p, 0, b2->top*Dbytes);
+ mpvecdigmulsub(b->p, b->top, d, b2->p);
+ }
+ if(loops > 1)
+ printf("%"PRId64" ns for a %d*%d vecdigmul\n", (nsec()-now)/loops, b->top*Dbits, Dbits);+ mpnorm(b2);
+ printf("0 - %s * %x = %s\n", MP(b), d, MP(b2));+
+ mpfree(b);
+ mpfree(b2);
+}
+
+static void
+testmul(char *str)
+{+ mpint *b, *b1, *b2;
+ int64_t now;
+ int i;
+
+ b = strtomp(str, nil, 16, nil);
+ b1 = mpcopy(b);
+ b2 = mpnew(0);
+
+ now = nsec();
+ for(i = 0; i < loops; i++)
+ mpmul(b, b1, b2);
+ if(loops > 1)
+ printf("%"PRId64" µs for a %d*%d mult\n", (nsec()-now)/(loops*1000),+ b->top*Dbits, b1->top*Dbits);
+ printf("%s * %s = %s\n", MP(b), MP(b1), MP(b2));+
+ mpfree(b);
+ mpfree(b1);
+ mpfree(b2);
+}
+
+static void
+testmul2(mpint *b, mpint *b1)
+{+ mpint *b2;
+ int64_t now;
+ int i;
+
+ b2 = mpnew(0);
+
+ now = nsec();
+ for(i = 0; i < loops; i++)
+ mpmul(b, b1, b2);
+ if(loops > 1)
+ printf("%"PRId64" µs for a %d*%d mult\n", (nsec()-now)/(loops*1000), b->top*Dbits, b1->top*Dbits);+ printf("%s * ", MP(b));+ printf("%s = ", MP(b1));+ printf("%s\n", MP(b2));+
+ mpfree(b2);
+}
+
+static void
+testdigdiv(char *str, mpdigit d)
+{+ mpint *b;
+ mpdigit q;
+ int i;
+ int64_t now;
+
+ b = strtomp(str, nil, 16, nil);
+ now = nsec();
+ for(i = 0; i < loops; i++)
+ mpdigdiv(b->p, d, &q);
+ if(loops > 1)
+ printf("%"PRId64" ns for a %d / %d div\n", (nsec()-now)/loops, 2*Dbits, Dbits);+ printf("%s / %x = %x\n", MP(b), d, q);+ mpfree(b);
+}
+
+static void
+testdiv(mpint *x, mpint *y)
+{+ mpint *b2, *b3;
+ int64_t now;
+ int i;
+
+ b2 = mpnew(0);
+ b3 = mpnew(0);
+ now = nsec();
+ for(i = 0; i < loops; i++)
+ mpdiv(x, y, b2, b3);
+ if(loops > 1)
+ printf("%"PRId64" µs for a %d/%d div\n", (nsec()-now)/(1000*loops),+ x->top*Dbits, y->top*Dbits);
+ printf("%s / %s = %s %s\n", MP(x), MP(y), MP(b2), MP(b3));+ mpfree(b2);
+ mpfree(b3);
+}
+
+static void
+testmod(mpint *x, mpint *y)
+{+ mpint *r;
+ int64_t now;
+ int i;
+
+ r = mpnew(0);
+ now = nsec();
+ for(i = 0; i < loops; i++)
+ mpmod(x, y, r);
+ if(loops > 1)
+ printf("%"PRId64" µs for a %d/%d mod\n", (nsec()-now)/(1000*loops),+ x->top*Dbits, y->top*Dbits);
+ printf("%s mod %s = %s\n", MP(x), MP(y), MP(r));+ mpfree(r);
+}
+
+static void
+testinvert(mpint *x, mpint *y)
+{+ mpint *r, *d1, *d2;
+ int64_t now;
+ int i;
+
+ r = mpnew(0);
+ d1 = mpnew(0);
+ d2 = mpnew(0);
+ now = nsec();
+ mpextendedgcd(x, y, r, d1, d2);
+ mpdiv(x, r, x, d1);
+ mpdiv(y, r, y, d1);
+ for(i = 0; i < loops; i++)
+ mpinvert(x, y, r);
+ if(loops > 1)
+ printf("%"PRId64" µs for a %d in %d invert\n", (nsec()-now)/(1000*loops),+ x->top*Dbits, y->top*Dbits);
+ printf("%s**-1 mod %s = %s\n", MP(x), MP(y), MP(r));+ mpmul(r, x, d1);
+ mpmod(d1, y, d2);
+ printf("%s*%s mod %s = %s\n", MP(x), MP(r), MP(y), MP(d2));+ mpfree(r);
+ mpfree(d1);
+ mpfree(d2);
+}
+
+static void
+testsub1(char *a, char *b)
+{+ mpint *b1, *b2, *b3;
+
+ b1 = strtomp(a, nil, 16, nil);
+ b2 = strtomp(b, nil, 16, nil);
+ b3 = mpnew(0);
+ mpsub(b1, b2, b3);
+ printf("%s - %s = %s\n", MP(b1), MP(b2), MP(b3));+}
+
+static void
+testmul1(char *a, char *b)
+{+ mpint *b1, *b2, *b3;
+
+ b1 = strtomp(a, nil, 16, nil);
+ b2 = strtomp(b, nil, 16, nil);
+ b3 = mpnew(0);
+ mpmul(b1, b2, b3);
+ printf("%s * %s = %s\n", MP(b1), MP(b2), MP(b3));+}
+
+static void
+testexp(char *base, char *exp, char *mod)
+{+ mpint *b, *e, *m, *res;
+ int i;
+ uint64_t now;
+
+ b = strtomp(base, nil, 16, nil);
+ e = strtomp(exp, nil, 16, nil);
+ res = mpnew(0);
+ if(mod != nil)
+ m = strtomp(mod, nil, 16, nil);
+ else
+ m = nil;
+ now = nsec();
+ for(i = 0; i < loops; i++)
+ mpexp(b, e, m, res);
+ if(loops > 1)
+ printf("%"PRIu64"µs for a %d to the %d bit exp\n", (nsec()-now)/(loops*1000),+ b->top*Dbits, e->top*Dbits);
+ if(m != nil)
+ printf("%s ^ %s mod %s == %s\n", MP(b), MP(e), MP(m), MP(res));+ else
+ printf("%s ^ %s == %s\n", MP(b), MP(e), MP(res));+ mpfree(b);
+ mpfree(e);
+ mpfree(res);
+ if(m != nil)
+ mpfree(m);
+}
+
+static void
+testgcd(void)
+{+ mpint *a, *b, *d, *x, *y, *t1, *t2;
+ int i;
+ uint64_t now, then;
+ uint64_t etime;
+
+ d = mpnew(0);
+ x = mpnew(0);
+ y = mpnew(0);
+ t1 = mpnew(0);
+ t2 = mpnew(0);
+
+ etime = 0;
+
+ a = strtomp("4EECAB3E04C4E6BC1F49D438731450396BF272B4D7B08F91C38E88ADCD281699889AFF872E2204C80CCAA8E460797103DE539D5DF8335A9B20C0B44886384F134C517287202FCA914D8A5096446B40CD861C641EF9C2730CB057D7B133F4C2B16BBD3D75FDDBD9151AAF0F9144AAA473AC93CF945DBFE0859FB685D5CBD0A8B3", nil, 16, nil);+ b = strtomp("C41CFBE4D4846F67A3DF7DE9921A49D3B42DC33728427AB159CEC8CBBDB12B5F0C244F1A734AEB9840804EA3C25036AD1B61AFF3ABBC247CD4B384224567A863A6F020E7EE9795554BCD08ABAD7321AF27E1E92E3DB1C6E7E94FAAE590AE9C48F96D93D178E809401ABE8A534A1EC44359733475A36A70C7B425125062B1142D", nil, 16, nil);+ mpextendedgcd(a, b, d, x, y);
+ printf("gcd %s*%s+%s*%s = %s?\n", MP(a), MP(x), MP(b), MP(y), MP(d));+ mpfree(a);
+ mpfree(b);
+
+ for(i = 0; i < loops; i++){+ a = mprand(2048, prng, nil);
+ b = mprand(2048, prng, nil);
+ then = nsec();
+ mpextendedgcd(a, b, d, x, y);
+ now = nsec();
+ etime += now-then;
+ mpmul(a, x, t1);
+ mpmul(b, y, t2);
+ mpadd(t1, t2, t2);
+ if(mpcmp(d, t2) != 0)
+ printf("%d gcd %s*%s+%s*%s != %s\n", i, MP(a), MP(x), MP(b), MP(y), MP(d));+// else
+// printf("%d euclid %s*%s+%s*%s == %s\n", i, MP(a), MP(x), MP(b), MP(y), MP(d));+ mpfree(a);
+ mpfree(b);
+ }
+
+ mpfree(x);
+ mpfree(y);
+ mpfree(d);
+ mpfree(t1);
+ mpfree(t2);
+
+ if(loops > 1)
+ printf("binary %"PRIu64"\n", etime);+}
+
+int
+main(int argc, char **argv)
+{+ mpint *x, *y;
+
+ if(argc == 3 && strcmp(argv[1], "-n") == 0)
+ loops = atoi(argv[2]);
+
+ memset(str, 0, sizeof(str));
+ D_PNAN = D_NNAN = strtod("+NaN", nil);+ D_PINF = D_NINF = strtod("+Inf", nil);+
+ union ieee754_double *d;
+ d = (union ieee754_double *)&D_NNAN;
+ d->ieee.negative = 1;
+ d = (union ieee754_double *)&D_NINF;
+ d->ieee.negative = 1;
+
+ srand(0);
+ mpsetminbits(2*Dbits);
+ testshift("1111111111111111");+ testaddsub("fffffffffffffffff");+ testdigdiv("1234567812345678", 0x76543218);+ testdigdiv("1ffff", 0xffff);+ testdigdiv("ffff", 0xffff);+ testdigdiv("fff", 0xffff);+ testdigdiv("effffffff", 0xffff);+ testdigdiv("ffffffff", 0x1);+ testdigdiv("ffffffff", 0);+ testdigdiv("200000000", 2);+ testdigdiv("ffffff00fffffff1", 0xfffffff1);+ testvecdigmuladd("fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff", 2);+ testconv("0");+ testconv("-abc0123456789abcedf");+ testconv("abc0123456789abcedf");+ testconv("ffffffff");+ testconv("aaaaaaaaaaaaaaaaa");+ testconv("1111111111111111");+ testconv("33333333333333333333333333333333");+
+ testvecdigmulsub("fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff", 2);+ testsub1("1FFFFFFFE00000000", "FFFFFFFE00000001");+ testmul1("ffffffff", "f");+ testmul("ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff");+ testmul1("100000000000000000000000000000000000000000000000000000002000000000000000000000000000000000000000000000000000000030000000000000000000000000000000000000000000000000000000400000000000000000000000000000000000000000000000000000004FFFFFFFFFFFFFFFE0000000200000000000000000000000000000003FFFFFFFFFFFFFFFE0000000200000000000000000000000000000002FFFFFFFFFFFFFFFE0000000200000000000000000000000000000001FFFFFFFFFFFFFFFE0000000200000000000000000000000000000000FFFFFFFFFFFFFFFE0000000200000000FFFFFFFE00000001", "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF");+ testmul1("1000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001", "1000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001");+ testmul1("1000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001", "1000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001000000000000000000000001");+ testmul1("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000", "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");+ x = mprand(256, prng, nil);
+ y = mprand(128, prng, nil);
+ testdiv(x, y);
+ x = mprand(2048, prng, nil);
+ y = mprand(1024, prng, nil);
+ testdiv(x, y);
+ x = mprand(4*1024, prng, nil);
+ y = mprand(4*1024, prng, nil);
+ testmul2(x, y);
+ testsub1("677132C9", "-A26559B6");+ testgcd();
+ x = mprand(512, prng, nil);
+ x->sign = -1;
+ y = mprand(256, prng, nil);
+ testdiv(x, y);
+ testmod(x, y);
+ x->sign = 1;
+ testinvert(y, x);
+ testexp("111111111", "222", "1000000000000000000000");+ testexp("ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff", "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff", "100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");+ return 0;
+}
--- a/meson.build
+++ b/meson.build
@@ -25,17 +25,22 @@
language: 'c',
)
-src = [
- '3rd/mt19937-64.c',
+mp = [
'3rd/mp/mpadd.c',
'3rd/mp/mpaux.c',
'3rd/mp/mpcmp.c',
'3rd/mp/mpdigdiv.c',
'3rd/mp/mpdiv.c',
+ '3rd/mp/mpextendedgcd.c',
+ '3rd/mp/mpexp.c',
'3rd/mp/mpfmt.c',
+ '3rd/mp/mpinvert.c',
'3rd/mp/mpleft.c',
'3rd/mp/mplogic.c',
+ '3rd/mp/mpmod.c',
+ '3rd/mp/mpmodop.c',
'3rd/mp/mpmul.c',
+ '3rd/mp/mprand.c',
'3rd/mp/mpright.c',
'3rd/mp/mpsub.c',
'3rd/mp/mptobe.c',
@@ -52,11 +57,18 @@
'3rd/mp/mpvectscmp.c',
'3rd/mp/strtomp.c',
'3rd/mp/u16.c',
- '3rd/spooky.c',
+]
+
+utf = [
'3rd/utf/rune.c',
'3rd/utf/runeistype.c',
'3rd/utf/runetotype.c',
'3rd/utf/utfnlen.c',
+]
+
+src = [
+ '3rd/mt19937-64.c',
+ '3rd/spooky.c',
'bitvector.c',
'builtins.c',
'cvalues.c',
@@ -145,6 +157,8 @@
'flisp',
sources: [
src,
+ mp,
+ utf,
boot,
builtins,
],
@@ -158,6 +172,20 @@
'posix',
),
)
+
+mptest = executable(
+ 'mptest',
+ sources: [
+ '3rd/mp/test.c',
+ mp,
+ ],
+ include_directories: include_directories(
+ '3rd',
+ '3rd/mp',
+ 'posix',
+ ),
+)
+test('mp', mptest)tests_dir = join_paths(meson.current_source_dir(), 'test')
--- a/posix/mp.h
+++ b/posix/mp.h
@@ -17,8 +17,6 @@
};
};
-#define sysfatal(s) do{ fprintf(stderr, "%s\n", s); abort(); }while(0)-
#define mpdighi (mpdigit)(1U<<(Dbits-1))
#define DIGITS(x) ((x) >= -(Dbits-1) ? ((Dbits - 1 + (x))/Dbits) : 0)