1 /* Copyright (C) 1999,2000,2001, 2003 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
4 * License as published by the Free Software Foundation; either
5 * version 2.1 of the License, or (at your option) any later version.
7 * This library 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 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 /* Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */
25 #include "libguile/_scm.h"
31 #include "libguile/smob.h"
32 #include "libguile/numbers.h"
33 #include "libguile/feature.h"
34 #include "libguile/strings.h"
35 #include "libguile/unif.h"
36 #include "libguile/vectors.h"
38 #include "libguile/validate.h"
39 #include "libguile/random.h"
43 * A plugin interface for RNGs
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.
51 * Look in random.h and how the default generator is "plugged in" in
55 scm_t_rng scm_the_rng
;
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).
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)
71 #define A 2131995753UL
74 #define M_PI 3.14159265359
77 #ifdef SCM_HAVE_T_INT64
80 scm_i_uniform32 (scm_t_i_rstate
*state
)
82 scm_t_int64 x
= (scm_t_int64
) A
* state
->w
+ state
->c
;
83 scm_t_int32 w
= x
& 0xffffffffUL
;
91 /* ww This is a portable version of the same RNG without 64 bit
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.
102 #define L(x) ((x) & 0xffff)
103 #define H(x) ((x) >> 16)
106 scm_i_uniform32 (scm_t_i_rstate
*state
)
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
);
114 state
->w
= w
= (L (m
) << 16) + L (w
);
115 state
->c
= H (x2
) + H (x3
) + x4
+ H (m
);
122 scm_i_init_rstate (scm_t_i_rstate
*state
, char *seed
, int n
)
127 for (i
= 0; i
< n
; ++i
)
131 w
+= seed
[i
] << (8 * m
);
133 c
+= seed
[i
] << (8 * (m
- 4));
135 if ((w
== 0 && c
== 0) || (w
== 0xffffffffUL
&& c
== A
- 1))
142 scm_i_copy_rstate (scm_t_i_rstate
*state
)
144 scm_t_rstate
*new_state
= scm_malloc (scm_the_rng
.rstate_size
);
146 scm_memory_error ("rstate");
147 return memcpy (new_state
, state
, scm_the_rng
.rstate_size
);
152 * Random number library functions
156 scm_c_make_rstate (char *seed
, int n
)
158 scm_t_rstate
*state
= scm_malloc (scm_the_rng
.rstate_size
);
160 scm_memory_error ("rstate");
161 state
->reserved0
= 0;
162 scm_the_rng
.init_rstate (state
, seed
, n
);
168 scm_c_default_rstate ()
169 #define FUNC_NAME "scm_c_default_rstate"
171 SCM state
= SCM_VARIABLE_REF (scm_var_random_state
);
172 if (!SCM_RSTATEP (state
))
173 SCM_MISC_ERROR ("*random-state* contains bogus random state", SCM_EOL
);
174 return SCM_RSTATE (state
);
180 scm_c_uniform01 (scm_t_rstate
*state
)
182 double x
= (double) scm_the_rng
.random_bits (state
) / (double) 0xffffffffUL
;
183 return ((x
+ (double) scm_the_rng
.random_bits (state
))
184 / (double) 0xffffffffUL
);
188 scm_c_normal01 (scm_t_rstate
*state
)
190 if (state
->reserved0
)
192 state
->reserved0
= 0;
193 return state
->reserved1
;
199 r
= sqrt (-2.0 * log (scm_c_uniform01 (state
)));
200 a
= 2.0 * M_PI
* scm_c_uniform01 (state
);
203 state
->reserved1
= r
* cos (a
);
204 state
->reserved0
= 1;
211 scm_c_exp1 (scm_t_rstate
*state
)
213 return - log (scm_c_uniform01 (state
));
216 unsigned char scm_masktab
[256];
219 scm_c_random (scm_t_rstate
*state
, unsigned long m
)
221 unsigned int r
, mask
;
225 ? scm_masktab
[m
>> 8] << 8 | 0xff
227 ? scm_masktab
[m
>> 16] << 16 | 0xffff
228 : scm_masktab
[m
>> 24] << 24 | 0xffffff)));
229 while ((r
= scm_the_rng
.random_bits (state
) & mask
) >= m
);
234 SCM scm_c_random_bignum (scm_t_rstate *state, SCM m)
236 Takes a random state (source of random bits) and a bignum m.
237 Returns a bignum b, 0 <= b < m.
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.)
249 scm_c_random_bignum (scm_t_rstate
*state
, SCM m
)
251 SCM result
= scm_i_mkbig ();
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? */
254 const size_t end_bits
= m_bits
% (sizeof (unsigned long) * SCM_CHAR_BIT
);
255 unsigned long *random_chunks
= NULL
;
256 const unsigned long num_full_chunks
=
257 m_bits
/ (sizeof (unsigned long) * SCM_CHAR_BIT
);
258 const unsigned long num_chunks
= num_full_chunks
+ ((end_bits
) ? 1 : 0);
260 /* we know the result will be this big */
261 mpz_realloc2 (SCM_I_BIG_MPZ (result
), m_bits
);
264 (unsigned long *) scm_gc_calloc (num_chunks
* sizeof (unsigned long),
265 "random bignum chunks");
269 unsigned long *current_chunk
= random_chunks
+ (num_chunks
- 1);
270 unsigned long chunks_left
= num_chunks
;
272 mpz_set_ui (SCM_I_BIG_MPZ (result
), 0);
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
);
279 int rshift
= (sizeof (unsigned long) * SCM_CHAR_BIT
) - end_bits
;
280 unsigned long mask
= ((unsigned long) ULONG_MAX
) >> rshift
;
281 unsigned long highest_bits
= rndbits
& mask
;
282 *current_chunk
-- = highest_bits
;
288 /* now fill in the remaining unsigned long sized chunks */
289 *current_chunk
-- = scm_the_rng
.random_bits (state
);
292 mpz_import (SCM_I_BIG_MPZ (result
),
295 sizeof (unsigned long),
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);
302 scm_gc_free (random_chunks
,
303 num_chunks
* sizeof (unsigned long),
304 "random bignum chunks");
305 return scm_i_normbig (result
);
309 * Scheme level representation of random states.
312 scm_t_bits scm_tc16_rstate
;
315 make_rstate (scm_t_rstate
*state
)
317 SCM_RETURN_NEWSMOB (scm_tc16_rstate
, state
);
321 rstate_free (SCM rstate
)
323 free (SCM_RSTATE (rstate
));
328 * Scheme level interface.
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")));
333 SCM_DEFINE (scm_random
, "random", 1, 1, 0,
335 "Return a number in [0, N).\n"
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"
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.")
347 #define FUNC_NAME s_scm_random
349 if (SCM_UNBNDP (state
))
350 state
= SCM_VARIABLE_REF (scm_var_random_state
);
351 SCM_VALIDATE_RSTATE (2, state
);
354 unsigned long m
= SCM_INUM (n
);
355 SCM_ASSERT_RANGE (1, n
, m
> 0);
356 return SCM_I_MAKINUM (scm_c_random (SCM_RSTATE (state
), m
));
358 SCM_VALIDATE_NIM (1, n
);
360 return scm_make_real (SCM_REAL_VALUE (n
)
361 * scm_c_uniform01 (SCM_RSTATE (state
)));
363 SCM_VALIDATE_BIGINT (1, n
);
364 return scm_c_random_bignum (SCM_RSTATE (state
), n
);
368 SCM_DEFINE (scm_copy_random_state
, "copy-random-state", 0, 1, 0,
370 "Return a copy of the random state @var{state}.")
371 #define FUNC_NAME s_scm_copy_random_state
373 if (SCM_UNBNDP (state
))
374 state
= SCM_VARIABLE_REF (scm_var_random_state
);
375 SCM_VALIDATE_RSTATE (1, state
);
376 return make_rstate (scm_the_rng
.copy_rstate (SCM_RSTATE (state
)));
380 SCM_DEFINE (scm_seed_to_random_state
, "seed->random-state", 1, 0, 0,
382 "Return a new random state using @var{seed}.")
383 #define FUNC_NAME s_scm_seed_to_random_state
385 if (SCM_NUMBERP (seed
))
386 seed
= scm_number_to_string (seed
, SCM_UNDEFINED
);
387 SCM_VALIDATE_STRING (1, seed
);
388 return make_rstate (scm_c_make_rstate (SCM_STRING_CHARS (seed
),
389 SCM_STRING_LENGTH (seed
)));
393 SCM_DEFINE (scm_random_uniform
, "random:uniform", 0, 1, 0,
395 "Return a uniformly distributed inexact real random number in\n"
397 #define FUNC_NAME s_scm_random_uniform
399 if (SCM_UNBNDP (state
))
400 state
= SCM_VARIABLE_REF (scm_var_random_state
);
401 SCM_VALIDATE_RSTATE (1, state
);
402 return scm_make_real (scm_c_uniform01 (SCM_RSTATE (state
)));
406 SCM_DEFINE (scm_random_normal
, "random:normal", 0, 1, 0,
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)))}.")
412 #define FUNC_NAME s_scm_random_normal
414 if (SCM_UNBNDP (state
))
415 state
= SCM_VARIABLE_REF (scm_var_random_state
);
416 SCM_VALIDATE_RSTATE (1, state
);
417 return scm_make_real (scm_c_normal01 (SCM_RSTATE (state
)));
424 vector_scale (SCM v
, double c
)
426 int n
= SCM_INUM (scm_uniform_vector_length (v
));
429 SCM_REAL_VALUE (SCM_VELTS (v
)[n
]) *= c
;
432 ((double *) SCM_VELTS (v
))[n
] *= c
;
436 vector_sum_squares (SCM v
)
439 int n
= SCM_INUM (scm_uniform_vector_length (v
));
443 x
= SCM_REAL_VALUE (SCM_VELTS (v
)[n
]);
449 x
= ((double *) SCM_VELTS (v
))[n
];
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).
460 SCM_DEFINE (scm_random_solid_sphere_x
, "random:solid-sphere!", 1, 1, 0,
462 "Fills vect with inexact real random numbers\n"
463 "the sum of whose squares is less than 1.0.\n"
464 "Thinking of vect as coordinates in space of\n"
465 "dimension n = (vector-length vect), the coordinates\n"
466 "are uniformly distributed within the unit n-sphere.\n"
467 "The sum of the squares of the numbers is returned.")
468 #define FUNC_NAME s_scm_random_solid_sphere_x
470 SCM_VALIDATE_VECTOR_OR_DVECTOR (1, v
);
471 if (SCM_UNBNDP (state
))
472 state
= SCM_VARIABLE_REF (scm_var_random_state
);
473 SCM_VALIDATE_RSTATE (2, state
);
474 scm_random_normal_vector_x (v
, state
);
476 pow (scm_c_uniform01 (SCM_RSTATE (state
)),
477 1.0 / SCM_INUM (scm_uniform_vector_length (v
)))
478 / sqrt (vector_sum_squares (v
)));
479 return SCM_UNSPECIFIED
;
483 SCM_DEFINE (scm_random_hollow_sphere_x
, "random:hollow-sphere!", 1, 1, 0,
485 "Fills vect with inexact real random numbers\n"
486 "the sum of whose squares is equal to 1.0.\n"
487 "Thinking of vect as coordinates in space of\n"
488 "dimension n = (vector-length vect), the coordinates\n"
489 "are uniformly distributed over the surface of the\n"
491 #define FUNC_NAME s_scm_random_hollow_sphere_x
493 SCM_VALIDATE_VECTOR_OR_DVECTOR (1, v
);
494 if (SCM_UNBNDP (state
))
495 state
= SCM_VARIABLE_REF (scm_var_random_state
);
496 SCM_VALIDATE_RSTATE (2, state
);
497 scm_random_normal_vector_x (v
, state
);
498 vector_scale (v
, 1 / sqrt (vector_sum_squares (v
)));
499 return SCM_UNSPECIFIED
;
504 SCM_DEFINE (scm_random_normal_vector_x
, "random:normal-vector!", 1, 1, 0,
506 "Fills vect with inexact real random numbers that are\n"
507 "independent and standard normally distributed\n"
508 "(i.e., with mean 0 and variance 1).")
509 #define FUNC_NAME s_scm_random_normal_vector_x
512 SCM_VALIDATE_VECTOR_OR_DVECTOR (1, v
);
513 if (SCM_UNBNDP (state
))
514 state
= SCM_VARIABLE_REF (scm_var_random_state
);
515 SCM_VALIDATE_RSTATE (2, state
);
516 n
= SCM_INUM (scm_uniform_vector_length (v
));
519 SCM_VECTOR_SET (v
, n
, scm_make_real (scm_c_normal01 (SCM_RSTATE (state
))));
522 ((double *) SCM_VELTS (v
))[n
] = scm_c_normal01 (SCM_RSTATE (state
));
523 return SCM_UNSPECIFIED
;
527 #endif /* SCM_HAVE_ARRAYS */
529 SCM_DEFINE (scm_random_exp
, "random:exp", 0, 1, 0,
531 "Return an inexact real in an exponential distribution with mean\n"
532 "1. For an exponential distribution with mean u use (* u\n"
534 #define FUNC_NAME s_scm_random_exp
536 if (SCM_UNBNDP (state
))
537 state
= SCM_VARIABLE_REF (scm_var_random_state
);
538 SCM_VALIDATE_RSTATE (1, state
);
539 return scm_make_real (scm_c_exp1 (SCM_RSTATE (state
)));
547 /* plug in default RNG */
550 sizeof (scm_t_i_rstate
),
551 (unsigned long (*)()) scm_i_uniform32
,
552 (void (*)()) scm_i_init_rstate
,
553 (scm_t_rstate
*(*)()) scm_i_copy_rstate
557 scm_tc16_rstate
= scm_make_smob_type ("random-state", 0);
558 scm_set_smob_free (scm_tc16_rstate
, rstate_free
);
560 for (m
= 1; m
<= 0x100; m
<<= 1)
561 for (i
= m
>> 1; i
< m
; ++i
)
562 scm_masktab
[i
] = m
- 1;
564 #include "libguile/random.x"
566 scm_add_feature ("random");