X-Git-Url: https://git.hcoop.net/bpt/guile.git/blobdiff_plain/99a0ee662050ad31e74acb3390d6901e3c916f57..37ae02ffa0d788f59c096cec7a3ac9744d87cf16:/libguile/random.c diff --git a/libguile/random.c b/libguile/random.c index a1191898c..1ee0459de 100644 --- a/libguile/random.c +++ b/libguile/random.c @@ -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 */ +/* Original Author: Mikael Djurfeldt */ #ifdef HAVE_CONFIG_H # include @@ -29,6 +31,9 @@ #include #include #include +#include +#include + #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 @@ -413,10 +451,8 @@ 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}.\n" - "\n" - "@var{datum} must be an external state representation obtained\n" - "from @code{random-state->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)); @@ -430,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 @@ -466,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) @@ -494,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) @@ -544,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" @@ -588,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) @@ -622,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 () {