fork(1) download
copy 
  1. // Random generator using jrand48 as bitstream. (3.01)
  2.  
  3. #include <stdlib.h>
  4. #include <stdint.h>
  5. #include <limits.h>
  6. #include <assert.h>
  7. #include <math.h>
  8. #include <float.h>
  9. #include <stdio.h>
  10. #include <time.h>
  11.  
  12. _Static_assert((-1 & 3) == 3, "Not 2's complement");
  13. _Static_assert(UINT_MAX == (unsigned int) INT_MAX - INT_MIN, "Bad integer");
  14.  
  15. #define BITS(x) (sizeof(x) * CHAR_BIT)
  16.  
  17. int bitlength(unsigned long x)
  18. {
  19.     return x ? BITS(x) - __builtin_clzl(x) : 0;
  20. }
  21.  
  22. // Private.
  23.  
  24. static _Thread_local unsigned short generator_state[3] = {
  25.     0x330E, 0xABCD, 0x1234};
  26.  
  27. static uint_fast32_t generator_next(void)
  28. {
  29.     return jrand48(generator_state) + (UINT32_C(1) << 31);
  30. }
  31.  
  32. static uint_fast32_t takebits(uint_fast32_t x, int available, int take)
  33. {
  34.     assert(take >= 1);
  35.     assert(take <= available);
  36.     return (x >> (available - take)) & (UINT32_C(-1) >> (BITS(x) - take));
  37. }
  38.  
  39. // Public.
  40.  
  41. void randseed(uint_least64_t seed)
  42. {
  43.     generator_state[0] = seed;
  44.     generator_state[1] = seed >> 16;
  45.     generator_state[2] = seed >> 32;
  46. }
  47.  
  48. uint_least64_t randbits(int n)
  49. {
  50.     assert(n <= 64);
  51.     static _Thread_local uint_fast32_t reservoir = 0;
  52.     static _Thread_local int available = 0;
  53.  
  54.     uint_fast64_t r = 0;
  55.  
  56.     while (n > 0)
  57.     {
  58.         if (available == 0)
  59.         {
  60.             reservoir = generator_next();
  61.             available = 32;
  62.         }
  63.  
  64.         int take = (available < n) ? available : n;
  65.         r = (r << take) | takebits(reservoir, available, take);
  66.         available -= take;
  67.         n -= take;
  68.     }
  69.     return r;
  70. }
  71.  
  72. uint_least32_t randuint(uint_least32_t max)
  73. {
  74.     max &= UINT32_C(0xFFFFFFFF);
  75.  
  76.     if (max == 0)
  77.         return 0;
  78.  
  79.     int bits = bitlength(max);
  80.     uint_fast32_t res_max = UINT32_C(-1) >> (BITS(max) - bits);
  81.     uint_fast32_t res = randbits(bits);
  82.  
  83.     if (max == res_max)
  84.         return res;
  85.  
  86.     uint_fast32_t mod = max + 1;
  87.     uint_fast32_t lim = res_max + 1;
  88.     uint_fast32_t min = (lim - mod) % mod;
  89.  
  90.     while (res < min)
  91.         res = randbits(bits);
  92.     res %= mod;
  93.     return res;
  94. }
  95.  
  96. double randreal(void)
  97. {
  98.     double r = randbits(DBL_MANT_DIG) / exp2(DBL_MANT_DIG);
  99.     assert(r < 1.0);
  100.     return r;
  101. }
  102.  
  103. _Bool randbool(double p)
  104. {
  105.     assert(0.0 <= p && p <= 1.0);
  106.     return randreal() < p;
  107. }
  108.  
  109. // Test.
  110.  
  111. int uniform_int_distribution_uint(int a, int b)
  112. {
  113.     assert(a <= b);
  114.     return randuint((unsigned int) b - a) + a;
  115. }
  116.  
  117. int uniform_int_distribution_real(int a, int b)
  118. {
  119.     assert(a <= b);
  120.     return randreal() * (1.0 + b - a) + a;
  121. }
  122.  
  123. int bernoulli_distribution_50_50(int a, int b)
  124. {
  125.     assert(a == 0 && b == 1);
  126.     return randbool(0.5);
  127. }
  128.  
  129. unsigned long random_device()
  130. {
  131.     unsigned long r;
  132.     size_t n = 0;
  133.     FILE *fp = fopen("/dev/urandom", "r");
  134.     if (fp)
  135.     {
  136.         n = fread(&r, sizeof r, 1, fp);
  137.         fclose(fp);
  138.     }
  139.     if (n == 1)
  140.         return r;
  141.     return time(0);
  142. }
  143.  
  144. double chi_square(int n, int k, int (*f)(int, int))
  145. {
  146.     int *h = calloc(k, sizeof *h);
  147.     assert(h != 0 || k == 0);
  148.     for (int j = 0; j < n; j++)
  149.     {
  150.         int i = f(0, k-1);
  151.         assert(0 <= i && i < k);
  152.         h[i] += 1;
  153.     }
  154.     double expect = (double) n / k;
  155.     double x2 = 0;
  156.     for (int i = 0; i < k; i++)
  157.         x2 += pow(h[i] - expect, 2) / expect;
  158.     free(h);
  159.     return x2;
  160. }
  161.  
  162. void test(int m, int n, int df, double x2, int (*f)(int, int),
  163.           const char *header)
  164. {
  165.     puts(header);
  166.     for (int i = 0; i < m; i++)
  167.     {
  168.         double r = chi_square(n, df, f);
  169.         printf(" significant: %-5s [%6.2f][%6.2f]\n",
  170.                r >= x2 ? "true" : "false", r, x2);
  171.     }
  172. }
  173.  
  174. int main(void)
  175. {
  176.     unsigned long seed = random_device();
  177.     printf("seed: %lu\n", seed);
  178.     randseed(seed);
  179.  
  180.     int m = 5;
  181.     int n = 10000;
  182.  
  183.     // PV=0.05 for all tests.
  184.  
  185.     test(m, n, 100, 124.34, uniform_int_distribution_uint,
  186.          "uniform_int_distribution_uint");
  187.     test(m, n, 100, 124.34, uniform_int_distribution_real,
  188.          "uniform_int_distribution_real");
  189.     test(m, n, 2, 5.99, bernoulli_distribution_50_50,
  190.          "bernoulli_distribution_50_50");
  191.     return 0;
  192. }
Success #stdin #stdout 0.01s 5288KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
seed: 2632378540274587558
uniform_int_distribution_uint
 significant: false [ 89.82][124.34]
 significant: false [120.12][124.34]
 significant: false [ 99.20][124.34]
 significant: false [ 99.50][124.34]
 significant: false [114.26][124.34]
uniform_int_distribution_real
 significant: false [ 91.08][124.34]
 significant: false [105.70][124.34]
 significant: false [102.44][124.34]
 significant: false [ 99.24][124.34]
 significant: false [108.32][124.34]
bernoulli_distribution_50_50
 significant: false [  0.58][  5.99]
 significant: false [  2.82][  5.99]
 significant: false [  0.49][  5.99]
 significant: false [  0.55][  5.99]
 significant: false [  0.71][  5.99]
https://ideone.com/Iuke69
language:
C (gcc 8.3)
created:
13 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!