fork download
copy 
  1. /**
  2.  * @file rand.004.c
  3.  * @ingroup experimental
  4.  * Random number generation using custom LCG as bitstream.
  5.  * @date 01/04/2025
  6.  */
  7.  
  8. #include <stdlib.h>
  9. #include <stdint.h>
  10. #include <limits.h>
  11. #include <assert.h>
  12. #include <math.h>
  13. #include <float.h>
  14. #include <stdio.h>
  15.  
  16. _Static_assert((-1 & 3) == 3, "Not 2's complement");
  17. _Static_assert(UINT_MAX == (unsigned int) INT_MAX - INT_MIN, "Bad integer");
  18.  
  19. //
  20. // Utility.
  21. //
  22.  
  23. #define BITS(x) (sizeof(x) * CHAR_BIT)
  24.  
  25. int bitlength(unsigned int x)
  26. {
  27.     return x ? BITS(x) - __builtin_clz(x) : 0;
  28. }
  29.  
  30. //
  31. // Private.
  32. //
  33.  
  34. #define GENERATOR_BITS 32
  35.  
  36. static _Thread_local uint_fast64_t generator_state = UINT64_C(0x1234ABCD330E);
  37.  
  38. static void generator_seed(uint_fast64_t seed)
  39. {
  40.     generator_state = seed;
  41.     generator_state %= (uint_fast64_t) 1 << 48;
  42. }
  43.  
  44. static uint_fast32_t generator_next(void)
  45. {
  46.     generator_state *= (uint_fast64_t) 25214903917;
  47.     generator_state += (uint_fast64_t) 11;
  48.     generator_state %= (uint_fast64_t) 1 << 48;
  49.     return generator_state >> 16;
  50. }
  51.  
  52. static uint_fast32_t takebits(uint_fast32_t x, int available, int take)
  53. {
  54.     assert(take >= 1);
  55.     assert(take <= available);
  56.     return (x >> (available - take)) & (UINT32_C(-1) >> (BITS(x) - take));
  57. }
  58.  
  59. //
  60. // Public.
  61. //
  62.  
  63. void randseed(uint_fast64_t seed)
  64. {
  65.     generator_seed(seed);
  66. }
  67.  
  68. uint_fast64_t randbits(int n)
  69. {
  70.     assert(n <= 64);
  71.     static _Thread_local uint_fast32_t reservoir = 0;
  72.     static _Thread_local int available = 0;
  73.  
  74.     uint_fast64_t r = 0;
  75.  
  76.     while (n > 0)
  77.     {
  78.         if (available == 0)
  79.         {
  80.             reservoir = generator_next();
  81.             available = GENERATOR_BITS;
  82.         }
  83.  
  84.         int take = (available < n) ? available : n;
  85.         r = (r << take) | takebits(reservoir, available, take);
  86.         available -= take;
  87.         n -= take;
  88.     }
  89.     return r;
  90. }
  91.  
  92. unsigned int randuint(unsigned int max)
  93. {
  94.     if (max == 0)
  95.         return 0;
  96.  
  97.     int n = bitlength(max);
  98.  
  99.     unsigned int bit_max = -1U >> (BITS(max) - n);
  100.     unsigned int r = randbits(n);
  101.  
  102.     if (max == bit_max)
  103.         return r;
  104.  
  105.     unsigned int mod = max + 1;
  106.     unsigned int min = (bit_max + 1 - mod) % mod;
  107.  
  108.     while (r < min)
  109.         r = randbits(n);
  110.     return r % mod;
  111. }
  112.  
  113. double randreal(void)
  114. {
  115.     double r = randbits(DBL_MANT_DIG) / exp2(DBL_MANT_DIG);
  116.     assert(r < 1.0);
  117.     return r;
  118. }
  119.  
  120. _Bool randbool(double p)
  121. {
  122.     assert(0.0 <= p && p <= 1.0);
  123.     return randreal() < p;
  124. }
  125.  
  126. //
  127. // Test.
  128. //
  129.  
  130. int uniform_int_distribution_uint(int a, int b)
  131. {
  132.     assert(a <= b);
  133.     return randuint((unsigned int) b - a) + a;
  134. }
  135.  
  136. int uniform_int_distribution_real(int a, int b)
  137. {
  138.     assert(a <= b);
  139.     return randreal() * (1.0 + b - a) + a;
  140. }
  141.  
  142. int bernoulli_distribution_50_50(int a, int b)
  143. {
  144.     assert(a == 0 && b == 1);
  145.     return randbool(0.5);
  146. }
  147.  
  148. unsigned long random_device()
  149. {
  150.     unsigned long r = -1;
  151.     FILE *fp = fopen("/dev/urandom", "r");
  152.     if (fp)
  153.     {
  154.         fread(&r, sizeof r, 1, fp);
  155.         fclose(fp);
  156.     }
  157.     return r;
  158. }
  159.  
  160. double chi_square(int n, int k, int (*f)(int, int))
  161. {
  162.     int *h = calloc(k, sizeof *h);
  163.     assert(h != 0 || k == 0);
  164.     for (int j = 0; j < n; j++)
  165.     {
  166.         int i = f(0, k-1);
  167.         assert(0 <= i && i < k);
  168.         h[i] += 1;
  169.     }
  170.     double expect = (double) n / k;
  171.     double x2 = 0;
  172.     for (int i = 0; i < k; i++)
  173.         x2 += pow(h[i] - expect, 2) / expect;
  174.     free(h);
  175.     return x2;
  176. }
  177.  
  178. void test(int m, int n, int df, double x2, int (*f)(int, int),
  179.           const char *header)
  180. {
  181.     puts(header);
  182.     for (int i = 0; i < m; i++)
  183.     {
  184.         double r = chi_square(n, df, f);
  185.         printf(" significant: %-5s [%6.2f][%6.2f]\n",
  186.                r >= x2 ? "true" : "false", r, x2);
  187.     }
  188. }
  189.  
  190. int main(void)
  191. {
  192.     unsigned long seed = random_device();
  193.     printf("seed: %lu\n", seed);
  194.     randseed(seed);
  195.  
  196.     int m = 5;
  197.     int n = 10000;
  198.  
  199.     // PV=0.05 for all tests.
  200.  
  201.     test(m, n, 100, 124.34, uniform_int_distribution_uint,
  202.          "uniform_int_distribution_uint");
  203.     test(m, n, 100, 124.34, uniform_int_distribution_real,
  204.          "uniform_int_distribution_real");
  205.     test(m, n, 2, 5.99, bernoulli_distribution_50_50,
  206.          "bernoulli_distribution_50_50");
  207.     return 0;
  208. }
Success #stdin #stdout 0.01s 5284KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
seed: 12977440896278311007
uniform_int_distribution_uint
 significant: false [108.74][124.34]
 significant: false [ 97.10][124.34]
 significant: false [103.78][124.34]
 significant: false [ 90.82][124.34]
 significant: false [ 97.66][124.34]
uniform_int_distribution_real
 significant: false [ 98.92][124.34]
 significant: false [ 77.58][124.34]
 significant: false [ 89.92][124.34]
 significant: false [100.56][124.34]
 significant: false [ 82.56][124.34]
bernoulli_distribution_50_50
 significant: false [  0.08][  5.99]
 significant: false [  0.05][  5.99]
 significant: false [  0.00][  5.99]
 significant: false [  0.08][  5.99]
 significant: false [  0.23][  5.99]
https://ideone.com/k1Zz8p
language:
C (gcc 8.3)
created:
11 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!