1 /* Copyright (C) 1999,2000,2001, 2003, 2005, 2006, 2009, 2010 Free Software Foundation, Inc.
2 * This library is free software; you can redistribute it and/or
3 * modify it under the terms of the GNU Lesser General Public License
4 * as published by the Free Software Foundation; either version 3 of
5 * the License, or (at your option) any later version.
7 * This library is distributed in the hope that it will be useful, but
8 * WITHOUT ANY WARRANTY; without even the implied warranty of
9 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
10 * Lesser General Public License for more details.
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., 51 Franklin Street, Fifth Floor, Boston, MA
20 /* Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */
26 #include "libguile/_scm.h"
32 #include "libguile/smob.h"
33 #include "libguile/numbers.h"
34 #include "libguile/feature.h"
35 #include "libguile/strings.h"
36 #include "libguile/arrays.h"
37 #include "libguile/srfi-4.h"
38 #include "libguile/vectors.h"
39 #include "libguile/generalized-vectors.h"
41 #include "libguile/validate.h"
42 #include "libguile/random.h"
46 * A plugin interface for RNGs
48 * Using this interface, it is possible for the application to tell
49 * libguile to use a different RNG. This is desirable if it is
50 * necessary to use the same RNG everywhere in the application in
51 * order to prevent interference, if the application uses RNG
52 * hardware, or if the application has special demands on the RNG.
54 * Look in random.h and how the default generator is "plugged in" in
58 scm_t_rng scm_the_rng
;
64 * This is the MWC (Multiply With Carry) random number generator
65 * described by George Marsaglia at the Department of Statistics and
66 * Supercomputer Computations Research Institute, The Florida State
67 * University (http://stat.fsu.edu/~geo).
69 * It uses 64 bits, has a period of 4578426017172946943 (4.6e18), and
70 * passes all tests in the DIEHARD test suite
71 * (http://stat.fsu.edu/~geo/diehard.html)
74 #define A 2131995753UL
77 #define M_PI 3.14159265359
81 scm_i_uniform32 (scm_t_i_rstate
*state
)
83 scm_t_uint64 x
= (scm_t_uint64
) A
* state
->w
+ state
->c
;
84 scm_t_uint32 w
= x
& 0xffffffffUL
;
91 scm_i_init_rstate (scm_t_i_rstate
*state
, const char *seed
, int n
)
96 for (i
= 0; i
< n
; ++i
)
100 w
+= seed
[i
] << (8 * m
);
102 c
+= seed
[i
] << (8 * (m
- 4));
104 if ((w
== 0 && c
== 0) || (w
== -1 && c
== A
- 1))
111 scm_i_copy_rstate (scm_t_i_rstate
*state
)
113 scm_t_rstate
*new_state
;
115 new_state
= scm_gc_malloc_pointerless (scm_the_rng
.rstate_size
,
117 return memcpy (new_state
, state
, scm_the_rng
.rstate_size
);
120 SCM_SYMBOL(scm_i_rstate_tag
, "multiply-with-carry");
123 scm_i_init_rstate_scm (scm_t_i_rstate
*state
, SCM value
)
124 #define FUNC_NAME "scm_i_init_rstate_scm"
129 SCM_VALIDATE_LIST_COPYLEN (SCM_ARG1
, value
, length
);
130 SCM_ASSERT (length
== 3, value
, SCM_ARG1
, FUNC_NAME
);
131 SCM_ASSERT (scm_is_eq (SCM_CAR (value
), scm_i_rstate_tag
),
132 value
, SCM_ARG1
, FUNC_NAME
);
133 SCM_VALIDATE_ULONG_COPY (SCM_ARG1
, SCM_CADR (value
), w
);
134 SCM_VALIDATE_ULONG_COPY (SCM_ARG1
, SCM_CADDR (value
), c
);
142 scm_i_expose_rstate (scm_t_i_rstate
*state
)
144 return scm_list_3 (scm_i_rstate_tag
,
145 scm_from_ulong (state
->w
),
146 scm_from_ulong (state
->c
));
151 * Random number library functions
155 scm_c_make_rstate (const char *seed
, int n
)
159 state
= scm_gc_malloc_pointerless (scm_the_rng
.rstate_size
,
161 state
->reserved0
= 0;
162 scm_the_rng
.init_rstate (state
, seed
, n
);
167 scm_c_make_rstate_scm (SCM external
)
171 state
= scm_gc_malloc_pointerless (scm_the_rng
.rstate_size
,
173 state
->reserved0
= 0;
174 scm_the_rng
.init_rstate_scm (state
, external
);
179 scm_c_default_rstate ()
180 #define FUNC_NAME "scm_c_default_rstate"
182 SCM state
= SCM_VARIABLE_REF (scm_var_random_state
);
183 if (!SCM_RSTATEP (state
))
184 SCM_MISC_ERROR ("*random-state* contains bogus random state", SCM_EOL
);
185 return SCM_RSTATE (state
);
191 scm_c_uniform01 (scm_t_rstate
*state
)
193 double x
= (double) scm_the_rng
.random_bits (state
) / (double) 0xffffffffUL
;
194 return ((x
+ (double) scm_the_rng
.random_bits (state
))
195 / (double) 0xffffffffUL
);
199 scm_c_normal01 (scm_t_rstate
*state
)
201 if (state
->reserved0
)
203 state
->reserved0
= 0;
204 return state
->reserved1
;
210 r
= sqrt (-2.0 * log (scm_c_uniform01 (state
)));
211 a
= 2.0 * M_PI
* scm_c_uniform01 (state
);
214 state
->reserved1
= r
* cos (a
);
215 state
->reserved0
= 1;
222 scm_c_exp1 (scm_t_rstate
*state
)
224 return - log (scm_c_uniform01 (state
));
227 unsigned char scm_masktab
[256];
230 scm_c_random (scm_t_rstate
*state
, unsigned long m
)
232 unsigned long r
, mask
;
233 #if SCM_SIZEOF_UNSIGNED_LONG == 4
237 ? scm_masktab
[m
>> 8] << 8 | 0xff
239 ? scm_masktab
[m
>> 16] << 16 | 0xffff
240 : scm_masktab
[m
>> 24] << 24 | 0xffffff)));
241 while ((r
= scm_the_rng
.random_bits (state
) & mask
) >= m
);
242 #elif SCM_SIZEOF_UNSIGNED_LONG == 8
246 ? scm_masktab
[m
>> 8] << 8 | 0xff
248 ? scm_masktab
[m
>> 16] << 16 | 0xffff
250 ? scm_masktab
[m
>> 24] << 24 | 0xffffff
252 ? ((unsigned long) scm_masktab
[m
>> 32] << 32
255 ? ((unsigned long) scm_masktab
[m
>> 40] << 40
258 ? ((unsigned long) scm_masktab
[m
>> 48] << 48
260 : ((unsigned long) scm_masktab
[m
>> 56] << 56
261 | 0xffffffffffffffUL
))))))));
262 while ((r
= ((scm_the_rng
.random_bits (state
) << 32
263 | scm_the_rng
.random_bits (state
))) & mask
) >= m
);
265 #error "Cannot deal with this platform's unsigned long size"
271 SCM scm_c_random_bignum (scm_t_rstate *state, SCM m)
273 Takes a random state (source of random bits) and a bignum m.
274 Returns a bignum b, 0 <= b < m.
276 It does this by allocating a bignum b with as many base 65536 digits
277 as m, filling b with random bits (in 32 bit chunks) up to the most
278 significant 1 in m, and, finally checking if the resultant b is too
279 large (>= m). If too large, we simply repeat the process again. (It
280 is important to throw away all generated random bits if b >= m,
281 otherwise we'll end up with a distorted distribution.)
286 scm_c_random_bignum (scm_t_rstate
*state
, SCM m
)
288 SCM result
= scm_i_mkbig ();
289 const size_t m_bits
= mpz_sizeinbase (SCM_I_BIG_MPZ (m
), 2);
290 /* how many bits would only partially fill the last unsigned long? */
291 const size_t end_bits
= m_bits
% (sizeof (unsigned long) * SCM_CHAR_BIT
);
292 unsigned long *random_chunks
= NULL
;
293 const unsigned long num_full_chunks
=
294 m_bits
/ (sizeof (unsigned long) * SCM_CHAR_BIT
);
295 const unsigned long num_chunks
= num_full_chunks
+ ((end_bits
) ? 1 : 0);
297 /* we know the result will be this big */
298 mpz_realloc2 (SCM_I_BIG_MPZ (result
), m_bits
);
301 (unsigned long *) scm_gc_calloc (num_chunks
* sizeof (unsigned long),
302 "random bignum chunks");
306 unsigned long *current_chunk
= random_chunks
+ (num_chunks
- 1);
307 unsigned long chunks_left
= num_chunks
;
309 mpz_set_ui (SCM_I_BIG_MPZ (result
), 0);
313 /* generate a mask with ones in the end_bits position, i.e. if
314 end_bits is 3, then we'd have a mask of ...0000000111 */
315 const unsigned long rndbits
= scm_the_rng
.random_bits (state
);
316 int rshift
= (sizeof (unsigned long) * SCM_CHAR_BIT
) - end_bits
;
317 unsigned long mask
= ((unsigned long) ULONG_MAX
) >> rshift
;
318 unsigned long highest_bits
= rndbits
& mask
;
319 *current_chunk
-- = highest_bits
;
325 /* now fill in the remaining unsigned long sized chunks */
326 *current_chunk
-- = scm_the_rng
.random_bits (state
);
329 mpz_import (SCM_I_BIG_MPZ (result
),
332 sizeof (unsigned long),
336 /* if result >= m, regenerate it (it is important to regenerate
337 all bits in order not to get a distorted distribution) */
338 } while (mpz_cmp (SCM_I_BIG_MPZ (result
), SCM_I_BIG_MPZ (m
)) >= 0);
339 scm_gc_free (random_chunks
,
340 num_chunks
* sizeof (unsigned long),
341 "random bignum chunks");
342 return scm_i_normbig (result
);
346 * Scheme level representation of random states.
349 scm_t_bits scm_tc16_rstate
;
352 make_rstate (scm_t_rstate
*state
)
354 SCM_RETURN_NEWSMOB (scm_tc16_rstate
, state
);
359 * Scheme level interface.
362 SCM_GLOBAL_VARIABLE_INIT (scm_var_random_state
, "*random-state*", scm_seed_to_random_state (scm_from_locale_string ("URL:http://stat.fsu.edu/~geo/diehard.html")));
364 SCM_DEFINE (scm_random
, "random", 1, 1, 0,
366 "Return a number in [0, N).\n"
368 "Accepts a positive integer or real n and returns a\n"
369 "number of the same type between zero (inclusive) and\n"
370 "N (exclusive). The values returned have a uniform\n"
373 "The optional argument @var{state} must be of the type produced\n"
374 "by @code{seed->random-state}. It defaults to the value of the\n"
375 "variable @var{*random-state*}. This object is used to maintain\n"
376 "the state of the pseudo-random-number generator and is altered\n"
377 "as a side effect of the random operation.")
378 #define FUNC_NAME s_scm_random
380 if (SCM_UNBNDP (state
))
381 state
= SCM_VARIABLE_REF (scm_var_random_state
);
382 SCM_VALIDATE_RSTATE (2, state
);
385 unsigned long m
= SCM_I_INUM (n
);
386 SCM_ASSERT_RANGE (1, n
, m
> 0);
387 return scm_from_ulong (scm_c_random (SCM_RSTATE (state
), m
));
389 SCM_VALIDATE_NIM (1, n
);
391 return scm_from_double (SCM_REAL_VALUE (n
)
392 * scm_c_uniform01 (SCM_RSTATE (state
)));
395 SCM_WRONG_TYPE_ARG (1, n
);
396 return scm_c_random_bignum (SCM_RSTATE (state
), n
);
400 SCM_DEFINE (scm_copy_random_state
, "copy-random-state", 0, 1, 0,
402 "Return a copy of the random state @var{state}.")
403 #define FUNC_NAME s_scm_copy_random_state
405 if (SCM_UNBNDP (state
))
406 state
= SCM_VARIABLE_REF (scm_var_random_state
);
407 SCM_VALIDATE_RSTATE (1, state
);
408 return make_rstate (scm_the_rng
.copy_rstate (SCM_RSTATE (state
)));
412 SCM_DEFINE (scm_seed_to_random_state
, "seed->random-state", 1, 0, 0,
414 "Return a new random state using @var{seed}.")
415 #define FUNC_NAME s_scm_seed_to_random_state
418 if (SCM_NUMBERP (seed
))
419 seed
= scm_number_to_string (seed
, SCM_UNDEFINED
);
420 SCM_VALIDATE_STRING (1, seed
);
421 res
= make_rstate (scm_c_make_rstate (scm_i_string_chars (seed
),
422 scm_i_string_length (seed
)));
423 scm_remember_upto_here_1 (seed
);
429 SCM_DEFINE (scm_external_to_random_state
, "external->random-state", 1, 0, 0,
431 "Return a new random state using @var{external}.\n"
433 "@var{external} must be an external state representation obtained\n"
434 "from @code{random-state->external}.")
435 #define FUNC_NAME s_scm_external_to_random_state
437 return make_rstate (scm_c_make_rstate_scm (external
));
441 SCM_DEFINE (scm_random_state_to_external
, "random-state->external", 1, 0, 0,
443 "Return an external representation of @var{state}.")
444 #define FUNC_NAME s_scm_random_state_to_external
446 SCM_VALIDATE_RSTATE (1, state
);
447 return scm_the_rng
.expose_rstate (SCM_RSTATE (state
));
451 SCM_DEFINE (scm_random_uniform
, "random:uniform", 0, 1, 0,
453 "Return a uniformly distributed inexact real random number in\n"
455 #define FUNC_NAME s_scm_random_uniform
457 if (SCM_UNBNDP (state
))
458 state
= SCM_VARIABLE_REF (scm_var_random_state
);
459 SCM_VALIDATE_RSTATE (1, state
);
460 return scm_from_double (scm_c_uniform01 (SCM_RSTATE (state
)));
464 SCM_DEFINE (scm_random_normal
, "random:normal", 0, 1, 0,
466 "Return an inexact real in a normal distribution. The\n"
467 "distribution used has mean 0 and standard deviation 1. For a\n"
468 "normal distribution with mean m and standard deviation d use\n"
469 "@code{(+ m (* d (random:normal)))}.")
470 #define FUNC_NAME s_scm_random_normal
472 if (SCM_UNBNDP (state
))
473 state
= SCM_VARIABLE_REF (scm_var_random_state
);
474 SCM_VALIDATE_RSTATE (1, state
);
475 return scm_from_double (scm_c_normal01 (SCM_RSTATE (state
)));
480 vector_scale_x (SCM v
, double c
)
483 if (scm_is_simple_vector (v
))
485 n
= SCM_SIMPLE_VECTOR_LENGTH (v
);
487 SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v
, n
)) *= c
;
491 /* must be a f64vector. */
492 scm_t_array_handle handle
;
497 elts
= scm_f64vector_writable_elements (v
, &handle
, &len
, &inc
);
499 for (i
= 0; i
< len
; i
++, elts
+= inc
)
502 scm_array_handle_release (&handle
);
507 vector_sum_squares (SCM v
)
511 if (scm_is_simple_vector (v
))
513 n
= SCM_SIMPLE_VECTOR_LENGTH (v
);
516 x
= SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v
, n
));
522 /* must be a f64vector. */
523 scm_t_array_handle handle
;
528 elts
= scm_f64vector_elements (v
, &handle
, &len
, &inc
);
530 for (i
= 0; i
< len
; i
++, elts
+= inc
)
536 scm_array_handle_release (&handle
);
541 /* For the uniform distribution on the solid sphere, note that in
542 * this distribution the length r of the vector has cumulative
543 * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be
544 * generated as r=u^(1/n).
546 SCM_DEFINE (scm_random_solid_sphere_x
, "random:solid-sphere!", 1, 1, 0,
548 "Fills @var{vect} with inexact real random numbers the sum of\n"
549 "whose squares is less than 1.0. Thinking of @var{vect} as\n"
550 "coordinates in space of dimension @var{n} @math{=}\n"
551 "@code{(vector-length @var{vect})}, the coordinates are\n"
552 "uniformly distributed within the unit @var{n}-sphere.")
553 #define FUNC_NAME s_scm_random_solid_sphere_x
555 if (SCM_UNBNDP (state
))
556 state
= SCM_VARIABLE_REF (scm_var_random_state
);
557 SCM_VALIDATE_RSTATE (2, state
);
558 scm_random_normal_vector_x (v
, state
);
560 pow (scm_c_uniform01 (SCM_RSTATE (state
)),
561 1.0 / scm_c_generalized_vector_length (v
))
562 / sqrt (vector_sum_squares (v
)));
563 return SCM_UNSPECIFIED
;
567 SCM_DEFINE (scm_random_hollow_sphere_x
, "random:hollow-sphere!", 1, 1, 0,
569 "Fills vect with inexact real random numbers\n"
570 "the sum of whose squares is equal to 1.0.\n"
571 "Thinking of vect as coordinates in space of\n"
572 "dimension n = (vector-length vect), the coordinates\n"
573 "are uniformly distributed over the surface of the\n"
575 #define FUNC_NAME s_scm_random_hollow_sphere_x
577 if (SCM_UNBNDP (state
))
578 state
= SCM_VARIABLE_REF (scm_var_random_state
);
579 SCM_VALIDATE_RSTATE (2, state
);
580 scm_random_normal_vector_x (v
, state
);
581 vector_scale_x (v
, 1 / sqrt (vector_sum_squares (v
)));
582 return SCM_UNSPECIFIED
;
587 SCM_DEFINE (scm_random_normal_vector_x
, "random:normal-vector!", 1, 1, 0,
589 "Fills vect with inexact real random numbers that are\n"
590 "independent and standard normally distributed\n"
591 "(i.e., with mean 0 and variance 1).")
592 #define FUNC_NAME s_scm_random_normal_vector_x
595 scm_t_array_handle handle
;
596 scm_t_array_dim
*dim
;
598 if (SCM_UNBNDP (state
))
599 state
= SCM_VARIABLE_REF (scm_var_random_state
);
600 SCM_VALIDATE_RSTATE (2, state
);
602 scm_generalized_vector_get_handle (v
, &handle
);
603 dim
= scm_array_handle_dims (&handle
);
605 if (scm_is_vector (v
))
607 SCM
*elts
= scm_array_handle_writable_elements (&handle
);
608 for (i
= dim
->lbnd
; i
<= dim
->ubnd
; i
++, elts
+= dim
->inc
)
609 *elts
= scm_from_double (scm_c_normal01 (SCM_RSTATE (state
)));
613 /* must be a f64vector. */
614 double *elts
= scm_array_handle_f64_writable_elements (&handle
);
615 for (i
= dim
->lbnd
; i
<= dim
->ubnd
; i
++, elts
+= dim
->inc
)
616 *elts
= scm_c_normal01 (SCM_RSTATE (state
));
619 scm_array_handle_release (&handle
);
621 return SCM_UNSPECIFIED
;
625 SCM_DEFINE (scm_random_exp
, "random:exp", 0, 1, 0,
627 "Return an inexact real in an exponential distribution with mean\n"
628 "1. For an exponential distribution with mean u use (* u\n"
630 #define FUNC_NAME s_scm_random_exp
632 if (SCM_UNBNDP (state
))
633 state
= SCM_VARIABLE_REF (scm_var_random_state
);
634 SCM_VALIDATE_RSTATE (1, state
);
635 return scm_from_double (scm_c_exp1 (SCM_RSTATE (state
)));
643 /* plug in default RNG */
646 sizeof (scm_t_i_rstate
),
647 (unsigned long (*)()) scm_i_uniform32
,
648 (void (*)()) scm_i_init_rstate
,
649 (scm_t_rstate
*(*)()) scm_i_copy_rstate
,
650 (void (*)(scm_t_rstate
*, SCM
)) scm_i_init_rstate_scm
,
651 (SCM (*)(scm_t_rstate
*)) scm_i_expose_rstate
655 scm_tc16_rstate
= scm_make_smob_type ("random-state", 0);
657 for (m
= 1; m
<= 0x100; m
<<= 1)
658 for (i
= m
>> 1; i
< m
; ++i
)
659 scm_masktab
[i
] = m
- 1;
661 #include "libguile/random.x"
663 scm_add_feature ("random");