guile18/libguile/random.c

   1 /* Copyright (C) 1999,2000,2001, 2003, 2005, 2006, 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
   4  * License as published by the Free Software Foundation; either
   5  * version 2.1 of the License, or (at your option) any later version.
   6  *
   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.
  11  *
  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 02110-1301 USA
  15  */
  16
  17
  18
  19 /* Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */
  20
  21 #ifdef HAVE_CONFIG_H
  22 #  include <config.h>
  23 #endif
  24
  25 #include "libguile/_scm.h"
  26
  27 #include <gmp.h>
  28 #include <stdio.h>
  29 #include <math.h>
  30 #include <string.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/srfi-4.h"
  37 #include "libguile/vectors.h"
  38
  39 #include "libguile/validate.h"
  40 #include "libguile/random.h"
  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
  56 scm_t_rng scm_the_rng;
  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
  74 #ifndef M_PI
  75 #define M_PI 3.14159265359
  76 #endif
  77
  78 unsigned long
  79 scm_i_uniform32 (scm_t_i_rstate *state)
  80 {
  81   scm_t_uint64 x = (scm_t_uint64) A * state->w + state->c;
  82   scm_t_uint32 w = x & 0xffffffffUL;
  83   state->w = w;
  84   state->c = x >> 32L;
  85   return w;
  86 }
  87
  88 void
  89 scm_i_init_rstate (scm_t_i_rstate *state, const char *seed, int n)
  90 {
  91   scm_t_uint32 w = 0L;
  92   scm_t_uint32 c = 0L;
  93   int i, m;
  94   for (i = 0; i < n; ++i)
  95     {
  96       m = i % 8;
  97       if (m < 4)
  98         w += seed[i] << (8 * m);
  99       else
 100         c += seed[i] << (8 * (m - 4));
 101     }
 102   if ((w == 0 && c == 0) || (w == -1 && c == A - 1))
 103     ++c;
 104   state->w = w;
 105   state->c = c;
 106 }
 107
 108 scm_t_i_rstate *
 109 scm_i_copy_rstate (scm_t_i_rstate *state)
 110 {
 111   scm_t_rstate *new_state = scm_malloc (scm_the_rng.rstate_size);
 112   return memcpy (new_state, state, scm_the_rng.rstate_size);
 113 }
 114
 115 \f
 116 /*
 117  * Random number library functions
 118  */
 119
 120 scm_t_rstate *
 121 scm_c_make_rstate (const char *seed, int n)
 122 {
 123   scm_t_rstate *state = scm_malloc (scm_the_rng.rstate_size);
 124   state->reserved0 = 0;
 125   scm_the_rng.init_rstate (state, seed, n);
 126   return state;
 127 }
 128
 129
 130 scm_t_rstate *
 131 scm_c_default_rstate ()
 132 #define FUNC_NAME "scm_c_default_rstate"
 133 {
 134   SCM state = SCM_VARIABLE_REF (scm_var_random_state);
 135   if (!SCM_RSTATEP (state))
 136     SCM_MISC_ERROR ("*random-state* contains bogus random state", SCM_EOL);
 137   return SCM_RSTATE (state);
 138 }
 139 #undef FUNC_NAME
 140
 141
 142 inline double
 143 scm_c_uniform01 (scm_t_rstate *state)
 144 {
 145   double x = (double) scm_the_rng.random_bits (state) / (double) 0xffffffffUL;
 146   return ((x + (double) scm_the_rng.random_bits (state))
 147           / (double) 0xffffffffUL);
 148 }
 149
 150 double
 151 scm_c_normal01 (scm_t_rstate *state)
 152 {
 153   if (state->reserved0)
 154     {
 155       state->reserved0 = 0;
 156       return state->reserved1;
 157     }
 158   else
 159     {
 160       double r, a, n;
 161
 162       r = sqrt (-2.0 * log (scm_c_uniform01 (state)));
 163       a = 2.0 * M_PI * scm_c_uniform01 (state);
 164
 165       n = r * sin (a);
 166       state->reserved1 = r * cos (a);
 167       state->reserved0 = 1;
 168
 169       return n;
 170     }
 171 }
 172
 173 double
 174 scm_c_exp1 (scm_t_rstate *state)
 175 {
 176   return - log (scm_c_uniform01 (state));
 177 }
 178
 179 unsigned char scm_masktab[256];
 180
 181 static inline scm_t_uint32
 182 scm_i_mask32 (scm_t_uint32 m)
 183 {
 184   return (m < 0x100
 185           ? scm_masktab[m]
 186           : (m < 0x10000
 187              ? scm_masktab[m >> 8] << 8 | 0xff
 188              : (m < 0x1000000
 189                 ? scm_masktab[m >> 16] << 16 | 0xffff
 190                 : scm_masktab[m >> 24] << 24 | 0xffffff)));
 191 }
 192
 193 static scm_t_uint32
 194 scm_c_random32 (scm_t_rstate *state, scm_t_uint32 m)
 195 {
 196   scm_t_uint32 r, mask = scm_i_mask32 (m);
 197   while ((r = scm_the_rng.random_bits (state) & mask) >= m);
 198   return r;
 199 }
 200
 201 /* Returns 32 random bits. */
 202 unsigned long
 203 scm_c_random (scm_t_rstate *state, unsigned long m)
 204 {
 205   return scm_c_random32 (state, (scm_t_uint32)m);
 206 }
 207
 208 scm_t_uint64
 209 scm_c_random64 (scm_t_rstate *state, scm_t_uint64 m)
 210 {
 211   scm_t_uint64 r;
 212   scm_t_uint32 mask;
 213
 214   if (m <= SCM_T_UINT32_MAX)
 215     return scm_c_random32 (state, (scm_t_uint32) m);
 216
 217   mask = scm_i_mask32 (m >> 32);
 218   while ((r = ((scm_t_uint64) (scm_the_rng.random_bits (state) & mask) << 32)
 219           | scm_the_rng.random_bits (state)) >= m)
 220     ;
 221   return r;
 222 }
 223
 224 /*
 225   SCM scm_c_random_bignum (scm_t_rstate *state, SCM m)
 226
 227   Takes a random state (source of random bits) and a bignum m.
 228   Returns a bignum b, 0 <= b < m.
 229
 230   It does this by allocating a bignum b with as many base 65536 digits
 231   as m, filling b with random bits (in 32 bit chunks) up to the most
 232   significant 1 in m, and, finally checking if the resultant b is too
 233   large (>= m).  If too large, we simply repeat the process again.  (It
 234   is important to throw away all generated random bits if b >= m,
 235   otherwise we'll end up with a distorted distribution.)
 236
 237 */
 238
 239 SCM
 240 scm_c_random_bignum (scm_t_rstate *state, SCM m)
 241 {
 242   SCM result = scm_i_mkbig ();
 243   const size_t m_bits = mpz_sizeinbase (SCM_I_BIG_MPZ (m), 2);
 244   /* how many bits would only partially fill the last u32? */
 245   const size_t end_bits = m_bits % (sizeof (scm_t_uint32) * SCM_CHAR_BIT);
 246   scm_t_uint32 *random_chunks = NULL;
 247   const scm_t_uint32 num_full_chunks =
 248     m_bits / (sizeof (scm_t_uint32) * SCM_CHAR_BIT);
 249   const scm_t_uint32 num_chunks = num_full_chunks + ((end_bits) ? 1 : 0);
 250
 251   /* we know the result will be this big */
 252   mpz_realloc2 (SCM_I_BIG_MPZ (result), m_bits);
 253
 254   random_chunks =
 255     (scm_t_uint32 *) scm_gc_calloc (num_chunks * sizeof (scm_t_uint32),
 256                                      "random bignum chunks");
 257
 258   do
 259     {
 260       scm_t_uint32 *current_chunk = random_chunks + (num_chunks - 1);
 261       scm_t_uint32 chunks_left = num_chunks;
 262
 263       mpz_set_ui (SCM_I_BIG_MPZ (result), 0);
 264
 265       if (end_bits)
 266         {
 267           /* generate a mask with ones in the end_bits position, i.e. if
 268              end_bits is 3, then we'd have a mask of ...0000000111 */
 269           const unsigned long rndbits = scm_the_rng.random_bits (state);
 270           int rshift = (sizeof (scm_t_uint32) * SCM_CHAR_BIT) - end_bits;
 271           scm_t_uint32 mask = 0xffffffff >> rshift;
 272           scm_t_uint32 highest_bits = ((scm_t_uint32) rndbits) & mask;
 273           *current_chunk-- = highest_bits;
 274           chunks_left--;
 275         }
 276
 277       while (chunks_left)
 278         {
 279           /* now fill in the remaining scm_t_uint32 sized chunks */
 280           *current_chunk-- = scm_the_rng.random_bits (state);
 281           chunks_left--;
 282         }
 283       mpz_import (SCM_I_BIG_MPZ (result),
 284                   num_chunks,
 285                   -1,
 286                   sizeof (scm_t_uint32),
 287                   0,
 288                   0,
 289                   random_chunks);
 290       /* if result >= m, regenerate it (it is important to regenerate
 291          all bits in order not to get a distorted distribution) */
 292     } while (mpz_cmp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (m)) >= 0);
 293   scm_gc_free (random_chunks,
 294                num_chunks * sizeof (scm_t_uint32),
 295                "random bignum chunks");
 296   return scm_i_normbig (result);
 297 }
 298
 299 /*
 300  * Scheme level representation of random states.
 301  */
 302
 303 scm_t_bits scm_tc16_rstate;
 304
 305 static SCM
 306 make_rstate (scm_t_rstate *state)
 307 {
 308   SCM_RETURN_NEWSMOB (scm_tc16_rstate, state);
 309 }
 310
 311 static size_t
 312 rstate_free (SCM rstate)
 313 {
 314   free (SCM_RSTATE (rstate));
 315   return 0;
 316 }
 317
 318 /*
 319  * Scheme level interface.
 320  */
 321
 322 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")));
 323
 324 SCM_DEFINE (scm_random, "random", 1, 1, 0,
 325             (SCM n, SCM state),
 326             "Return a number in [0, N).\n"
 327             "\n"
 328             "Accepts a positive integer or real n and returns a\n"
 329             "number of the same type between zero (inclusive) and\n"
 330             "N (exclusive). The values returned have a uniform\n"
 331             "distribution.\n"
 332             "\n"
 333             "The optional argument @var{state} must be of the type produced\n"
 334             "by @code{seed->random-state}. It defaults to the value of the\n"
 335             "variable @var{*random-state*}. This object is used to maintain\n"
 336             "the state of the pseudo-random-number generator and is altered\n"
 337             "as a side effect of the random operation.")
 338 #define FUNC_NAME s_scm_random
 339 {
 340   if (SCM_UNBNDP (state))
 341     state = SCM_VARIABLE_REF (scm_var_random_state);
 342   SCM_VALIDATE_RSTATE (2, state);
 343   if (SCM_I_INUMP (n))
 344     {
 345       unsigned long m = (unsigned long) SCM_I_INUM (n);
 346       SCM_ASSERT_RANGE (1, n, SCM_I_INUM (n) > 0);
 347 #if SCM_SIZEOF_UNSIGNED_LONG <= 4
 348       return scm_from_uint32 (scm_c_random (SCM_RSTATE (state),
 349                                             (scm_t_uint32) m));
 350 #elif SCM_SIZEOF_UNSIGNED_LONG <= 8
 351       return scm_from_uint64 (scm_c_random64 (SCM_RSTATE (state),
 352                                               (scm_t_uint64) m));
 353 #else
 354 #error "Cannot deal with this platform's unsigned long size"
 355 #endif
 356     }
 357   SCM_VALIDATE_NIM (1, n);
 358   if (SCM_REALP (n))
 359     return scm_from_double (SCM_REAL_VALUE (n)
 360                             * scm_c_uniform01 (SCM_RSTATE (state)));
 361
 362   if (!SCM_BIGP (n))
 363     SCM_WRONG_TYPE_ARG (1, n);
 364   return scm_c_random_bignum (SCM_RSTATE (state), n);
 365 }
 366 #undef FUNC_NAME
 367
 368 SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0,
 369             (SCM state),
 370             "Return a copy of the random state @var{state}.")
 371 #define FUNC_NAME s_scm_copy_random_state
 372 {
 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)));
 377 }
 378 #undef FUNC_NAME
 379
 380 SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0,
 381             (SCM seed),
 382             "Return a new random state using @var{seed}.")
 383 #define FUNC_NAME s_scm_seed_to_random_state
 384 {
 385   SCM res;
 386   if (SCM_NUMBERP (seed))
 387     seed = scm_number_to_string (seed, SCM_UNDEFINED);
 388   SCM_VALIDATE_STRING (1, seed);
 389   res = make_rstate (scm_c_make_rstate (scm_i_string_chars (seed),
 390                                         scm_i_string_length (seed)));
 391   scm_remember_upto_here_1 (seed);
 392   return res;
 393
 394 }
 395 #undef FUNC_NAME
 396
 397 SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0,
 398             (SCM state),
 399             "Return a uniformly distributed inexact real random number in\n"
 400             "[0,1).")
 401 #define FUNC_NAME s_scm_random_uniform
 402 {
 403   if (SCM_UNBNDP (state))
 404     state = SCM_VARIABLE_REF (scm_var_random_state);
 405   SCM_VALIDATE_RSTATE (1, state);
 406   return scm_from_double (scm_c_uniform01 (SCM_RSTATE (state)));
 407 }
 408 #undef FUNC_NAME
 409
 410 SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0,
 411             (SCM state),
 412             "Return an inexact real in a normal distribution.  The\n"
 413             "distribution used has mean 0 and standard deviation 1.  For a\n"
 414             "normal distribution with mean m and standard deviation d use\n"
 415             "@code{(+ m (* d (random:normal)))}.")
 416 #define FUNC_NAME s_scm_random_normal
 417 {
 418   if (SCM_UNBNDP (state))
 419     state = SCM_VARIABLE_REF (scm_var_random_state);
 420   SCM_VALIDATE_RSTATE (1, state);
 421   return scm_from_double (scm_c_normal01 (SCM_RSTATE (state)));
 422 }
 423 #undef FUNC_NAME
 424
 425 static void
 426 vector_scale_x (SCM v, double c)
 427 {
 428   size_t n;
 429   if (scm_is_simple_vector (v))
 430     {
 431       n = SCM_SIMPLE_VECTOR_LENGTH (v);
 432       while (n-- > 0)
 433         SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n)) *= c;
 434     }
 435   else
 436     {
 437       /* must be a f64vector. */
 438       scm_t_array_handle handle;
 439       size_t i, len;
 440       ssize_t inc;
 441       double *elts;
 442
 443       elts = scm_f64vector_writable_elements (v, &handle, &len, &inc);
 444
 445       for (i = 0; i < len; i++, elts += inc)
 446         *elts *= c;
 447
 448       scm_array_handle_release (&handle);
 449     }
 450 }
 451
 452 static double
 453 vector_sum_squares (SCM v)
 454 {
 455   double x, sum = 0.0;
 456   size_t n;
 457   if (scm_is_simple_vector (v))
 458     {
 459       n = SCM_SIMPLE_VECTOR_LENGTH (v);
 460       while (n-- > 0)
 461         {
 462           x = SCM_REAL_VALUE (SCM_SIMPLE_VECTOR_REF (v, n));
 463           sum += x * x;
 464         }
 465     }
 466   else
 467     {
 468       /* must be a f64vector. */
 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)
 477         {
 478           x = *elts;
 479           sum += x * x;
 480         }
 481
 482       scm_array_handle_release (&handle);
 483     }
 484   return sum;
 485 }
 486
 487 /* For the uniform distribution on the solid sphere, note that in
 488  * this distribution the length r of the vector has cumulative
 489  * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be
 490  * generated as r=u^(1/n).
 491  */
 492 SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0,
 493             (SCM v, SCM state),
 494             "Fills @var{vect} with inexact real random numbers the sum of\n"
 495             "whose squares is less than 1.0.  Thinking of @var{vect} as\n"
 496             "coordinates in space of dimension @var{n} @math{=}\n"
 497             "@code{(vector-length @var{vect})}, the coordinates are\n"
 498             "uniformly distributed within the unit @var{n}-sphere.")
 499 #define FUNC_NAME s_scm_random_solid_sphere_x
 500 {
 501   if (SCM_UNBNDP (state))
 502     state = SCM_VARIABLE_REF (scm_var_random_state);
 503   SCM_VALIDATE_RSTATE (2, state);
 504   scm_random_normal_vector_x (v, state);
 505   vector_scale_x (v,
 506                   pow (scm_c_uniform01 (SCM_RSTATE (state)),
 507                        1.0 / scm_c_generalized_vector_length (v))
 508                   / sqrt (vector_sum_squares (v)));
 509   return SCM_UNSPECIFIED;
 510 }
 511 #undef FUNC_NAME
 512
 513 SCM_DEFINE (scm_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0,
 514             (SCM v, SCM state),
 515             "Fills vect with inexact real random numbers\n"
 516             "the sum of whose squares is equal to 1.0.\n"
 517             "Thinking of vect as coordinates in space of\n"
 518             "dimension n = (vector-length vect), the coordinates\n"
 519             "are uniformly distributed over the surface of the\n"
 520             "unit n-sphere.")
 521 #define FUNC_NAME s_scm_random_hollow_sphere_x
 522 {
 523   if (SCM_UNBNDP (state))
 524     state = SCM_VARIABLE_REF (scm_var_random_state);
 525   SCM_VALIDATE_RSTATE (2, state);
 526   scm_random_normal_vector_x (v, state);
 527   vector_scale_x (v, 1 / sqrt (vector_sum_squares (v)));
 528   return SCM_UNSPECIFIED;
 529 }
 530 #undef FUNC_NAME
 531
 532
 533 SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0,
 534             (SCM v, SCM state),
 535             "Fills vect with inexact real random numbers that are\n"
 536             "independent and standard normally distributed\n"
 537             "(i.e., with mean 0 and variance 1).")
 538 #define FUNC_NAME s_scm_random_normal_vector_x
 539 {
 540   long i;
 541   scm_t_array_handle handle;
 542   scm_t_array_dim *dim;
 543
 544   if (SCM_UNBNDP (state))
 545     state = SCM_VARIABLE_REF (scm_var_random_state);
 546   SCM_VALIDATE_RSTATE (2, state);
 547
 548   scm_generalized_vector_get_handle (v, &handle);
 549   dim = scm_array_handle_dims (&handle);
 550
 551   if (scm_is_vector (v))
 552     {
 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)));
 556     }
 557   else
 558     {
 559       /* must be a f64vector. */
 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));
 563     }
 564
 565   scm_array_handle_release (&handle);
 566
 567   return SCM_UNSPECIFIED;
 568 }
 569 #undef FUNC_NAME
 570
 571 SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0,
 572             (SCM state),
 573             "Return an inexact real in an exponential distribution with mean\n"
 574             "1.  For an exponential distribution with mean u use (* u\n"
 575             "(random:exp)).")
 576 #define FUNC_NAME s_scm_random_exp
 577 {
 578   if (SCM_UNBNDP (state))
 579     state = SCM_VARIABLE_REF (scm_var_random_state);
 580   SCM_VALIDATE_RSTATE (1, state);
 581   return scm_from_double (scm_c_exp1 (SCM_RSTATE (state)));
 582 }
 583 #undef FUNC_NAME
 584
 585 void
 586 scm_init_random ()
 587 {
 588   int i, m;
 589   /* plug in default RNG */
 590   scm_t_rng rng =
 591   {
 592     sizeof (scm_t_i_rstate),
 593     (unsigned long (*)()) scm_i_uniform32,
 594     (void (*)())          scm_i_init_rstate,
 595     (scm_t_rstate *(*)())    scm_i_copy_rstate
 596   };
 597   scm_the_rng = rng;
 598
 599   scm_tc16_rstate = scm_make_smob_type ("random-state", 0);
 600   scm_set_smob_free (scm_tc16_rstate, rstate_free);
 601
 602   for (m = 1; m <= 0x100; m <<= 1)
 603     for (i = m >> 1; i < m; ++i)
 604       scm_masktab[i] = m - 1;
 605
 606 #include "libguile/random.x"
 607
 608   scm_add_feature ("random");
 609 }
 610
 611 /*
 612   Local Variables:
 613   c-file-style: "gnu"
 614   End:
 615 */