htswanalysis/MACS/lib/gsl/gsl-1.11/rng/mt.c

   1 /* This program is free software; you can redistribute it and/or
   2    modify it under the terms of the GNU General Public License as
   3    published by the Free Software Foundation; either version 3 of the
   4    License, or (at your option) any later version.
   5
   6    This program is distributed in the hope that it will be useful, but
   7    WITHOUT ANY WARRANTY; without even the implied warranty of
   8    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   9    General Public License for more details.  You should have received
  10    a copy of the GNU General Public License along with this program;
  11    if not, write to the Free Foundation, Inc., 59 Temple Place, Suite
  12    330, Boston, MA 02111-1307 USA
  13
  14    Original implementation was copyright (C) 1997 Makoto Matsumoto and
  15    Takuji Nishimura. Coded by Takuji Nishimura, considering the
  16    suggestions by Topher Cooper and Marc Rieffel in July-Aug. 1997, "A
  17    C-program for MT19937: Integer version (1998/4/6)"
  18
  19    This implementation copyright (C) 1998 Brian Gough. I reorganized
  20    the code to use the module framework of GSL.  The license on this
  21    implementation was changed from LGPL to GPL, following paragraph 3
  22    of the LGPL, version 2.
  23
  24    Update:
  25
  26    The seeding procedure has been updated to match the 10/99 release
  27    of MT19937.
  28
  29    Update:
  30
  31    The seeding procedure has been updated again to match the 2002
  32    release of MT19937
  33
  34    The original code included the comment: "When you use this, send an
  35    email to: matumoto@math.keio.ac.jp with an appropriate reference to
  36    your work".
  37
  38    Makoto Matsumoto has a web page with more information about the
  39    generator, http://www.math.keio.ac.jp/~matumoto/emt.html.
  40
  41    The paper below has details of the algorithm.
  42
  43    From: Makoto Matsumoto and Takuji Nishimura, "Mersenne Twister: A
  44    623-dimensionally equidistributerd uniform pseudorandom number
  45    generator". ACM Transactions on Modeling and Computer Simulation,
  46    Vol. 8, No. 1 (Jan. 1998), Pages 3-30
  47
  48    You can obtain the paper directly from Makoto Matsumoto's web page.
  49
  50    The period of this generator is 2^{19937} - 1.
  51
  52 */
  53
  54 #include <config.h>
  55 #include <stdlib.h>
  56 #include <gsl/gsl_rng.h>
  57
  58 static inline unsigned long int mt_get (void *vstate);
  59 static double mt_get_double (void *vstate);
  60 static void mt_set (void *state, unsigned long int s);
  61
  62 #define N 624   /* Period parameters */
  63 #define M 397
  64
  65 /* most significant w-r bits */
  66 static const unsigned long UPPER_MASK = 0x80000000UL;
  67
  68 /* least significant r bits */
  69 static const unsigned long LOWER_MASK = 0x7fffffffUL;
  70
  71 typedef struct
  72   {
  73     unsigned long mt[N];
  74     int mti;
  75   }
  76 mt_state_t;
  77
  78 static inline unsigned long
  79 mt_get (void *vstate)
  80 {
  81   mt_state_t *state = (mt_state_t *) vstate;
  82
  83   unsigned long k ;
  84   unsigned long int *const mt = state->mt;
  85
  86 #define MAGIC(y) (((y)&0x1) ? 0x9908b0dfUL : 0)
  87
  88   if (state->mti >= N)
  89     {   /* generate N words at one time */
  90       int kk;
  91
  92       for (kk = 0; kk < N - M; kk++)
  93         {
  94           unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
  95           mt[kk] = mt[kk + M] ^ (y >> 1) ^ MAGIC(y);
  96         }
  97       for (; kk < N - 1; kk++)
  98         {
  99           unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
 100           mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ MAGIC(y);
 101         }
 102
 103       {
 104         unsigned long y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
 105         mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ MAGIC(y);
 106       }
 107
 108       state->mti = 0;
 109     }
 110
 111   /* Tempering */
 112
 113   k = mt[state->mti];
 114   k ^= (k >> 11);
 115   k ^= (k << 7) & 0x9d2c5680UL;
 116   k ^= (k << 15) & 0xefc60000UL;
 117   k ^= (k >> 18);
 118
 119   state->mti++;
 120
 121   return k;
 122 }
 123
 124 static double
 125 mt_get_double (void * vstate)
 126 {
 127   return mt_get (vstate) / 4294967296.0 ;
 128 }
 129
 130 static void
 131 mt_set (void *vstate, unsigned long int s)
 132 {
 133   mt_state_t *state = (mt_state_t *) vstate;
 134   int i;
 135
 136   if (s == 0)
 137     s = 4357;   /* the default seed is 4357 */
 138
 139   state->mt[0]= s & 0xffffffffUL;
 140
 141   for (i = 1; i < N; i++)
 142     {
 143       /* See Knuth's "Art of Computer Programming" Vol. 2, 3rd
 144          Ed. p.106 for multiplier. */
 145
 146       state->mt[i] =
 147         (1812433253UL * (state->mt[i-1] ^ (state->mt[i-1] >> 30)) + i);
 148
 149       state->mt[i] &= 0xffffffffUL;
 150     }
 151
 152   state->mti = i;
 153 }
 154
 155 static void
 156 mt_1999_set (void *vstate, unsigned long int s)
 157 {
 158   mt_state_t *state = (mt_state_t *) vstate;
 159   int i;
 160
 161   if (s == 0)
 162     s = 4357;   /* the default seed is 4357 */
 163
 164   /* This is the October 1999 version of the seeding procedure. It
 165      was updated by the original developers to avoid the periodicity
 166      in the simple congruence originally used.
 167
 168      Note that an ANSI-C unsigned long integer arithmetic is
 169      automatically modulo 2^32 (or a higher power of two), so we can
 170      safely ignore overflow. */
 171
 172 #define LCG(x) ((69069 * x) + 1) &0xffffffffUL
 173
 174   for (i = 0; i < N; i++)
 175     {
 176       state->mt[i] = s & 0xffff0000UL;
 177       s = LCG(s);
 178       state->mt[i] |= (s &0xffff0000UL) >> 16;
 179       s = LCG(s);
 180     }
 181
 182   state->mti = i;
 183 }
 184
 185 /* This is the original version of the seeding procedure, no longer
 186    used but available for compatibility with the original MT19937. */
 187
 188 static void
 189 mt_1998_set (void *vstate, unsigned long int s)
 190 {
 191   mt_state_t *state = (mt_state_t *) vstate;
 192   int i;
 193
 194   if (s == 0)
 195     s = 4357;   /* the default seed is 4357 */
 196
 197   state->mt[0] = s & 0xffffffffUL;
 198
 199 #define LCG1998(n) ((69069 * n) & 0xffffffffUL)
 200
 201   for (i = 1; i < N; i++)
 202     state->mt[i] = LCG1998 (state->mt[i - 1]);
 203
 204   state->mti = i;
 205 }
 206
 207 static const gsl_rng_type mt_type =
 208 {"mt19937",                     /* name */
 209  0xffffffffUL,                  /* RAND_MAX  */
 210  0,                             /* RAND_MIN  */
 211  sizeof (mt_state_t),
 212  &mt_set,
 213  &mt_get,
 214  &mt_get_double};
 215
 216 static const gsl_rng_type mt_1999_type =
 217 {"mt19937_1999",                /* name */
 218  0xffffffffUL,                  /* RAND_MAX  */
 219  0,                             /* RAND_MIN  */
 220  sizeof (mt_state_t),
 221  &mt_1999_set,
 222  &mt_get,
 223  &mt_get_double};
 224
 225 static const gsl_rng_type mt_1998_type =
 226 {"mt19937_1998",                /* name */
 227  0xffffffffUL,                  /* RAND_MAX  */
 228  0,                             /* RAND_MIN  */
 229  sizeof (mt_state_t),
 230  &mt_1998_set,
 231  &mt_get,
 232  &mt_get_double};
 233
 234 const gsl_rng_type *gsl_rng_mt19937 = &mt_type;
 235 const gsl_rng_type *gsl_rng_mt19937_1999 = &mt_1999_type;
 236 const gsl_rng_type *gsl_rng_mt19937_1998 = &mt_1998_type;
 237
 238 /* MT19937 is the default generator, so define that here too */
 239
 240 const gsl_rng_type *gsl_rng_default = &mt_type;
 241 unsigned long int gsl_rng_default_seed = 0;