fork download
copy 
  1. // Pseudorandom number generation. (1.10)
  2.  
  3. #include <stdlib.h>
  4. #include <limits.h>
  5. #include <assert.h>
  6. #include <stdio.h>
  7. #include <time.h>
  8.  
  9. //
  10. // Utility.
  11. //
  12.  
  13. #define BITS(x) \
  14.     ((int) sizeof(x) * CHAR_BIT)
  15.  
  16. #define NEW_T(T, n) \
  17.     ((T *) allocate(n * sizeof(T)))
  18.  
  19. void *allocate(size_t n)
  20. {
  21.     void *r = malloc(n);
  22.     assert(r != 0 || n == 0);
  23.     return r;
  24. }
  25.  
  26. unsigned long uniform_ulong(unsigned short *state)
  27. {
  28.     unsigned long r = 0;
  29.     for (int i = 0; i < BITS(r); i += 32)
  30.         r += (unsigned long) jrand48(state) << i;
  31.     return r;
  32. }
  33.  
  34. //
  35. // Engine Base.
  36. //
  37.  
  38. typedef struct {
  39.     struct EngineBaseMethod *method;
  40. } EngineBase;
  41.  
  42. struct EngineBaseMethod {
  43.     long (*next)    (EngineBase *);
  44.     void (*discard) (EngineBase *, long);
  45.     void (*delete)  (EngineBase *);
  46. };
  47.  
  48. long engine_next(EngineBase *self)
  49. {
  50.     return self->method->next(self);
  51. }
  52.  
  53. void engine_discard(EngineBase *self, long n)
  54. {
  55.     self->method->discard(self, n);
  56. }
  57.  
  58. void engine_delete(EngineBase *self)
  59. {
  60.     self->method->delete(self);
  61. }
  62.  
  63. //
  64. // Linear Congruential.
  65. //
  66.  
  67. typedef struct {
  68.     EngineBase base;
  69.     long x, a, c, m; int d;
  70. } LinearCongruentialEngine;
  71.  
  72. static long linearCongruentialEngine_next(EngineBase *base)
  73. {
  74.     LinearCongruentialEngine *self = (LinearCongruentialEngine *) base;
  75.  
  76.     self->x = ((unsigned long long) self->x * self->a + self->c) % self->m;
  77.     return self->x >> self->d;
  78. }
  79.  
  80. static void linearCongruentialEngine_discard(EngineBase *base, long n)
  81. {
  82.     for (long i = 0; i < n; i++)
  83.         linearCongruentialEngine_next(base);
  84. }
  85.  
  86. static void linearCongruentialEngine_delete(EngineBase *base)
  87. {
  88.     free(base);
  89. }
  90.  
  91. static struct EngineBaseMethod linearCongruentialEngine_method = {
  92.     linearCongruentialEngine_next,
  93.     linearCongruentialEngine_discard,
  94.     linearCongruentialEngine_delete
  95. };
  96.  
  97. EngineBase *linearCongruentialEngine_new(long a, long c, long m, int d,
  98.         unsigned long seed)
  99. {
  100.     assert(m != 0);
  101.     assert(0 <= d && d < BITS(long));
  102.  
  103.     LinearCongruentialEngine *self = NEW_T(LinearCongruentialEngine, 1);
  104.  
  105.     self->base.method = &linearCongruentialEngine_method;
  106.     self->x = seed % m;
  107.     self->a = a;
  108.     self->c = c;
  109.     self->m = m;
  110.     self->d = d;
  111.     return (EngineBase *) self;
  112. }
  113.  
  114. //
  115. // Subtract With Carry.
  116. //
  117.  
  118. typedef struct {
  119.     EngineBase base;
  120.     long *x, m; int c, i, r, s;
  121. } SubtractWithCarryEngine;
  122.  
  123. static inline long next(long *x, long mask, int i_s, int i_r, int *carry)
  124. {
  125.     long y = x[i_s] - x[i_r] - *carry;
  126.     *carry = -(y >> (BITS(y) - 1));
  127.     return y & mask;
  128. }
  129.  
  130. static long subtractWithCarryEngine_next(EngineBase *base)
  131. {
  132.     SubtractWithCarryEngine *self = (SubtractWithCarryEngine *) base;
  133.  
  134.     self->i += 1;
  135.     long *x = self->x;
  136.     int   i = self->i, r = self->r;
  137.  
  138.     if (i == r)
  139.     {
  140.         long m = self->m;
  141.         int  c = self->c, s = self->s, t = r - s;
  142.  
  143.         for (i = 0; i < s; i++)
  144.             x[i] = next(x, m, i + t, i, &c);
  145.         for (i = s; i < r; i++)
  146.             x[i] = next(x, m, i - s, i, &c);
  147.         self->c = c;
  148.         self->i = i = 0;
  149.     }
  150.     return x[i];
  151. }
  152.  
  153. static void subtractWithCarryEngine_discard(EngineBase *base, long n)
  154. {
  155.     for (long i = 0; i < n; i++)
  156.         subtractWithCarryEngine_next(base);
  157. }
  158.  
  159. static void subtractWithCarryEngine_delete(EngineBase *base)
  160. {
  161.     SubtractWithCarryEngine *self = (SubtractWithCarryEngine *) base;
  162.  
  163.     if (self != 0)
  164.     {
  165.         free(self->x);
  166.         free(self);
  167.     }
  168. }
  169.  
  170. static struct EngineBaseMethod subtractWithCarryEngine_method = {
  171.     subtractWithCarryEngine_next,
  172.     subtractWithCarryEngine_discard,
  173.     subtractWithCarryEngine_delete
  174. };
  175.  
  176. EngineBase *subtractWithCarryEngine_new(int w, int s, int r, unsigned long seed)
  177. {
  178.     assert(0 < w && w < BITS(long));
  179.     assert(0 < s && s < r);
  180.  
  181.     SubtractWithCarryEngine *self = NEW_T(SubtractWithCarryEngine, 1);
  182.  
  183.     self->base.method = &subtractWithCarryEngine_method;
  184.     self->x = NEW_T(long, r);
  185.     self->m = -1UL >> (BITS(long) - w);
  186.     self->c = 0;
  187.     self->i = r-1;
  188.     self->r = r;
  189.     self->s = s;
  190.  
  191.     unsigned short state[3] = {0x330E, seed, seed>>16};
  192.  
  193.     for (int i = 0; i < r; i++)
  194.         self->x[i] = uniform_ulong(state) & self->m;
  195.     return (EngineBase *) self;
  196. }
  197.  
  198. //
  199. // Discard Block.
  200. //
  201.  
  202. typedef struct {
  203.     EngineBase base, *owned;
  204.     int p, r, i;
  205. } DiscardBlockEngine;
  206.  
  207. static long discardBlockEngine_next(EngineBase *base)
  208. {
  209.     DiscardBlockEngine *self = (DiscardBlockEngine *) base;
  210.  
  211.     if (self->i == 0)
  212.     {
  213.         engine_discard(self->owned, self->p - self->r);
  214.         self->i = self->r;
  215.     }
  216.     self->i -= 1;
  217.     return engine_next(self->owned);
  218. }
  219.  
  220. static void discardBlockEngine_discard(EngineBase *base, long n)
  221. {
  222.     for (long i = 0; i < n; i++)
  223.         discardBlockEngine_next(base);
  224. }
  225.  
  226. static void discardBlockEngine_delete(EngineBase *base)
  227. {
  228.     DiscardBlockEngine *self = (DiscardBlockEngine *) base;
  229.  
  230.     if (self != 0)
  231.     {
  232.         engine_delete(self->owned);
  233.         free(self);
  234.     }
  235. }
  236.  
  237. static struct EngineBaseMethod discardBlockEngine_method = {
  238.     discardBlockEngine_next,
  239.     discardBlockEngine_discard,
  240.     discardBlockEngine_delete,
  241. };
  242.  
  243. EngineBase *discardBlockEngine_new(EngineBase **unique, int p, int r)
  244. {
  245.     assert(0 < r && r <= p);
  246.  
  247.     DiscardBlockEngine *self = NEW_T(DiscardBlockEngine, 1);
  248.  
  249.     self->base.method = &discardBlockEngine_method;
  250.     self->owned = *unique; *unique = 0; // Transfer ownership.
  251.     self->p = p;
  252.     self->r = r;
  253.     self->i = r;
  254.     return (EngineBase *) self;
  255. }
  256.  
  257. //
  258. // Predefined Engines.
  259. //
  260.  
  261. EngineBase *default_rand_new(unsigned long seed)
  262. {
  263.     assert(BITS(long) >= 50);
  264.  
  265.     // 32-bit with a period of 2^48.
  266.     seed = seed * 0x10000L + 0x330E;
  267.     return linearCongruentialEngine_new(25214903917L, 11, 1L << 48, 16, seed);
  268. }
  269.  
  270. EngineBase *minstd_rand_new(unsigned long seed)
  271. {
  272.     // 31-bit on the closed interval [1, 2^31-2].
  273.     if (seed % 2147483647L == 0)
  274.         seed = 1;
  275.     return linearCongruentialEngine_new(48271L, 0, 2147483647L, 0, seed);
  276. }
  277.  
  278. EngineBase *ranlux24_new(unsigned long seed)
  279. {
  280.     // 24-bit chaotic with long period.
  281.     EngineBase *ranlux24_base = subtractWithCarryEngine_new(24, 10, 24, seed);
  282.     return discardBlockEngine_new(&ranlux24_base, 223, 23);
  283. }
  284.  
  285. EngineBase *ranlux48_new(unsigned long seed)
  286. {
  287.     // 48-bit chaotic with long period.
  288.     EngineBase *ranlux48_base = subtractWithCarryEngine_new(48, 5, 12, seed);
  289.     return discardBlockEngine_new(&ranlux48_base, 389, 11);
  290. }
  291.  
  292. //
  293. // Main.
  294. //
  295.  
  296. double clock_now(void)
  297. {
  298.     struct timespec now;
  299.     clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &now);
  300.     return now.tv_sec + now.tv_nsec / 1.0E+09;
  301. }
  302.  
  303. void show(int n, EngineBase *(*new_)(unsigned long), unsigned long seed)
  304. {
  305.     EngineBase *engine = new_(seed);
  306.  
  307.     for (int i = 0; i < n; i++)
  308.         printf("%ld\n", engine_next(engine));
  309.  
  310.     double t = clock_now();
  311.     engine_discard(engine, 1000000L);
  312.     printf("Elapsed: %.9fs\n", clock_now() - t);
  313.  
  314.     engine_delete(engine);
  315. }
  316.  
  317. int main(void)
  318. {
  319.     unsigned long seed = time(0);
  320.     show(8, default_rand_new, seed);
  321.     show(8, minstd_rand_new, seed);
  322.     show(8, ranlux24_new, seed);
  323.     show(8, ranlux48_new, seed);
  324.     return 0;
  325. }
Success #stdin #stdout 0.16s 5288KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
657968727
1541920992
3842421640
4058446105
420737484
2769223732
2322098806
2843797956
Elapsed: 0.009451811s
520552454
2028837134
153057526
896091866
634845812
50548362
478559110
45208031
Elapsed: 0.009451343s
8269242
15861073
16576166
11246857
14803049
2374360
6144976
2803864
Elapsed: 0.033663982s
162947252066944
67107160573369
98596998364663
32702569952025
109376477829477
232982971992861
243134774958870
260142457476383
Elapsed: 0.107014445s
https://ideone.com/IEOFd0
language:
C (gcc 8.3)
created:
46 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!