Commit | Line | Data |
---|---|---|
be6a36a0 | 1 | /* Copyright (C) 1999, 2000, 2001, 2003, 2005, 2006, 2009, 2010, |
bc8e6d7d MW |
2 | * 2012, 2013, 2014 Free Software Foundation, Inc. |
3 | * | |
73be1d9e | 4 | * This library is free software; you can redistribute it and/or |
53befeb7 NJ |
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. | |
e7a72986 | 8 | * |
53befeb7 NJ |
9 | * This library is distributed in the hope that it will be useful, but |
10 | * WITHOUT ANY WARRANTY; without even the implied warranty of | |
73be1d9e MV |
11 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
12 | * Lesser General Public License for more details. | |
e7a72986 | 13 | * |
73be1d9e MV |
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 | |
53befeb7 NJ |
16 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
17 | * 02110-1301 USA | |
73be1d9e | 18 | */ |
e7a72986 | 19 | |
1bbd0b84 GB |
20 | |
21 | ||
587f4edd | 22 | /* Original Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */ |
e7a72986 | 23 | |
dbb605f5 | 24 | #ifdef HAVE_CONFIG_H |
2fb8013c RB |
25 | # include <config.h> |
26 | #endif | |
27 | ||
a0599745 | 28 | #include "libguile/_scm.h" |
e7a72986 | 29 | |
8db9cc6c | 30 | #include <gmp.h> |
7f146094 | 31 | #include <stdio.h> |
e7a72986 | 32 | #include <math.h> |
f34d19c7 | 33 | #include <string.h> |
587f4edd | 34 | #include <sys/types.h> |
587f4edd | 35 | #include <unistd.h> |
587f4edd | 36 | |
a0599745 MD |
37 | #include "libguile/smob.h" |
38 | #include "libguile/numbers.h" | |
39 | #include "libguile/feature.h" | |
40 | #include "libguile/strings.h" | |
2fa901a5 | 41 | #include "libguile/arrays.h" |
46d25cff | 42 | #include "libguile/srfi-4.h" |
a0599745 | 43 | #include "libguile/vectors.h" |
f332e957 | 44 | #include "libguile/generalized-vectors.h" |
e7a72986 | 45 | |
a0599745 MD |
46 | #include "libguile/validate.h" |
47 | #include "libguile/random.h" | |
e7a72986 MD |
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 | ||
92c2555f | 63 | scm_t_rng scm_the_rng; |
e7a72986 MD |
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 | ||
99a0ee66 AW |
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 | ||
e7a72986 MD |
86 | #define A 2131995753UL |
87 | ||
82893676 MG |
88 | #ifndef M_PI |
89 | #define M_PI 3.14159265359 | |
90 | #endif | |
91 | ||
99a0ee66 AW |
92 | static scm_t_uint32 |
93 | scm_i_uniform32 (scm_t_rstate *state) | |
e7a72986 | 94 | { |
99a0ee66 AW |
95 | scm_t_i_rstate *istate = (scm_t_i_rstate*) state; |
96 | scm_t_uint64 x = (scm_t_uint64) A * istate->w + istate->c; | |
4a9f83ff | 97 | scm_t_uint32 w = x & 0xffffffffUL; |
99a0ee66 AW |
98 | istate->w = w; |
99 | istate->c = x >> 32L; | |
e7a72986 MD |
100 | return w; |
101 | } | |
102 | ||
99a0ee66 AW |
103 | static void |
104 | scm_i_init_rstate (scm_t_rstate *state, const char *seed, int n) | |
e7a72986 | 105 | { |
99a0ee66 | 106 | scm_t_i_rstate *istate = (scm_t_i_rstate*) state; |
4a9f83ff MD |
107 | scm_t_uint32 w = 0L; |
108 | scm_t_uint32 c = 0L; | |
e7a72986 MD |
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 | } | |
8b3747f9 | 118 | if ((w == 0 && c == 0) || (w == -1 && c == A - 1)) |
e7a72986 | 119 | ++c; |
99a0ee66 AW |
120 | istate->w = w; |
121 | istate->c = c; | |
e7a72986 MD |
122 | } |
123 | ||
99a0ee66 AW |
124 | static scm_t_rstate * |
125 | scm_i_copy_rstate (scm_t_rstate *state) | |
e7a72986 | 126 | { |
92d8fd32 LC |
127 | scm_t_rstate *new_state; |
128 | ||
a2a95453 | 129 | new_state = scm_gc_malloc_pointerless (state->rng->rstate_size, |
92d8fd32 | 130 | "random-state"); |
a2a95453 | 131 | return memcpy (new_state, state, state->rng->rstate_size); |
e7a72986 MD |
132 | } |
133 | ||
77b13912 AR |
134 | SCM_SYMBOL(scm_i_rstate_tag, "multiply-with-carry"); |
135 | ||
99a0ee66 AW |
136 | static void |
137 | scm_i_rstate_from_datum (scm_t_rstate *state, SCM value) | |
138 | #define FUNC_NAME "scm_i_rstate_from_datum" | |
77b13912 | 139 | { |
99a0ee66 | 140 | scm_t_i_rstate *istate = (scm_t_i_rstate*) state; |
b606ff6a | 141 | scm_t_uint32 w, c; |
77b13912 AR |
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); | |
b606ff6a AW |
148 | SCM_VALIDATE_UINT_COPY (SCM_ARG1, SCM_CADR (value), w); |
149 | SCM_VALIDATE_UINT_COPY (SCM_ARG1, SCM_CADDR (value), c); | |
77b13912 | 150 | |
99a0ee66 AW |
151 | istate->w = w; |
152 | istate->c = c; | |
77b13912 AR |
153 | } |
154 | #undef FUNC_NAME | |
155 | ||
99a0ee66 AW |
156 | static SCM |
157 | scm_i_rstate_to_datum (scm_t_rstate *state) | |
77b13912 | 158 | { |
99a0ee66 | 159 | scm_t_i_rstate *istate = (scm_t_i_rstate*) state; |
77b13912 | 160 | return scm_list_3 (scm_i_rstate_tag, |
99a0ee66 AW |
161 | scm_from_uint32 (istate->w), |
162 | scm_from_uint32 (istate->c)); | |
77b13912 AR |
163 | } |
164 | ||
e7a72986 MD |
165 | \f |
166 | /* | |
167 | * Random number library functions | |
168 | */ | |
169 | ||
92c2555f | 170 | scm_t_rstate * |
cc95e00a | 171 | scm_c_make_rstate (const char *seed, int n) |
5ee11b7c | 172 | { |
92d8fd32 LC |
173 | scm_t_rstate *state; |
174 | ||
175 | state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size, | |
176 | "random-state"); | |
a2a95453 AW |
177 | state->rng = &scm_the_rng; |
178 | state->normal_next = 0.0; | |
179 | state->rng->init_rstate (state, seed, n); | |
5ee11b7c MD |
180 | return state; |
181 | } | |
182 | ||
77b13912 | 183 | scm_t_rstate * |
99a0ee66 | 184 | scm_c_rstate_from_datum (SCM datum) |
77b13912 AR |
185 | { |
186 | scm_t_rstate *state; | |
187 | ||
188 | state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size, | |
189 | "random-state"); | |
a2a95453 AW |
190 | state->rng = &scm_the_rng; |
191 | state->normal_next = 0.0; | |
192 | state->rng->from_datum (state, datum); | |
77b13912 AR |
193 | return state; |
194 | } | |
2ade72d7 | 195 | |
92c2555f | 196 | scm_t_rstate * |
9b741bb6 | 197 | scm_c_default_rstate () |
2ade72d7 | 198 | #define FUNC_NAME "scm_c_default_rstate" |
9b741bb6 | 199 | { |
c5f268f8 | 200 | SCM state = SCM_VARIABLE_REF (scm_var_random_state); |
2ade72d7 DH |
201 | if (!SCM_RSTATEP (state)) |
202 | SCM_MISC_ERROR ("*random-state* contains bogus random state", SCM_EOL); | |
9b741bb6 MD |
203 | return SCM_RSTATE (state); |
204 | } | |
2ade72d7 DH |
205 | #undef FUNC_NAME |
206 | ||
9b741bb6 | 207 | |
a2a95453 | 208 | double |
92c2555f | 209 | scm_c_uniform01 (scm_t_rstate *state) |
e7a72986 | 210 | { |
a2a95453 AW |
211 | double x = (double) state->rng->random_bits (state) / (double) 0xffffffffUL; |
212 | return ((x + (double) state->rng->random_bits (state)) | |
5a92ddfd | 213 | / (double) 0xffffffffUL); |
e7a72986 MD |
214 | } |
215 | ||
216 | double | |
92c2555f | 217 | scm_c_normal01 (scm_t_rstate *state) |
e7a72986 | 218 | { |
a2a95453 | 219 | if (state->normal_next != 0.0) |
e7a72986 | 220 | { |
a2a95453 AW |
221 | double ret = state->normal_next; |
222 | ||
223 | state->normal_next = 0.0; | |
224 | ||
225 | return ret; | |
e7a72986 MD |
226 | } |
227 | else | |
228 | { | |
229 | double r, a, n; | |
e7a72986 | 230 | |
9b741bb6 MD |
231 | r = sqrt (-2.0 * log (scm_c_uniform01 (state))); |
232 | a = 2.0 * M_PI * scm_c_uniform01 (state); | |
e7a72986 MD |
233 | |
234 | n = r * sin (a); | |
a2a95453 | 235 | state->normal_next = r * cos (a); |
e7a72986 MD |
236 | |
237 | return n; | |
238 | } | |
239 | } | |
240 | ||
241 | double | |
92c2555f | 242 | scm_c_exp1 (scm_t_rstate *state) |
e7a72986 | 243 | { |
9b741bb6 | 244 | return - log (scm_c_uniform01 (state)); |
e7a72986 MD |
245 | } |
246 | ||
247 | unsigned char scm_masktab[256]; | |
248 | ||
2af6e135 AR |
249 | static inline scm_t_uint32 |
250 | scm_i_mask32 (scm_t_uint32 m) | |
e7a72986 | 251 | { |
2af6e135 | 252 | return (m < 0x100 |
e7a72986 MD |
253 | ? scm_masktab[m] |
254 | : (m < 0x10000 | |
5a92ddfd | 255 | ? scm_masktab[m >> 8] << 8 | 0xff |
e7a72986 | 256 | : (m < 0x1000000 |
5a92ddfd MD |
257 | ? scm_masktab[m >> 16] << 16 | 0xffff |
258 | : scm_masktab[m >> 24] << 24 | 0xffffff))); | |
2af6e135 AR |
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); | |
a2a95453 | 265 | while ((r = state->rng->random_bits (state) & mask) >= m); |
e7a72986 MD |
266 | return r; |
267 | } | |
268 | ||
2af6e135 AR |
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 | ||
ea7f6344 MD |
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 | ||
e7a72986 | 300 | SCM |
92c2555f | 301 | scm_c_random_bignum (scm_t_rstate *state, SCM m) |
e7a72986 | 302 | { |
969d3bd0 | 303 | SCM result = scm_i_mkbig (); |
969d3bd0 | 304 | const size_t m_bits = mpz_sizeinbase (SCM_I_BIG_MPZ (m), 2); |
b606ff6a AW |
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); | |
969d3bd0 RB |
311 | |
312 | /* we know the result will be this big */ | |
313 | mpz_realloc2 (SCM_I_BIG_MPZ (result), m_bits); | |
314 | ||
315 | random_chunks = | |
b606ff6a | 316 | (scm_t_uint32 *) scm_gc_calloc (num_chunks * sizeof (scm_t_uint32), |
969d3bd0 RB |
317 | "random bignum chunks"); |
318 | ||
372691d8 | 319 | do |
e7a72986 | 320 | { |
b606ff6a AW |
321 | scm_t_uint32 *current_chunk = random_chunks + (num_chunks - 1); |
322 | scm_t_uint32 chunks_left = num_chunks; | |
969d3bd0 RB |
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 */ | |
a2a95453 | 330 | const scm_t_uint32 rndbits = state->rng->random_bits (state); |
b606ff6a AW |
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; | |
969d3bd0 RB |
334 | *current_chunk-- = highest_bits; |
335 | chunks_left--; | |
336 | } | |
337 | ||
338 | while (chunks_left) | |
339 | { | |
b606ff6a | 340 | /* now fill in the remaining scm_t_uint32 sized chunks */ |
a2a95453 | 341 | *current_chunk-- = state->rng->random_bits (state); |
969d3bd0 RB |
342 | chunks_left--; |
343 | } | |
344 | mpz_import (SCM_I_BIG_MPZ (result), | |
345 | num_chunks, | |
346 | -1, | |
b606ff6a | 347 | sizeof (scm_t_uint32), |
969d3bd0 RB |
348 | 0, |
349 | 0, | |
350 | random_chunks); | |
372691d8 MD |
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); | |
969d3bd0 | 354 | scm_gc_free (random_chunks, |
b606ff6a | 355 | num_chunks * sizeof (scm_t_uint32), |
969d3bd0 | 356 | "random bignum chunks"); |
8db9cc6c | 357 | return scm_i_normbig (result); |
e7a72986 MD |
358 | } |
359 | ||
360 | /* | |
361 | * Scheme level representation of random states. | |
362 | */ | |
363 | ||
92c2555f | 364 | scm_t_bits scm_tc16_rstate; |
e7a72986 MD |
365 | |
366 | static SCM | |
92c2555f | 367 | make_rstate (scm_t_rstate *state) |
e7a72986 | 368 | { |
23a62151 | 369 | SCM_RETURN_NEWSMOB (scm_tc16_rstate, state); |
e7a72986 MD |
370 | } |
371 | ||
e7a72986 | 372 | |
e7a72986 MD |
373 | /* |
374 | * Scheme level interface. | |
375 | */ | |
376 | ||
cc95e00a | 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"))); |
e7a72986 | 378 | |
a1ec6916 | 379 | SCM_DEFINE (scm_random, "random", 1, 1, 0, |
1bbd0b84 | 380 | (SCM n, SCM state), |
34d19ef6 | 381 | "Return a number in [0, N).\n" |
d928e0b4 | 382 | "\n" |
9401323e NJ |
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" | |
d928e0b4 GB |
386 | "distribution.\n" |
387 | "\n" | |
3b644514 MG |
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.") | |
1bbd0b84 | 393 | #define FUNC_NAME s_scm_random |
e7a72986 MD |
394 | { |
395 | if (SCM_UNBNDP (state)) | |
86d31dfe | 396 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 397 | SCM_VALIDATE_RSTATE (2, state); |
e11e83f3 | 398 | if (SCM_I_INUMP (n)) |
e7a72986 | 399 | { |
e25f3727 | 400 | scm_t_bits m = (scm_t_bits) SCM_I_INUM (n); |
9defb641 | 401 | SCM_ASSERT_RANGE (1, n, SCM_I_INUM (n) > 0); |
e25f3727 | 402 | #if SCM_SIZEOF_UINTPTR_T <= 4 |
9defb641 AW |
403 | return scm_from_uint32 (scm_c_random (SCM_RSTATE (state), |
404 | (scm_t_uint32) m)); | |
e25f3727 | 405 | #elif SCM_SIZEOF_UINTPTR_T <= 8 |
2af6e135 AR |
406 | return scm_from_uint64 (scm_c_random64 (SCM_RSTATE (state), |
407 | (scm_t_uint64) m)); | |
9defb641 | 408 | #else |
e25f3727 | 409 | #error "Cannot deal with this platform's scm_t_bits size" |
9defb641 | 410 | #endif |
e7a72986 | 411 | } |
34d19ef6 | 412 | SCM_VALIDATE_NIM (1, n); |
e7a72986 | 413 | if (SCM_REALP (n)) |
d9a67fc4 MV |
414 | return scm_from_double (SCM_REAL_VALUE (n) |
415 | * scm_c_uniform01 (SCM_RSTATE (state))); | |
969d3bd0 | 416 | |
a55c2b68 MV |
417 | if (!SCM_BIGP (n)) |
418 | SCM_WRONG_TYPE_ARG (1, n); | |
9b741bb6 | 419 | return scm_c_random_bignum (SCM_RSTATE (state), n); |
e7a72986 | 420 | } |
1bbd0b84 | 421 | #undef FUNC_NAME |
e7a72986 | 422 | |
a1ec6916 | 423 | SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0, |
1bbd0b84 | 424 | (SCM state), |
3b644514 | 425 | "Return a copy of the random state @var{state}.") |
1bbd0b84 | 426 | #define FUNC_NAME s_scm_copy_random_state |
e7a72986 MD |
427 | { |
428 | if (SCM_UNBNDP (state)) | |
86d31dfe | 429 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 430 | SCM_VALIDATE_RSTATE (1, state); |
a2a95453 | 431 | return make_rstate (SCM_RSTATE (state)->rng->copy_rstate (SCM_RSTATE (state))); |
e7a72986 | 432 | } |
1bbd0b84 | 433 | #undef FUNC_NAME |
e7a72986 | 434 | |
a1ec6916 | 435 | SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0, |
1bbd0b84 | 436 | (SCM seed), |
3b644514 | 437 | "Return a new random state using @var{seed}.") |
1bbd0b84 | 438 | #define FUNC_NAME s_scm_seed_to_random_state |
5ee11b7c | 439 | { |
8824ac88 | 440 | SCM res; |
5ee11b7c MD |
441 | if (SCM_NUMBERP (seed)) |
442 | seed = scm_number_to_string (seed, SCM_UNDEFINED); | |
34d19ef6 | 443 | SCM_VALIDATE_STRING (1, seed); |
cc95e00a MV |
444 | res = make_rstate (scm_c_make_rstate (scm_i_string_chars (seed), |
445 | scm_i_string_length (seed))); | |
8824ac88 MV |
446 | scm_remember_upto_here_1 (seed); |
447 | return res; | |
448 | ||
5ee11b7c | 449 | } |
1bbd0b84 | 450 | #undef FUNC_NAME |
5ee11b7c | 451 | |
99a0ee66 AW |
452 | SCM_DEFINE (scm_datum_to_random_state, "datum->random-state", 1, 0, 0, |
453 | (SCM datum), | |
1d454874 | 454 | "Return a new random state using @var{datum}, which should have\n" |
ffb62a43 | 455 | "been obtained from @code{random-state->datum}.") |
99a0ee66 | 456 | #define FUNC_NAME s_scm_datum_to_random_state |
77b13912 | 457 | { |
99a0ee66 | 458 | return make_rstate (scm_c_rstate_from_datum (datum)); |
77b13912 AR |
459 | } |
460 | #undef FUNC_NAME | |
461 | ||
99a0ee66 | 462 | SCM_DEFINE (scm_random_state_to_datum, "random-state->datum", 1, 0, 0, |
77b13912 | 463 | (SCM state), |
99a0ee66 AW |
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 | |
77b13912 AR |
467 | { |
468 | SCM_VALIDATE_RSTATE (1, state); | |
a2a95453 | 469 | return SCM_RSTATE (state)->rng->to_datum (SCM_RSTATE (state)); |
77b13912 AR |
470 | } |
471 | #undef FUNC_NAME | |
472 | ||
a1ec6916 | 473 | SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0, |
1bbd0b84 | 474 | (SCM state), |
1e6808ea MG |
475 | "Return a uniformly distributed inexact real random number in\n" |
476 | "[0,1).") | |
1bbd0b84 | 477 | #define FUNC_NAME s_scm_random_uniform |
e7a72986 MD |
478 | { |
479 | if (SCM_UNBNDP (state)) | |
86d31dfe | 480 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 481 | SCM_VALIDATE_RSTATE (1, state); |
d9a67fc4 | 482 | return scm_from_double (scm_c_uniform01 (SCM_RSTATE (state))); |
e7a72986 | 483 | } |
1bbd0b84 | 484 | #undef FUNC_NAME |
e7a72986 | 485 | |
a1ec6916 | 486 | SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0, |
1bbd0b84 | 487 | (SCM state), |
1e6808ea MG |
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)))}.") | |
1bbd0b84 | 492 | #define FUNC_NAME s_scm_random_normal |
afe5177e GH |
493 | { |
494 | if (SCM_UNBNDP (state)) | |
86d31dfe | 495 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 496 | SCM_VALIDATE_RSTATE (1, state); |
d9a67fc4 | 497 | return scm_from_double (scm_c_normal01 (SCM_RSTATE (state))); |
afe5177e | 498 | } |
1bbd0b84 | 499 | #undef FUNC_NAME |
afe5177e | 500 | |
e7a72986 | 501 | static void |
46d25cff | 502 | vector_scale_x (SCM v, double c) |
e7a72986 | 503 | { |
46d25cff | 504 | size_t n; |
4057a3e0 | 505 | if (scm_is_simple_vector (v)) |
46d25cff | 506 | { |
4057a3e0 | 507 | n = SCM_SIMPLE_VECTOR_LENGTH (v); |
46d25cff | 508 | while (n-- > 0) |
4057a3e0 | 509 | SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n)) *= c; |
46d25cff | 510 | } |
e7a72986 | 511 | else |
46d25cff MV |
512 | { |
513 | /* must be a f64vector. */ | |
4057a3e0 MV |
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; | |
c8857a4d MV |
523 | |
524 | scm_array_handle_release (&handle); | |
46d25cff | 525 | } |
e7a72986 MD |
526 | } |
527 | ||
528 | static double | |
529 | vector_sum_squares (SCM v) | |
530 | { | |
531 | double x, sum = 0.0; | |
46d25cff | 532 | size_t n; |
4057a3e0 | 533 | if (scm_is_simple_vector (v)) |
46d25cff | 534 | { |
4057a3e0 | 535 | n = SCM_SIMPLE_VECTOR_LENGTH (v); |
46d25cff MV |
536 | while (n-- > 0) |
537 | { | |
4057a3e0 | 538 | x = SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n)); |
46d25cff MV |
539 | sum += x * x; |
540 | } | |
541 | } | |
e7a72986 | 542 | else |
46d25cff MV |
543 | { |
544 | /* must be a f64vector. */ | |
4057a3e0 MV |
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) | |
46d25cff | 553 | { |
4057a3e0 | 554 | x = *elts; |
46d25cff MV |
555 | sum += x * x; |
556 | } | |
4057a3e0 | 557 | |
c8857a4d | 558 | scm_array_handle_release (&handle); |
46d25cff | 559 | } |
e7a72986 MD |
560 | return sum; |
561 | } | |
562 | ||
e7a72986 MD |
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 | */ | |
a1ec6916 | 568 | SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0, |
1bbd0b84 | 569 | (SCM v, SCM state), |
6efaeb35 KR |
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.") | |
1bbd0b84 | 575 | #define FUNC_NAME s_scm_random_solid_sphere_x |
e7a72986 | 576 | { |
e7a72986 | 577 | if (SCM_UNBNDP (state)) |
86d31dfe | 578 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 579 | SCM_VALIDATE_RSTATE (2, state); |
e7a72986 | 580 | scm_random_normal_vector_x (v, state); |
46d25cff MV |
581 | vector_scale_x (v, |
582 | pow (scm_c_uniform01 (SCM_RSTATE (state)), | |
f160e709 | 583 | 1.0 / scm_c_generalized_vector_length (v)) |
46d25cff | 584 | / sqrt (vector_sum_squares (v))); |
e7a72986 MD |
585 | return SCM_UNSPECIFIED; |
586 | } | |
1bbd0b84 | 587 | #undef FUNC_NAME |
e7a72986 | 588 | |
a1ec6916 | 589 | SCM_DEFINE (scm_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0, |
1bbd0b84 | 590 | (SCM v, SCM state), |
d928e0b4 GB |
591 | "Fills vect with inexact real random numbers\n" |
592 | "the sum of whose squares is equal to 1.0.\n" | |
9401323e | 593 | "Thinking of vect as coordinates in space of\n" |
d928e0b4 | 594 | "dimension n = (vector-length vect), the coordinates\n" |
9401323e | 595 | "are uniformly distributed over the surface of the\n" |
72dd0a03 | 596 | "unit n-sphere.") |
1bbd0b84 | 597 | #define FUNC_NAME s_scm_random_hollow_sphere_x |
e7a72986 | 598 | { |
e7a72986 | 599 | if (SCM_UNBNDP (state)) |
86d31dfe | 600 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 601 | SCM_VALIDATE_RSTATE (2, state); |
e7a72986 | 602 | scm_random_normal_vector_x (v, state); |
46d25cff | 603 | vector_scale_x (v, 1 / sqrt (vector_sum_squares (v))); |
e7a72986 MD |
604 | return SCM_UNSPECIFIED; |
605 | } | |
1bbd0b84 | 606 | #undef FUNC_NAME |
e7a72986 | 607 | |
1bbd0b84 | 608 | |
a1ec6916 | 609 | SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0, |
1bbd0b84 | 610 | (SCM v, SCM state), |
d928e0b4 GB |
611 | "Fills vect with inexact real random numbers that are\n" |
612 | "independent and standard normally distributed\n" | |
64ba8e85 | 613 | "(i.e., with mean 0 and variance 1).") |
1bbd0b84 | 614 | #define FUNC_NAME s_scm_random_normal_vector_x |
e7a72986 | 615 | { |
4057a3e0 MV |
616 | long i; |
617 | scm_t_array_handle handle; | |
618 | scm_t_array_dim *dim; | |
46d25cff | 619 | |
e7a72986 | 620 | if (SCM_UNBNDP (state)) |
86d31dfe | 621 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 622 | SCM_VALIDATE_RSTATE (2, state); |
4057a3e0 | 623 | |
c8857a4d | 624 | scm_generalized_vector_get_handle (v, &handle); |
4057a3e0 MV |
625 | dim = scm_array_handle_dims (&handle); |
626 | ||
627 | if (scm_is_vector (v)) | |
46d25cff | 628 | { |
4057a3e0 MV |
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))); | |
46d25cff | 632 | } |
e7a72986 | 633 | else |
46d25cff MV |
634 | { |
635 | /* must be a f64vector. */ | |
4057a3e0 MV |
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)); | |
46d25cff | 639 | } |
4057a3e0 | 640 | |
c8857a4d MV |
641 | scm_array_handle_release (&handle); |
642 | ||
e7a72986 MD |
643 | return SCM_UNSPECIFIED; |
644 | } | |
1bbd0b84 | 645 | #undef FUNC_NAME |
e7a72986 | 646 | |
a1ec6916 | 647 | SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0, |
1bbd0b84 | 648 | (SCM state), |
1e6808ea MG |
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)).") | |
1bbd0b84 | 652 | #define FUNC_NAME s_scm_random_exp |
e7a72986 MD |
653 | { |
654 | if (SCM_UNBNDP (state)) | |
86d31dfe | 655 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 656 | SCM_VALIDATE_RSTATE (1, state); |
d9a67fc4 | 657 | return scm_from_double (scm_c_exp1 (SCM_RSTATE (state))); |
e7a72986 | 658 | } |
1bbd0b84 | 659 | #undef FUNC_NAME |
e7a72986 | 660 | |
444b26f7 MW |
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. */ | |
d47db067 MW |
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 */ | |
587f4edd MW |
673 | /* Avoid scm_getpid, since it depends on HAVE_POSIX. */ |
674 | scm_from_unsigned_integer (getpid ()), /* process ID */ | |
d47db067 MW |
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); | |
444b26f7 | 680 | |
d47db067 | 681 | /* Concatenate the sources bitwise to form the seed */ |
3ec19a78 | 682 | SCM seed = SCM_INUM0; |
d47db067 MW |
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 | ||
e7a72986 MD |
763 | void |
764 | scm_init_random () | |
765 | { | |
766 | int i, m; | |
767 | /* plug in default RNG */ | |
92c2555f | 768 | scm_t_rng rng = |
e7a72986 | 769 | { |
92c2555f | 770 | sizeof (scm_t_i_rstate), |
99a0ee66 AW |
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 | |
e7a72986 MD |
776 | }; |
777 | scm_the_rng = rng; | |
778 | ||
e841c3e0 | 779 | scm_tc16_rstate = scm_make_smob_type ("random-state", 0); |
e7a72986 MD |
780 | |
781 | for (m = 1; m <= 0x100; m <<= 1) | |
782 | for (i = m >> 1; i < m; ++i) | |
783 | scm_masktab[i] = m - 1; | |
784 | ||
a0599745 | 785 | #include "libguile/random.x" |
e7a72986 MD |
786 | |
787 | scm_add_feature ("random"); | |
788 | } | |
89e00824 ML |
789 | |
790 | /* | |
791 | Local Variables: | |
792 | c-file-style: "gnu" | |
793 | End: | |
794 | */ |