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

   1 /* rng/ran2.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 the algorithm used in Numerical
  25    Recipe's ran2 generator.  It is a L'Ecuyer combined recursive
  26    generator with a 32-element shuffle-box.
  27
  28    As far as I can tell, in general the effects of adding a shuffle
  29    box cannot be proven theoretically, so the period of this generator
  30    is unknown.
  31
  32    The period of the underlying combined generator is O(2^60). */
  33
  34 static inline unsigned long int ran2_get (void *vstate);
  35 static double ran2_get_double (void *vstate);
  36 static void ran2_set (void *state, unsigned long int s);
  37
  38 static const long int m1 = 2147483563, a1 = 40014, q1 = 53668, r1 = 12211;
  39 static const long int m2 = 2147483399, a2 = 40692, q2 = 52774, r2 = 3791;
  40
  41 #define N_SHUFFLE 32
  42 #define N_DIV (1 + 2147483562/N_SHUFFLE)
  43
  44 typedef struct
  45   {
  46     unsigned long int x;
  47     unsigned long int y;
  48     unsigned long int n;
  49     unsigned long int shuffle[N_SHUFFLE];
  50   }
  51 ran2_state_t;
  52
  53 static inline unsigned long int
  54 ran2_get (void *vstate)
  55 {
  56   ran2_state_t *state = (ran2_state_t *) vstate;
  57
  58   const unsigned long int x = state->x;
  59   const unsigned long int y = state->y;
  60
  61   long int h1 = x / q1;
  62   long int t1 = a1 * (x - h1 * q1) - h1 * r1;
  63
  64   long int h2 = y / q2;
  65   long int t2 = a2 * (y - h2 * q2) - h2 * r2;
  66
  67   if (t1 < 0)
  68     t1 += m1;
  69
  70   if (t2 < 0)
  71     t2 += m2;
  72
  73   state->x = t1;
  74   state->y = t2;
  75
  76   {
  77     unsigned long int j = state->n / N_DIV;
  78     long int delta = state->shuffle[j] - t2;
  79     if (delta < 1)
  80       delta += m1 - 1;
  81     state->n = delta;
  82     state->shuffle[j] = t1;
  83   }
  84
  85   return state->n;
  86 }
  87
  88 static double
  89 ran2_get_double (void *vstate)
  90 {
  91   float x_max = 1 - 1.2e-7f ; /* Numerical Recipes version of 1-FLT_EPS */
  92
  93   float x = ran2_get (vstate) / 2147483563.0f ;
  94
  95   if (x > x_max)
  96     return x_max ;
  97
  98   return x ;
  99 }
 100
 101 static void
 102 ran2_set (void *vstate, unsigned long int s)
 103 {
 104   ran2_state_t *state = (ran2_state_t *) vstate;
 105   int i;
 106
 107   if (s == 0)
 108     s = 1;      /* default seed is 1 */
 109
 110   state->y = s;
 111
 112   for (i = 0; i < 8; i++)
 113     {
 114       long int h = s / q1;
 115       long int t = a1 * (s - h * q1) - h * r1;
 116       if (t < 0)
 117         t += m1;
 118       s = t;
 119     }
 120
 121   for (i = N_SHUFFLE - 1; i >= 0; i--)
 122     {
 123       long int h = s / q1;
 124       long int t = a1 * (s - h * q1) - h * r1;
 125       if (t < 0)
 126         t += m1;
 127       s = t;
 128       state->shuffle[i] = s;
 129     }
 130
 131   state->x = s;
 132   state->n = s;
 133
 134   return;
 135 }
 136
 137 static const gsl_rng_type ran2_type =
 138 {"ran2",                        /* name */
 139  2147483562,                    /* RAND_MAX */
 140  1,                             /* RAND_MIN */
 141  sizeof (ran2_state_t),
 142  &ran2_set,
 143  &ran2_get,
 144  &ran2_get_double};
 145
 146 const gsl_rng_type *gsl_rng_ran2 = &ran2_type;