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

   1 /* rng/ran1.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 ran1 generator.  It is MINSTD with a 32-element
  26    shuffle-box. */
  27
  28 static inline unsigned long int ran1_get (void *vstate);
  29 static double ran1_get_double (void *vstate);
  30 static void ran1_set (void *state, unsigned long int s);
  31
  32 static const long int m = 2147483647, a = 16807, q = 127773, r = 2836;
  33
  34 #define N_SHUFFLE 32
  35 #define N_DIV (1 + 2147483646/N_SHUFFLE)
  36
  37 typedef struct
  38   {
  39     unsigned long int x;
  40     unsigned long int n;
  41     unsigned long int shuffle[N_SHUFFLE];
  42   }
  43 ran1_state_t;
  44
  45 static inline unsigned long int
  46 ran1_get (void *vstate)
  47 {
  48   ran1_state_t *state = (ran1_state_t *) vstate;
  49
  50   const unsigned long int x = state->x;
  51
  52   const long int h = x / q;
  53   const long int t = a * (x - h * q) - h * r;
  54
  55   if (t < 0)
  56     {
  57       state->x = t + m;
  58     }
  59   else
  60     {
  61       state->x = t;
  62     }
  63
  64   {
  65     unsigned long int j = state->n / N_DIV;
  66     state->n = state->shuffle[j];
  67     state->shuffle[j] = state->x;
  68   }
  69
  70   return state->n;
  71 }
  72
  73 static double
  74 ran1_get_double (void *vstate)
  75 {
  76   float x_max = 1 - 1.2e-7f ; /* Numerical Recipes version of 1-FLT_EPS */
  77
  78   float x = ran1_get (vstate) / 2147483647.0f ;
  79
  80   if (x > x_max)
  81     return x_max ;
  82
  83   return x ;
  84 }
  85
  86
  87 static void
  88 ran1_set (void *vstate, unsigned long int s)
  89 {
  90   ran1_state_t *state = (ran1_state_t *) vstate;
  91   int i;
  92
  93   if (s == 0)
  94     s = 1;      /* default seed is 1 */
  95
  96   for (i = 0; i < 8; i++)
  97     {
  98       long int h = s / q;
  99       long int t = a * (s - h * q) - h * r;
 100       if (t < 0)
 101         t += m;
 102       s = t;
 103     }
 104
 105   for (i = N_SHUFFLE - 1; i >= 0; i--)
 106     {
 107       long int h = s / q;
 108       long int t = a * (s - h * q) - h * r;
 109       if (t < 0)
 110         t += m;
 111       s = t;
 112       state->shuffle[i] = s;
 113     }
 114
 115   state->x = s;
 116   state->n = s;
 117
 118   return;
 119 }
 120
 121 static const gsl_rng_type ran1_type =
 122 {"ran1",                        /* name */
 123  2147483646,                    /* RAND_MAX */
 124  1,                             /* RAND_MIN */
 125  sizeof (ran1_state_t),
 126  &ran1_set,
 127  &ran1_get,
 128  &ran1_get_double};
 129
 130 const gsl_rng_type *gsl_rng_ran1 = &ran1_type;