Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  1. #include <float.h>
  2.  
  3. // TODO: use large period prng
  4. static uint64_t seed = -1;
  5. static uint32_t rand32(void)
  6. {
  7.         seed = 6364136223846793005ULL*seed + 1;
  8.         return seed >> 32;
  9. }
  10. static uint64_t rand64(void)
  11. {
  12.         uint64_t u = rand32();
  13.         return u<<32 | rand32();
  14. }
  15. static double frand()
  16. {
  17.         return rand64() * 0x1p-64;
  18. }
  19. static float frandf()
  20. {
  21.         return rand32() * 0x1p-32f;
  22. }
  23. static long double frandl()
  24. {
  25.         return rand64() * 0x1p-64L
  26. #if LDBL_MANT_DIG > 64
  27. + rand64() * 0x1p-128L
  28. #endif
  29. ;
  30. }
  31.  
  32. void t_randseed(uint64_t s)
  33. {
  34.         seed = s;
  35. }
  36.  
  37. /* uniform random in [0,n), n > 0 must hold */
  38. uint64_t t_randn(uint64_t n)
  39. {
  40.         uint64_t r, m;
  41.  
  42.         /* m is the largest multiple of n */
  43.         m = -1;
  44.         m -= m%n;
  45.         while ((r = rand64()) >= m);
  46.         return r%n;
  47. }
  48.  
  49. /* uniform on [a,b], a <= b must hold */
  50. uint64_t t_randint(uint64_t a, uint64_t b)
  51. {
  52.         uint64_t n = b - a + 1;
  53.         if (n)
  54.                 return a + t_randn(n);
  55.         return rand64();
  56. }
  57.  
  58. /* shuffle the elements of p and q until the elements in p are well shuffled */
  59. static void shuffle2(uint64_t *p, uint64_t *q, size_t np, size_t nq)
  60. {
  61.         size_t r;
  62.         uint64_t t;
  63.  
  64.         while (np) {
  65.                 r = t_randn(nq+np--);
  66.                 t = p[np];
  67.                 if (r < nq) {
  68.                         p[np] = q[r];
  69.                         q[r] = t;
  70.                 } else {
  71.                         p[np] = p[r-nq];
  72.                         p[r-nq] = t;
  73.                 }
  74.         }
  75. }
  76.  
  77. /* shuffle the elements of p */
  78. void t_shuffle(uint64_t *p, size_t n)
  79. {
  80.         shuffle2(p,0,n,0);
  81. }
  82.  
  83. void t_randrange(uint64_t *p, size_t n)
  84. {
  85.         size_t i;
  86.         for (i = 0; i < n; i++)
  87.                 p[i] = i;
  88.         t_shuffle(p, n);
  89. }
  90.  
  91. /* hash table insert, 0 means empty, v > 0 must hold, len is power-of-2 */
  92. static int insert(uint64_t *tab, size_t len, uint64_t v)
  93. {
  94.         size_t i = v & (len-1);
  95.         size_t j = 1;
  96.  
  97.         while (tab[i]) {
  98.                 if (tab[i] == v)
  99.                         return -1;
  100.                 i += j++;
  101.                 i &= len-1;
  102.         }
  103.         tab[i] = v;
  104.         return 0;
  105. }
  106.  
  107. /* choose k unique numbers from [0,n), k <= n */
  108. int t_choose(uint64_t n, size_t k, uint64_t *p)
  109. {
  110.         uint64_t *tab;
  111.         size_t i, j, len;
  112.  
  113.         if (n < k)
  114.                 return -1;
  115.  
  116.         if (n < 16) {
  117.                 /* no alloc */
  118.                 while (k)
  119.                         if (t_randn(n--) < k)
  120.                                 p[--k] = n;
  121.                 return 0;
  122.         }
  123.  
  124.         if (k < 8) {
  125.                 /* no alloc, n > 15 > 2*k */
  126.                 for (i = 0; i < k;) {
  127.                         p[i] = t_randn(n);
  128.                         for (j = 0; p[j] != p[i]; j++);
  129.                         if (j == i)
  130.                                 i++;
  131.                 }
  132.                 return 0;
  133.         }
  134.  
  135.         // TODO: if k < n/k use k*log(k) solution without alloc
  136.  
  137.         if (n < 5*k && (n-k)*sizeof *tab < (size_t)-1) {
  138.                 /* allocation is n-k < 4*k */
  139.                 tab = malloc((n-k) * sizeof *tab);
  140.                 if (!tab)
  141.                         return -1;
  142.                 for (i = 0; i < k; i++)
  143.                         p[i] = i;
  144.                 for (; i < n; i++)
  145.                         tab[i-k] = i;
  146.                 if (k < n-k)
  147.                         shuffle2(p, tab, k, n-k);
  148.                 else
  149.                         shuffle2(tab, p, n-k, k);
  150.                 free(tab);
  151.                 return 0;
  152.         }
  153.  
  154.         /* allocation is 2*k <= len < 4*k */
  155.         for (len = 16; len < 2*k; len *= 2);
  156.         tab = calloc(len, sizeof *tab);
  157.         if (!tab)
  158.                 return -1;
  159.         for (i = 0; i < k; i++)
  160.                 while (insert(tab, len, t_randn(n)+1));
  161.         for (i = 0; i < len; i++)
  162.                 if (tab[i])
  163.                         *p++ = tab[i]-1;
  164.         free(tab);
  165.         return 0;
  166. }
  167.  
  168.