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

   1 /* rng/knuthran.c
   2  *
   3  * Copyright (C) 2001, 2007 Brian Gough, Carlo Perassi
   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 /*
  21  * This generator is taken from
  22  *
  23  * Donald E. Knuth
  24  * The Art of Computer Programming
  25  * Volume 2
  26  * Third Edition
  27  * Addison-Wesley
  28  * Section 3.6
  29  *
  30  * The comments are taken from the book
  31  * Our comments are signed
  32  */
  33
  34 #include <config.h>
  35 #include <stdlib.h>
  36 #include <gsl/gsl_rng.h>
  37
  38 #define BUFLEN 2009             /* [Brian]: length of the buffer aa[] */
  39 #define KK 100                  /* the long lag */
  40 #define LL 37                   /* the short lag */
  41 #define MM (1L << 30)           /* the modulus */
  42 #define TT 70                   /* guaranteed separation between streams */
  43
  44 #define evenize(x) ((x) & (MM - 2))     /* make x even */
  45 #define is_odd(x) ((x) & 1)     /* the units bit of x */
  46 #define mod_diff(x, y) (((x) - (y)) & (MM - 1)) /* (x - y) mod MM */
  47
  48 static inline void ran_array (unsigned long int aa[], unsigned int n,
  49                               unsigned long int ran_x[]);
  50 static inline unsigned long int ran_get (void *vstate);
  51 static double ran_get_double (void *vstate);
  52 static void ran_set (void *state, unsigned long int s);
  53
  54 typedef struct
  55 {
  56   unsigned int i;
  57   unsigned long int aa[BUFLEN]; /* [Carlo]: I can't pass n to ran_array like
  58                                    Knuth does */
  59   unsigned long int ran_x[KK];  /* the generator state */
  60 }
  61 ran_state_t;
  62
  63 static inline void
  64 ran_array (unsigned long int aa[], unsigned int n, unsigned long int ran_x[])
  65 {
  66   unsigned int i;
  67   unsigned int j;
  68
  69   for (j = 0; j < KK; j++)
  70     aa[j] = ran_x[j];
  71
  72   for (; j < n; j++)
  73     aa[j] = mod_diff (aa[j - KK], aa[j - LL]);
  74
  75   for (i = 0; i < LL; i++, j++)
  76     ran_x[i] = mod_diff (aa[j - KK], aa[j - LL]);
  77
  78   for (; i < KK; i++, j++)
  79     ran_x[i] = mod_diff (aa[j - KK], ran_x[i - LL]);
  80 }
  81
  82 static inline unsigned long int
  83 ran_get (void *vstate)
  84 {
  85   ran_state_t *state = (ran_state_t *) vstate;
  86
  87   unsigned int i = state->i;
  88
  89   if (i == 0)
  90     {
  91       /* fill buffer with new random numbers */
  92       ran_array (state->aa, BUFLEN, state->ran_x);
  93     }
  94
  95   state->i = (i + 1) % BUFLEN;
  96
  97   return state->aa[i];
  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;        /* [Carlo]: 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 = evenize (s + 2);
 118
 119   for (j = 0; j < KK; j++)
 120     {
 121       x[j] = ss;                /* bootstrap the buffer */
 122       ss <<= 1;
 123       if (ss >= MM)             /* cyclic shift 29 bits */
 124         ss -= MM - 2;
 125     }
 126   for (; j < KK + KK - 1; j++)
 127     x[j] = 0;
 128   x[1]++;                       /* make x[1] (and only x[1]) odd */
 129   ss = s & (MM - 1);
 130   t = TT - 1;
 131   while (t)
 132     {
 133       for (j = KK - 1; j > 0; j--)      /* square */
 134         x[j + j] = x[j];
 135       for (j = KK + KK - 2; j > KK - LL; j -= 2)
 136         x[KK + KK - 1 - j] = evenize (x[j]);
 137       for (j = KK + KK - 2; j >= KK; j--)
 138         if (is_odd (x[j]))
 139           {
 140             x[j - (KK - LL)] = mod_diff (x[j - (KK - LL)], x[j]);
 141             x[j - KK] = mod_diff (x[j - KK], x[j]);
 142           }
 143       if (is_odd (ss))
 144         {                       /* multiply by "z" */
 145           for (j = KK; j > 0; j--)
 146             x[j] = x[j - 1];
 147           x[0] = x[KK];         /* shift the buffer cyclically */
 148           if (is_odd (x[KK]))
 149             x[LL] = mod_diff (x[LL], x[KK]);
 150         }
 151       if (ss)
 152         ss >>= 1;
 153       else
 154         t--;
 155     }
 156
 157   state->i = 0;
 158
 159   for (j = 0; j < LL; j++)
 160     state->ran_x[j + KK - LL] = x[j];
 161   for (; j < KK; j++)
 162     state->ran_x[j - LL] = x[j];
 163
 164   return;
 165 }
 166
 167 static const gsl_rng_type ran_type = {
 168   "knuthran",                   /* name */
 169   0x3fffffffUL,                 /* RAND_MAX *//* [Carlo]: (2 ^ 30) - 1 */
 170   0,                            /* RAND_MIN */
 171   sizeof (ran_state_t),
 172   &ran_set,
 173   &ran_get,
 174   &ran_get_double
 175 };
 176
 177 const gsl_rng_type *gsl_rng_knuthran = &ran_type;