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)
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.
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
17 * As a special exception, the Free Software Foundation gives permission
18 * for additional uses of the text contained in its release of GUILE.
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.
26 * This exception does not however invalidate any other reasons why
27 * the executable file might be covered by the GNU General Public License.
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.
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. */
41 /* Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */
56 * A plugin interface for RNGs
58 * Using this interface, it is possible for the application to tell
59 * libguile to use a different RNG. This is desirable if it is
60 * necessary to use the same RNG everywhere in the application in
61 * order to prevent interference, if the application uses RNG
62 * hardware, or if the application has special demands on the RNG.
64 * Look in random.h and how the default generator is "plugged in" in
74 * This is the MWC (Multiply With Carry) random number generator
75 * described by George Marsaglia at the Department of Statistics and
76 * Supercomputer Computations Research Institute, The Florida State
77 * University (http://stat.fsu.edu/~geo).
79 * It uses 64 bits, has a period of 4578426017172946943 (4.6e18), and
80 * passes all tests in the DIEHARD test suite
81 * (http://stat.fsu.edu/~geo/diehard.html)
84 #define A 2131995753UL
88 #define LONG32 unsigned short
90 #define LONG32 unsigned int
92 #define LONG64 unsigned long
94 #define LONG32 unsigned long
95 #define LONG64 unsigned long long
98 #if SIZEOF_LONG > 4 || defined (HAVE_LONG_LONGS)
101 scm_i_uniform32 (scm_i_rstate
*state
)
103 LONG64 x
= (LONG64
) A
* state
->w
+ state
->c
;
104 LONG32 w
= x
& 0xffffffffUL
;
112 /* ww This is a portable version of the same RNG without 64 bit
115 * xx It is only intended to provide identical behaviour on
116 * xx platforms without 8 byte longs or long longs until
117 * xx someone has implemented the routine in assembler code.
123 #define L(x) ((x) & 0xffff)
124 #define H(x) ((x) >> 16)
127 scm_i_uniform32 (scm_i_rstate
*state
)
129 LONG32 x1
= L (A
) * L (state
->w
);
130 LONG32 x2
= L (A
) * H (state
->w
);
131 LONG32 x3
= H (A
) * L (state
->w
);
132 LONG32 w
= L (x1
) + L (state
->c
);
133 LONG32 m
= H (x1
) + L (x2
) + L (x3
) + H (state
->c
) + H (w
);
134 LONG32 x4
= H (A
) * H (state
->w
);
135 state
->w
= w
= (L (m
) << 16) + L (w
);
136 state
->c
= H (x2
) + H (x3
) + x4
+ H (m
);
143 scm_i_init_rstate (scm_i_rstate
*state
, char *seed
, int n
)
148 for (i
= 0; i
< n
; ++i
)
152 w
+= seed
[i
] << (8 * m
);
154 c
+= seed
[i
] << (8 * (m
- 4));
156 if ((w
== 0 && c
== 0) || (w
== 0xffffffffUL
&& c
== A
- 1))
163 scm_i_copy_rstate (scm_i_rstate
*state
)
165 scm_rstate
*new_state
= malloc (scm_the_rng
.rstate_size
);
167 scm_wta (SCM_MAKINUM (scm_the_rng
.rstate_size
),
168 (char *) SCM_NALLOC
, "rstate");
169 return memcpy (new_state
, state
, scm_the_rng
.rstate_size
);
174 * Random number library functions
178 scm_i_make_rstate (char *seed
, int n
)
180 scm_rstate
*state
= malloc (scm_the_rng
.rstate_size
);
182 scm_wta (SCM_MAKINUM (scm_the_rng
.rstate_size
),
185 state
->reserved0
= 0;
186 scm_the_rng
.init_rstate (state
, seed
, n
);
191 scm_i_uniform01 (scm_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_i_normal01 (scm_rstate
*state
)
201 if (state
->reserved0
)
203 state
->reserved0
= 0;
204 return state
->reserved1
;
210 r
= sqrt (-2.0 * log (scm_i_uniform01 (state
)));
211 a
= 2.0 * M_PI
* scm_i_uniform01 (state
);
214 state
->reserved1
= r
* cos (a
);
215 state
->reserved0
= 1;
222 scm_i_exp1 (scm_rstate
*state
)
224 return - log (scm_i_uniform01 (state
));
227 unsigned char scm_masktab
[256];
230 scm_i_random (unsigned long m
, scm_rstate
*state
)
232 unsigned int r
, mask
;
236 ? scm_masktab
[m
>> 8] << 8 | 0xff
238 ? scm_masktab
[m
>> 16] << 16 | 0xffff
239 : scm_masktab
[m
>> 24] << 24 | 0xffffff)));
240 while ((r
= scm_the_rng
.random_bits (state
) & mask
) >= m
);
245 scm_i_random_bignum (SCM m
, scm_rstate
*state
)
249 LONG32
*bits
, mask
, w
;
250 nd
= SCM_NUMDIGS (m
);
251 /* calculate mask for most significant digit */
256 /* fix most significant 16 bits */
257 unsigned short s
= SCM_BDIGITS (m
)[nd
- 1];
258 mask
= s
< 0x100 ? scm_masktab
[s
] : scm_masktab
[s
>> 8] << 8 | 0xff;
263 /* fix most significant 32 bits */
265 w
= SCM_BDIGITS (m
)[nd
- 1] << 16 | SCM_BDIGITS (m
)[nd
- 2];
267 w
= SCM_BDIGITS (m
)[nd
- 1];
272 : scm_masktab
[w
>> 8] << 8 | 0xff)
274 ? scm_masktab
[w
>> 16] << 16 | 0xffff
275 : scm_masktab
[w
>> 24] << 24 | 0xffffff));
277 b
= scm_mkbig (nd
, 0);
278 bits
= (LONG32
*) SCM_BDIGITS (b
);
282 /* treat most significant digit specially */
287 ((SCM_BIGDIG
*) bits
)[i
- 1] = scm_the_rng
.random_bits (state
) & mask
;
293 /* fix most significant 32 bits */
295 w
= scm_the_rng
.random_bits (state
) & mask
;
296 ((SCM_BIGDIG
*) bits
)[i
- 2] = w
& 0xffff;
297 ((SCM_BIGDIG
*) bits
)[i
- 1] = w
>> 16;
301 bits
[--i
] = scm_the_rng
.random_bits (state
) & mask
;
304 /* now fill up the rest of the bignum */
306 bits
[--i
] = scm_the_rng
.random_bits (state
);
310 } while (scm_bigcomp (b
, m
) <= 0);
315 * Scheme level representation of random states.
318 long scm_tc16_rstate
;
321 make_rstate (scm_rstate
*state
)
326 SCM_SETCDR (cell
, state
);
327 SCM_SETCAR (cell
, scm_tc16_rstate
);
333 print_rstate (SCM rstate
, SCM port
, scm_print_state
*pstate
)
335 scm_puts ("#<random-state ", port
);
336 scm_intprint ((long) SCM_RSTATE (rstate
), 16, port
);
337 scm_putc ('>', port
);
342 free_rstate (SCM rstate
)
344 free (SCM_RSTATE (rstate
));
345 return scm_the_rng
.rstate_size
;
348 static scm_smobfuns rstate_smob
= { 0, free_rstate
, print_rstate
, 0};
351 * Scheme level interface.
354 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")));
356 SCM_PROC (s_random
, "random", 1, 1, 0, scm_random
);
359 scm_random (SCM n
, SCM state
)
361 if (SCM_UNBNDP (state
))
362 state
= SCM_CDR (scm_var_random_state
);
363 SCM_ASSERT (SCM_NIMP (state
) && SCM_RSTATEP (state
),
364 state
, SCM_ARG2
, s_random
);
367 unsigned long m
= SCM_INUM (n
);
368 SCM_ASSERT (m
> 0, n
, SCM_ARG1
, s_random
);
369 return SCM_MAKINUM (scm_i_random (m
, SCM_RSTATE (state
)));
371 SCM_ASSERT (SCM_NIMP (n
), n
, SCM_ARG1
, s_random
);
373 return scm_makdbl (SCM_REALPART (n
) * scm_i_uniform01 (SCM_RSTATE (state
)),
375 SCM_ASSERT (SCM_TYP16 (n
) == scm_tc16_bigpos
, n
, SCM_ARG1
, s_random
);
376 return scm_i_random_bignum (n
, SCM_RSTATE (state
));
379 SCM_PROC (s_copy_random_state
, "copy-random-state", 0, 1, 0, scm_copy_random_state
);
382 scm_copy_random_state (SCM state
)
384 if (SCM_UNBNDP (state
))
385 state
= SCM_CDR (scm_var_random_state
);
386 SCM_ASSERT (SCM_NIMP (state
) && SCM_RSTATEP (state
),
389 s_copy_random_state
);
390 return make_rstate (scm_the_rng
.copy_rstate (SCM_RSTATE (state
)));
393 SCM_PROC (s_seed_to_random_state
, "seed->random-state", 1, 0, 0, scm_seed_to_random_state
);
396 scm_seed_to_random_state (SCM seed
)
398 if (SCM_NUMBERP (seed
))
399 seed
= scm_number_to_string (seed
, SCM_UNDEFINED
);
400 SCM_ASSERT (SCM_NIMP (seed
) && SCM_STRINGP (seed
),
403 s_seed_to_random_state
);
404 return make_rstate (scm_i_make_rstate (SCM_ROCHARS (seed
),
408 SCM_PROC (s_random_uniform
, "random:uniform", 0, 1, 0, scm_random_uniform
);
411 scm_random_uniform (SCM state
)
413 if (SCM_UNBNDP (state
))
414 state
= SCM_CDR (scm_var_random_state
);
415 SCM_ASSERT (SCM_NIMP (state
) && SCM_RSTATEP (state
),
419 return scm_makdbl (scm_i_uniform01 (SCM_RSTATE (state
)), 0.0);
423 vector_scale (SCM v
, double c
)
425 int n
= SCM_LENGTH (v
);
428 SCM_REAL (SCM_VELTS (v
)[n
]) *= c
;
431 ((double *) SCM_VELTS (v
))[n
] *= c
;
435 vector_sum_squares (SCM v
)
438 int n
= SCM_LENGTH (v
);
442 x
= SCM_REAL (SCM_VELTS (v
)[n
]);
448 x
= ((double *) SCM_VELTS (v
))[n
];
454 SCM_PROC (s_random_solid_sphere_x
, "random:solid-sphere!", 1, 1, 0, scm_random_solid_sphere_x
);
456 /* For the uniform distribution on the solid sphere, note that in
457 * this distribution the length r of the vector has cumulative
458 * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be
459 * generated as r=u^(1/n).
462 scm_random_solid_sphere_x (SCM v
, SCM state
)
464 SCM_ASSERT (SCM_NIMP (v
)
465 && (SCM_VECTORP (v
) || SCM_TYP7 (v
) == scm_tc7_dvect
),
466 v
, SCM_ARG1
, s_random_solid_sphere_x
);
467 if (SCM_UNBNDP (state
))
468 state
= SCM_CDR (scm_var_random_state
);
469 SCM_ASSERT (SCM_NIMP (state
) && SCM_RSTATEP (state
),
472 s_random_solid_sphere_x
);
473 scm_random_normal_vector_x (v
, state
);
475 pow (scm_i_uniform01 (SCM_RSTATE (state
)),
476 1.0 / SCM_LENGTH (v
))
477 / sqrt (vector_sum_squares (v
)));
478 return SCM_UNSPECIFIED
;
481 SCM_PROC (s_random_hollow_sphere_x
, "random:hollow-sphere!", 1, 1, 0, scm_random_hollow_sphere_x
);
484 scm_random_hollow_sphere_x (SCM v
, SCM state
)
486 SCM_ASSERT (SCM_NIMP (v
)
487 && (SCM_VECTORP (v
) || SCM_TYP7 (v
) == scm_tc7_dvect
),
488 v
, SCM_ARG1
, s_random_solid_sphere_x
);
489 if (SCM_UNBNDP (state
))
490 state
= SCM_CDR (scm_var_random_state
);
491 SCM_ASSERT (SCM_NIMP (state
) && SCM_RSTATEP (state
),
494 s_random_hollow_sphere_x
);
495 scm_random_normal_vector_x (v
, state
);
496 vector_scale (v
, 1 / sqrt (vector_sum_squares (v
)));
497 return SCM_UNSPECIFIED
;
500 SCM_PROC (s_random_normal
, "random:normal", 0, 1, 0, scm_random_normal
);
503 scm_random_normal (SCM state
)
505 if (SCM_UNBNDP (state
))
506 state
= SCM_CDR (scm_var_random_state
);
507 SCM_ASSERT (SCM_NIMP (state
) && SCM_RSTATEP (state
),
511 return scm_makdbl (scm_i_normal01 (SCM_RSTATE (state
)), 0.0);
514 SCM_PROC (s_random_normal_vector_x
, "random:normal-vector!", 1, 1, 0, scm_random_normal_vector_x
);
517 scm_random_normal_vector_x (SCM v
, SCM state
)
520 SCM_ASSERT (SCM_NIMP (v
)
521 && (SCM_VECTORP (v
) || SCM_TYP7 (v
) == scm_tc7_dvect
),
522 v
, SCM_ARG1
, s_random_solid_sphere_x
);
523 if (SCM_UNBNDP (state
))
524 state
= SCM_CDR (scm_var_random_state
);
525 SCM_ASSERT (SCM_NIMP (state
) && SCM_RSTATEP (state
),
528 s_random_normal_vector_x
);
532 SCM_VELTS (v
)[n
] = scm_makdbl (scm_i_normal01 (SCM_RSTATE (state
)), 0.0);
535 ((double *) SCM_VELTS (v
))[n
] = scm_i_normal01 (SCM_RSTATE (state
));
536 return SCM_UNSPECIFIED
;
539 SCM_PROC (s_random_exp
, "random:exp", 0, 1, 0, scm_random_exp
);
542 scm_random_exp (SCM state
)
544 if (SCM_UNBNDP (state
))
545 state
= SCM_CDR (scm_var_random_state
);
546 SCM_ASSERT (SCM_NIMP (state
) && SCM_RSTATEP (state
),
550 return scm_makdbl (scm_i_exp1 (SCM_RSTATE (state
)), 0.0);
557 /* plug in default RNG */
560 sizeof (scm_i_rstate
),
561 (unsigned long (*)()) scm_i_uniform32
,
562 (void (*)()) scm_i_init_rstate
,
563 (scm_rstate
*(*)()) scm_i_copy_rstate
567 scm_tc16_rstate
= scm_newsmob (&rstate_smob
);
569 for (m
= 1; m
<= 0x100; m
<<= 1)
570 for (i
= m
>> 1; i
< m
; ++i
)
571 scm_masktab
[i
] = m
- 1;
575 /* Check that the assumptions about bits per bignum digit are correct. */
581 if (m
!= SCM_BITSPERDIG
)
583 fprintf (stderr
, "Internal inconsistency: Confused about bignum digit size in random.c\n");
587 scm_add_feature ("random");