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

   1 /* rng/ranlxd.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 Martin Luescher's second generation
  25    double-precision (48-bit) version of the RANLUX generator.
  26
  27    Thanks to Martin Luescher for providing information on this
  28    generator.
  29
  30 */
  31
  32 static inline unsigned long int ranlxd_get (void *vstate);
  33 static double ranlxd_get_double (void *vstate);
  34 static void ranlxd_set_lux (void *state, unsigned long int s, unsigned int luxury);
  35 static void ranlxd1_set (void *state, unsigned long int s);
  36 static void ranlxd2_set (void *state, unsigned long int s);
  37
  38 static const int next[12] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0};
  39
  40 static const double one_bit = 1.0 / 281474976710656.0;  /* 1/2^48 */
  41
  42 #define RANLUX_STEP(x1,x2,i1,i2,i3)      \
  43           x1=xdbl[i1] - xdbl[i2];        \
  44           if (x2 < 0)                    \
  45           {                              \
  46             x1-=one_bit;                 \
  47             x2+=1;                       \
  48           }                              \
  49          xdbl[i3]=x2
  50
  51 typedef struct
  52   {
  53     double xdbl[12];
  54     double carry;
  55     unsigned int ir;
  56     unsigned int jr;
  57     unsigned int ir_old;
  58     unsigned int pr;
  59   }
  60 ranlxd_state_t;
  61
  62 static inline void increment_state (ranlxd_state_t * state);
  63
  64 static inline void
  65 increment_state (ranlxd_state_t * state)
  66 {
  67   int k, kmax;
  68   double y1, y2, y3;
  69
  70   double *xdbl = state->xdbl;
  71   double carry = state->carry;
  72   unsigned int ir = state->ir;
  73   unsigned int jr = state->jr;
  74
  75   for (k = 0; ir > 0; ++k)
  76     {
  77       y1 = xdbl[jr] - xdbl[ir];
  78       y2 = y1 - carry;
  79       if (y2 < 0)
  80         {
  81           carry = one_bit;
  82           y2 += 1;
  83         }
  84       else
  85         {
  86           carry = 0;
  87         }
  88       xdbl[ir] = y2;
  89       ir = next[ir];
  90       jr = next[jr];
  91     }
  92
  93   kmax = state->pr - 12;
  94
  95   for (; k <= kmax; k += 12)
  96     {
  97       y1 = xdbl[7] - xdbl[0];
  98       y1 -= carry;
  99
 100       RANLUX_STEP (y2, y1, 8, 1, 0);
 101       RANLUX_STEP (y3, y2, 9, 2, 1);
 102       RANLUX_STEP (y1, y3, 10, 3, 2);
 103       RANLUX_STEP (y2, y1, 11, 4, 3);
 104       RANLUX_STEP (y3, y2, 0, 5, 4);
 105       RANLUX_STEP (y1, y3, 1, 6, 5);
 106       RANLUX_STEP (y2, y1, 2, 7, 6);
 107       RANLUX_STEP (y3, y2, 3, 8, 7);
 108       RANLUX_STEP (y1, y3, 4, 9, 8);
 109       RANLUX_STEP (y2, y1, 5, 10, 9);
 110       RANLUX_STEP (y3, y2, 6, 11, 10);
 111
 112       if (y3 < 0)
 113         {
 114           carry = one_bit;
 115           y3 += 1;
 116         }
 117       else
 118         {
 119           carry = 0;
 120         }
 121       xdbl[11] = y3;
 122     }
 123
 124   kmax = state->pr;
 125
 126   for (; k < kmax; ++k)
 127     {
 128       y1 = xdbl[jr] - xdbl[ir];
 129       y2 = y1 - carry;
 130       if (y2 < 0)
 131         {
 132           carry = one_bit;
 133           y2 += 1;
 134         }
 135       else
 136         {
 137           carry = 0;
 138         }
 139       xdbl[ir] = y2;
 140       ir = next[ir];
 141       jr = next[jr];
 142     }
 143   state->ir = ir;
 144   state->ir_old = ir;
 145   state->jr = jr;
 146   state->carry = carry;
 147 }
 148
 149 static inline unsigned long int
 150 ranlxd_get (void *vstate)
 151 {
 152   return ranlxd_get_double (vstate) * 4294967296.0;     /* 2^32 */
 153 }
 154
 155 static double
 156 ranlxd_get_double (void *vstate)
 157 {
 158   ranlxd_state_t *state = (ranlxd_state_t *) vstate;
 159
 160   int ir = state->ir;
 161
 162   state->ir = next[ir];
 163
 164   if (state->ir == state->ir_old)
 165     increment_state (state);
 166
 167   return state->xdbl[state->ir];
 168 }
 169
 170 static void
 171 ranlxd_set_lux (void *vstate, unsigned long int s, unsigned int luxury)
 172 {
 173   ranlxd_state_t *state = (ranlxd_state_t *) vstate;
 174
 175   int ibit, jbit, i, k, l, xbit[31];
 176   double x, y;
 177
 178   long int seed;
 179
 180   if (s == 0)
 181     s = 1;                      /* default seed is 1 */
 182
 183   seed = s;
 184
 185   i = seed & 0xFFFFFFFFUL;
 186
 187   for (k = 0; k < 31; ++k)
 188     {
 189       xbit[k] = i % 2;
 190       i /= 2;
 191     }
 192
 193   ibit = 0;
 194   jbit = 18;
 195
 196   for (k = 0; k < 12; ++k)
 197     {
 198       x = 0;
 199
 200       for (l = 1; l <= 48; ++l)
 201         {
 202           y = (double) ((xbit[ibit] + 1) % 2);
 203           x += x + y;
 204           xbit[ibit] = (xbit[ibit] + xbit[jbit]) % 2;
 205           ibit = (ibit + 1) % 31;
 206           jbit = (jbit + 1) % 31;
 207         }
 208       state->xdbl[k] = one_bit * x;
 209     }
 210
 211   state->carry = 0;
 212   state->ir = 11;
 213   state->jr = 7;
 214   state->ir_old = 0;
 215   state->pr = luxury;
 216 }
 217
 218 static void
 219 ranlxd1_set (void *vstate, unsigned long int s)
 220 {
 221   ranlxd_set_lux (vstate, s, 202);
 222 }
 223
 224 static void
 225 ranlxd2_set (void *vstate, unsigned long int s)
 226 {
 227   ranlxd_set_lux (vstate, s, 397);
 228 }
 229
 230 static const gsl_rng_type ranlxd1_type =
 231 {"ranlxd1",                     /* name */
 232  0xffffffffUL,                  /* RAND_MAX */
 233  0,                             /* RAND_MIN */
 234  sizeof (ranlxd_state_t),
 235  &ranlxd1_set,
 236  &ranlxd_get,
 237  &ranlxd_get_double};
 238
 239 static const gsl_rng_type ranlxd2_type =
 240 {"ranlxd2",                     /* name */
 241  0xffffffffUL,                  /* RAND_MAX */
 242  0,                             /* RAND_MIN */
 243  sizeof (ranlxd_state_t),
 244  &ranlxd2_set,
 245  &ranlxd_get,
 246  &ranlxd_get_double};
 247
 248 const gsl_rng_type *gsl_rng_ranlxd1 = &ranlxd1_type;
 249 const gsl_rng_type *gsl_rng_ranlxd2 = &ranlxd2_type;