Blame src/math/gen/rnd.c

Szabolcs Nagy 7fd3b8
#include <stdint.h>
Szabolcs Nagy 7fd3b8
#include <stdlib.h>
Szabolcs Nagy 7fd3b8
#include <stdio.h>
Szabolcs Nagy 7fd3b8
#include <errno.h>
Szabolcs Nagy 7fd3b8
#include <unistd.h>
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
static uint64_t seed = -1;
Szabolcs Nagy 7fd3b8
static uint32_t rand32(void)
Szabolcs Nagy 7fd3b8
{
Szabolcs Nagy 7fd3b8
	seed = 6364136223846793005ull*seed + 1;
Szabolcs Nagy 7fd3b8
	return seed >> 32;
Szabolcs Nagy 7fd3b8
}
Szabolcs Nagy 7fd3b8
static uint64_t rand64(void)
Szabolcs Nagy 7fd3b8
{
Szabolcs Nagy 7fd3b8
	uint64_t u = rand32();
Szabolcs Nagy 7fd3b8
	return u<<32 | rand32();
Szabolcs Nagy 7fd3b8
}
Szabolcs Nagy 7fd3b8
static double frand(){return rand64() * 0x1p-64;}
Szabolcs Nagy 7fd3b8
static float frandf(){return rand32() * 0x1p-32f;}
Szabolcs Nagy 7fd3b8
static long double frandl(){return rand64() * 0x1p-64L;}
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
/* uniform random in [0,n), n > 0 must hold */
Szabolcs Nagy 7fd3b8
uint64_t randn(uint64_t n)
Szabolcs Nagy 7fd3b8
{
Szabolcs Nagy 7fd3b8
	uint64_t r, m;
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
	/* m is the largest multiple of n */
Szabolcs Nagy 7fd3b8
	m = -1;
Szabolcs Nagy 7fd3b8
	m -= m%n;
Szabolcs Nagy 7fd3b8
	while ((r = rand64()) >= m);
Szabolcs Nagy 7fd3b8
	return r%n;
Szabolcs Nagy 7fd3b8
}
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
/* uniform on [a,b] */
Szabolcs Nagy 7fd3b8
uint64_t randint(uint64_t a, uint64_t b)
Szabolcs Nagy 7fd3b8
{
Szabolcs Nagy 7fd3b8
	if (b < a) {
Szabolcs Nagy 7fd3b8
		uint64_t t = b;
Szabolcs Nagy 7fd3b8
		b = a;
Szabolcs Nagy 7fd3b8
		a = t;
Szabolcs Nagy 7fd3b8
	}
Szabolcs Nagy 7fd3b8
	return a + randn(b - a + 1);
Szabolcs Nagy 7fd3b8
}
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
int insert(uint64_t *tab, size_t len, uint64_t v)
Szabolcs Nagy 7fd3b8
{
Szabolcs Nagy 7fd3b8
	size_t i = v & (len-1);
Szabolcs Nagy 7fd3b8
	size_t j = 1;
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
	/* 0 means empty, v > 0 must hold */
Szabolcs Nagy 7fd3b8
	while (tab[i]) {
Szabolcs Nagy 7fd3b8
		if (tab[i] == v)
Szabolcs Nagy 7fd3b8
			return -1;
Szabolcs Nagy 7fd3b8
		i += j++;
Szabolcs Nagy 7fd3b8
		i &= len-1;
Szabolcs Nagy 7fd3b8
	}
Szabolcs Nagy 7fd3b8
	tab[i] = v;
Szabolcs Nagy 7fd3b8
	return 0;
Szabolcs Nagy 7fd3b8
}
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
static void shuffle2(uint64_t *p, uint64_t *q, size_t np, size_t nq)
Szabolcs Nagy 7fd3b8
{
Szabolcs Nagy 7fd3b8
	size_t i,r,t;
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
	i = np+nq;
Szabolcs Nagy 7fd3b8
	while (i > np) {
Szabolcs Nagy 7fd3b8
		r = randn(i);
Szabolcs Nagy 7fd3b8
		i--;
Szabolcs Nagy 7fd3b8
		t = q[i-np];
Szabolcs Nagy 7fd3b8
		if (r < np) {
Szabolcs Nagy 7fd3b8
			q[i-np] = p[r];
Szabolcs Nagy 7fd3b8
			p[r] = t;
Szabolcs Nagy 7fd3b8
		} else {
Szabolcs Nagy 7fd3b8
			q[i-np] = q[r-np];
Szabolcs Nagy 7fd3b8
			q[r-np] = t;
Szabolcs Nagy 7fd3b8
		}
Szabolcs Nagy 7fd3b8
	}
Szabolcs Nagy 7fd3b8
}
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
/* choose k unique numbers from [0,n), k <= n */
Szabolcs Nagy 7fd3b8
int choose(uint64_t n, size_t k, uint64_t *p)
Szabolcs Nagy 7fd3b8
{
Szabolcs Nagy 7fd3b8
	uint64_t *tab;
Szabolcs Nagy 7fd3b8
	size_t i, j, len;
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
	if (n < k)
Szabolcs Nagy 7fd3b8
		return -1;
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
	if (n < 16) {
Szabolcs Nagy 7fd3b8
		/* no alloc */
Szabolcs Nagy 7fd3b8
		while (k)
Szabolcs Nagy 7fd3b8
			if (randn(n--) < k)
Szabolcs Nagy 7fd3b8
				p[--k] = n;
Szabolcs Nagy 7fd3b8
		return 0;
Szabolcs Nagy 7fd3b8
	}
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
	if (k < 8) {
Szabolcs Nagy 7fd3b8
		/* no alloc, n > 15 > 2*k */
Szabolcs Nagy 7fd3b8
		for (i = 0; i < k;) {
Szabolcs Nagy 7fd3b8
			p[i] = randn(n);
Szabolcs Nagy 7fd3b8
			for (j = 0; p[j] != p[i]; j++);
Szabolcs Nagy 7fd3b8
			if (j == i)
Szabolcs Nagy 7fd3b8
				i++;
Szabolcs Nagy 7fd3b8
		}
Szabolcs Nagy 7fd3b8
		return 0;
Szabolcs Nagy 7fd3b8
	}
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
	if (n < 5*k && (n-k)*sizeof *tab < (size_t)-1) {
Szabolcs Nagy 7fd3b8
		/* allocation is < 4*k */
Szabolcs Nagy 7fd3b8
		tab = malloc((n-k) * sizeof *tab);
Szabolcs Nagy 7fd3b8
		if (!tab)
Szabolcs Nagy 7fd3b8
			return -1;
Szabolcs Nagy 7fd3b8
		for (i = 0; i < k; i++)
Szabolcs Nagy 7fd3b8
			p[i] = i;
Szabolcs Nagy 7fd3b8
		for (; i < n; i++)
Szabolcs Nagy 7fd3b8
			tab[i-k] = i;
Szabolcs Nagy 7fd3b8
		if (n-k < k)
Szabolcs Nagy 7fd3b8
			shuffle2(p, tab, k, n-k);
Szabolcs Nagy 7fd3b8
		else
Szabolcs Nagy 7fd3b8
			shuffle2(tab, p, n-k, k);
Szabolcs Nagy 7fd3b8
		free(tab);
Szabolcs Nagy 7fd3b8
		return 0;
Szabolcs Nagy 7fd3b8
	}
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
	/* allocation is < 4*k */
Szabolcs Nagy 7fd3b8
	for (len = 16; len <= 2*k; len *= 2);
Szabolcs Nagy 7fd3b8
	tab = calloc(len, sizeof *tab);
Szabolcs Nagy 7fd3b8
	if (!tab)
Szabolcs Nagy 7fd3b8
		return -1;
Szabolcs Nagy 7fd3b8
	for (i = 0; i < k; i++)
Szabolcs Nagy 7fd3b8
		while (insert(tab, len, randn(n)+1));
Szabolcs Nagy 7fd3b8
	for (i = 0; i < len; i++)
Szabolcs Nagy 7fd3b8
		if (tab[i])
Szabolcs Nagy 7fd3b8
			*p++ = tab[i]-1;
Szabolcs Nagy 7fd3b8
	free(tab);
Szabolcs Nagy 7fd3b8
	return 0;
Szabolcs Nagy 7fd3b8
}
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
static int cmp64(const void *a, const void *b)
Szabolcs Nagy 7fd3b8
{
Szabolcs Nagy 7fd3b8
	const uint64_t *ua = a, *ub = b;
Szabolcs Nagy 7fd3b8
	return *ua < *ub ? -1 : (*ua > *ub ? 1 : 0);
Szabolcs Nagy 7fd3b8
}
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
// todo: in place flip problem
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
/* choose k unique uint64_t numbers */
Szabolcs Nagy 7fd3b8
int choose64(size_t k, uint64_t *p)
Szabolcs Nagy 7fd3b8
{
Szabolcs Nagy 7fd3b8
	size_t i, c;
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
	/* no alloc, collisions should be very rare */
Szabolcs Nagy 7fd3b8
	for (i = 0; i < k; i++)
Szabolcs Nagy 7fd3b8
		p[i] = rand64();
Szabolcs Nagy 7fd3b8
	do {
Szabolcs Nagy 7fd3b8
		c = 0;
Szabolcs Nagy 7fd3b8
		qsort(p, k, sizeof *p, cmp64);
Szabolcs Nagy 7fd3b8
		for (i = 1; i < k; i++)
Szabolcs Nagy 7fd3b8
			if (p[i] == p[i-1]) {
Szabolcs Nagy 7fd3b8
				p[i-1] = rand64();
Szabolcs Nagy 7fd3b8
				c = 1;
Szabolcs Nagy 7fd3b8
			}
Szabolcs Nagy 7fd3b8
	} while (c);
Szabolcs Nagy 7fd3b8
	return 0;
Szabolcs Nagy 7fd3b8
}
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
/* equidistant sampling with some randomness */
Szabolcs Nagy 7fd3b8
int sample(uint64_t n, size_t k, uint64_t *p)
Szabolcs Nagy 7fd3b8
{
Szabolcs Nagy 7fd3b8
	uint64_t a = 0;
Szabolcs Nagy 7fd3b8
	uint64_t d = n/k;
Szabolcs Nagy 7fd3b8
	size_t m = n%k;
Szabolcs Nagy 7fd3b8
	size_t i, j;
Szabolcs Nagy 7fd3b8
	uint64_t *q;
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
	if (!d)
Szabolcs Nagy 7fd3b8
		return -1;
Szabolcs Nagy 7fd3b8
	q = malloc((m+1) * sizeof *q);
Szabolcs Nagy 7fd3b8
	if (!q)
Szabolcs Nagy 7fd3b8
		return -1;
Szabolcs Nagy 7fd3b8
	if (choose(k, m, q))
Szabolcs Nagy 7fd3b8
		return -1;
Szabolcs Nagy 7fd3b8
	qsort(q, m, sizeof *q, cmp64);
Szabolcs Nagy 7fd3b8
	q[m] = k;
Szabolcs Nagy 7fd3b8
	for (i = j = 0; i < k; i++) {
Szabolcs Nagy 7fd3b8
		uint64_t t;
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
		while (q[j] < i)
Szabolcs Nagy 7fd3b8
			j++;
Szabolcs Nagy 7fd3b8
		if (q[j] == i)
Szabolcs Nagy 7fd3b8
			t = d+1;
Szabolcs Nagy 7fd3b8
		else
Szabolcs Nagy 7fd3b8
			t = d;
Szabolcs Nagy 7fd3b8
		p[i] = a + randn(t);
Szabolcs Nagy 7fd3b8
		a += t;
Szabolcs Nagy 7fd3b8
	}
Szabolcs Nagy 7fd3b8
	free(q);
Szabolcs Nagy 7fd3b8
	return 0;
Szabolcs Nagy 7fd3b8
}
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
/* [-inf,inf] uniform on representation */
Szabolcs Nagy 7fd3b8
int genall(size_t k, uint64_t *p)
Szabolcs Nagy 7fd3b8
{
Szabolcs Nagy 7fd3b8
	size_t i;
Szabolcs Nagy 7fd3b8
	uint64_t n, d;
Szabolcs Nagy 7fd3b8
	d = 1;
Szabolcs Nagy 7fd3b8
	d <<= 52;
Szabolcs Nagy 7fd3b8
	if (sample(-2*d, k, p))
Szabolcs Nagy 7fd3b8
		return -1;
Szabolcs Nagy 7fd3b8
	n = 0x7ff;
Szabolcs Nagy 7fd3b8
	n <<= 52;
Szabolcs Nagy 7fd3b8
	for (i = 0; i < k; i++)
Szabolcs Nagy 7fd3b8
		if (p[i] > n)
Szabolcs Nagy 7fd3b8
			p[i] += d-1;
Szabolcs Nagy 7fd3b8
	return 0;
Szabolcs Nagy 7fd3b8
}
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
/* [a,b) uniform on representation, 0 <= a <= b */
Szabolcs Nagy 7fd3b8
int genab(size_t k, uint64_t a, uint64_t b, uint64_t *p)
Szabolcs Nagy 7fd3b8
{
Szabolcs Nagy 7fd3b8
	size_t i;
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
	if (sample(b-a, k, p))
Szabolcs Nagy 7fd3b8
		return -1;
Szabolcs Nagy 7fd3b8
	for (i = 0; i < k; i++)
Szabolcs Nagy 7fd3b8
		p[i] += a;
Szabolcs Nagy 7fd3b8
	return 0;
Szabolcs Nagy 7fd3b8
}
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
#define asfloat(x) ((union{uint64_t af_i; double af_f;}){.af_i=x}.af_f)
Szabolcs Nagy 7fd3b8
#define asint(x)   ((union{uint64_t af_i; double af_f;}){.af_f=x}.af_i)
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
int main(int argc, char *argv[])
Szabolcs Nagy 7fd3b8
{
Szabolcs Nagy 7fd3b8
	uint64_t k, i;
Szabolcs Nagy 7fd3b8
	uint64_t *p;
Szabolcs Nagy 7fd3b8
	double a,b,m;
Szabolcs Nagy 7fd3b8
	char *e;
Szabolcs Nagy 7fd3b8
	int opt;
Szabolcs Nagy 7fd3b8
Szabolcs Nagy 7fd3b8
	k = 1000;
Szabolcs Nagy 7fd3b8
	a = 0;
Szabolcs Nagy 7fd3b8
	b = 1;
Szabolcs Nagy 7fd3b8
	m = 1;
Szabolcs Nagy 7fd3b8
	while ((opt = getopt(argc, argv, "n:a:b:m:s:")) != -1) {
Szabolcs Nagy 7fd3b8
		switch(opt) {
Szabolcs Nagy 7fd3b8
		case 'n':
Szabolcs Nagy 7fd3b8
			k = strtoull(optarg,&e,0);
Szabolcs Nagy 7fd3b8
			break;
Szabolcs Nagy 7fd3b8
		case 'a':
Szabolcs Nagy 7fd3b8
			a = strtod(optarg,&e);
Szabolcs Nagy 7fd3b8
			if (a < 0)
Szabolcs Nagy 7fd3b8
				goto usage;
Szabolcs Nagy 7fd3b8
			break;
Szabolcs Nagy 7fd3b8
		case 'b':
Szabolcs Nagy 7fd3b8
			b = strtod(optarg,&e);
Szabolcs Nagy 7fd3b8
			if (b < 0)
Szabolcs Nagy 7fd3b8
				goto usage;
Szabolcs Nagy 7fd3b8
			break;
Szabolcs Nagy 7fd3b8
		case 'm':
Szabolcs Nagy 7fd3b8
			m = strtod(optarg,&e);
Szabolcs Nagy 7fd3b8
			break;
Szabolcs Nagy 7fd3b8
		case 's':
Szabolcs Nagy 7fd3b8
			seed = strtoull(optarg,&e,0);
Szabolcs Nagy 7fd3b8
			break;
Szabolcs Nagy 7fd3b8
		default:
Szabolcs Nagy 7fd3b8
usage:
Szabolcs Nagy 7fd3b8
			fprintf(stderr, "usage: %s -n num -a absmin -b absmax -m mult -s seed\n", argv[0]);
Szabolcs Nagy 7fd3b8
			return -1;
Szabolcs Nagy 7fd3b8
		}
Szabolcs Nagy 7fd3b8
		if (*e || errno)
Szabolcs Nagy 7fd3b8
			goto usage;
Szabolcs Nagy 7fd3b8
	}
Szabolcs Nagy 7fd3b8
	if (!(a <= b))
Szabolcs Nagy 7fd3b8
		goto usage;
Szabolcs Nagy 7fd3b8
	p = malloc(k * sizeof *p);
Szabolcs Nagy 7fd3b8
	if (!p)
Szabolcs Nagy 7fd3b8
		return -1;
Szabolcs Nagy 7fd3b8
	if (genab(k, asint(a), asint(b), p))
Szabolcs Nagy 7fd3b8
//	if (genall(k,p))
Szabolcs Nagy 7fd3b8
		return -1;
Szabolcs Nagy 7fd3b8
	for (i = 0; i < k; i++)
Szabolcs Nagy 7fd3b8
//		printf("0x%016llx\n", p[i]);
Szabolcs Nagy 7fd3b8
		printf("%a\n", m*asfloat(p[i]));
Szabolcs Nagy 7fd3b8
	return 0;
Szabolcs Nagy 7fd3b8
}