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

   1 /* rng/r250.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 shift-register random number generator. The sequence is
  25
  26    x_n = x_{n-103} ^ x_{n-250}        ("^" means XOR)
  27
  28    defined on 32-bit words.
  29
  30    BJG: Note that this implementation actually uses the sequence, x_n
  31    = x_{n-147} ^ x_{n-250} which generates the outputs in
  32    time-reversed order but is otherwise completely equivalent.
  33
  34    The first 250 elements x_1 .. x_250 are first initialized as x_n =
  35    s_n, where s_n = (69069*s_{n-1}) mod 2^32 and s_0=s is the
  36    user-supplied seed. To ensure that the sequence does not lie on a
  37    subspace we force 32 of the entries to be linearly independent.  We
  38    take the 32 elements x[3], x[10], x[17], x[24], ..., 213 and apply
  39    the following operations,
  40
  41    x[3]   &= 11111111111111111111111111111111
  42    x[3]   |= 10000000000000000000000000000000
  43    x[10]  &= 01111111111111111111111111111111
  44    x[10]  |= 01000000000000000000000000000000
  45    x[17]  &= 00111111111111111111111111111111
  46    x[17]  |= 00100000000000000000000000000000
  47    ....      ...
  48    x[206] &= 00000000000000000000000000000111
  49    x[206] |= 00000000000000000000000000000100
  50    x[213] &= 00000000000000000000000000000011
  51    x[213] |= 00000000000000000000000000000010
  52    x[220] &= 00000000000000000000000000000001
  53    x[220] |= 00000000000000000000000000000001
  54
  55    i.e. if we consider the bits of the 32 elements as forming a 32x32
  56    array then we are setting the diagonal bits of the array to one and
  57    masking the lower triangle below the diagonal to zero.
  58
  59    With this initialization procedure the theoretical value of
  60    x_{10001} is 1100653588 for s = 1 (Actually I got this by running
  61    the original code). The subscript 10001 means (1) seed the
  62    generator with s = 1 and then do 10000 actual iterations.
  63
  64    The period of this generator is about 2^250.
  65
  66    The algorithm works for any number of bits. It is implemented here
  67    for 32 bits.
  68
  69    From: S. Kirkpatrick and E. Stoll, "A very fast shift-register
  70    sequence random number generator", Journal of Computational Physics,
  71    40, 517-526 (1981). */
  72
  73 static inline unsigned long int r250_get (void *vstate);
  74 static double r250_get_double (void *vstate);
  75 static void r250_set (void *state, unsigned long int s);
  76
  77 typedef struct
  78   {
  79     int i;
  80     unsigned long x[250];
  81   }
  82 r250_state_t;
  83
  84 static inline unsigned long int
  85 r250_get (void *vstate)
  86 {
  87   r250_state_t *state = (r250_state_t *) vstate;
  88   unsigned long int k;
  89   int j;
  90
  91   int i = state->i;
  92
  93   if (i >= 147)
  94     {
  95       j = i - 147;
  96     }
  97   else
  98     {
  99       j = i + 103;
 100     }
 101
 102   k = state->x[i] ^ state->x[j];
 103   state->x[i] = k;
 104
 105   if (i >= 249)
 106     {
 107       state->i = 0;
 108     }
 109   else
 110     {
 111       state->i = i + 1;
 112     }
 113
 114   return k;
 115 }
 116
 117 static double
 118 r250_get_double (void *vstate)
 119 {
 120   return r250_get (vstate) /  4294967296.0 ;
 121 }
 122
 123 static void
 124 r250_set (void *vstate, unsigned long int s)
 125 {
 126   r250_state_t *state = (r250_state_t *) vstate;
 127
 128   int i;
 129
 130   if (s == 0)
 131     s = 1;      /* default seed is 1 */
 132
 133   state->i = 0;
 134
 135 #define LCG(n) ((69069 * n) & 0xffffffffUL)
 136
 137   for (i = 0; i < 250; i++)     /* Fill the buffer  */
 138     {
 139       s = LCG (s);
 140       state->x[i] = s;
 141     }
 142
 143   {
 144     /* Masks for turning on the diagonal bit and turning off the
 145        leftmost bits */
 146
 147     unsigned long int msb = 0x80000000UL;
 148     unsigned long int mask = 0xffffffffUL;
 149
 150     for (i = 0; i < 32; i++)
 151       {
 152         int k = 7 * i + 3;      /* Select a word to operate on        */
 153         state->x[k] &= mask;    /* Turn off bits left of the diagonal */
 154         state->x[k] |= msb;     /* Turn on the diagonal bit           */
 155         mask >>= 1;
 156         msb >>= 1;
 157       }
 158   }
 159
 160   return;
 161 }
 162
 163 static const gsl_rng_type r250_type =
 164 {"r250",                        /* name */
 165  0xffffffffUL,                  /* RAND_MAX */
 166  0,                             /* RAND_MIN */
 167  sizeof (r250_state_t),
 168  &r250_set,
 169  &r250_get,
 170  &r250_get_double};
 171
 172 const gsl_rng_type *gsl_rng_r250 = &r250_type;