Commit | Line | Data |
---|---|---|
2ade72d7 | 1 | /* Copyright (C) 1999,2000,2001 Free Software Foundation, Inc. |
e7a72986 MD |
2 | * This program is free software; you can redistribute it and/or modify |
3 | * it under the terms of the GNU General Public License as published by | |
4 | * the Free Software Foundation; either version 2, or (at your option) | |
5 | * any later version. | |
6 | * | |
7 | * This program is distributed in the hope that it will be useful, | |
8 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
9 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
10 | * GNU General Public License for more details. | |
11 | * | |
12 | * You should have received a copy of the GNU General Public License | |
13 | * along with this software; see the file COPYING. If not, write to | |
14 | * the Free Software Foundation, Inc., 59 Temple Place, Suite 330, | |
15 | * Boston, MA 02111-1307 USA | |
16 | * | |
17 | * As a special exception, the Free Software Foundation gives permission | |
18 | * for additional uses of the text contained in its release of GUILE. | |
19 | * | |
20 | * The exception is that, if you link the GUILE library with other files | |
21 | * to produce an executable, this does not by itself cause the | |
22 | * resulting executable to be covered by the GNU General Public License. | |
23 | * Your use of that executable is in no way restricted on account of | |
24 | * linking the GUILE library code into it. | |
25 | * | |
26 | * This exception does not however invalidate any other reasons why | |
27 | * the executable file might be covered by the GNU General Public License. | |
28 | * | |
29 | * This exception applies only to the code released by the | |
30 | * Free Software Foundation under the name GUILE. If you copy | |
31 | * code from other Free Software Foundation releases into a copy of | |
32 | * GUILE, as the General Public License permits, the exception does | |
33 | * not apply to the code that you add in this way. To avoid misleading | |
34 | * anyone as to the status of such modified files, you must delete | |
35 | * this exception notice from them. | |
36 | * | |
37 | * If you write modifications of your own for GUILE, it is your choice | |
38 | * whether to permit this exception to apply to your modifications. | |
39 | * If you do not wish that, delete this exception notice. */ | |
40 | ||
1bbd0b84 GB |
41 | |
42 | ||
e7a72986 MD |
43 | /* Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */ |
44 | ||
a0599745 | 45 | #include "libguile/_scm.h" |
e7a72986 | 46 | |
7f146094 | 47 | #include <stdio.h> |
e7a72986 | 48 | #include <math.h> |
f34d19c7 | 49 | #include <string.h> |
a0599745 MD |
50 | #include "libguile/smob.h" |
51 | #include "libguile/numbers.h" | |
52 | #include "libguile/feature.h" | |
53 | #include "libguile/strings.h" | |
e9bfab50 | 54 | #include "libguile/unif.h" |
a0599745 | 55 | #include "libguile/vectors.h" |
e7a72986 | 56 | |
a0599745 MD |
57 | #include "libguile/validate.h" |
58 | #include "libguile/random.h" | |
e7a72986 MD |
59 | |
60 | \f | |
61 | /* | |
62 | * A plugin interface for RNGs | |
63 | * | |
64 | * Using this interface, it is possible for the application to tell | |
65 | * libguile to use a different RNG. This is desirable if it is | |
66 | * necessary to use the same RNG everywhere in the application in | |
67 | * order to prevent interference, if the application uses RNG | |
68 | * hardware, or if the application has special demands on the RNG. | |
69 | * | |
70 | * Look in random.h and how the default generator is "plugged in" in | |
71 | * scm_init_random(). | |
72 | */ | |
73 | ||
92c2555f | 74 | scm_t_rng scm_the_rng; |
e7a72986 MD |
75 | |
76 | \f | |
77 | /* | |
78 | * The prepackaged RNG | |
79 | * | |
80 | * This is the MWC (Multiply With Carry) random number generator | |
81 | * described by George Marsaglia at the Department of Statistics and | |
82 | * Supercomputer Computations Research Institute, The Florida State | |
83 | * University (http://stat.fsu.edu/~geo). | |
84 | * | |
85 | * It uses 64 bits, has a period of 4578426017172946943 (4.6e18), and | |
86 | * passes all tests in the DIEHARD test suite | |
87 | * (http://stat.fsu.edu/~geo/diehard.html) | |
88 | */ | |
89 | ||
90 | #define A 2131995753UL | |
91 | ||
82893676 MG |
92 | #ifndef M_PI |
93 | #define M_PI 3.14159265359 | |
94 | #endif | |
95 | ||
e7a72986 MD |
96 | #if SIZEOF_LONG > 4 |
97 | #if SIZEOF_INT > 4 | |
98 | #define LONG32 unsigned short | |
99 | #else | |
100 | #define LONG32 unsigned int | |
101 | #endif | |
102 | #define LONG64 unsigned long | |
103 | #else | |
104 | #define LONG32 unsigned long | |
82893676 MG |
105 | #ifdef __MINGW32__ |
106 | #define LONG64 unsigned __int64 | |
107 | #else | |
e7a72986 MD |
108 | #define LONG64 unsigned long long |
109 | #endif | |
82893676 | 110 | #endif |
e7a72986 MD |
111 | |
112 | #if SIZEOF_LONG > 4 || defined (HAVE_LONG_LONGS) | |
113 | ||
114 | unsigned long | |
92c2555f | 115 | scm_i_uniform32 (scm_t_i_rstate *state) |
e7a72986 MD |
116 | { |
117 | LONG64 x = (LONG64) A * state->w + state->c; | |
118 | LONG32 w = x & 0xffffffffUL; | |
119 | state->w = w; | |
120 | state->c = x >> 32L; | |
121 | return w; | |
122 | } | |
123 | ||
124 | #else | |
125 | ||
126 | /* ww This is a portable version of the same RNG without 64 bit | |
127 | * * aa arithmetic. | |
128 | * ---- | |
129 | * xx It is only intended to provide identical behaviour on | |
130 | * xx platforms without 8 byte longs or long longs until | |
131 | * xx someone has implemented the routine in assembler code. | |
132 | * xxcc | |
133 | * ---- | |
134 | * ccww | |
135 | */ | |
136 | ||
137 | #define L(x) ((x) & 0xffff) | |
138 | #define H(x) ((x) >> 16) | |
139 | ||
140 | unsigned long | |
92c2555f | 141 | scm_i_uniform32 (scm_t_i_rstate *state) |
e7a72986 MD |
142 | { |
143 | LONG32 x1 = L (A) * L (state->w); | |
144 | LONG32 x2 = L (A) * H (state->w); | |
145 | LONG32 x3 = H (A) * L (state->w); | |
146 | LONG32 w = L (x1) + L (state->c); | |
147 | LONG32 m = H (x1) + L (x2) + L (x3) + H (state->c) + H (w); | |
148 | LONG32 x4 = H (A) * H (state->w); | |
149 | state->w = w = (L (m) << 16) + L (w); | |
150 | state->c = H (x2) + H (x3) + x4 + H (m); | |
151 | return w; | |
152 | } | |
153 | ||
154 | #endif | |
155 | ||
156 | void | |
92c2555f | 157 | scm_i_init_rstate (scm_t_i_rstate *state, char *seed, int n) |
e7a72986 MD |
158 | { |
159 | LONG32 w = 0L; | |
160 | LONG32 c = 0L; | |
161 | int i, m; | |
162 | for (i = 0; i < n; ++i) | |
163 | { | |
164 | m = i % 8; | |
165 | if (m < 4) | |
166 | w += seed[i] << (8 * m); | |
167 | else | |
168 | c += seed[i] << (8 * (m - 4)); | |
169 | } | |
170 | if ((w == 0 && c == 0) || (w == 0xffffffffUL && c == A - 1)) | |
171 | ++c; | |
172 | state->w = w; | |
173 | state->c = c; | |
174 | } | |
175 | ||
92c2555f MV |
176 | scm_t_i_rstate * |
177 | scm_i_copy_rstate (scm_t_i_rstate *state) | |
e7a72986 | 178 | { |
92c2555f | 179 | scm_t_rstate *new_state = malloc (scm_the_rng.rstate_size); |
e7a72986 | 180 | if (new_state == 0) |
2500356c | 181 | scm_memory_error ("rstate"); |
e7a72986 MD |
182 | return memcpy (new_state, state, scm_the_rng.rstate_size); |
183 | } | |
184 | ||
185 | \f | |
186 | /* | |
187 | * Random number library functions | |
188 | */ | |
189 | ||
92c2555f | 190 | scm_t_rstate * |
9b741bb6 | 191 | scm_c_make_rstate (char *seed, int n) |
5ee11b7c | 192 | { |
92c2555f | 193 | scm_t_rstate *state = malloc (scm_the_rng.rstate_size); |
5ee11b7c | 194 | if (state == 0) |
2500356c | 195 | scm_memory_error ("rstate"); |
5ee11b7c MD |
196 | state->reserved0 = 0; |
197 | scm_the_rng.init_rstate (state, seed, n); | |
198 | return state; | |
199 | } | |
200 | ||
2ade72d7 | 201 | |
92c2555f | 202 | scm_t_rstate * |
9b741bb6 | 203 | scm_c_default_rstate () |
2ade72d7 | 204 | #define FUNC_NAME "scm_c_default_rstate" |
9b741bb6 MD |
205 | { |
206 | SCM state = SCM_CDR (scm_var_random_state); | |
2ade72d7 DH |
207 | if (!SCM_RSTATEP (state)) |
208 | SCM_MISC_ERROR ("*random-state* contains bogus random state", SCM_EOL); | |
9b741bb6 MD |
209 | return SCM_RSTATE (state); |
210 | } | |
2ade72d7 DH |
211 | #undef FUNC_NAME |
212 | ||
9b741bb6 | 213 | |
e7a72986 | 214 | inline double |
92c2555f | 215 | scm_c_uniform01 (scm_t_rstate *state) |
e7a72986 | 216 | { |
5a92ddfd | 217 | double x = (double) scm_the_rng.random_bits (state) / (double) 0xffffffffUL; |
e7a72986 | 218 | return ((x + (double) scm_the_rng.random_bits (state)) |
5a92ddfd | 219 | / (double) 0xffffffffUL); |
e7a72986 MD |
220 | } |
221 | ||
222 | double | |
92c2555f | 223 | scm_c_normal01 (scm_t_rstate *state) |
e7a72986 MD |
224 | { |
225 | if (state->reserved0) | |
226 | { | |
227 | state->reserved0 = 0; | |
228 | return state->reserved1; | |
229 | } | |
230 | else | |
231 | { | |
232 | double r, a, n; | |
e7a72986 | 233 | |
9b741bb6 MD |
234 | r = sqrt (-2.0 * log (scm_c_uniform01 (state))); |
235 | a = 2.0 * M_PI * scm_c_uniform01 (state); | |
e7a72986 MD |
236 | |
237 | n = r * sin (a); | |
238 | state->reserved1 = r * cos (a); | |
5a92ddfd | 239 | state->reserved0 = 1; |
e7a72986 MD |
240 | |
241 | return n; | |
242 | } | |
243 | } | |
244 | ||
245 | double | |
92c2555f | 246 | scm_c_exp1 (scm_t_rstate *state) |
e7a72986 | 247 | { |
9b741bb6 | 248 | return - log (scm_c_uniform01 (state)); |
e7a72986 MD |
249 | } |
250 | ||
251 | unsigned char scm_masktab[256]; | |
252 | ||
253 | unsigned long | |
92c2555f | 254 | scm_c_random (scm_t_rstate *state, unsigned long m) |
e7a72986 MD |
255 | { |
256 | unsigned int r, mask; | |
257 | mask = (m < 0x100 | |
258 | ? scm_masktab[m] | |
259 | : (m < 0x10000 | |
5a92ddfd | 260 | ? scm_masktab[m >> 8] << 8 | 0xff |
e7a72986 | 261 | : (m < 0x1000000 |
5a92ddfd MD |
262 | ? scm_masktab[m >> 16] << 16 | 0xffff |
263 | : scm_masktab[m >> 24] << 24 | 0xffffff))); | |
e7a72986 MD |
264 | while ((r = scm_the_rng.random_bits (state) & mask) >= m); |
265 | return r; | |
266 | } | |
267 | ||
268 | SCM | |
92c2555f | 269 | scm_c_random_bignum (scm_t_rstate *state, SCM m) |
e7a72986 MD |
270 | { |
271 | SCM b; | |
a7e7ea3e MD |
272 | int i, nd; |
273 | LONG32 *bits, mask, w; | |
274 | nd = SCM_NUMDIGS (m); | |
275 | /* calculate mask for most significant digit */ | |
e7a72986 MD |
276 | #if SIZEOF_INT == 4 |
277 | /* 16 bit digits */ | |
a7e7ea3e | 278 | if (nd & 1) |
e7a72986 MD |
279 | { |
280 | /* fix most significant 16 bits */ | |
a7e7ea3e | 281 | unsigned short s = SCM_BDIGITS (m)[nd - 1]; |
5a92ddfd | 282 | mask = s < 0x100 ? scm_masktab[s] : scm_masktab[s >> 8] << 8 | 0xff; |
e7a72986 MD |
283 | } |
284 | else | |
285 | #endif | |
286 | { | |
287 | /* fix most significant 32 bits */ | |
96e263d6 | 288 | #if SIZEOF_INT == 4 |
2a0279c9 MD |
289 | w = SCM_BDIGITS (m)[nd - 1] << 16 | SCM_BDIGITS (m)[nd - 2]; |
290 | #else | |
96e263d6 | 291 | w = SCM_BDIGITS (m)[nd - 1]; |
2a0279c9 | 292 | #endif |
e7a72986 MD |
293 | mask = (w < 0x10000 |
294 | ? (w < 0x100 | |
295 | ? scm_masktab[w] | |
5a92ddfd | 296 | : scm_masktab[w >> 8] << 8 | 0xff) |
e7a72986 | 297 | : (w < 0x1000000 |
5a92ddfd MD |
298 | ? scm_masktab[w >> 16] << 16 | 0xffff |
299 | : scm_masktab[w >> 24] << 24 | 0xffffff)); | |
e7a72986 | 300 | } |
1be6b49c | 301 | b = scm_i_mkbig (nd, 0); |
a7e7ea3e MD |
302 | bits = (LONG32 *) SCM_BDIGITS (b); |
303 | do | |
e7a72986 | 304 | { |
a7e7ea3e MD |
305 | i = nd; |
306 | /* treat most significant digit specially */ | |
307 | #if SIZEOF_INT == 4 | |
308 | /* 16 bit digits */ | |
309 | if (i & 1) | |
310 | { | |
311 | ((SCM_BIGDIG*) bits)[i - 1] = scm_the_rng.random_bits (state) & mask; | |
312 | i /= 2; | |
313 | } | |
314 | else | |
315 | #endif | |
316 | { | |
317 | /* fix most significant 32 bits */ | |
96e263d6 | 318 | #if SIZEOF_INT == 4 |
2a0279c9 | 319 | w = scm_the_rng.random_bits (state) & mask; |
96e263d6 MD |
320 | ((SCM_BIGDIG*) bits)[i - 2] = w & 0xffff; |
321 | ((SCM_BIGDIG*) bits)[i - 1] = w >> 16; | |
322 | i = i / 2 - 1; | |
2a0279c9 | 323 | #else |
96e263d6 | 324 | i /= 2; |
a7e7ea3e | 325 | bits[--i] = scm_the_rng.random_bits (state) & mask; |
2a0279c9 | 326 | #endif |
a7e7ea3e MD |
327 | } |
328 | /* now fill up the rest of the bignum */ | |
329 | while (i) | |
330 | bits[--i] = scm_the_rng.random_bits (state); | |
1be6b49c | 331 | b = scm_i_normbig (b); |
a7e7ea3e MD |
332 | if (SCM_INUMP (b)) |
333 | return b; | |
334 | } while (scm_bigcomp (b, m) <= 0); | |
335 | return b; | |
e7a72986 MD |
336 | } |
337 | ||
338 | /* | |
339 | * Scheme level representation of random states. | |
340 | */ | |
341 | ||
92c2555f | 342 | scm_t_bits scm_tc16_rstate; |
e7a72986 MD |
343 | |
344 | static SCM | |
92c2555f | 345 | make_rstate (scm_t_rstate *state) |
e7a72986 | 346 | { |
23a62151 | 347 | SCM_RETURN_NEWSMOB (scm_tc16_rstate, state); |
e7a72986 MD |
348 | } |
349 | ||
1be6b49c | 350 | static size_t |
e841c3e0 | 351 | rstate_free (SCM rstate) |
e7a72986 MD |
352 | { |
353 | free (SCM_RSTATE (rstate)); | |
354 | return scm_the_rng.rstate_size; | |
355 | } | |
356 | ||
e7a72986 MD |
357 | /* |
358 | * Scheme level interface. | |
359 | */ | |
360 | ||
86d31dfe | 361 | 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 | 362 | |
a1ec6916 | 363 | SCM_DEFINE (scm_random, "random", 1, 1, 0, |
1bbd0b84 | 364 | (SCM n, SCM state), |
d928e0b4 GB |
365 | "Return a number in [0,N).\n" |
366 | "\n" | |
367 | "Accepts a positive integer or real n and returns a \n" | |
368 | "number of the same type between zero (inclusive) and \n" | |
369 | "N (exclusive). The values returned have a uniform \n" | |
370 | "distribution.\n" | |
371 | "\n" | |
3b644514 MG |
372 | "The optional argument @var{state} must be of the type produced\n" |
373 | "by @code{seed->random-state}. It defaults to the value of the\n" | |
374 | "variable @var{*random-state*}. This object is used to maintain\n" | |
375 | "the state of the pseudo-random-number generator and is altered\n" | |
376 | "as a side effect of the random operation.") | |
1bbd0b84 | 377 | #define FUNC_NAME s_scm_random |
e7a72986 MD |
378 | { |
379 | if (SCM_UNBNDP (state)) | |
86d31dfe | 380 | state = SCM_VARIABLE_REF (scm_var_random_state); |
3b3b36dd | 381 | SCM_VALIDATE_RSTATE (2,state); |
e7a72986 MD |
382 | if (SCM_INUMP (n)) |
383 | { | |
384 | unsigned long m = SCM_INUM (n); | |
1bbd0b84 | 385 | SCM_ASSERT_RANGE (1,n,m > 0); |
9b741bb6 | 386 | return SCM_MAKINUM (scm_c_random (SCM_RSTATE (state), m)); |
e7a72986 | 387 | } |
6b5a304f | 388 | SCM_VALIDATE_NIM (1,n); |
e7a72986 | 389 | if (SCM_REALP (n)) |
f8de44c1 DH |
390 | return scm_make_real (SCM_REAL_VALUE (n) |
391 | * scm_c_uniform01 (SCM_RSTATE (state))); | |
950cc72b | 392 | SCM_VALIDATE_SMOB (1, n, big); |
9b741bb6 | 393 | return scm_c_random_bignum (SCM_RSTATE (state), n); |
e7a72986 | 394 | } |
1bbd0b84 | 395 | #undef FUNC_NAME |
e7a72986 | 396 | |
a1ec6916 | 397 | SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0, |
1bbd0b84 | 398 | (SCM state), |
3b644514 | 399 | "Return a copy of the random state @var{state}.") |
1bbd0b84 | 400 | #define FUNC_NAME s_scm_copy_random_state |
e7a72986 MD |
401 | { |
402 | if (SCM_UNBNDP (state)) | |
86d31dfe | 403 | state = SCM_VARIABLE_REF (scm_var_random_state); |
3b3b36dd | 404 | SCM_VALIDATE_RSTATE (1,state); |
e7a72986 MD |
405 | return make_rstate (scm_the_rng.copy_rstate (SCM_RSTATE (state))); |
406 | } | |
1bbd0b84 | 407 | #undef FUNC_NAME |
e7a72986 | 408 | |
a1ec6916 | 409 | SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0, |
1bbd0b84 | 410 | (SCM seed), |
3b644514 | 411 | "Return a new random state using @var{seed}.") |
1bbd0b84 | 412 | #define FUNC_NAME s_scm_seed_to_random_state |
5ee11b7c MD |
413 | { |
414 | if (SCM_NUMBERP (seed)) | |
415 | seed = scm_number_to_string (seed, SCM_UNDEFINED); | |
3b3b36dd | 416 | SCM_VALIDATE_STRING (1,seed); |
34f0f2b8 | 417 | return make_rstate (scm_c_make_rstate (SCM_STRING_CHARS (seed), |
b7ead2ae | 418 | SCM_STRING_LENGTH (seed))); |
5ee11b7c | 419 | } |
1bbd0b84 | 420 | #undef FUNC_NAME |
5ee11b7c | 421 | |
a1ec6916 | 422 | SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0, |
1bbd0b84 | 423 | (SCM state), |
1e6808ea MG |
424 | "Return a uniformly distributed inexact real random number in\n" |
425 | "[0,1).") | |
1bbd0b84 | 426 | #define FUNC_NAME s_scm_random_uniform |
e7a72986 MD |
427 | { |
428 | if (SCM_UNBNDP (state)) | |
86d31dfe | 429 | state = SCM_VARIABLE_REF (scm_var_random_state); |
3b3b36dd | 430 | SCM_VALIDATE_RSTATE (1,state); |
f8de44c1 | 431 | return scm_make_real (scm_c_uniform01 (SCM_RSTATE (state))); |
e7a72986 | 432 | } |
1bbd0b84 | 433 | #undef FUNC_NAME |
e7a72986 | 434 | |
a1ec6916 | 435 | SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0, |
1bbd0b84 | 436 | (SCM state), |
1e6808ea MG |
437 | "Return an inexact real in a normal distribution. The\n" |
438 | "distribution used has mean 0 and standard deviation 1. For a\n" | |
439 | "normal distribution with mean m and standard deviation d use\n" | |
440 | "@code{(+ m (* d (random:normal)))}.") | |
1bbd0b84 | 441 | #define FUNC_NAME s_scm_random_normal |
afe5177e GH |
442 | { |
443 | if (SCM_UNBNDP (state)) | |
86d31dfe | 444 | state = SCM_VARIABLE_REF (scm_var_random_state); |
3b3b36dd | 445 | SCM_VALIDATE_RSTATE (1,state); |
f8de44c1 | 446 | return scm_make_real (scm_c_normal01 (SCM_RSTATE (state))); |
afe5177e | 447 | } |
1bbd0b84 | 448 | #undef FUNC_NAME |
afe5177e GH |
449 | |
450 | #ifdef HAVE_ARRAYS | |
451 | ||
e7a72986 MD |
452 | static void |
453 | vector_scale (SCM v, double c) | |
454 | { | |
b7ead2ae | 455 | int n = SCM_INUM (scm_uniform_vector_length (v)); |
e7a72986 MD |
456 | if (SCM_VECTORP (v)) |
457 | while (--n >= 0) | |
eb42e2f0 | 458 | SCM_REAL_VALUE (SCM_VELTS (v)[n]) *= c; |
e7a72986 MD |
459 | else |
460 | while (--n >= 0) | |
461 | ((double *) SCM_VELTS (v))[n] *= c; | |
462 | } | |
463 | ||
464 | static double | |
465 | vector_sum_squares (SCM v) | |
466 | { | |
467 | double x, sum = 0.0; | |
b7ead2ae | 468 | int n = SCM_INUM (scm_uniform_vector_length (v)); |
e7a72986 MD |
469 | if (SCM_VECTORP (v)) |
470 | while (--n >= 0) | |
471 | { | |
eb42e2f0 | 472 | x = SCM_REAL_VALUE (SCM_VELTS (v)[n]); |
e7a72986 MD |
473 | sum += x * x; |
474 | } | |
475 | else | |
476 | while (--n >= 0) | |
477 | { | |
478 | x = ((double *) SCM_VELTS (v))[n]; | |
479 | sum += x * x; | |
480 | } | |
481 | return sum; | |
482 | } | |
483 | ||
e7a72986 MD |
484 | /* For the uniform distribution on the solid sphere, note that in |
485 | * this distribution the length r of the vector has cumulative | |
486 | * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be | |
487 | * generated as r=u^(1/n). | |
488 | */ | |
a1ec6916 | 489 | SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0, |
1bbd0b84 | 490 | (SCM v, SCM state), |
d928e0b4 GB |
491 | "Fills vect with inexact real random numbers\n" |
492 | "the sum of whose squares is less than 1.0.\n" | |
493 | "Thinking of vect as coordinates in space of \n" | |
494 | "dimension n = (vector-length vect), the coordinates \n" | |
495 | "are uniformly distributed within the unit n-shere.\n" | |
64ba8e85 | 496 | "The sum of the squares of the numbers is returned.") |
1bbd0b84 | 497 | #define FUNC_NAME s_scm_random_solid_sphere_x |
e7a72986 | 498 | { |
3b3b36dd | 499 | SCM_VALIDATE_VECTOR_OR_DVECTOR (1,v); |
e7a72986 | 500 | if (SCM_UNBNDP (state)) |
86d31dfe | 501 | state = SCM_VARIABLE_REF (scm_var_random_state); |
3b3b36dd | 502 | SCM_VALIDATE_RSTATE (2,state); |
e7a72986 MD |
503 | scm_random_normal_vector_x (v, state); |
504 | vector_scale (v, | |
9b741bb6 | 505 | pow (scm_c_uniform01 (SCM_RSTATE (state)), |
b7ead2ae | 506 | 1.0 / SCM_INUM (scm_uniform_vector_length (v))) |
e7a72986 MD |
507 | / sqrt (vector_sum_squares (v))); |
508 | return SCM_UNSPECIFIED; | |
509 | } | |
1bbd0b84 | 510 | #undef FUNC_NAME |
e7a72986 | 511 | |
a1ec6916 | 512 | SCM_DEFINE (scm_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0, |
1bbd0b84 | 513 | (SCM v, SCM state), |
d928e0b4 GB |
514 | "Fills vect with inexact real random numbers\n" |
515 | "the sum of whose squares is equal to 1.0.\n" | |
516 | "Thinking of vect as coordinates in space of \n" | |
517 | "dimension n = (vector-length vect), the coordinates\n" | |
518 | "are uniformly distributed over the surface of the \n" | |
64ba8e85 | 519 | "unit n-shere.") |
1bbd0b84 | 520 | #define FUNC_NAME s_scm_random_hollow_sphere_x |
e7a72986 | 521 | { |
3b3b36dd | 522 | SCM_VALIDATE_VECTOR_OR_DVECTOR (1,v); |
e7a72986 | 523 | if (SCM_UNBNDP (state)) |
86d31dfe | 524 | state = SCM_VARIABLE_REF (scm_var_random_state); |
3b3b36dd | 525 | SCM_VALIDATE_RSTATE (2,state); |
e7a72986 MD |
526 | scm_random_normal_vector_x (v, state); |
527 | vector_scale (v, 1 / sqrt (vector_sum_squares (v))); | |
528 | return SCM_UNSPECIFIED; | |
529 | } | |
1bbd0b84 | 530 | #undef FUNC_NAME |
e7a72986 | 531 | |
1bbd0b84 | 532 | |
a1ec6916 | 533 | SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0, |
1bbd0b84 | 534 | (SCM v, SCM state), |
d928e0b4 GB |
535 | "Fills vect with inexact real random numbers that are\n" |
536 | "independent and standard normally distributed\n" | |
64ba8e85 | 537 | "(i.e., with mean 0 and variance 1).") |
1bbd0b84 | 538 | #define FUNC_NAME s_scm_random_normal_vector_x |
e7a72986 MD |
539 | { |
540 | int n; | |
3b3b36dd | 541 | SCM_VALIDATE_VECTOR_OR_DVECTOR (1,v); |
e7a72986 | 542 | if (SCM_UNBNDP (state)) |
86d31dfe | 543 | state = SCM_VARIABLE_REF (scm_var_random_state); |
3b3b36dd | 544 | SCM_VALIDATE_RSTATE (2,state); |
b7ead2ae | 545 | n = SCM_INUM (scm_uniform_vector_length (v)); |
e7a72986 MD |
546 | if (SCM_VECTORP (v)) |
547 | while (--n >= 0) | |
f8de44c1 | 548 | SCM_VELTS (v)[n] = scm_make_real (scm_c_normal01 (SCM_RSTATE (state))); |
e7a72986 MD |
549 | else |
550 | while (--n >= 0) | |
9b741bb6 | 551 | ((double *) SCM_VELTS (v))[n] = scm_c_normal01 (SCM_RSTATE (state)); |
e7a72986 MD |
552 | return SCM_UNSPECIFIED; |
553 | } | |
1bbd0b84 | 554 | #undef FUNC_NAME |
e7a72986 | 555 | |
afe5177e GH |
556 | #endif /* HAVE_ARRAYS */ |
557 | ||
a1ec6916 | 558 | SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0, |
1bbd0b84 | 559 | (SCM state), |
1e6808ea MG |
560 | "Return an inexact real in an exponential distribution with mean\n" |
561 | "1. For an exponential distribution with mean u use (* u\n" | |
562 | "(random:exp)).") | |
1bbd0b84 | 563 | #define FUNC_NAME s_scm_random_exp |
e7a72986 MD |
564 | { |
565 | if (SCM_UNBNDP (state)) | |
86d31dfe | 566 | state = SCM_VARIABLE_REF (scm_var_random_state); |
3b3b36dd | 567 | SCM_VALIDATE_RSTATE (1,state); |
f8de44c1 | 568 | return scm_make_real (scm_c_exp1 (SCM_RSTATE (state))); |
e7a72986 | 569 | } |
1bbd0b84 | 570 | #undef FUNC_NAME |
e7a72986 MD |
571 | |
572 | void | |
573 | scm_init_random () | |
574 | { | |
575 | int i, m; | |
576 | /* plug in default RNG */ | |
92c2555f | 577 | scm_t_rng rng = |
e7a72986 | 578 | { |
92c2555f | 579 | sizeof (scm_t_i_rstate), |
e7a72986 MD |
580 | (unsigned long (*)()) scm_i_uniform32, |
581 | (void (*)()) scm_i_init_rstate, | |
92c2555f | 582 | (scm_t_rstate *(*)()) scm_i_copy_rstate |
e7a72986 MD |
583 | }; |
584 | scm_the_rng = rng; | |
585 | ||
e841c3e0 KN |
586 | scm_tc16_rstate = scm_make_smob_type ("random-state", 0); |
587 | scm_set_smob_free (scm_tc16_rstate, rstate_free); | |
e7a72986 MD |
588 | |
589 | for (m = 1; m <= 0x100; m <<= 1) | |
590 | for (i = m >> 1; i < m; ++i) | |
591 | scm_masktab[i] = m - 1; | |
592 | ||
8dc9439f | 593 | #ifndef SCM_MAGIC_SNARFER |
a0599745 | 594 | #include "libguile/random.x" |
8dc9439f | 595 | #endif |
e7a72986 | 596 | |
2a0279c9 MD |
597 | /* Check that the assumptions about bits per bignum digit are correct. */ |
598 | #if SIZEOF_INT == 4 | |
599 | m = 16; | |
600 | #else | |
601 | m = 32; | |
602 | #endif | |
603 | if (m != SCM_BITSPERDIG) | |
604 | { | |
605 | fprintf (stderr, "Internal inconsistency: Confused about bignum digit size in random.c\n"); | |
606 | exit (1); | |
607 | } | |
608 | ||
e7a72986 MD |
609 | scm_add_feature ("random"); |
610 | } | |
89e00824 ML |
611 | |
612 | /* | |
613 | Local Variables: | |
614 | c-file-style: "gnu" | |
615 | End: | |
616 | */ |