libguile/random.c

   1 /* Copyright (C) 1999,2000,2001, 2003 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  15  */
  16
  17
  18
  19 /* Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */
  20
  21 #if 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 #ifdef SCM_HAVE_T_INT64
  79
  80 unsigned long
  81 scm_i_uniform32 (scm_t_i_rstate *state)
  82 {
  83   scm_t_int64 x = (scm_t_int64) A * state->w + state->c;
  84   scm_t_int32 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_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);
 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_int32 w = 0L;
 126   scm_t_int32 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 == 0xffffffffUL && 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   if (new_state == 0)
 147     scm_memory_error ("rstate");
 148   return memcpy (new_state, state, scm_the_rng.rstate_size);
 149 }
 150
 151 \f
 152 /*
 153  * Random number library functions
 154  */
 155
 156 scm_t_rstate *
 157 scm_c_make_rstate (const char *seed, int n)
 158 {
 159   scm_t_rstate *state = scm_malloc (scm_the_rng.rstate_size);
 160   if (state == 0)
 161     scm_memory_error ("rstate");
 162   state->reserved0 = 0;
 163   scm_the_rng.init_rstate (state, seed, n);
 164   return state;
 165 }
 166
 167
 168 scm_t_rstate *
 169 scm_c_default_rstate ()
 170 #define FUNC_NAME "scm_c_default_rstate"
 171 {
 172   SCM state = SCM_VARIABLE_REF (scm_var_random_state);
 173   if (!SCM_RSTATEP (state))
 174     SCM_MISC_ERROR ("*random-state* contains bogus random state", SCM_EOL);
 175   return SCM_RSTATE (state);
 176 }
 177 #undef FUNC_NAME
 178
 179
 180 inline double
 181 scm_c_uniform01 (scm_t_rstate *state)
 182 {
 183   double x = (double) scm_the_rng.random_bits (state) / (double) 0xffffffffUL;
 184   return ((x + (double) scm_the_rng.random_bits (state))
 185           / (double) 0xffffffffUL);
 186 }
 187
 188 double
 189 scm_c_normal01 (scm_t_rstate *state)
 190 {
 191   if (state->reserved0)
 192     {
 193       state->reserved0 = 0;
 194       return state->reserved1;
 195     }
 196   else
 197     {
 198       double r, a, n;
 199
 200       r = sqrt (-2.0 * log (scm_c_uniform01 (state)));
 201       a = 2.0 * M_PI * scm_c_uniform01 (state);
 202
 203       n = r * sin (a);
 204       state->reserved1 = r * cos (a);
 205       state->reserved0 = 1;
 206
 207       return n;
 208     }
 209 }
 210
 211 double
 212 scm_c_exp1 (scm_t_rstate *state)
 213 {
 214   return - log (scm_c_uniform01 (state));
 215 }
 216
 217 unsigned char scm_masktab[256];
 218
 219 unsigned long
 220 scm_c_random (scm_t_rstate *state, unsigned long m)
 221 {
 222   unsigned int r, mask;
 223   mask = (m < 0x100
 224           ? scm_masktab[m]
 225           : (m < 0x10000
 226              ? scm_masktab[m >> 8] << 8 | 0xff
 227              : (m < 0x1000000
 228                 ? scm_masktab[m >> 16] << 16 | 0xffff
 229                 : scm_masktab[m >> 24] << 24 | 0xffffff)));
 230   while ((r = scm_the_rng.random_bits (state) & mask) >= m);
 231   return r;
 232 }
 233
 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
 249 SCM
 250 scm_c_random_bignum (scm_t_rstate *state, SCM m)
 251 {
 252   SCM result = scm_i_mkbig ();
 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? */
 255   const size_t end_bits = m_bits % (sizeof (unsigned long) * SCM_CHAR_BIT);
 256   unsigned long *random_chunks = NULL;
 257   const unsigned long num_full_chunks =
 258     m_bits / (sizeof (unsigned long) * SCM_CHAR_BIT);
 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
 268   do
 269     {
 270       unsigned long *current_chunk = random_chunks + (num_chunks - 1);
 271       unsigned long chunks_left = num_chunks;
 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);
 280           int rshift = (sizeof (unsigned long) * SCM_CHAR_BIT) - end_bits;
 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);
 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);
 303   scm_gc_free (random_chunks,
 304                num_chunks * sizeof (unsigned long),
 305                "random bignum chunks");
 306   return scm_i_normbig (result);
 307 }
 308
 309 /*
 310  * Scheme level representation of random states.
 311  */
 312
 313 scm_t_bits scm_tc16_rstate;
 314
 315 static SCM
 316 make_rstate (scm_t_rstate *state)
 317 {
 318   SCM_RETURN_NEWSMOB (scm_tc16_rstate, state);
 319 }
 320
 321 static size_t
 322 rstate_free (SCM rstate)
 323 {
 324   free (SCM_RSTATE (rstate));
 325   return 0;
 326 }
 327
 328 /*
 329  * Scheme level interface.
 330  */
 331
 332 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")));
 333
 334 SCM_DEFINE (scm_random, "random", 1, 1, 0,
 335             (SCM n, SCM state),
 336             "Return a number in [0, N).\n"
 337             "\n"
 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"
 341             "distribution.\n"
 342             "\n"
 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.")
 348 #define FUNC_NAME s_scm_random
 349 {
 350   if (SCM_UNBNDP (state))
 351     state = SCM_VARIABLE_REF (scm_var_random_state);
 352   SCM_VALIDATE_RSTATE (2, state);
 353   if (SCM_I_INUMP (n))
 354     {
 355       unsigned long m = SCM_I_INUM (n);
 356       SCM_ASSERT_RANGE (1, n, m > 0);
 357       return scm_from_ulong (scm_c_random (SCM_RSTATE (state), m));
 358     }
 359   SCM_VALIDATE_NIM (1, n);
 360   if (SCM_REALP (n))
 361     return scm_from_double (SCM_REAL_VALUE (n)
 362                             * scm_c_uniform01 (SCM_RSTATE (state)));
 363
 364   if (!SCM_BIGP (n))
 365     SCM_WRONG_TYPE_ARG (1, n);
 366   return scm_c_random_bignum (SCM_RSTATE (state), n);
 367 }
 368 #undef FUNC_NAME
 369
 370 SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0,
 371             (SCM state),
 372             "Return a copy of the random state @var{state}.")
 373 #define FUNC_NAME s_scm_copy_random_state
 374 {
 375   if (SCM_UNBNDP (state))
 376     state = SCM_VARIABLE_REF (scm_var_random_state);
 377   SCM_VALIDATE_RSTATE (1, state);
 378   return make_rstate (scm_the_rng.copy_rstate (SCM_RSTATE (state)));
 379 }
 380 #undef FUNC_NAME
 381
 382 SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0,
 383             (SCM seed),
 384             "Return a new random state using @var{seed}.")
 385 #define FUNC_NAME s_scm_seed_to_random_state
 386 {
 387   SCM res;
 388   if (SCM_NUMBERP (seed))
 389     seed = scm_number_to_string (seed, SCM_UNDEFINED);
 390   SCM_VALIDATE_STRING (1, seed);
 391   res = make_rstate (scm_c_make_rstate (scm_i_string_chars (seed),
 392                                         scm_i_string_length (seed)));
 393   scm_remember_upto_here_1 (seed);
 394   return res;
 395
 396 }
 397 #undef FUNC_NAME
 398
 399 SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0,
 400             (SCM state),
 401             "Return a uniformly distributed inexact real random number in\n"
 402             "[0,1).")
 403 #define FUNC_NAME s_scm_random_uniform
 404 {
 405   if (SCM_UNBNDP (state))
 406     state = SCM_VARIABLE_REF (scm_var_random_state);
 407   SCM_VALIDATE_RSTATE (1, state);
 408   return scm_from_double (scm_c_uniform01 (SCM_RSTATE (state)));
 409 }
 410 #undef FUNC_NAME
 411
 412 SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0,
 413             (SCM state),
 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)))}.")
 418 #define FUNC_NAME s_scm_random_normal
 419 {
 420   if (SCM_UNBNDP (state))
 421     state = SCM_VARIABLE_REF (scm_var_random_state);
 422   SCM_VALIDATE_RSTATE (1, state);
 423   return scm_from_double (scm_c_normal01 (SCM_RSTATE (state)));
 424 }
 425 #undef FUNC_NAME
 426
 427 #if SCM_HAVE_ARRAYS
 428
 429 static void
 430 vector_scale_x (SCM v, double c)
 431 {
 432   size_t n;
 433   if (SCM_VECTORP (v))
 434     {
 435       n = SCM_VECTOR_LENGTH (v);
 436       while (n-- > 0)
 437         SCM_REAL_VALUE (SCM_VELTS (v)[n]) *= c;
 438     }
 439   else
 440     {
 441       /* must be a f64vector. */
 442       double *elts = scm_f64vector_elements (v);
 443       n = scm_c_uniform_vector_length (v);
 444       while (n-- > 0)
 445         elts[n] *= c;
 446       scm_uniform_vector_release (v);
 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_VECTORP (v))
 456     {
 457       n = SCM_VECTOR_LENGTH (v);
 458       while (n-- > 0)
 459         {
 460           x = SCM_REAL_VALUE (SCM_VELTS (v)[n]);
 461           sum += x * x;
 462         }
 463     }
 464   else
 465     {
 466       /* must be a f64vector. */
 467       double *elts = scm_f64vector_elements (v);
 468       n = scm_c_uniform_vector_length (v);
 469       while (n-- > 0)
 470         {
 471           x = elts[n];
 472           sum += x * x;
 473         }
 474       scm_uniform_vector_release (v);
 475     }
 476   return sum;
 477 }
 478
 479 /* For the uniform distribution on the solid sphere, note that in
 480  * this distribution the length r of the vector has cumulative
 481  * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be
 482  * generated as r=u^(1/n).
 483  */
 484 SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0,
 485             (SCM v, SCM state),
 486             "Fills vect with inexact real random numbers\n"
 487             "the sum of whose squares is less than 1.0.\n"
 488             "Thinking of vect as coordinates in space of\n"
 489             "dimension n = (vector-length vect), the coordinates\n"
 490             "are uniformly distributed within the unit n-sphere.\n"
 491             "The sum of the squares of the numbers is returned.")
 492 #define FUNC_NAME s_scm_random_solid_sphere_x
 493 {
 494   SCM_VALIDATE_VECTOR_OR_DVECTOR (1, v);
 495   if (SCM_UNBNDP (state))
 496     state = SCM_VARIABLE_REF (scm_var_random_state);
 497   SCM_VALIDATE_RSTATE (2, state);
 498   scm_random_normal_vector_x (v, state);
 499   vector_scale_x (v,
 500                   pow (scm_c_uniform01 (SCM_RSTATE (state)),
 501                        1.0 / scm_to_int (scm_uniform_vector_length (v)))
 502                   / sqrt (vector_sum_squares (v)));
 503   return SCM_UNSPECIFIED;
 504 }
 505 #undef FUNC_NAME
 506
 507 SCM_DEFINE (scm_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0,
 508             (SCM v, SCM state),
 509             "Fills vect with inexact real random numbers\n"
 510             "the sum of whose squares is equal to 1.0.\n"
 511             "Thinking of vect as coordinates in space of\n"
 512             "dimension n = (vector-length vect), the coordinates\n"
 513             "are uniformly distributed over the surface of the\n"
 514             "unit n-sphere.")
 515 #define FUNC_NAME s_scm_random_hollow_sphere_x
 516 {
 517   SCM_VALIDATE_VECTOR_OR_DVECTOR (1, v);
 518   if (SCM_UNBNDP (state))
 519     state = SCM_VARIABLE_REF (scm_var_random_state);
 520   SCM_VALIDATE_RSTATE (2, state);
 521   scm_random_normal_vector_x (v, state);
 522   vector_scale_x (v, 1 / sqrt (vector_sum_squares (v)));
 523   return SCM_UNSPECIFIED;
 524 }
 525 #undef FUNC_NAME
 526
 527
 528 SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0,
 529             (SCM v, SCM state),
 530             "Fills vect with inexact real random numbers that are\n"
 531             "independent and standard normally distributed\n"
 532             "(i.e., with mean 0 and variance 1).")
 533 #define FUNC_NAME s_scm_random_normal_vector_x
 534 {
 535   size_t n;
 536
 537   if (SCM_UNBNDP (state))
 538     state = SCM_VARIABLE_REF (scm_var_random_state);
 539   SCM_VALIDATE_RSTATE (2, state);
 540   if (SCM_VECTORP (v))
 541     {
 542       n = SCM_VECTOR_LENGTH (v);
 543       while (n-- > 0)
 544         SCM_VECTOR_SET (v, n,
 545                         scm_from_double (scm_c_normal01 (SCM_RSTATE (state))));
 546     }
 547   else
 548     {
 549       /* must be a f64vector. */
 550       double *elts = scm_f64vector_elements (v);
 551       n = scm_c_uniform_vector_length (v);
 552       while (n-- > 0)
 553         elts[n] = scm_c_normal01 (SCM_RSTATE (state));
 554       scm_uniform_vector_release (v);
 555     }
 556   return SCM_UNSPECIFIED;
 557 }
 558 #undef FUNC_NAME
 559
 560 #endif /* SCM_HAVE_ARRAYS */
 561
 562 SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0,
 563             (SCM state),
 564             "Return an inexact real in an exponential distribution with mean\n"
 565             "1.  For an exponential distribution with mean u use (* u\n"
 566             "(random:exp)).")
 567 #define FUNC_NAME s_scm_random_exp
 568 {
 569   if (SCM_UNBNDP (state))
 570     state = SCM_VARIABLE_REF (scm_var_random_state);
 571   SCM_VALIDATE_RSTATE (1, state);
 572   return scm_from_double (scm_c_exp1 (SCM_RSTATE (state)));
 573 }
 574 #undef FUNC_NAME
 575
 576 void
 577 scm_init_random ()
 578 {
 579   int i, m;
 580   /* plug in default RNG */
 581   scm_t_rng rng =
 582   {
 583     sizeof (scm_t_i_rstate),
 584     (unsigned long (*)()) scm_i_uniform32,
 585     (void (*)())          scm_i_init_rstate,
 586     (scm_t_rstate *(*)())    scm_i_copy_rstate
 587   };
 588   scm_the_rng = rng;
 589
 590   scm_tc16_rstate = scm_make_smob_type ("random-state", 0);
 591   scm_set_smob_free (scm_tc16_rstate, rstate_free);
 592
 593   for (m = 1; m <= 0x100; m <<= 1)
 594     for (i = m >> 1; i < m; ++i)
 595       scm_masktab[i] = m - 1;
 596
 597 #include "libguile/random.x"
 598
 599   scm_add_feature ("random");
 600 }
 601
 602 /*
 603   Local Variables:
 604   c-file-style: "gnu"
 605   End:
 606 */