Back to index

plt-scheme  4.2.1
newrandom.inc
Go to the documentation of this file.
00001 /*
00002   Based on
00003 
00004    Implementation of SRFI-27 core generator in C for MzScheme.
00005    dvanhorn@cs.uvm.edu
00006 
00007   and
00008 
00009    54-BIT (double) IMPLEMENTATION IN C OF THE "MRG32K3A" GENERATOR
00010    ===============================================================
00011 
00012    Sebastian.Egner@philips.com, Mar-2002, in ANSI-C and Scheme 48 0.57
00013 
00014    This code is a C-implementation of Pierre L'Ecuyer's MRG32k3a generator.
00015    The code uses (double)-arithmetics, assuming that it covers the range
00016    {-2^53..2^53-1} exactly (!). The code of the generator is based on the
00017    L'Ecuyer's own implementation of the generator. Please refer to the
00018    file 'mrg32k3a.scm' for more information about the method.
00019 */
00020 
00021 /* maximum value for random_integer: min(S48_MAX_FIXNUM_VALUE, m1) */
00022 #define m_max (((long)1 << 31) - 1)
00023 
00024 /* The Generator
00025    =============
00026 */
00027 
00028 /* moduli of the components */
00029 #define Im1 0xffffff2f
00030 #define Im2 0xffffa6bb
00031 #define m1 4294967087.0
00032 #define m2 4294944443.0
00033 
00034 /* recursion coefficients of the components */
00035 #define a12  1403580.0
00036 #define a13n  810728.0
00037 #define a21   527612.0
00038 #define a23n 1370589.0
00039 
00040 /* normalization factor 1/(m1 + 1) */
00041 #define norm 2.328306549295728e-10
00042 
00043 /* the actual generator */
00044 
00045 static double mrg32k3a(Scheme_Random_State *s) { /* (double), in {0..m1-1} */
00046   double x10, x20, y;
00047   long   k10, k20;
00048 
00049   /* component 1 */
00050   x10  = a12*(s->x11) - a13n*(s->x12);
00051   k10  = (long)(x10 / m1);
00052   x10 -= k10 * m1;
00053   if (x10 < 0.0)
00054     x10 += m1;
00055   s->x12 = s->x11;
00056   s->x11 = s->x10;
00057   s->x10 = x10;
00058 
00059   /* component 2 */
00060   x20  = a21*(s->x20) - a23n*(s->x22);
00061   k20  = (long)(x20 / m2);
00062   x20 -= k20 * m2;
00063   if (x20 < 0.0)
00064     x20 += m2;
00065   s->x22 = s->x21;
00066   s->x21 = s->x20;
00067   s->x20 = x20;
00068 
00069   /* combination of component */
00070   y = x10 - x20;
00071   if (y < 0.0)
00072     y += m1;
00073   return y;
00074 }
00075 
00076 /**************************************************/
00077 
00078 static Scheme_Object *pack_rand_state(Scheme_Object *vec, Scheme_Random_State *s)
00079 {
00080   if (!s) {
00081     s = (Scheme_Random_State *)scheme_malloc_atomic_tagged(sizeof(Scheme_Random_State));
00082     s->so.type = scheme_random_state_type;
00083 
00084   }
00085 #define REF(r, i, top) \
00086   { \
00087     unsigned long l; \
00088     if (!scheme_get_unsigned_int_val(SCHEME_VEC_ELS(vec)[i],  &l)) \
00089       return NULL; \
00090     if (l > top - 1) \
00091       return NULL; \
00092     r = (double)l; \
00093   }
00094 
00095   /* copy the numbers from state into s */
00096   REF(s->x10, 0, Im1)
00097   REF(s->x11, 1, Im1)
00098   REF(s->x12, 2, Im1)
00099   REF(s->x20, 3, Im2)
00100   REF(s->x21, 4, Im2)
00101   REF(s->x22, 5, Im2)
00102 
00103 #undef REF
00104 
00105   if (!s->x10 && !s->x11 && !s->x12)
00106     return NULL;
00107   if (!s->x20 && !s->x21 && !s->x22)
00108     return NULL;
00109 
00110   return (Scheme_Object *)s;
00111 }
00112 
00113 static Scheme_Object *unpack_rand_state(Scheme_Random_State *s)
00114 {
00115   Scheme_Object *result;
00116 
00117   /* make and fill a Scheme vector with the numbers */
00118   result = scheme_make_vector((long)6, NULL);
00119 
00120 #define SET(i, x) \
00121   { \
00122     Scheme_Object *o; \
00123     o = scheme_make_integer_value_from_unsigned((unsigned long)x); \
00124     SCHEME_VEC_ELS(result)[i] = o; \
00125   }
00126 
00127   SET(0, s->x10)
00128   SET(1, s->x11)
00129   SET(2, s->x12)
00130   SET(3, s->x20)
00131   SET(4, s->x21)
00132   SET(5, s->x22)
00133 
00134 #undef SET
00135 
00136   return result;
00137 }
00138 
00139 /**************************************************/
00140 
00141 static unsigned int _random_m(unsigned int *_x)
00142 {
00143   unsigned int x, y;
00144   x = *_x;
00145   y = x & 0xFFFF;
00146   x = (30903 * y) + (x >> 16);
00147   *_x = x;
00148   return y;
00149 }
00150 
00151 static int _random_n(unsigned int *_x, int n)
00152 {
00153   return ((_random_m(_x) << 16) + _random_m(_x)) % n;
00154 }
00155 
00156 static void sch_srand_half(unsigned int x, Scheme_Random_State *s)
00157 {
00158   /* Due to integer overflow, this doesn't match the Scheme implementation!
00159      We use "int" instead of "long" to make the overflow consistent
00160      across platforms (since "long" is sometimes 64 bits). */
00161   unsigned int z;
00162   z = _random_n(&x, Im1-1);
00163   s->x10 = (double)(1 + (((unsigned int)s->x10 + z) % (Im1 - 1)));
00164   z = _random_n(&x, Im1);
00165   s->x11 = (double)(((unsigned int)s->x11 + z) % Im1);
00166   z = _random_n(&x, Im1);
00167   s->x12 = (double)(((unsigned int)s->x12 + z) % Im1);
00168   z = _random_n(&x, Im2-1);
00169   s->x20 = (double)(1 + (((unsigned int)s->x20 + z) % (Im2 - 1)));
00170   z = _random_n(&x, Im2);
00171   s->x21 = (double)(((unsigned int)s->x21 + z) % Im2);
00172   z = _random_n(&x, Im2);
00173   s->x22 = (double)(((unsigned int)s->x22 + z) % Im2);
00174 
00175   /* Due to the mismatch, maybe it's possible that we can hit a degeneracy?
00176      Double-check, just in case... */
00177   if (!s->x10 && !s->x11 && !s->x12)
00178     s->x10 = 1;
00179   if (!s->x20 && !s->x21 && !s->x22)
00180     s->x20 = 1;
00181 }
00182 
00183 static void sch_srand(unsigned int x, Scheme_Random_State *s)
00184 {
00185   /* Initial values are from Sebastian Egner's implementation: */
00186   s->x10 = 1062452522.0;
00187   s->x11 = 2961816100.0;
00188   s->x12 = 342112271.0;
00189   s->x20 = 2854655037.0;
00190   s->x21 = 3321940838.0;
00191   s->x22 = 3542344109.0;
00192 
00193   sch_srand_half(x & 0xFFFF, s);
00194   sch_srand_half((x >> 16) & 0xFFFF, s);
00195 }
00196 
00197 /**************************************************/
00198 
00199 static unsigned long sch_int_rand(unsigned long n, Scheme_Random_State *rs)
00200 {
00201   double  x, q, qn, xq;
00202 
00203   /* generate result in {0..n-1} using the rejection method */
00204   q  = (double)( (unsigned long)(m1 / (double)n) );
00205   qn = q * n;
00206   do {
00207     x = mrg32k3a(rs);
00208   } while (x >= qn);
00209   xq = x / q;
00210 
00211   /* return result */
00212   return (unsigned long)xq;
00213 }
00214 
00215 static double sch_double_rand(Scheme_Random_State *rs)
00216 {
00217   double  x;
00218   x = mrg32k3a(rs);
00219   return (x + 1.0) * norm;
00220 }
00221 
00222 Scheme_Object *scheme_make_random_state(long seed)
00223 {
00224   Scheme_Random_State *s;
00225 
00226   s = (Scheme_Random_State *)scheme_malloc_atomic_tagged(sizeof(Scheme_Random_State));
00227   s->so.type = scheme_random_state_type;
00228 
00229   sch_srand(seed, s);
00230 
00231   return (Scheme_Object *)s;
00232 }