diff options
Diffstat (limited to 'libguile/random.c')
-rw-r--r-- | libguile/random.c | 809 |
1 files changed, 809 insertions, 0 deletions
diff --git a/libguile/random.c b/libguile/random.c new file mode 100644 index 000000000..df592a972 --- /dev/null +++ b/libguile/random.c @@ -0,0 +1,809 @@ +/* Copyright 1999-2001,2003,2005-2006,2009-2010,2012-2014,2017-2018 + Free Software Foundation, Inc. + + This file is part of Guile. + + Guile 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. + + Guile 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 Guile. If not, see + <https://www.gnu.org/licenses/>. */ + + +/* Original Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */ + +#ifdef HAVE_CONFIG_H +# include <config.h> +#endif + +#include <math.h> +#include <stdio.h> +#include <string.h> +#include <sys/types.h> +#include <unistd.h> + +#include <gmp.h> + +#include "arrays.h" +#include "feature.h" +#include "generalized-arrays.h" +#include "generalized-vectors.h" +#include "gsubr.h" +#include "list.h" +#include "modules.h" +#include "numbers.h" +#include "numbers.h" +#include "pairs.h" +#include "smob.h" +#include "srfi-4.h" +#include "stime.h" +#include "strings.h" +#include "symbols.h" +#include "variable.h" +#include "vectors.h" + +#include "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) + */ + +typedef struct scm_t_i_rstate { + scm_t_rstate rstate; + uint32_t w; + uint32_t c; +} scm_t_i_rstate; + + +#define A 2131995753UL + +#ifndef M_PI +#define M_PI 3.14159265359 +#endif + +static uint32_t +scm_i_uniform32 (scm_t_rstate *state) +{ + scm_t_i_rstate *istate = (scm_t_i_rstate*) state; + uint64_t x = (uint64_t) A * istate->w + istate->c; + uint32_t w = x & 0xffffffffUL; + istate->w = w; + istate->c = x >> 32L; + return w; +} + +static void +scm_i_init_rstate (scm_t_rstate *state, const char *seed, int n) +{ + scm_t_i_rstate *istate = (scm_t_i_rstate*) state; + uint32_t w = 0L; + uint32_t 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; + istate->w = w; + istate->c = c; +} + +static scm_t_rstate * +scm_i_copy_rstate (scm_t_rstate *state) +{ + scm_t_rstate *new_state; + + new_state = scm_gc_malloc_pointerless (state->rng->rstate_size, + "random-state"); + return memcpy (new_state, state, state->rng->rstate_size); +} + +SCM_SYMBOL(scm_i_rstate_tag, "multiply-with-carry"); + +static void +scm_i_rstate_from_datum (scm_t_rstate *state, SCM value) +#define FUNC_NAME "scm_i_rstate_from_datum" +{ + scm_t_i_rstate *istate = (scm_t_i_rstate*) state; + uint32_t w, c; + long length; + + SCM_VALIDATE_LIST_COPYLEN (SCM_ARG1, value, length); + SCM_ASSERT (length == 3, value, SCM_ARG1, FUNC_NAME); + SCM_ASSERT (scm_is_eq (SCM_CAR (value), scm_i_rstate_tag), + value, SCM_ARG1, FUNC_NAME); + SCM_VALIDATE_UINT_COPY (SCM_ARG1, SCM_CADR (value), w); + SCM_VALIDATE_UINT_COPY (SCM_ARG1, SCM_CADDR (value), c); + + istate->w = w; + istate->c = c; +} +#undef FUNC_NAME + +static SCM +scm_i_rstate_to_datum (scm_t_rstate *state) +{ + scm_t_i_rstate *istate = (scm_t_i_rstate*) state; + return scm_list_3 (scm_i_rstate_tag, + scm_from_uint32 (istate->w), + scm_from_uint32 (istate->c)); +} + + +/* + * Random number library functions + */ + +scm_t_rstate * +scm_c_make_rstate (const char *seed, int n) +{ + scm_t_rstate *state; + + state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size, + "random-state"); + state->rng = &scm_the_rng; + state->normal_next = 0.0; + state->rng->init_rstate (state, seed, n); + return state; +} + +scm_t_rstate * +scm_c_rstate_from_datum (SCM datum) +{ + scm_t_rstate *state; + + state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size, + "random-state"); + state->rng = &scm_the_rng; + state->normal_next = 0.0; + state->rng->from_datum (state, datum); + 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 + + +double +scm_c_uniform01 (scm_t_rstate *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->normal_next != 0.0) + { + double ret = state->normal_next; + + state->normal_next = 0.0; + + return ret; + } + 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->normal_next = r * cos (a); + + return n; + } +} + +double +scm_c_exp1 (scm_t_rstate *state) +{ + return - log (scm_c_uniform01 (state)); +} + +unsigned char scm_masktab[256]; + +static inline uint32_t +scm_i_mask32 (uint32_t m) +{ + return (m < 0x100 + ? scm_masktab[m] + : (m < 0x10000 + ? scm_masktab[m >> 8] << 8 | 0xff + : (m < 0x1000000 + ? scm_masktab[m >> 16] << 16 | 0xffff + : ((uint32_t) scm_masktab[m >> 24]) << 24 | 0xffffff))); +} + +uint32_t +scm_c_random (scm_t_rstate *state, uint32_t m) +{ + uint32_t r, mask = scm_i_mask32 (m); + while ((r = state->rng->random_bits (state) & mask) >= m); + return r; +} + +uint64_t +scm_c_random64 (scm_t_rstate *state, uint64_t m) +{ + uint64_t r; + uint32_t mask; + + if (m <= UINT32_MAX) + return scm_c_random (state, (uint32_t) m); + + mask = scm_i_mask32 (m >> 32); + while ((r = ((uint64_t) (state->rng->random_bits (state) & mask) << 32) + | state->rng->random_bits (state)) >= 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 uint32_t? */ + const size_t end_bits = m_bits % (sizeof (uint32_t) * SCM_CHAR_BIT); + uint32_t *random_chunks = NULL; + const uint32_t num_full_chunks = + m_bits / (sizeof (uint32_t) * SCM_CHAR_BIT); + const uint32_t 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 = + (uint32_t *) scm_gc_calloc (num_chunks * sizeof (uint32_t), + "random bignum chunks"); + + do + { + uint32_t *current_chunk = random_chunks + (num_chunks - 1); + uint32_t 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 uint32_t rndbits = state->rng->random_bits (state); + int rshift = (sizeof (uint32_t) * SCM_CHAR_BIT) - end_bits; + uint32_t mask = ((uint32_t)-1) >> rshift; + uint32_t highest_bits = rndbits & mask; + *current_chunk-- = highest_bits; + chunks_left--; + } + + while (chunks_left) + { + /* now fill in the remaining uint32_t sized chunks */ + *current_chunk-- = state->rng->random_bits (state); + chunks_left--; + } + mpz_import (SCM_I_BIG_MPZ (result), + num_chunks, + -1, + sizeof (uint32_t), + 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 (uint32_t), + "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); +} + + +/* + * 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)) + { + 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), + (uint32_t) m)); +#elif SCM_SIZEOF_UINTPTR_T <= 8 + return scm_from_uint64 (scm_c_random64 (SCM_RSTATE (state), + (uint64_t) m)); +#else +#error "Cannot deal with this platform's scm_t_bits size" +#endif + } + 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_RSTATE (state)->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_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 obtained from @code{random-state->datum}.") +#define FUNC_NAME s_scm_datum_to_random_state +{ + return make_rstate (scm_c_rstate_from_datum (datum)); +} +#undef FUNC_NAME + +SCM_DEFINE (scm_random_state_to_datum, "random-state->datum", 1, 0, 0, + (SCM state), + "Return a datum representation of @var{state} that may be\n" + "written out and read back with the Scheme reader.") +#define FUNC_NAME s_scm_random_state_to_datum +{ + SCM_VALIDATE_RSTATE (1, state); + return SCM_RSTATE (state)->rng->to_datum (SCM_RSTATE (state)); +} +#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) +{ + scm_t_array_handle handle; + scm_t_array_dim const * dims; + ssize_t i, inc, ubnd; + + scm_array_get_handle (v, &handle); + dims = scm_array_handle_dims (&handle); + if (1 == scm_array_handle_rank (&handle)) + { + ubnd = dims[0].ubnd; + inc = dims[0].inc; + + if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_F64) + { + double *elts = (double *)(handle.writable_elements) + handle.base; + for (i = dims[0].lbnd; i <= ubnd; ++i, elts += inc) + *elts *= c; + return; + } + else if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_SCM) + { + SCM *elts = (SCM *)(handle.writable_elements) + handle.base; + for (i = dims[0].lbnd; i <= ubnd; ++i, elts += inc) + SCM_REAL_VALUE (*elts) *= c; + return; + } + } + scm_array_handle_release (&handle); + scm_misc_error (NULL, "must be a rank-1 array of type #t or 'f64", scm_list_1 (v)); +} + +static double +vector_sum_squares (SCM v) +{ + double x, sum = 0.0; + scm_t_array_handle handle; + scm_t_array_dim const * dims; + ssize_t i, inc, ubnd; + + scm_array_get_handle (v, &handle); + dims = scm_array_handle_dims (&handle); + if (1 == scm_array_handle_rank (&handle)) + { + ubnd = dims[0].ubnd; + inc = dims[0].inc; + if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_F64) + { + const double *elts = (const double *)(handle.elements) + handle.base; + for (i = dims[0].lbnd; i <= ubnd; ++i, elts += inc) + { + x = *elts; + sum += x * x; + } + return sum; + } + else if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_SCM) + { + const SCM *elts = (const SCM *)(handle.elements) + handle.base; + for (i = dims[0].lbnd; i <= ubnd; ++i, elts += inc) + { + x = SCM_REAL_VALUE (*elts); + sum += x * x; + } + return sum; + } + } + scm_array_handle_release (&handle); + scm_misc_error (NULL, "must be an array of type #t or 'f64", scm_list_1 (v)); +} + +/* 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_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 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 +{ + scm_t_array_handle handle; + scm_t_array_dim const * dims; + ssize_t i; + + if (SCM_UNBNDP (state)) + state = SCM_VARIABLE_REF (scm_var_random_state); + SCM_VALIDATE_RSTATE (2, state); + + scm_array_get_handle (v, &handle); + if (1 != scm_array_handle_rank (&handle)) + { + scm_array_handle_release (&handle); + scm_wrong_type_arg_msg (NULL, 0, v, "rank 1 array"); + } + + dims = scm_array_handle_dims (&handle); + + if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_SCM) + { + SCM *elts = scm_array_handle_writable_elements (&handle); + for (i = dims->lbnd; i <= dims->ubnd; i++, elts += dims->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 = dims->lbnd; i <= dims->ubnd; i++, elts += dims->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 + +/* 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 () +{ + int i, m; + /* plug in default RNG */ + scm_t_rng rng = + { + sizeof (scm_t_i_rstate), + scm_i_uniform32, + scm_i_init_rstate, + scm_i_copy_rstate, + scm_i_rstate_from_datum, + scm_i_rstate_to_datum + }; + scm_the_rng = rng; + + scm_tc16_rstate = scm_make_smob_type ("random-state", 0); + + for (m = 1; m <= 0x100; m <<= 1) + for (i = m >> 1; i < m; ++i) + scm_masktab[i] = m - 1; + +#include "random.x" + + scm_add_feature ("random"); +} |