Commit | Line | Data |
---|---|---|
77b13912 | 1 | /* Copyright (C) 1999,2000,2001, 2003, 2005, 2006, 2009, 2010 Free Software Foundation, Inc. |
73be1d9e | 2 | * This library is free software; you can redistribute it and/or |
53befeb7 NJ |
3 | * modify it under the terms of the GNU Lesser General Public License |
4 | * as published by the Free Software Foundation; either version 3 of | |
5 | * the License, or (at your option) any later version. | |
e7a72986 | 6 | * |
53befeb7 NJ |
7 | * This library is distributed in the hope that it will be useful, but |
8 | * WITHOUT ANY WARRANTY; without even the implied warranty of | |
73be1d9e MV |
9 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
10 | * Lesser General Public License for more details. | |
e7a72986 | 11 | * |
73be1d9e MV |
12 | * You should have received a copy of the GNU Lesser General Public |
13 | * License along with this library; if not, write to the Free Software | |
53befeb7 NJ |
14 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
15 | * 02110-1301 USA | |
73be1d9e | 16 | */ |
e7a72986 | 17 | |
1bbd0b84 GB |
18 | |
19 | ||
e7a72986 MD |
20 | /* Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */ |
21 | ||
dbb605f5 | 22 | #ifdef HAVE_CONFIG_H |
2fb8013c RB |
23 | # include <config.h> |
24 | #endif | |
25 | ||
a0599745 | 26 | #include "libguile/_scm.h" |
e7a72986 | 27 | |
8db9cc6c | 28 | #include <gmp.h> |
7f146094 | 29 | #include <stdio.h> |
e7a72986 | 30 | #include <math.h> |
f34d19c7 | 31 | #include <string.h> |
a0599745 MD |
32 | #include "libguile/smob.h" |
33 | #include "libguile/numbers.h" | |
34 | #include "libguile/feature.h" | |
35 | #include "libguile/strings.h" | |
2fa901a5 | 36 | #include "libguile/arrays.h" |
46d25cff | 37 | #include "libguile/srfi-4.h" |
a0599745 | 38 | #include "libguile/vectors.h" |
f332e957 | 39 | #include "libguile/generalized-vectors.h" |
e7a72986 | 40 | |
a0599745 MD |
41 | #include "libguile/validate.h" |
42 | #include "libguile/random.h" | |
e7a72986 MD |
43 | |
44 | \f | |
45 | /* | |
46 | * A plugin interface for RNGs | |
47 | * | |
48 | * Using this interface, it is possible for the application to tell | |
49 | * libguile to use a different RNG. This is desirable if it is | |
50 | * necessary to use the same RNG everywhere in the application in | |
51 | * order to prevent interference, if the application uses RNG | |
52 | * hardware, or if the application has special demands on the RNG. | |
53 | * | |
54 | * Look in random.h and how the default generator is "plugged in" in | |
55 | * scm_init_random(). | |
56 | */ | |
57 | ||
92c2555f | 58 | scm_t_rng scm_the_rng; |
e7a72986 MD |
59 | |
60 | \f | |
61 | /* | |
62 | * The prepackaged RNG | |
63 | * | |
64 | * This is the MWC (Multiply With Carry) random number generator | |
65 | * described by George Marsaglia at the Department of Statistics and | |
66 | * Supercomputer Computations Research Institute, The Florida State | |
67 | * University (http://stat.fsu.edu/~geo). | |
68 | * | |
69 | * It uses 64 bits, has a period of 4578426017172946943 (4.6e18), and | |
70 | * passes all tests in the DIEHARD test suite | |
71 | * (http://stat.fsu.edu/~geo/diehard.html) | |
72 | */ | |
73 | ||
99a0ee66 AW |
74 | typedef struct scm_t_i_rstate { |
75 | scm_t_rstate rstate; | |
76 | scm_t_uint32 w; | |
77 | scm_t_uint32 c; | |
78 | } scm_t_i_rstate; | |
79 | ||
80 | ||
e7a72986 MD |
81 | #define A 2131995753UL |
82 | ||
82893676 MG |
83 | #ifndef M_PI |
84 | #define M_PI 3.14159265359 | |
85 | #endif | |
86 | ||
99a0ee66 AW |
87 | static scm_t_uint32 |
88 | scm_i_uniform32 (scm_t_rstate *state) | |
e7a72986 | 89 | { |
99a0ee66 AW |
90 | scm_t_i_rstate *istate = (scm_t_i_rstate*) state; |
91 | scm_t_uint64 x = (scm_t_uint64) A * istate->w + istate->c; | |
4a9f83ff | 92 | scm_t_uint32 w = x & 0xffffffffUL; |
99a0ee66 AW |
93 | istate->w = w; |
94 | istate->c = x >> 32L; | |
e7a72986 MD |
95 | return w; |
96 | } | |
97 | ||
99a0ee66 AW |
98 | static void |
99 | scm_i_init_rstate (scm_t_rstate *state, const char *seed, int n) | |
e7a72986 | 100 | { |
99a0ee66 | 101 | scm_t_i_rstate *istate = (scm_t_i_rstate*) state; |
4a9f83ff MD |
102 | scm_t_uint32 w = 0L; |
103 | scm_t_uint32 c = 0L; | |
e7a72986 MD |
104 | int i, m; |
105 | for (i = 0; i < n; ++i) | |
106 | { | |
107 | m = i % 8; | |
108 | if (m < 4) | |
109 | w += seed[i] << (8 * m); | |
110 | else | |
111 | c += seed[i] << (8 * (m - 4)); | |
112 | } | |
8b3747f9 | 113 | if ((w == 0 && c == 0) || (w == -1 && c == A - 1)) |
e7a72986 | 114 | ++c; |
99a0ee66 AW |
115 | istate->w = w; |
116 | istate->c = c; | |
e7a72986 MD |
117 | } |
118 | ||
99a0ee66 AW |
119 | static scm_t_rstate * |
120 | scm_i_copy_rstate (scm_t_rstate *state) | |
e7a72986 | 121 | { |
92d8fd32 LC |
122 | scm_t_rstate *new_state; |
123 | ||
124 | new_state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size, | |
125 | "random-state"); | |
e7a72986 MD |
126 | return memcpy (new_state, state, scm_the_rng.rstate_size); |
127 | } | |
128 | ||
77b13912 AR |
129 | SCM_SYMBOL(scm_i_rstate_tag, "multiply-with-carry"); |
130 | ||
99a0ee66 AW |
131 | static void |
132 | scm_i_rstate_from_datum (scm_t_rstate *state, SCM value) | |
133 | #define FUNC_NAME "scm_i_rstate_from_datum" | |
77b13912 | 134 | { |
99a0ee66 | 135 | scm_t_i_rstate *istate = (scm_t_i_rstate*) state; |
b606ff6a | 136 | scm_t_uint32 w, c; |
77b13912 AR |
137 | long length; |
138 | ||
139 | SCM_VALIDATE_LIST_COPYLEN (SCM_ARG1, value, length); | |
140 | SCM_ASSERT (length == 3, value, SCM_ARG1, FUNC_NAME); | |
141 | SCM_ASSERT (scm_is_eq (SCM_CAR (value), scm_i_rstate_tag), | |
142 | value, SCM_ARG1, FUNC_NAME); | |
b606ff6a AW |
143 | SCM_VALIDATE_UINT_COPY (SCM_ARG1, SCM_CADR (value), w); |
144 | SCM_VALIDATE_UINT_COPY (SCM_ARG1, SCM_CADDR (value), c); | |
77b13912 | 145 | |
99a0ee66 AW |
146 | istate->w = w; |
147 | istate->c = c; | |
77b13912 AR |
148 | } |
149 | #undef FUNC_NAME | |
150 | ||
99a0ee66 AW |
151 | static SCM |
152 | scm_i_rstate_to_datum (scm_t_rstate *state) | |
77b13912 | 153 | { |
99a0ee66 | 154 | scm_t_i_rstate *istate = (scm_t_i_rstate*) state; |
77b13912 | 155 | return scm_list_3 (scm_i_rstate_tag, |
99a0ee66 AW |
156 | scm_from_uint32 (istate->w), |
157 | scm_from_uint32 (istate->c)); | |
77b13912 AR |
158 | } |
159 | ||
e7a72986 MD |
160 | \f |
161 | /* | |
162 | * Random number library functions | |
163 | */ | |
164 | ||
92c2555f | 165 | scm_t_rstate * |
cc95e00a | 166 | scm_c_make_rstate (const char *seed, int n) |
5ee11b7c | 167 | { |
92d8fd32 LC |
168 | scm_t_rstate *state; |
169 | ||
170 | state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size, | |
171 | "random-state"); | |
5ee11b7c MD |
172 | state->reserved0 = 0; |
173 | scm_the_rng.init_rstate (state, seed, n); | |
174 | return state; | |
175 | } | |
176 | ||
77b13912 | 177 | scm_t_rstate * |
99a0ee66 | 178 | scm_c_rstate_from_datum (SCM datum) |
77b13912 AR |
179 | { |
180 | scm_t_rstate *state; | |
181 | ||
182 | state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size, | |
183 | "random-state"); | |
184 | state->reserved0 = 0; | |
99a0ee66 | 185 | scm_the_rng.from_datum (state, datum); |
77b13912 AR |
186 | return state; |
187 | } | |
2ade72d7 | 188 | |
92c2555f | 189 | scm_t_rstate * |
9b741bb6 | 190 | scm_c_default_rstate () |
2ade72d7 | 191 | #define FUNC_NAME "scm_c_default_rstate" |
9b741bb6 | 192 | { |
c5f268f8 | 193 | SCM state = SCM_VARIABLE_REF (scm_var_random_state); |
2ade72d7 DH |
194 | if (!SCM_RSTATEP (state)) |
195 | SCM_MISC_ERROR ("*random-state* contains bogus random state", SCM_EOL); | |
9b741bb6 MD |
196 | return SCM_RSTATE (state); |
197 | } | |
2ade72d7 DH |
198 | #undef FUNC_NAME |
199 | ||
9b741bb6 | 200 | |
e7a72986 | 201 | inline double |
92c2555f | 202 | scm_c_uniform01 (scm_t_rstate *state) |
e7a72986 | 203 | { |
5a92ddfd | 204 | double x = (double) scm_the_rng.random_bits (state) / (double) 0xffffffffUL; |
e7a72986 | 205 | return ((x + (double) scm_the_rng.random_bits (state)) |
5a92ddfd | 206 | / (double) 0xffffffffUL); |
e7a72986 MD |
207 | } |
208 | ||
209 | double | |
92c2555f | 210 | scm_c_normal01 (scm_t_rstate *state) |
e7a72986 MD |
211 | { |
212 | if (state->reserved0) | |
213 | { | |
214 | state->reserved0 = 0; | |
215 | return state->reserved1; | |
216 | } | |
217 | else | |
218 | { | |
219 | double r, a, n; | |
e7a72986 | 220 | |
9b741bb6 MD |
221 | r = sqrt (-2.0 * log (scm_c_uniform01 (state))); |
222 | a = 2.0 * M_PI * scm_c_uniform01 (state); | |
e7a72986 MD |
223 | |
224 | n = r * sin (a); | |
225 | state->reserved1 = r * cos (a); | |
5a92ddfd | 226 | state->reserved0 = 1; |
e7a72986 MD |
227 | |
228 | return n; | |
229 | } | |
230 | } | |
231 | ||
232 | double | |
92c2555f | 233 | scm_c_exp1 (scm_t_rstate *state) |
e7a72986 | 234 | { |
9b741bb6 | 235 | return - log (scm_c_uniform01 (state)); |
e7a72986 MD |
236 | } |
237 | ||
238 | unsigned char scm_masktab[256]; | |
239 | ||
b606ff6a AW |
240 | scm_t_uint32 |
241 | scm_c_random (scm_t_rstate *state, scm_t_uint32 m) | |
e7a72986 | 242 | { |
b606ff6a | 243 | scm_t_uint32 r, mask; |
e7a72986 MD |
244 | mask = (m < 0x100 |
245 | ? scm_masktab[m] | |
246 | : (m < 0x10000 | |
5a92ddfd | 247 | ? scm_masktab[m >> 8] << 8 | 0xff |
e7a72986 | 248 | : (m < 0x1000000 |
5a92ddfd MD |
249 | ? scm_masktab[m >> 16] << 16 | 0xffff |
250 | : scm_masktab[m >> 24] << 24 | 0xffffff))); | |
e7a72986 MD |
251 | while ((r = scm_the_rng.random_bits (state) & mask) >= m); |
252 | return r; | |
253 | } | |
254 | ||
ea7f6344 MD |
255 | /* |
256 | SCM scm_c_random_bignum (scm_t_rstate *state, SCM m) | |
257 | ||
258 | Takes a random state (source of random bits) and a bignum m. | |
259 | Returns a bignum b, 0 <= b < m. | |
260 | ||
261 | It does this by allocating a bignum b with as many base 65536 digits | |
262 | as m, filling b with random bits (in 32 bit chunks) up to the most | |
263 | significant 1 in m, and, finally checking if the resultant b is too | |
264 | large (>= m). If too large, we simply repeat the process again. (It | |
265 | is important to throw away all generated random bits if b >= m, | |
266 | otherwise we'll end up with a distorted distribution.) | |
267 | ||
268 | */ | |
269 | ||
e7a72986 | 270 | SCM |
92c2555f | 271 | scm_c_random_bignum (scm_t_rstate *state, SCM m) |
e7a72986 | 272 | { |
969d3bd0 | 273 | SCM result = scm_i_mkbig (); |
969d3bd0 | 274 | const size_t m_bits = mpz_sizeinbase (SCM_I_BIG_MPZ (m), 2); |
b606ff6a AW |
275 | /* how many bits would only partially fill the last scm_t_uint32? */ |
276 | const size_t end_bits = m_bits % (sizeof (scm_t_uint32) * SCM_CHAR_BIT); | |
277 | scm_t_uint32 *random_chunks = NULL; | |
278 | const scm_t_uint32 num_full_chunks = | |
279 | m_bits / (sizeof (scm_t_uint32) * SCM_CHAR_BIT); | |
280 | const scm_t_uint32 num_chunks = num_full_chunks + ((end_bits) ? 1 : 0); | |
969d3bd0 RB |
281 | |
282 | /* we know the result will be this big */ | |
283 | mpz_realloc2 (SCM_I_BIG_MPZ (result), m_bits); | |
284 | ||
285 | random_chunks = | |
b606ff6a | 286 | (scm_t_uint32 *) scm_gc_calloc (num_chunks * sizeof (scm_t_uint32), |
969d3bd0 RB |
287 | "random bignum chunks"); |
288 | ||
372691d8 | 289 | do |
e7a72986 | 290 | { |
b606ff6a AW |
291 | scm_t_uint32 *current_chunk = random_chunks + (num_chunks - 1); |
292 | scm_t_uint32 chunks_left = num_chunks; | |
969d3bd0 RB |
293 | |
294 | mpz_set_ui (SCM_I_BIG_MPZ (result), 0); | |
295 | ||
296 | if (end_bits) | |
297 | { | |
298 | /* generate a mask with ones in the end_bits position, i.e. if | |
299 | end_bits is 3, then we'd have a mask of ...0000000111 */ | |
b606ff6a AW |
300 | const scm_t_uint32 rndbits = scm_the_rng.random_bits (state); |
301 | int rshift = (sizeof (scm_t_uint32) * SCM_CHAR_BIT) - end_bits; | |
302 | scm_t_uint32 mask = ((scm_t_uint32)-1) >> rshift; | |
303 | scm_t_uint32 highest_bits = rndbits & mask; | |
969d3bd0 RB |
304 | *current_chunk-- = highest_bits; |
305 | chunks_left--; | |
306 | } | |
307 | ||
308 | while (chunks_left) | |
309 | { | |
b606ff6a | 310 | /* now fill in the remaining scm_t_uint32 sized chunks */ |
969d3bd0 RB |
311 | *current_chunk-- = scm_the_rng.random_bits (state); |
312 | chunks_left--; | |
313 | } | |
314 | mpz_import (SCM_I_BIG_MPZ (result), | |
315 | num_chunks, | |
316 | -1, | |
b606ff6a | 317 | sizeof (scm_t_uint32), |
969d3bd0 RB |
318 | 0, |
319 | 0, | |
320 | random_chunks); | |
372691d8 MD |
321 | /* if result >= m, regenerate it (it is important to regenerate |
322 | all bits in order not to get a distorted distribution) */ | |
323 | } while (mpz_cmp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (m)) >= 0); | |
969d3bd0 | 324 | scm_gc_free (random_chunks, |
b606ff6a | 325 | num_chunks * sizeof (scm_t_uint32), |
969d3bd0 | 326 | "random bignum chunks"); |
8db9cc6c | 327 | return scm_i_normbig (result); |
e7a72986 MD |
328 | } |
329 | ||
330 | /* | |
331 | * Scheme level representation of random states. | |
332 | */ | |
333 | ||
92c2555f | 334 | scm_t_bits scm_tc16_rstate; |
e7a72986 MD |
335 | |
336 | static SCM | |
92c2555f | 337 | make_rstate (scm_t_rstate *state) |
e7a72986 | 338 | { |
23a62151 | 339 | SCM_RETURN_NEWSMOB (scm_tc16_rstate, state); |
e7a72986 MD |
340 | } |
341 | ||
e7a72986 | 342 | |
e7a72986 MD |
343 | /* |
344 | * Scheme level interface. | |
345 | */ | |
346 | ||
cc95e00a | 347 | 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 | 348 | |
a1ec6916 | 349 | SCM_DEFINE (scm_random, "random", 1, 1, 0, |
1bbd0b84 | 350 | (SCM n, SCM state), |
34d19ef6 | 351 | "Return a number in [0, N).\n" |
d928e0b4 | 352 | "\n" |
9401323e NJ |
353 | "Accepts a positive integer or real n and returns a\n" |
354 | "number of the same type between zero (inclusive) and\n" | |
355 | "N (exclusive). The values returned have a uniform\n" | |
d928e0b4 GB |
356 | "distribution.\n" |
357 | "\n" | |
3b644514 MG |
358 | "The optional argument @var{state} must be of the type produced\n" |
359 | "by @code{seed->random-state}. It defaults to the value of the\n" | |
360 | "variable @var{*random-state*}. This object is used to maintain\n" | |
361 | "the state of the pseudo-random-number generator and is altered\n" | |
362 | "as a side effect of the random operation.") | |
1bbd0b84 | 363 | #define FUNC_NAME s_scm_random |
e7a72986 MD |
364 | { |
365 | if (SCM_UNBNDP (state)) | |
86d31dfe | 366 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 367 | SCM_VALIDATE_RSTATE (2, state); |
e11e83f3 | 368 | if (SCM_I_INUMP (n)) |
e7a72986 | 369 | { |
b606ff6a | 370 | scm_t_uint32 m = SCM_I_INUM (n); |
34d19ef6 | 371 | SCM_ASSERT_RANGE (1, n, m > 0); |
b606ff6a | 372 | return scm_from_uint32 (scm_c_random (SCM_RSTATE (state), m)); |
e7a72986 | 373 | } |
34d19ef6 | 374 | SCM_VALIDATE_NIM (1, n); |
e7a72986 | 375 | if (SCM_REALP (n)) |
d9a67fc4 MV |
376 | return scm_from_double (SCM_REAL_VALUE (n) |
377 | * scm_c_uniform01 (SCM_RSTATE (state))); | |
969d3bd0 | 378 | |
a55c2b68 MV |
379 | if (!SCM_BIGP (n)) |
380 | SCM_WRONG_TYPE_ARG (1, n); | |
9b741bb6 | 381 | return scm_c_random_bignum (SCM_RSTATE (state), n); |
e7a72986 | 382 | } |
1bbd0b84 | 383 | #undef FUNC_NAME |
e7a72986 | 384 | |
a1ec6916 | 385 | SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0, |
1bbd0b84 | 386 | (SCM state), |
3b644514 | 387 | "Return a copy of the random state @var{state}.") |
1bbd0b84 | 388 | #define FUNC_NAME s_scm_copy_random_state |
e7a72986 MD |
389 | { |
390 | if (SCM_UNBNDP (state)) | |
86d31dfe | 391 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 392 | SCM_VALIDATE_RSTATE (1, state); |
e7a72986 MD |
393 | return make_rstate (scm_the_rng.copy_rstate (SCM_RSTATE (state))); |
394 | } | |
1bbd0b84 | 395 | #undef FUNC_NAME |
e7a72986 | 396 | |
a1ec6916 | 397 | SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0, |
1bbd0b84 | 398 | (SCM seed), |
3b644514 | 399 | "Return a new random state using @var{seed}.") |
1bbd0b84 | 400 | #define FUNC_NAME s_scm_seed_to_random_state |
5ee11b7c | 401 | { |
8824ac88 | 402 | SCM res; |
5ee11b7c MD |
403 | if (SCM_NUMBERP (seed)) |
404 | seed = scm_number_to_string (seed, SCM_UNDEFINED); | |
34d19ef6 | 405 | SCM_VALIDATE_STRING (1, seed); |
cc95e00a MV |
406 | res = make_rstate (scm_c_make_rstate (scm_i_string_chars (seed), |
407 | scm_i_string_length (seed))); | |
8824ac88 MV |
408 | scm_remember_upto_here_1 (seed); |
409 | return res; | |
410 | ||
5ee11b7c | 411 | } |
1bbd0b84 | 412 | #undef FUNC_NAME |
5ee11b7c | 413 | |
99a0ee66 AW |
414 | SCM_DEFINE (scm_datum_to_random_state, "datum->random-state", 1, 0, 0, |
415 | (SCM datum), | |
1d454874 AW |
416 | "Return a new random state using @var{datum}, which should have\n" |
417 | "been obtailed from @code{random-state->datum}.") | |
99a0ee66 | 418 | #define FUNC_NAME s_scm_datum_to_random_state |
77b13912 | 419 | { |
99a0ee66 | 420 | return make_rstate (scm_c_rstate_from_datum (datum)); |
77b13912 AR |
421 | } |
422 | #undef FUNC_NAME | |
423 | ||
99a0ee66 | 424 | SCM_DEFINE (scm_random_state_to_datum, "random-state->datum", 1, 0, 0, |
77b13912 | 425 | (SCM state), |
99a0ee66 AW |
426 | "Return a datum representation of @var{state} that may be\n" |
427 | "written out and read back with the Scheme reader.") | |
428 | #define FUNC_NAME s_scm_random_state_to_datum | |
77b13912 AR |
429 | { |
430 | SCM_VALIDATE_RSTATE (1, state); | |
99a0ee66 | 431 | return scm_the_rng.to_datum (SCM_RSTATE (state)); |
77b13912 AR |
432 | } |
433 | #undef FUNC_NAME | |
434 | ||
a1ec6916 | 435 | SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0, |
1bbd0b84 | 436 | (SCM state), |
1e6808ea MG |
437 | "Return a uniformly distributed inexact real random number in\n" |
438 | "[0,1).") | |
1bbd0b84 | 439 | #define FUNC_NAME s_scm_random_uniform |
e7a72986 MD |
440 | { |
441 | if (SCM_UNBNDP (state)) | |
86d31dfe | 442 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 443 | SCM_VALIDATE_RSTATE (1, state); |
d9a67fc4 | 444 | return scm_from_double (scm_c_uniform01 (SCM_RSTATE (state))); |
e7a72986 | 445 | } |
1bbd0b84 | 446 | #undef FUNC_NAME |
e7a72986 | 447 | |
a1ec6916 | 448 | SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0, |
1bbd0b84 | 449 | (SCM state), |
1e6808ea MG |
450 | "Return an inexact real in a normal distribution. The\n" |
451 | "distribution used has mean 0 and standard deviation 1. For a\n" | |
452 | "normal distribution with mean m and standard deviation d use\n" | |
453 | "@code{(+ m (* d (random:normal)))}.") | |
1bbd0b84 | 454 | #define FUNC_NAME s_scm_random_normal |
afe5177e GH |
455 | { |
456 | if (SCM_UNBNDP (state)) | |
86d31dfe | 457 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 458 | SCM_VALIDATE_RSTATE (1, state); |
d9a67fc4 | 459 | return scm_from_double (scm_c_normal01 (SCM_RSTATE (state))); |
afe5177e | 460 | } |
1bbd0b84 | 461 | #undef FUNC_NAME |
afe5177e | 462 | |
e7a72986 | 463 | static void |
46d25cff | 464 | vector_scale_x (SCM v, double c) |
e7a72986 | 465 | { |
46d25cff | 466 | size_t n; |
4057a3e0 | 467 | if (scm_is_simple_vector (v)) |
46d25cff | 468 | { |
4057a3e0 | 469 | n = SCM_SIMPLE_VECTOR_LENGTH (v); |
46d25cff | 470 | while (n-- > 0) |
4057a3e0 | 471 | SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n)) *= c; |
46d25cff | 472 | } |
e7a72986 | 473 | else |
46d25cff MV |
474 | { |
475 | /* must be a f64vector. */ | |
4057a3e0 MV |
476 | scm_t_array_handle handle; |
477 | size_t i, len; | |
478 | ssize_t inc; | |
479 | double *elts; | |
480 | ||
481 | elts = scm_f64vector_writable_elements (v, &handle, &len, &inc); | |
482 | ||
483 | for (i = 0; i < len; i++, elts += inc) | |
484 | *elts *= c; | |
c8857a4d MV |
485 | |
486 | scm_array_handle_release (&handle); | |
46d25cff | 487 | } |
e7a72986 MD |
488 | } |
489 | ||
490 | static double | |
491 | vector_sum_squares (SCM v) | |
492 | { | |
493 | double x, sum = 0.0; | |
46d25cff | 494 | size_t n; |
4057a3e0 | 495 | if (scm_is_simple_vector (v)) |
46d25cff | 496 | { |
4057a3e0 | 497 | n = SCM_SIMPLE_VECTOR_LENGTH (v); |
46d25cff MV |
498 | while (n-- > 0) |
499 | { | |
4057a3e0 | 500 | x = SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n)); |
46d25cff MV |
501 | sum += x * x; |
502 | } | |
503 | } | |
e7a72986 | 504 | else |
46d25cff MV |
505 | { |
506 | /* must be a f64vector. */ | |
4057a3e0 MV |
507 | scm_t_array_handle handle; |
508 | size_t i, len; | |
509 | ssize_t inc; | |
510 | const double *elts; | |
511 | ||
512 | elts = scm_f64vector_elements (v, &handle, &len, &inc); | |
513 | ||
514 | for (i = 0; i < len; i++, elts += inc) | |
46d25cff | 515 | { |
4057a3e0 | 516 | x = *elts; |
46d25cff MV |
517 | sum += x * x; |
518 | } | |
4057a3e0 | 519 | |
c8857a4d | 520 | scm_array_handle_release (&handle); |
46d25cff | 521 | } |
e7a72986 MD |
522 | return sum; |
523 | } | |
524 | ||
e7a72986 MD |
525 | /* For the uniform distribution on the solid sphere, note that in |
526 | * this distribution the length r of the vector has cumulative | |
527 | * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be | |
528 | * generated as r=u^(1/n). | |
529 | */ | |
a1ec6916 | 530 | SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0, |
1bbd0b84 | 531 | (SCM v, SCM state), |
6efaeb35 KR |
532 | "Fills @var{vect} with inexact real random numbers the sum of\n" |
533 | "whose squares is less than 1.0. Thinking of @var{vect} as\n" | |
534 | "coordinates in space of dimension @var{n} @math{=}\n" | |
535 | "@code{(vector-length @var{vect})}, the coordinates are\n" | |
536 | "uniformly distributed within the unit @var{n}-sphere.") | |
1bbd0b84 | 537 | #define FUNC_NAME s_scm_random_solid_sphere_x |
e7a72986 | 538 | { |
e7a72986 | 539 | if (SCM_UNBNDP (state)) |
86d31dfe | 540 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 541 | SCM_VALIDATE_RSTATE (2, state); |
e7a72986 | 542 | scm_random_normal_vector_x (v, state); |
46d25cff MV |
543 | vector_scale_x (v, |
544 | pow (scm_c_uniform01 (SCM_RSTATE (state)), | |
f160e709 | 545 | 1.0 / scm_c_generalized_vector_length (v)) |
46d25cff | 546 | / sqrt (vector_sum_squares (v))); |
e7a72986 MD |
547 | return SCM_UNSPECIFIED; |
548 | } | |
1bbd0b84 | 549 | #undef FUNC_NAME |
e7a72986 | 550 | |
a1ec6916 | 551 | SCM_DEFINE (scm_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0, |
1bbd0b84 | 552 | (SCM v, SCM state), |
d928e0b4 GB |
553 | "Fills vect with inexact real random numbers\n" |
554 | "the sum of whose squares is equal to 1.0.\n" | |
9401323e | 555 | "Thinking of vect as coordinates in space of\n" |
d928e0b4 | 556 | "dimension n = (vector-length vect), the coordinates\n" |
9401323e | 557 | "are uniformly distributed over the surface of the\n" |
72dd0a03 | 558 | "unit n-sphere.") |
1bbd0b84 | 559 | #define FUNC_NAME s_scm_random_hollow_sphere_x |
e7a72986 | 560 | { |
e7a72986 | 561 | if (SCM_UNBNDP (state)) |
86d31dfe | 562 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 563 | SCM_VALIDATE_RSTATE (2, state); |
e7a72986 | 564 | scm_random_normal_vector_x (v, state); |
46d25cff | 565 | vector_scale_x (v, 1 / sqrt (vector_sum_squares (v))); |
e7a72986 MD |
566 | return SCM_UNSPECIFIED; |
567 | } | |
1bbd0b84 | 568 | #undef FUNC_NAME |
e7a72986 | 569 | |
1bbd0b84 | 570 | |
a1ec6916 | 571 | SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0, |
1bbd0b84 | 572 | (SCM v, SCM state), |
d928e0b4 GB |
573 | "Fills vect with inexact real random numbers that are\n" |
574 | "independent and standard normally distributed\n" | |
64ba8e85 | 575 | "(i.e., with mean 0 and variance 1).") |
1bbd0b84 | 576 | #define FUNC_NAME s_scm_random_normal_vector_x |
e7a72986 | 577 | { |
4057a3e0 MV |
578 | long i; |
579 | scm_t_array_handle handle; | |
580 | scm_t_array_dim *dim; | |
46d25cff | 581 | |
e7a72986 | 582 | if (SCM_UNBNDP (state)) |
86d31dfe | 583 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 584 | SCM_VALIDATE_RSTATE (2, state); |
4057a3e0 | 585 | |
c8857a4d | 586 | scm_generalized_vector_get_handle (v, &handle); |
4057a3e0 MV |
587 | dim = scm_array_handle_dims (&handle); |
588 | ||
589 | if (scm_is_vector (v)) | |
46d25cff | 590 | { |
4057a3e0 MV |
591 | SCM *elts = scm_array_handle_writable_elements (&handle); |
592 | for (i = dim->lbnd; i <= dim->ubnd; i++, elts += dim->inc) | |
593 | *elts = scm_from_double (scm_c_normal01 (SCM_RSTATE (state))); | |
46d25cff | 594 | } |
e7a72986 | 595 | else |
46d25cff MV |
596 | { |
597 | /* must be a f64vector. */ | |
4057a3e0 MV |
598 | double *elts = scm_array_handle_f64_writable_elements (&handle); |
599 | for (i = dim->lbnd; i <= dim->ubnd; i++, elts += dim->inc) | |
600 | *elts = scm_c_normal01 (SCM_RSTATE (state)); | |
46d25cff | 601 | } |
4057a3e0 | 602 | |
c8857a4d MV |
603 | scm_array_handle_release (&handle); |
604 | ||
e7a72986 MD |
605 | return SCM_UNSPECIFIED; |
606 | } | |
1bbd0b84 | 607 | #undef FUNC_NAME |
e7a72986 | 608 | |
a1ec6916 | 609 | SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0, |
1bbd0b84 | 610 | (SCM state), |
1e6808ea MG |
611 | "Return an inexact real in an exponential distribution with mean\n" |
612 | "1. For an exponential distribution with mean u use (* u\n" | |
613 | "(random:exp)).") | |
1bbd0b84 | 614 | #define FUNC_NAME s_scm_random_exp |
e7a72986 MD |
615 | { |
616 | if (SCM_UNBNDP (state)) | |
86d31dfe | 617 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 618 | SCM_VALIDATE_RSTATE (1, state); |
d9a67fc4 | 619 | return scm_from_double (scm_c_exp1 (SCM_RSTATE (state))); |
e7a72986 | 620 | } |
1bbd0b84 | 621 | #undef FUNC_NAME |
e7a72986 MD |
622 | |
623 | void | |
624 | scm_init_random () | |
625 | { | |
626 | int i, m; | |
627 | /* plug in default RNG */ | |
92c2555f | 628 | scm_t_rng rng = |
e7a72986 | 629 | { |
92c2555f | 630 | sizeof (scm_t_i_rstate), |
99a0ee66 AW |
631 | scm_i_uniform32, |
632 | scm_i_init_rstate, | |
633 | scm_i_copy_rstate, | |
634 | scm_i_rstate_from_datum, | |
635 | scm_i_rstate_to_datum | |
e7a72986 MD |
636 | }; |
637 | scm_the_rng = rng; | |
638 | ||
e841c3e0 | 639 | scm_tc16_rstate = scm_make_smob_type ("random-state", 0); |
e7a72986 MD |
640 | |
641 | for (m = 1; m <= 0x100; m <<= 1) | |
642 | for (i = m >> 1; i < m; ++i) | |
643 | scm_masktab[i] = m - 1; | |
644 | ||
a0599745 | 645 | #include "libguile/random.x" |
e7a72986 MD |
646 | |
647 | scm_add_feature ("random"); | |
648 | } | |
89e00824 ML |
649 | |
650 | /* | |
651 | Local Variables: | |
652 | c-file-style: "gnu" | |
653 | End: | |
654 | */ |