xref: /third_party/musl/libc-test/src/math/gen/rnd.c

#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;
}