From: Szabolcs Nagy <nsz@port70.net>
Date: Fri, 14 Dec 2012 17:39:30 +0000 (+0100)
Subject: math: add random float generator to gen, make check work for ulp tests
X-Git-Url: http://nsz.repo.hu/git/?a=commitdiff_plain;h=7fd3b84625cfa955207526eda74d9261c5324549;p=libc-test

math: add random float generator to gen, make check work for ulp tests
---

diff --git a/src/math/gen/Makefile b/src/math/gen/Makefile
index 9c1a452..f3c3429 100644
--- a/src/math/gen/Makefile
+++ b/src/math/gen/Makefile
@@ -2,7 +2,7 @@ CFLAGS=-I. -Wall -fno-builtin -ffloat-store -D_GNU_SOURCE
 U=mpfr
 T=$(wildcard t*.c)
 
-all: gen check mgen tof toa toe next prev
+all: gen check mgen tof toa toe next prev rnd
 
 #tx: tx.c $(U)/lib/libmpfr.a $(U)/lib/libgmp.a
 #	musl-gcc -o $@ $(CFLAGS) -I$(U)/include $^
diff --git a/src/math/gen/functions.h b/src/math/gen/functions.h
index f528de3..0bca4a8 100644
--- a/src/math/gen/functions.h
+++ b/src/math/gen/functions.h
@@ -1,7 +1,14 @@
 T(sinpi,       d_d)
+
 T(add,         dd_d)
+T(addf,        ff_f)
+T(addl,        ll_l)
 T(mul,         dd_d)
+T(mulf,        ff_f)
+T(mull,        ll_l)
 T(div,         dd_d)
+T(divf,        ff_f)
+T(divl,        ll_l)
 
 T(acos,        d_d)
 T(acosf,       f_f)
diff --git a/src/math/gen/gen.c b/src/math/gen/gen.c
index 6bc3c96..54acbed 100644
--- a/src/math/gen/gen.c
+++ b/src/math/gen/gen.c
@@ -36,7 +36,7 @@ static int print(const char *fmt, struct t *t, char *buf, int n);
 
 // TODO: many output, fmt->ulp
 struct fun;
-static int check(struct t *want, struct t *got, struct fun *f, float ulpthres);
+static int check(struct t *want, struct t *got, struct fun *f, float ulpthres, float *abserr);
 
 struct fun {
 	char *name;
@@ -57,17 +57,28 @@ int main(int argc, char *argv[])
 	struct t t;
 	struct t tread;
 	struct fun *f = 0;
+	double ulpthres = 1.0;
+	float maxerr = 0;
+	float abserr;
+	struct t terr;
 
-	if (argc < 2) {
-		fprintf(stderr, "%s func\n", argv[0]);
-		return 1;
-	}
 	p = strrchr(argv[0], '/');
 	if (!p)
 		p = argv[0];
 	else
 		p++;
 	checkmode = strcmp(p, "check") == 0;
+	if (argc < 2) {
+		fprintf(stderr, "%s func%s\n", argv[0], checkmode ? " ulpthres" : "");
+		return 1;
+	}
+	if (argc > 2 && checkmode) {
+		ulpthres = strtod(argv[2], &p);
+		if (*p) {
+			fprintf(stderr, "invalid ulperr %s\n", argv[2]);
+			return 1;
+		}
+	}
 	for (i = 0; i < sizeof fun/sizeof *fun; i++)
 		if (strcmp(fun[i].name, argv[1]) == 0) {
 			f = fun + i;
@@ -88,11 +99,15 @@ int main(int argc, char *argv[])
 		if (f->mpf(&t))
 			fprintf(stderr, "error mpf %s, line %d\n", f->name, i);
 		if (checkmode) {
-			if (check(&tread, &t, f, 1.0)) {
+			if (check(&tread, &t, f, ulpthres, &abserr)) {
 				print(f->fmt, &tread, buf, sizeof buf);
 				fputs(buf, stdout);
-				print(f->fmt, &t, buf, sizeof buf);
-				fputs(buf, stdout);
+//				print(f->fmt, &t, buf, sizeof buf);
+//				fputs(buf, stdout);
+			}
+			if (abserr > maxerr) {
+				maxerr = abserr;
+				terr = tread;
 			}
 		} else {
 			if (print(f->fmt, &t, buf, sizeof buf))
@@ -100,16 +115,21 @@ int main(int argc, char *argv[])
 			fputs(buf, stdout);
 		}
 	}
+	if (checkmode && maxerr) {
+		printf("// maxerr: %f, ", maxerr);
+		print(f->fmt, &terr, buf, sizeof buf);
+		fputs(buf, stdout);
+	}
 	return 0;
 }
 
-static int check(struct t *want, struct t *got, struct fun *f, float ulpthres)
+static int check(struct t *want, struct t *got, struct fun *f, float ulpthres, float *abserr)
 {
 	int err = 0;
 	int m = INEXACT|UNDERFLOW; // TODO: dont check inexact and underflow for now
 
 	if ((got->e|m) != (want->e|m)) {
-		fprintf(stdout, "%s %s(%La,%La)==%La except: want %s",
+		fprintf(stdout, "//%s %s(%La,%La)==%La except: want %s",
 			rstr(want->r), f->name, want->x, want->x2, want->y, estr(want->e));
 		fprintf(stdout, " got %s\n", estr(got->e));
 		err++;
@@ -135,9 +155,10 @@ static int check(struct t *want, struct t *got, struct fun *f, float ulpthres)
 			return -1;
 
 		d = scalbnl(got->y - want->y, -n);
-		if (fabsf(d + want->dy) <= ulpthres)
+		*abserr = fabsf(d + want->dy);
+		if (*abserr <= ulpthres)
 			return err;
-		fprintf(stdout, "%s %s(%La,%La) want %La got %La ulperr %.3f = %a + %a\n",
+		fprintf(stdout, "//%s %s(%La,%La) want %La got %La ulperr %.3f = %a + %a\n",
 			rstr(want->r), f->name, want->x, want->x2, want->y, got->y, d + want->dy, d, want->dy);
 		err++;
 	}
diff --git a/src/math/gen/mp.c b/src/math/gen/mp.c
index 34b4f65..eca0b69 100644
--- a/src/math/gen/mp.c
+++ b/src/math/gen/mp.c
@@ -350,8 +350,14 @@ static int wrap_sinpi(mpfr_t my, const mpfr_t mx, mpfr_rnd_t r)
 int mpsinpi(struct t *t) { return mpd1(t, wrap_sinpi); }
 
 int mpadd(struct t *t) { return mpd2(t, mpfr_add); }
+int mpaddf(struct t *t) { return mpf2(t, mpfr_add); }
+int mpaddl(struct t *t) { return mpl2(t, mpfr_add); }
 int mpmul(struct t *t) { return mpd2(t, mpfr_mul); }
+int mpmulf(struct t *t) { return mpf2(t, mpfr_mul); }
+int mpmull(struct t *t) { return mpl2(t, mpfr_mul); }
 int mpdiv(struct t *t) { return mpd2(t, mpfr_div); }
+int mpdivf(struct t *t) { return mpf2(t, mpfr_div); }
+int mpdivl(struct t *t) { return mpl2(t, mpfr_div); }
 
 int mpacos(struct t *t) { return mpd1(t, mpfr_acos); }
 int mpacosf(struct t *t) { return mpf1(t, mpfr_acos); }
diff --git a/src/math/gen/mplibm.c b/src/math/gen/mplibm.c
index 2151259..141fa85 100644
--- a/src/math/gen/mplibm.c
+++ b/src/math/gen/mplibm.c
@@ -57,12 +57,26 @@ static int mpl2(struct t *s, long double (*f)(long double, long double))
 static double sinpi(double x) { return sin(3.141592653589793238*x); }
 int mpsinpi(struct t *t) { return mpd1(t, sinpi); }
 
-static double add(double x, double y) { double z = x + y; return z; }
+
+#define OP(n,op,t) static t n(t x, t y) { t z = x op y; return z; }
+OP(add,+,double)
+OP(addf,+,float)
+OP(addl,+,long double)
+OP(mul,*,double)
+OP(mulf,*,float)
+OP(mull,*,long double)
+OP(div,/,double)
+OP(divf,/,float)
+OP(divl,/,long double)
 int mpadd(struct t *t) { return mpd2(t, add); }
-static double mul(double x, double y) { double z = x * y; return z; }
+int mpaddf(struct t *t) { return mpf2(t, addf); }
+int mpaddl(struct t *t) { return mpl2(t, addl); }
 int mpmul(struct t *t) { return mpd2(t, mul); }
-static double div(double x, double y) { double z = x / y; return z; }
+int mpmulf(struct t *t) { return mpf2(t, mulf); }
+int mpmull(struct t *t) { return mpl2(t, mull); }
 int mpdiv(struct t *t) { return mpd2(t, div); }
+int mpdivf(struct t *t) { return mpf2(t, divf); }
+int mpdivl(struct t *t) { return mpl2(t, divl); }
 
 int mpacos(struct t *t) { return mpd1(t, acos); }
 int mpacosf(struct t *t) { return mpf1(t, acosf); }
diff --git a/src/math/gen/rnd.c b/src/math/gen/rnd.c
new file mode 100644
index 0000000..21ca038
--- /dev/null
+++ b/src/math/gen/rnd.c
@@ -0,0 +1,286 @@
+#include <stdint.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <errno.h>
+#include <unistd.h>
+
+static uint64_t seed = -1;
+static uint32_t rand32(void)
+{
+	seed = 6364136223846793005ull*seed + 1;
+	return seed >> 32;
+}
+static uint64_t rand64(void)
+{
+	uint64_t u = rand32();
+	return u<<32 | rand32();
+}
+static double frand(){return rand64() * 0x1p-64;}
+static float frandf(){return rand32() * 0x1p-32f;}
+static long double frandl(){return rand64() * 0x1p-64L;}
+
+/* uniform random in [0,n), n > 0 must hold */
+uint64_t randn(uint64_t n)
+{
+	uint64_t r, m;
+
+	/* m is the largest multiple of n */
+	m = -1;
+	m -= m%n;
+	while ((r = rand64()) >= m);
+	return r%n;
+}
+
+/* uniform on [a,b] */
+uint64_t randint(uint64_t a, uint64_t b)
+{
+	if (b < a) {
+		uint64_t t = b;
+		b = a;
+		a = t;
+	}
+	return a + randn(b - a + 1);
+}
+
+int insert(uint64_t *tab, size_t len, uint64_t v)
+{
+	size_t i = v & (len-1);
+	size_t j = 1;
+
+	/* 0 means empty, v > 0 must hold */
+	while (tab[i]) {
+		if (tab[i] == v)
+			return -1;
+		i += j++;
+		i &= len-1;
+	}
+	tab[i] = v;
+	return 0;
+}
+
+static void shuffle2(uint64_t *p, uint64_t *q, size_t np, size_t nq)
+{
+	size_t i,r,t;
+
+	i = np+nq;
+	while (i > np) {
+		r = randn(i);
+		i--;
+		t = q[i-np];
+		if (r < np) {
+			q[i-np] = p[r];
+			p[r] = t;
+		} else {
+			q[i-np] = q[r-np];
+			q[r-np] = t;
+		}
+	}
+}
+
+/* choose k unique numbers from [0,n), k <= n */
+int choose(uint64_t n, size_t k, uint64_t *p)
+{
+	uint64_t *tab;
+	size_t i, j, len;
+
+	if (n < k)
+		return -1;
+
+	if (n < 16) {
+		/* no alloc */
+		while (k)
+			if (randn(n--) < k)
+				p[--k] = n;
+		return 0;
+	}
+
+	if (k < 8) {
+		/* no alloc, n > 15 > 2*k */
+		for (i = 0; i < k;) {
+			p[i] = randn(n);
+			for (j = 0; p[j] != p[i]; j++);
+			if (j == i)
+				i++;
+		}
+		return 0;
+	}
+
+	if (n < 5*k && (n-k)*sizeof *tab < (size_t)-1) {
+		/* allocation is < 4*k */
+		tab = malloc((n-k) * sizeof *tab);
+		if (!tab)
+			return -1;
+		for (i = 0; i < k; i++)
+			p[i] = i;
+		for (; i < n; i++)
+			tab[i-k] = i;
+		if (n-k < k)
+			shuffle2(p, tab, k, n-k);
+		else
+			shuffle2(tab, p, n-k, k);
+		free(tab);
+		return 0;
+	}
+
+	/* allocation is < 4*k */
+	for (len = 16; len <= 2*k; len *= 2);
+	tab = calloc(len, sizeof *tab);
+	if (!tab)
+		return -1;
+	for (i = 0; i < k; i++)
+		while (insert(tab, len, randn(n)+1));
+	for (i = 0; i < len; i++)
+		if (tab[i])
+			*p++ = tab[i]-1;
+	free(tab);
+	return 0;
+}
+
+static int cmp64(const void *a, const void *b)
+{
+	const uint64_t *ua = a, *ub = b;
+	return *ua < *ub ? -1 : (*ua > *ub ? 1 : 0);
+}
+
+// todo: in place flip problem
+
+/* choose k unique uint64_t numbers */
+int choose64(size_t k, uint64_t *p)
+{
+	size_t i, c;
+
+	/* no alloc, collisions should be very rare */
+	for (i = 0; i < k; i++)
+		p[i] = rand64();
+	do {
+		c = 0;
+		qsort(p, k, sizeof *p, cmp64);
+		for (i = 1; i < k; i++)
+			if (p[i] == p[i-1]) {
+				p[i-1] = rand64();
+				c = 1;
+			}
+	} while (c);
+	return 0;
+}
+
+/* equidistant sampling with some randomness */
+int sample(uint64_t n, size_t k, uint64_t *p)
+{
+	uint64_t a = 0;
+	uint64_t d = n/k;
+	size_t m = n%k;
+	size_t i, j;
+	uint64_t *q;
+
+	if (!d)
+		return -1;
+	q = malloc((m+1) * sizeof *q);
+	if (!q)
+		return -1;
+	if (choose(k, m, q))
+		return -1;
+	qsort(q, m, sizeof *q, cmp64);
+	q[m] = k;
+	for (i = j = 0; i < k; i++) {
+		uint64_t t;
+
+		while (q[j] < i)
+			j++;
+		if (q[j] == i)
+			t = d+1;
+		else
+			t = d;
+		p[i] = a + randn(t);
+		a += t;
+	}
+	free(q);
+	return 0;
+}
+
+/* [-inf,inf] uniform on representation */
+int genall(size_t k, uint64_t *p)
+{
+	size_t i;
+	uint64_t n, d;
+	d = 1;
+	d <<= 52;
+	if (sample(-2*d, k, p))
+		return -1;
+	n = 0x7ff;
+	n <<= 52;
+	for (i = 0; i < k; i++)
+		if (p[i] > n)
+			p[i] += d-1;
+	return 0;
+}
+
+/* [a,b) uniform on representation, 0 <= a <= b */
+int genab(size_t k, uint64_t a, uint64_t b, uint64_t *p)
+{
+	size_t i;
+
+	if (sample(b-a, k, p))
+		return -1;
+	for (i = 0; i < k; i++)
+		p[i] += a;
+	return 0;
+}
+
+#define asfloat(x) ((union{uint64_t af_i; double af_f;}){.af_i=x}.af_f)
+#define asint(x)   ((union{uint64_t af_i; double af_f;}){.af_f=x}.af_i)
+
+int main(int argc, char *argv[])
+{
+	uint64_t k, i;
+	uint64_t *p;
+	double a,b,m;
+	char *e;
+	int opt;
+
+	k = 1000;
+	a = 0;
+	b = 1;
+	m = 1;
+	while ((opt = getopt(argc, argv, "n:a:b:m:s:")) != -1) {
+		switch(opt) {
+		case 'n':
+			k = strtoull(optarg,&e,0);
+			break;
+		case 'a':
+			a = strtod(optarg,&e);
+			if (a < 0)
+				goto usage;
+			break;
+		case 'b':
+			b = strtod(optarg,&e);
+			if (b < 0)
+				goto usage;
+			break;
+		case 'm':
+			m = strtod(optarg,&e);
+			break;
+		case 's':
+			seed = strtoull(optarg,&e,0);
+			break;
+		default:
+usage:
+			fprintf(stderr, "usage: %s -n num -a absmin -b absmax -m mult -s seed\n", argv[0]);
+			return -1;
+		}
+		if (*e || errno)
+			goto usage;
+	}
+	if (!(a <= b))
+		goto usage;
+	p = malloc(k * sizeof *p);
+	if (!p)
+		return -1;
+	if (genab(k, asint(a), asint(b), p))
+//	if (genall(k,p))
+		return -1;
+	for (i = 0; i < k; i++)
+//		printf("0x%016llx\n", p[i]);
+		printf("%a\n", m*asfloat(p[i]));
+	return 0;
+}