src/random.c

   1 /* PSPP - computes sample statistics.
   2    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
   3    Written by Ben Pfaff <blp@gnu.org>.
   4
   5    This program is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU General Public License as
   7    published by the Free Software Foundation; either version 2 of the
   8    License, or (at 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., 59 Temple Place - Suite 330, Boston, MA
  18    02111-1307, USA. */
  19
  20 #include <config.h>
  21 #include "random.h"
  22 #include <assert.h>
  23 #include <inttypes.h>
  24 #include <limits.h>
  25 #include <math.h>
  26 #include <stdlib.h>
  27 #include <time.h>
  28 #include "alloc.h"
  29 #include "magic.h"
  30 #include "settings.h"
  31
  32 /* Random number generator. */
  33 struct rng
  34   {
  35     /* RC4-based random bytes. */
  36     uint8_t s[256];
  37     uint8_t i, j;
  38
  39     /* Normal distribution. */
  40     double next_normal;
  41   };
  42
  43
  44 /* Return a `random' seed by using the real time clock */
  45 unsigned long
  46 random_seed(void)
  47 {
  48   time_t t;
  49
  50   time(&t);
  51
  52   return (unsigned long) t;
  53 }
  54
  55 /* Creates a new random number generator, seeds it based on
  56    the current time, and returns it. */
  57 struct rng *
  58 rng_create (void)
  59 {
  60   struct rng *rng;
  61   static unsigned long seed=0;
  62   unsigned long s;
  63
  64   rng = xmalloc (sizeof *rng);
  65
  66
  67   if ( seed_is_set(&s) )
  68     {
  69       seed = s;
  70     }
  71   else if ( seed == 0 )
  72     {
  73       seed = random_seed();
  74     }
  75   assert(seed);
  76   /*
  77   if (t == 0 || set_seed_used)
  78   {
  79     if (set_seed == NOT_LONG)
  80       time (&t);
  81     else
  82       t = set_seed;
  83     set_seed_used=0;
  84   }
  85   else
  86     t++;
  87   */
  88   rng_seed (rng, &seed, sizeof seed);
  89   rng->next_normal = NOT_DOUBLE;
  90   return rng;
  91 }
  92
  93 /* Destroys RNG. */
  94 void
  95 rng_destroy (struct rng *rng)
  96 {
  97   free (rng);
  98 }
  99
 100 /* Swap bytes. */
 101 static void
 102 swap_byte (uint8_t *a, uint8_t *b)
 103 {
 104   uint8_t t = *a;
 105   *a = *b;
 106   *b = t;
 107 }
 108
 109 /* Seeds RNG based on the SIZE bytes in BUF.
 110    At most the first 256 bytes of BUF are used. */
 111 void
 112 rng_seed (struct rng *rng, const void *key_, size_t size)
 113 {
 114   const uint8_t *key = key_;
 115   size_t key_idx;
 116   uint8_t *s;
 117   int i, j;
 118
 119   assert (rng != NULL);
 120
 121   s = rng->s;
 122   rng->i = rng->j = 0;
 123   for (i = 0; i < 256; i++)
 124     s[i] = i;
 125   for (key_idx = 0, i = j = 0; i < 256; i++)
 126     {
 127       j = (j + s[i] + key[key_idx]) & 255;
 128       swap_byte (s + i, s + j);
 129       if (++key_idx >= size)
 130         key_idx = 0;
 131     }
 132 }
 133
 134 /* Reads SIZE random bytes from RNG into BUF. */
 135 void
 136 rng_get_bytes (struct rng *rng, void *buf_, size_t size)
 137 {
 138   uint8_t *buf = buf_;
 139   uint8_t *s;
 140   uint8_t i, j;
 141
 142   assert (rng != 0);
 143
 144   s = rng->s;
 145   i = rng->i;
 146   j = rng->j;
 147   while (size-- > 0)
 148     {
 149       i += 1;
 150       j += s[i];
 151       swap_byte (s + i, s + j);
 152       *buf++ = s[(s[i] + s[j]) & 255];
 153     }
 154   rng->i = i;
 155   rng->j = j;
 156 }
 157
 158 /* Returns a random int in the range [0, INT_MAX]. */
 159 int
 160 rng_get_int (struct rng *rng)
 161 {
 162   int value;
 163
 164   do
 165     {
 166       rng_get_bytes (rng, &value, sizeof value);
 167       value = abs (value);
 168     }
 169   while (value < 0);
 170
 171   return value;
 172 }
 173
 174 /* Returns a random unsigned in the range [0, UINT_MAX]. */
 175 unsigned
 176 rng_get_unsigned (struct rng *rng)
 177 {
 178   unsigned value;
 179
 180   rng_get_bytes (rng, &value, sizeof value);
 181   return value;
 182 }
 183
 184 /* Returns a random number from the uniform distribution with
 185    range [0,1). */
 186 double
 187 rng_get_double (struct rng *rng)
 188 {
 189   for (;;)
 190     {
 191       unsigned long ulng;
 192       double dbl;
 193
 194       rng_get_bytes (rng, &ulng, sizeof ulng);
 195       dbl = ulng / (ULONG_MAX + 1.0);
 196       if (dbl >= 0 && dbl < 1)
 197         return dbl;
 198     }
 199 }
 200
 201 /* Returns a random number from the distribution with mean 0 and
 202    standard deviation 1.  (Multiply the result by the desired
 203    standard deviation, then add the desired mean.) */
 204 double
 205 rng_get_double_normal (struct rng *rng)
 206 {
 207   /* Knuth, _The Art of Computer Programming_, Vol. 2, 3.4.1C,
 208      Algorithm P. */
 209   double this_normal;
 210
 211   if (rng->next_normal != NOT_DOUBLE)
 212     {
 213       this_normal = rng->next_normal;
 214       rng->next_normal = NOT_DOUBLE;
 215     }
 216   else
 217     {
 218       double v1, v2, s;
 219
 220       do
 221         {
 222           double u1 = rng_get_double (rng);
 223           double u2 = rng_get_double (rng);
 224           v1 = 2.0 * u1 - 1.0;
 225           v2 = 2.0 * u2 - 1.0;
 226           s = v1 * v1 + v2 * v2;
 227         }
 228       while (s >= 1);
 229
 230       this_normal = v1 * sqrt (-2. * log (s) / s);
 231       rng->next_normal = v2 * sqrt (-2. * log (s) / s);
 232     }
 233
 234   return this_normal;
 235 }
 236
 237 /* Gets an initialized RNG for use in PSPP transformations and
 238    procedures. */
 239 struct rng *
 240 pspp_rng (void)
 241 {
 242   static struct rng *rng;
 243
 244   if (rng == NULL)
 245     rng = rng_create ();
 246   return rng;
 247 }