Fix infinite loop in expander
[bpt/guile.git] / libguile / random.c
index 4586b27..1ee0459 100644 (file)
@@ -1,4 +1,6 @@
-/* Copyright (C) 1999,2000,2001, 2003, 2005, 2006, 2009, 2010 Free Software Foundation, Inc.
+/* Copyright (C) 1999, 2000, 2001, 2003, 2005, 2006, 2009, 2010,
+ *    2012, 2013, 2014 Free Software Foundation, Inc.
+ *
  * This library is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public License
  * as published by the Free Software Foundation; either version 3 of
@@ -17,7 +19,7 @@
 
 
 
-/* Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */
+/* Original Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */
 
 #ifdef HAVE_CONFIG_H
 #  include <config.h>
@@ -29,6 +31,9 @@
 #include <stdio.h>
 #include <math.h>
 #include <string.h>
+#include <sys/types.h>
+#include <unistd.h>
+
 #include "libguile/smob.h"
 #include "libguile/numbers.h"
 #include "libguile/feature.h"
@@ -121,9 +126,9 @@ scm_i_copy_rstate (scm_t_rstate *state)
 {
   scm_t_rstate *new_state;
 
-  new_state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size,
+  new_state = scm_gc_malloc_pointerless (state->rng->rstate_size,
                                         "random-state");
-  return memcpy (new_state, state, scm_the_rng.rstate_size);
+  return memcpy (new_state, state, state->rng->rstate_size);
 }
 
 SCM_SYMBOL(scm_i_rstate_tag, "multiply-with-carry");
@@ -169,8 +174,9 @@ scm_c_make_rstate (const char *seed, int n)
 
   state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size,
                                     "random-state");
-  state->reserved0 = 0;
-  scm_the_rng.init_rstate (state, seed, n);
+  state->rng = &scm_the_rng;
+  state->normal_next = 0.0;
+  state->rng->init_rstate (state, seed, n);
   return state;
 }
 
@@ -181,8 +187,9 @@ scm_c_rstate_from_datum (SCM datum)
 
   state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size,
                                     "random-state");
-  state->reserved0 = 0;
-  scm_the_rng.from_datum (state, datum);
+  state->rng = &scm_the_rng;
+  state->normal_next = 0.0;
+  state->rng->from_datum (state, datum);
   return state;
 }
 
@@ -198,21 +205,24 @@ scm_c_default_rstate ()
 #undef FUNC_NAME
 
 
-inline double
+double
 scm_c_uniform01 (scm_t_rstate *state)
 {
-  double x = (double) scm_the_rng.random_bits (state) / (double) 0xffffffffUL;
-  return ((x + (double) scm_the_rng.random_bits (state))
+  double x = (double) state->rng->random_bits (state) / (double) 0xffffffffUL;
+  return ((x + (double) state->rng->random_bits (state))
          / (double) 0xffffffffUL);
 }
 
 double
 scm_c_normal01 (scm_t_rstate *state)
 {
-  if (state->reserved0)
+  if (state->normal_next != 0.0)
     {
-      state->reserved0 = 0;
-      return state->reserved1;
+      double ret = state->normal_next;
+
+      state->normal_next = 0.0;
+
+      return ret;
     }
   else
     {
@@ -222,8 +232,7 @@ scm_c_normal01 (scm_t_rstate *state)
       a = 2.0 * M_PI * scm_c_uniform01 (state);
       
       n = r * sin (a);
-      state->reserved1 = r * cos (a);
-      state->reserved0 = 1;
+      state->normal_next = r * cos (a);
       
       return n;
     }
@@ -237,18 +246,39 @@ scm_c_exp1 (scm_t_rstate *state)
 
 unsigned char scm_masktab[256];
 
-scm_t_uint32
-scm_c_random (scm_t_rstate *state, scm_t_uint32 m)
+static inline scm_t_uint32
+scm_i_mask32 (scm_t_uint32 m)
 {
-  scm_t_uint32 r, mask;
-  mask = (m < 0x100
+  return (m < 0x100
          ? scm_masktab[m]
          : (m < 0x10000
             ? scm_masktab[m >> 8] << 8 | 0xff
             : (m < 0x1000000
                ? scm_masktab[m >> 16] << 16 | 0xffff
-               : scm_masktab[m >> 24] << 24 | 0xffffff)));
-  while ((r = scm_the_rng.random_bits (state) & mask) >= m);
+               : ((scm_t_uint32) scm_masktab[m >> 24]) << 24 | 0xffffff)));
+}
+
+scm_t_uint32
+scm_c_random (scm_t_rstate *state, scm_t_uint32 m)
+{
+  scm_t_uint32 r, mask = scm_i_mask32 (m);
+  while ((r = state->rng->random_bits (state) & mask) >= m);
+  return r;
+}
+
+scm_t_uint64
+scm_c_random64 (scm_t_rstate *state, scm_t_uint64 m)
+{
+  scm_t_uint64 r;
+  scm_t_uint32 mask;
+
+  if (m <= SCM_T_UINT32_MAX)
+    return scm_c_random (state, (scm_t_uint32) m);
+  
+  mask = scm_i_mask32 (m >> 32);
+  while ((r = ((scm_t_uint64) (state->rng->random_bits (state) & mask) << 32)
+          | state->rng->random_bits (state)) >= m)
+    ;
   return r;
 }
 
@@ -297,7 +327,7 @@ scm_c_random_bignum (scm_t_rstate *state, SCM m)
         {
           /* generate a mask with ones in the end_bits position, i.e. if
              end_bits is 3, then we'd have a mask of ...0000000111 */
-          const scm_t_uint32 rndbits = scm_the_rng.random_bits (state);
+          const scm_t_uint32 rndbits = state->rng->random_bits (state);
           int rshift = (sizeof (scm_t_uint32) * SCM_CHAR_BIT) - end_bits;
           scm_t_uint32 mask = ((scm_t_uint32)-1) >> rshift;
           scm_t_uint32 highest_bits = rndbits & mask;
@@ -308,7 +338,7 @@ scm_c_random_bignum (scm_t_rstate *state, SCM m)
       while (chunks_left)
         {
           /* now fill in the remaining scm_t_uint32 sized chunks */
-          *current_chunk-- = scm_the_rng.random_bits (state);
+          *current_chunk-- = state->rng->random_bits (state);
           chunks_left--;
         }
       mpz_import (SCM_I_BIG_MPZ (result),
@@ -367,9 +397,17 @@ SCM_DEFINE (scm_random, "random", 1, 1, 0,
   SCM_VALIDATE_RSTATE (2, state);
   if (SCM_I_INUMP (n))
     {
-      scm_t_uint32 m = SCM_I_INUM (n);
-      SCM_ASSERT_RANGE (1, n, m > 0);
-      return scm_from_uint32 (scm_c_random (SCM_RSTATE (state), m));
+      scm_t_bits m = (scm_t_bits) SCM_I_INUM (n);
+      SCM_ASSERT_RANGE (1, n, SCM_I_INUM (n) > 0);
+#if SCM_SIZEOF_UINTPTR_T <= 4
+      return scm_from_uint32 (scm_c_random (SCM_RSTATE (state),
+                                            (scm_t_uint32) m));
+#elif SCM_SIZEOF_UINTPTR_T <= 8
+      return scm_from_uint64 (scm_c_random64 (SCM_RSTATE (state),
+                                              (scm_t_uint64) m));
+#else
+#error "Cannot deal with this platform's scm_t_bits size"
+#endif
     }
   SCM_VALIDATE_NIM (1, n);
   if (SCM_REALP (n))
@@ -390,7 +428,7 @@ SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0,
   if (SCM_UNBNDP (state))
     state = SCM_VARIABLE_REF (scm_var_random_state);
   SCM_VALIDATE_RSTATE (1, state);
-  return make_rstate (scm_the_rng.copy_rstate (SCM_RSTATE (state)));
+  return make_rstate (SCM_RSTATE (state)->rng->copy_rstate (SCM_RSTATE (state)));
 }
 #undef FUNC_NAME
 
@@ -414,7 +452,7 @@ SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0,
 SCM_DEFINE (scm_datum_to_random_state, "datum->random-state", 1, 0, 0, 
             (SCM datum),
             "Return a new random state using @var{datum}, which should have\n"
-            "been obtailed from @code{random-state->datum}.")
+            "been obtained from @code{random-state->datum}.")
 #define FUNC_NAME s_scm_datum_to_random_state
 {
   return make_rstate (scm_c_rstate_from_datum (datum));
@@ -428,7 +466,7 @@ SCM_DEFINE (scm_random_state_to_datum, "random-state->datum", 1, 0, 0,
 #define FUNC_NAME s_scm_random_state_to_datum
 {
   SCM_VALIDATE_RSTATE (1, state);
-  return scm_the_rng.to_datum (SCM_RSTATE (state));
+  return SCM_RSTATE (state)->rng->to_datum (SCM_RSTATE (state));
 }
 #undef FUNC_NAME
 
@@ -464,7 +502,7 @@ static void
 vector_scale_x (SCM v, double c)
 {
   size_t n;
-  if (scm_is_simple_vector (v))
+  if (scm_is_vector (v))
     {
       n = SCM_SIMPLE_VECTOR_LENGTH (v);
       while (n-- > 0)
@@ -492,7 +530,7 @@ vector_sum_squares (SCM v)
 {
   double x, sum = 0.0;
   size_t n;
-  if (scm_is_simple_vector (v))
+  if (scm_is_vector (v))
     {
       n = SCM_SIMPLE_VECTOR_LENGTH (v);
       while (n-- > 0)
@@ -542,13 +580,13 @@ SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0,
   scm_random_normal_vector_x (v, state);
   vector_scale_x (v,
                  pow (scm_c_uniform01 (SCM_RSTATE (state)),
-                      1.0 / scm_c_generalized_vector_length (v))
+                      1.0 / scm_c_array_length (v))
                  / sqrt (vector_sum_squares (v)));
   return SCM_UNSPECIFIED;
 }
 #undef FUNC_NAME
 
-SCM_DEFINE (scm_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0, 
+SCM_DEFINE (scm_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0,
             (SCM v, SCM state),
             "Fills vect with inexact real random numbers\n"
             "the sum of whose squares is equal to 1.0.\n"
@@ -586,7 +624,7 @@ SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0,
   scm_generalized_vector_get_handle (v, &handle);
   dim = scm_array_handle_dims (&handle);
 
-  if (scm_is_vector (v))
+  if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_SCM)
     {
       SCM *elts = scm_array_handle_writable_elements (&handle);
       for (i = dim->lbnd; i <= dim->ubnd; i++, elts += dim->inc)
@@ -620,6 +658,108 @@ SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0,
 }
 #undef FUNC_NAME
 
+/* Return a new random-state seeded from the time, date, process ID, an
+   address from a freshly allocated heap cell, an address from the local
+   stack frame, and a high-resolution timer if available.  This is only
+   to be used as a last resort, when no better source of entropy is
+   available. */
+static SCM
+random_state_of_last_resort (void)
+{
+  SCM state;
+  SCM time_of_day = scm_gettimeofday ();
+  SCM sources = scm_list_n
+    (scm_from_unsigned_integer (SCM_UNPACK (time_of_day)),  /* heap addr */
+     /* Avoid scm_getpid, since it depends on HAVE_POSIX. */
+     scm_from_unsigned_integer (getpid ()),                 /* process ID */
+     scm_get_internal_real_time (), /* high-resolution process timer */
+     scm_from_unsigned_integer ((scm_t_bits) &time_of_day), /* stack addr */
+     scm_car (time_of_day), /* seconds since midnight 1970-01-01 UTC */
+     scm_cdr (time_of_day), /* microsecond component of the above clock */
+     SCM_UNDEFINED);
+
+  /* Concatenate the sources bitwise to form the seed */
+  SCM seed = SCM_INUM0;
+  while (scm_is_pair (sources))
+    {
+      seed = scm_logxor (seed, scm_ash (scm_car (sources),
+                                        scm_integer_length (seed)));
+      sources = scm_cdr (sources);
+    }
+
+  /* FIXME The following code belongs in `scm_seed_to_random_state',
+     and here we should simply do:
+
+       return scm_seed_to_random_state (seed);
+
+     Unfortunately, `scm_seed_to_random_state' only preserves around 32
+     bits of entropy from the provided seed.  I don't know if it's okay
+     to fix that in 2.0, so for now we have this workaround. */
+  {
+    int i, len;
+    unsigned char *buf;
+    len = scm_to_int (scm_ceiling_quotient (scm_integer_length (seed),
+                                            SCM_I_MAKINUM (8)));
+    buf = (unsigned char *) malloc (len);
+    for (i = len-1; i >= 0; --i)
+      {
+        buf[i] = scm_to_int (scm_logand (seed, SCM_I_MAKINUM (255)));
+        seed = scm_ash (seed, SCM_I_MAKINUM (-8));
+      }
+    state = make_rstate (scm_c_make_rstate ((char *) buf, len));
+    free (buf);
+  }
+  return state;
+}
+
+/* Attempt to fill buffer with random bytes from /dev/urandom.
+   Return 1 if successful, else return 0. */
+static int
+read_dev_urandom (unsigned char *buf, size_t len)
+{
+  size_t res = 0;
+  FILE *f = fopen ("/dev/urandom", "r");
+  if (f)
+    {
+      res = fread(buf, 1, len, f);
+      fclose (f);
+    }
+  return (res == len);
+}
+
+/* Fill a buffer with random bytes seeded from a platform-specific
+   source of entropy.  /dev/urandom is used if available.  Note that
+   this function provides no guarantees about the amount of entropy
+   present in the returned bytes. */
+void
+scm_i_random_bytes_from_platform (unsigned char *buf, size_t len)
+{
+  if (read_dev_urandom (buf, len))
+    return;
+  else  /* FIXME: support other platform sources */
+    {
+      /* When all else fails, use this (rather weak) fallback */
+      SCM random_state = random_state_of_last_resort ();
+      int i;
+      for (i = len-1; i >= 0; --i)
+        buf[i] = scm_to_int (scm_random (SCM_I_MAKINUM (256), random_state));
+    }
+}
+
+SCM_DEFINE (scm_random_state_from_platform, "random-state-from-platform", 0, 0, 0,
+            (void),
+            "Construct a new random state seeded from a platform-specific\n\
+source of entropy, appropriate for use in non-security-critical applications.")
+#define FUNC_NAME s_scm_random_state_from_platform
+{
+  unsigned char buf[32];
+  if (read_dev_urandom (buf, sizeof(buf)))
+    return make_rstate (scm_c_make_rstate ((char *) buf, sizeof(buf)));
+  else
+    return random_state_of_last_resort ();
+}
+#undef FUNC_NAME
+
 void
 scm_init_random ()
 {