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