fork download
copy 
  1. // Ranlux24 PRNG. (1.05)
  2. // Chaotic with long period but slow due to discard.
  3.  
  4. #include <stdlib.h>
  5. #include <limits.h>
  6. #include <assert.h>
  7. #include <stdio.h>
  8. #include <time.h>
  9.  
  10. // Utility.
  11.  
  12. #define BITS(x) \
  13.     (sizeof(x) * CHAR_BIT)
  14.  
  15. #define UNLIKELY(cond) \
  16.     __builtin_expect(!!(cond), 0)
  17.  
  18. #define NEW_T(T, n) \
  19.     ((T *) allocate(n * sizeof(T)))
  20.  
  21. void *allocate(size_t n)
  22. {
  23.     void *r = malloc(n);
  24.     assert(r != 0 || n == 0);
  25.     return r;
  26. }
  27.  
  28. // Engine Base.
  29.  
  30. typedef struct EngineBase {
  31.     int  (*next)    (struct EngineBase*);
  32.     void (*discard) (struct EngineBase*, int);
  33.     void (*delete)  (struct EngineBase*);
  34. } EngineBase;
  35.  
  36. int engineBase_next(EngineBase *self)
  37. {
  38.     return self->next(self);
  39. }
  40.  
  41. void engineBase_discard(EngineBase *self, int n)
  42. {
  43.     return self->discard(self, n);
  44. }
  45.  
  46. void engineBase_delete(EngineBase *self)
  47. {
  48.     return self->delete(self);
  49. }
  50.  
  51. // Subtract With Carry.
  52.  
  53. typedef struct {
  54.     EngineBase base;
  55.     int *x, c, i, m, r, s;
  56. } SubtractWithCarryEngine;
  57.  
  58. static inline int next(int *x, int mask, int i_s, int i_r, int *carry)
  59. {
  60.     int y = x[i_s] - x[i_r] - *carry;
  61.     *carry = -(y >> (BITS(int) - 1));
  62.     return y & mask;
  63. }
  64.  
  65. int subtractWithCarryEngine_next(EngineBase *base)
  66. {
  67.     SubtractWithCarryEngine *self = (SubtractWithCarryEngine *) base;
  68.  
  69.     self->i += 1;
  70.     int *x = self->x, i = self->i, r = self->r;
  71.  
  72.     if (UNLIKELY(i >= r))
  73.     {
  74.         int c = self->c, m = self->m, s = self->s, t = r - s;
  75.  
  76.         for (i = 0; i < s; i++)
  77.             x[i] = next(x, m, i + t, i, &c);
  78.         for (i = s; i < r; i++)
  79.             x[i] = next(x, m, i - s, i, &c);
  80.         self->c = c;
  81.         self->i = i = 0;
  82.     }
  83.     return x[i];
  84. }
  85.  
  86. void subtractWithCarryEngine_discard(EngineBase *base, int n)
  87. {
  88.     for (int i = 0; i < n; i++)
  89.         subtractWithCarryEngine_next(base);
  90. }
  91.  
  92. void subtractWithCarryEngine_delete(EngineBase *base)
  93. {
  94.     SubtractWithCarryEngine *self = (SubtractWithCarryEngine *) base;
  95.  
  96.     if (self != 0)
  97.     {
  98.         free(self->x);
  99.         free(self);
  100.     }
  101. }
  102.  
  103. EngineBase *subtractWithCarryEngine_new(int w, int s, int r, unsigned int seed)
  104. {
  105.     assert(0 < w && w < BITS(int));
  106.     assert(0 < s && s < r);
  107.  
  108.     SubtractWithCarryEngine *self = NEW_T(SubtractWithCarryEngine, 1);
  109.  
  110.     self->base.next = subtractWithCarryEngine_next;
  111.     self->base.discard = subtractWithCarryEngine_discard;
  112.     self->base.delete = subtractWithCarryEngine_delete;
  113.     self->x = NEW_T(int, r);
  114.     self->c = 0;
  115.     self->i = r-1;
  116.     self->m = -1U >> (BITS(int) - w);
  117.     self->r = r;
  118.     self->s = s;
  119.  
  120.     unsigned short state[] = {0x330E, seed, seed>>16};
  121.     for (int i = 0; i < r; i++)
  122.         self->x[i] = nrand48(state) & self->m;
  123.     return (EngineBase *) self;
  124. }
  125.  
  126. // Discard Block.
  127.  
  128. typedef struct {
  129.     EngineBase base, *owned;
  130.     int p, r, i;
  131. } DiscardBlockEngine;
  132.  
  133. int discardBlockEngine_next(EngineBase *base)
  134. {
  135.     DiscardBlockEngine *self = (DiscardBlockEngine *) base;
  136.  
  137.     if (self->i == 0)
  138.     {
  139.         engineBase_discard(self->owned, self->p - self->r);
  140.         self->i = self->r;
  141.     }
  142.     self->i -= 1;
  143.     return engineBase_next(self->owned);
  144. }
  145.  
  146. void discardBlockEngine_discard(EngineBase *base, int n)
  147. {
  148.     for (int i = 0; i < n; i++)
  149.         discardBlockEngine_next(base);
  150. }
  151.  
  152. void discardBlockEngine_delete(EngineBase *base)
  153. {
  154.     DiscardBlockEngine *self = (DiscardBlockEngine *) base;
  155.  
  156.     if (self != 0)
  157.     {
  158.         engineBase_delete(self->owned);
  159.         free(self);
  160.     }
  161. }
  162.  
  163. EngineBase *discardBlockEngine_new(EngineBase **unique, int p, int r)
  164. {
  165.     assert(0 < r && r <= p);
  166.  
  167.     DiscardBlockEngine *self = NEW_T(DiscardBlockEngine, 1);
  168.  
  169.     self->base.next = discardBlockEngine_next;
  170.     self->base.discard = discardBlockEngine_discard;
  171.     self->base.delete = discardBlockEngine_delete;
  172.     self->owned = *unique; *unique = 0; // Transfer ownership.
  173.     self->p = p;
  174.     self->r = r;
  175.     self->i = r;
  176.     return (EngineBase *) self;
  177. }
  178.  
  179. // Ranlux24.
  180.  
  181. EngineBase *ranlux24_new(unsigned int seed)
  182. {
  183.     EngineBase *ranlux24_base = subtractWithCarryEngine_new(24, 10, 24, seed);
  184.     return discardBlockEngine_new(&ranlux24_base, 223, 23);
  185. }
  186.  
  187. // Main.
  188.  
  189. double clock_now(void)
  190. {
  191.     struct timespec now;
  192.     clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &now);
  193.     return now.tv_sec + now.tv_nsec / 1.0E+09;
  194. }
  195.  
  196. int main(void)
  197. {
  198.     unsigned int seed = time(0);
  199.     EngineBase *engine = ranlux24_new(seed);
  200.  
  201.     for (int i = 0; i < 24; i++)
  202.     {
  203.         int r = engineBase_next(engine);
  204.         printf("%d\n", r);
  205.     }
  206.  
  207.     int n = 10000000;
  208.     double t = clock_now();
  209.     engineBase_discard(engine, n);
  210.     printf("Elapsed: %.9fs\n", clock_now() - t);
  211.  
  212.     engineBase_delete(engine);
  213.     return 0;
  214. }
Success #stdin #stdout 0.33s 5272KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
12440664
3989329
14269365
7504508
12484499
13983117
10217443
16333398
10657482
2176501
6693054
14669557
5769682
3883724
13257135
8343678
12244641
14577848
132831
13856560
14070853
5990614
2043614
2205617
Elapsed: 0.322357405s
https://ideone.com/JppEU1
language:
C (gcc 8.3)
created:
10 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!