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

   1 /* rng/ran3.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 Knuths's
  25    subtractive generator, with the Numerical Recipe's ran3 paramters.
  26    It is a subtractive lagged fibonnaci generator. */
  27
  28 static inline unsigned long int ran3_get (void *vstate);
  29 static double ran3_get_double (void *vstate);
  30 static void ran3_set (void *state, unsigned long int s);
  31
  32 #define M_BIG 1000000000
  33 #define M_SEED 161803398
  34
  35 typedef struct
  36   {
  37     unsigned int x;
  38     unsigned int y;
  39     unsigned long int buffer[56];
  40   }
  41 ran3_state_t;
  42
  43 static inline unsigned long int
  44 ran3_get (void *vstate)
  45 {
  46   ran3_state_t *state = (ran3_state_t *) vstate;
  47   long int j;
  48
  49   state->x++;
  50
  51   if (state->x == 56)
  52     state->x = 1;
  53
  54   state->y++;
  55
  56   if (state->y == 56)
  57     state->y = 1;
  58
  59   j = state->buffer[state->x] - state->buffer[state->y];
  60
  61   if (j < 0)
  62     j += M_BIG;
  63
  64   state->buffer[state->x] = j;
  65
  66   return j;
  67 }
  68
  69 static double
  70 ran3_get_double (void *vstate)
  71 {
  72   return ran3_get (vstate) / (double) M_BIG ;
  73 }
  74
  75 static void
  76 ran3_set (void *vstate, unsigned long int s)
  77 {
  78   ran3_state_t *state = (ran3_state_t *) vstate;
  79   int i, i1;
  80   long int j, k;
  81
  82   if (s == 0)
  83     s = 1;      /* default seed is 1 */
  84
  85   j = (M_SEED - s) % M_BIG;
  86
  87   /* the zeroth element is never used, but we initialize it for
  88      consistency between states */
  89
  90   state->buffer[0] = 0;
  91
  92   state->buffer[55] = j;
  93
  94   k = 1;
  95   for (i = 1; i < 55; i++)
  96     {
  97       int n = (21 * i) % 55;
  98       state->buffer[n] = k;
  99       k = j - k;
 100       if (k < 0)
 101         k += M_BIG;
 102       j = state->buffer[n];
 103
 104     }
 105
 106   for (i1 = 0; i1 < 4; i1++)
 107     {
 108       for (i = 1; i < 56; i++)
 109         {
 110           long int t = state->buffer[i] - state->buffer[1 + (i + 30) % 55];
 111           if (t < 0)
 112             t += M_BIG;
 113           state->buffer[i] = t;
 114         }
 115     }
 116
 117   state->x = 0;
 118   state->y = 31;
 119
 120   return;
 121 }
 122
 123 static const gsl_rng_type ran3_type =
 124 {"ran3",                        /* name */
 125  M_BIG,                         /* RAND_MAX */
 126  0,                             /* RAND_MIN */
 127  sizeof (ran3_state_t),
 128  &ran3_set,
 129  &ran3_get,
 130  &ran3_get_double};
 131
 132 const gsl_rng_type *gsl_rng_ran3 = &ran3_type;