fork download
copy 
  1. // Ranlux24 PRNG. (2.02)
  2. // Chaotic with long period but slow due to discard.
  3.  
  4. #include <stdlib.h>
  5. #include <stdint.h>
  6. #include <limits.h>
  7. #include <assert.h>
  8. #include <stdio.h>
  9. #include <time.h>
  10.  
  11. // Utility.
  12.  
  13. #define BITS(x) \
  14.     (sizeof(x) * CHAR_BIT)
  15.  
  16. #define UNLIKELY(cond) \
  17.     __builtin_expect(!!(cond), 0)
  18.  
  19. #define NEW_T(T, n) \
  20.     ((T *) allocate(n * sizeof(T)))
  21.  
  22. void *allocate(size_t n)
  23. {
  24.     void *r = malloc(n);
  25.     assert(r != 0 || n == 0);
  26.     return r;
  27. }
  28.  
  29. // Private.
  30.  
  31. #define SWC_M 0xFFFFFFL
  32. #define SWC_S 10
  33. #define SWC_R 24
  34. #define BLK_P 223
  35. #define BLK_R 23
  36.  
  37. typedef struct {
  38.     int_least32_t *x; int c, i, j;
  39. } Ranlux24;
  40.  
  41. static inline long next(int_least32_t *x, int i_s, int i_r, int *carry)
  42. {
  43.     int_least32_t y = x[i_s] - x[i_r] - *carry;
  44.     *carry = -(y >> (BITS(y) - 1));
  45.     return y & SWC_M;
  46. }
  47.  
  48. static long subtract_with_carry_next(Ranlux24 *self)
  49. {
  50.     self->i += 1;
  51.     int_least32_t *x = self->x;
  52.     int i = self->i, r = SWC_R;
  53.  
  54.     if (UNLIKELY(i == r))
  55.     {
  56.         int c = self->c, s = SWC_S, t = r - s;
  57.  
  58.         for (i = 0; i < s; i++)
  59.             x[i] = next(x, i + t, i, &c);
  60.         for (i = s; i < r; i++)
  61.             x[i] = next(x, i - s, i, &c);
  62.         self->c = c;
  63.         self->i = i = 0;
  64.     }
  65.     return x[i];
  66. }
  67.  
  68. static void subtract_with_carry_discard(Ranlux24 *self, long n)
  69. {
  70.     for (long i = 0; i < n; i++)
  71.         subtract_with_carry_next(self);
  72. }
  73.  
  74. // Public.
  75.  
  76. long ranlux24_next(Ranlux24 *self)
  77. {
  78.     if (UNLIKELY(self->j == 0))
  79.     {
  80.         subtract_with_carry_discard(self, BLK_P - BLK_R);
  81.         self->j = BLK_R;
  82.     }
  83.     self->j -= 1;
  84.     return subtract_with_carry_next(self);
  85. }
  86.  
  87. void ranlux24_discard(Ranlux24 *self, long n)
  88. {
  89.     for (long i = 0; i < n; i++)
  90.         ranlux24_next(self);
  91. }
  92.  
  93. void ranlux24_delete(Ranlux24 *self)
  94. {
  95.     if (self != 0)
  96.     {
  97.         free(self->x);
  98.         free(self);
  99.     }
  100. }
  101.  
  102. Ranlux24 *ranlux24_new(unsigned long seed)
  103. {
  104.     Ranlux24 *self = NEW_T(Ranlux24, 1);
  105.     self->x = NEW_T(int_least32_t, SWC_R);
  106.     self->c = 0;
  107.     self->i = SWC_R-1;
  108.     self->j = BLK_R;
  109.  
  110.     unsigned short rs[3] = {0x330E, seed, seed>>16};
  111.  
  112.     for (int i = 0; i < SWC_R; i++)
  113.         self->x[i] = nrand48(rs) & SWC_M;
  114.     return self;
  115. }
  116.  
  117. // Main.
  118.  
  119. double clock_now(void)
  120. {
  121.     struct timespec now;
  122.     clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &now);
  123.     return now.tv_sec + now.tv_nsec / 1.0E+09;
  124. }
  125.  
  126. int main(void)
  127. {
  128.     unsigned long seed = time(0);
  129.     Ranlux24 *engine = ranlux24_new(seed);
  130.  
  131.     for (int i = 0; i < 24; i++)
  132.     {
  133.         long r = ranlux24_next(engine);
  134.         printf("%ld\n", r);
  135.     }
  136.  
  137.     long n = 10000000;
  138.     double t = clock_now();
  139.     ranlux24_discard(engine, n);
  140.     printf("Elapsed: %.9fs\n", clock_now() - t);
  141.  
  142.     ranlux24_delete(engine);
  143.     return 0;
  144. }
Success #stdin #stdout 0.25s 5288KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
4178616
3776049
7618836
9538652
5567476
15936109
14508867
3342902
5241131
10678486
15286458
12639552
14680109
4575223
11398508
14086026
16668349
8725668
3320411
12219133
9300157
2133998
5268501
11861005
Elapsed: 0.245708888s
https://ideone.com/APtTN3
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!