fork download
copy 
  1. /*
  2.  * @file    main.c
  3.  * @author  Piotr Gregor <piotrek.gregor gmail com>
  4.  * @brief   Bernoulli distribution.
  5.  * @details Calls random device kernel module to get
  6.  *          sample of random bits and draws a Bernoulli
  7.  *          distribution from it.
  8. */
  9.  
  10.  
  11. #include <stdio.h>
  12. #include <stdlib.h>
  13. #include <stdint.h>
  14. #include <unistd.h>
  15. #include <fcntl.h> /* O_RDONLY */
  16.  
  17.  
  18. #define BITS_PER_BYTE 0x8u
  19. #define BIT_0 0x1u
  20. #define BIT_1 0x2u
  21. #define BIT_2 0x4u
  22. #define BIT_3 0x8u
  23. #define BIT_4 0x10u
  24. #define BIT_5 0x20u
  25. #define BIT_6 0x40u
  26. #define BIT_7 0x80u
  27.  
  28.  
  29. /* if defined use /dev/urandom (will not block),
  30.  * if not defined use /dev/random (may block)*/
  31. #define URANDOM_DEVICE 1
  32.  
  33. int
  34. log_int(uint32_t x)
  35. {
  36.     int8_t y; /* = log(x) */
  37.     switch (x)
  38.     {
  39.         case 0:
  40.             y = -1;
  41.             break;
  42.         case 1:
  43.             y = 0;
  44.             break;
  45.         case 2:
  46.             y = 1;
  47.             break;
  48.         case 4:
  49.             y = 2;
  50.             break;
  51.         case 8:
  52.             y = 3;
  53.             break;
  54.         case 16:
  55.             y = 4;
  56.             break;
  57.         case 32:
  58.             y = 5;
  59.             break;
  60.         case 64:
  61.             y = 6;
  62.             break;
  63.         case 128:
  64.             y = 7;
  65.             break;
  66.         case 256:
  67.             y = 8;
  68.             break;
  69.         case 512:
  70.             y = 9;
  71.             break;
  72.         case 1024:
  73.             y = 10;
  74.             break;
  75.         case 2048:
  76.             y = 11;
  77.             break;
  78.         case 4096:
  79.             y = 12;
  80.             break;
  81.         case 8192:
  82.             y = 13;
  83.             break;
  84.         case 16384:
  85.             y = 14;
  86.             break;
  87.         case 32768:
  88.             y = 15;
  89.             break;
  90.         case 65536:
  91.             y = 16;
  92.             break;
  93.         case 131072:
  94.             y = 17;
  95.             break;
  96.         case 262144:
  97.             y = 18;
  98.             break;
  99.         case 524288:
  100.             y = 19;
  101.             break;
  102.         case 1048576:
  103.             y = 20;
  104.             break;
  105.         case 2097152:
  106.             y = 21;
  107.             break;
  108.         case 4194304:
  109.             y = 22;
  110.             break;
  111.         case 8388608:
  112.             y = 23;
  113.             break;
  114.         case 16777216:
  115.             y = 24;
  116.             break;
  117.     }
  118.     return y;
  119. }
  120. /*
  121.  * @brief   Read @outlen bytes from random device
  122.  *          to array @out.
  123.  */
  124. int
  125. get_random_samples(unsigned char *out, size_t outlen)
  126. {
  127.     ssize_t res;
  128.  
  129. #ifdef URANDOM_DEVICE
  130.     int fd = open("/dev/urandom", O_RDONLY);
  131.     if (fd == -1)
  132.         return -1;
  133.        res = read(fd, out, outlen);
  134.        if (res < 0)
  135.        {
  136.            // error, unable to read /dev/urandom
  137.            close(fd);
  138.            return -2;
  139.        }
  140. #else
  141.     size_t read_n;
  142.     int fd = open("/dev/random", O_RDONLY);
  143.     if (fd == -1)
  144.         return -1;
  145.     read_n = 0;
  146.     while (read_n < outlen)
  147.     {
  148.        res = read(fd, out + read_n, outlen - read_n);
  149.        if (res < 0)
  150.        {
  151.            // error, unable to read /dev/random
  152.            close(fd);
  153.            return -3;
  154.        }
  155.        read_n += res;
  156.     }
  157. #endif /* URANDOM_DEVICE */
  158.     close(fd);
  159.     return 0;
  160. }
  161.  
  162. /*
  163.  * @brief   Draw vector of Bernoulli samples.
  164.  * @details @x and @resolution determines probability
  165.  *          of success in Bernoulli distribution
  166.  *          and accuracy of results: p = x/resolution.
  167.  * @param   resolution: number of segments per sample of output array 
  168.  *          as power of 2: max resolution supported is 2^24=16777216
  169.  * @param   x: determines used probability, x = [0, resolution - 1]
  170.  * @param   n: number of samples in result vector
  171.  */
  172. int
  173. get_bernoulli_samples(char *out, uint32_t n, uint32_t resolution, uint32_t x)
  174. {
  175.     int res;
  176.     size_t i, j;
  177.     uint32_t bytes_per_byte, word;
  178.     unsigned char *rnd_bytes;
  179.     uint32_t uniform_byte;
  180.     uint8_t bits_per_byte;
  181.  
  182.     if (out == NULL || n == 0 || resolution == 0 || x > (resolution -1))
  183.         return -1;
  184.  
  185.     bits_per_byte = log_int(resolution);
  186.     bytes_per_byte = bits_per_byte / BITS_PER_BYTE + 
  187.                         (bits_per_byte % BITS_PER_BYTE ? 1 : 0);
  188.     rnd_bytes = malloc(n * bytes_per_byte);
  189.     if (rnd_bytes == NULL)
  190.         return -2;
  191.     res = get_random_samples(rnd_bytes, n * bytes_per_byte);
  192.     if (res < 0)
  193.     {
  194.         free(rnd_bytes);
  195.         return -3;
  196.     }
  197.  
  198.     i = 0;
  199.     while (i < n)
  200.     {
  201.         /* get Bernoulli sample */
  202.         /* read byte */
  203.         j = 0;
  204.         word = 0;
  205.         while (j < bytes_per_byte)
  206.         {
  207.             word |= (rnd_bytes[i * bytes_per_byte + j] << (BITS_PER_BYTE * j));
  208.             ++j;
  209.         }
  210.         uniform_byte = word & ((1u << bits_per_byte) - 1);
  211.         /* decision */
  212.         if (uniform_byte < x)
  213.             out[i] = 1;
  214.         else
  215.             out[i] = 0;
  216.         ++i;
  217.     }
  218.  
  219.     free(rnd_bytes);    
  220.     return 0;
  221. }
  222.  
  223. void
  224. print_c(char *c, int n)
  225. {
  226.     int i = 0;
  227.     while (i < n)
  228.     {
  229.         printf( "[%d]", (int)c[i]);
  230.         ++i;
  231.     }
  232.     printf("\n");
  233. }
  234.  
  235. int
  236. main(void)
  237. {
  238.     int res;
  239.     char c[256];
  240.  
  241.     res = get_bernoulli_samples(c, 256, 256, 1); /* 1/256 = 0.0039 */
  242.     if (res < 0) return -1;
  243.     print_c(c, 40);
  244.  
  245.     res = get_bernoulli_samples(c, sizeof(c), 256, 2); /* 2/256 = 0.0078 */ 
  246.     if (res < 0) return -1;
  247.     print_c(c, 40);
  248.  
  249.     res = get_bernoulli_samples(c, sizeof(c), 256, 13); /* 13/256 = 0.0508 */
  250.     if (res < 0) return -1;
  251.     print_c(c, 40);
  252.  
  253.     res = get_bernoulli_samples(c, sizeof(c), 256*256, 328); /* 328/256^2 = 0.0050 */
  254.     if (res < 0) return -1;
  255.     print_c(c, 40);
  256.  
  257.     res = get_bernoulli_samples(c, sizeof(c), 256, 200); /* 200/256 = 0.7813 */
  258.     if (res < 0) return -1;
  259.     print_c(c, 40);
  260.  
  261.     res = get_bernoulli_samples(c, sizeof(c), 256*256, 200); /* 200/256^2 = 0.0031 */
  262.     if (res < 0) return -1;
  263.     print_c(c, 40);
  264.  
  265.     return 0;
  266. }
Success #stdin #stdout 0s 2292KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
[0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0]
[0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0]
[0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][1][0][0][0][1][0][0][0][0][0][0][0][0][0][0][0][0][0][1][0]
[0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0]
[1][0][0][1][1][0][1][0][0][0][1][1][1][0][1][0][1][0][1][1][1][0][1][1][1][1][1][1][1][0][1][1][1][1][1][1][1][1][0][1]
[0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0][0]
https://ideone.com/Xfy6Po
language:
C (gcc 8.3)
created:
8 years 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!