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

   1 /* This program is free software; you can redistribute it and/or
   2    modify it under the terms of the GNU General Public License as
   3    published by the Free Software Foundation; either version 3 of the
   4    License, or (at your option) any later version.
   5
   6    This program is distributed in the hope that it will be useful, but
   7    WITHOUT ANY WARRANTY; without even the implied warranty of
   8    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   9    General Public License for more details.  You should have received
  10    a copy of the GNU General Public License along with this program;
  11    if not, write to the Free Foundation, Inc., 59 Temple Place, Suite
  12    330, Boston, MA 02111-1307 USA
  13
  14    From Robert M. Ziff, "Four-tap shift-register-sequence
  15    random-number generators," Computers in Physics 12(4), Jul/Aug
  16    1998, pp 385-392.  A generalized feedback shift-register (GFSR)
  17    is basically an xor-sum of particular past lagged values.  A
  18    four-tap register looks like:
  19       ra[nd] = ra[nd-A] ^ ra[nd-B] ^ ra[nd-C] ^ ra[nd-D]
  20
  21    Ziff notes that "it is now widely known" that two-tap registers
  22    have serious flaws, the most obvious one being the three-point
  23    correlation that comes from the defn of the generator.  Nice
  24    mathematical properties can be derived for GFSR's, and numerics
  25    bears out the claim that 4-tap GFSR's with appropriately chosen
  26    offsets are as random as can be measured, using the author's test.
  27
  28    This implementation uses the values suggested the the author's
  29    example on p392, but altered to fit the GSL framework.  The "state"
  30    is 2^14 longs, or 64Kbytes; 2^14 is the smallest power of two that
  31    is larger than D, the largest offset.  We really only need a state
  32    with the last D values, but by going to a power of two, we can do a
  33    masking operation instead of a modulo, and this is presumably
  34    faster, though I haven't actually tried it.  The article actually
  35    suggested a short/fast hack:
  36
  37    #define RandomInteger (++nd, ra[nd&M]=ra[(nd-A)&M]\
  38                           ^ra[(nd-B)&M]^ra[(nd-C)&M]^ra[(nd-D)&M])
  39
  40    so that (as long as you've defined nd,ra[M+1]), then you ca do things
  41    like: 'if (RandomInteger < p) {...}'.
  42
  43    Note that n&M varies from 0 to M, *including* M, so that the
  44    array has to be of size M+1.  Since M+1 is a power of two, n&M
  45    is a potentially quicker implementation of the equivalent n%(M+1).
  46
  47    This implementation copyright (C) 1998 James Theiler, based on
  48    the example mt.c in the GSL, as implemented by Brian Gough.
  49 */
  50
  51 #include <config.h>
  52 #include <stdlib.h>
  53 #include <gsl/gsl_rng.h>
  54
  55 static inline unsigned long int gfsr4_get (void *vstate);
  56 static double gfsr4_get_double (void *vstate);
  57 static void gfsr4_set (void *state, unsigned long int s);
  58
  59 /* Magic numbers */
  60 #define A 471
  61 #define B 1586
  62 #define C 6988
  63 #define D 9689
  64 #define M 16383 /* = 2^14-1 */
  65 /* #define M 0x0003fff */
  66
  67 typedef struct
  68   {
  69     int nd;
  70     unsigned long ra[M+1];
  71   }
  72 gfsr4_state_t;
  73
  74 static inline unsigned long
  75 gfsr4_get (void *vstate)
  76 {
  77   gfsr4_state_t *state = (gfsr4_state_t *) vstate;
  78
  79   state->nd = ((state->nd)+1)&M;
  80   return state->ra[(state->nd)] =
  81       state->ra[((state->nd)+(M+1-A))&M]^
  82       state->ra[((state->nd)+(M+1-B))&M]^
  83       state->ra[((state->nd)+(M+1-C))&M]^
  84       state->ra[((state->nd)+(M+1-D))&M];
  85
  86 }
  87
  88 static double
  89 gfsr4_get_double (void * vstate)
  90 {
  91   return gfsr4_get (vstate) / 4294967296.0 ;
  92 }
  93
  94 static void
  95 gfsr4_set (void *vstate, unsigned long int s)
  96 {
  97   gfsr4_state_t *state = (gfsr4_state_t *) vstate;
  98   int i, j;
  99   /* Masks for turning on the diagonal bit and turning off the
 100      leftmost bits */
 101   unsigned long int msb = 0x80000000UL;
 102   unsigned long int mask = 0xffffffffUL;
 103
 104   if (s == 0)
 105     s = 4357;   /* the default seed is 4357 */
 106
 107   /* We use the congruence s_{n+1} = (69069*s_n) mod 2^32 to
 108      initialize the state. This works because ANSI-C unsigned long
 109      integer arithmetic is automatically modulo 2^32 (or a higher
 110      power of two), so we can safely ignore overflow. */
 111
 112 #define LCG(n) ((69069 * n) & 0xffffffffUL)
 113
 114   /* Brian Gough suggests this to avoid low-order bit correlations */
 115   for (i = 0; i <= M; i++)
 116     {
 117       unsigned long t = 0 ;
 118       unsigned long bit = msb ;
 119       for (j = 0; j < 32; j++)
 120         {
 121           s = LCG(s) ;
 122           if (s & msb)
 123             t |= bit ;
 124           bit >>= 1 ;
 125         }
 126       state->ra[i] = t ;
 127     }
 128
 129   /* Perform the "orthogonalization" of the matrix */
 130   /* Based on the orthogonalization used in r250, as suggested initially
 131    * by Kirkpatrick and Stoll, and pointed out to me by Brian Gough
 132    */
 133
 134   /* BJG: note that this orthogonalisation doesn't have any effect
 135      here because the the initial 6695 elements do not participate in
 136      the calculation.  For practical purposes this orthogonalisation
 137      is somewhat irrelevant, because the probability of the original
 138      sequence being degenerate should be exponentially small. */
 139
 140   for (i=0; i<32; ++i) {
 141       int k=7+i*3;
 142       state->ra[k] &= mask;     /* Turn off bits left of the diagonal */
 143       state->ra[k] |= msb;      /* Turn on the diagonal bit           */
 144       mask >>= 1;
 145       msb >>= 1;
 146   }
 147
 148   state->nd = i;
 149 }
 150
 151 static const gsl_rng_type gfsr4_type =
 152 {"gfsr4",                       /* name */
 153  0xffffffffUL,                  /* RAND_MAX  */
 154  0,                             /* RAND_MIN  */
 155  sizeof (gfsr4_state_t),
 156  &gfsr4_set,
 157  &gfsr4_get,
 158  &gfsr4_get_double};
 159
 160 const gsl_rng_type *gsl_rng_gfsr4 = &gfsr4_type;
 161
 162
 163
 164
 165