(scm_c_uniform_vector_element_size,
[bpt/guile.git] / libguile / random.c
CommitLineData
ea7f6344 1/* Copyright (C) 1999,2000,2001, 2003 Free Software Foundation, Inc.
73be1d9e
MV
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.
e7a72986 6 *
73be1d9e 7 * This library is distributed in the hope that it will be useful,
e7a72986 8 * but WITHOUT ANY WARRANTY; without even the implied warranty of
73be1d9e
MV
9 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
10 * Lesser General Public License for more details.
e7a72986 11 *
73be1d9e
MV
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
15 */
e7a72986 16
1bbd0b84
GB
17
18
e7a72986
MD
19/* Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */
20
2fb8013c
RB
21#if HAVE_CONFIG_H
22# include <config.h>
23#endif
24
a0599745 25#include "libguile/_scm.h"
e7a72986 26
8db9cc6c 27#include <gmp.h>
7f146094 28#include <stdio.h>
e7a72986 29#include <math.h>
f34d19c7 30#include <string.h>
a0599745
MD
31#include "libguile/smob.h"
32#include "libguile/numbers.h"
33#include "libguile/feature.h"
34#include "libguile/strings.h"
e9bfab50 35#include "libguile/unif.h"
46d25cff 36#include "libguile/srfi-4.h"
a0599745 37#include "libguile/vectors.h"
e7a72986 38
a0599745
MD
39#include "libguile/validate.h"
40#include "libguile/random.h"
e7a72986
MD
41
42\f
43/*
44 * A plugin interface for RNGs
45 *
46 * Using this interface, it is possible for the application to tell
47 * libguile to use a different RNG. This is desirable if it is
48 * necessary to use the same RNG everywhere in the application in
49 * order to prevent interference, if the application uses RNG
50 * hardware, or if the application has special demands on the RNG.
51 *
52 * Look in random.h and how the default generator is "plugged in" in
53 * scm_init_random().
54 */
55
92c2555f 56scm_t_rng scm_the_rng;
e7a72986
MD
57
58\f
59/*
60 * The prepackaged RNG
61 *
62 * This is the MWC (Multiply With Carry) random number generator
63 * described by George Marsaglia at the Department of Statistics and
64 * Supercomputer Computations Research Institute, The Florida State
65 * University (http://stat.fsu.edu/~geo).
66 *
67 * It uses 64 bits, has a period of 4578426017172946943 (4.6e18), and
68 * passes all tests in the DIEHARD test suite
69 * (http://stat.fsu.edu/~geo/diehard.html)
70 */
71
72#define A 2131995753UL
73
82893676
MG
74#ifndef M_PI
75#define M_PI 3.14159265359
76#endif
77
2fb8013c 78#ifdef SCM_HAVE_T_INT64
e7a72986
MD
79
80unsigned long
92c2555f 81scm_i_uniform32 (scm_t_i_rstate *state)
e7a72986 82{
2fb8013c
RB
83 scm_t_int64 x = (scm_t_int64) A * state->w + state->c;
84 scm_t_int32 w = x & 0xffffffffUL;
e7a72986
MD
85 state->w = w;
86 state->c = x >> 32L;
87 return w;
88}
89
90#else
91
92/* ww This is a portable version of the same RNG without 64 bit
93 * * aa arithmetic.
94 * ----
95 * xx It is only intended to provide identical behaviour on
96 * xx platforms without 8 byte longs or long longs until
97 * xx someone has implemented the routine in assembler code.
98 * xxcc
99 * ----
100 * ccww
101 */
102
103#define L(x) ((x) & 0xffff)
104#define H(x) ((x) >> 16)
105
106unsigned long
92c2555f 107scm_i_uniform32 (scm_t_i_rstate *state)
e7a72986 108{
2fb8013c
RB
109 scm_t_int32 x1 = L (A) * L (state->w);
110 scm_t_int32 x2 = L (A) * H (state->w);
111 scm_t_int32 x3 = H (A) * L (state->w);
112 scm_t_int32 w = L (x1) + L (state->c);
113 scm_t_int32 m = H (x1) + L (x2) + L (x3) + H (state->c) + H (w);
114 scm_t_int32 x4 = H (A) * H (state->w);
e7a72986
MD
115 state->w = w = (L (m) << 16) + L (w);
116 state->c = H (x2) + H (x3) + x4 + H (m);
117 return w;
118}
119
120#endif
121
122void
cc95e00a 123scm_i_init_rstate (scm_t_i_rstate *state, const char *seed, int n)
e7a72986 124{
2fb8013c
RB
125 scm_t_int32 w = 0L;
126 scm_t_int32 c = 0L;
e7a72986
MD
127 int i, m;
128 for (i = 0; i < n; ++i)
129 {
130 m = i % 8;
131 if (m < 4)
132 w += seed[i] << (8 * m);
133 else
134 c += seed[i] << (8 * (m - 4));
135 }
136 if ((w == 0 && c == 0) || (w == 0xffffffffUL && c == A - 1))
137 ++c;
138 state->w = w;
139 state->c = c;
140}
141
92c2555f
MV
142scm_t_i_rstate *
143scm_i_copy_rstate (scm_t_i_rstate *state)
e7a72986 144{
67329a9e 145 scm_t_rstate *new_state = scm_malloc (scm_the_rng.rstate_size);
e7a72986 146 if (new_state == 0)
2500356c 147 scm_memory_error ("rstate");
e7a72986
MD
148 return memcpy (new_state, state, scm_the_rng.rstate_size);
149}
150
151\f
152/*
153 * Random number library functions
154 */
155
92c2555f 156scm_t_rstate *
cc95e00a 157scm_c_make_rstate (const char *seed, int n)
5ee11b7c 158{
67329a9e 159 scm_t_rstate *state = scm_malloc (scm_the_rng.rstate_size);
5ee11b7c 160 if (state == 0)
2500356c 161 scm_memory_error ("rstate");
5ee11b7c
MD
162 state->reserved0 = 0;
163 scm_the_rng.init_rstate (state, seed, n);
164 return state;
165}
166
2ade72d7 167
92c2555f 168scm_t_rstate *
9b741bb6 169scm_c_default_rstate ()
2ade72d7 170#define FUNC_NAME "scm_c_default_rstate"
9b741bb6 171{
c5f268f8 172 SCM state = SCM_VARIABLE_REF (scm_var_random_state);
2ade72d7
DH
173 if (!SCM_RSTATEP (state))
174 SCM_MISC_ERROR ("*random-state* contains bogus random state", SCM_EOL);
9b741bb6
MD
175 return SCM_RSTATE (state);
176}
2ade72d7
DH
177#undef FUNC_NAME
178
9b741bb6 179
e7a72986 180inline double
92c2555f 181scm_c_uniform01 (scm_t_rstate *state)
e7a72986 182{
5a92ddfd 183 double x = (double) scm_the_rng.random_bits (state) / (double) 0xffffffffUL;
e7a72986 184 return ((x + (double) scm_the_rng.random_bits (state))
5a92ddfd 185 / (double) 0xffffffffUL);
e7a72986
MD
186}
187
188double
92c2555f 189scm_c_normal01 (scm_t_rstate *state)
e7a72986
MD
190{
191 if (state->reserved0)
192 {
193 state->reserved0 = 0;
194 return state->reserved1;
195 }
196 else
197 {
198 double r, a, n;
e7a72986 199
9b741bb6
MD
200 r = sqrt (-2.0 * log (scm_c_uniform01 (state)));
201 a = 2.0 * M_PI * scm_c_uniform01 (state);
e7a72986
MD
202
203 n = r * sin (a);
204 state->reserved1 = r * cos (a);
5a92ddfd 205 state->reserved0 = 1;
e7a72986
MD
206
207 return n;
208 }
209}
210
211double
92c2555f 212scm_c_exp1 (scm_t_rstate *state)
e7a72986 213{
9b741bb6 214 return - log (scm_c_uniform01 (state));
e7a72986
MD
215}
216
217unsigned char scm_masktab[256];
218
219unsigned long
92c2555f 220scm_c_random (scm_t_rstate *state, unsigned long m)
e7a72986
MD
221{
222 unsigned int r, mask;
223 mask = (m < 0x100
224 ? scm_masktab[m]
225 : (m < 0x10000
5a92ddfd 226 ? scm_masktab[m >> 8] << 8 | 0xff
e7a72986 227 : (m < 0x1000000
5a92ddfd
MD
228 ? scm_masktab[m >> 16] << 16 | 0xffff
229 : scm_masktab[m >> 24] << 24 | 0xffffff)));
e7a72986
MD
230 while ((r = scm_the_rng.random_bits (state) & mask) >= m);
231 return r;
232}
233
ea7f6344
MD
234/*
235 SCM scm_c_random_bignum (scm_t_rstate *state, SCM m)
236
237 Takes a random state (source of random bits) and a bignum m.
238 Returns a bignum b, 0 <= b < m.
239
240 It does this by allocating a bignum b with as many base 65536 digits
241 as m, filling b with random bits (in 32 bit chunks) up to the most
242 significant 1 in m, and, finally checking if the resultant b is too
243 large (>= m). If too large, we simply repeat the process again. (It
244 is important to throw away all generated random bits if b >= m,
245 otherwise we'll end up with a distorted distribution.)
246
247*/
248
e7a72986 249SCM
92c2555f 250scm_c_random_bignum (scm_t_rstate *state, SCM m)
e7a72986 251{
969d3bd0 252 SCM result = scm_i_mkbig ();
969d3bd0
RB
253 const size_t m_bits = mpz_sizeinbase (SCM_I_BIG_MPZ (m), 2);
254 /* how many bits would only partially fill the last unsigned long? */
0d79003d 255 const size_t end_bits = m_bits % (sizeof (unsigned long) * SCM_CHAR_BIT);
969d3bd0 256 unsigned long *random_chunks = NULL;
0d79003d
RB
257 const unsigned long num_full_chunks =
258 m_bits / (sizeof (unsigned long) * SCM_CHAR_BIT);
969d3bd0
RB
259 const unsigned long num_chunks = num_full_chunks + ((end_bits) ? 1 : 0);
260
261 /* we know the result will be this big */
262 mpz_realloc2 (SCM_I_BIG_MPZ (result), m_bits);
263
264 random_chunks =
265 (unsigned long *) scm_gc_calloc (num_chunks * sizeof (unsigned long),
266 "random bignum chunks");
267
372691d8 268 do
e7a72986 269 {
969d3bd0
RB
270 unsigned long *current_chunk = random_chunks + (num_chunks - 1);
271 unsigned long chunks_left = num_chunks;
969d3bd0
RB
272
273 mpz_set_ui (SCM_I_BIG_MPZ (result), 0);
274
275 if (end_bits)
276 {
277 /* generate a mask with ones in the end_bits position, i.e. if
278 end_bits is 3, then we'd have a mask of ...0000000111 */
279 const unsigned long rndbits = scm_the_rng.random_bits (state);
0d79003d 280 int rshift = (sizeof (unsigned long) * SCM_CHAR_BIT) - end_bits;
969d3bd0
RB
281 unsigned long mask = ((unsigned long) ULONG_MAX) >> rshift;
282 unsigned long highest_bits = rndbits & mask;
283 *current_chunk-- = highest_bits;
284 chunks_left--;
285 }
286
287 while (chunks_left)
288 {
289 /* now fill in the remaining unsigned long sized chunks */
290 *current_chunk-- = scm_the_rng.random_bits (state);
291 chunks_left--;
292 }
293 mpz_import (SCM_I_BIG_MPZ (result),
294 num_chunks,
295 -1,
296 sizeof (unsigned long),
297 0,
298 0,
299 random_chunks);
372691d8
MD
300 /* if result >= m, regenerate it (it is important to regenerate
301 all bits in order not to get a distorted distribution) */
302 } while (mpz_cmp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (m)) >= 0);
969d3bd0
RB
303 scm_gc_free (random_chunks,
304 num_chunks * sizeof (unsigned long),
305 "random bignum chunks");
8db9cc6c 306 return scm_i_normbig (result);
e7a72986
MD
307}
308
309/*
310 * Scheme level representation of random states.
311 */
312
92c2555f 313scm_t_bits scm_tc16_rstate;
e7a72986
MD
314
315static SCM
92c2555f 316make_rstate (scm_t_rstate *state)
e7a72986 317{
23a62151 318 SCM_RETURN_NEWSMOB (scm_tc16_rstate, state);
e7a72986
MD
319}
320
1be6b49c 321static size_t
e841c3e0 322rstate_free (SCM rstate)
e7a72986
MD
323{
324 free (SCM_RSTATE (rstate));
dfd71aba 325 return 0;
e7a72986
MD
326}
327
e7a72986
MD
328/*
329 * Scheme level interface.
330 */
331
cc95e00a 332SCM_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")));
e7a72986 333
a1ec6916 334SCM_DEFINE (scm_random, "random", 1, 1, 0,
1bbd0b84 335 (SCM n, SCM state),
34d19ef6 336 "Return a number in [0, N).\n"
d928e0b4 337 "\n"
9401323e
NJ
338 "Accepts a positive integer or real n and returns a\n"
339 "number of the same type between zero (inclusive) and\n"
340 "N (exclusive). The values returned have a uniform\n"
d928e0b4
GB
341 "distribution.\n"
342 "\n"
3b644514
MG
343 "The optional argument @var{state} must be of the type produced\n"
344 "by @code{seed->random-state}. It defaults to the value of the\n"
345 "variable @var{*random-state*}. This object is used to maintain\n"
346 "the state of the pseudo-random-number generator and is altered\n"
347 "as a side effect of the random operation.")
1bbd0b84 348#define FUNC_NAME s_scm_random
e7a72986
MD
349{
350 if (SCM_UNBNDP (state))
86d31dfe 351 state = SCM_VARIABLE_REF (scm_var_random_state);
34d19ef6 352 SCM_VALIDATE_RSTATE (2, state);
e11e83f3 353 if (SCM_I_INUMP (n))
e7a72986 354 {
e11e83f3 355 unsigned long m = SCM_I_INUM (n);
34d19ef6 356 SCM_ASSERT_RANGE (1, n, m > 0);
e11e83f3 357 return scm_from_ulong (scm_c_random (SCM_RSTATE (state), m));
e7a72986 358 }
34d19ef6 359 SCM_VALIDATE_NIM (1, n);
e7a72986 360 if (SCM_REALP (n))
d9a67fc4
MV
361 return scm_from_double (SCM_REAL_VALUE (n)
362 * scm_c_uniform01 (SCM_RSTATE (state)));
969d3bd0 363
a55c2b68
MV
364 if (!SCM_BIGP (n))
365 SCM_WRONG_TYPE_ARG (1, n);
9b741bb6 366 return scm_c_random_bignum (SCM_RSTATE (state), n);
e7a72986 367}
1bbd0b84 368#undef FUNC_NAME
e7a72986 369
a1ec6916 370SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0,
1bbd0b84 371 (SCM state),
3b644514 372 "Return a copy of the random state @var{state}.")
1bbd0b84 373#define FUNC_NAME s_scm_copy_random_state
e7a72986
MD
374{
375 if (SCM_UNBNDP (state))
86d31dfe 376 state = SCM_VARIABLE_REF (scm_var_random_state);
34d19ef6 377 SCM_VALIDATE_RSTATE (1, state);
e7a72986
MD
378 return make_rstate (scm_the_rng.copy_rstate (SCM_RSTATE (state)));
379}
1bbd0b84 380#undef FUNC_NAME
e7a72986 381
a1ec6916 382SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0,
1bbd0b84 383 (SCM seed),
3b644514 384 "Return a new random state using @var{seed}.")
1bbd0b84 385#define FUNC_NAME s_scm_seed_to_random_state
5ee11b7c 386{
8824ac88 387 SCM res;
5ee11b7c
MD
388 if (SCM_NUMBERP (seed))
389 seed = scm_number_to_string (seed, SCM_UNDEFINED);
34d19ef6 390 SCM_VALIDATE_STRING (1, seed);
cc95e00a
MV
391 res = make_rstate (scm_c_make_rstate (scm_i_string_chars (seed),
392 scm_i_string_length (seed)));
8824ac88
MV
393 scm_remember_upto_here_1 (seed);
394 return res;
395
5ee11b7c 396}
1bbd0b84 397#undef FUNC_NAME
5ee11b7c 398
a1ec6916 399SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0,
1bbd0b84 400 (SCM state),
1e6808ea
MG
401 "Return a uniformly distributed inexact real random number in\n"
402 "[0,1).")
1bbd0b84 403#define FUNC_NAME s_scm_random_uniform
e7a72986
MD
404{
405 if (SCM_UNBNDP (state))
86d31dfe 406 state = SCM_VARIABLE_REF (scm_var_random_state);
34d19ef6 407 SCM_VALIDATE_RSTATE (1, state);
d9a67fc4 408 return scm_from_double (scm_c_uniform01 (SCM_RSTATE (state)));
e7a72986 409}
1bbd0b84 410#undef FUNC_NAME
e7a72986 411
a1ec6916 412SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0,
1bbd0b84 413 (SCM state),
1e6808ea
MG
414 "Return an inexact real in a normal distribution. The\n"
415 "distribution used has mean 0 and standard deviation 1. For a\n"
416 "normal distribution with mean m and standard deviation d use\n"
417 "@code{(+ m (* d (random:normal)))}.")
1bbd0b84 418#define FUNC_NAME s_scm_random_normal
afe5177e
GH
419{
420 if (SCM_UNBNDP (state))
86d31dfe 421 state = SCM_VARIABLE_REF (scm_var_random_state);
34d19ef6 422 SCM_VALIDATE_RSTATE (1, state);
d9a67fc4 423 return scm_from_double (scm_c_normal01 (SCM_RSTATE (state)));
afe5177e 424}
1bbd0b84 425#undef FUNC_NAME
afe5177e 426
e7a72986 427static void
46d25cff 428vector_scale_x (SCM v, double c)
e7a72986 429{
46d25cff 430 size_t n;
4057a3e0 431 if (scm_is_simple_vector (v))
46d25cff 432 {
4057a3e0 433 n = SCM_SIMPLE_VECTOR_LENGTH (v);
46d25cff 434 while (n-- > 0)
4057a3e0 435 SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n)) *= c;
46d25cff 436 }
e7a72986 437 else
46d25cff
MV
438 {
439 /* must be a f64vector. */
4057a3e0
MV
440 scm_t_array_handle handle;
441 size_t i, len;
442 ssize_t inc;
443 double *elts;
444
445 elts = scm_f64vector_writable_elements (v, &handle, &len, &inc);
446
447 for (i = 0; i < len; i++, elts += inc)
448 *elts *= c;
46d25cff 449 }
e7a72986
MD
450}
451
452static double
453vector_sum_squares (SCM v)
454{
455 double x, sum = 0.0;
46d25cff 456 size_t n;
4057a3e0 457 if (scm_is_simple_vector (v))
46d25cff 458 {
4057a3e0 459 n = SCM_SIMPLE_VECTOR_LENGTH (v);
46d25cff
MV
460 while (n-- > 0)
461 {
4057a3e0 462 x = SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n));
46d25cff
MV
463 sum += x * x;
464 }
465 }
e7a72986 466 else
46d25cff
MV
467 {
468 /* must be a f64vector. */
4057a3e0
MV
469 scm_t_array_handle handle;
470 size_t i, len;
471 ssize_t inc;
472 const double *elts;
473
474 elts = scm_f64vector_elements (v, &handle, &len, &inc);
475
476 for (i = 0; i < len; i++, elts += inc)
46d25cff 477 {
4057a3e0 478 x = *elts;
46d25cff
MV
479 sum += x * x;
480 }
4057a3e0 481
46d25cff 482 }
e7a72986
MD
483 return sum;
484}
485
e7a72986
MD
486/* For the uniform distribution on the solid sphere, note that in
487 * this distribution the length r of the vector has cumulative
488 * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be
489 * generated as r=u^(1/n).
490 */
a1ec6916 491SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0,
1bbd0b84 492 (SCM v, SCM state),
d928e0b4
GB
493 "Fills vect with inexact real random numbers\n"
494 "the sum of whose squares is less than 1.0.\n"
9401323e
NJ
495 "Thinking of vect as coordinates in space of\n"
496 "dimension n = (vector-length vect), the coordinates\n"
72dd0a03 497 "are uniformly distributed within the unit n-sphere.\n"
64ba8e85 498 "The sum of the squares of the numbers is returned.")
1bbd0b84 499#define FUNC_NAME s_scm_random_solid_sphere_x
e7a72986 500{
e7a72986 501 if (SCM_UNBNDP (state))
86d31dfe 502 state = SCM_VARIABLE_REF (scm_var_random_state);
34d19ef6 503 SCM_VALIDATE_RSTATE (2, state);
e7a72986 504 scm_random_normal_vector_x (v, state);
46d25cff
MV
505 vector_scale_x (v,
506 pow (scm_c_uniform01 (SCM_RSTATE (state)),
507 1.0 / scm_to_int (scm_uniform_vector_length (v)))
508 / sqrt (vector_sum_squares (v)));
e7a72986
MD
509 return SCM_UNSPECIFIED;
510}
1bbd0b84 511#undef FUNC_NAME
e7a72986 512
a1ec6916 513SCM_DEFINE (scm_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0,
1bbd0b84 514 (SCM v, SCM state),
d928e0b4
GB
515 "Fills vect with inexact real random numbers\n"
516 "the sum of whose squares is equal to 1.0.\n"
9401323e 517 "Thinking of vect as coordinates in space of\n"
d928e0b4 518 "dimension n = (vector-length vect), the coordinates\n"
9401323e 519 "are uniformly distributed over the surface of the\n"
72dd0a03 520 "unit n-sphere.")
1bbd0b84 521#define FUNC_NAME s_scm_random_hollow_sphere_x
e7a72986 522{
e7a72986 523 if (SCM_UNBNDP (state))
86d31dfe 524 state = SCM_VARIABLE_REF (scm_var_random_state);
34d19ef6 525 SCM_VALIDATE_RSTATE (2, state);
e7a72986 526 scm_random_normal_vector_x (v, state);
46d25cff 527 vector_scale_x (v, 1 / sqrt (vector_sum_squares (v)));
e7a72986
MD
528 return SCM_UNSPECIFIED;
529}
1bbd0b84 530#undef FUNC_NAME
e7a72986 531
1bbd0b84 532
a1ec6916 533SCM_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 539{
4057a3e0
MV
540 long i;
541 scm_t_array_handle handle;
542 scm_t_array_dim *dim;
46d25cff 543
e7a72986 544 if (SCM_UNBNDP (state))
86d31dfe 545 state = SCM_VARIABLE_REF (scm_var_random_state);
34d19ef6 546 SCM_VALIDATE_RSTATE (2, state);
4057a3e0
MV
547
548 scm_vector_get_handle (v, &handle);
549 dim = scm_array_handle_dims (&handle);
550
551 if (scm_is_vector (v))
46d25cff 552 {
4057a3e0
MV
553 SCM *elts = scm_array_handle_writable_elements (&handle);
554 for (i = dim->lbnd; i <= dim->ubnd; i++, elts += dim->inc)
555 *elts = scm_from_double (scm_c_normal01 (SCM_RSTATE (state)));
46d25cff 556 }
e7a72986 557 else
46d25cff
MV
558 {
559 /* must be a f64vector. */
4057a3e0
MV
560 double *elts = scm_array_handle_f64_writable_elements (&handle);
561 for (i = dim->lbnd; i <= dim->ubnd; i++, elts += dim->inc)
562 *elts = scm_c_normal01 (SCM_RSTATE (state));
46d25cff 563 }
4057a3e0 564
e7a72986
MD
565 return SCM_UNSPECIFIED;
566}
1bbd0b84 567#undef FUNC_NAME
e7a72986 568
a1ec6916 569SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0,
1bbd0b84 570 (SCM state),
1e6808ea
MG
571 "Return an inexact real in an exponential distribution with mean\n"
572 "1. For an exponential distribution with mean u use (* u\n"
573 "(random:exp)).")
1bbd0b84 574#define FUNC_NAME s_scm_random_exp
e7a72986
MD
575{
576 if (SCM_UNBNDP (state))
86d31dfe 577 state = SCM_VARIABLE_REF (scm_var_random_state);
34d19ef6 578 SCM_VALIDATE_RSTATE (1, state);
d9a67fc4 579 return scm_from_double (scm_c_exp1 (SCM_RSTATE (state)));
e7a72986 580}
1bbd0b84 581#undef FUNC_NAME
e7a72986
MD
582
583void
584scm_init_random ()
585{
586 int i, m;
587 /* plug in default RNG */
92c2555f 588 scm_t_rng rng =
e7a72986 589 {
92c2555f 590 sizeof (scm_t_i_rstate),
e7a72986
MD
591 (unsigned long (*)()) scm_i_uniform32,
592 (void (*)()) scm_i_init_rstate,
92c2555f 593 (scm_t_rstate *(*)()) scm_i_copy_rstate
e7a72986
MD
594 };
595 scm_the_rng = rng;
596
e841c3e0
KN
597 scm_tc16_rstate = scm_make_smob_type ("random-state", 0);
598 scm_set_smob_free (scm_tc16_rstate, rstate_free);
e7a72986
MD
599
600 for (m = 1; m <= 0x100; m <<= 1)
601 for (i = m >> 1; i < m; ++i)
602 scm_masktab[i] = m - 1;
603
a0599745 604#include "libguile/random.x"
e7a72986
MD
605
606 scm_add_feature ("random");
607}
89e00824
ML
608
609/*
610 Local Variables:
611 c-file-style: "gnu"
612 End:
613*/