diff --git a/libguile/random.c b/libguile/random.c index a8ad07567..6dae7f3ed 100644 --- a/libguile/random.c +++ b/libguile/random.c @@ -1,5 +1,5 @@ /* Copyright (C) 1999, 2000, 2001, 2003, 2005, 2006, 2009, 2010, - * 2012, 2013, 2014 Free Software Foundation, Inc. + * 2012, 2013, 2014, 2017 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 @@ -501,63 +501,73 @@ SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0, static void vector_scale_x (SCM v, double c) { - size_t n; - if (scm_is_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; + scm_t_array_handle handle; + scm_t_array_dim const * dims; + ssize_t i, inc, ubnd; - 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) - *elts *= c; - - scm_array_handle_release (&handle); + 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; - size_t n; - if (scm_is_vector (v)) + 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)) { - n = SCM_SIMPLE_VECTOR_LENGTH (v); - while (n-- > 0) - { - x = SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n)); - sum += x * x; - } + 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; + } } - 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; + 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 @@ -606,16 +616,16 @@ SCM_DEFINE (scm_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0, #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), "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; + scm_t_array_dim const * dims; + ssize_t i; if (SCM_UNBNDP (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_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) { 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))); + 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 = dim->lbnd; i <= dim->ubnd; i++, elts += dim->inc) - *elts = scm_c_normal01 (SCM_RSTATE (state)); + 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_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" @@ -781,13 +790,13 @@ scm_init_random () 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 "libguile/random.x" scm_add_feature ("random"); diff --git a/test-suite/tests/random.test b/test-suite/tests/random.test index ab20b581d..678bd88b7 100644 --- a/test-suite/tests/random.test +++ b/test-suite/tests/random.test @@ -20,7 +20,8 @@ #:use-module ((system base compile) #:select (compile)) #:use-module (test-suite lib) #: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. (define exception:wrong-type-arg @@ -52,4 +53,46 @@ (begin (random:normal-vector! b (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))