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

   1 /* rng/ranlxs.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 an implementation of M. Luescher's second generation
  25    version of the RANLUX generator.
  26
  27    Thanks to Martin Luescher for providing information on this
  28    generator.
  29
  30  */
  31
  32 static unsigned long int ranlxs_get (void *vstate);
  33 static inline double ranlxs_get_double (void *vstate);
  34 static void ranlxs_set_lux (void *state, unsigned long int s, unsigned int luxury);
  35 static void ranlxs0_set (void *state, unsigned long int s);
  36 static void ranlxs1_set (void *state, unsigned long int s);
  37 static void ranlxs2_set (void *state, unsigned long int s);
  38
  39 static const int next[12] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0};
  40 static const int snext[24] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
  41                               14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0};
  42
  43 static const double sbase = 16777216.0;         /* 2^24 */
  44 static const double sone_bit = 1.0 / 16777216.0;        /* 1/2^24 */
  45 static const double one_bit = 1.0 / 281474976710656.0;  /* 1/2^48 */
  46
  47 static const double shift = 268435456.0;        /* 2^28 */
  48
  49 #define RANLUX_STEP(x1,x2,i1,i2,i3)      \
  50           x1=xdbl[i1] - xdbl[i2];        \
  51           if (x2 < 0)                    \
  52           {                              \
  53             x1-=one_bit;                 \
  54             x2+=1;                       \
  55           }                              \
  56           xdbl[i3]=x2
  57
  58 typedef struct
  59   {
  60     double xdbl[12], ydbl[12];  /* doubles first so they are 8-byte aligned */
  61     double carry;
  62     float xflt[24];
  63     unsigned int ir;
  64     unsigned int jr;
  65     unsigned int is;
  66     unsigned int is_old;
  67     unsigned int pr;
  68   }
  69 ranlxs_state_t;
  70
  71 static void increment_state (ranlxs_state_t * state);
  72
  73 static void
  74 increment_state (ranlxs_state_t * state)
  75 {
  76   int k, kmax, m;
  77   double x, y1, y2, y3;
  78
  79   float *xflt = state->xflt;
  80   double *xdbl = state->xdbl;
  81   double *ydbl = state->ydbl;
  82   double carry = state->carry;
  83   unsigned int ir = state->ir;
  84   unsigned int jr = state->jr;
  85
  86   for (k = 0; ir > 0; ++k)
  87     {
  88       y1 = xdbl[jr] - xdbl[ir];
  89       y2 = y1 - carry;
  90       if (y2 < 0)
  91         {
  92           carry = one_bit;
  93           y2 += 1;
  94         }
  95       else
  96         {
  97           carry = 0;
  98         }
  99       xdbl[ir] = y2;
 100       ir = next[ir];
 101       jr = next[jr];
 102     }
 103
 104   kmax = state->pr - 12;
 105
 106   for (; k <= kmax; k += 12)
 107     {
 108       y1 = xdbl[7] - xdbl[0];
 109       y1 -= carry;
 110
 111       RANLUX_STEP (y2, y1, 8, 1, 0);
 112       RANLUX_STEP (y3, y2, 9, 2, 1);
 113       RANLUX_STEP (y1, y3, 10, 3, 2);
 114       RANLUX_STEP (y2, y1, 11, 4, 3);
 115       RANLUX_STEP (y3, y2, 0, 5, 4);
 116       RANLUX_STEP (y1, y3, 1, 6, 5);
 117       RANLUX_STEP (y2, y1, 2, 7, 6);
 118       RANLUX_STEP (y3, y2, 3, 8, 7);
 119       RANLUX_STEP (y1, y3, 4, 9, 8);
 120       RANLUX_STEP (y2, y1, 5, 10, 9);
 121       RANLUX_STEP (y3, y2, 6, 11, 10);
 122
 123       if (y3 < 0)
 124         {
 125           carry = one_bit;
 126           y3 += 1;
 127         }
 128       else
 129         {
 130           carry = 0;
 131         }
 132       xdbl[11] = y3;
 133     }
 134
 135   kmax = state->pr;
 136
 137   for (; k < kmax; ++k)
 138     {
 139       y1 = xdbl[jr] - xdbl[ir];
 140       y2 = y1 - carry;
 141       if (y2 < 0)
 142         {
 143           carry = one_bit;
 144           y2 += 1;
 145         }
 146       else
 147         {
 148           carry = 0;
 149         }
 150       xdbl[ir] = y2;
 151       ydbl[ir] = y2 + shift;
 152       ir = next[ir];
 153       jr = next[jr];
 154     }
 155
 156   ydbl[ir] = xdbl[ir] + shift;
 157
 158   for (k = next[ir]; k > 0;)
 159     {
 160       ydbl[k] = xdbl[k] + shift;
 161       k = next[k];
 162     }
 163
 164   for (k = 0, m = 0; k < 12; ++k)
 165     {
 166       x = xdbl[k];
 167       y2 = ydbl[k] - shift;
 168       if (y2 > x)
 169         y2 -= sone_bit;
 170       y1 = (x - y2) * sbase;
 171
 172       xflt[m++] = (float) y1;
 173       xflt[m++] = (float) y2;
 174     }
 175
 176   state->ir = ir;
 177   state->is = 2 * ir;
 178   state->is_old = 2 * ir;
 179   state->jr = jr;
 180   state->carry = carry;
 181 }
 182
 183
 184 static inline double
 185 ranlxs_get_double (void *vstate)
 186 {
 187   ranlxs_state_t *state = (ranlxs_state_t *) vstate;
 188
 189   const unsigned int is = snext[state->is];
 190
 191   state->is = is;
 192
 193   if (is == state->is_old)
 194     increment_state (state);
 195
 196   return state->xflt[state->is];
 197 }
 198
 199 static unsigned long int
 200 ranlxs_get (void *vstate)
 201 {
 202   return ranlxs_get_double (vstate) * 16777216.0;       /* 2^24 */
 203 }
 204
 205 static void
 206 ranlxs_set_lux (void *vstate, unsigned long int s, unsigned int luxury)
 207 {
 208   ranlxs_state_t *state = (ranlxs_state_t *) vstate;
 209
 210   int ibit, jbit, i, k, m, xbit[31];
 211   double x, y;
 212
 213   long int seed;
 214
 215   if (s == 0)
 216     s = 1;                      /* default seed is 1 */
 217
 218   seed = s;
 219
 220   i = seed & 0xFFFFFFFFUL;
 221
 222   for (k = 0; k < 31; ++k)
 223     {
 224       xbit[k] = i % 2;
 225       i /= 2;
 226     }
 227
 228   ibit = 0;
 229   jbit = 18;
 230
 231   for (k = 0; k < 12; ++k)
 232     {
 233       x = 0;
 234
 235       for (m = 1; m <= 48; ++m)
 236         {
 237           y = (double) xbit[ibit];
 238           x += x + y;
 239           xbit[ibit] = (xbit[ibit] + xbit[jbit]) % 2;
 240           ibit = (ibit + 1) % 31;
 241           jbit = (jbit + 1) % 31;
 242         }
 243       state->xdbl[k] = one_bit * x;
 244     }
 245
 246   state->carry = 0;
 247   state->ir = 0;
 248   state->jr = 7;
 249   state->is = 23;
 250   state->is_old = 0;
 251   state->pr = luxury;
 252 }
 253
 254 static void
 255 ranlxs0_set (void *vstate, unsigned long int s)
 256 {
 257   ranlxs_set_lux (vstate, s, 109);
 258 }
 259
 260 void
 261 ranlxs1_set (void *vstate, unsigned long int s)
 262 {
 263   ranlxs_set_lux (vstate, s, 202);
 264 }
 265
 266 static void
 267 ranlxs2_set (void *vstate, unsigned long int s)
 268 {
 269   ranlxs_set_lux (vstate, s, 397);
 270 }
 271
 272
 273 static const gsl_rng_type ranlxs0_type =
 274 {"ranlxs0",                     /* name */
 275  0x00ffffffUL,                  /* RAND_MAX */
 276  0,                             /* RAND_MIN */
 277  sizeof (ranlxs_state_t),
 278  &ranlxs0_set,
 279  &ranlxs_get,
 280  &ranlxs_get_double};
 281
 282 static const gsl_rng_type ranlxs1_type =
 283 {"ranlxs1",                     /* name */
 284  0x00ffffffUL,                  /* RAND_MAX */
 285  0,                             /* RAND_MIN */
 286  sizeof (ranlxs_state_t),
 287  &ranlxs1_set,
 288  &ranlxs_get,
 289  &ranlxs_get_double};
 290
 291 static const gsl_rng_type ranlxs2_type =
 292 {"ranlxs2",                     /* name */
 293  0x00ffffffUL,                  /* RAND_MAX */
 294  0,                             /* RAND_MIN */
 295  sizeof (ranlxs_state_t),
 296  &ranlxs2_set,
 297  &ranlxs_get,
 298  &ranlxs_get_double};
 299
 300 const gsl_rng_type *gsl_rng_ranlxs0 = &ranlxs0_type;
 301 const gsl_rng_type *gsl_rng_ranlxs1 = &ranlxs1_type;
 302 const gsl_rng_type *gsl_rng_ranlxs2 = &ranlxs2_type;