-/* Copyright (C) 1999,2000,2001, 2003 Free Software Foundation, Inc.
+/* Copyright (C) 1999,2000,2001, 2003, 2005, 2006 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
*
* 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
/* Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */
-#if HAVE_CONFIG_H
+#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/feature.h"
#include "libguile/strings.h"
#include "libguile/unif.h"
+#include "libguile/srfi-4.h"
#include "libguile/vectors.h"
#include "libguile/validate.h"
#define M_PI 3.14159265359
#endif
-#ifdef SCM_HAVE_T_INT64
+#if SCM_HAVE_T_UINT64
unsigned long
scm_i_uniform32 (scm_t_i_rstate *state)
{
- scm_t_int64 x = (scm_t_int64) A * state->w + state->c;
- scm_t_int32 w = x & 0xffffffffUL;
+ 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;
unsigned long
scm_i_uniform32 (scm_t_i_rstate *state)
{
- scm_t_int32 x1 = L (A) * L (state->w);
- scm_t_int32 x2 = L (A) * H (state->w);
- scm_t_int32 x3 = H (A) * L (state->w);
- scm_t_int32 w = L (x1) + L (state->c);
- scm_t_int32 m = H (x1) + L (x2) + L (x3) + H (state->c) + H (w);
- scm_t_int32 x4 = H (A) * H (state->w);
+ 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, char *seed, int n)
+scm_i_init_rstate (scm_t_i_rstate *state, const char *seed, int n)
{
- scm_t_int32 w = 0L;
- scm_t_int32 c = 0L;
+ scm_t_uint32 w = 0L;
+ scm_t_uint32 c = 0L;
int i, m;
for (i = 0; i < n; ++i)
{
else
c += seed[i] << (8 * (m - 4));
}
- if ((w == 0 && c == 0) || (w == 0xffffffffUL && c == A - 1))
+ 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);
- if (new_state == 0)
- scm_memory_error ("rstate");
+ scm_t_rstate *new_state;
+
+ new_state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size,
+ "random-state");
return memcpy (new_state, state, scm_the_rng.rstate_size);
}
*/
scm_t_rstate *
-scm_c_make_rstate (char *seed, int n)
+scm_c_make_rstate (const char *seed, int n)
{
- scm_t_rstate *state = scm_malloc (scm_the_rng.rstate_size);
- if (state == 0)
- scm_memory_error ("rstate");
+ scm_t_rstate *state;
+
+ state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size,
+ "random-state");
state->reserved0 = 0;
scm_the_rng.init_rstate (state, seed, n);
return state;
scm_c_default_rstate ()
#define FUNC_NAME "scm_c_default_rstate"
{
- SCM state = SCM_CDR (scm_var_random_state);
+ 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);
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) * 8);
+ 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) * 8);
+ 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 */
(unsigned long *) scm_gc_calloc (num_chunks * sizeof (unsigned long),
"random bignum chunks");
- /* FIXME: what about chance that bignums end up with zeroes at the
- front? -- do we need to normalize? */
-
do
{
unsigned long *current_chunk = random_chunks + (num_chunks - 1);
/* 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) * 8) - end_bits;
+ 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;
scm_gc_free (random_chunks,
num_chunks * sizeof (unsigned long),
"random bignum chunks");
- return result;
+ return scm_i_normbig (result);
}
/*
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_makfrom0str ("URL:http://stat.fsu.edu/~geo/diehard.html")));
+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),
if (SCM_UNBNDP (state))
state = SCM_VARIABLE_REF (scm_var_random_state);
SCM_VALIDATE_RSTATE (2, state);
- if (SCM_INUMP (n))
+ if (SCM_I_INUMP (n))
{
- unsigned long m = SCM_INUM (n);
+ unsigned long m = SCM_I_INUM (n);
SCM_ASSERT_RANGE (1, n, m > 0);
- return SCM_MAKINUM (scm_c_random (SCM_RSTATE (state), m));
+ return scm_from_ulong (scm_c_random (SCM_RSTATE (state), m));
}
SCM_VALIDATE_NIM (1, n);
if (SCM_REALP (n))
- return scm_make_real (SCM_REAL_VALUE (n)
- * scm_c_uniform01 (SCM_RSTATE (state)));
+ return scm_from_double (SCM_REAL_VALUE (n)
+ * scm_c_uniform01 (SCM_RSTATE (state)));
- SCM_VALIDATE_BIGINT (1, n);
+ if (!SCM_BIGP (n))
+ SCM_WRONG_TYPE_ARG (1, n);
return scm_c_random_bignum (SCM_RSTATE (state), n);
}
#undef FUNC_NAME
"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);
- return make_rstate (scm_c_make_rstate (SCM_STRING_CHARS (seed),
- SCM_STRING_LENGTH (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
if (SCM_UNBNDP (state))
state = SCM_VARIABLE_REF (scm_var_random_state);
SCM_VALIDATE_RSTATE (1, state);
- return scm_make_real (scm_c_uniform01 (SCM_RSTATE (state)));
+ return scm_from_double (scm_c_uniform01 (SCM_RSTATE (state)));
}
#undef FUNC_NAME
if (SCM_UNBNDP (state))
state = SCM_VARIABLE_REF (scm_var_random_state);
SCM_VALIDATE_RSTATE (1, state);
- return scm_make_real (scm_c_normal01 (SCM_RSTATE (state)));
+ return scm_from_double (scm_c_normal01 (SCM_RSTATE (state)));
}
#undef FUNC_NAME
-#if SCM_HAVE_ARRAYS
-
static void
-vector_scale (SCM v, double c)
+vector_scale_x (SCM v, double c)
{
- int n = SCM_INUM (scm_uniform_vector_length (v));
- if (SCM_VECTORP (v))
- while (--n >= 0)
- SCM_REAL_VALUE (SCM_VELTS (v)[n]) *= 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
- while (--n >= 0)
- ((double *) SCM_VELTS (v))[n] *= c;
+ {
+ /* 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;
- int n = SCM_INUM (scm_uniform_vector_length (v));
- if (SCM_VECTORP (v))
- while (--n >= 0)
- {
- x = SCM_REAL_VALUE (SCM_VELTS (v)[n]);
- sum += x * x;
- }
+ 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
- while (--n >= 0)
- {
- x = ((double *) SCM_VELTS (v))[n];
- sum += x * x;
- }
+ {
+ /* 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;
}
*/
SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0,
(SCM v, SCM state),
- "Fills vect with inexact real random numbers\n"
- "the sum of whose squares is less than 1.0.\n"
- "Thinking of vect as coordinates in space of\n"
- "dimension n = (vector-length vect), the coordinates\n"
- "are uniformly distributed within the unit n-sphere.\n"
- "The sum of the squares of the numbers is returned.")
+ "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
{
- SCM_VALIDATE_VECTOR_OR_DVECTOR (1, v);
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 (v,
- pow (scm_c_uniform01 (SCM_RSTATE (state)),
- 1.0 / SCM_INUM (scm_uniform_vector_length (v)))
- / sqrt (vector_sum_squares (v)));
+ 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
"unit n-sphere.")
#define FUNC_NAME s_scm_random_hollow_sphere_x
{
- SCM_VALIDATE_VECTOR_OR_DVECTOR (1, v);
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 (v, 1 / sqrt (vector_sum_squares (v)));
+ vector_scale_x (v, 1 / sqrt (vector_sum_squares (v)));
return SCM_UNSPECIFIED;
}
#undef FUNC_NAME
"(i.e., with mean 0 and variance 1).")
#define FUNC_NAME s_scm_random_normal_vector_x
{
- int n;
- SCM_VALIDATE_VECTOR_OR_DVECTOR (1, v);
+ 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);
- n = SCM_INUM (scm_uniform_vector_length (v));
- if (SCM_VECTORP (v))
- while (--n >= 0)
- SCM_VECTOR_SET (v, n, scm_make_real (scm_c_normal01 (SCM_RSTATE (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
- while (--n >= 0)
- ((double *) SCM_VELTS (v))[n] = scm_c_normal01 (SCM_RSTATE (state));
+ {
+ /* 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
-#endif /* SCM_HAVE_ARRAYS */
-
SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0,
(SCM state),
"Return an inexact real in an exponential distribution with mean\n"
if (SCM_UNBNDP (state))
state = SCM_VARIABLE_REF (scm_var_random_state);
SCM_VALIDATE_RSTATE (1, state);
- return scm_make_real (scm_c_exp1 (SCM_RSTATE (state)));
+ return scm_from_double (scm_c_exp1 (SCM_RSTATE (state)));
}
#undef FUNC_NAME
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)
HCoop / bpt / guile.git / blobdiff