fork download
copy 
  1. // Ranlux24 PRNG. (1.06)
  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.     long (*next)    (struct EngineBase *);
  32.     void (*discard) (struct EngineBase *, long);
  33.     void (*delete)  (struct EngineBase *);
  34. } EngineBase;
  35.  
  36. long engineBase_next(EngineBase *self)
  37. {
  38.     return self->next(self);
  39. }
  40.  
  41. void engineBase_discard(EngineBase *self, long 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.     long *x, m; int c, i, r, s;
  56. } SubtractWithCarryEngine;
  57.  
  58. static inline long next(long *x, long mask, int i_s, int i_r, int *carry)
  59. {
  60.     long y = x[i_s] - x[i_r] - *carry;
  61.     *carry = -(y >> (BITS(y) - 1));
  62.     return y & mask;
  63. }
  64.  
  65. long subtractWithCarryEngine_next(EngineBase *base)
  66. {
  67.     SubtractWithCarryEngine *self = (SubtractWithCarryEngine *) base;
  68.  
  69.     self->i += 1;
  70.     long *x = self->x;
  71.     int   i = self->i, r = self->r;
  72.  
  73.     if (UNLIKELY(i == r))
  74.     {
  75.         long m = self->m;
  76.         int  c = self->c, s = self->s, t = r - s;
  77.  
  78.         for (i = 0; i < s; i++)
  79.             x[i] = next(x, m, i + t, i, &c);
  80.         for (i = s; i < r; i++)
  81.             x[i] = next(x, m, i - s, i, &c);
  82.         self->c = c;
  83.         self->i = i = 0;
  84.     }
  85.     return x[i];
  86. }
  87.  
  88. void subtractWithCarryEngine_discard(EngineBase *base, long n)
  89. {
  90.     for (long i = 0; i < n; i++)
  91.         subtractWithCarryEngine_next(base);
  92. }
  93.  
  94. void subtractWithCarryEngine_delete(EngineBase *base)
  95. {
  96.     SubtractWithCarryEngine *self = (SubtractWithCarryEngine *) base;
  97.  
  98.     if (self != 0)
  99.     {
  100.         free(self->x);
  101.         free(self);
  102.     }
  103. }
  104.  
  105. EngineBase *subtractWithCarryEngine_new(int w, int s, int r, unsigned long seed)
  106. {
  107.     assert(0 < w && w < BITS(long));
  108.     assert(0 < s && s < r);
  109.  
  110.     SubtractWithCarryEngine *self = NEW_T(SubtractWithCarryEngine, 1);
  111.  
  112.     self->base.next = subtractWithCarryEngine_next;
  113.     self->base.discard = subtractWithCarryEngine_discard;
  114.     self->base.delete = subtractWithCarryEngine_delete;
  115.     self->x = NEW_T(long, r);
  116.     self->m = -1UL >> (BITS(long) - w);
  117.     self->c = 0;
  118.     self->i = r-1;
  119.     self->r = r;
  120.     self->s = s;
  121.  
  122.     unsigned short state[3] = {0x330E, seed, seed>>16};
  123.  
  124.     for (int i = 0; i < r; i++)
  125.         self->x[i] = nrand48(state) & self->m;
  126.     return (EngineBase *) self;
  127. }
  128.  
  129. // Discard Block.
  130.  
  131. typedef struct {
  132.     EngineBase base, *owned;
  133.     int p, r, i;
  134. } DiscardBlockEngine;
  135.  
  136. long discardBlockEngine_next(EngineBase *base)
  137. {
  138.     DiscardBlockEngine *self = (DiscardBlockEngine *) base;
  139.  
  140.     if (self->i == 0)
  141.     {
  142.         engineBase_discard(self->owned, self->p - self->r);
  143.         self->i = self->r;
  144.     }
  145.     self->i -= 1;
  146.     return engineBase_next(self->owned);
  147. }
  148.  
  149. void discardBlockEngine_discard(EngineBase *base, long n)
  150. {
  151.     for (long i = 0; i < n; i++)
  152.         discardBlockEngine_next(base);
  153. }
  154.  
  155. void discardBlockEngine_delete(EngineBase *base)
  156. {
  157.     DiscardBlockEngine *self = (DiscardBlockEngine *) base;
  158.  
  159.     if (self != 0)
  160.     {
  161.         engineBase_delete(self->owned);
  162.         free(self);
  163.     }
  164. }
  165.  
  166. EngineBase *discardBlockEngine_new(EngineBase **unique, int p, int r)
  167. {
  168.     assert(0 < r && r <= p);
  169.  
  170.     DiscardBlockEngine *self = NEW_T(DiscardBlockEngine, 1);
  171.  
  172.     self->base.next = discardBlockEngine_next;
  173.     self->base.discard = discardBlockEngine_discard;
  174.     self->base.delete = discardBlockEngine_delete;
  175.     self->owned = *unique; *unique = 0; // Transfer ownership.
  176.     self->p = p;
  177.     self->r = r;
  178.     self->i = r;
  179.     return (EngineBase *) self;
  180. }
  181.  
  182. // Ranlux24.
  183.  
  184. EngineBase *ranlux24_new(unsigned long seed)
  185. {
  186.     EngineBase *ranlux24_base = subtractWithCarryEngine_new(24, 10, 24, seed);
  187.     return discardBlockEngine_new(&ranlux24_base, 223, 23);
  188. }
  189.  
  190. // Main.
  191.  
  192. double clock_now(void)
  193. {
  194.     struct timespec now;
  195.     clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &now);
  196.     return now.tv_sec + now.tv_nsec / 1.0E+09;
  197. }
  198.  
  199. int main(void)
  200. {
  201.     unsigned long seed = time(0);
  202.     EngineBase *engine = ranlux24_new(seed);
  203.  
  204.     for (int i = 0; i < 24; i++)
  205.     {
  206.         long r = engineBase_next(engine);
  207.         printf("%ld\n", r);
  208.     }
  209.  
  210.     long n = 10000000;
  211.     double t = clock_now();
  212.     engineBase_discard(engine, n);
  213.     printf("Elapsed: %.9fs\n", clock_now() - t);
  214.  
  215.     engineBase_delete(engine);
  216.     return 0;
  217. }
Success #stdin #stdout 0.35s 5288KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
14462952
6536096
1403332
16645964
361507
15445726
15361779
5092390
13327451
8997957
8354423
12592934
16718619
1973069
12403977
12530129
16085643
10424154
14884313
13030305
13295224
2889977
6981169
14555912
Elapsed: 0.344467723s
https://ideone.com/GLuACu
language:
C (gcc 8.3)
created:
9 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!