X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Frandom.c;h=e6db94c1501f06f606d8fc1af86c4897dcc8d43f;hb=f98e66bead58bd287261d7bc14dab07515498f55;hp=7f63f477edcf3a28525f2bdb1912aade623e63bd;hpb=4944c86a9318bc5b5578ab145a95c116ffd2c9fd;p=pspp-builds.git

diff --git a/src/random.c b/src/random.c
index 7f63f477..e6db94c1 100644
--- a/src/random.c
+++ b/src/random.c
@@ -18,132 +18,230 @@
    02111-1307, USA. */
 
 #include <config.h>
+#include "random.h"
+#include "error.h"
+#include <inttypes.h>
+#include <limits.h>
 #include <math.h>
 #include <stdlib.h>
 #include <time.h>
+#include "alloc.h"
 #include "magic.h"
-#include "random.h"
 #include "settings.h"
 
-/* Deal with broken system random number generator. */
-#if HAVE_GOOD_RANDOM
-#define real_rand rand
-#define real_srand srand
-#define REAL_RAND_MAX RAND_MAX
-#else /* !HAVE_GOOD_RANDOM */
-#define REAL_RAND_MAX 32767
+/* Random number generator. */
+struct rng 
+  {
+    /* RC4-based random bytes. */
+    uint8_t s[256];
+    uint8_t i, j;
 
-/* Some systems are so broken that they do not supply a value for
-   RAND_MAX.  There is absolutely no reliable way to determine this
-   value, either.  So we must supply our own.  This one is the one
-   presented in the ANSI C standard as strictly compliant. */
-static unsigned long int next = 1;
+    /* Normal distribution. */
+    double next_normal;
+  };
 
-int
-real_rand (void)
+
+/* Return a `random' seed by using the real time clock */
+unsigned long
+random_seed(void)
 {
-  next = next * 1103515245 + 12345;
-  return (unsigned int)(next / 65536) % 32768;
+  time_t t;
+  
+  time(&t);
+
+  return (unsigned long) t;
 }
 
-void
-real_srand (unsigned int seed)
+/* Creates a new random number generator, seeds it based on
+   the current time, and returns it. */
+struct rng *
+rng_create (void) 
 {
-  next = seed;
+  struct rng *rng;
+  static unsigned long seed=0;
+  unsigned long s;
+
+  rng = xmalloc (sizeof *rng);
+
+
+  if ( seed_is_set(&s) ) 
+    {
+      seed = s;
+    }
+  else if ( seed == 0 ) 
+    {
+      seed = random_seed();
+    }
+  assert(seed);
+  /* 
+  if (t == 0 || set_seed_used)
+  {
+    if (set_seed == NOT_LONG) 
+      time (&t);
+    else
+      t = set_seed;
+    set_seed_used=0;
+  }
+  else
+    t++;
+  */
+  rng_seed (rng, &seed, sizeof seed);
+  rng->next_normal = NOT_DOUBLE;
+  return rng;
 }
-#endif /* !HAVE_GOOD_RANDOM */
 
-/* The random number generator here is an implementation in C of
-   Knuth's Algorithm 3.2.2B (Randomizing by Shuffling) in _The Art of
-   Computer Programming_, Vol. 2. */
+/* Destroys RNG. */
+void
+rng_destroy (struct rng *rng) 
+{
+  free (rng);
+}
 
-#define k 13
-static int V[k];
-static int Y;
+/* Swap bytes. */
+static void
+swap_byte (uint8_t *a, uint8_t *b) 
+{
+  uint8_t t = *a;
+  *a = *b;
+  *b = t;
+}
 
-static double X2;
+/* Seeds RNG based on the SIZE bytes in BUF.
+   At most the first 256 bytes of BUF are used. */
+void
+rng_seed (struct rng *rng, const void *key_, size_t size) 
+{
+  const uint8_t *key = key_;
+  size_t key_idx;
+  uint8_t *s;
+  int i, j;
+
+  assert (rng != NULL);
+
+  s = rng->s;
+  rng->i = rng->j = 0;
+  for (i = 0; i < 256; i++) 
+    s[i] = i;
+  for (key_idx = 0, i = j = 0; i < 256; i++) 
+    {
+      j = (j + s[i] + key[key_idx]) & 255;
+      swap_byte (s + i, s + j);
+      if (++key_idx >= size)
+        key_idx = 0;
+    }
+}
 
-/* Initializes the random number generator.  Should be called once by
-   every cmd_*() that uses random numbers.  Note that this includes
-   all procedures that use expressions since they may generate random
-   numbers. */
+/* Reads SIZE random bytes from RNG into BUF. */
 void
-setup_randomize (void)
+rng_get_bytes (struct rng *rng, void *buf_, size_t size) 
 {
-  static time_t curtime;
-  int i;
+  uint8_t *buf = buf_;
+  uint8_t *s;
+  uint8_t i, j;
+
+  assert (rng != 0);
 
-  if (set_seed == NOT_LONG)
+  s = rng->s;
+  i = rng->i;
+  j = rng->j;
+  while (size-- > 0) 
     {
-      if (!curtime)
-	time (&curtime);
-      real_srand (curtime++);
+      i += 1;
+      j += s[i];
+      swap_byte (s + i, s + j);
+      *buf++ = s[(s[i] + s[j]) & 255];
     }
-  else
-    real_srand (set_seed);
+  rng->i = i;
+  rng->j = j;
+}
 
-  set_seed_used = 1;
+/* Returns a random int in the range [0, INT_MAX]. */
+int
+rng_get_int (struct rng *rng) 
+{
+  int value;
 
-  for (i = 0; i < k; i++)
-    V[i] = real_rand ();
-  Y = real_rand ();
-  X2 = NOT_DOUBLE;
+  do 
+    {
+      rng_get_bytes (rng, &value, sizeof value);
+      value = abs (value);
+    }
+  while (value < 0);
+   
+  return value;
 }
 
-/* Standard shuffling procedure for increasing randomness of the ANSI
-   C random number generator. Returns a random number R where 0 <= R
-   <= RAND_MAX. */
-inline int
-shuffle (void)
+/* Returns a random unsigned in the range [0, UINT_MAX]. */
+unsigned
+rng_get_unsigned (struct rng *rng) 
 {
-  int j = k * Y / RAND_MAX;
-  Y = V[j];
-  V[j] = real_rand ();
-  return Y;
+  unsigned value;
+
+  rng_get_bytes (rng, &value, sizeof value);
+  return value;
 }
 
-/* Returns a random number R where 0 <= R <= X. */
-double 
-rand_uniform (double x)
+/* Returns a random number from the uniform distribution with
+   range [0,1). */
+double
+rng_get_double (struct rng *rng) 
 {
-  return ((double) shuffle ()) / (((double) RAND_MAX) / x);
+  for (;;) 
+    {
+      unsigned long ulng;
+      double dbl;
+  
+      rng_get_bytes (rng, &ulng, sizeof ulng);
+      dbl = ulng / (ULONG_MAX + 1.0);
+      if (dbl >= 0 && dbl < 1)
+        return dbl;
+    }
 }
 
 /* Returns a random number from the distribution with mean 0 and
-   standard deviation X.  This uses algorithm P in section 3.4.1C of
-   Knuth's _Art of Computer Programming_, Vol 2. */
+   standard deviation 1.  (Multiply the result by the desired
+   standard deviation, then add the desired mean.) */
 double 
-rand_normal (double x)
+rng_get_double_normal (struct rng *rng)
 {
-  double U1, U2;
-  double V1, V2;
-  double S;
-  double X1;
-
-  if (X2 != NOT_DOUBLE)
+  /* Knuth, _The Art of Computer Programming_, Vol. 2, 3.4.1C,
+     Algorithm P. */
+  double this_normal;
+  
+  if (rng->next_normal != NOT_DOUBLE)
     {
-      double t = X2;
-      X2 = NOT_DOUBLE;
-      return t * x;
+      this_normal = rng->next_normal;
+      rng->next_normal = NOT_DOUBLE;
     }
-  do
+  else 
     {
-      U1 = ((double) shuffle ()) / RAND_MAX;
-      U2 = ((double) shuffle ()) / RAND_MAX;
-      V1 = 2 * U1 - 1;
-      V2 = 2 * U2 - 1;
-      S = V1 * V1 + V2 * V2;
+      double v1, v2, s;
+      
+      do
+        {
+          double u1 = rng_get_double (rng);
+          double u2 = rng_get_double (rng);
+          v1 = 2.0 * u1 - 1.0;
+          v2 = 2.0 * u2 - 1.0;
+          s = v1 * v1 + v2 * v2;
+        }
+      while (s >= 1);
+
+      this_normal = v1 * sqrt (-2. * log (s) / s);
+      rng->next_normal = v2 * sqrt (-2. * log (s) / s); 
     }
-  while (S >= 1);
-  X1 = V1 * sqrt (-2. * log (S) / S);
-  X2 = V2 * sqrt (-2. * log (S) / S);
-  return X1 * x;
+  
+  return this_normal;
 }
 
-/* Returns a random integer R, where 0 <= R < X. */
-int
-rand_simple (int x)
+/* Gets an initialized RNG for use in PSPP transformations and
+   procedures. */
+struct rng *
+pspp_rng (void)
 {
-  return shuffle () % x;
-}
+  static struct rng *rng;
 
+  if (rng == NULL)
+    rng = rng_create ();
+  return rng;
+}