7 static uint64_t seed = -1;
8 static uint32_t rand32(void)
10 seed = 6364136223846793005ull*seed + 1;
13 static uint64_t rand64(void)
15 uint64_t u = rand32();
16 return u<<32 | rand32();
18 static double frand(){return rand64() * 0x1p-64;}
19 static float frandf(){return rand32() * 0x1p-32f;}
20 static long double frandl(){return rand64() * 0x1p-64L;}
22 /* uniform random in [0,n), n > 0 must hold */
23 uint64_t randn(uint64_t n)
27 /* m is the largest multiple of n */
30 while ((r = rand64()) >= m);
34 /* uniform on [a,b] */
35 uint64_t randint(uint64_t a, uint64_t b)
42 return a + randn(b - a + 1);
45 int insert(uint64_t *tab, size_t len, uint64_t v)
47 size_t i = v & (len-1);
50 /* 0 means empty, v > 0 must hold */
61 static void shuffle2(uint64_t *p, uint64_t *q, size_t np, size_t nq)
80 /* choose k unique numbers from [0,n), k <= n */
81 int choose(uint64_t n, size_t k, uint64_t *p)
98 /* no alloc, n > 15 > 2*k */
101 for (j = 0; p[j] != p[i]; j++);
108 if (n < 5*k && (n-k)*sizeof *tab < (size_t)-1) {
109 /* allocation is < 4*k */
110 tab = malloc((n-k) * sizeof *tab);
113 for (i = 0; i < k; i++)
118 shuffle2(p, tab, k, n-k);
120 shuffle2(tab, p, n-k, k);
125 /* allocation is < 4*k */
126 for (len = 16; len <= 2*k; len *= 2);
127 tab = calloc(len, sizeof *tab);
130 for (i = 0; i < k; i++)
131 while (insert(tab, len, randn(n)+1));
132 for (i = 0; i < len; i++)
139 static int cmp64(const void *a, const void *b)
141 const uint64_t *ua = a, *ub = b;
142 return *ua < *ub ? -1 : (*ua > *ub ? 1 : 0);
145 // todo: in place flip problem
147 /* choose k unique uint64_t numbers */
148 int choose64(size_t k, uint64_t *p)
152 /* no alloc, collisions should be very rare */
153 for (i = 0; i < k; i++)
157 qsort(p, k, sizeof *p, cmp64);
158 for (i = 1; i < k; i++)
159 if (p[i] == p[i-1]) {
167 /* equidistant sampling with some randomness */
168 int sample(uint64_t n, size_t k, uint64_t *p)
178 q = malloc((m+1) * sizeof *q);
183 qsort(q, m, sizeof *q, cmp64);
185 for (i = j = 0; i < k; i++) {
201 /* [-inf,inf] uniform on representation */
202 int genall(size_t k, uint64_t *p)
208 if (sample(-2*d, k, p))
212 for (i = 0; i < k; i++)
218 /* [a,b) uniform on representation, 0 <= a <= b */
219 int genab(size_t k, uint64_t a, uint64_t b, uint64_t *p)
223 if (sample(b-a, k, p))
225 for (i = 0; i < k; i++)
230 #define asfloat(x) ((union{uint64_t af_i; double af_f;}){.af_i=x}.af_f)
231 #define asint(x) ((union{uint64_t af_i; double af_f;}){.af_f=x}.af_i)
233 int main(int argc, char *argv[])
245 while ((opt = getopt(argc, argv, "n:a:b:m:s:")) != -1) {
248 k = strtoull(optarg,&e,0);
251 a = strtod(optarg,&e);
256 b = strtod(optarg,&e);
261 m = strtod(optarg,&e);
264 seed = strtoull(optarg,&e,0);
268 fprintf(stderr, "usage: %s -n num -a absmin -b absmax -m mult -s seed\n", argv[0]);
276 p = malloc(k * sizeof *p);
279 if (genab(k, asint(a), asint(b), p))
282 for (i = 0; i < k; i++)
283 // printf("0x%016llx\n", p[i]);
284 printf("%a\n", m*asfloat(p[i]));