libguile/random.c

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