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

   1 /* rng/zuf.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 /* It is crucial that m == n-273 mod 607 at all times;
  25    For speed of execution, however, this is never enforced.
  26    Instead is is set in the initializer: note 607-273=334
  27    Note also that the state.u[607] is not initialized */
  28
  29 static inline unsigned long int zuf_get (void *vstate);
  30 static double zuf_get_double (void *vstate);
  31 static void zuf_set (void *state, unsigned long int s);
  32
  33 static const unsigned long int zuf_randmax = 16777216;  /* 2^24 */
  34
  35 typedef struct
  36   {
  37     int n;
  38     unsigned long int u[607];
  39   }
  40 zuf_state_t;
  41
  42 /* The zufall package was implemented with float's, which is to say 24
  43    bits of precision.  Since I'm using long's instead, my RANDMAX
  44    reflects this. */
  45
  46 static inline unsigned long int
  47 zuf_get (void *vstate)
  48 {
  49   zuf_state_t *state = (zuf_state_t *) vstate;
  50   const int n = state->n;
  51   const int m = (n - 273 + 607) % 607;
  52   unsigned long int t = state->u[n] + state->u[m];
  53
  54   while (t > zuf_randmax)
  55     t -= zuf_randmax;
  56
  57   state->u[n] = t;
  58
  59   if (n == 606)
  60     {
  61       state->n = 0;
  62     }
  63   else
  64     {
  65       state->n = n + 1;
  66     }
  67
  68   return t;
  69 }
  70
  71 static double
  72 zuf_get_double (void *vstate)
  73 {
  74   return zuf_get (vstate) / 16777216.0 ;
  75 }
  76
  77 static void
  78 zuf_set (void *vstate, unsigned long int s)
  79 {
  80   /* A very elaborate seeding procedure is provided with the
  81      zufall package; this is virtually a copy of that procedure */
  82
  83   /* Initialized data */
  84
  85   long int kl = 9373;
  86   long int ij = 1802;
  87
  88   /* Local variables */
  89   long int i, j, k, l, m;
  90   double x, y;
  91   long int ii, jj;
  92
  93   zuf_state_t *state = (zuf_state_t *) vstate;
  94
  95   state->n = 0;
  96
  97 /*  generates initial seed buffer by linear congruential */
  98 /*  method. Taken from Marsaglia, FSU report FSU-SCRI-87-50 */
  99 /*  variable seed should be 0 < seed <31328 */
 100
 101   if (s == 0)
 102     s = 1802;   /* default seed is 1802 */
 103
 104   ij = s;
 105
 106   i = ij / 177 % 177 + 2;
 107   j = ij % 177 + 2;
 108   k = kl / 169 % 178 + 1;
 109   l = kl % 169;
 110   for (ii = 0; ii < 607; ++ii)
 111     {
 112       x = 0.0;
 113       y = 0.5;
 114       /* 24 bits?? */
 115       for (jj = 1; jj <= 24; ++jj)
 116         {
 117           m = i * j % 179 * k % 179;
 118           i = j;
 119           j = k;
 120           k = m;
 121           l = (l * 53 + 1) % 169;
 122           if (l * m % 64 >= 32)
 123             {
 124               x += y;
 125             }
 126           y *= 0.5;
 127         }
 128       state->u[ii] = (unsigned long int) (x * zuf_randmax);
 129     }
 130 }
 131
 132 static const gsl_rng_type zuf_type =
 133 {"zuf",                         /* name */
 134  0x00ffffffUL,                  /* RAND_MAX */
 135  0,                             /* RAND_MIN */
 136  sizeof (zuf_state_t),
 137  &zuf_set,
 138  &zuf_get,
 139  &zuf_get_double};
 140
 141 const gsl_rng_type *gsl_rng_zuf = &zuf_type;