nsz Git - libc-test/blob - src/math/gen/rnd.c

   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;
   8 static uint32_t rand32(void)
   9 {
  10         seed = 6364136223846793005ull*seed + 1;
  11         return seed >> 32;
  12 }
  13 static uint64_t rand64(void)
  14 {
  15         uint64_t u = rand32();
  16         return u<<32 | rand32();
  17 }
  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;}
  21
  22 /* uniform random in [0,n), n > 0 must hold */
  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] */
  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
  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
  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 */
  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
 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 */
 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 */
 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 */
 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 */
 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
 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 }