pintos-os.org Git - pspp/blob - 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 /* Creates a new random number generator, seeds it based on
  44    the current time, and returns it. */
  45 struct rng *
  46 rng_create (void)
  47 {
  48   struct rng *rng;
  49   static time_t t=0;
  50
  51   rng = xmalloc (sizeof *rng);
  52   if (t == 0 || set_seed == NOT_LONG)
  53   {
  54     if (set_seed == NOT_LONG)
  55       time (&t);
  56     else
  57       t = set_seed;
  58   }
  59   else
  60     t++;
  61   rng_seed (rng, &t, sizeof t);
  62   rng->next_normal = NOT_DOUBLE;
  63   return rng;
  64 }
  65
  66 /* Destroys RNG. */
  67 void
  68 rng_destroy (struct rng *rng)
  69 {
  70   free (rng);
  71 }
  72
  73 /* Swap bytes. */
  74 static void
  75 swap_byte (uint8_t *a, uint8_t *b)
  76 {
  77   uint8_t t = *a;
  78   *a = *b;
  79   *b = t;
  80 }
  81
  82 /* Seeds RNG based on the SIZE bytes in BUF.
  83    At most the first 256 bytes of BUF are used. */
  84 void
  85 rng_seed (struct rng *rng, const void *key_, size_t size)
  86 {
  87   const uint8_t *key = key_;
  88   size_t key_idx;
  89   uint8_t *s;
  90   int i, j;
  91
  92   assert (rng != NULL);
  93
  94   s = rng->s;
  95   rng->i = rng->j = 0;
  96   for (i = 0; i < 256; i++)
  97     s[i] = i;
  98   for (key_idx = 0, i = j = 0; i < 256; i++)
  99     {
 100       j = (j + s[i] + key[key_idx]) & 255;
 101       swap_byte (s + i, s + j);
 102       if (++key_idx >= size)
 103         key_idx = 0;
 104     }
 105 }
 106
 107 /* Reads SIZE random bytes from RNG into BUF. */
 108 void
 109 rng_get_bytes (struct rng *rng, void *buf_, size_t size)
 110 {
 111   uint8_t *buf = buf_;
 112   uint8_t *s;
 113   uint8_t i, j;
 114
 115   assert (rng != 0);
 116
 117   s = rng->s;
 118   i = rng->i;
 119   j = rng->j;
 120   while (size-- > 0)
 121     {
 122       i += 1;
 123       j += s[i];
 124       swap_byte (s + i, s + j);
 125       *buf++ = s[(s[i] + s[j]) & 255];
 126     }
 127   rng->i = i;
 128   rng->j = j;
 129 }
 130
 131 /* Returns a random int in the range [0, INT_MAX]. */
 132 int
 133 rng_get_int (struct rng *rng)
 134 {
 135   int value;
 136
 137   do
 138     {
 139       rng_get_bytes (rng, &value, sizeof value);
 140       value = abs (value);
 141     }
 142   while (value < 0);
 143
 144   return value;
 145 }
 146
 147 /* Returns a random unsigned in the range [0, UINT_MAX]. */
 148 unsigned
 149 rng_get_unsigned (struct rng *rng)
 150 {
 151   unsigned value;
 152
 153   rng_get_bytes (rng, &value, sizeof value);
 154   return value;
 155 }
 156
 157 /* Returns a random number from the uniform distribution with
 158    range [0,1). */
 159 double
 160 rng_get_double (struct rng *rng)
 161 {
 162   for (;;)
 163     {
 164       unsigned long ulng;
 165       double dbl;
 166
 167       rng_get_bytes (rng, &ulng, sizeof ulng);
 168       dbl = ulng / (ULONG_MAX + 1.0);
 169       if (dbl >= 0 && dbl < 1)
 170         return dbl;
 171     }
 172 }
 173
 174 /* Returns a random number from the distribution with mean 0 and
 175    standard deviation 1.  (Multiply the result by the desired
 176    standard deviation, then add the desired mean.) */
 177 double
 178 rng_get_double_normal (struct rng *rng)
 179 {
 180   /* Knuth, _The Art of Computer Programming_, Vol. 2, 3.4.1C,
 181      Algorithm P. */
 182   double this_normal;
 183
 184   if (rng->next_normal != NOT_DOUBLE)
 185     {
 186       this_normal = rng->next_normal;
 187       rng->next_normal = NOT_DOUBLE;
 188     }
 189   else
 190     {
 191       double v1, v2, s;
 192
 193       do
 194         {
 195           double u1 = rng_get_double (rng);
 196           double u2 = rng_get_double (rng);
 197           v1 = 2.0 * u1 - 1.0;
 198           v2 = 2.0 * u2 - 1.0;
 199           s = v1 * v1 + v2 * v2;
 200         }
 201       while (s >= 1);
 202
 203       this_normal = v1 * sqrt (-2. * log (s) / s);
 204       rng->next_normal = v2 * sqrt (-2. * log (s) / s);
 205     }
 206
 207   return this_normal;
 208 }
 209
 210 /* Gets an initialized RNG for use in PSPP transformations and
 211    procedures. */
 212 struct rng *
 213 pspp_rng (void)
 214 {
 215   static struct rng *rng;
 216
 217   if (rng == NULL)
 218     rng = rng_create ();
 219   return rng;
 220 }