1 files changed, 241 insertions, 0 deletions
diff --git a/gsl-1.9/rng/mt.c b/gsl-1.9/rng/mt.c
new file mode 100644
index 0000000..46de635
--- /dev/null
+++ b/gsl-1.9/rng/mt.c
@@ -0,0 +1,241 @@
+/* This program is free software; you can redistribute it and/or
+   modify it under the terms of the GNU General Public License as
+   published by the Free Software Foundation; either version 2 of the
+   License, or (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful, but
+   WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   General Public License for more details.  You should have received
+   a copy of the GNU General Public License along with this program;
+   if not, write to the Free Foundation, Inc., 59 Temple Place, Suite
+   330, Boston, MA 02111-1307 USA
+
+   Original implementation was copyright (C) 1997 Makoto Matsumoto and
+   Takuji Nishimura. Coded by Takuji Nishimura, considering the
+   suggestions by Topher Cooper and Marc Rieffel in July-Aug. 1997, "A
+   C-program for MT19937: Integer version (1998/4/6)"
+
+   This implementation copyright (C) 1998 Brian Gough. I reorganized
+   the code to use the module framework of GSL.  The license on this
+   implementation was changed from LGPL to GPL, following paragraph 3
+   of the LGPL, version 2.
+
+   Update:
+
+   The seeding procedure has been updated to match the 10/99 release
+   of MT19937.
+
+   Update:
+
+   The seeding procedure has been updated again to match the 2002
+   release of MT19937
+
+   The original code included the comment: "When you use this, send an
+   email to: matumoto@math.keio.ac.jp with an appropriate reference to
+   your work".
+
+   Makoto Matsumoto has a web page with more information about the
+   generator, http://www.math.keio.ac.jp/~matumoto/emt.html. 
+
+   The paper below has details of the algorithm.
+
+   From: Makoto Matsumoto and Takuji Nishimura, "Mersenne Twister: A
+   623-dimensionally equidistributerd uniform pseudorandom number
+   generator". ACM Transactions on Modeling and Computer Simulation,
+   Vol. 8, No. 1 (Jan. 1998), Pages 3-30
+
+   You can obtain the paper directly from Makoto Matsumoto's web page.
+
+   The period of this generator is 2^{19937} - 1.
+
+*/
+
+#include <config.h>
+#include <stdlib.h>
+#include <gsl/gsl_rng.h>
+
+static inline unsigned long int mt_get (void *vstate);
+static double mt_get_double (void *vstate);
+static void mt_set (void *state, unsigned long int s);
+
+#define N 624   /* Period parameters */
+#define M 397
+
+/* most significant w-r bits */
+static const unsigned long UPPER_MASK = 0x80000000UL;   
+
+/* least significant r bits */
+static const unsigned long LOWER_MASK = 0x7fffffffUL;   
+
+typedef struct
+  {
+    unsigned long mt[N];
+    int mti;
+  }
+mt_state_t;
+
+static inline unsigned long
+mt_get (void *vstate)
+{
+  mt_state_t *state = (mt_state_t *) vstate;
+
+  unsigned long k ;
+  unsigned long int *const mt = state->mt;
+
+#define MAGIC(y) (((y)&0x1) ? 0x9908b0dfUL : 0)
+
+  if (state->mti >= N)
+    {   /* generate N words at one time */
+      int kk;
+
+      for (kk = 0; kk < N - M; kk++)
+        {
+          unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
+          mt[kk] = mt[kk + M] ^ (y >> 1) ^ MAGIC(y);
+        }
+      for (; kk < N - 1; kk++)
+        {
+          unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
+          mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ MAGIC(y);
+        }
+
+      {
+        unsigned long y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
+        mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ MAGIC(y);
+      }
+
+      state->mti = 0;
+    }
+
+  /* Tempering */
+  
+  k = mt[state->mti];
+  k ^= (k >> 11);
+  k ^= (k << 7) & 0x9d2c5680UL;
+  k ^= (k << 15) & 0xefc60000UL;
+  k ^= (k >> 18);
+
+  state->mti++;
+
+  return k;
+}
+
+static double
+mt_get_double (void * vstate)
+{
+  return mt_get (vstate) / 4294967296.0 ;
+}
+
+static void
+mt_set (void *vstate, unsigned long int s)
+{
+  mt_state_t *state = (mt_state_t *) vstate;
+  int i;
+
+  if (s == 0)
+    s = 4357;   /* the default seed is 4357 */
+
+  state->mt[0]= s & 0xffffffffUL;
+
+  for (i = 1; i < N; i++)
+    {
+      /* See Knuth's "Art of Computer Programming" Vol. 2, 3rd
+         Ed. p.106 for multiplier. */
+
+      state->mt[i] =
+        (1812433253UL * (state->mt[i-1] ^ (state->mt[i-1] >> 30)) + i);
+      
+      state->mt[i] &= 0xffffffffUL;
+    }
+
+  state->mti = i;
+}
+
+static void
+mt_1999_set (void *vstate, unsigned long int s)
+{
+  mt_state_t *state = (mt_state_t *) vstate;
+  int i;
+
+  if (s == 0)
+    s = 4357;   /* the default seed is 4357 */
+
+  /* This is the October 1999 version of the seeding procedure. It
+     was updated by the original developers to avoid the periodicity
+     in the simple congruence originally used.
+
+     Note that an ANSI-C unsigned long integer arithmetic is
+     automatically modulo 2^32 (or a higher power of two), so we can
+     safely ignore overflow. */
+
+#define LCG(x) ((69069 * x) + 1) &0xffffffffUL
+
+  for (i = 0; i < N; i++)
+    {
+      state->mt[i] = s & 0xffff0000UL;
+      s = LCG(s);
+      state->mt[i] |= (s &0xffff0000UL) >> 16;
+      s = LCG(s);
+    }
+
+  state->mti = i;
+}
+
+/* This is the original version of the seeding procedure, no longer
+   used but available for compatibility with the original MT19937. */
+
+static void
+mt_1998_set (void *vstate, unsigned long int s)
+{
+  mt_state_t *state = (mt_state_t *) vstate;
+  int i;
+
+  if (s == 0)
+    s = 4357;   /* the default seed is 4357 */
+
+  state->mt[0] = s & 0xffffffffUL;
+
+#define LCG1998(n) ((69069 * n) & 0xffffffffUL)
+
+  for (i = 1; i < N; i++)
+    state->mt[i] = LCG1998 (state->mt[i - 1]);
+
+  state->mti = i;
+}
+
+static const gsl_rng_type mt_type =
+{"mt19937",                     /* name */
+ 0xffffffffUL,                  /* RAND_MAX  */
+ 0,                             /* RAND_MIN  */
+ sizeof (mt_state_t),
+ &mt_set,
+ &mt_get,
+ &mt_get_double};
+
+static const gsl_rng_type mt_1999_type =
+{"mt19937_1999",                /* name */
+ 0xffffffffUL,                  /* RAND_MAX  */
+ 0,                             /* RAND_MIN  */
+ sizeof (mt_state_t),
+ &mt_1999_set,
+ &mt_get,
+ &mt_get_double};
+
+static const gsl_rng_type mt_1998_type =
+{"mt19937_1998",                /* name */
+ 0xffffffffUL,                  /* RAND_MAX  */
+ 0,                             /* RAND_MIN  */
+ sizeof (mt_state_t),
+ &mt_1998_set,
+ &mt_get,
+ &mt_get_double};
+
+const gsl_rng_type *gsl_rng_mt19937 = &mt_type;
+const gsl_rng_type *gsl_rng_mt19937_1999 = &mt_1999_type;
+const gsl_rng_type *gsl_rng_mt19937_1998 = &mt_1998_type;
+
+/* MT19937 is the default generator, so define that here too */
+
+const gsl_rng_type *gsl_rng_default = &mt_type;
+unsigned long int gsl_rng_default_seed = 0;