src/rng.c

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