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

   1 /* rng/taus.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 a maximally equidistributed combined Tausworthe
  25    generator. The sequence is,
  26
  27    x_n = (s1_n ^ s2_n ^ s3_n)
  28
  29    s1_{n+1} = (((s1_n & 4294967294) <<12) ^ (((s1_n <<13) ^ s1_n) >>19))
  30    s2_{n+1} = (((s2_n & 4294967288) << 4) ^ (((s2_n << 2) ^ s2_n) >>25))
  31    s3_{n+1} = (((s3_n & 4294967280) <<17) ^ (((s3_n << 3) ^ s3_n) >>11))
  32
  33    computed modulo 2^32. In the three formulas above '^' means
  34    exclusive-or (C-notation), not exponentiation. Note that the
  35    algorithm relies on the properties of 32-bit unsigned integers (it
  36    is formally defined on bit-vectors of length 32). I have added a
  37    bitmask to make it work on 64 bit machines.
  38
  39    We initialize the generator with s1_1 .. s3_1 = s_n MOD m, where
  40    s_n = (69069 * s_{n-1}) mod 2^32, and s_0 = s is the user-supplied
  41    seed.
  42
  43    The theoretical value of x_{10007} is 2733957125. The subscript
  44    10007 means (1) seed the generator with s=1 (2) do six warm-up
  45    iterations, (3) then do 10000 actual iterations.
  46
  47    The period of this generator is about 2^88.
  48
  49    From: P. L'Ecuyer, "Maximally Equidistributed Combined Tausworthe
  50    Generators", Mathematics of Computation, 65, 213 (1996), 203--213.
  51
  52    This is available on the net from L'Ecuyer's home page,
  53
  54    http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme.ps
  55    ftp://ftp.iro.umontreal.ca/pub/simulation/lecuyer/papers/tausme.ps
  56
  57    Update: April 2002
  58
  59    There is an erratum in the paper "Tables of Maximally
  60    Equidistributed Combined LFSR Generators", Mathematics of
  61    Computation, 68, 225 (1999), 261--269:
  62    http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme2.ps
  63
  64         ... the k_j most significant bits of z_j must be non-
  65         zero, for each j. (Note: this restriction also applies to the
  66         computer code given in [4], but was mistakenly not mentioned in
  67         that paper.)
  68
  69    This affects the seeding procedure by imposing the requirement
  70    s1 > 1, s2 > 7, s3 > 15.
  71
  72    The generator taus2 has been added to satisfy this requirement.
  73    The original taus generator is unchanged.
  74
  75    Update: November 2002
  76
  77    There was a bug in the correction to the seeding procedure for s2.
  78    It affected the following seeds 254679140 1264751179 1519430319
  79    2274823218 2529502358 3284895257 3539574397 (s2 < 8).
  80
  81 */
  82
  83 static inline unsigned long int taus_get (void *vstate);
  84 static double taus_get_double (void *vstate);
  85 static void taus_set (void *state, unsigned long int s);
  86
  87 typedef struct
  88   {
  89     unsigned long int s1, s2, s3;
  90   }
  91 taus_state_t;
  92
  93 static inline unsigned long
  94 taus_get (void *vstate)
  95 {
  96   taus_state_t *state = (taus_state_t *) vstate;
  97
  98 #define MASK 0xffffffffUL
  99 #define TAUSWORTHE(s,a,b,c,d) (((s &c) <<d) &MASK) ^ ((((s <<a) &MASK)^s) >>b)
 100
 101   state->s1 = TAUSWORTHE (state->s1, 13, 19, 4294967294UL, 12);
 102   state->s2 = TAUSWORTHE (state->s2, 2, 25, 4294967288UL, 4);
 103   state->s3 = TAUSWORTHE (state->s3, 3, 11, 4294967280UL, 17);
 104
 105   return (state->s1 ^ state->s2 ^ state->s3);
 106 }
 107
 108 static double
 109 taus_get_double (void *vstate)
 110 {
 111   return taus_get (vstate) / 4294967296.0 ;
 112 }
 113
 114 static void
 115 taus_set (void *vstate, unsigned long int s)
 116 {
 117   taus_state_t *state = (taus_state_t *) vstate;
 118
 119   if (s == 0)
 120     s = 1;      /* default seed is 1 */
 121
 122 #define LCG(n) ((69069 * n) & 0xffffffffUL)
 123   state->s1 = LCG (s);
 124   state->s2 = LCG (state->s1);
 125   state->s3 = LCG (state->s2);
 126
 127   /* "warm it up" */
 128   taus_get (state);
 129   taus_get (state);
 130   taus_get (state);
 131   taus_get (state);
 132   taus_get (state);
 133   taus_get (state);
 134   return;
 135 }
 136
 137 static void
 138 taus2_set (void *vstate, unsigned long int s)
 139 {
 140   taus_state_t *state = (taus_state_t *) vstate;
 141
 142   if (s == 0)
 143     s = 1;      /* default seed is 1 */
 144
 145 #define LCG(n) ((69069 * n) & 0xffffffffUL)
 146   state->s1 = LCG (s);
 147   if (state->s1 < 2) state->s1 += 2UL;
 148   state->s2 = LCG (state->s1);
 149   if (state->s2 < 8) state->s2 += 8UL;
 150   state->s3 = LCG (state->s2);
 151   if (state->s3 < 16) state->s3 += 16UL;
 152
 153   /* "warm it up" */
 154   taus_get (state);
 155   taus_get (state);
 156   taus_get (state);
 157   taus_get (state);
 158   taus_get (state);
 159   taus_get (state);
 160   return;
 161 }
 162
 163
 164 static const gsl_rng_type taus_type =
 165 {"taus",                        /* name */
 166  0xffffffffUL,                  /* RAND_MAX */
 167  0,                             /* RAND_MIN */
 168  sizeof (taus_state_t),
 169  &taus_set,
 170  &taus_get,
 171  &taus_get_double};
 172
 173 const gsl_rng_type *gsl_rng_taus = &taus_type;
 174
 175 static const gsl_rng_type taus2_type =
 176 {"taus2",                       /* name */
 177  0xffffffffUL,                  /* RAND_MAX */
 178  0,                             /* RAND_MIN */
 179  sizeof (taus_state_t),
 180  &taus2_set,
 181  &taus_get,
 182  &taus_get_double};
 183
 184 const gsl_rng_type *gsl_rng_taus2 = &taus2_type;