From: Mikael Djurfeldt Date: Sun, 10 Jan 1999 07:57:27 +0000 (+0000) Subject: * random.c, random.h: New files: Random number support. X-Git-Url: https://git.hcoop.net/bpt/guile.git/commitdiff_plain/e7a72986684114de55c73799e99c4fc9897f1cb8 * random.c, random.h: New files: Random number support. Interface same as slib's. --- diff --git a/libguile/random.c b/libguile/random.c new file mode 100644 index 000000000..83496aa76 --- /dev/null +++ b/libguile/random.c @@ -0,0 +1,577 @@ +/* Copyright (C) 1999 Free Software Foundation, Inc. + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this software; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 59 Temple Place, Suite 330, + * Boston, MA 02111-1307 USA + * + * As a special exception, the Free Software Foundation gives permission + * for additional uses of the text contained in its release of GUILE. + * + * The exception is that, if you link the GUILE library with other files + * to produce an executable, this does not by itself cause the + * resulting executable to be covered by the GNU General Public License. + * Your use of that executable is in no way restricted on account of + * linking the GUILE library code into it. + * + * This exception does not however invalidate any other reasons why + * the executable file might be covered by the GNU General Public License. + * + * This exception applies only to the code released by the + * Free Software Foundation under the name GUILE. If you copy + * code from other Free Software Foundation releases into a copy of + * GUILE, as the General Public License permits, the exception does + * not apply to the code that you add in this way. To avoid misleading + * anyone as to the status of such modified files, you must delete + * this exception notice from them. + * + * If you write modifications of your own for GUILE, it is your choice + * whether to permit this exception to apply to your modifications. + * If you do not wish that, delete this exception notice. */ + +/* Author: Mikael Djurfeldt */ + +/* We need this to get the definitions for HAVE_ALLOCA_H, etc. */ +#include "scmconfig.h" + +/* AIX requires this to be the first thing in the file. The #pragma + directive is indented so pre-ANSI compilers will ignore it, rather + than choke on it. */ +#ifndef __GNUC__ +# if HAVE_ALLOCA_H +# include +# else +# ifdef _AIX + #pragma alloca +# else +# ifndef alloca /* predefined by HP cc +Olibcalls */ +char *alloca (); +# endif +# endif +# endif +#endif + +#include "_scm.h" + +#include +#include "genio.h" +#include "smob.h" +#include "numbers.h" +#include "feature.h" + +#include "random.h" + + +/* + * A plugin interface for RNGs + * + * Using this interface, it is possible for the application to tell + * libguile to use a different RNG. This is desirable if it is + * necessary to use the same RNG everywhere in the application in + * order to prevent interference, if the application uses RNG + * hardware, or if the application has special demands on the RNG. + * + * Look in random.h and how the default generator is "plugged in" in + * scm_init_random(). + */ + +scm_rng scm_the_rng; + + +/* + * The prepackaged RNG + * + * This is the MWC (Multiply With Carry) random number generator + * described by George Marsaglia at the Department of Statistics and + * Supercomputer Computations Research Institute, The Florida State + * University (http://stat.fsu.edu/~geo). + * + * It uses 64 bits, has a period of 4578426017172946943 (4.6e18), and + * passes all tests in the DIEHARD test suite + * (http://stat.fsu.edu/~geo/diehard.html) + */ + +#define A 2131995753UL + +#if SIZEOF_LONG > 4 +#if SIZEOF_INT > 4 +#define LONG32 unsigned short +#else +#define LONG32 unsigned int +#endif +#define LONG64 unsigned long +#else +#define LONG32 unsigned long +#define LONG64 unsigned long long +#endif + +#if SIZEOF_LONG > 4 || defined (HAVE_LONG_LONGS) + +unsigned long +scm_i_uniform32 (scm_i_rstate *state) +{ + LONG64 x = (LONG64) A * state->w + state->c; + LONG32 w = x & 0xffffffffUL; + state->w = w; + state->c = x >> 32L; + return w; +} + +#else + +/* ww This is a portable version of the same RNG without 64 bit + * * aa arithmetic. + * ---- + * xx It is only intended to provide identical behaviour on + * xx platforms without 8 byte longs or long longs until + * xx someone has implemented the routine in assembler code. + * xxcc + * ---- + * ccww + */ + +#define L(x) ((x) & 0xffff) +#define H(x) ((x) >> 16) + +unsigned long +scm_i_uniform32 (scm_i_rstate *state) +{ + LONG32 x1 = L (A) * L (state->w); + LONG32 x2 = L (A) * H (state->w); + LONG32 x3 = H (A) * L (state->w); + LONG32 w = L (x1) + L (state->c); + LONG32 m = H (x1) + L (x2) + L (x3) + H (state->c) + H (w); + LONG32 x4 = H (A) * H (state->w); + state->w = w = (L (m) << 16) + L (w); + state->c = H (x2) + H (x3) + x4 + H (m); + return w; +} + +#endif + +void +scm_i_init_rstate (scm_i_rstate *state, char *seed, int n) +{ + LONG32 w = 0L; + LONG32 c = 0L; + int i, m; + for (i = 0; i < n; ++i) + { + m = i % 8; + if (m < 4) + w += seed[i] << (8 * m); + else + c += seed[i] << (8 * (m - 4)); + } + if ((w == 0 && c == 0) || (w == 0xffffffffUL && c == A - 1)) + ++c; + state->w = w; + state->c = c; +} + +scm_i_rstate * +scm_i_copy_rstate (scm_i_rstate *state) +{ + scm_rstate *new_state = malloc (scm_the_rng.rstate_size); + if (new_state == 0) + scm_wta (SCM_MAKINUM (scm_the_rng.rstate_size), + (char *) SCM_NALLOC, "rstate"); + return memcpy (new_state, state, scm_the_rng.rstate_size); +} + + +/* + * Random number library functions + */ + +inline double +scm_i_uniform01 (scm_rstate *state) +{ + double x = (double) scm_the_rng.random_bits (state) / (double) 0xFFFFFFFFUL; + return ((x + (double) scm_the_rng.random_bits (state)) + / (double) 0xFFFFFFFFUL); +} + +double +scm_i_normal01 (scm_rstate *state) +{ + if (state->reserved0) + { + state->reserved0 = 0; + return state->reserved1; + } + else + { + double r, a, n; + state->reserved0 = 1; + + r = sqrt (-2.0 * log (scm_i_uniform01 (state))); + a = 2.0 * M_PI * scm_i_uniform01 (state); + + n = r * sin (a); + state->reserved1 = r * cos (a); + + return n; + } +} + +double +scm_i_exp1 (scm_rstate *state) +{ + return - log (scm_i_uniform01 (state)); +} + +unsigned char scm_masktab[256]; + +unsigned long +scm_i_random (unsigned long m, scm_rstate *state) +{ + unsigned int r, mask; + mask = (m < 0x100 + ? scm_masktab[m] + : (m < 0x10000 + ? scm_masktab[m >> 8] << 8 | 0xFF + : (m < 0x1000000 + ? scm_masktab[m >> 16] << 16 | 0xFFFF + : scm_masktab[m >> 24] << 24 | 0xFFFFFF))); + while ((r = scm_the_rng.random_bits (state) & mask) >= m); + return r; +} + +SCM +scm_i_random_bignum (SCM m, scm_rstate *state) +{ + SCM b; + int i; + LONG32 *bits, w, r; + i = SCM_NUMDIGS (m); + b = scm_mkbig (i, 0); + bits = (LONG32 *) SCM_BDIGITS (b); + /* treat most significant digit specially */ +#if SIZEOF_INT == 4 + /* 16 bit digits */ + if (i & 1) + { + /* fix most significant 16 bits */ + unsigned short mask, r, s = SCM_BDIGITS (m)[--i]; + mask = s < 0x100 ? scm_masktab[s] : scm_masktab[s >> 8] << 8 | 0xFF; + while ((r = scm_the_rng.random_bits (state) & mask) > s); + + ((SCM_BIGDIG *) bits)[i] = r; + i /= 2; + if (r < s) + goto rest; + } + else +#endif + { + /* fix most significant 32 bits */ + LONG32 mask; + i /= 2; + w = ((LONG32 *) SCM_BDIGITS (m))[--i]; + mask = (w < 0x10000 + ? (w < 0x100 + ? scm_masktab[w] + : scm_masktab[w >> 8] << 8 | 0xFF) + : (w < 0x1000000 + ? scm_masktab[w >> 16] << 16 | 0xFFFF + : scm_masktab[w >> 24] << 24 | 0xFFFFFF)); + while ((r = scm_the_rng.random_bits (state) & mask) > w); + bits[i] = r; + if (r < w) + goto rest; + } + /* fill and compare */ + while (i) + { + w = ((LONG32 *) SCM_BDIGITS (m))[--i]; + while ((r = scm_the_rng.random_bits (state)) > w); + bits[i] = r; + if (r < w) + goto rest; + } + /* now fill up the rest of the bignum */ + rest: + while (i) + bits[--i] = scm_the_rng.random_bits (state); + return scm_normbig (b); +} + +/* + * Scheme level representation of random states. + */ + +long scm_tc16_rstate; + +static SCM +make_rstate (scm_rstate *state) +{ + SCM cell; + SCM_NEWCELL (cell); + SCM_ENTER_A_SECTION; + SCM_SETCDR (cell, state); + SCM_SETCAR (cell, scm_tc16_rstate); + SCM_EXIT_A_SECTION; + return cell; +} + +static int +print_rstate (SCM rstate, SCM port, scm_print_state *pstate) +{ + scm_puts ("#', port); + return 1; +} + +static scm_sizet +free_rstate (SCM rstate) +{ + free (SCM_RSTATE (rstate)); + return scm_the_rng.rstate_size; +} + +static scm_smobfuns rstate_smob = { 0, free_rstate, print_rstate, 0}; + +/* + * Scheme level interface. + */ + +SCM_GLOBAL_VCELL_INIT (scm_var_random_state, "*random-state*", scm_make_random_state (scm_makfrom0str ("URL:http://stat.fsu.edu/~geo/diehard.html"))); + +SCM_PROC (s_random, "random", 1, 1, 0, scm_random); + +SCM +scm_random (SCM n, SCM state) +{ + if (SCM_UNBNDP (state)) + state = SCM_CDR (scm_var_random_state); + SCM_ASSERT (SCM_NIMP (state) && SCM_RSTATEP (state), + state, SCM_ARG2, s_random); + if (SCM_INUMP (n)) + { + unsigned long m = SCM_INUM (n); + SCM_ASSERT (m > 0, n, SCM_ARG1, s_random); + return SCM_MAKINUM (scm_i_random (m, SCM_RSTATE (state))); + } + SCM_ASSERT (SCM_NIMP (n), n, SCM_ARG1, s_random); + if (SCM_REALP (n)) + return scm_makdbl (SCM_REALPART (n) * scm_i_uniform01 (SCM_RSTATE (state)), + 0.0); + SCM_ASSERT (SCM_TYP16 (n) == scm_tc16_bigpos, n, SCM_ARG1, s_random); + return scm_i_random_bignum (n, SCM_RSTATE (state)); +} + +SCM_PROC (s_make_random_state, "make-random-state", 0, 1, 0, scm_make_random_state); + +SCM +scm_make_random_state (SCM state) +{ + if (SCM_UNBNDP (state)) + { + state = SCM_CDR (scm_var_random_state); + goto copy_state; + } + else if (SCM_NUMBERP (state)) + { + state = scm_number_to_string (state, SCM_UNDEFINED); + goto seed; + } + else if (SCM_NIMP (state) && SCM_STRINGP (state)) + seed: + { + scm_rstate *nstate = malloc (scm_the_rng.rstate_size); + if (nstate == 0) + scm_wta (SCM_MAKINUM (scm_the_rng.rstate_size), + (char *) SCM_NALLOC, + "rstate"); + nstate->reserved0 = 0; + scm_the_rng.init_rstate (nstate, SCM_ROCHARS (state), SCM_LENGTH (state)); + return make_rstate (nstate); + } + copy_state: + SCM_ASSERT (SCM_NIMP (state) && SCM_RSTATEP (state), + state, + SCM_ARG1, + s_make_random_state); + return make_rstate (scm_the_rng.copy_rstate (SCM_RSTATE (state))); +} + +SCM_PROC (s_random_uniform, "random:uniform", 0, 1, 0, scm_random_uniform); + +SCM +scm_random_uniform (SCM state) +{ + if (SCM_UNBNDP (state)) + state = SCM_CDR (scm_var_random_state); + SCM_ASSERT (SCM_NIMP (state) && SCM_RSTATEP (state), + state, + SCM_ARG1, + s_random_uniform); + return scm_makdbl (scm_i_uniform01 (SCM_RSTATE (state)), 0.0); +} + +static void +vector_scale (SCM v, double c) +{ + int n = SCM_LENGTH (v); + if (SCM_VECTORP (v)) + while (--n >= 0) + SCM_REAL (SCM_VELTS (v)[n]) *= c; + else + while (--n >= 0) + ((double *) SCM_VELTS (v))[n] *= c; +} + +static double +vector_sum_squares (SCM v) +{ + double x, sum = 0.0; + int n = SCM_LENGTH (v); + if (SCM_VECTORP (v)) + while (--n >= 0) + { + x = SCM_REAL (SCM_VELTS (v)[n]); + sum += x * x; + } + else + while (--n >= 0) + { + x = ((double *) SCM_VELTS (v))[n]; + sum += x * x; + } + return sum; +} + +SCM_PROC (s_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0, scm_random_solid_sphere_x); + +/* For the uniform distribution on the solid sphere, note that in + * this distribution the length r of the vector has cumulative + * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be + * generated as r=u^(1/n). + */ +SCM +scm_random_solid_sphere_x (SCM v, SCM state) +{ + SCM_ASSERT (SCM_NIMP (v) + && (SCM_VECTORP (v) || SCM_TYP7 (v) == scm_tc7_dvect), + v, SCM_ARG1, s_random_solid_sphere_x); + if (SCM_UNBNDP (state)) + state = SCM_CDR (scm_var_random_state); + SCM_ASSERT (SCM_NIMP (state) && SCM_RSTATEP (state), + state, + SCM_ARG2, + s_random_solid_sphere_x); + scm_random_normal_vector_x (v, state); + vector_scale (v, + pow (scm_i_uniform01 (SCM_RSTATE (state)), + 1.0 / SCM_LENGTH (v)) + / sqrt (vector_sum_squares (v))); + return SCM_UNSPECIFIED; +} + +SCM_PROC (s_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0, scm_random_hollow_sphere_x); + +SCM +scm_random_hollow_sphere_x (SCM v, SCM state) +{ + SCM_ASSERT (SCM_NIMP (v) + && (SCM_VECTORP (v) || SCM_TYP7 (v) == scm_tc7_dvect), + v, SCM_ARG1, s_random_solid_sphere_x); + if (SCM_UNBNDP (state)) + state = SCM_CDR (scm_var_random_state); + SCM_ASSERT (SCM_NIMP (state) && SCM_RSTATEP (state), + state, + SCM_ARG2, + s_random_hollow_sphere_x); + scm_random_normal_vector_x (v, state); + vector_scale (v, 1 / sqrt (vector_sum_squares (v))); + return SCM_UNSPECIFIED; +} + +SCM_PROC (s_random_normal, "random:normal", 0, 1, 0, scm_random_normal); + +SCM +scm_random_normal (SCM state) +{ + if (SCM_UNBNDP (state)) + state = SCM_CDR (scm_var_random_state); + SCM_ASSERT (SCM_NIMP (state) && SCM_RSTATEP (state), + state, + SCM_ARG1, + s_random_normal); + return scm_makdbl (scm_i_normal01 (SCM_RSTATE (state)), 0.0); +} + +SCM_PROC (s_random_normal_vector_x, "random:normal-vector!", 1, 1, 0, scm_random_normal_vector_x); + +SCM +scm_random_normal_vector_x (SCM v, SCM state) +{ + int n; + SCM_ASSERT (SCM_NIMP (v) + && (SCM_VECTORP (v) || SCM_TYP7 (v) == scm_tc7_dvect), + v, SCM_ARG1, s_random_solid_sphere_x); + if (SCM_UNBNDP (state)) + state = SCM_CDR (scm_var_random_state); + SCM_ASSERT (SCM_NIMP (state) && SCM_RSTATEP (state), + state, + SCM_ARG2, + s_random_normal_vector_x); + n = SCM_LENGTH (v); + if (SCM_VECTORP (v)) + while (--n >= 0) + SCM_VELTS (v)[n] = scm_makdbl (scm_i_normal01 (SCM_RSTATE (state)), 0.0); + else + while (--n >= 0) + ((double *) SCM_VELTS (v))[n] = scm_i_normal01 (SCM_RSTATE (state)); + return SCM_UNSPECIFIED; +} + +SCM_PROC (s_random_exp, "random:exp", 0, 1, 0, scm_random_exp); + +SCM +scm_random_exp (SCM state) +{ + if (SCM_UNBNDP (state)) + state = SCM_CDR (scm_var_random_state); + SCM_ASSERT (SCM_NIMP (state) && SCM_RSTATEP (state), + state, + SCM_ARG1, + s_random_exp); + return scm_makdbl (scm_i_exp1 (SCM_RSTATE (state)), 0.0); +} + +void +scm_init_random () +{ + int i, m; + /* plug in default RNG */ + scm_rng rng = + { + sizeof (scm_i_rstate), + (unsigned long (*)()) scm_i_uniform32, + (void (*)()) scm_i_init_rstate, + (scm_rstate *(*)()) scm_i_copy_rstate + }; + scm_the_rng = rng; + + scm_tc16_rstate = scm_newsmob (&rstate_smob); + + for (m = 1; m <= 0x100; m <<= 1) + for (i = m >> 1; i < m; ++i) + scm_masktab[i] = m - 1; + +#include "random.x" + + scm_add_feature ("random"); +} diff --git a/libguile/random.h b/libguile/random.h new file mode 100644 index 000000000..9aa769270 --- /dev/null +++ b/libguile/random.h @@ -0,0 +1,122 @@ +/* classes: h_files */ + +#ifndef RANDOMH +#define RANDOMH +/* Copyright (C) 1999 Free Software Foundation, Inc. + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2, or (at your option) + * any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this software; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 59 Temple Place, Suite 330, + * Boston, MA 02111-1307 USA + * + * As a special exception, the Free Software Foundation gives permission + * for additional uses of the text contained in its release of GUILE. + * + * The exception is that, if you link the GUILE library with other files + * to produce an executable, this does not by itself cause the + * resulting executable to be covered by the GNU General Public License. + * Your use of that executable is in no way restricted on account of + * linking the GUILE library code into it. + * + * This exception does not however invalidate any other reasons why + * the executable file might be covered by the GNU General Public License. + * + * This exception applies only to the code released by the + * Free Software Foundation under the name GUILE. If you copy + * code from other Free Software Foundation releases into a copy of + * GUILE, as the General Public License permits, the exception does + * not apply to the code that you add in this way. To avoid misleading + * anyone as to the status of such modified files, you must delete + * this exception notice from them. + * + * If you write modifications of your own for GUILE, it is your choice + * whether to permit this exception to apply to your modifications. + * If you do not wish that, delete this exception notice. */ + + +#include "libguile/__scm.h" + + +/* + * A plugin interface for RNGs + * + * Using this interface, it is possible for the application to tell + * libguile to use a different RNG. This is desirable if it is + * necessary to use the same RNG everywhere in the application in + * order to prevent interference, if the application uses RNG + * hardware, or if the application has special demands on the RNG. + * + * Look how the default generator is "plugged in" in scm_init_random(). + */ + +typedef struct scm_rstate { + int reserved0; + double reserved1; + /* Custom fields follow here */ +} scm_rstate; + +typedef struct scm_rng { + size_t rstate_size; /* size of random state */ + unsigned long (*random_bits) (scm_rstate *state); /* gives 32 random bits */ + void (*init_rstate) (scm_rstate *state, char *seed, int n); + scm_rstate *(*copy_rstate) (scm_rstate *state); +} scm_rng; + +extern scm_rng scm_the_rng; + + +/* + * Default RNG + */ +typedef struct scm_i_rstate { + scm_rstate rstate; + unsigned long w; + unsigned long c; +} scm_i_rstate; + +extern unsigned long scm_i_uniform32 (scm_i_rstate *); +extern void scm_i_init_rstate (scm_i_rstate *, char *seed, int n); +extern scm_i_rstate *scm_i_copy_rstate (scm_i_rstate *); + + +/* + * Random number library functions + */ +extern double scm_i_uniform01 (scm_rstate *); +extern double scm_i_normal01 (scm_rstate *); +extern double scm_i_exp1 (scm_rstate *); +extern unsigned long scm_i_random (unsigned long m, scm_rstate *); +extern SCM scm_i_random_bignum (SCM m, scm_rstate *); + + +/* + * Scheme level interface + */ +extern long scm_tc16_rstate; +#define SCM_RSTATE(obj) ((scm_rstate *) SCM_CDR (obj)) +#define SCM_RSTATEP(obj) (SCM_TYP16 (obj) == scm_tc16_rstate) + +extern unsigned char scm_masktab[256]; + +extern SCM scm_var_random_state; +extern SCM scm_random (SCM n, SCM state); +extern SCM scm_make_random_state (SCM state); +extern SCM scm_random_uniform (SCM state); +extern SCM scm_random_solid_sphere_x (SCM v, SCM state); +extern SCM scm_random_hollow_sphere_x (SCM v, SCM state); +extern SCM scm_random_normal (SCM state); +extern SCM scm_random_normal_vector_x (SCM v, SCM state); +extern SCM scm_random_exp (SCM state); +extern void scm_init_random (void); + +#endif /* RANDOMH */