mirror/guile
mirror/guile
1
Fork 0
mirror of https://git.savannah.gnu.org/git/guile.git synced 2025-04-30 03:40:34 +02:00
Code Activity

Support general arrays in random:hollow-sphere!

Browse source 
* libguile/random.c (vector_scale_x, vector_sum_squares): Handle general
  rank-1 #t or 'f64 arrays.
* test-suite/tests/random.test: Add tests for random:hollow-sphere!.
This commit is contained in:
Daniel Llorens 2017-03-30 14:29:59 +02:00
parent 46ed57a177
commit 6311ea3a38
2 changed files with 132 additions and 78 deletions
Download patch file Download diff file Expand all files Collapse all files

161
libguile/random.c
View file

@ -140,7 +140,7 @@ scm_i_rstate_from_datum (scm_t_rstate *state, SCM value)
scm_t_i_rstate *istate = (scm_t_i_rstate*) state; scm_t_i_rstate *istate = (scm_t_i_rstate*) state;
scm_t_uint32 w, c; scm_t_uint32 w, c;
long length; long length;
SCM_VALIDATE_LIST_COPYLEN (SCM_ARG1, value, length); SCM_VALIDATE_LIST_COPYLEN (SCM_ARG1, value, length);
SCM_ASSERT (length == 3, value, SCM_ARG1, FUNC_NAME); SCM_ASSERT (length == 3, value, SCM_ARG1, FUNC_NAME);
SCM_ASSERT (scm_is_eq (SCM_CAR (value), scm_i_rstate_tag), SCM_ASSERT (scm_is_eq (SCM_CAR (value), scm_i_rstate_tag),
@ -227,13 +227,13 @@ scm_c_normal01 (scm_t_rstate *state)
else else
{ {
double r, a, n; double r, a, n;
r = sqrt (-2.0 * log (scm_c_uniform01 (state))); r = sqrt (-2.0 * log (scm_c_uniform01 (state)));
a = 2.0 * M_PI * scm_c_uniform01 (state); a = 2.0 * M_PI * scm_c_uniform01 (state);
n = r * sin (a); n = r * sin (a);
state->normal_next = r * cos (a); state->normal_next = r * cos (a);
return n; return n;
} }
} }
@ -274,7 +274,7 @@ scm_c_random64 (scm_t_rstate *state, scm_t_uint64 m)
if (m <= SCM_T_UINT32_MAX) if (m <= SCM_T_UINT32_MAX)
return scm_c_random (state, (scm_t_uint32) m); return scm_c_random (state, (scm_t_uint32) m);
mask = scm_i_mask32 (m >> 32); mask = scm_i_mask32 (m >> 32);
while ((r = ((scm_t_uint64) (state->rng->random_bits (state) & mask) << 32) while ((r = ((scm_t_uint64) (state->rng->random_bits (state) & mask) << 32)
| state->rng->random_bits (state)) >= m) | state->rng->random_bits (state)) >= m)
@ -322,7 +322,7 @@ scm_c_random_bignum (scm_t_rstate *state, SCM m)
scm_t_uint32 chunks_left = num_chunks; scm_t_uint32 chunks_left = num_chunks;
mpz_set_ui (SCM_I_BIG_MPZ (result), 0); mpz_set_ui (SCM_I_BIG_MPZ (result), 0);
if (end_bits) if (end_bits)
{ {
/* generate a mask with ones in the end_bits position, i.e. if /* generate a mask with ones in the end_bits position, i.e. if
@ -334,7 +334,7 @@ scm_c_random_bignum (scm_t_rstate *state, SCM m)
*current_chunk-- = highest_bits; *current_chunk-- = highest_bits;
chunks_left--; chunks_left--;
} }
while (chunks_left) while (chunks_left)
{ {
/* now fill in the remaining scm_t_uint32 sized chunks */ /* now fill in the remaining scm_t_uint32 sized chunks */
@ -360,7 +360,7 @@ scm_c_random_bignum (scm_t_rstate *state, SCM m)
/* /*
* Scheme level representation of random states. * Scheme level representation of random states.
*/ */
scm_t_bits scm_tc16_rstate; scm_t_bits scm_tc16_rstate;
static SCM static SCM
@ -376,7 +376,7 @@ make_rstate (scm_t_rstate *state)
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_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_DEFINE (scm_random, "random", 1, 1, 0,
(SCM n, SCM state), (SCM n, SCM state),
"Return a number in [0, N).\n" "Return a number in [0, N).\n"
"\n" "\n"
@ -420,7 +420,7 @@ SCM_DEFINE (scm_random, "random", 1, 1, 0,
} }
#undef FUNC_NAME #undef FUNC_NAME
SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0, SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0,
(SCM state), (SCM state),
"Return a copy of the random state @var{state}.") "Return a copy of the random state @var{state}.")
#define FUNC_NAME s_scm_copy_random_state #define FUNC_NAME s_scm_copy_random_state
@ -432,7 +432,7 @@ SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0,
} }
#undef FUNC_NAME #undef FUNC_NAME
SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0, SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0,
(SCM seed), (SCM seed),
"Return a new random state using @var{seed}.") "Return a new random state using @var{seed}.")
#define FUNC_NAME s_scm_seed_to_random_state #define FUNC_NAME s_scm_seed_to_random_state
@ -445,11 +445,11 @@ SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0,
scm_i_string_length (seed))); scm_i_string_length (seed)));
scm_remember_upto_here_1 (seed); scm_remember_upto_here_1 (seed);
return res; return res;
} }
#undef FUNC_NAME #undef FUNC_NAME
SCM_DEFINE (scm_datum_to_random_state, "datum->random-state", 1, 0, 0, SCM_DEFINE (scm_datum_to_random_state, "datum->random-state", 1, 0, 0,
(SCM datum), (SCM datum),
"Return a new random state using @var{datum}, which should have\n" "Return a new random state using @var{datum}, which should have\n"
"been obtained from @code{random-state->datum}.") "been obtained from @code{random-state->datum}.")
@ -459,7 +459,7 @@ SCM_DEFINE (scm_datum_to_random_state, "datum->random-state", 1, 0, 0,
} }
#undef FUNC_NAME #undef FUNC_NAME
SCM_DEFINE (scm_random_state_to_datum, "random-state->datum", 1, 0, 0, SCM_DEFINE (scm_random_state_to_datum, "random-state->datum", 1, 0, 0,
(SCM state), (SCM state),
"Return a datum representation of @var{state} that may be\n" "Return a datum representation of @var{state} that may be\n"
"written out and read back with the Scheme reader.") "written out and read back with the Scheme reader.")
@ -470,7 +470,7 @@ SCM_DEFINE (scm_random_state_to_datum, "random-state->datum", 1, 0, 0,
} }
#undef FUNC_NAME #undef FUNC_NAME
SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0, SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0,
(SCM state), (SCM state),
"Return a uniformly distributed inexact real random number in\n" "Return a uniformly distributed inexact real random number in\n"
"[0,1).") "[0,1).")
@ -483,7 +483,7 @@ SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0,
} }
#undef FUNC_NAME #undef FUNC_NAME
SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0, SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0,
(SCM state), (SCM state),
"Return an inexact real in a normal distribution. The\n" "Return an inexact real in a normal distribution. The\n"
"distribution used has mean 0 and standard deviation 1. For a\n" "distribution used has mean 0 and standard deviation 1. For a\n"
@ -501,63 +501,73 @@ SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0,
static void static void
vector_scale_x (SCM v, double c) vector_scale_x (SCM v, double c)
{ {
size_t n; scm_t_array_handle handle;
if (scm_is_vector (v)) scm_t_array_dim const * dims;
{ ssize_t i, inc, ubnd;
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); 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;
for (i = 0; i < len; i++, elts += inc) if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_F64)
*elts *= c; {
double *elts = (double *)(handle.writable_elements) + handle.base;
scm_array_handle_release (&handle); 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 static double
vector_sum_squares (SCM v) vector_sum_squares (SCM v)
{ {
double x, sum = 0.0; double x, sum = 0.0;
size_t n; scm_t_array_handle handle;
if (scm_is_vector (v)) 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))
{ {
n = SCM_SIMPLE_VECTOR_LENGTH (v); ubnd = dims[0].ubnd;
while (n-- > 0) inc = dims[0].inc;
{ if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_F64)
x = SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n)); {
sum += x * x; 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;
}
} }
else scm_array_handle_release (&handle);
{ scm_misc_error (NULL, "must be an array of type #t or 'f64", scm_list_1 (v));
/* 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 /* For the uniform distribution on the solid sphere, note that in
@ -565,7 +575,7 @@ vector_sum_squares (SCM v)
* distribution r^n; i.e., u=r^n is uniform [0,1], so r can be * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be
* generated as r=u^(1/n). * generated as r=u^(1/n).
*/ */
SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0, SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0,
(SCM v, SCM state), (SCM v, SCM state),
"Fills @var{vect} with inexact real random numbers the sum of\n" "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" "whose squares is less than 1.0. Thinking of @var{vect} as\n"
@ -606,16 +616,16 @@ SCM_DEFINE (scm_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0,
#undef FUNC_NAME #undef FUNC_NAME
SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0, SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0,
(SCM v, SCM state), (SCM v, SCM state),
"Fills vect with inexact real random numbers that are\n" "Fills vect with inexact real random numbers that are\n"
"independent and standard normally distributed\n" "independent and standard normally distributed\n"
"(i.e., with mean 0 and variance 1).") "(i.e., with mean 0 and variance 1).")
#define FUNC_NAME s_scm_random_normal_vector_x #define FUNC_NAME s_scm_random_normal_vector_x
{ {
long i;
scm_t_array_handle handle; scm_t_array_handle handle;
scm_t_array_dim *dim; scm_t_array_dim const * dims;
ssize_t i;
if (SCM_UNBNDP (state)) if (SCM_UNBNDP (state))
state = SCM_VARIABLE_REF (scm_var_random_state); state = SCM_VARIABLE_REF (scm_var_random_state);
@ -627,30 +637,29 @@ SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0,
scm_array_handle_release (&handle); scm_array_handle_release (&handle);
scm_wrong_type_arg_msg (NULL, 0, v, "rank 1 array"); scm_wrong_type_arg_msg (NULL, 0, v, "rank 1 array");
} }
dim = scm_array_handle_dims (&handle); dims = scm_array_handle_dims (&handle);
if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_SCM) if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_SCM)
{ {
SCM *elts = scm_array_handle_writable_elements (&handle); SCM *elts = scm_array_handle_writable_elements (&handle);
for (i = dim->lbnd; i <= dim->ubnd; i++, elts += dim->inc) for (i = dims->lbnd; i <= dims->ubnd; i++, elts += dims->inc)
*elts = scm_from_double (scm_c_normal01 (SCM_RSTATE (state))); *elts = scm_from_double (scm_c_normal01 (SCM_RSTATE (state)));
} }
else else
{ {
/* must be a f64vector. */ /* must be a f64vector. */
double *elts = scm_array_handle_f64_writable_elements (&handle); double *elts = scm_array_handle_f64_writable_elements (&handle);
for (i = dim->lbnd; i <= dim->ubnd; i++, elts += dim->inc) for (i = dims->lbnd; i <= dims->ubnd; i++, elts += dims->inc)
*elts = scm_c_normal01 (SCM_RSTATE (state)); *elts = scm_c_normal01 (SCM_RSTATE (state));
} }
scm_array_handle_release (&handle); scm_array_handle_release (&handle);
return SCM_UNSPECIFIED; return SCM_UNSPECIFIED;
} }
#undef FUNC_NAME #undef FUNC_NAME
SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0, SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0,
(SCM state), (SCM state),
"Return an inexact real in an exponential distribution with mean\n" "Return an inexact real in an exponential distribution with mean\n"
"1. For an exponential distribution with mean u use (* u\n" "1. For an exponential distribution with mean u use (* u\n"
@ -781,13 +790,13 @@ scm_init_random ()
scm_i_rstate_to_datum scm_i_rstate_to_datum
}; };
scm_the_rng = rng; scm_the_rng = rng;
scm_tc16_rstate = scm_make_smob_type ("random-state", 0); scm_tc16_rstate = scm_make_smob_type ("random-state", 0);
for (m = 1; m <= 0x100; m <<= 1) for (m = 1; m <= 0x100; m <<= 1)
for (i = m >> 1; i < m; ++i) for (i = m >> 1; i < m; ++i)
scm_masktab[i] = m - 1; scm_masktab[i] = m - 1;
#include "libguile/random.x" #include "libguile/random.x"
scm_add_feature ("random"); scm_add_feature ("random");

49
test-suite/tests/random.test
View file

@ -20,7 +20,8 @@
#:use-module ((system base compile) #:select (compile)) #:use-module ((system base compile) #:select (compile))
#:use-module (test-suite lib) #:use-module (test-suite lib)
#:use-module (srfi srfi-4) #:use-module (srfi srfi-4)
#:use-module (srfi srfi-4 gnu)) #:use-module (srfi srfi-4 gnu)
#:use-module ((ice-9 control) #:select (let/ec)))
; see strings.test, arrays.test. ; see strings.test, arrays.test.
(define exception:wrong-type-arg (define exception:wrong-type-arg
@ -52,4 +53,48 @@
(begin (begin
(random:normal-vector! b (random-state-from-platform)) (random:normal-vector! b (random-state-from-platform))
(random:normal-vector! c (random-state-from-platform)) (random:normal-vector! c (random-state-from-platform))
(and (not (equal? a b)) (not (equal? a c))))))) (and (not (equal? a b)) (not (equal? a c))))))
(pass-if "empty argument"
(random:normal-vector! (vector) (random-state-from-platform))
(random:normal-vector! (f64vector) (random-state-from-platform))
#t))
;;;
;;; random:hollow-sphere!
;;;
(with-test-prefix "random:hollow-sphere!"
(define (sqr a)
(* a a))
(define (norm a)
(sqrt (+ (sqr (array-ref a 0)) (sqr (array-ref a 1)) (sqr (array-ref a 2)))))
(define double-eps 1e-15)
(pass-if "non uniform"
(let ((a (transpose-array (make-array 0. 3 10) 1 0)))
(let/ec exit
(array-slice-for-each 1
(lambda (a)
(random:hollow-sphere! a)
(if (> (magnitude (- 1 (norm a))) double-eps) (exit #f)))
a)
#t)))
(pass-if "uniform (f64)"
(let ((a (transpose-array (make-array 0. 3 10) 1 0)))
(let/ec exit
(array-slice-for-each 1
(lambda (a)
(random:hollow-sphere! a)
(if (> (magnitude (- 1 (norm a))) double-eps) (exit #f)))
a)
#t)))
(pass-if "empty argument"
(random:hollow-sphere! (vector))
(random:hollow-sphere! (f64vector))
#t))