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