1 /* Copyright 1999-2001,2003,2005-2006,2009-2010,2012-2014,2017-2019
2 Free Software Foundation, Inc.
3
4 This file is part of Guile.
5
6 Guile is free software: you can redistribute it and/or modify it
7 under the terms of the GNU Lesser General Public License as published
8 by the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 Guile is distributed in the hope that it will be useful, but WITHOUT
12 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14 License for more details.
15
16 You should have received a copy of the GNU Lesser General Public
17 License along with Guile. If not, see
18 <https://www.gnu.org/licenses/>. */
19
20
21 /* Original Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */
22
23 #ifdef HAVE_CONFIG_H
24 # include <config.h>
25 #endif
26
27 #include <math.h>
28 #include <stdio.h>
29 #include <string.h>
30 #include <sys/types.h>
31 #include <unistd.h>
32
33 #include "scm.h"
34 #if SCM_ENABLE_MINI_GMP
35 #include "mini-gmp.h"
36 #else
37 #include <gmp.h>
38 #endif
39
40 #include "arrays.h"
41 #include "feature.h"
42 #include "generalized-arrays.h"
43 #include "generalized-vectors.h"
44 #include "gsubr.h"
45 #include "list.h"
46 #include "modules.h"
47 #include "numbers.h"
48 #include "numbers.h"
49 #include "pairs.h"
50 #include "smob.h"
51 #include "srfi-4.h"
52 #include "stime.h"
53 #include "strings.h"
54 #include "symbols.h"
55 #include "variable.h"
56 #include "vectors.h"
57
58 #include "random.h"
59
60
61
62 /*
63 * A plugin interface for RNGs
64 *
65 * Using this interface, it is possible for the application to tell
66 * libguile to use a different RNG. This is desirable if it is
67 * necessary to use the same RNG everywhere in the application in
68 * order to prevent interference, if the application uses RNG
69 * hardware, or if the application has special demands on the RNG.
70 *
71 * Look in random.h and how the default generator is "plugged in" in
72 * scm_init_random().
73 */
74
75 scm_t_rng scm_the_rng;
76
77
78 /*
79 * The prepackaged RNG
80 *
81 * This is the MWC (Multiply With Carry) random number generator
82 * described by George Marsaglia at the Department of Statistics and
83 * Supercomputer Computations Research Institute, The Florida State
84 * University (http://stat.fsu.edu/~geo).
85 *
86 * It uses 64 bits, has a period of 4578426017172946943 (4.6e18), and
87 * passes all tests in the DIEHARD test suite
88 * (http://stat.fsu.edu/~geo/diehard.html)
89 */
90
91 typedef struct scm_t_i_rstate {
92 scm_t_rstate rstate;
93 uint32_t w;
94 uint32_t c;
95 } scm_t_i_rstate;
96
97
98 #define A 2131995753UL
99
100 #ifndef M_PI
101 #define M_PI 3.14159265359
102 #endif
103
104 static uint32_t
scm_i_uniform32(scm_t_rstate * state)105 scm_i_uniform32 (scm_t_rstate *state)
106 {
107 scm_t_i_rstate *istate = (scm_t_i_rstate*) state;
108 uint64_t x = (uint64_t) A * istate->w + istate->c;
109 uint32_t w = x & 0xffffffffUL;
110 istate->w = w;
111 istate->c = x >> 32L;
112 return w;
113 }
114
115 static void
scm_i_init_rstate(scm_t_rstate * state,const char * seed,int n)116 scm_i_init_rstate (scm_t_rstate *state, const char *seed, int n)
117 {
118 scm_t_i_rstate *istate = (scm_t_i_rstate*) state;
119 uint32_t w = 0L;
120 uint32_t c = 0L;
121 int i, m;
122 for (i = 0; i < n; ++i)
123 {
124 m = i % 8;
125 if (m < 4)
126 w += seed[i] << (8 * m);
127 else
128 c += seed[i] << (8 * (m - 4));
129 }
130 if ((w == 0 && c == 0) || (w == -1 && c == A - 1))
131 ++c;
132 istate->w = w;
133 istate->c = c;
134 }
135
136 static scm_t_rstate *
scm_i_copy_rstate(scm_t_rstate * state)137 scm_i_copy_rstate (scm_t_rstate *state)
138 {
139 scm_t_rstate *new_state;
140
141 new_state = scm_gc_malloc_pointerless (state->rng->rstate_size,
142 "random-state");
143 return memcpy (new_state, state, state->rng->rstate_size);
144 }
145
146 SCM_SYMBOL(scm_i_rstate_tag, "multiply-with-carry");
147
148 static void
scm_i_rstate_from_datum(scm_t_rstate * state,SCM value)149 scm_i_rstate_from_datum (scm_t_rstate *state, SCM value)
150 #define FUNC_NAME "scm_i_rstate_from_datum"
151 {
152 scm_t_i_rstate *istate = (scm_t_i_rstate*) state;
153 uint32_t w, c;
154 long length;
155
156 SCM_VALIDATE_LIST_COPYLEN (SCM_ARG1, value, length);
157 SCM_ASSERT (length == 3, value, SCM_ARG1, FUNC_NAME);
158 SCM_ASSERT (scm_is_eq (SCM_CAR (value), scm_i_rstate_tag),
159 value, SCM_ARG1, FUNC_NAME);
160 SCM_VALIDATE_UINT_COPY (SCM_ARG1, SCM_CADR (value), w);
161 SCM_VALIDATE_UINT_COPY (SCM_ARG1, SCM_CADDR (value), c);
162
163 istate->w = w;
164 istate->c = c;
165 }
166 #undef FUNC_NAME
167
168 static SCM
scm_i_rstate_to_datum(scm_t_rstate * state)169 scm_i_rstate_to_datum (scm_t_rstate *state)
170 {
171 scm_t_i_rstate *istate = (scm_t_i_rstate*) state;
172 return scm_list_3 (scm_i_rstate_tag,
173 scm_from_uint32 (istate->w),
174 scm_from_uint32 (istate->c));
175 }
176
177
178 /*
179 * Random number library functions
180 */
181
182 scm_t_rstate *
scm_c_make_rstate(const char * seed,int n)183 scm_c_make_rstate (const char *seed, int n)
184 {
185 scm_t_rstate *state;
186
187 state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size,
188 "random-state");
189 state->rng = &scm_the_rng;
190 state->normal_next = 0.0;
191 state->rng->init_rstate (state, seed, n);
192 return state;
193 }
194
195 scm_t_rstate *
scm_c_rstate_from_datum(SCM datum)196 scm_c_rstate_from_datum (SCM datum)
197 {
198 scm_t_rstate *state;
199
200 state = scm_gc_malloc_pointerless (scm_the_rng.rstate_size,
201 "random-state");
202 state->rng = &scm_the_rng;
203 state->normal_next = 0.0;
204 state->rng->from_datum (state, datum);
205 return state;
206 }
207
208 scm_t_rstate *
scm_c_default_rstate()209 scm_c_default_rstate ()
210 #define FUNC_NAME "scm_c_default_rstate"
211 {
212 SCM state = SCM_VARIABLE_REF (scm_var_random_state);
213 if (!SCM_RSTATEP (state))
214 SCM_MISC_ERROR ("*random-state* contains bogus random state", SCM_EOL);
215 return SCM_RSTATE (state);
216 }
217 #undef FUNC_NAME
218
219
220 double
scm_c_uniform01(scm_t_rstate * state)221 scm_c_uniform01 (scm_t_rstate *state)
222 {
223 double x = (double) state->rng->random_bits (state) / (double) 0xffffffffUL;
224 return ((x + (double) state->rng->random_bits (state))
225 / (double) 0xffffffffUL);
226 }
227
228 double
scm_c_normal01(scm_t_rstate * state)229 scm_c_normal01 (scm_t_rstate *state)
230 {
231 if (state->normal_next != 0.0)
232 {
233 double ret = state->normal_next;
234
235 state->normal_next = 0.0;
236
237 return ret;
238 }
239 else
240 {
241 double r, a, n;
242
243 r = sqrt (-2.0 * log (scm_c_uniform01 (state)));
244 a = 2.0 * M_PI * scm_c_uniform01 (state);
245
246 n = r * sin (a);
247 state->normal_next = r * cos (a);
248
249 return n;
250 }
251 }
252
253 double
scm_c_exp1(scm_t_rstate * state)254 scm_c_exp1 (scm_t_rstate *state)
255 {
256 return - log (scm_c_uniform01 (state));
257 }
258
259 unsigned char scm_masktab[256];
260
261 static inline uint32_t
scm_i_mask32(uint32_t m)262 scm_i_mask32 (uint32_t m)
263 {
264 return (m < 0x100
265 ? scm_masktab[m]
266 : (m < 0x10000
267 ? scm_masktab[m >> 8] << 8 | 0xff
268 : (m < 0x1000000
269 ? scm_masktab[m >> 16] << 16 | 0xffff
270 : ((uint32_t) scm_masktab[m >> 24]) << 24 | 0xffffff)));
271 }
272
273 uint32_t
scm_c_random(scm_t_rstate * state,uint32_t m)274 scm_c_random (scm_t_rstate *state, uint32_t m)
275 {
276 uint32_t r, mask = scm_i_mask32 (m);
277 while ((r = state->rng->random_bits (state) & mask) >= m);
278 return r;
279 }
280
281 uint64_t
scm_c_random64(scm_t_rstate * state,uint64_t m)282 scm_c_random64 (scm_t_rstate *state, uint64_t m)
283 {
284 uint64_t r;
285 uint32_t mask;
286
287 if (m <= UINT32_MAX)
288 return scm_c_random (state, (uint32_t) m);
289
290 mask = scm_i_mask32 (m >> 32);
291 while ((r = ((uint64_t) (state->rng->random_bits (state) & mask) << 32)
292 | state->rng->random_bits (state)) >= m)
293 ;
294 return r;
295 }
296
297 /*
298 SCM scm_c_random_bignum (scm_t_rstate *state, SCM m)
299
300 Takes a random state (source of random bits) and a bignum m.
301 Returns a bignum b, 0 <= b < m.
302
303 It does this by allocating a bignum b with as many base 65536 digits
304 as m, filling b with random bits (in 32 bit chunks) up to the most
305 significant 1 in m, and, finally checking if the resultant b is too
306 large (>= m). If too large, we simply repeat the process again. (It
307 is important to throw away all generated random bits if b >= m,
308 otherwise we'll end up with a distorted distribution.)
309
310 */
311
312 SCM
scm_c_random_bignum(scm_t_rstate * state,SCM m)313 scm_c_random_bignum (scm_t_rstate *state, SCM m)
314 {
315 SCM result = scm_i_mkbig ();
316 const size_t m_bits = mpz_sizeinbase (SCM_I_BIG_MPZ (m), 2);
317 /* how many bits would only partially fill the last uint32_t? */
318 const size_t end_bits = m_bits % (sizeof (uint32_t) * SCM_CHAR_BIT);
319 uint32_t *random_chunks = NULL;
320 const uint32_t num_full_chunks =
321 m_bits / (sizeof (uint32_t) * SCM_CHAR_BIT);
322 const uint32_t num_chunks = num_full_chunks + ((end_bits) ? 1 : 0);
323
324 /* we know the result will be this big */
325 mpz_realloc2 (SCM_I_BIG_MPZ (result), m_bits);
326
327 random_chunks =
328 (uint32_t *) scm_gc_calloc (num_chunks * sizeof (uint32_t),
329 "random bignum chunks");
330
331 do
332 {
333 uint32_t *current_chunk = random_chunks + (num_chunks - 1);
334 uint32_t chunks_left = num_chunks;
335
336 mpz_set_ui (SCM_I_BIG_MPZ (result), 0);
337
338 if (end_bits)
339 {
340 /* generate a mask with ones in the end_bits position, i.e. if
341 end_bits is 3, then we'd have a mask of ...0000000111 */
342 const uint32_t rndbits = state->rng->random_bits (state);
343 int rshift = (sizeof (uint32_t) * SCM_CHAR_BIT) - end_bits;
344 uint32_t mask = ((uint32_t)-1) >> rshift;
345 uint32_t highest_bits = rndbits & mask;
346 *current_chunk-- = highest_bits;
347 chunks_left--;
348 }
349
350 while (chunks_left)
351 {
352 /* now fill in the remaining uint32_t sized chunks */
353 *current_chunk-- = state->rng->random_bits (state);
354 chunks_left--;
355 }
356 mpz_import (SCM_I_BIG_MPZ (result),
357 num_chunks,
358 -1,
359 sizeof (uint32_t),
360 0,
361 0,
362 random_chunks);
363 /* if result >= m, regenerate it (it is important to regenerate
364 all bits in order not to get a distorted distribution) */
365 } while (mpz_cmp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (m)) >= 0);
366 scm_gc_free (random_chunks,
367 num_chunks * sizeof (uint32_t),
368 "random bignum chunks");
369 return scm_i_normbig (result);
370 }
371
372 /*
373 * Scheme level representation of random states.
374 */
375
376 scm_t_bits scm_tc16_rstate;
377
378 static SCM
make_rstate(scm_t_rstate * state)379 make_rstate (scm_t_rstate *state)
380 {
381 SCM_RETURN_NEWSMOB (scm_tc16_rstate, state);
382 }
383
384
385 /*
386 * Scheme level interface.
387 */
388
389 SCM_GLOBAL_VARIABLE_INIT (scm_var_random_state, "*random-state*",
390 scm_seed_to_random_state
391 (scm_from_utf8_string
392 ("URL:http://stat.fsu.edu/~geo/diehard.html")));
393
394 SCM_DEFINE (scm_random, "random", 1, 1, 0,
395 (SCM n, SCM state),
396 "Return a number in [0, N).\n"
397 "\n"
398 "Accepts a positive integer or real n and returns a\n"
399 "number of the same type between zero (inclusive) and\n"
400 "N (exclusive). The values returned have a uniform\n"
401 "distribution.\n"
402 "\n"
403 "The optional argument @var{state} must be of the type produced\n"
404 "by @code{seed->random-state}. It defaults to the value of the\n"
405 "variable @var{*random-state*}. This object is used to maintain\n"
406 "the state of the pseudo-random-number generator and is altered\n"
407 "as a side effect of the random operation.")
408 #define FUNC_NAME s_scm_random
409 {
410 if (SCM_UNBNDP (state))
411 state = SCM_VARIABLE_REF (scm_var_random_state);
412 SCM_VALIDATE_RSTATE (2, state);
413 if (SCM_I_INUMP (n))
414 {
415 scm_t_bits m = (scm_t_bits) SCM_I_INUM (n);
416 SCM_ASSERT_RANGE (1, n, SCM_I_INUM (n) > 0);
417 #if SCM_SIZEOF_UINTPTR_T <= 4
418 return scm_from_uint32 (scm_c_random (SCM_RSTATE (state),
419 (uint32_t) m));
420 #elif SCM_SIZEOF_UINTPTR_T <= 8
421 return scm_from_uint64 (scm_c_random64 (SCM_RSTATE (state),
422 (uint64_t) m));
423 #else
424 #error "Cannot deal with this platform's scm_t_bits size"
425 #endif
426 }
427 if (SCM_REALP (n))
428 return scm_from_double (SCM_REAL_VALUE (n)
429 * scm_c_uniform01 (SCM_RSTATE (state)));
430
431 if (!SCM_BIGP (n))
432 SCM_WRONG_TYPE_ARG (1, n);
433 return scm_c_random_bignum (SCM_RSTATE (state), n);
434 }
435 #undef FUNC_NAME
436
437 SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0,
438 (SCM state),
439 "Return a copy of the random state @var{state}.")
440 #define FUNC_NAME s_scm_copy_random_state
441 {
442 if (SCM_UNBNDP (state))
443 state = SCM_VARIABLE_REF (scm_var_random_state);
444 SCM_VALIDATE_RSTATE (1, state);
445 return make_rstate (SCM_RSTATE (state)->rng->copy_rstate (SCM_RSTATE (state)));
446 }
447 #undef FUNC_NAME
448
449 SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0,
450 (SCM seed),
451 "Return a new random state using @var{seed}.")
452 #define FUNC_NAME s_scm_seed_to_random_state
453 {
454 SCM res;
455 char *c_str;
456 size_t len;
457
458 if (SCM_NUMBERP (seed))
459 seed = scm_number_to_string (seed, SCM_UNDEFINED);
460 SCM_VALIDATE_STRING (1, seed);
461
462 if (scm_i_is_narrow_string (seed))
463 /* This special case of a narrow string, where latin1 is used, is
464 for backward compatibility during the 2.2 stable series. In
465 future major releases, we should use UTF-8 uniformly. */
466 c_str = scm_to_latin1_stringn (seed, &len);
467 else
468 c_str = scm_to_utf8_stringn (seed, &len);
469
470 /* 'scm_to_*_stringn' returns a 'size_t' for the length in bytes, but
471 'scm_c_make_rstate' accepts an 'int'. Make sure the length fits in
472 an 'int'. */
473 if (len > INT_MAX)
474 {
475 free (c_str);
476 SCM_OUT_OF_RANGE (1, seed);
477 }
478
479 res = make_rstate (scm_c_make_rstate (c_str, len));
480 free (c_str);
481
482 scm_remember_upto_here_1 (seed);
483 return res;
484
485 }
486 #undef FUNC_NAME
487
488 SCM_DEFINE (scm_datum_to_random_state, "datum->random-state", 1, 0, 0,
489 (SCM datum),
490 "Return a new random state using @var{datum}, which should have\n"
491 "been obtained from @code{random-state->datum}.")
492 #define FUNC_NAME s_scm_datum_to_random_state
493 {
494 return make_rstate (scm_c_rstate_from_datum (datum));
495 }
496 #undef FUNC_NAME
497
498 SCM_DEFINE (scm_random_state_to_datum, "random-state->datum", 1, 0, 0,
499 (SCM state),
500 "Return a datum representation of @var{state} that may be\n"
501 "written out and read back with the Scheme reader.")
502 #define FUNC_NAME s_scm_random_state_to_datum
503 {
504 SCM_VALIDATE_RSTATE (1, state);
505 return SCM_RSTATE (state)->rng->to_datum (SCM_RSTATE (state));
506 }
507 #undef FUNC_NAME
508
509 SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0,
510 (SCM state),
511 "Return a uniformly distributed inexact real random number in\n"
512 "[0,1).")
513 #define FUNC_NAME s_scm_random_uniform
514 {
515 if (SCM_UNBNDP (state))
516 state = SCM_VARIABLE_REF (scm_var_random_state);
517 SCM_VALIDATE_RSTATE (1, state);
518 return scm_from_double (scm_c_uniform01 (SCM_RSTATE (state)));
519 }
520 #undef FUNC_NAME
521
522 SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0,
523 (SCM state),
524 "Return an inexact real in a normal distribution. The\n"
525 "distribution used has mean 0 and standard deviation 1. For a\n"
526 "normal distribution with mean m and standard deviation d use\n"
527 "@code{(+ m (* d (random:normal)))}.")
528 #define FUNC_NAME s_scm_random_normal
529 {
530 if (SCM_UNBNDP (state))
531 state = SCM_VARIABLE_REF (scm_var_random_state);
532 SCM_VALIDATE_RSTATE (1, state);
533 return scm_from_double (scm_c_normal01 (SCM_RSTATE (state)));
534 }
535 #undef FUNC_NAME
536
537 static void
vector_scale_x(SCM v,double c)538 vector_scale_x (SCM v, double c)
539 {
540 scm_t_array_handle handle;
541 scm_t_array_dim const * dims;
542 ssize_t i, inc, ubnd;
543
544 scm_array_get_handle (v, &handle);
545 dims = scm_array_handle_dims (&handle);
546 if (1 == scm_array_handle_rank (&handle))
547 {
548 ubnd = dims[0].ubnd;
549 inc = dims[0].inc;
550
551 if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_F64)
552 {
553 double *elts = (double *)(handle.writable_elements) + handle.base;
554 for (i = dims[0].lbnd; i <= ubnd; ++i, elts += inc)
555 *elts *= c;
556 return;
557 }
558 else if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_SCM)
559 {
560 SCM *elts = (SCM *)(handle.writable_elements) + handle.base;
561 for (i = dims[0].lbnd; i <= ubnd; ++i, elts += inc)
562 SCM_REAL_VALUE (*elts) *= c;
563 return;
564 }
565 }
566 scm_array_handle_release (&handle);
567 scm_misc_error (NULL, "must be a rank-1 array of type #t or 'f64", scm_list_1 (v));
568 }
569
570 static double
vector_sum_squares(SCM v)571 vector_sum_squares (SCM v)
572 {
573 double x, sum = 0.0;
574 scm_t_array_handle handle;
575 scm_t_array_dim const * dims;
576 ssize_t i, inc, ubnd;
577
578 scm_array_get_handle (v, &handle);
579 dims = scm_array_handle_dims (&handle);
580 if (1 == scm_array_handle_rank (&handle))
581 {
582 ubnd = dims[0].ubnd;
583 inc = dims[0].inc;
584 if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_F64)
585 {
586 const double *elts = (const double *)(handle.elements) + handle.base;
587 for (i = dims[0].lbnd; i <= ubnd; ++i, elts += inc)
588 {
589 x = *elts;
590 sum += x * x;
591 }
592 return sum;
593 }
594 else if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_SCM)
595 {
596 const SCM *elts = (const SCM *)(handle.elements) + handle.base;
597 for (i = dims[0].lbnd; i <= ubnd; ++i, elts += inc)
598 {
599 x = SCM_REAL_VALUE (*elts);
600 sum += x * x;
601 }
602 return sum;
603 }
604 }
605 scm_array_handle_release (&handle);
606 scm_misc_error (NULL, "must be an array of type #t or 'f64", scm_list_1 (v));
607 }
608
609 /* For the uniform distribution on the solid sphere, note that in
610 * this distribution the length r of the vector has cumulative
611 * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be
612 * generated as r=u^(1/n).
613 */
614 SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0,
615 (SCM v, SCM state),
616 "Fills @var{vect} with inexact real random numbers the sum of\n"
617 "whose squares is less than 1.0. Thinking of @var{vect} as\n"
618 "coordinates in space of dimension @var{n} @math{=}\n"
619 "@code{(vector-length @var{vect})}, the coordinates are\n"
620 "uniformly distributed within the unit @var{n}-sphere.")
621 #define FUNC_NAME s_scm_random_solid_sphere_x
622 {
623 if (SCM_UNBNDP (state))
624 state = SCM_VARIABLE_REF (scm_var_random_state);
625 SCM_VALIDATE_RSTATE (2, state);
626 scm_random_normal_vector_x (v, state);
627 vector_scale_x (v,
628 pow (scm_c_uniform01 (SCM_RSTATE (state)),
629 1.0 / scm_c_array_length (v))
630 / sqrt (vector_sum_squares (v)));
631 return SCM_UNSPECIFIED;
632 }
633 #undef FUNC_NAME
634
635 SCM_DEFINE (scm_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0,
636 (SCM v, SCM state),
637 "Fills vect with inexact real random numbers\n"
638 "the sum of whose squares is equal to 1.0.\n"
639 "Thinking of vect as coordinates in space of\n"
640 "dimension n = (vector-length vect), the coordinates\n"
641 "are uniformly distributed over the surface of the\n"
642 "unit n-sphere.")
643 #define FUNC_NAME s_scm_random_hollow_sphere_x
644 {
645 if (SCM_UNBNDP (state))
646 state = SCM_VARIABLE_REF (scm_var_random_state);
647 SCM_VALIDATE_RSTATE (2, state);
648 scm_random_normal_vector_x (v, state);
649 vector_scale_x (v, 1 / sqrt (vector_sum_squares (v)));
650 return SCM_UNSPECIFIED;
651 }
652 #undef FUNC_NAME
653
654
655 SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0,
656 (SCM v, SCM state),
657 "Fills vect with inexact real random numbers that are\n"
658 "independent and standard normally distributed\n"
659 "(i.e., with mean 0 and variance 1).")
660 #define FUNC_NAME s_scm_random_normal_vector_x
661 {
662 scm_t_array_handle handle;
663 scm_t_array_dim const * dims;
664 ssize_t i;
665
666 if (SCM_UNBNDP (state))
667 state = SCM_VARIABLE_REF (scm_var_random_state);
668 SCM_VALIDATE_RSTATE (2, state);
669
670 scm_array_get_handle (v, &handle);
671 if (1 != scm_array_handle_rank (&handle))
672 {
673 scm_array_handle_release (&handle);
674 scm_wrong_type_arg_msg (NULL, 0, v, "rank 1 array");
675 }
676
677 dims = scm_array_handle_dims (&handle);
678
679 if (handle.element_type == SCM_ARRAY_ELEMENT_TYPE_SCM)
680 {
681 SCM *elts = scm_array_handle_writable_elements (&handle);
682 for (i = dims->lbnd; i <= dims->ubnd; i++, elts += dims->inc)
683 *elts = scm_from_double (scm_c_normal01 (SCM_RSTATE (state)));
684 }
685 else
686 {
687 /* must be a f64vector. */
688 double *elts = scm_array_handle_f64_writable_elements (&handle);
689 for (i = dims->lbnd; i <= dims->ubnd; i++, elts += dims->inc)
690 *elts = scm_c_normal01 (SCM_RSTATE (state));
691 }
692
693 scm_array_handle_release (&handle);
694 return SCM_UNSPECIFIED;
695 }
696 #undef FUNC_NAME
697
698 SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0,
699 (SCM state),
700 "Return an inexact real in an exponential distribution with mean\n"
701 "1. For an exponential distribution with mean u use (* u\n"
702 "(random:exp)).")
703 #define FUNC_NAME s_scm_random_exp
704 {
705 if (SCM_UNBNDP (state))
706 state = SCM_VARIABLE_REF (scm_var_random_state);
707 SCM_VALIDATE_RSTATE (1, state);
708 return scm_from_double (scm_c_exp1 (SCM_RSTATE (state)));
709 }
710 #undef FUNC_NAME
711
712 /* Return a new random-state seeded from the time, date, process ID, an
713 address from a freshly allocated heap cell, an address from the local
714 stack frame, and a high-resolution timer if available. This is only
715 to be used as a last resort, when no better source of entropy is
716 available. */
717 static SCM
random_state_of_last_resort(void)718 random_state_of_last_resort (void)
719 {
720 SCM state;
721 SCM time_of_day = scm_gettimeofday ();
722 SCM sources = scm_list_n
723 (scm_from_unsigned_integer (SCM_UNPACK (time_of_day)), /* heap addr */
724 /* Avoid scm_getpid, since it depends on HAVE_POSIX. */
725 scm_from_unsigned_integer (getpid ()), /* process ID */
726 scm_get_internal_real_time (), /* high-resolution process timer */
727 scm_from_unsigned_integer ((scm_t_bits) &time_of_day), /* stack addr */
728 scm_car (time_of_day), /* seconds since midnight 1970-01-01 UTC */
729 scm_cdr (time_of_day), /* microsecond component of the above clock */
730 SCM_UNDEFINED);
731
732 /* Concatenate the sources bitwise to form the seed */
733 SCM seed = SCM_INUM0;
734 while (scm_is_pair (sources))
735 {
736 seed = scm_logxor (seed, scm_ash (scm_car (sources),
737 scm_integer_length (seed)));
738 sources = scm_cdr (sources);
739 }
740
741 /* FIXME The following code belongs in `scm_seed_to_random_state',
742 and here we should simply do:
743
744 return scm_seed_to_random_state (seed);
745
746 Unfortunately, `scm_seed_to_random_state' only preserves around 32
747 bits of entropy from the provided seed. I don't know if it's okay
748 to fix that in 2.0, so for now we have this workaround. */
749 {
750 int i, len;
751 unsigned char *buf;
752 len = scm_to_int (scm_ceiling_quotient (scm_integer_length (seed),
753 SCM_I_MAKINUM (8)));
754 buf = (unsigned char *) malloc (len);
755 for (i = len-1; i >= 0; --i)
756 {
757 buf[i] = scm_to_int (scm_logand (seed, SCM_I_MAKINUM (255)));
758 seed = scm_ash (seed, SCM_I_MAKINUM (-8));
759 }
760 state = make_rstate (scm_c_make_rstate ((char *) buf, len));
761 free (buf);
762 }
763 return state;
764 }
765
766 /* Attempt to fill buffer with random bytes from /dev/urandom.
767 Return 1 if successful, else return 0. */
768 static int
read_dev_urandom(unsigned char * buf,size_t len)769 read_dev_urandom (unsigned char *buf, size_t len)
770 {
771 size_t res = 0;
772 FILE *f = fopen ("/dev/urandom", "r");
773 if (f)
774 {
775 res = fread(buf, 1, len, f);
776 fclose (f);
777 }
778 return (res == len);
779 }
780
781 /* Fill a buffer with random bytes seeded from a platform-specific
782 source of entropy. /dev/urandom is used if available. Note that
783 this function provides no guarantees about the amount of entropy
784 present in the returned bytes. */
785 void
scm_i_random_bytes_from_platform(unsigned char * buf,size_t len)786 scm_i_random_bytes_from_platform (unsigned char *buf, size_t len)
787 {
788 if (read_dev_urandom (buf, len))
789 return;
790 else /* FIXME: support other platform sources */
791 {
792 /* When all else fails, use this (rather weak) fallback */
793 SCM random_state = random_state_of_last_resort ();
794 int i;
795 for (i = len-1; i >= 0; --i)
796 buf[i] = scm_to_int (scm_random (SCM_I_MAKINUM (256), random_state));
797 }
798 }
799
800 SCM_DEFINE (scm_random_state_from_platform, "random-state-from-platform", 0, 0, 0,
801 (void),
802 "Construct a new random state seeded from a platform-specific\n\
803 source of entropy, appropriate for use in non-security-critical applications.")
804 #define FUNC_NAME s_scm_random_state_from_platform
805 {
806 unsigned char buf[32];
807 if (read_dev_urandom (buf, sizeof(buf)))
808 return make_rstate (scm_c_make_rstate ((char *) buf, sizeof(buf)));
809 else
810 return random_state_of_last_resort ();
811 }
812 #undef FUNC_NAME
813
814 void
scm_init_random()815 scm_init_random ()
816 {
817 int i, m;
818 /* plug in default RNG */
819 scm_t_rng rng =
820 {
821 sizeof (scm_t_i_rstate),
822 scm_i_uniform32,
823 scm_i_init_rstate,
824 scm_i_copy_rstate,
825 scm_i_rstate_from_datum,
826 scm_i_rstate_to_datum
827 };
828 scm_the_rng = rng;
829
830 scm_tc16_rstate = scm_make_smob_type ("random-state", 0);
831
832 for (m = 1; m <= 0x100; m <<= 1)
833 for (i = m >> 1; i < m; ++i)
834 scm_masktab[i] = m - 1;
835
836 #include "random.x"
837
838 scm_add_feature ("random");
839 }
840