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

   1 /* rng/knuthran2002.c
   2  *
   3  * Copyright (C) 2007 Brian Gough
   4  * Copyright (C) 2001, 2007 Brian Gough, Carlo Perassi
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 3 of the License, or (at
   9  * your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  19  */
  20
  21 /*
  22  * This generator is taken from
  23  *
  24  * Donald E. Knuth, The Art of Computer Programming, Volume 2, Section 3.6
  25  * Third Edition, Addison-Wesley,
  26  *
  27  * The modifications introduced in the 9th printing (2002) are
  28  * included here; there's no backwards compatibility with the
  29  * original.  [ see http://www-cs-faculty.stanford.edu/~knuth/taocp.html ]
  30  *
  31  */
  32
  33 #include <config.h>
  34 #include <stdlib.h>
  35 #include <gsl/gsl_rng.h>
  36
  37 #define BUFLEN 1009             /* length of the buffer aa[] */
  38 #define KK 100                  /* the long lag */
  39 #define LL 37                   /* the short lag */
  40 #define MM (1L << 30)           /* the modulus */
  41 #define TT 70                   /* guaranteed separation between streams */
  42
  43 #define is_odd(x) ((x) & 1)                       /* the units bit of x */
  44 #define mod_diff(x, y) (((x) - (y)) & (MM - 1))   /* (x - y) mod MM */
  45
  46 static inline void ran_array (long int aa[], unsigned int n,
  47                               long int ran_x[]);
  48 static inline unsigned long int ran_get (void *vstate);
  49 static double ran_get_double (void *vstate);
  50 static void ran_set (void *state, unsigned long int s);
  51
  52 typedef struct
  53 {
  54   unsigned int i;
  55   long int aa[BUFLEN];
  56   long int ran_x[KK];  /* the generator state */
  57 }
  58 ran_state_t;
  59
  60 static inline void
  61 ran_array (long int aa[], unsigned int n, long int ran_x[])
  62 {
  63   unsigned int i;
  64   unsigned int j;
  65
  66   for (j = 0; j < KK; j++)
  67     aa[j] = ran_x[j];
  68
  69   for (; j < n; j++)
  70     aa[j] = mod_diff (aa[j - KK], aa[j - LL]);
  71
  72   for (i = 0; i < LL; i++, j++)
  73     ran_x[i] = mod_diff (aa[j - KK], aa[j - LL]);
  74
  75   for (; i < KK; i++, j++)
  76     ran_x[i] = mod_diff (aa[j - KK], ran_x[i - LL]);
  77 }
  78
  79 static inline unsigned long int
  80 ran_get (void *vstate)
  81 {
  82   ran_state_t *state = (ran_state_t *) vstate;
  83
  84   unsigned int i = state->i;
  85   unsigned long int v;
  86
  87   if (i  == 0)
  88     {
  89       /* fill buffer with new random numbers */
  90       ran_array (state->aa, BUFLEN, state->ran_x);
  91     }
  92
  93   v = state->aa[i];
  94
  95   state->i = (i + 1) % KK;
  96
  97   return v;
  98 }
  99
 100 static double
 101 ran_get_double (void *vstate)
 102 {
 103   ran_state_t *state = (ran_state_t *) vstate;
 104
 105   return ran_get (state) / 1073741824.0;        /* RAND_MAX + 1 */
 106 }
 107
 108 static void
 109 ran_set (void *vstate, unsigned long int s)
 110 {
 111   ran_state_t *state = (ran_state_t *) vstate;
 112
 113   long x[KK + KK - 1];          /* the preparation buffer */
 114
 115   register int j;
 116   register int t;
 117   register long ss;
 118
 119   if (s == 0 )
 120     s = 314159;                 /* default seed used by Knuth */
 121
 122   ss = (s + 2)&(MM-2);
 123
 124   for (j = 0; j < KK; j++)
 125     {
 126       x[j] = ss;                /* bootstrap the buffer */
 127       ss <<= 1;
 128       if (ss >= MM)             /* cyclic shift 29 bits */
 129         ss -= MM - 2;
 130     }
 131   x[1]++;                       /* make x[1] (and only x[1]) odd */
 132
 133   ss = s & (MM - 1);
 134   t = TT - 1;
 135   while (t)
 136     {
 137       for (j = KK - 1; j > 0; j--)      /* square */
 138         {
 139           x[j + j] = x[j];
 140           x[j + j - 1] = 0;
 141         }
 142
 143       for (j = KK + KK - 2; j >= KK; j--)
 144         {
 145           x[j - (KK - LL)] = mod_diff (x[j - (KK - LL)], x[j]);
 146           x[j - KK] = mod_diff (x[j - KK], x[j]);
 147         }
 148
 149       if (is_odd (ss))
 150         {                       /* multiply by "z" */
 151           for (j = KK; j > 0; j--)
 152             {
 153               x[j] = x[j - 1];
 154             }
 155           x[0] = x[KK];         /* shift the buffer cyclically */
 156           x[LL] = mod_diff (x[LL], x[KK]);
 157         }
 158
 159       if (ss)
 160         ss >>= 1;
 161       else
 162         t--;
 163     }
 164
 165   for (j = 0; j < LL; j++)
 166     state->ran_x[j + KK - LL] = x[j];
 167   for (; j < KK; j++)
 168     state->ran_x[j - LL] = x[j];
 169
 170
 171   for (j = 0; j< 10; j++)
 172     ran_array(x, KK+KK-1, state->ran_x);  /* warm things up */
 173
 174   state->i = 0;
 175
 176   return;
 177 }
 178
 179 static const gsl_rng_type ran_type = {
 180   "knuthran2002",                /* name */
 181   0x3fffffffUL,                 /* RAND_MAX = (2 ^ 30) - 1 */
 182   0,                            /* RAND_MIN */
 183   sizeof (ran_state_t),
 184   &ran_set,
 185   &ran_get,
 186   &ran_get_double
 187 };
 188
 189 const gsl_rng_type *gsl_rng_knuthran2002 = &ran_type;