fork(1) download
  1. // Ranlux24 PRNG. (1.03)
  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)
  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. for (int i = 0; i < r; i++)
  121. self->x[i] = rand() & self->m;
  122. return (EngineBase *) self;
  123. }
  124.  
  125. // Discard Block.
  126.  
  127. typedef struct {
  128. EngineBase base, *owned;
  129. int p, r, i;
  130. } DiscardBlockEngine;
  131.  
  132. int DiscardBlockEngine_next(EngineBase *base)
  133. {
  134. DiscardBlockEngine *self = (DiscardBlockEngine *) base;
  135.  
  136. if (self->i == 0)
  137. {
  138. EngineBase_discard(self->owned, self->p - self->r);
  139. self->i = self->r;
  140. }
  141. self->i -= 1;
  142. return EngineBase_next(self->owned);
  143. }
  144.  
  145. void DiscardBlockEngine_discard(EngineBase *base, int n)
  146. {
  147. for (int i = 0; i < n; i++)
  148. DiscardBlockEngine_next(base);
  149. }
  150.  
  151. void DiscardBlockEngine_delete(EngineBase *base)
  152. {
  153. DiscardBlockEngine *self = (DiscardBlockEngine *) base;
  154.  
  155. if (self != 0)
  156. {
  157. EngineBase_delete(self->owned);
  158. free(self);
  159. }
  160. }
  161.  
  162. EngineBase *DiscardBlockEngine_new(EngineBase **unique, int p, int r)
  163. {
  164. assert(0 < r && r <= p);
  165.  
  166. DiscardBlockEngine *self = NEW_T(DiscardBlockEngine, 1);
  167.  
  168. self->base.next = DiscardBlockEngine_next;
  169. self->base.discard = DiscardBlockEngine_discard;
  170. self->base.delete = DiscardBlockEngine_delete;
  171. self->owned = *unique; *unique = 0; // Transfer ownership.
  172. self->p = p;
  173. self->r = r;
  174. self->i = r;
  175. return (EngineBase *) self;
  176. }
  177.  
  178. // Ranlux24.
  179.  
  180. EngineBase *Ranlux24_new(void)
  181. {
  182. EngineBase *ranlux24_base = SubtractWithCarryEngine_new(24, 10, 24);
  183. return DiscardBlockEngine_new(&ranlux24_base, 223, 23);
  184. }
  185.  
  186. // Main.
  187.  
  188. double clock_now(void)
  189. {
  190. struct timespec now;
  191. clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &now);
  192. return now.tv_sec + now.tv_nsec / 1.0E+09;
  193. }
  194.  
  195. int main(void)
  196. {
  197. srand(time(0));
  198. EngineBase *ranlux24 = Ranlux24_new();
  199.  
  200. for (int i = 0; i < 24; i++)
  201. {
  202. int r = EngineBase_next(ranlux24);
  203. printf("%d\n", r);
  204. }
  205.  
  206. int n = 10000000;
  207. double t = clock_now();
  208. EngineBase_discard(ranlux24, n);
  209. printf("Elapsed: %.9fs\n", clock_now() - t);
  210.  
  211. EngineBase_delete(ranlux24);
  212. return 0;
  213. }
Success #stdin #stdout 0.26s 5284KB
stdin
Standard input is empty
stdout
3352764
5920169
10445798
12617574
3820924
3539462
8716122
5897964
11840758
3846825
4906109
6528382
10633070
4488847
8338194
8792474
9950155
2915225
14011476
6341764
12752839
5359972
11594245
14435949
Elapsed: 0.256149355s