| 1 | /* Copyright (C) 1999, 2000, 2001, 2003, 2005, 2006, 2009, 2010, |
| 2 | * 2012, 2013, 2014 Free Software Foundation, Inc. |
| 3 | * |
| 4 | * This library is free software; you can redistribute it and/or |
| 5 | * modify it under the terms of the GNU Lesser General Public License |
| 6 | * as published by the Free Software Foundation; either version 3 of |
| 7 | * the License, or (at your option) any later version. |
| 8 | * |
| 9 | * This library is distributed in the hope that it will be useful, but |
| 10 | * WITHOUT ANY WARRANTY; without even the implied warranty of |
| 11 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| 12 | * Lesser General Public License for more details. |
| 13 | * |
| 14 | * You should have received a copy of the GNU Lesser General Public |
| 15 | * License along with this library; if not, write to the Free Software |
| 16 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
| 17 | * 02110-1301 USA |
| 18 | */ |
| 19 | |
| 20 | |
| 21 | |
| 22 | /* Original Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */ |
| 23 | |
| 24 | #ifdef HAVE_CONFIG_H |
| 25 | # include <config.h> |
| 26 | #endif |
| 27 | |
| 28 | #include "libguile/_scm.h" |
| 29 | |
| 30 | #include <gmp.h> |
| 31 | #include <stdio.h> |
| 32 | #include <math.h> |
| 33 | #include <string.h> |
| 34 | #include <sys/types.h> |
| 35 | #include <unistd.h> |
| 36 | |
| 37 | #include "libguile/smob.h" |
| 38 | #include "libguile/numbers.h" |
| 39 | #include "libguile/feature.h" |
| 40 | #include "libguile/strings.h" |
| 41 | #include "libguile/arrays.h" |
| 42 | #include "libguile/srfi-4.h" |
| 43 | #include "libguile/vectors.h" |
| 44 | #include "libguile/generalized-vectors.h" |
| 45 | |
| 46 | #include "libguile/validate.h" |
| 47 | #include "libguile/random.h" |
| 48 | |
| 49 | \f |
| 50 | /* |
| 51 | * A plugin interface for RNGs |
| 52 | * |
| 53 | * Using this interface, it is possible for the application to tell |
| 54 | * libguile to use a different RNG. This is desirable if it is |
| 55 | * necessary to use the same RNG everywhere in the application in |
| 56 | * order to prevent interference, if the application uses RNG |
| 57 | * hardware, or if the application has special demands on the RNG. |
| 58 | * |
| 59 | * Look in random.h and how the default generator is "plugged in" in |
| 60 | * scm_init_random(). |
| 61 | */ |
| 62 | |
| 63 | scm_t_rng scm_the_rng; |
| 64 | |
| 65 | \f |
| 66 | /* |
| 67 | * The prepackaged RNG |
| 68 | * |
| 69 | * This is the MWC (Multiply With Carry) random number generator |
| 70 | * described by George Marsaglia at the Department of Statistics and |
| 71 | * Supercomputer Computations Research Institute, The Florida State |
| 72 | * University (http://stat.fsu.edu/~geo). |
| 73 | * |
| 74 | * It uses 64 bits, has a period of 4578426017172946943 (4.6e18), and |
| 75 | * passes all tests in the DIEHARD test suite |
| 76 | * (http://stat.fsu.edu/~geo/diehard.html) |
| 77 | */ |
| 78 | |
| 79 | typedef struct scm_t_i_rstate { |
| 80 | scm_t_rstate rstate; |
| 81 | scm_t_uint32 w; |
| 82 | scm_t_uint32 c; |
| 83 | } scm_t_i_rstate; |
| 84 | |
| 85 | |
| 86 | #define A 2131995753UL |
| 87 | |
| 88 | #ifndef M_PI |
| 89 | #define M_PI 3.14159265359 |
| 90 | #endif |
| 91 | |
| 92 | static scm_t_uint32 |
| 93 | scm_i_uniform32 (scm_t_rstate *state) |
| 94 | { |
| 95 | scm_t_i_rstate *istate = (scm_t_i_rstate*) state; |
| 96 | scm_t_uint64 x = (scm_t_uint64) A * istate->w + istate->c; |
| 97 | scm_t_uint32 w = x & 0xffffffffUL; |
| 98 | istate->w = w; |
| 99 | istate->c = x >> 32L; |
| 100 | return w; |
| 101 | } |
| 102 | |
| 103 | static void |
| 104 | scm_i_init_rstate (scm_t_rstate *state, const char *seed, int n) |
| 105 | { |
| 106 | scm_t_i_rstate *istate = (scm_t_i_rstate*) state; |
| 107 | scm_t_uint32 w = 0L; |
| 108 | scm_t_uint32 c = 0L; |
| 109 | int i, m; |
| 110 | for (i = 0; i < n; ++i) |
| 111 | { |
| 112 | m = i % 8; |
| 113 | if (m < 4) |
| 114 | w += seed[i] << (8 * m); |
| 115 | else |
| 116 | c += seed[i] << (8 * (m - 4)); |
| 117 | } |
| 118 | if ((w == 0 && c == 0) || (w == -1 && c == A - 1)) |
| 119 | ++c; |
| 120 | istate->w = w; |
| 121 | istate->c = c; |
| 122 | } |
| 123 | |
| 124 | static scm_t_rstate * |
| 125 | scm_i_copy_rstate (scm_t_rstate *state) |
| 126 | { |
| 127 | scm_t_rstate *new_state; |
| 128 | |
| 129 | new_state = scm_gc_malloc_pointerless (state->rng->rstate_size, |
| 130 | "random-state"); |
| 131 | return memcpy (new_state, state, state->rng->rstate_size); |
| 132 | } |
| 133 | |
| 134 | SCM_SYMBOL(scm_i_rstate_tag, "multiply-with-carry"); |
| 135 | |
| 136 | static void |
| 137 | scm_i_rstate_from_datum (scm_t_rstate *state, SCM value) |
| 138 | #define FUNC_NAME "scm_i_rstate_from_datum" |
| 139 | { |
| 140 | scm_t_i_rstate *istate = (scm_t_i_rstate*) state; |
| 141 | scm_t_uint32 w, c; |
| 142 | long length; |
| 143 | |
| 144 | SCM_VALIDATE_LIST_COPYLEN (SCM_ARG1, value, length); |
| 145 | SCM_ASSERT (length == 3, value, SCM_ARG1, FUNC_NAME); |
| 146 | SCM_ASSERT (scm_is_eq (SCM_CAR (value), scm_i_rstate_tag), |
| 147 | value, SCM_ARG1, FUNC_NAME); |
| 148 | SCM_VALIDATE_UINT_COPY (SCM_ARG1, SCM_CADR (value), w); |
| 149 | SCM_VALIDATE_UINT_COPY (SCM_ARG1, SCM_CADDR (value), c); |
| 150 | |
| 151 | istate->w = w; |
| 152 | istate->c = c; |
| 153 | } |
| 154 | #undef FUNC_NAME |
| 155 | |
| 156 | static SCM |
| 157 | scm_i_rstate_to_datum (scm_t_rstate *state) |
| 158 | { |
| 159 | scm_t_i_rstate *istate = (scm_t_i_rstate*) state; |
| 160 | return scm_list_3 (scm_i_rstate_tag, |
| 161 | scm_from_uint32 (istate->w), |
| 162 | scm_from_uint32 (istate->c)); |
| 163 | } |
| 164 | |
| 165 | \f |
| 166 | /* |
| 167 | * Random number library functions |
| 168 | */ |
| 169 | |
| 170 | scm_t_rstate * |
| 171 | scm_c_make_rstate (const char *seed, int n) |
| 172 | { |
| 173 | scm_t_rstate *state; |
| 174 | |
| 175 | state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size, |
| 176 | "random-state"); |
| 177 | state->rng = &scm_the_rng; |
| 178 | state->normal_next = 0.0; |
| 179 | state->rng->init_rstate (state, seed, n); |
| 180 | return state; |
| 181 | } |
| 182 | |
| 183 | scm_t_rstate * |
| 184 | scm_c_rstate_from_datum (SCM datum) |
| 185 | { |
| 186 | scm_t_rstate *state; |
| 187 | |
| 188 | state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size, |
| 189 | "random-state"); |
| 190 | state->rng = &scm_the_rng; |
| 191 | state->normal_next = 0.0; |
| 192 | state->rng->from_datum (state, datum); |
| 193 | return state; |
| 194 | } |
| 195 | |
| 196 | scm_t_rstate * |
| 197 | scm_c_default_rstate () |
| 198 | #define FUNC_NAME "scm_c_default_rstate" |
| 199 | { |
| 200 | SCM state = SCM_VARIABLE_REF (scm_var_random_state); |
| 201 | if (!SCM_RSTATEP (state)) |
| 202 | SCM_MISC_ERROR ("*random-state* contains bogus random state", SCM_EOL); |
| 203 | return SCM_RSTATE (state); |
| 204 | } |
| 205 | #undef FUNC_NAME |
| 206 | |
| 207 | |
| 208 | double |
| 209 | scm_c_uniform01 (scm_t_rstate *state) |
| 210 | { |
| 211 | double x = (double) state->rng->random_bits (state) / (double) 0xffffffffUL; |
| 212 | return ((x + (double) state->rng->random_bits (state)) |
| 213 | / (double) 0xffffffffUL); |
| 214 | } |
| 215 | |
| 216 | double |
| 217 | scm_c_normal01 (scm_t_rstate *state) |
| 218 | { |
| 219 | if (state->normal_next != 0.0) |
| 220 | { |
| 221 | double ret = state->normal_next; |
| 222 | |
| 223 | state->normal_next = 0.0; |
| 224 | |
| 225 | return ret; |
| 226 | } |
| 227 | else |
| 228 | { |
| 229 | double r, a, n; |
| 230 | |
| 231 | r = sqrt (-2.0 * log (scm_c_uniform01 (state))); |
| 232 | a = 2.0 * M_PI * scm_c_uniform01 (state); |
| 233 | |
| 234 | n = r * sin (a); |
| 235 | state->normal_next = r * cos (a); |
| 236 | |
| 237 | return n; |
| 238 | } |
| 239 | } |
| 240 | |
| 241 | double |
| 242 | scm_c_exp1 (scm_t_rstate *state) |
| 243 | { |
| 244 | return - log (scm_c_uniform01 (state)); |
| 245 | } |
| 246 | |
| 247 | unsigned char scm_masktab[256]; |
| 248 | |
| 249 | static inline scm_t_uint32 |
| 250 | scm_i_mask32 (scm_t_uint32 m) |
| 251 | { |
| 252 | return (m < 0x100 |
| 253 | ? scm_masktab[m] |
| 254 | : (m < 0x10000 |
| 255 | ? scm_masktab[m >> 8] << 8 | 0xff |
| 256 | : (m < 0x1000000 |
| 257 | ? scm_masktab[m >> 16] << 16 | 0xffff |
| 258 | : ((scm_t_uint32) scm_masktab[m >> 24]) << 24 | 0xffffff))); |
| 259 | } |
| 260 | |
| 261 | scm_t_uint32 |
| 262 | scm_c_random (scm_t_rstate *state, scm_t_uint32 m) |
| 263 | { |
| 264 | scm_t_uint32 r, mask = scm_i_mask32 (m); |
| 265 | while ((r = state->rng->random_bits (state) & mask) >= m); |
| 266 | return r; |
| 267 | } |
| 268 | |
| 269 | scm_t_uint64 |
| 270 | scm_c_random64 (scm_t_rstate *state, scm_t_uint64 m) |
| 271 | { |
| 272 | scm_t_uint64 r; |
| 273 | scm_t_uint32 mask; |
| 274 | |
| 275 | if (m <= SCM_T_UINT32_MAX) |
| 276 | return scm_c_random (state, (scm_t_uint32) m); |
| 277 | |
| 278 | mask = scm_i_mask32 (m >> 32); |
| 279 | while ((r = ((scm_t_uint64) (state->rng->random_bits (state) & mask) << 32) |
| 280 | | state->rng->random_bits (state)) >= m) |
| 281 | ; |
| 282 | return r; |
| 283 | } |
| 284 | |
| 285 | /* |
| 286 | SCM scm_c_random_bignum (scm_t_rstate *state, SCM m) |
| 287 | |
| 288 | Takes a random state (source of random bits) and a bignum m. |
| 289 | Returns a bignum b, 0 <= b < m. |
| 290 | |
| 291 | It does this by allocating a bignum b with as many base 65536 digits |
| 292 | as m, filling b with random bits (in 32 bit chunks) up to the most |
| 293 | significant 1 in m, and, finally checking if the resultant b is too |
| 294 | large (>= m). If too large, we simply repeat the process again. (It |
| 295 | is important to throw away all generated random bits if b >= m, |
| 296 | otherwise we'll end up with a distorted distribution.) |
| 297 | |
| 298 | */ |
| 299 | |
| 300 | SCM |
| 301 | scm_c_random_bignum (scm_t_rstate *state, SCM m) |
| 302 | { |
| 303 | SCM result = scm_i_mkbig (); |
| 304 | const size_t m_bits = mpz_sizeinbase (SCM_I_BIG_MPZ (m), 2); |
| 305 | /* how many bits would only partially fill the last scm_t_uint32? */ |
| 306 | const size_t end_bits = m_bits % (sizeof (scm_t_uint32) * SCM_CHAR_BIT); |
| 307 | scm_t_uint32 *random_chunks = NULL; |
| 308 | const scm_t_uint32 num_full_chunks = |
| 309 | m_bits / (sizeof (scm_t_uint32) * SCM_CHAR_BIT); |
| 310 | const scm_t_uint32 num_chunks = num_full_chunks + ((end_bits) ? 1 : 0); |
| 311 | |
| 312 | /* we know the result will be this big */ |
| 313 | mpz_realloc2 (SCM_I_BIG_MPZ (result), m_bits); |
| 314 | |
| 315 | random_chunks = |
| 316 | (scm_t_uint32 *) scm_gc_calloc (num_chunks * sizeof (scm_t_uint32), |
| 317 | "random bignum chunks"); |
| 318 | |
| 319 | do |
| 320 | { |
| 321 | scm_t_uint32 *current_chunk = random_chunks + (num_chunks - 1); |
| 322 | scm_t_uint32 chunks_left = num_chunks; |
| 323 | |
| 324 | mpz_set_ui (SCM_I_BIG_MPZ (result), 0); |
| 325 | |
| 326 | if (end_bits) |
| 327 | { |
| 328 | /* generate a mask with ones in the end_bits position, i.e. if |
| 329 | end_bits is 3, then we'd have a mask of ...0000000111 */ |
| 330 | const scm_t_uint32 rndbits = state->rng->random_bits (state); |
| 331 | int rshift = (sizeof (scm_t_uint32) * SCM_CHAR_BIT) - end_bits; |
| 332 | scm_t_uint32 mask = ((scm_t_uint32)-1) >> rshift; |
| 333 | scm_t_uint32 highest_bits = rndbits & mask; |
| 334 | *current_chunk-- = highest_bits; |
| 335 | chunks_left--; |
| 336 | } |
| 337 | |
| 338 | while (chunks_left) |
| 339 | { |
| 340 | /* now fill in the remaining scm_t_uint32 sized chunks */ |
| 341 | *current_chunk-- = state->rng->random_bits (state); |
| 342 | chunks_left--; |
| 343 | } |
| 344 | mpz_import (SCM_I_BIG_MPZ (result), |
| 345 | num_chunks, |
| 346 | -1, |
| 347 | sizeof (scm_t_uint32), |
| 348 | 0, |
| 349 | 0, |
| 350 | random_chunks); |
| 351 | /* if result >= m, regenerate it (it is important to regenerate |
| 352 | all bits in order not to get a distorted distribution) */ |
| 353 | } while (mpz_cmp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (m)) >= 0); |
| 354 | scm_gc_free (random_chunks, |
| 355 | num_chunks * sizeof (scm_t_uint32), |
| 356 | "random bignum chunks"); |
| 357 | return scm_i_normbig (result); |
| 358 | } |
| 359 | |
| 360 | /* |
| 361 | * Scheme level representation of random states. |
| 362 | */ |
| 363 | |
| 364 | scm_t_bits scm_tc16_rstate; |
| 365 | |
| 366 | static SCM |
| 367 | make_rstate (scm_t_rstate *state) |
| 368 | { |
| 369 | SCM_RETURN_NEWSMOB (scm_tc16_rstate, state); |
| 370 | } |
| 371 | |
| 372 | |
| 373 | /* |
| 374 | * Scheme level interface. |
| 375 | */ |
| 376 | |
| 377 | 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"))); |
| 378 | |
| 379 | SCM_DEFINE (scm_random, "random", 1, 1, 0, |
| 380 | (SCM n, SCM state), |
| 381 | "Return a number in [0, N).\n" |
| 382 | "\n" |
| 383 | "Accepts a positive integer or real n and returns a\n" |
| 384 | "number of the same type between zero (inclusive) and\n" |
| 385 | "N (exclusive). The values returned have a uniform\n" |
| 386 | "distribution.\n" |
| 387 | "\n" |
| 388 | "The optional argument @var{state} must be of the type produced\n" |
| 389 | "by @code{seed->random-state}. It defaults to the value of the\n" |
| 390 | "variable @var{*random-state*}. This object is used to maintain\n" |
| 391 | "the state of the pseudo-random-number generator and is altered\n" |
| 392 | "as a side effect of the random operation.") |
| 393 | #define FUNC_NAME s_scm_random |
| 394 | { |
| 395 | if (SCM_UNBNDP (state)) |
| 396 | state = SCM_VARIABLE_REF (scm_var_random_state); |
| 397 | SCM_VALIDATE_RSTATE (2, state); |
| 398 | if (SCM_I_INUMP (n)) |
| 399 | { |
| 400 | scm_t_bits m = (scm_t_bits) SCM_I_INUM (n); |
| 401 | SCM_ASSERT_RANGE (1, n, SCM_I_INUM (n) > 0); |
| 402 | #if SCM_SIZEOF_UINTPTR_T <= 4 |
| 403 | return scm_from_uint32 (scm_c_random (SCM_RSTATE (state), |
| 404 | (scm_t_uint32) m)); |
| 405 | #elif SCM_SIZEOF_UINTPTR_T <= 8 |
| 406 | return scm_from_uint64 (scm_c_random64 (SCM_RSTATE (state), |
| 407 | (scm_t_uint64) m)); |
| 408 | #else |
| 409 | #error "Cannot deal with this platform's scm_t_bits size" |
| 410 | #endif |
| 411 | } |
| 412 | SCM_VALIDATE_NIM (1, n); |
| 413 | if (SCM_REALP (n)) |
| 414 | return scm_from_double (SCM_REAL_VALUE (n) |
| 415 | * scm_c_uniform01 (SCM_RSTATE (state))); |
| 416 | |
| 417 | if (!SCM_BIGP (n)) |
| 418 | SCM_WRONG_TYPE_ARG (1, n); |
| 419 | return scm_c_random_bignum (SCM_RSTATE (state), n); |
| 420 | } |
| 421 | #undef FUNC_NAME |
| 422 | |
| 423 | SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0, |
| 424 | (SCM state), |
| 425 | "Return a copy of the random state @var{state}.") |
| 426 | #define FUNC_NAME s_scm_copy_random_state |
| 427 | { |
| 428 | if (SCM_UNBNDP (state)) |
| 429 | state = SCM_VARIABLE_REF (scm_var_random_state); |
| 430 | SCM_VALIDATE_RSTATE (1, state); |
| 431 | return make_rstate (SCM_RSTATE (state)->rng->copy_rstate (SCM_RSTATE (state))); |
| 432 | } |
| 433 | #undef FUNC_NAME |
| 434 | |
| 435 | SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0, |
| 436 | (SCM seed), |
| 437 | "Return a new random state using @var{seed}.") |
| 438 | #define FUNC_NAME s_scm_seed_to_random_state |
| 439 | { |
| 440 | SCM res; |
| 441 | if (SCM_NUMBERP (seed)) |
| 442 | seed = scm_number_to_string (seed, SCM_UNDEFINED); |
| 443 | SCM_VALIDATE_STRING (1, seed); |
| 444 | res = make_rstate (scm_c_make_rstate (scm_i_string_chars (seed), |
| 445 | scm_i_string_length (seed))); |
| 446 | scm_remember_upto_here_1 (seed); |
| 447 | return res; |
| 448 | |
| 449 | } |
| 450 | #undef FUNC_NAME |
| 451 | |
| 452 | SCM_DEFINE (scm_datum_to_random_state, "datum->random-state", 1, 0, 0, |
| 453 | (SCM datum), |
| 454 | "Return a new random state using @var{datum}, which should have\n" |
| 455 | "been obtained from @code{random-state->datum}.") |
| 456 | #define FUNC_NAME s_scm_datum_to_random_state |
| 457 | { |
| 458 | return make_rstate (scm_c_rstate_from_datum (datum)); |
| 459 | } |
| 460 | #undef FUNC_NAME |
| 461 | |
| 462 | SCM_DEFINE (scm_random_state_to_datum, "random-state->datum", 1, 0, 0, |
| 463 | (SCM state), |
| 464 | "Return a datum representation of @var{state} that may be\n" |
| 465 | "written out and read back with the Scheme reader.") |
| 466 | #define FUNC_NAME s_scm_random_state_to_datum |
| 467 | { |
| 468 | SCM_VALIDATE_RSTATE (1, state); |
| 469 | return SCM_RSTATE (state)->rng->to_datum (SCM_RSTATE (state)); |
| 470 | } |
| 471 | #undef FUNC_NAME |
| 472 | |
| 473 | SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0, |
| 474 | (SCM state), |
| 475 | "Return a uniformly distributed inexact real random number in\n" |
| 476 | "[0,1).") |
| 477 | #define FUNC_NAME s_scm_random_uniform |
| 478 | { |
| 479 | if (SCM_UNBNDP (state)) |
| 480 | state = SCM_VARIABLE_REF (scm_var_random_state); |
| 481 | SCM_VALIDATE_RSTATE (1, state); |
| 482 | return scm_from_double (scm_c_uniform01 (SCM_RSTATE (state))); |
| 483 | } |
| 484 | #undef FUNC_NAME |
| 485 | |
| 486 | SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0, |
| 487 | (SCM state), |
| 488 | "Return an inexact real in a normal distribution. The\n" |
| 489 | "distribution used has mean 0 and standard deviation 1. For a\n" |
| 490 | "normal distribution with mean m and standard deviation d use\n" |
| 491 | "@code{(+ m (* d (random:normal)))}.") |
| 492 | #define FUNC_NAME s_scm_random_normal |
| 493 | { |
| 494 | if (SCM_UNBNDP (state)) |
| 495 | state = SCM_VARIABLE_REF (scm_var_random_state); |
| 496 | SCM_VALIDATE_RSTATE (1, state); |
| 497 | return scm_from_double (scm_c_normal01 (SCM_RSTATE (state))); |
| 498 | } |
| 499 | #undef FUNC_NAME |
| 500 | |
| 501 | static void |
| 502 | vector_scale_x (SCM v, double c) |
| 503 | { |
| 504 | size_t n; |
| 505 | if (scm_is_vector (v)) |
| 506 | { |
| 507 | n = SCM_SIMPLE_VECTOR_LENGTH (v); |
| 508 | while (n-- > 0) |
| 509 | SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n)) *= c; |
| 510 | } |
| 511 | else |
| 512 | { |
| 513 | /* must be a f64vector. */ |
| 514 | scm_t_array_handle handle; |
| 515 | size_t i, len; |
| 516 | ssize_t inc; |
| 517 | double *elts; |
| 518 | |
| 519 | elts = scm_f64vector_writable_elements (v, &handle, &len, &inc); |
| 520 | |
| 521 | for (i = 0; i < len; i++, elts += inc) |
| 522 | *elts *= c; |
| 523 | |
| 524 | scm_array_handle_release (&handle); |
| 525 | } |
| 526 | } |
| 527 | |
| 528 | static double |
| 529 | vector_sum_squares (SCM v) |
| 530 | { |
| 531 | double x, sum = 0.0; |
| 532 | size_t n; |
| 533 | if (scm_is_vector (v)) |
| 534 | { |
| 535 | n = SCM_SIMPLE_VECTOR_LENGTH (v); |
| 536 | while (n-- > 0) |
| 537 | { |
| 538 | x = SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n)); |
| 539 | sum += x * x; |
| 540 | } |
| 541 | } |
| 542 | else |
| 543 | { |
| 544 | /* must be a f64vector. */ |
| 545 | scm_t_array_handle handle; |
| 546 | size_t i, len; |
| 547 | ssize_t inc; |
| 548 | const double *elts; |
| 549 | |
| 550 | elts = scm_f64vector_elements (v, &handle, &len, &inc); |
| 551 | |
| 552 | for (i = 0; i < len; i++, elts += inc) |
| 553 | { |
| 554 | x = *elts; |
| 555 | sum += x * x; |
| 556 | } |
| 557 | |
| 558 | scm_array_handle_release (&handle); |
| 559 | } |
| 560 | return sum; |
| 561 | } |
| 562 | |
| 563 | /* For the uniform distribution on the solid sphere, note that in |
| 564 | * this distribution the length r of the vector has cumulative |
| 565 | * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be |
| 566 | * generated as r=u^(1/n). |
| 567 | */ |
| 568 | SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0, |
| 569 | (SCM v, SCM state), |
| 570 | "Fills @var{vect} with inexact real random numbers the sum of\n" |
| 571 | "whose squares is less than 1.0. Thinking of @var{vect} as\n" |
| 572 | "coordinates in space of dimension @var{n} @math{=}\n" |
| 573 | "@code{(vector-length @var{vect})}, the coordinates are\n" |
| 574 | "uniformly distributed within the unit @var{n}-sphere.") |
| 575 | #define FUNC_NAME s_scm_random_solid_sphere_x |
| 576 | { |
| 577 | if (SCM_UNBNDP (state)) |
| 578 | state = SCM_VARIABLE_REF (scm_var_random_state); |
| 579 | SCM_VALIDATE_RSTATE (2, state); |
| 580 | scm_random_normal_vector_x (v, state); |
| 581 | vector_scale_x (v, |
| 582 | pow (scm_c_uniform01 (SCM_RSTATE (state)), |
| 583 | 1.0 / scm_c_array_length (v)) |
| 584 | / sqrt (vector_sum_squares (v))); |
| 585 | return SCM_UNSPECIFIED; |
| 586 | } |
| 587 | #undef FUNC_NAME |
| 588 | |
| 589 | SCM_DEFINE (scm_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0, |
| 590 | (SCM v, SCM state), |
| 591 | "Fills vect with inexact real random numbers\n" |
| 592 | "the sum of whose squares is equal to 1.0.\n" |
| 593 | "Thinking of vect as coordinates in space of\n" |
| 594 | "dimension n = (vector-length vect), the coordinates\n" |
| 595 | "are uniformly distributed over the surface of the\n" |
| 596 | "unit n-sphere.") |
| 597 | #define FUNC_NAME s_scm_random_hollow_sphere_x |
| 598 | { |
| 599 | if (SCM_UNBNDP (state)) |
| 600 | state = SCM_VARIABLE_REF (scm_var_random_state); |
| 601 | SCM_VALIDATE_RSTATE (2, state); |
| 602 | scm_random_normal_vector_x (v, state); |
| 603 | vector_scale_x (v, 1 / sqrt (vector_sum_squares (v))); |
| 604 | return SCM_UNSPECIFIED; |
| 605 | } |
| 606 | #undef FUNC_NAME |
| 607 | |
| 608 | |
| 609 | SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0, |
| 610 | (SCM v, SCM state), |
| 611 | "Fills vect with inexact real random numbers that are\n" |
| 612 | "independent and standard normally distributed\n" |
| 613 | "(i.e., with mean 0 and variance 1).") |
| 614 | #define FUNC_NAME s_scm_random_normal_vector_x |
| 615 | { |
| 616 | long i; |
| 617 | scm_t_array_handle handle; |
| 618 | scm_t_array_dim *dim; |
| 619 | |
| 620 | if (SCM_UNBNDP (state)) |
| 621 | state = SCM_VARIABLE_REF (scm_var_random_state); |
| 622 | SCM_VALIDATE_RSTATE (2, state); |
| 623 | |
| 624 | scm_generalized_vector_get_handle (v, &handle); |
| 625 | dim = scm_array_handle_dims (&handle); |
| 626 | |
| 627 | if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_SCM) |
| 628 | { |
| 629 | SCM *elts = scm_array_handle_writable_elements (&handle); |
| 630 | for (i = dim->lbnd; i <= dim->ubnd; i++, elts += dim->inc) |
| 631 | *elts = scm_from_double (scm_c_normal01 (SCM_RSTATE (state))); |
| 632 | } |
| 633 | else |
| 634 | { |
| 635 | /* must be a f64vector. */ |
| 636 | double *elts = scm_array_handle_f64_writable_elements (&handle); |
| 637 | for (i = dim->lbnd; i <= dim->ubnd; i++, elts += dim->inc) |
| 638 | *elts = scm_c_normal01 (SCM_RSTATE (state)); |
| 639 | } |
| 640 | |
| 641 | scm_array_handle_release (&handle); |
| 642 | |
| 643 | return SCM_UNSPECIFIED; |
| 644 | } |
| 645 | #undef FUNC_NAME |
| 646 | |
| 647 | SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0, |
| 648 | (SCM state), |
| 649 | "Return an inexact real in an exponential distribution with mean\n" |
| 650 | "1. For an exponential distribution with mean u use (* u\n" |
| 651 | "(random:exp)).") |
| 652 | #define FUNC_NAME s_scm_random_exp |
| 653 | { |
| 654 | if (SCM_UNBNDP (state)) |
| 655 | state = SCM_VARIABLE_REF (scm_var_random_state); |
| 656 | SCM_VALIDATE_RSTATE (1, state); |
| 657 | return scm_from_double (scm_c_exp1 (SCM_RSTATE (state))); |
| 658 | } |
| 659 | #undef FUNC_NAME |
| 660 | |
| 661 | /* Return a new random-state seeded from the time, date, process ID, an |
| 662 | address from a freshly allocated heap cell, an address from the local |
| 663 | stack frame, and a high-resolution timer if available. This is only |
| 664 | to be used as a last resort, when no better source of entropy is |
| 665 | available. */ |
| 666 | static SCM |
| 667 | random_state_of_last_resort (void) |
| 668 | { |
| 669 | SCM state; |
| 670 | SCM time_of_day = scm_gettimeofday (); |
| 671 | SCM sources = scm_list_n |
| 672 | (scm_from_unsigned_integer (SCM_UNPACK (time_of_day)), /* heap addr */ |
| 673 | /* Avoid scm_getpid, since it depends on HAVE_POSIX. */ |
| 674 | scm_from_unsigned_integer (getpid ()), /* process ID */ |
| 675 | scm_get_internal_real_time (), /* high-resolution process timer */ |
| 676 | scm_from_unsigned_integer ((scm_t_bits) &time_of_day), /* stack addr */ |
| 677 | scm_car (time_of_day), /* seconds since midnight 1970-01-01 UTC */ |
| 678 | scm_cdr (time_of_day), /* microsecond component of the above clock */ |
| 679 | SCM_UNDEFINED); |
| 680 | |
| 681 | /* Concatenate the sources bitwise to form the seed */ |
| 682 | SCM seed = SCM_INUM0; |
| 683 | while (scm_is_pair (sources)) |
| 684 | { |
| 685 | seed = scm_logxor (seed, scm_ash (scm_car (sources), |
| 686 | scm_integer_length (seed))); |
| 687 | sources = scm_cdr (sources); |
| 688 | } |
| 689 | |
| 690 | /* FIXME The following code belongs in `scm_seed_to_random_state', |
| 691 | and here we should simply do: |
| 692 | |
| 693 | return scm_seed_to_random_state (seed); |
| 694 | |
| 695 | Unfortunately, `scm_seed_to_random_state' only preserves around 32 |
| 696 | bits of entropy from the provided seed. I don't know if it's okay |
| 697 | to fix that in 2.0, so for now we have this workaround. */ |
| 698 | { |
| 699 | int i, len; |
| 700 | unsigned char *buf; |
| 701 | len = scm_to_int (scm_ceiling_quotient (scm_integer_length (seed), |
| 702 | SCM_I_MAKINUM (8))); |
| 703 | buf = (unsigned char *) malloc (len); |
| 704 | for (i = len-1; i >= 0; --i) |
| 705 | { |
| 706 | buf[i] = scm_to_int (scm_logand (seed, SCM_I_MAKINUM (255))); |
| 707 | seed = scm_ash (seed, SCM_I_MAKINUM (-8)); |
| 708 | } |
| 709 | state = make_rstate (scm_c_make_rstate ((char *) buf, len)); |
| 710 | free (buf); |
| 711 | } |
| 712 | return state; |
| 713 | } |
| 714 | |
| 715 | /* Attempt to fill buffer with random bytes from /dev/urandom. |
| 716 | Return 1 if successful, else return 0. */ |
| 717 | static int |
| 718 | read_dev_urandom (unsigned char *buf, size_t len) |
| 719 | { |
| 720 | size_t res = 0; |
| 721 | FILE *f = fopen ("/dev/urandom", "r"); |
| 722 | if (f) |
| 723 | { |
| 724 | res = fread(buf, 1, len, f); |
| 725 | fclose (f); |
| 726 | } |
| 727 | return (res == len); |
| 728 | } |
| 729 | |
| 730 | /* Fill a buffer with random bytes seeded from a platform-specific |
| 731 | source of entropy. /dev/urandom is used if available. Note that |
| 732 | this function provides no guarantees about the amount of entropy |
| 733 | present in the returned bytes. */ |
| 734 | void |
| 735 | scm_i_random_bytes_from_platform (unsigned char *buf, size_t len) |
| 736 | { |
| 737 | if (read_dev_urandom (buf, len)) |
| 738 | return; |
| 739 | else /* FIXME: support other platform sources */ |
| 740 | { |
| 741 | /* When all else fails, use this (rather weak) fallback */ |
| 742 | SCM random_state = random_state_of_last_resort (); |
| 743 | int i; |
| 744 | for (i = len-1; i >= 0; --i) |
| 745 | buf[i] = scm_to_int (scm_random (SCM_I_MAKINUM (256), random_state)); |
| 746 | } |
| 747 | } |
| 748 | |
| 749 | SCM_DEFINE (scm_random_state_from_platform, "random-state-from-platform", 0, 0, 0, |
| 750 | (void), |
| 751 | "Construct a new random state seeded from a platform-specific\n\ |
| 752 | source of entropy, appropriate for use in non-security-critical applications.") |
| 753 | #define FUNC_NAME s_scm_random_state_from_platform |
| 754 | { |
| 755 | unsigned char buf[32]; |
| 756 | if (read_dev_urandom (buf, sizeof(buf))) |
| 757 | return make_rstate (scm_c_make_rstate ((char *) buf, sizeof(buf))); |
| 758 | else |
| 759 | return random_state_of_last_resort (); |
| 760 | } |
| 761 | #undef FUNC_NAME |
| 762 | |
| 763 | void |
| 764 | scm_init_random () |
| 765 | { |
| 766 | int i, m; |
| 767 | /* plug in default RNG */ |
| 768 | scm_t_rng rng = |
| 769 | { |
| 770 | sizeof (scm_t_i_rstate), |
| 771 | scm_i_uniform32, |
| 772 | scm_i_init_rstate, |
| 773 | scm_i_copy_rstate, |
| 774 | scm_i_rstate_from_datum, |
| 775 | scm_i_rstate_to_datum |
| 776 | }; |
| 777 | scm_the_rng = rng; |
| 778 | |
| 779 | scm_tc16_rstate = scm_make_smob_type ("random-state", 0); |
| 780 | |
| 781 | for (m = 1; m <= 0x100; m <<= 1) |
| 782 | for (i = m >> 1; i < m; ++i) |
| 783 | scm_masktab[i] = m - 1; |
| 784 | |
| 785 | #include "libguile/random.x" |
| 786 | |
| 787 | scm_add_feature ("random"); |
| 788 | } |
| 789 | |
| 790 | /* |
| 791 | Local Variables: |
| 792 | c-file-style: "gnu" |
| 793 | End: |
| 794 | */ |