1/*
2  Based on
3
4   Implementation of SRFI-27 core generator in C for Racket.
5   dvanhorn@cs.uvm.edu
6
7  and
8
9   54-BIT (double) IMPLEMENTATION IN C OF THE "MRG32K3A" GENERATOR
10   ===============================================================
11
12   Sebastian.Egner@philips.com, Mar-2002, in ANSI-C and Scheme 48 0.57
13
14   This code is a C-implementation of Pierre L'Ecuyer's MRG32k3a generator.
15   The code uses (double)-arithmetics, assuming that it covers the range
16   {-2^53..2^53-1} exactly (!). The code of the generator is based on the
17   L'Ecuyer's own implementation of the generator. Please refer to the
18   file 'mrg32k3a.scm' for more information about the method.
19*/
20
21/* maximum value for random_integer: min(S48_MAX_FIXNUM_VALUE, m1) */
22#define m_max (((intptr_t)1 << 31) - 1)
23
24/* The Generator
25   =============
26*/
27
28/* moduli of the components */
29#define Im1 0xffffff2f
30#define Im2 0xffffa6bb
31#define m1 4294967087.0
32#define m2 4294944443.0
33
34/* recursion coefficients of the components */
35#define a12  1403580.0
36#define a13n  810728.0
37#define a21   527612.0
38#define a23n 1370589.0
39
40/* normalization factor 1/(m1 + 1) */
41#define norm 2.328306549295728e-10
42
43/* the actual generator */
44
45XFORM_NONGCING static double mrg32k3a(Scheme_Random_State *s) { /* (double), in {0..m1-1} */
46  double x10, x20, y;
47  intptr_t   k10, k20;
48
49  /* component 1 */
50  x10  = a12*(s->x11) - a13n*(s->x12);
51  k10  = (intptr_t)(x10 / m1);
52  x10 -= k10 * m1;
53  if (x10 < 0.0)
54    x10 += m1;
55  s->x12 = s->x11;
56  s->x11 = s->x10;
57  s->x10 = x10;
58
59  /* component 2 */
60  x20  = a21*(s->x20) - a23n*(s->x22);
61  k20  = (intptr_t)(x20 / m2);
62  x20 -= k20 * m2;
63  if (x20 < 0.0)
64    x20 += m2;
65  s->x22 = s->x21;
66  s->x21 = s->x20;
67  s->x20 = x20;
68
69  /* combination of component */
70  y = x10 - x20;
71  if (y < 0.0)
72    y += m1;
73  return y;
74}
75
76/**************************************************/
77
78static Scheme_Object *pack_rand_state(Scheme_Object *vec, Scheme_Random_State *s)
79{
80  if (!s) {
81    s = (Scheme_Random_State *)scheme_malloc_atomic_tagged(sizeof(Scheme_Random_State));
82    s->so.type = scheme_random_state_type;
83
84  }
85#define REF(r, i, top) \
86  { \
87    uintptr_t l; \
88    if (!scheme_get_unsigned_int_val(scheme_chaperone_vector_ref(vec, i),  &l)) \
89      return NULL; \
90    if (l > top - 1) \
91      return NULL; \
92    r = (double)l; \
93  }
94
95  /* copy the numbers from state into s */
96  REF(s->x10, 0, Im1)
97  REF(s->x11, 1, Im1)
98  REF(s->x12, 2, Im1)
99  REF(s->x20, 3, Im2)
100  REF(s->x21, 4, Im2)
101  REF(s->x22, 5, Im2)
102
103#undef REF
104
105  if (!s->x10 && !s->x11 && !s->x12)
106    return NULL;
107  if (!s->x20 && !s->x21 && !s->x22)
108    return NULL;
109
110  return (Scheme_Object *)s;
111}
112
113static Scheme_Object *unpack_rand_state(Scheme_Random_State *s)
114{
115  Scheme_Object *result;
116
117  /* make and fill a Scheme vector with the numbers */
118  result = scheme_make_vector((intptr_t)6, NULL);
119
120#define SET(i, x) \
121  { \
122    Scheme_Object *o; \
123    o = scheme_make_integer_value_from_unsigned((uintptr_t)x); \
124    SCHEME_VEC_ELS(result)[i] = o; \
125  }
126
127  SET(0, s->x10)
128  SET(1, s->x11)
129  SET(2, s->x12)
130  SET(3, s->x20)
131  SET(4, s->x21)
132  SET(5, s->x22)
133
134#undef SET
135
136  return result;
137}
138
139/**************************************************/
140
141static unsigned int _random_m(unsigned int *_x)
142{
143  unsigned int x, y;
144  x = *_x;
145  y = x & 0xFFFF;
146  x = (30903 * y) + (x >> 16);
147  *_x = x;
148  return y;
149}
150
151static int _random_n(unsigned int *_x, int n)
152{
153  return ((_random_m(_x) << 16) + _random_m(_x)) % n;
154}
155
156static void sch_srand_half(unsigned int x, Scheme_Random_State *s)
157{
158  /* Due to integer overflow, this doesn't match the Scheme implementation!
159     We use "int" instead of "long" to make the overflow consistent
160     across platforms (since "long" is sometimes 64 bits). */
161  unsigned int z;
162  z = _random_n(&x, Im1-1);
163  s->x10 = (double)(1 + (((unsigned int)s->x10 + z) % (Im1 - 1)));
164  z = _random_n(&x, Im1);
165  s->x11 = (double)(((unsigned int)s->x11 + z) % Im1);
166  z = _random_n(&x, Im1);
167  s->x12 = (double)(((unsigned int)s->x12 + z) % Im1);
168  z = _random_n(&x, Im2-1);
169  s->x20 = (double)(1 + (((unsigned int)s->x20 + z) % (Im2 - 1)));
170  z = _random_n(&x, Im2);
171  s->x21 = (double)(((unsigned int)s->x21 + z) % Im2);
172  z = _random_n(&x, Im2);
173  s->x22 = (double)(((unsigned int)s->x22 + z) % Im2);
174
175  /* Due to the mismatch, maybe it's possible that we can hit a degeneracy?
176     Double-check, just in case... */
177  if (!s->x10 && !s->x11 && !s->x12)
178    s->x10 = 1;
179  if (!s->x20 && !s->x21 && !s->x22)
180    s->x20 = 1;
181}
182
183static void sch_srand(unsigned int x, Scheme_Random_State *s)
184{
185  /* Initial values are from Sebastian Egner's implementation: */
186  s->x10 = 1062452522.0;
187  s->x11 = 2961816100.0;
188  s->x12 = 342112271.0;
189  s->x20 = 2854655037.0;
190  s->x21 = 3321940838.0;
191  s->x22 = 3542344109.0;
192
193  sch_srand_half(x & 0xFFFF, s);
194  sch_srand_half((x >> 16) & 0xFFFF, s);
195}
196
197/**************************************************/
198
199static uintptr_t sch_int_rand(uintptr_t n, Scheme_Random_State *rs)
200{
201  double  x, q, qn, xq;
202
203  /* generate result in {0..n-1} using the rejection method */
204  q  = (double)( (uintptr_t)(m1 / (double)n) );
205  qn = q * n;
206  do {
207    x = mrg32k3a(rs);
208  } while (x >= qn);
209  xq = x / q;
210
211  /* return result */
212  return (uintptr_t)xq;
213}
214
215XFORM_NONGCING static double sch_double_rand(Scheme_Random_State *rs)
216{
217  double  x;
218  x = mrg32k3a(rs);
219  return (x + 1.0) * norm;
220}
221
222Scheme_Object *scheme_make_random_state(intptr_t seed)
223{
224  Scheme_Random_State *s;
225
226  s = (Scheme_Random_State *)scheme_malloc_atomic_tagged(sizeof(Scheme_Random_State));
227  s->so.type = scheme_random_state_type;
228
229  sch_srand(seed, s);
230
231  return (Scheme_Object *)s;
232}
233