Commit | Line | Data |
---|---|---|
e7a72986 MD |
1 | /* Copyright (C) 1999 Free Software Foundation, Inc. |
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 | /* Software engineering face-lift by Greg J. Badros, 11-Dec-1999, |
42 | gjb@cs.washington.edu, http://www.cs.washington.edu/homes/gjb */ | |
43 | ||
44 | ||
e7a72986 MD |
45 | /* Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */ |
46 | ||
e7a72986 MD |
47 | #include "_scm.h" |
48 | ||
7f146094 | 49 | #include <stdio.h> |
e7a72986 MD |
50 | #include <math.h> |
51 | #include "genio.h" | |
52 | #include "smob.h" | |
53 | #include "numbers.h" | |
54 | #include "feature.h" | |
55 | ||
1bbd0b84 | 56 | #include "scm_validate.h" |
e7a72986 MD |
57 | #include "random.h" |
58 | ||
59 | \f | |
60 | /* | |
61 | * A plugin interface for RNGs | |
62 | * | |
63 | * Using this interface, it is possible for the application to tell | |
64 | * libguile to use a different RNG. This is desirable if it is | |
65 | * necessary to use the same RNG everywhere in the application in | |
66 | * order to prevent interference, if the application uses RNG | |
67 | * hardware, or if the application has special demands on the RNG. | |
68 | * | |
69 | * Look in random.h and how the default generator is "plugged in" in | |
70 | * scm_init_random(). | |
71 | */ | |
72 | ||
73 | scm_rng scm_the_rng; | |
74 | ||
75 | \f | |
76 | /* | |
77 | * The prepackaged RNG | |
78 | * | |
79 | * This is the MWC (Multiply With Carry) random number generator | |
80 | * described by George Marsaglia at the Department of Statistics and | |
81 | * Supercomputer Computations Research Institute, The Florida State | |
82 | * University (http://stat.fsu.edu/~geo). | |
83 | * | |
84 | * It uses 64 bits, has a period of 4578426017172946943 (4.6e18), and | |
85 | * passes all tests in the DIEHARD test suite | |
86 | * (http://stat.fsu.edu/~geo/diehard.html) | |
87 | */ | |
88 | ||
89 | #define A 2131995753UL | |
90 | ||
91 | #if SIZEOF_LONG > 4 | |
92 | #if SIZEOF_INT > 4 | |
93 | #define LONG32 unsigned short | |
94 | #else | |
95 | #define LONG32 unsigned int | |
96 | #endif | |
97 | #define LONG64 unsigned long | |
98 | #else | |
99 | #define LONG32 unsigned long | |
100 | #define LONG64 unsigned long long | |
101 | #endif | |
102 | ||
103 | #if SIZEOF_LONG > 4 || defined (HAVE_LONG_LONGS) | |
104 | ||
105 | unsigned long | |
106 | scm_i_uniform32 (scm_i_rstate *state) | |
107 | { | |
108 | LONG64 x = (LONG64) A * state->w + state->c; | |
109 | LONG32 w = x & 0xffffffffUL; | |
110 | state->w = w; | |
111 | state->c = x >> 32L; | |
112 | return w; | |
113 | } | |
114 | ||
115 | #else | |
116 | ||
117 | /* ww This is a portable version of the same RNG without 64 bit | |
118 | * * aa arithmetic. | |
119 | * ---- | |
120 | * xx It is only intended to provide identical behaviour on | |
121 | * xx platforms without 8 byte longs or long longs until | |
122 | * xx someone has implemented the routine in assembler code. | |
123 | * xxcc | |
124 | * ---- | |
125 | * ccww | |
126 | */ | |
127 | ||
128 | #define L(x) ((x) & 0xffff) | |
129 | #define H(x) ((x) >> 16) | |
130 | ||
131 | unsigned long | |
132 | scm_i_uniform32 (scm_i_rstate *state) | |
133 | { | |
134 | LONG32 x1 = L (A) * L (state->w); | |
135 | LONG32 x2 = L (A) * H (state->w); | |
136 | LONG32 x3 = H (A) * L (state->w); | |
137 | LONG32 w = L (x1) + L (state->c); | |
138 | LONG32 m = H (x1) + L (x2) + L (x3) + H (state->c) + H (w); | |
139 | LONG32 x4 = H (A) * H (state->w); | |
140 | state->w = w = (L (m) << 16) + L (w); | |
141 | state->c = H (x2) + H (x3) + x4 + H (m); | |
142 | return w; | |
143 | } | |
144 | ||
145 | #endif | |
146 | ||
147 | void | |
148 | scm_i_init_rstate (scm_i_rstate *state, char *seed, int n) | |
149 | { | |
150 | LONG32 w = 0L; | |
151 | LONG32 c = 0L; | |
152 | int i, m; | |
153 | for (i = 0; i < n; ++i) | |
154 | { | |
155 | m = i % 8; | |
156 | if (m < 4) | |
157 | w += seed[i] << (8 * m); | |
158 | else | |
159 | c += seed[i] << (8 * (m - 4)); | |
160 | } | |
161 | if ((w == 0 && c == 0) || (w == 0xffffffffUL && c == A - 1)) | |
162 | ++c; | |
163 | state->w = w; | |
164 | state->c = c; | |
165 | } | |
166 | ||
167 | scm_i_rstate * | |
168 | scm_i_copy_rstate (scm_i_rstate *state) | |
169 | { | |
170 | scm_rstate *new_state = malloc (scm_the_rng.rstate_size); | |
171 | if (new_state == 0) | |
172 | scm_wta (SCM_MAKINUM (scm_the_rng.rstate_size), | |
173 | (char *) SCM_NALLOC, "rstate"); | |
174 | return memcpy (new_state, state, scm_the_rng.rstate_size); | |
175 | } | |
176 | ||
177 | \f | |
178 | /* | |
179 | * Random number library functions | |
180 | */ | |
181 | ||
5ee11b7c | 182 | scm_rstate * |
9b741bb6 | 183 | scm_c_make_rstate (char *seed, int n) |
5ee11b7c MD |
184 | { |
185 | scm_rstate *state = malloc (scm_the_rng.rstate_size); | |
186 | if (state == 0) | |
187 | scm_wta (SCM_MAKINUM (scm_the_rng.rstate_size), | |
188 | (char *) SCM_NALLOC, | |
189 | "rstate"); | |
190 | state->reserved0 = 0; | |
191 | scm_the_rng.init_rstate (state, seed, n); | |
192 | return state; | |
193 | } | |
194 | ||
9b741bb6 MD |
195 | scm_rstate * |
196 | scm_c_default_rstate () | |
197 | { | |
198 | SCM state = SCM_CDR (scm_var_random_state); | |
0c95b57d | 199 | SCM_ASSERT (SCM_RSTATEP (state), |
9b741bb6 MD |
200 | state, "*random-state* contains bogus random state", 0); |
201 | return SCM_RSTATE (state); | |
202 | } | |
203 | ||
e7a72986 | 204 | inline double |
9b741bb6 | 205 | scm_c_uniform01 (scm_rstate *state) |
e7a72986 | 206 | { |
5a92ddfd | 207 | double x = (double) scm_the_rng.random_bits (state) / (double) 0xffffffffUL; |
e7a72986 | 208 | return ((x + (double) scm_the_rng.random_bits (state)) |
5a92ddfd | 209 | / (double) 0xffffffffUL); |
e7a72986 MD |
210 | } |
211 | ||
212 | double | |
9b741bb6 | 213 | scm_c_normal01 (scm_rstate *state) |
e7a72986 MD |
214 | { |
215 | if (state->reserved0) | |
216 | { | |
217 | state->reserved0 = 0; | |
218 | return state->reserved1; | |
219 | } | |
220 | else | |
221 | { | |
222 | double r, a, n; | |
e7a72986 | 223 | |
9b741bb6 MD |
224 | r = sqrt (-2.0 * log (scm_c_uniform01 (state))); |
225 | a = 2.0 * M_PI * scm_c_uniform01 (state); | |
e7a72986 MD |
226 | |
227 | n = r * sin (a); | |
228 | state->reserved1 = r * cos (a); | |
5a92ddfd | 229 | state->reserved0 = 1; |
e7a72986 MD |
230 | |
231 | return n; | |
232 | } | |
233 | } | |
234 | ||
235 | double | |
9b741bb6 | 236 | scm_c_exp1 (scm_rstate *state) |
e7a72986 | 237 | { |
9b741bb6 | 238 | return - log (scm_c_uniform01 (state)); |
e7a72986 MD |
239 | } |
240 | ||
241 | unsigned char scm_masktab[256]; | |
242 | ||
243 | unsigned long | |
9b741bb6 | 244 | scm_c_random (scm_rstate *state, unsigned long m) |
e7a72986 MD |
245 | { |
246 | unsigned int r, mask; | |
247 | mask = (m < 0x100 | |
248 | ? scm_masktab[m] | |
249 | : (m < 0x10000 | |
5a92ddfd | 250 | ? scm_masktab[m >> 8] << 8 | 0xff |
e7a72986 | 251 | : (m < 0x1000000 |
5a92ddfd MD |
252 | ? scm_masktab[m >> 16] << 16 | 0xffff |
253 | : scm_masktab[m >> 24] << 24 | 0xffffff))); | |
e7a72986 MD |
254 | while ((r = scm_the_rng.random_bits (state) & mask) >= m); |
255 | return r; | |
256 | } | |
257 | ||
258 | SCM | |
9b741bb6 | 259 | scm_c_random_bignum (scm_rstate *state, SCM m) |
e7a72986 MD |
260 | { |
261 | SCM b; | |
a7e7ea3e MD |
262 | int i, nd; |
263 | LONG32 *bits, mask, w; | |
264 | nd = SCM_NUMDIGS (m); | |
265 | /* calculate mask for most significant digit */ | |
e7a72986 MD |
266 | #if SIZEOF_INT == 4 |
267 | /* 16 bit digits */ | |
a7e7ea3e | 268 | if (nd & 1) |
e7a72986 MD |
269 | { |
270 | /* fix most significant 16 bits */ | |
a7e7ea3e | 271 | unsigned short s = SCM_BDIGITS (m)[nd - 1]; |
5a92ddfd | 272 | mask = s < 0x100 ? scm_masktab[s] : scm_masktab[s >> 8] << 8 | 0xff; |
e7a72986 MD |
273 | } |
274 | else | |
275 | #endif | |
276 | { | |
277 | /* fix most significant 32 bits */ | |
96e263d6 | 278 | #if SIZEOF_INT == 4 |
2a0279c9 MD |
279 | w = SCM_BDIGITS (m)[nd - 1] << 16 | SCM_BDIGITS (m)[nd - 2]; |
280 | #else | |
96e263d6 | 281 | w = SCM_BDIGITS (m)[nd - 1]; |
2a0279c9 | 282 | #endif |
e7a72986 MD |
283 | mask = (w < 0x10000 |
284 | ? (w < 0x100 | |
285 | ? scm_masktab[w] | |
5a92ddfd | 286 | : scm_masktab[w >> 8] << 8 | 0xff) |
e7a72986 | 287 | : (w < 0x1000000 |
5a92ddfd MD |
288 | ? scm_masktab[w >> 16] << 16 | 0xffff |
289 | : scm_masktab[w >> 24] << 24 | 0xffffff)); | |
e7a72986 | 290 | } |
a7e7ea3e MD |
291 | b = scm_mkbig (nd, 0); |
292 | bits = (LONG32 *) SCM_BDIGITS (b); | |
293 | do | |
e7a72986 | 294 | { |
a7e7ea3e MD |
295 | i = nd; |
296 | /* treat most significant digit specially */ | |
297 | #if SIZEOF_INT == 4 | |
298 | /* 16 bit digits */ | |
299 | if (i & 1) | |
300 | { | |
301 | ((SCM_BIGDIG*) bits)[i - 1] = scm_the_rng.random_bits (state) & mask; | |
302 | i /= 2; | |
303 | } | |
304 | else | |
305 | #endif | |
306 | { | |
307 | /* fix most significant 32 bits */ | |
96e263d6 | 308 | #if SIZEOF_INT == 4 |
2a0279c9 | 309 | w = scm_the_rng.random_bits (state) & mask; |
96e263d6 MD |
310 | ((SCM_BIGDIG*) bits)[i - 2] = w & 0xffff; |
311 | ((SCM_BIGDIG*) bits)[i - 1] = w >> 16; | |
312 | i = i / 2 - 1; | |
2a0279c9 | 313 | #else |
96e263d6 | 314 | i /= 2; |
a7e7ea3e | 315 | bits[--i] = scm_the_rng.random_bits (state) & mask; |
2a0279c9 | 316 | #endif |
a7e7ea3e MD |
317 | } |
318 | /* now fill up the rest of the bignum */ | |
319 | while (i) | |
320 | bits[--i] = scm_the_rng.random_bits (state); | |
321 | b = scm_normbig (b); | |
322 | if (SCM_INUMP (b)) | |
323 | return b; | |
324 | } while (scm_bigcomp (b, m) <= 0); | |
325 | return b; | |
e7a72986 MD |
326 | } |
327 | ||
328 | /* | |
329 | * Scheme level representation of random states. | |
330 | */ | |
331 | ||
332 | long scm_tc16_rstate; | |
333 | ||
334 | static SCM | |
335 | make_rstate (scm_rstate *state) | |
336 | { | |
23a62151 | 337 | SCM_RETURN_NEWSMOB (scm_tc16_rstate, state); |
e7a72986 MD |
338 | } |
339 | ||
340 | static scm_sizet | |
341 | free_rstate (SCM rstate) | |
342 | { | |
343 | free (SCM_RSTATE (rstate)); | |
344 | return scm_the_rng.rstate_size; | |
345 | } | |
346 | ||
e7a72986 MD |
347 | /* |
348 | * Scheme level interface. | |
349 | */ | |
350 | ||
5ee11b7c | 351 | SCM_GLOBAL_VCELL_INIT (scm_var_random_state, "*random-state*", scm_seed_to_random_state (scm_makfrom0str ("URL:http://stat.fsu.edu/~geo/diehard.html"))); |
e7a72986 | 352 | |
a1ec6916 | 353 | SCM_DEFINE (scm_random, "random", 1, 1, 0, |
1bbd0b84 | 354 | (SCM n, SCM state), |
717050c8 | 355 | "") |
1bbd0b84 | 356 | #define FUNC_NAME s_scm_random |
e7a72986 MD |
357 | { |
358 | if (SCM_UNBNDP (state)) | |
359 | state = SCM_CDR (scm_var_random_state); | |
3b3b36dd | 360 | SCM_VALIDATE_RSTATE (2,state); |
e7a72986 MD |
361 | if (SCM_INUMP (n)) |
362 | { | |
363 | unsigned long m = SCM_INUM (n); | |
1bbd0b84 | 364 | SCM_ASSERT_RANGE (1,n,m > 0); |
9b741bb6 | 365 | return SCM_MAKINUM (scm_c_random (SCM_RSTATE (state), m)); |
e7a72986 | 366 | } |
6b5a304f | 367 | SCM_VALIDATE_NIM (1,n); |
e7a72986 | 368 | if (SCM_REALP (n)) |
9b741bb6 | 369 | return scm_makdbl (SCM_REALPART (n) * scm_c_uniform01 (SCM_RSTATE (state)), |
e7a72986 | 370 | 0.0); |
1bbd0b84 | 371 | SCM_VALIDATE_SMOB (1,n,bigpos); |
9b741bb6 | 372 | return scm_c_random_bignum (SCM_RSTATE (state), n); |
e7a72986 | 373 | } |
1bbd0b84 | 374 | #undef FUNC_NAME |
e7a72986 | 375 | |
a1ec6916 | 376 | SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0, |
1bbd0b84 | 377 | (SCM state), |
717050c8 | 378 | "") |
1bbd0b84 | 379 | #define FUNC_NAME s_scm_copy_random_state |
e7a72986 MD |
380 | { |
381 | if (SCM_UNBNDP (state)) | |
5ee11b7c | 382 | state = SCM_CDR (scm_var_random_state); |
3b3b36dd | 383 | SCM_VALIDATE_RSTATE (1,state); |
e7a72986 MD |
384 | return make_rstate (scm_the_rng.copy_rstate (SCM_RSTATE (state))); |
385 | } | |
1bbd0b84 | 386 | #undef FUNC_NAME |
e7a72986 | 387 | |
a1ec6916 | 388 | SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0, |
1bbd0b84 | 389 | (SCM seed), |
717050c8 | 390 | "") |
1bbd0b84 | 391 | #define FUNC_NAME s_scm_seed_to_random_state |
5ee11b7c MD |
392 | { |
393 | if (SCM_NUMBERP (seed)) | |
394 | seed = scm_number_to_string (seed, SCM_UNDEFINED); | |
3b3b36dd | 395 | SCM_VALIDATE_STRING (1,seed); |
9b741bb6 | 396 | return make_rstate (scm_c_make_rstate (SCM_ROCHARS (seed), |
5ee11b7c MD |
397 | SCM_LENGTH (seed))); |
398 | } | |
1bbd0b84 | 399 | #undef FUNC_NAME |
5ee11b7c | 400 | |
a1ec6916 | 401 | SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0, |
1bbd0b84 | 402 | (SCM state), |
717050c8 | 403 | "") |
1bbd0b84 | 404 | #define FUNC_NAME s_scm_random_uniform |
e7a72986 MD |
405 | { |
406 | if (SCM_UNBNDP (state)) | |
407 | state = SCM_CDR (scm_var_random_state); | |
3b3b36dd | 408 | SCM_VALIDATE_RSTATE (1,state); |
9b741bb6 | 409 | return scm_makdbl (scm_c_uniform01 (SCM_RSTATE (state)), 0.0); |
e7a72986 | 410 | } |
1bbd0b84 | 411 | #undef FUNC_NAME |
e7a72986 | 412 | |
a1ec6916 | 413 | SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0, |
1bbd0b84 | 414 | (SCM state), |
717050c8 | 415 | "") |
1bbd0b84 | 416 | #define FUNC_NAME s_scm_random_normal |
afe5177e GH |
417 | { |
418 | if (SCM_UNBNDP (state)) | |
419 | state = SCM_CDR (scm_var_random_state); | |
3b3b36dd | 420 | SCM_VALIDATE_RSTATE (1,state); |
afe5177e GH |
421 | return scm_makdbl (scm_c_normal01 (SCM_RSTATE (state)), 0.0); |
422 | } | |
1bbd0b84 | 423 | #undef FUNC_NAME |
afe5177e GH |
424 | |
425 | #ifdef HAVE_ARRAYS | |
426 | ||
e7a72986 MD |
427 | static void |
428 | vector_scale (SCM v, double c) | |
429 | { | |
430 | int n = SCM_LENGTH (v); | |
431 | if (SCM_VECTORP (v)) | |
432 | while (--n >= 0) | |
433 | SCM_REAL (SCM_VELTS (v)[n]) *= c; | |
434 | else | |
435 | while (--n >= 0) | |
436 | ((double *) SCM_VELTS (v))[n] *= c; | |
437 | } | |
438 | ||
439 | static double | |
440 | vector_sum_squares (SCM v) | |
441 | { | |
442 | double x, sum = 0.0; | |
443 | int n = SCM_LENGTH (v); | |
444 | if (SCM_VECTORP (v)) | |
445 | while (--n >= 0) | |
446 | { | |
447 | x = SCM_REAL (SCM_VELTS (v)[n]); | |
448 | sum += x * x; | |
449 | } | |
450 | else | |
451 | while (--n >= 0) | |
452 | { | |
453 | x = ((double *) SCM_VELTS (v))[n]; | |
454 | sum += x * x; | |
455 | } | |
456 | return sum; | |
457 | } | |
458 | ||
e7a72986 MD |
459 | /* For the uniform distribution on the solid sphere, note that in |
460 | * this distribution the length r of the vector has cumulative | |
461 | * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be | |
462 | * generated as r=u^(1/n). | |
463 | */ | |
a1ec6916 | 464 | SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0, |
1bbd0b84 | 465 | (SCM v, SCM state), |
717050c8 | 466 | "") |
1bbd0b84 | 467 | #define FUNC_NAME s_scm_random_solid_sphere_x |
e7a72986 | 468 | { |
3b3b36dd | 469 | SCM_VALIDATE_VECTOR_OR_DVECTOR (1,v); |
e7a72986 MD |
470 | if (SCM_UNBNDP (state)) |
471 | state = SCM_CDR (scm_var_random_state); | |
3b3b36dd | 472 | SCM_VALIDATE_RSTATE (2,state); |
e7a72986 MD |
473 | scm_random_normal_vector_x (v, state); |
474 | vector_scale (v, | |
9b741bb6 | 475 | pow (scm_c_uniform01 (SCM_RSTATE (state)), |
e7a72986 MD |
476 | 1.0 / SCM_LENGTH (v)) |
477 | / sqrt (vector_sum_squares (v))); | |
478 | return SCM_UNSPECIFIED; | |
479 | } | |
1bbd0b84 | 480 | #undef FUNC_NAME |
e7a72986 | 481 | |
a1ec6916 | 482 | SCM_DEFINE (scm_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0, |
1bbd0b84 | 483 | (SCM v, SCM state), |
717050c8 | 484 | "") |
1bbd0b84 | 485 | #define FUNC_NAME s_scm_random_hollow_sphere_x |
e7a72986 | 486 | { |
3b3b36dd | 487 | SCM_VALIDATE_VECTOR_OR_DVECTOR (1,v); |
e7a72986 MD |
488 | if (SCM_UNBNDP (state)) |
489 | state = SCM_CDR (scm_var_random_state); | |
3b3b36dd | 490 | SCM_VALIDATE_RSTATE (2,state); |
e7a72986 MD |
491 | scm_random_normal_vector_x (v, state); |
492 | vector_scale (v, 1 / sqrt (vector_sum_squares (v))); | |
493 | return SCM_UNSPECIFIED; | |
494 | } | |
1bbd0b84 | 495 | #undef FUNC_NAME |
e7a72986 | 496 | |
1bbd0b84 | 497 | |
a1ec6916 | 498 | SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0, |
1bbd0b84 GB |
499 | (SCM v, SCM state), |
500 | "") | |
501 | #define FUNC_NAME s_scm_random_normal_vector_x | |
e7a72986 MD |
502 | { |
503 | int n; | |
3b3b36dd | 504 | SCM_VALIDATE_VECTOR_OR_DVECTOR (1,v); |
e7a72986 MD |
505 | if (SCM_UNBNDP (state)) |
506 | state = SCM_CDR (scm_var_random_state); | |
3b3b36dd | 507 | SCM_VALIDATE_RSTATE (2,state); |
e7a72986 MD |
508 | n = SCM_LENGTH (v); |
509 | if (SCM_VECTORP (v)) | |
510 | while (--n >= 0) | |
9b741bb6 | 511 | SCM_VELTS (v)[n] = scm_makdbl (scm_c_normal01 (SCM_RSTATE (state)), 0.0); |
e7a72986 MD |
512 | else |
513 | while (--n >= 0) | |
9b741bb6 | 514 | ((double *) SCM_VELTS (v))[n] = scm_c_normal01 (SCM_RSTATE (state)); |
e7a72986 MD |
515 | return SCM_UNSPECIFIED; |
516 | } | |
1bbd0b84 | 517 | #undef FUNC_NAME |
e7a72986 | 518 | |
afe5177e GH |
519 | #endif /* HAVE_ARRAYS */ |
520 | ||
a1ec6916 | 521 | SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0, |
1bbd0b84 | 522 | (SCM state), |
717050c8 | 523 | "") |
1bbd0b84 | 524 | #define FUNC_NAME s_scm_random_exp |
e7a72986 MD |
525 | { |
526 | if (SCM_UNBNDP (state)) | |
527 | state = SCM_CDR (scm_var_random_state); | |
3b3b36dd | 528 | SCM_VALIDATE_RSTATE (1,state); |
9b741bb6 | 529 | return scm_makdbl (scm_c_exp1 (SCM_RSTATE (state)), 0.0); |
e7a72986 | 530 | } |
1bbd0b84 | 531 | #undef FUNC_NAME |
e7a72986 MD |
532 | |
533 | void | |
534 | scm_init_random () | |
535 | { | |
536 | int i, m; | |
537 | /* plug in default RNG */ | |
538 | scm_rng rng = | |
539 | { | |
540 | sizeof (scm_i_rstate), | |
541 | (unsigned long (*)()) scm_i_uniform32, | |
542 | (void (*)()) scm_i_init_rstate, | |
543 | (scm_rstate *(*)()) scm_i_copy_rstate | |
544 | }; | |
545 | scm_the_rng = rng; | |
546 | ||
23a62151 MD |
547 | scm_tc16_rstate = scm_make_smob_type_mfpe ("random-state", 0, |
548 | NULL, free_rstate, NULL, NULL); | |
e7a72986 MD |
549 | |
550 | for (m = 1; m <= 0x100; m <<= 1) | |
551 | for (i = m >> 1; i < m; ++i) | |
552 | scm_masktab[i] = m - 1; | |
553 | ||
554 | #include "random.x" | |
555 | ||
2a0279c9 MD |
556 | /* Check that the assumptions about bits per bignum digit are correct. */ |
557 | #if SIZEOF_INT == 4 | |
558 | m = 16; | |
559 | #else | |
560 | m = 32; | |
561 | #endif | |
562 | if (m != SCM_BITSPERDIG) | |
563 | { | |
564 | fprintf (stderr, "Internal inconsistency: Confused about bignum digit size in random.c\n"); | |
565 | exit (1); | |
566 | } | |
567 | ||
e7a72986 MD |
568 | scm_add_feature ("random"); |
569 | } |