Commit | Line | Data |
---|---|---|
ea7f6344 | 1 | /* Copyright (C) 1999,2000,2001, 2003 Free Software Foundation, Inc. |
73be1d9e MV |
2 | * This library is free software; you can redistribute it and/or |
3 | * modify it under the terms of the GNU Lesser General Public | |
4 | * License as published by the Free Software Foundation; either | |
5 | * version 2.1 of the License, or (at your option) any later version. | |
e7a72986 | 6 | * |
73be1d9e | 7 | * This library is distributed in the hope that it will be useful, |
e7a72986 | 8 | * but 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 | |
14 | * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA | |
15 | */ | |
e7a72986 | 16 | |
1bbd0b84 GB |
17 | |
18 | ||
e7a72986 MD |
19 | /* Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */ |
20 | ||
2fb8013c RB |
21 | #if HAVE_CONFIG_H |
22 | # include <config.h> | |
23 | #endif | |
24 | ||
a0599745 | 25 | #include "libguile/_scm.h" |
e7a72986 | 26 | |
8db9cc6c | 27 | #include <gmp.h> |
7f146094 | 28 | #include <stdio.h> |
e7a72986 | 29 | #include <math.h> |
f34d19c7 | 30 | #include <string.h> |
a0599745 MD |
31 | #include "libguile/smob.h" |
32 | #include "libguile/numbers.h" | |
33 | #include "libguile/feature.h" | |
34 | #include "libguile/strings.h" | |
e9bfab50 | 35 | #include "libguile/unif.h" |
a0599745 | 36 | #include "libguile/vectors.h" |
e7a72986 | 37 | |
a0599745 MD |
38 | #include "libguile/validate.h" |
39 | #include "libguile/random.h" | |
e7a72986 MD |
40 | |
41 | \f | |
42 | /* | |
43 | * A plugin interface for RNGs | |
44 | * | |
45 | * Using this interface, it is possible for the application to tell | |
46 | * libguile to use a different RNG. This is desirable if it is | |
47 | * necessary to use the same RNG everywhere in the application in | |
48 | * order to prevent interference, if the application uses RNG | |
49 | * hardware, or if the application has special demands on the RNG. | |
50 | * | |
51 | * Look in random.h and how the default generator is "plugged in" in | |
52 | * scm_init_random(). | |
53 | */ | |
54 | ||
92c2555f | 55 | scm_t_rng scm_the_rng; |
e7a72986 MD |
56 | |
57 | \f | |
58 | /* | |
59 | * The prepackaged RNG | |
60 | * | |
61 | * This is the MWC (Multiply With Carry) random number generator | |
62 | * described by George Marsaglia at the Department of Statistics and | |
63 | * Supercomputer Computations Research Institute, The Florida State | |
64 | * University (http://stat.fsu.edu/~geo). | |
65 | * | |
66 | * It uses 64 bits, has a period of 4578426017172946943 (4.6e18), and | |
67 | * passes all tests in the DIEHARD test suite | |
68 | * (http://stat.fsu.edu/~geo/diehard.html) | |
69 | */ | |
70 | ||
71 | #define A 2131995753UL | |
72 | ||
82893676 MG |
73 | #ifndef M_PI |
74 | #define M_PI 3.14159265359 | |
75 | #endif | |
76 | ||
2fb8013c | 77 | #ifdef SCM_HAVE_T_INT64 |
e7a72986 MD |
78 | |
79 | unsigned long | |
92c2555f | 80 | scm_i_uniform32 (scm_t_i_rstate *state) |
e7a72986 | 81 | { |
2fb8013c RB |
82 | scm_t_int64 x = (scm_t_int64) A * state->w + state->c; |
83 | scm_t_int32 w = x & 0xffffffffUL; | |
e7a72986 MD |
84 | state->w = w; |
85 | state->c = x >> 32L; | |
86 | return w; | |
87 | } | |
88 | ||
89 | #else | |
90 | ||
91 | /* ww This is a portable version of the same RNG without 64 bit | |
92 | * * aa arithmetic. | |
93 | * ---- | |
94 | * xx It is only intended to provide identical behaviour on | |
95 | * xx platforms without 8 byte longs or long longs until | |
96 | * xx someone has implemented the routine in assembler code. | |
97 | * xxcc | |
98 | * ---- | |
99 | * ccww | |
100 | */ | |
101 | ||
102 | #define L(x) ((x) & 0xffff) | |
103 | #define H(x) ((x) >> 16) | |
104 | ||
105 | unsigned long | |
92c2555f | 106 | scm_i_uniform32 (scm_t_i_rstate *state) |
e7a72986 | 107 | { |
2fb8013c RB |
108 | scm_t_int32 x1 = L (A) * L (state->w); |
109 | scm_t_int32 x2 = L (A) * H (state->w); | |
110 | scm_t_int32 x3 = H (A) * L (state->w); | |
111 | scm_t_int32 w = L (x1) + L (state->c); | |
112 | scm_t_int32 m = H (x1) + L (x2) + L (x3) + H (state->c) + H (w); | |
113 | scm_t_int32 x4 = H (A) * H (state->w); | |
e7a72986 MD |
114 | state->w = w = (L (m) << 16) + L (w); |
115 | state->c = H (x2) + H (x3) + x4 + H (m); | |
116 | return w; | |
117 | } | |
118 | ||
119 | #endif | |
120 | ||
121 | void | |
92c2555f | 122 | scm_i_init_rstate (scm_t_i_rstate *state, char *seed, int n) |
e7a72986 | 123 | { |
2fb8013c RB |
124 | scm_t_int32 w = 0L; |
125 | scm_t_int32 c = 0L; | |
e7a72986 MD |
126 | int i, m; |
127 | for (i = 0; i < n; ++i) | |
128 | { | |
129 | m = i % 8; | |
130 | if (m < 4) | |
131 | w += seed[i] << (8 * m); | |
132 | else | |
133 | c += seed[i] << (8 * (m - 4)); | |
134 | } | |
135 | if ((w == 0 && c == 0) || (w == 0xffffffffUL && c == A - 1)) | |
136 | ++c; | |
137 | state->w = w; | |
138 | state->c = c; | |
139 | } | |
140 | ||
92c2555f MV |
141 | scm_t_i_rstate * |
142 | scm_i_copy_rstate (scm_t_i_rstate *state) | |
e7a72986 | 143 | { |
67329a9e | 144 | scm_t_rstate *new_state = scm_malloc (scm_the_rng.rstate_size); |
e7a72986 | 145 | if (new_state == 0) |
2500356c | 146 | scm_memory_error ("rstate"); |
e7a72986 MD |
147 | return memcpy (new_state, state, scm_the_rng.rstate_size); |
148 | } | |
149 | ||
150 | \f | |
151 | /* | |
152 | * Random number library functions | |
153 | */ | |
154 | ||
92c2555f | 155 | scm_t_rstate * |
9b741bb6 | 156 | scm_c_make_rstate (char *seed, int n) |
5ee11b7c | 157 | { |
67329a9e | 158 | scm_t_rstate *state = scm_malloc (scm_the_rng.rstate_size); |
5ee11b7c | 159 | if (state == 0) |
2500356c | 160 | scm_memory_error ("rstate"); |
5ee11b7c MD |
161 | state->reserved0 = 0; |
162 | scm_the_rng.init_rstate (state, seed, n); | |
163 | return state; | |
164 | } | |
165 | ||
2ade72d7 | 166 | |
92c2555f | 167 | scm_t_rstate * |
9b741bb6 | 168 | scm_c_default_rstate () |
2ade72d7 | 169 | #define FUNC_NAME "scm_c_default_rstate" |
9b741bb6 | 170 | { |
c5f268f8 | 171 | SCM state = SCM_VARIABLE_REF (scm_var_random_state); |
2ade72d7 DH |
172 | if (!SCM_RSTATEP (state)) |
173 | SCM_MISC_ERROR ("*random-state* contains bogus random state", SCM_EOL); | |
9b741bb6 MD |
174 | return SCM_RSTATE (state); |
175 | } | |
2ade72d7 DH |
176 | #undef FUNC_NAME |
177 | ||
9b741bb6 | 178 | |
e7a72986 | 179 | inline double |
92c2555f | 180 | scm_c_uniform01 (scm_t_rstate *state) |
e7a72986 | 181 | { |
5a92ddfd | 182 | double x = (double) scm_the_rng.random_bits (state) / (double) 0xffffffffUL; |
e7a72986 | 183 | return ((x + (double) scm_the_rng.random_bits (state)) |
5a92ddfd | 184 | / (double) 0xffffffffUL); |
e7a72986 MD |
185 | } |
186 | ||
187 | double | |
92c2555f | 188 | scm_c_normal01 (scm_t_rstate *state) |
e7a72986 MD |
189 | { |
190 | if (state->reserved0) | |
191 | { | |
192 | state->reserved0 = 0; | |
193 | return state->reserved1; | |
194 | } | |
195 | else | |
196 | { | |
197 | double r, a, n; | |
e7a72986 | 198 | |
9b741bb6 MD |
199 | r = sqrt (-2.0 * log (scm_c_uniform01 (state))); |
200 | a = 2.0 * M_PI * scm_c_uniform01 (state); | |
e7a72986 MD |
201 | |
202 | n = r * sin (a); | |
203 | state->reserved1 = r * cos (a); | |
5a92ddfd | 204 | state->reserved0 = 1; |
e7a72986 MD |
205 | |
206 | return n; | |
207 | } | |
208 | } | |
209 | ||
210 | double | |
92c2555f | 211 | scm_c_exp1 (scm_t_rstate *state) |
e7a72986 | 212 | { |
9b741bb6 | 213 | return - log (scm_c_uniform01 (state)); |
e7a72986 MD |
214 | } |
215 | ||
216 | unsigned char scm_masktab[256]; | |
217 | ||
218 | unsigned long | |
92c2555f | 219 | scm_c_random (scm_t_rstate *state, unsigned long m) |
e7a72986 MD |
220 | { |
221 | unsigned int r, mask; | |
222 | mask = (m < 0x100 | |
223 | ? scm_masktab[m] | |
224 | : (m < 0x10000 | |
5a92ddfd | 225 | ? scm_masktab[m >> 8] << 8 | 0xff |
e7a72986 | 226 | : (m < 0x1000000 |
5a92ddfd MD |
227 | ? scm_masktab[m >> 16] << 16 | 0xffff |
228 | : scm_masktab[m >> 24] << 24 | 0xffffff))); | |
e7a72986 MD |
229 | while ((r = scm_the_rng.random_bits (state) & mask) >= m); |
230 | return r; | |
231 | } | |
232 | ||
ea7f6344 MD |
233 | /* |
234 | SCM scm_c_random_bignum (scm_t_rstate *state, SCM m) | |
235 | ||
236 | Takes a random state (source of random bits) and a bignum m. | |
237 | Returns a bignum b, 0 <= b < m. | |
238 | ||
239 | It does this by allocating a bignum b with as many base 65536 digits | |
240 | as m, filling b with random bits (in 32 bit chunks) up to the most | |
241 | significant 1 in m, and, finally checking if the resultant b is too | |
242 | large (>= m). If too large, we simply repeat the process again. (It | |
243 | is important to throw away all generated random bits if b >= m, | |
244 | otherwise we'll end up with a distorted distribution.) | |
245 | ||
246 | */ | |
247 | ||
e7a72986 | 248 | SCM |
92c2555f | 249 | scm_c_random_bignum (scm_t_rstate *state, SCM m) |
e7a72986 | 250 | { |
969d3bd0 | 251 | SCM result = scm_i_mkbig (); |
969d3bd0 RB |
252 | const size_t m_bits = mpz_sizeinbase (SCM_I_BIG_MPZ (m), 2); |
253 | /* how many bits would only partially fill the last unsigned long? */ | |
0d79003d | 254 | const size_t end_bits = m_bits % (sizeof (unsigned long) * SCM_CHAR_BIT); |
969d3bd0 | 255 | unsigned long *random_chunks = NULL; |
0d79003d RB |
256 | const unsigned long num_full_chunks = |
257 | m_bits / (sizeof (unsigned long) * SCM_CHAR_BIT); | |
969d3bd0 RB |
258 | const unsigned long num_chunks = num_full_chunks + ((end_bits) ? 1 : 0); |
259 | ||
260 | /* we know the result will be this big */ | |
261 | mpz_realloc2 (SCM_I_BIG_MPZ (result), m_bits); | |
262 | ||
263 | random_chunks = | |
264 | (unsigned long *) scm_gc_calloc (num_chunks * sizeof (unsigned long), | |
265 | "random bignum chunks"); | |
266 | ||
372691d8 | 267 | do |
e7a72986 | 268 | { |
969d3bd0 RB |
269 | unsigned long *current_chunk = random_chunks + (num_chunks - 1); |
270 | unsigned long chunks_left = num_chunks; | |
969d3bd0 RB |
271 | |
272 | mpz_set_ui (SCM_I_BIG_MPZ (result), 0); | |
273 | ||
274 | if (end_bits) | |
275 | { | |
276 | /* generate a mask with ones in the end_bits position, i.e. if | |
277 | end_bits is 3, then we'd have a mask of ...0000000111 */ | |
278 | const unsigned long rndbits = scm_the_rng.random_bits (state); | |
0d79003d | 279 | int rshift = (sizeof (unsigned long) * SCM_CHAR_BIT) - end_bits; |
969d3bd0 RB |
280 | unsigned long mask = ((unsigned long) ULONG_MAX) >> rshift; |
281 | unsigned long highest_bits = rndbits & mask; | |
282 | *current_chunk-- = highest_bits; | |
283 | chunks_left--; | |
284 | } | |
285 | ||
286 | while (chunks_left) | |
287 | { | |
288 | /* now fill in the remaining unsigned long sized chunks */ | |
289 | *current_chunk-- = scm_the_rng.random_bits (state); | |
290 | chunks_left--; | |
291 | } | |
292 | mpz_import (SCM_I_BIG_MPZ (result), | |
293 | num_chunks, | |
294 | -1, | |
295 | sizeof (unsigned long), | |
296 | 0, | |
297 | 0, | |
298 | random_chunks); | |
372691d8 MD |
299 | /* if result >= m, regenerate it (it is important to regenerate |
300 | all bits in order not to get a distorted distribution) */ | |
301 | } while (mpz_cmp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (m)) >= 0); | |
969d3bd0 RB |
302 | scm_gc_free (random_chunks, |
303 | num_chunks * sizeof (unsigned long), | |
304 | "random bignum chunks"); | |
8db9cc6c | 305 | return scm_i_normbig (result); |
e7a72986 MD |
306 | } |
307 | ||
308 | /* | |
309 | * Scheme level representation of random states. | |
310 | */ | |
311 | ||
92c2555f | 312 | scm_t_bits scm_tc16_rstate; |
e7a72986 MD |
313 | |
314 | static SCM | |
92c2555f | 315 | make_rstate (scm_t_rstate *state) |
e7a72986 | 316 | { |
23a62151 | 317 | SCM_RETURN_NEWSMOB (scm_tc16_rstate, state); |
e7a72986 MD |
318 | } |
319 | ||
1be6b49c | 320 | static size_t |
e841c3e0 | 321 | rstate_free (SCM rstate) |
e7a72986 MD |
322 | { |
323 | free (SCM_RSTATE (rstate)); | |
dfd71aba | 324 | return 0; |
e7a72986 MD |
325 | } |
326 | ||
e7a72986 MD |
327 | /* |
328 | * Scheme level interface. | |
329 | */ | |
330 | ||
86d31dfe | 331 | SCM_GLOBAL_VARIABLE_INIT (scm_var_random_state, "*random-state*", scm_seed_to_random_state (scm_makfrom0str ("URL:http://stat.fsu.edu/~geo/diehard.html"))); |
e7a72986 | 332 | |
a1ec6916 | 333 | SCM_DEFINE (scm_random, "random", 1, 1, 0, |
1bbd0b84 | 334 | (SCM n, SCM state), |
34d19ef6 | 335 | "Return a number in [0, N).\n" |
d928e0b4 | 336 | "\n" |
9401323e NJ |
337 | "Accepts a positive integer or real n and returns a\n" |
338 | "number of the same type between zero (inclusive) and\n" | |
339 | "N (exclusive). The values returned have a uniform\n" | |
d928e0b4 GB |
340 | "distribution.\n" |
341 | "\n" | |
3b644514 MG |
342 | "The optional argument @var{state} must be of the type produced\n" |
343 | "by @code{seed->random-state}. It defaults to the value of the\n" | |
344 | "variable @var{*random-state*}. This object is used to maintain\n" | |
345 | "the state of the pseudo-random-number generator and is altered\n" | |
346 | "as a side effect of the random operation.") | |
1bbd0b84 | 347 | #define FUNC_NAME s_scm_random |
e7a72986 MD |
348 | { |
349 | if (SCM_UNBNDP (state)) | |
86d31dfe | 350 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 351 | SCM_VALIDATE_RSTATE (2, state); |
e7a72986 MD |
352 | if (SCM_INUMP (n)) |
353 | { | |
354 | unsigned long m = SCM_INUM (n); | |
34d19ef6 | 355 | SCM_ASSERT_RANGE (1, n, m > 0); |
93ccaef0 | 356 | return SCM_I_MAKINUM (scm_c_random (SCM_RSTATE (state), m)); |
e7a72986 | 357 | } |
34d19ef6 | 358 | SCM_VALIDATE_NIM (1, n); |
e7a72986 | 359 | if (SCM_REALP (n)) |
f8de44c1 DH |
360 | return scm_make_real (SCM_REAL_VALUE (n) |
361 | * scm_c_uniform01 (SCM_RSTATE (state))); | |
969d3bd0 RB |
362 | |
363 | SCM_VALIDATE_BIGINT (1, n); | |
9b741bb6 | 364 | return scm_c_random_bignum (SCM_RSTATE (state), n); |
e7a72986 | 365 | } |
1bbd0b84 | 366 | #undef FUNC_NAME |
e7a72986 | 367 | |
a1ec6916 | 368 | SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0, |
1bbd0b84 | 369 | (SCM state), |
3b644514 | 370 | "Return a copy of the random state @var{state}.") |
1bbd0b84 | 371 | #define FUNC_NAME s_scm_copy_random_state |
e7a72986 MD |
372 | { |
373 | if (SCM_UNBNDP (state)) | |
86d31dfe | 374 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 375 | SCM_VALIDATE_RSTATE (1, state); |
e7a72986 MD |
376 | return make_rstate (scm_the_rng.copy_rstate (SCM_RSTATE (state))); |
377 | } | |
1bbd0b84 | 378 | #undef FUNC_NAME |
e7a72986 | 379 | |
a1ec6916 | 380 | SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0, |
1bbd0b84 | 381 | (SCM seed), |
3b644514 | 382 | "Return a new random state using @var{seed}.") |
1bbd0b84 | 383 | #define FUNC_NAME s_scm_seed_to_random_state |
5ee11b7c MD |
384 | { |
385 | if (SCM_NUMBERP (seed)) | |
386 | seed = scm_number_to_string (seed, SCM_UNDEFINED); | |
34d19ef6 | 387 | SCM_VALIDATE_STRING (1, seed); |
34f0f2b8 | 388 | return make_rstate (scm_c_make_rstate (SCM_STRING_CHARS (seed), |
b7ead2ae | 389 | SCM_STRING_LENGTH (seed))); |
5ee11b7c | 390 | } |
1bbd0b84 | 391 | #undef FUNC_NAME |
5ee11b7c | 392 | |
a1ec6916 | 393 | SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0, |
1bbd0b84 | 394 | (SCM state), |
1e6808ea MG |
395 | "Return a uniformly distributed inexact real random number in\n" |
396 | "[0,1).") | |
1bbd0b84 | 397 | #define FUNC_NAME s_scm_random_uniform |
e7a72986 MD |
398 | { |
399 | if (SCM_UNBNDP (state)) | |
86d31dfe | 400 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 401 | SCM_VALIDATE_RSTATE (1, state); |
f8de44c1 | 402 | return scm_make_real (scm_c_uniform01 (SCM_RSTATE (state))); |
e7a72986 | 403 | } |
1bbd0b84 | 404 | #undef FUNC_NAME |
e7a72986 | 405 | |
a1ec6916 | 406 | SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0, |
1bbd0b84 | 407 | (SCM state), |
1e6808ea MG |
408 | "Return an inexact real in a normal distribution. The\n" |
409 | "distribution used has mean 0 and standard deviation 1. For a\n" | |
410 | "normal distribution with mean m and standard deviation d use\n" | |
411 | "@code{(+ m (* d (random:normal)))}.") | |
1bbd0b84 | 412 | #define FUNC_NAME s_scm_random_normal |
afe5177e GH |
413 | { |
414 | if (SCM_UNBNDP (state)) | |
86d31dfe | 415 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 416 | SCM_VALIDATE_RSTATE (1, state); |
f8de44c1 | 417 | return scm_make_real (scm_c_normal01 (SCM_RSTATE (state))); |
afe5177e | 418 | } |
1bbd0b84 | 419 | #undef FUNC_NAME |
afe5177e | 420 | |
c7ad1f89 | 421 | #if SCM_HAVE_ARRAYS |
afe5177e | 422 | |
e7a72986 MD |
423 | static void |
424 | vector_scale (SCM v, double c) | |
425 | { | |
b7ead2ae | 426 | int n = SCM_INUM (scm_uniform_vector_length (v)); |
e7a72986 MD |
427 | if (SCM_VECTORP (v)) |
428 | while (--n >= 0) | |
eb42e2f0 | 429 | SCM_REAL_VALUE (SCM_VELTS (v)[n]) *= c; |
e7a72986 MD |
430 | else |
431 | while (--n >= 0) | |
432 | ((double *) SCM_VELTS (v))[n] *= c; | |
433 | } | |
434 | ||
435 | static double | |
436 | vector_sum_squares (SCM v) | |
437 | { | |
438 | double x, sum = 0.0; | |
b7ead2ae | 439 | int n = SCM_INUM (scm_uniform_vector_length (v)); |
e7a72986 MD |
440 | if (SCM_VECTORP (v)) |
441 | while (--n >= 0) | |
442 | { | |
eb42e2f0 | 443 | x = SCM_REAL_VALUE (SCM_VELTS (v)[n]); |
e7a72986 MD |
444 | sum += x * x; |
445 | } | |
446 | else | |
447 | while (--n >= 0) | |
448 | { | |
449 | x = ((double *) SCM_VELTS (v))[n]; | |
450 | sum += x * x; | |
451 | } | |
452 | return sum; | |
453 | } | |
454 | ||
e7a72986 MD |
455 | /* For the uniform distribution on the solid sphere, note that in |
456 | * this distribution the length r of the vector has cumulative | |
457 | * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be | |
458 | * generated as r=u^(1/n). | |
459 | */ | |
a1ec6916 | 460 | SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0, |
1bbd0b84 | 461 | (SCM v, SCM state), |
d928e0b4 GB |
462 | "Fills vect with inexact real random numbers\n" |
463 | "the sum of whose squares is less than 1.0.\n" | |
9401323e NJ |
464 | "Thinking of vect as coordinates in space of\n" |
465 | "dimension n = (vector-length vect), the coordinates\n" | |
72dd0a03 | 466 | "are uniformly distributed within the unit n-sphere.\n" |
64ba8e85 | 467 | "The sum of the squares of the numbers is returned.") |
1bbd0b84 | 468 | #define FUNC_NAME s_scm_random_solid_sphere_x |
e7a72986 | 469 | { |
34d19ef6 | 470 | SCM_VALIDATE_VECTOR_OR_DVECTOR (1, v); |
e7a72986 | 471 | if (SCM_UNBNDP (state)) |
86d31dfe | 472 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 473 | SCM_VALIDATE_RSTATE (2, state); |
e7a72986 MD |
474 | scm_random_normal_vector_x (v, state); |
475 | vector_scale (v, | |
9b741bb6 | 476 | pow (scm_c_uniform01 (SCM_RSTATE (state)), |
b7ead2ae | 477 | 1.0 / SCM_INUM (scm_uniform_vector_length (v))) |
e7a72986 MD |
478 | / sqrt (vector_sum_squares (v))); |
479 | return SCM_UNSPECIFIED; | |
480 | } | |
1bbd0b84 | 481 | #undef FUNC_NAME |
e7a72986 | 482 | |
a1ec6916 | 483 | SCM_DEFINE (scm_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0, |
1bbd0b84 | 484 | (SCM v, SCM state), |
d928e0b4 GB |
485 | "Fills vect with inexact real random numbers\n" |
486 | "the sum of whose squares is equal to 1.0.\n" | |
9401323e | 487 | "Thinking of vect as coordinates in space of\n" |
d928e0b4 | 488 | "dimension n = (vector-length vect), the coordinates\n" |
9401323e | 489 | "are uniformly distributed over the surface of the\n" |
72dd0a03 | 490 | "unit n-sphere.") |
1bbd0b84 | 491 | #define FUNC_NAME s_scm_random_hollow_sphere_x |
e7a72986 | 492 | { |
34d19ef6 | 493 | SCM_VALIDATE_VECTOR_OR_DVECTOR (1, v); |
e7a72986 | 494 | if (SCM_UNBNDP (state)) |
86d31dfe | 495 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 496 | SCM_VALIDATE_RSTATE (2, state); |
e7a72986 MD |
497 | scm_random_normal_vector_x (v, state); |
498 | vector_scale (v, 1 / sqrt (vector_sum_squares (v))); | |
499 | return SCM_UNSPECIFIED; | |
500 | } | |
1bbd0b84 | 501 | #undef FUNC_NAME |
e7a72986 | 502 | |
1bbd0b84 | 503 | |
a1ec6916 | 504 | SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0, |
1bbd0b84 | 505 | (SCM v, SCM state), |
d928e0b4 GB |
506 | "Fills vect with inexact real random numbers that are\n" |
507 | "independent and standard normally distributed\n" | |
64ba8e85 | 508 | "(i.e., with mean 0 and variance 1).") |
1bbd0b84 | 509 | #define FUNC_NAME s_scm_random_normal_vector_x |
e7a72986 MD |
510 | { |
511 | int n; | |
34d19ef6 | 512 | SCM_VALIDATE_VECTOR_OR_DVECTOR (1, v); |
e7a72986 | 513 | if (SCM_UNBNDP (state)) |
86d31dfe | 514 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 515 | SCM_VALIDATE_RSTATE (2, state); |
b7ead2ae | 516 | n = SCM_INUM (scm_uniform_vector_length (v)); |
e7a72986 MD |
517 | if (SCM_VECTORP (v)) |
518 | while (--n >= 0) | |
34d19ef6 | 519 | SCM_VECTOR_SET (v, n, scm_make_real (scm_c_normal01 (SCM_RSTATE (state)))); |
e7a72986 MD |
520 | else |
521 | while (--n >= 0) | |
9b741bb6 | 522 | ((double *) SCM_VELTS (v))[n] = scm_c_normal01 (SCM_RSTATE (state)); |
e7a72986 MD |
523 | return SCM_UNSPECIFIED; |
524 | } | |
1bbd0b84 | 525 | #undef FUNC_NAME |
e7a72986 | 526 | |
2fb8013c | 527 | #endif /* SCM_HAVE_ARRAYS */ |
afe5177e | 528 | |
a1ec6916 | 529 | SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0, |
1bbd0b84 | 530 | (SCM state), |
1e6808ea MG |
531 | "Return an inexact real in an exponential distribution with mean\n" |
532 | "1. For an exponential distribution with mean u use (* u\n" | |
533 | "(random:exp)).") | |
1bbd0b84 | 534 | #define FUNC_NAME s_scm_random_exp |
e7a72986 MD |
535 | { |
536 | if (SCM_UNBNDP (state)) | |
86d31dfe | 537 | state = SCM_VARIABLE_REF (scm_var_random_state); |
34d19ef6 | 538 | SCM_VALIDATE_RSTATE (1, state); |
f8de44c1 | 539 | return scm_make_real (scm_c_exp1 (SCM_RSTATE (state))); |
e7a72986 | 540 | } |
1bbd0b84 | 541 | #undef FUNC_NAME |
e7a72986 MD |
542 | |
543 | void | |
544 | scm_init_random () | |
545 | { | |
546 | int i, m; | |
547 | /* plug in default RNG */ | |
92c2555f | 548 | scm_t_rng rng = |
e7a72986 | 549 | { |
92c2555f | 550 | sizeof (scm_t_i_rstate), |
e7a72986 MD |
551 | (unsigned long (*)()) scm_i_uniform32, |
552 | (void (*)()) scm_i_init_rstate, | |
92c2555f | 553 | (scm_t_rstate *(*)()) scm_i_copy_rstate |
e7a72986 MD |
554 | }; |
555 | scm_the_rng = rng; | |
556 | ||
e841c3e0 KN |
557 | scm_tc16_rstate = scm_make_smob_type ("random-state", 0); |
558 | scm_set_smob_free (scm_tc16_rstate, rstate_free); | |
e7a72986 MD |
559 | |
560 | for (m = 1; m <= 0x100; m <<= 1) | |
561 | for (i = m >> 1; i < m; ++i) | |
562 | scm_masktab[i] = m - 1; | |
563 | ||
a0599745 | 564 | #include "libguile/random.x" |
e7a72986 MD |
565 | |
566 | scm_add_feature ("random"); | |
567 | } | |
89e00824 ML |
568 | |
569 | /* | |
570 | Local Variables: | |
571 | c-file-style: "gnu" | |
572 | End: | |
573 | */ |