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

   1 /* rng/tt.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <stdlib.h>
  22 #include <gsl/gsl_rng.h>
  23
  24 /* This is the TT800 twisted GSFR generator for 32 bit integers. It
  25    has been superceded by MT19937 (mt.c). The period is 2^800.
  26
  27    This implementation is based on tt800.c, July 8th 1996 version by
  28    M. Matsumoto, email: matumoto@math.keio.ac.jp
  29
  30    From: Makoto Matsumoto and Yoshiharu Kurita, "Twisted GFSR
  31    Generators II", ACM Transactions on Modelling and Computer
  32    Simulation, Vol. 4, No. 3, 1994, pages 254-266. */
  33
  34 static inline unsigned long int tt_get (void *vstate);
  35 static double tt_get_double (void *vstate);
  36 static void tt_set (void *state, unsigned long int s);
  37
  38 #define N 25
  39 #define M 7
  40
  41 typedef struct
  42   {
  43     int n;
  44     unsigned long int x[N];
  45   }
  46 tt_state_t;
  47
  48 static inline unsigned long int
  49 tt_get (void *vstate)
  50 {
  51   tt_state_t *state = (tt_state_t *) vstate;
  52
  53   /* this is the magic vector, a */
  54
  55   const unsigned long mag01[2] =
  56   {0x00000000, 0x8ebfd028UL};
  57   unsigned long int y;
  58   unsigned long int *const x = state->x;
  59   int n = state->n;
  60
  61   if (n >= N)
  62     {
  63       int i;
  64       for (i = 0; i < N - M; i++)
  65         {
  66           x[i] = x[i + M] ^ (x[i] >> 1) ^ mag01[x[i] % 2];
  67         }
  68       for (; i < N; i++)
  69         {
  70           x[i] = x[i + (M - N)] ^ (x[i] >> 1) ^ mag01[x[i] % 2];
  71         };
  72       n = 0;
  73     }
  74
  75   y = x[n];
  76   y ^= (y << 7) & 0x2b5b2500UL;         /* s and b, magic vectors */
  77   y ^= (y << 15) & 0xdb8b0000UL;        /* t and c, magic vectors */
  78   y &= 0xffffffffUL;    /* you may delete this line if word size = 32 */
  79
  80   /* The following line was added by Makoto Matsumoto in the 1996
  81      version to improve lower bit's correlation.  Delete this line
  82      to use the code published in 1994.  */
  83
  84   y ^= (y >> 16);       /* added to the 1994 version */
  85
  86   state->n = n + 1;
  87
  88   return y;
  89 }
  90
  91 static double
  92 tt_get_double (void * vstate)
  93 {
  94   return tt_get (vstate) / 4294967296.0 ;
  95 }
  96
  97 static void
  98 tt_set (void *vstate, unsigned long int s)
  99 {
 100   tt_state_t *state = (tt_state_t *) vstate;
 101
 102   const tt_state_t init_state =
 103   {0,
 104    {0x95f24dabUL, 0x0b685215UL, 0xe76ccae7UL,
 105     0xaf3ec239UL, 0x715fad23UL, 0x24a590adUL,
 106     0x69e4b5efUL, 0xbf456141UL, 0x96bc1b7bUL,
 107     0xa7bdf825UL, 0xc1de75b7UL, 0x8858a9c9UL,
 108     0x2da87693UL, 0xb657f9ddUL, 0xffdc8a9fUL,
 109     0x8121da71UL, 0x8b823ecbUL, 0x885d05f5UL,
 110     0x4e20cd47UL, 0x5a9ad5d9UL, 0x512c0c03UL,
 111     0xea857ccdUL, 0x4cc1d30fUL, 0x8891a8a1UL,
 112     0xa6b7aadbUL}};
 113
 114
 115   if (s == 0)   /* default seed is given explicitly in the original code */
 116     {
 117       *state = init_state;
 118     }
 119   else
 120     {
 121       int i;
 122
 123       state->n = 0;
 124
 125       state->x[0] = s & 0xffffffffUL;
 126
 127       for (i = 1; i < N; i++)
 128         state->x[i] = (69069 * state->x[i - 1]) & 0xffffffffUL;
 129     }
 130
 131   return;
 132 }
 133
 134 static const gsl_rng_type tt_type =
 135 {"tt800",                       /* name */
 136  0xffffffffUL,                  /* RAND_MAX */
 137  0,                             /* RAND_MIN */
 138  sizeof (tt_state_t),
 139  &tt_set,
 140  &tt_get,
 141  &tt_get_double};
 142
 143 const gsl_rng_type *gsl_rng_tt800 = &tt_type;