1 #include <stdint.h>
2 #include <stdlib.h>
3 #include <stdio.h>
4 #include <errno.h>
5 #include <unistd.h>
6
7 static uint64_t seed = -1;
rand32(void)8 static uint32_t rand32(void)
9 {
10 seed = 6364136223846793005ull*seed + 1;
11 return seed >> 32;
12 }
rand64(void)13 static uint64_t rand64(void)
14 {
15 uint64_t u = rand32();
16 return u<<32 | rand32();
17 }
frand()18 static double frand(){return rand64() * 0x1p-64;}
frandf()19 static float frandf(){return rand32() * 0x1p-32f;}
frandl()20 static long double frandl(){return rand64() * 0x1p-64L;}
21
22 /* uniform random in [0,n), n > 0 must hold */
randn(uint64_t n)23 uint64_t randn(uint64_t n)
24 {
25 uint64_t r, m;
26
27 /* m is the largest multiple of n */
28 m = -1;
29 m -= m%n;
30 while ((r = rand64()) >= m);
31 return r%n;
32 }
33
34 /* uniform on [a,b] */
randint(uint64_t a,uint64_t b)35 uint64_t randint(uint64_t a, uint64_t b)
36 {
37 if (b < a) {
38 uint64_t t = b;
39 b = a;
40 a = t;
41 }
42 return a + randn(b - a + 1);
43 }
44
insert(uint64_t * tab,size_t len,uint64_t v)45 int insert(uint64_t *tab, size_t len, uint64_t v)
46 {
47 size_t i = v & (len-1);
48 size_t j = 1;
49
50 /* 0 means empty, v > 0 must hold */
51 while (tab[i]) {
52 if (tab[i] == v)
53 return -1;
54 i += j++;
55 i &= len-1;
56 }
57 tab[i] = v;
58 return 0;
59 }
60
shuffle2(uint64_t * p,uint64_t * q,size_t np,size_t nq)61 static void shuffle2(uint64_t *p, uint64_t *q, size_t np, size_t nq)
62 {
63 size_t i,r,t;
64
65 i = np+nq;
66 while (i > np) {
67 r = randn(i);
68 i--;
69 t = q[i-np];
70 if (r < np) {
71 q[i-np] = p[r];
72 p[r] = t;
73 } else {
74 q[i-np] = q[r-np];
75 q[r-np] = t;
76 }
77 }
78 }
79
80 /* choose k unique numbers from [0,n), k <= n */
choose(uint64_t n,size_t k,uint64_t * p)81 int choose(uint64_t n, size_t k, uint64_t *p)
82 {
83 uint64_t *tab;
84 size_t i, j, len;
85
86 if (n < k)
87 return -1;
88
89 if (n < 16) {
90 /* no alloc */
91 while (k)
92 if (randn(n--) < k)
93 p[--k] = n;
94 return 0;
95 }
96
97 if (k < 8) {
98 /* no alloc, n > 15 > 2*k */
99 for (i = 0; i < k;) {
100 p[i] = randn(n);
101 for (j = 0; p[j] != p[i]; j++);
102 if (j == i)
103 i++;
104 }
105 return 0;
106 }
107
108 if (n < 5*k && (n-k)*sizeof *tab < (size_t)-1) {
109 /* allocation is < 4*k */
110 tab = malloc((n-k) * sizeof *tab);
111 if (!tab)
112 return -1;
113 for (i = 0; i < k; i++)
114 p[i] = i;
115 for (; i < n; i++)
116 tab[i-k] = i;
117 if (n-k < k)
118 shuffle2(p, tab, k, n-k);
119 else
120 shuffle2(tab, p, n-k, k);
121 free(tab);
122 return 0;
123 }
124
125 /* allocation is < 4*k */
126 for (len = 16; len <= 2*k; len *= 2);
127 tab = calloc(len, sizeof *tab);
128 if (!tab)
129 return -1;
130 for (i = 0; i < k; i++)
131 while (insert(tab, len, randn(n)+1));
132 for (i = 0; i < len; i++)
133 if (tab[i])
134 *p++ = tab[i]-1;
135 free(tab);
136 return 0;
137 }
138
cmp64(const void * a,const void * b)139 static int cmp64(const void *a, const void *b)
140 {
141 const uint64_t *ua = a, *ub = b;
142 return *ua < *ub ? -1 : (*ua > *ub ? 1 : 0);
143 }
144
145 // todo: in place flip problem
146
147 /* choose k unique uint64_t numbers */
choose64(size_t k,uint64_t * p)148 int choose64(size_t k, uint64_t *p)
149 {
150 size_t i, c;
151
152 /* no alloc, collisions should be very rare */
153 for (i = 0; i < k; i++)
154 p[i] = rand64();
155 do {
156 c = 0;
157 qsort(p, k, sizeof *p, cmp64);
158 for (i = 1; i < k; i++)
159 if (p[i] == p[i-1]) {
160 p[i-1] = rand64();
161 c = 1;
162 }
163 } while (c);
164 return 0;
165 }
166
167 /* equidistant sampling with some randomness */
sample(uint64_t n,size_t k,uint64_t * p)168 int sample(uint64_t n, size_t k, uint64_t *p)
169 {
170 uint64_t a = 0;
171 uint64_t d = n/k;
172 size_t m = n%k;
173 size_t i, j;
174 uint64_t *q;
175
176 if (!d)
177 return -1;
178 q = malloc((m+1) * sizeof *q);
179 if (!q)
180 return -1;
181 if (choose(k, m, q))
182 return -1;
183 qsort(q, m, sizeof *q, cmp64);
184 q[m] = k;
185 for (i = j = 0; i < k; i++) {
186 uint64_t t;
187
188 while (q[j] < i)
189 j++;
190 if (q[j] == i)
191 t = d+1;
192 else
193 t = d;
194 p[i] = a + randn(t);
195 a += t;
196 }
197 free(q);
198 return 0;
199 }
200
201 /* [-inf,inf] uniform on representation */
genall(size_t k,uint64_t * p)202 int genall(size_t k, uint64_t *p)
203 {
204 size_t i;
205 uint64_t n, d;
206 d = 1;
207 d <<= 52;
208 if (sample(-2*d, k, p))
209 return -1;
210 n = 0x7ff;
211 n <<= 52;
212 for (i = 0; i < k; i++)
213 if (p[i] > n)
214 p[i] += d-1;
215 return 0;
216 }
217
218 /* [a,b) uniform on representation, 0 <= a <= b */
genab(size_t k,uint64_t a,uint64_t b,uint64_t * p)219 int genab(size_t k, uint64_t a, uint64_t b, uint64_t *p)
220 {
221 size_t i;
222
223 if (sample(b-a, k, p))
224 return -1;
225 for (i = 0; i < k; i++)
226 p[i] += a;
227 return 0;
228 }
229
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)
232
main(int argc,char * argv[])233 int main(int argc, char *argv[])
234 {
235 uint64_t k, i;
236 uint64_t *p;
237 double a,b,m;
238 char *e;
239 int opt;
240
241 k = 1000;
242 a = 0;
243 b = 1;
244 m = 1;
245 while ((opt = getopt(argc, argv, "n:a:b:m:s:")) != -1) {
246 switch(opt) {
247 case 'n':
248 k = strtoull(optarg,&e,0);
249 break;
250 case 'a':
251 a = strtod(optarg,&e);
252 if (a < 0)
253 goto usage;
254 break;
255 case 'b':
256 b = strtod(optarg,&e);
257 if (b < 0)
258 goto usage;
259 break;
260 case 'm':
261 m = strtod(optarg,&e);
262 break;
263 case 's':
264 seed = strtoull(optarg,&e,0);
265 break;
266 default:
267 usage:
268 fprintf(stderr, "usage: %s -n num -a absmin -b absmax -m mult -s seed\n", argv[0]);
269 return -1;
270 }
271 if (*e || errno)
272 goto usage;
273 }
274 if (!(a <= b))
275 goto usage;
276 p = malloc(k * sizeof *p);
277 if (!p)
278 return -1;
279 if (genab(k, asint(a), asint(b), p))
280 // if (genall(k,p))
281 return -1;
282 for (i = 0; i < k; i++)
283 // printf("0x%016llx\n", p[i]);
284 printf("%a\n", m*asfloat(p[i]));
285 return 0;
286 }
287