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

   1 /* rng/cmrg.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 combined multiple recursive generator. The sequence is,
  25
  26    z_n = (x_n - y_n) mod m1
  27
  28    where the two underlying generators x and y are,
  29
  30    x_n = (a_{1} x_{n-1} + a_{2} x_{n-2} + a_{3} x_{n-3}) mod m1
  31    y_n = (b_{1} y_{n-1} + b_{2} y_{n-2} + b_{3} y_{n-3}) mod m2
  32
  33    with coefficients a11 ... a23,
  34
  35    a_{1} = 0,     a_{2} = 63308, a_{3} = -183326
  36    b_{1} = 86098, b_{2} = 0,     b_{3} = -539608
  37
  38    and moduli m1, m2,
  39
  40    m1 = 2^31 - 1 = 2147483647
  41    m2 = 2^31 - 2000169 = 2145483479
  42
  43    We initialize the generator with
  44
  45    x_1 = s_1 MOD m1, x_2 = s_2 MOD m1, x_3 = s_3 MOD m1
  46    y_1 = s_4 MOD m2, y_2 = s_5 MOD m2, y_3 = s_6 MOD m2
  47
  48    where s_n = (69069 * s_{n-1}) mod 2^32 and s_0 = s is the
  49    user-supplied seed.
  50
  51    NOTE: According to the paper the initial values for x_n must lie in
  52    the range 0 <= x_n <= (m1 - 1) and the initial values for y_n must
  53    lie in the range 0 <= y_n <= (m2 - 1), with at least one non-zero
  54    value -- our seeding procedure satisfies these constraints.
  55
  56    We then use 7 iterations of the generator to "warm up" the internal
  57    state.
  58
  59    The theoretical value of z_{10008} is 719452880. The subscript 10008
  60    means (1) seed the generator with s=1, (2) do the seven warm-up
  61    iterations that are part of the seeding process, (3) then do 10000
  62    actual iterations.
  63
  64    The period of this generator is about 2^205.
  65
  66    From: P. L'Ecuyer, "Combined Multiple Recursive Random Number
  67    Generators," Operations Research, 44, 5 (1996), 816--822.
  68
  69    This is available on the net from L'Ecuyer's home page,
  70
  71    http://www.iro.umontreal.ca/~lecuyer/myftp/papers/combmrg.ps
  72    ftp://ftp.iro.umontreal.ca/pub/simulation/lecuyer/papers/combmrg.ps */
  73
  74 static inline unsigned long int cmrg_get (void *vstate);
  75 static double cmrg_get_double (void *vstate);
  76 static void cmrg_set (void *state, unsigned long int s);
  77
  78 static const long int m1 = 2147483647, m2 = 2145483479;
  79
  80 static const long int a2 = 63308, qa2 = 33921, ra2 = 12979;
  81 static const long int a3 = -183326, qa3 = 11714, ra3 = 2883;
  82 static const long int b1 = 86098, qb1 = 24919, rb1 = 7417;
  83 static const long int b3 = -539608, qb3 = 3976, rb3 = 2071;
  84
  85 typedef struct
  86   {
  87     long int x1, x2, x3;        /* first component */
  88     long int y1, y2, y3;        /* second component */
  89   }
  90 cmrg_state_t;
  91
  92 static inline unsigned long int
  93 cmrg_get (void *vstate)
  94 {
  95   cmrg_state_t *state = (cmrg_state_t *) vstate;
  96
  97   /* Component 1 */
  98
  99   {
 100     long int h3 = state->x3 / qa3;
 101     long int p3 = -a3 * (state->x3 - h3 * qa3) - h3 * ra3;
 102
 103     long int h2 = state->x2 / qa2;
 104     long int p2 = a2 * (state->x2 - h2 * qa2) - h2 * ra2;
 105
 106     if (p3 < 0)
 107       p3 += m1;
 108     if (p2 < 0)
 109       p2 += m1;
 110
 111     state->x3 = state->x2;
 112     state->x2 = state->x1;
 113     state->x1 = p2 - p3;
 114     if (state->x1 < 0)
 115       state->x1 += m1;
 116   }
 117
 118   /* Component 2 */
 119
 120   {
 121     long int h3 = state->y3 / qb3;
 122     long int p3 = -b3 * (state->y3 - h3 * qb3) - h3 * rb3;
 123
 124     long int h1 = state->y1 / qb1;
 125     long int p1 = b1 * (state->y1 - h1 * qb1) - h1 * rb1;
 126
 127     if (p3 < 0)
 128       p3 += m2;
 129     if (p1 < 0)
 130       p1 += m2;
 131
 132     state->y3 = state->y2;
 133     state->y2 = state->y1;
 134     state->y1 = p1 - p3;
 135     if (state->y1 < 0)
 136       state->y1 += m2;
 137   }
 138
 139   if (state->x1 < state->y1)
 140     return (state->x1 - state->y1 + m1);
 141   else
 142     return (state->x1 - state->y1);
 143 }
 144
 145 static double
 146 cmrg_get_double (void *vstate)
 147 {
 148   return cmrg_get (vstate) / 2147483647.0 ;
 149 }
 150
 151
 152 static void
 153 cmrg_set (void *vstate, unsigned long int s)
 154 {
 155   /* An entirely adhoc way of seeding! This does **not** come from
 156      L'Ecuyer et al */
 157
 158   cmrg_state_t *state = (cmrg_state_t *) vstate;
 159
 160   if (s == 0)
 161     s = 1;      /* default seed is 1 */
 162
 163 #define LCG(n) ((69069 * n) & 0xffffffffUL)
 164   s = LCG (s);
 165   state->x1 = s % m1;
 166   s = LCG (s);
 167   state->x2 = s % m1;
 168   s = LCG (s);
 169   state->x3 = s % m1;
 170
 171   s = LCG (s);
 172   state->y1 = s % m2;
 173   s = LCG (s);
 174   state->y2 = s % m2;
 175   s = LCG (s);
 176   state->y3 = s % m2;
 177
 178   /* "warm it up" */
 179   cmrg_get (state);
 180   cmrg_get (state);
 181   cmrg_get (state);
 182   cmrg_get (state);
 183   cmrg_get (state);
 184   cmrg_get (state);
 185   cmrg_get (state);
 186 }
 187
 188 static const gsl_rng_type cmrg_type =
 189 {"cmrg",                        /* name */
 190  2147483646,                    /* RAND_MAX */
 191  0,                             /* RAND_MIN */
 192  sizeof (cmrg_state_t),
 193  &cmrg_set,
 194  &cmrg_get,
 195  &cmrg_get_double};
 196
 197 const gsl_rng_type *gsl_rng_cmrg = &cmrg_type;