guile/libguile/random.c

/* Copyright (C) 1999,2000,2001, 2003, 2005, 2006, 2009 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
 * the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 * 02110-1301 USA
 */



/* Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */

#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif

#include "libguile/_scm.h"

#include <gmp.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "libguile/smob.h"
#include "libguile/numbers.h"
#include "libguile/feature.h"
#include "libguile/strings.h"
#include "libguile/arrays.h"
#include "libguile/srfi-4.h"
#include "libguile/vectors.h"
#include "libguile/generalized-vectors.h"

#include "libguile/validate.h"
#include "libguile/random.h"


/*
 * A plugin interface for RNGs
 *
 * Using this interface, it is possible for the application to tell
 * libguile to use a different RNG.  This is desirable if it is
 * necessary to use the same RNG everywhere in the application in
 * order to prevent interference, if the application uses RNG
 * hardware, or if the application has special demands on the RNG.
 *
 * Look in random.h and how the default generator is "plugged in" in
 * scm_init_random().
 */

scm_t_rng scm_the_rng;


/*
 * The prepackaged RNG
 *
 * This is the MWC (Multiply With Carry) random number generator
 * described by George Marsaglia at the Department of Statistics and
 * Supercomputer Computations Research Institute, The Florida State
 * University (http://stat.fsu.edu/~geo).
 *
 * It uses 64 bits, has a period of 4578426017172946943 (4.6e18), and
 * passes all tests in the DIEHARD test suite
 * (http://stat.fsu.edu/~geo/diehard.html)
 */

#define A 2131995753UL

#ifndef M_PI
#define M_PI 3.14159265359
#endif

#if SCM_HAVE_T_UINT64

unsigned long
scm_i_uniform32 (scm_t_i_rstate *state)
{
  scm_t_uint64 x = (scm_t_uint64) A * state->w + state->c;
  scm_t_uint32 w = x & 0xffffffffUL;
  state->w = w;
  state->c = x >> 32L;
  return w;
}

#else

/*     ww  This is a portable version of the same RNG without 64 bit
 *   * aa  arithmetic.
 *   ----
 *     xx  It is only intended to provide identical behaviour on
 *    xx   platforms without 8 byte longs or long longs until
 *    xx   someone has implemented the routine in assembler code.
 *   xxcc
 *   ----
 *   ccww
 */

#define L(x) ((x) & 0xffff)
#define H(x) ((x) >> 16)

unsigned long
scm_i_uniform32 (scm_t_i_rstate *state)
{
  scm_t_uint32 x1 = L (A) * L (state->w);
  scm_t_uint32 x2 = L (A) * H (state->w);
  scm_t_uint32 x3 = H (A) * L (state->w);
  scm_t_uint32 w = L (x1) + L (state->c);
  scm_t_uint32 m = H (x1) + L (x2) + L (x3) + H (state->c) + H (w);
  scm_t_uint32 x4 = H (A) * H (state->w);
  state->w = w = (L (m) << 16) + L (w);
  state->c = H (x2) + H (x3) + x4 + H (m);
  return w;
}

#endif

void
scm_i_init_rstate (scm_t_i_rstate *state, const char *seed, int n)
{
  scm_t_uint32 w = 0L;
  scm_t_uint32 c = 0L;
  int i, m;
  for (i = 0; i < n; ++i)
    {
      m = i % 8;
      if (m < 4)
	w += seed[i] << (8 * m);
      else
        c += seed[i] << (8 * (m - 4));
    }
  if ((w == 0 && c == 0) || (w == -1 && c == A - 1))
    ++c;
  state->w = w;
  state->c = c;
}

scm_t_i_rstate *
scm_i_copy_rstate (scm_t_i_rstate *state)
{
  scm_t_rstate *new_state = scm_malloc (scm_the_rng.rstate_size);
  return memcpy (new_state, state, scm_the_rng.rstate_size);
}


/*
 * Random number library functions
 */

scm_t_rstate *
scm_c_make_rstate (const char *seed, int n)
{
  scm_t_rstate *state = scm_malloc (scm_the_rng.rstate_size);
  state->reserved0 = 0;
  scm_the_rng.init_rstate (state, seed, n);
  return state;
}


scm_t_rstate *
scm_c_default_rstate ()
#define FUNC_NAME "scm_c_default_rstate"
{
  SCM state = SCM_VARIABLE_REF (scm_var_random_state);
  if (!SCM_RSTATEP (state))
    SCM_MISC_ERROR ("*random-state* contains bogus random state", SCM_EOL);
  return SCM_RSTATE (state);
}
#undef FUNC_NAME


inline 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) 0xffffffffUL);
}

double
scm_c_normal01 (scm_t_rstate *state)
{
  if (state->reserved0)
    {
      state->reserved0 = 0;
      return state->reserved1;
    }
  else
    {
      double r, a, n;
      
      r = sqrt (-2.0 * log (scm_c_uniform01 (state)));
      a = 2.0 * M_PI * scm_c_uniform01 (state);
      
      n = r * sin (a);
      state->reserved1 = r * cos (a);
      state->reserved0 = 1;
      
      return n;
    }
}

double
scm_c_exp1 (scm_t_rstate *state)
{
  return - log (scm_c_uniform01 (state));
}

unsigned char scm_masktab[256];

unsigned long
scm_c_random (scm_t_rstate *state, unsigned long m)
{
  unsigned int r, mask;
  mask = (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);
  return r;
}

/*
  SCM scm_c_random_bignum (scm_t_rstate *state, SCM m)

  Takes a random state (source of random bits) and a bignum m.
  Returns a bignum b, 0 <= b < m.

  It does this by allocating a bignum b with as many base 65536 digits
  as m, filling b with random bits (in 32 bit chunks) up to the most
  significant 1 in m, and, finally checking if the resultant b is too
  large (>= m).  If too large, we simply repeat the process again.  (It
  is important to throw away all generated random bits if b >= m,
  otherwise we'll end up with a distorted distribution.)

*/

SCM
scm_c_random_bignum (scm_t_rstate *state, SCM m)
{
  SCM result = scm_i_mkbig ();
  const size_t m_bits = mpz_sizeinbase (SCM_I_BIG_MPZ (m), 2);
  /* how many bits would only partially fill the last unsigned long? */
  const size_t end_bits = m_bits % (sizeof (unsigned long) * SCM_CHAR_BIT);
  unsigned long *random_chunks = NULL;
  const unsigned long num_full_chunks =
    m_bits / (sizeof (unsigned long) * SCM_CHAR_BIT);
  const unsigned long num_chunks = num_full_chunks + ((end_bits) ? 1 : 0);

  /* we know the result will be this big */
  mpz_realloc2 (SCM_I_BIG_MPZ (result), m_bits);

  random_chunks =
    (unsigned long *) scm_gc_calloc (num_chunks * sizeof (unsigned long),
                                     "random bignum chunks");

  do
    {
      unsigned long *current_chunk = random_chunks + (num_chunks - 1);
      unsigned long chunks_left = num_chunks;

      mpz_set_ui (SCM_I_BIG_MPZ (result), 0);
      
      if (end_bits)
        {
          /* 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 unsigned long rndbits = scm_the_rng.random_bits (state);
          int rshift = (sizeof (unsigned long) * SCM_CHAR_BIT) - end_bits;
          unsigned long mask = ((unsigned long) ULONG_MAX) >> rshift;
          unsigned long highest_bits = rndbits & mask;
          *current_chunk-- = highest_bits;
          chunks_left--;
        }
      
      while (chunks_left)
        {
          /* now fill in the remaining unsigned long sized chunks */
          *current_chunk-- = scm_the_rng.random_bits (state);
          chunks_left--;
        }
      mpz_import (SCM_I_BIG_MPZ (result),
                  num_chunks,
                  -1,
                  sizeof (unsigned long),
                  0,
                  0,
                  random_chunks);
      /* if result >= m, regenerate it (it is important to regenerate
	 all bits in order not to get a distorted distribution) */
    } while (mpz_cmp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (m)) >= 0);
  scm_gc_free (random_chunks,
               num_chunks * sizeof (unsigned long),
               "random bignum chunks");
  return scm_i_normbig (result);
}

/*
 * Scheme level representation of random states.
 */
 
scm_t_bits scm_tc16_rstate;

static SCM
make_rstate (scm_t_rstate *state)
{
  SCM_RETURN_NEWSMOB (scm_tc16_rstate, state);
}

static size_t
rstate_free (SCM rstate)
{
  free (SCM_RSTATE (rstate));
  return 0;
}

/*
 * Scheme level interface.
 */

SCM_GLOBAL_VARIABLE_INIT (scm_var_random_state, "*random-state*", scm_seed_to_random_state (scm_from_locale_string ("URL:http://stat.fsu.edu/~geo/diehard.html")));

SCM_DEFINE (scm_random, "random", 1, 1, 0, 
            (SCM n, SCM state),
            "Return a number in [0, N).\n"
            "\n"
            "Accepts a positive integer or real n and returns a\n"
            "number of the same type between zero (inclusive) and\n"
            "N (exclusive). The values returned have a uniform\n"
            "distribution.\n"
            "\n"
            "The optional argument @var{state} must be of the type produced\n"
	    "by @code{seed->random-state}. It defaults to the value of the\n"
	    "variable @var{*random-state*}. This object is used to maintain\n"
	    "the state of the pseudo-random-number generator and is altered\n"
	    "as a side effect of the random operation.")
#define FUNC_NAME s_scm_random
{
  if (SCM_UNBNDP (state))
    state = SCM_VARIABLE_REF (scm_var_random_state);
  SCM_VALIDATE_RSTATE (2, state);
  if (SCM_I_INUMP (n))
    {
      unsigned long m = SCM_I_INUM (n);
      SCM_ASSERT_RANGE (1, n, m > 0);
      return scm_from_ulong (scm_c_random (SCM_RSTATE (state), m));
    }
  SCM_VALIDATE_NIM (1, n);
  if (SCM_REALP (n))
    return scm_from_double (SCM_REAL_VALUE (n)
			    * scm_c_uniform01 (SCM_RSTATE (state)));

  if (!SCM_BIGP (n))
    SCM_WRONG_TYPE_ARG (1, n);
  return scm_c_random_bignum (SCM_RSTATE (state), n);
}
#undef FUNC_NAME

SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0, 
            (SCM state),
            "Return a copy of the random state @var{state}.")
#define FUNC_NAME s_scm_copy_random_state
{
  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)));
}
#undef FUNC_NAME

SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0, 
            (SCM seed),
            "Return a new random state using @var{seed}.")
#define FUNC_NAME s_scm_seed_to_random_state
{
  SCM res;
  if (SCM_NUMBERP (seed))
    seed = scm_number_to_string (seed, SCM_UNDEFINED);
  SCM_VALIDATE_STRING (1, seed);
  res = make_rstate (scm_c_make_rstate (scm_i_string_chars (seed),
					scm_i_string_length (seed)));
  scm_remember_upto_here_1 (seed);
  return res;
  
}
#undef FUNC_NAME

SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0, 
            (SCM state),
	    "Return a uniformly distributed inexact real random number in\n"
	    "[0,1).")
#define FUNC_NAME s_scm_random_uniform
{
  if (SCM_UNBNDP (state))
    state = SCM_VARIABLE_REF (scm_var_random_state);
  SCM_VALIDATE_RSTATE (1, state);
  return scm_from_double (scm_c_uniform01 (SCM_RSTATE (state)));
}
#undef FUNC_NAME

SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0, 
            (SCM state),
	    "Return an inexact real in a normal distribution.  The\n"
	    "distribution used has mean 0 and standard deviation 1.  For a\n"
	    "normal distribution with mean m and standard deviation d use\n"
	    "@code{(+ m (* d (random:normal)))}.")
#define FUNC_NAME s_scm_random_normal
{
  if (SCM_UNBNDP (state))
    state = SCM_VARIABLE_REF (scm_var_random_state);
  SCM_VALIDATE_RSTATE (1, state);
  return scm_from_double (scm_c_normal01 (SCM_RSTATE (state)));
}
#undef FUNC_NAME

static void
vector_scale_x (SCM v, double c)
{
  size_t n;
  if (scm_is_simple_vector (v))
    {
      n = SCM_SIMPLE_VECTOR_LENGTH (v);
      while (n-- > 0)
	SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n)) *= c;
    }
  else
    {
      /* must be a f64vector. */
      scm_t_array_handle handle;
      size_t i, len;
      ssize_t inc;
      double *elts;

      elts = scm_f64vector_writable_elements (v, &handle, &len, &inc);

      for (i = 0; i < len; i++, elts += inc)
	*elts *= c;
      
      scm_array_handle_release (&handle);
    }
}

static double
vector_sum_squares (SCM v)
{
  double x, sum = 0.0;
  size_t n;
  if (scm_is_simple_vector (v))
    {
      n = SCM_SIMPLE_VECTOR_LENGTH (v);
      while (n-- > 0)
	{
	  x = SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n));
	  sum += x * x;
	}
    }
  else
    {
      /* must be a f64vector. */
      scm_t_array_handle handle;
      size_t i, len;
      ssize_t inc;
      const double *elts;

      elts = scm_f64vector_elements (v, &handle, &len, &inc);

      for (i = 0; i < len; i++, elts += inc)
	{
	  x = *elts;
	  sum += x * x;
	}

      scm_array_handle_release (&handle);
    }
  return sum;
}

/* For the uniform distribution on the solid sphere, note that in
 * this distribution the length r of the vector has cumulative
 * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be
 * generated as r=u^(1/n).
 */
SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0, 
            (SCM v, SCM state),
	    "Fills @var{vect} with inexact real random numbers the sum of\n"
	    "whose squares is less than 1.0.  Thinking of @var{vect} as\n"
	    "coordinates in space of dimension @var{n} @math{=}\n"
	    "@code{(vector-length @var{vect})}, the coordinates are\n"
	    "uniformly distributed within the unit @var{n}-sphere.")
#define FUNC_NAME s_scm_random_solid_sphere_x
{
  if (SCM_UNBNDP (state))
    state = SCM_VARIABLE_REF (scm_var_random_state);
  SCM_VALIDATE_RSTATE (2, state);
  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))
		  / 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 v, SCM state),
            "Fills vect with inexact real random numbers\n"
            "the sum of whose squares is equal to 1.0.\n"
            "Thinking of vect as coordinates in space of\n"
            "dimension n = (vector-length vect), the coordinates\n"
            "are uniformly distributed over the surface of the\n"
            "unit n-sphere.")
#define FUNC_NAME s_scm_random_hollow_sphere_x
{
  if (SCM_UNBNDP (state))
    state = SCM_VARIABLE_REF (scm_var_random_state);
  SCM_VALIDATE_RSTATE (2, state);
  scm_random_normal_vector_x (v, state);
  vector_scale_x (v, 1 / sqrt (vector_sum_squares (v)));
  return SCM_UNSPECIFIED;
}
#undef FUNC_NAME


SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0, 
            (SCM v, SCM state),
            "Fills vect with inexact real random numbers that are\n"
            "independent and standard normally distributed\n"
            "(i.e., with mean 0 and variance 1).")
#define FUNC_NAME s_scm_random_normal_vector_x
{
  long i;
  scm_t_array_handle handle;
  scm_t_array_dim *dim;

  if (SCM_UNBNDP (state))
    state = SCM_VARIABLE_REF (scm_var_random_state);
  SCM_VALIDATE_RSTATE (2, state);

  scm_generalized_vector_get_handle (v, &handle);
  dim = scm_array_handle_dims (&handle);

  if (scm_is_vector (v))
    {
      SCM *elts = scm_array_handle_writable_elements (&handle);
      for (i = dim->lbnd; i <= dim->ubnd; i++, elts += dim->inc)
	*elts = scm_from_double (scm_c_normal01 (SCM_RSTATE (state)));
    }
  else
    {
      /* must be a f64vector. */
      double *elts = scm_array_handle_f64_writable_elements (&handle);
      for (i = dim->lbnd; i <= dim->ubnd; i++, elts += dim->inc)
	*elts = scm_c_normal01 (SCM_RSTATE (state));
    }

  scm_array_handle_release (&handle);

  return SCM_UNSPECIFIED;
}
#undef FUNC_NAME

SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0, 
            (SCM state),
	    "Return an inexact real in an exponential distribution with mean\n"
	    "1.  For an exponential distribution with mean u use (* u\n"
	    "(random:exp)).")
#define FUNC_NAME s_scm_random_exp
{
  if (SCM_UNBNDP (state))
    state = SCM_VARIABLE_REF (scm_var_random_state);
  SCM_VALIDATE_RSTATE (1, state);
  return scm_from_double (scm_c_exp1 (SCM_RSTATE (state)));
}
#undef FUNC_NAME

void
scm_init_random ()
{
  int i, m;
  /* plug in default RNG */
  scm_t_rng rng =
  {
    sizeof (scm_t_i_rstate),
    (unsigned long (*)()) scm_i_uniform32,
    (void (*)())          scm_i_init_rstate,
    (scm_t_rstate *(*)())    scm_i_copy_rstate
  };
  scm_the_rng = rng;
  
  scm_tc16_rstate = scm_make_smob_type ("random-state", 0);
  scm_set_smob_free (scm_tc16_rstate, rstate_free);

  for (m = 1; m <= 0x100; m <<= 1)
    for (i = m >> 1; i < m; ++i)
      scm_masktab[i] = m - 1;
  
#include "libguile/random.x"

  scm_add_feature ("random");
}

/*
  Local Variables:
  c-file-style: "gnu"
  End:
*/