libguile/random.c

   1 /* Copyright (C) 1999,2000,2001 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)
   5  * any later version.
   6  *
   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.
  11  *
  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
  16  *
  17  * As a special exception, the Free Software Foundation gives permission
  18  * for additional uses of the text contained in its release of GUILE.
  19  *
  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.
  25  *
  26  * This exception does not however invalidate any other reasons why
  27  * the executable file might be covered by the GNU General Public License.
  28  *
  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.
  36  *
  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.  */
  40
  41 /* Software engineering face-lift by Greg J. Badros, 11-Dec-1999,
  42    gjb@cs.washington.edu, http://www.cs.washington.edu/homes/gjb */
  43
  44
  45 /* Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */
  46
  47 #include "libguile/_scm.h"
  48
  49 #include <stdio.h>
  50 #include <math.h>
  51 #include <string.h>
  52 #include "libguile/smob.h"
  53 #include "libguile/numbers.h"
  54 #include "libguile/feature.h"
  55 #include "libguile/strings.h"
  56 #include "libguile/unif.h"
  57 #include "libguile/vectors.h"
  58
  59 #include "libguile/validate.h"
  60 #include "libguile/random.h"
  61
  62 \f
  63 /*
  64  * A plugin interface for RNGs
  65  *
  66  * Using this interface, it is possible for the application to tell
  67  * libguile to use a different RNG.  This is desirable if it is
  68  * necessary to use the same RNG everywhere in the application in
  69  * order to prevent interference, if the application uses RNG
  70  * hardware, or if the application has special demands on the RNG.
  71  *
  72  * Look in random.h and how the default generator is "plugged in" in
  73  * scm_init_random().
  74  */
  75
  76 scm_rng scm_the_rng;
  77
  78 \f
  79 /*
  80  * The prepackaged RNG
  81  *
  82  * This is the MWC (Multiply With Carry) random number generator
  83  * described by George Marsaglia at the Department of Statistics and
  84  * Supercomputer Computations Research Institute, The Florida State
  85  * University (http://stat.fsu.edu/~geo).
  86  *
  87  * It uses 64 bits, has a period of 4578426017172946943 (4.6e18), and
  88  * passes all tests in the DIEHARD test suite
  89  * (http://stat.fsu.edu/~geo/diehard.html)
  90  */
  91
  92 #define A 2131995753UL
  93
  94 #if SIZEOF_LONG > 4
  95 #if SIZEOF_INT > 4
  96 #define LONG32 unsigned short
  97 #else
  98 #define LONG32 unsigned int
  99 #endif
 100 #define LONG64 unsigned long
 101 #else
 102 #define LONG32 unsigned long
 103 #define LONG64 unsigned long long
 104 #endif
 105
 106 #if SIZEOF_LONG > 4 || defined (HAVE_LONG_LONGS)
 107
 108 unsigned long
 109 scm_i_uniform32 (scm_i_rstate *state)
 110 {
 111   LONG64 x = (LONG64) A * state->w + state->c;
 112   LONG32 w = x & 0xffffffffUL;
 113   state->w = w;
 114   state->c = x >> 32L;
 115   return w;
 116 }
 117
 118 #else
 119
 120 /*     ww  This is a portable version of the same RNG without 64 bit
 121  *   * aa  arithmetic.
 122  *   ----
 123  *     xx  It is only intended to provide identical behaviour on
 124  *    xx   platforms without 8 byte longs or long longs until
 125  *    xx   someone has implemented the routine in assembler code.
 126  *   xxcc
 127  *   ----
 128  *   ccww
 129  */
 130
 131 #define L(x) ((x) & 0xffff)
 132 #define H(x) ((x) >> 16)
 133
 134 unsigned long
 135 scm_i_uniform32 (scm_i_rstate *state)
 136 {
 137   LONG32 x1 = L (A) * L (state->w);
 138   LONG32 x2 = L (A) * H (state->w);
 139   LONG32 x3 = H (A) * L (state->w);
 140   LONG32 w = L (x1) + L (state->c);
 141   LONG32 m = H (x1) + L (x2) + L (x3) + H (state->c) + H (w);
 142   LONG32 x4 = H (A) * H (state->w);
 143   state->w = w = (L (m) << 16) + L (w);
 144   state->c = H (x2) + H (x3) + x4 + H (m);
 145   return w;
 146 }
 147
 148 #endif
 149
 150 void
 151 scm_i_init_rstate (scm_i_rstate *state, char *seed, int n)
 152 {
 153   LONG32 w = 0L;
 154   LONG32 c = 0L;
 155   int i, m;
 156   for (i = 0; i < n; ++i)
 157     {
 158       m = i % 8;
 159       if (m < 4)
 160         w += seed[i] << (8 * m);
 161       else
 162         c += seed[i] << (8 * (m - 4));
 163     }
 164   if ((w == 0 && c == 0) || (w == 0xffffffffUL && c == A - 1))
 165     ++c;
 166   state->w = w;
 167   state->c = c;
 168 }
 169
 170 scm_i_rstate *
 171 scm_i_copy_rstate (scm_i_rstate *state)
 172 {
 173   scm_rstate *new_state = malloc (scm_the_rng.rstate_size);
 174   if (new_state == 0)
 175     scm_memory_error ("rstate");
 176   return memcpy (new_state, state, scm_the_rng.rstate_size);
 177 }
 178
 179 \f
 180 /*
 181  * Random number library functions
 182  */
 183
 184 scm_rstate *
 185 scm_c_make_rstate (char *seed, int n)
 186 {
 187   scm_rstate *state = malloc (scm_the_rng.rstate_size);
 188   if (state == 0)
 189     scm_memory_error ("rstate");
 190   state->reserved0 = 0;
 191   scm_the_rng.init_rstate (state, seed, n);
 192   return state;
 193 }
 194
 195
 196 scm_rstate *
 197 scm_c_default_rstate ()
 198 #define FUNC_NAME "scm_c_default_rstate"
 199 {
 200   SCM state = SCM_CDR (scm_var_random_state);
 201   if (!SCM_RSTATEP (state))
 202     SCM_MISC_ERROR ("*random-state* contains bogus random state", SCM_EOL);
 203   return SCM_RSTATE (state);
 204 }
 205 #undef FUNC_NAME
 206
 207
 208 inline double
 209 scm_c_uniform01 (scm_rstate *state)
 210 {
 211   double x = (double) scm_the_rng.random_bits (state) / (double) 0xffffffffUL;
 212   return ((x + (double) scm_the_rng.random_bits (state))
 213           / (double) 0xffffffffUL);
 214 }
 215
 216 double
 217 scm_c_normal01 (scm_rstate *state)
 218 {
 219   if (state->reserved0)
 220     {
 221       state->reserved0 = 0;
 222       return state->reserved1;
 223     }
 224   else
 225     {
 226       double r, a, n;
 227
 228       r = sqrt (-2.0 * log (scm_c_uniform01 (state)));
 229       a = 2.0 * M_PI * scm_c_uniform01 (state);
 230
 231       n = r * sin (a);
 232       state->reserved1 = r * cos (a);
 233       state->reserved0 = 1;
 234
 235       return n;
 236     }
 237 }
 238
 239 double
 240 scm_c_exp1 (scm_rstate *state)
 241 {
 242   return - log (scm_c_uniform01 (state));
 243 }
 244
 245 unsigned char scm_masktab[256];
 246
 247 unsigned long
 248 scm_c_random (scm_rstate *state, unsigned long m)
 249 {
 250   unsigned int r, mask;
 251   mask = (m < 0x100
 252           ? scm_masktab[m]
 253           : (m < 0x10000
 254              ? scm_masktab[m >> 8] << 8 | 0xff
 255              : (m < 0x1000000
 256                 ? scm_masktab[m >> 16] << 16 | 0xffff
 257                 : scm_masktab[m >> 24] << 24 | 0xffffff)));
 258   while ((r = scm_the_rng.random_bits (state) & mask) >= m);
 259   return r;
 260 }
 261
 262 SCM
 263 scm_c_random_bignum (scm_rstate *state, SCM m)
 264 {
 265   SCM b;
 266   int i, nd;
 267   LONG32 *bits, mask, w;
 268   nd = SCM_NUMDIGS (m);
 269   /* calculate mask for most significant digit */
 270 #if SIZEOF_INT == 4
 271   /* 16 bit digits */
 272   if (nd & 1)
 273     {
 274       /* fix most significant 16 bits */
 275       unsigned short s = SCM_BDIGITS (m)[nd - 1];
 276       mask = s < 0x100 ? scm_masktab[s] : scm_masktab[s >> 8] << 8 | 0xff;
 277     }
 278   else
 279 #endif
 280     {
 281       /* fix most significant 32 bits */
 282 #if SIZEOF_INT == 4
 283       w = SCM_BDIGITS (m)[nd - 1] << 16 | SCM_BDIGITS (m)[nd - 2];
 284 #else
 285       w = SCM_BDIGITS (m)[nd - 1];
 286 #endif
 287       mask = (w < 0x10000
 288               ? (w < 0x100
 289                  ? scm_masktab[w]
 290                  : scm_masktab[w >> 8] << 8 | 0xff)
 291               : (w < 0x1000000
 292                  ? scm_masktab[w >> 16] << 16 | 0xffff
 293                  : scm_masktab[w >> 24] << 24 | 0xffffff));
 294     }
 295   b = scm_mkbig (nd, 0);
 296   bits = (LONG32 *) SCM_BDIGITS (b);
 297   do
 298     {
 299       i = nd;
 300       /* treat most significant digit specially */
 301 #if SIZEOF_INT == 4
 302       /* 16 bit digits */
 303       if (i & 1)
 304         {
 305           ((SCM_BIGDIG*) bits)[i - 1] = scm_the_rng.random_bits (state) & mask;
 306           i /= 2;
 307         }
 308       else
 309 #endif
 310         {
 311           /* fix most significant 32 bits */
 312 #if SIZEOF_INT == 4
 313           w = scm_the_rng.random_bits (state) & mask;
 314           ((SCM_BIGDIG*) bits)[i - 2] = w & 0xffff;
 315           ((SCM_BIGDIG*) bits)[i - 1] = w >> 16;
 316           i = i / 2 - 1;
 317 #else
 318           i /= 2;
 319           bits[--i] = scm_the_rng.random_bits (state) & mask;
 320 #endif
 321         }
 322       /* now fill up the rest of the bignum */
 323       while (i)
 324         bits[--i] = scm_the_rng.random_bits (state);
 325       b = scm_normbig (b);
 326       if (SCM_INUMP (b))
 327         return b;
 328     } while (scm_bigcomp (b, m) <= 0);
 329   return b;
 330 }
 331
 332 /*
 333  * Scheme level representation of random states.
 334  */
 335
 336 scm_bits_t scm_tc16_rstate;
 337
 338 static SCM
 339 make_rstate (scm_rstate *state)
 340 {
 341   SCM_RETURN_NEWSMOB (scm_tc16_rstate, state);
 342 }
 343
 344 static scm_sizet
 345 rstate_free (SCM rstate)
 346 {
 347   free (SCM_RSTATE (rstate));
 348   return scm_the_rng.rstate_size;
 349 }
 350
 351 /*
 352  * Scheme level interface.
 353  */
 354
 355 SCM_GLOBAL_VARIABLE_INIT (scm_var_random_state, "*random-state*", scm_seed_to_random_state (scm_makfrom0str ("URL:http://stat.fsu.edu/~geo/diehard.html")));
 356
 357 SCM_DEFINE (scm_random, "random", 1, 1, 0,
 358             (SCM n, SCM state),
 359             "Return a number in [0,N).\n"
 360             "\n"
 361             "Accepts a positive integer or real n and returns a \n"
 362             "number of the same type between zero (inclusive) and \n"
 363             "N (exclusive). The values returned have a uniform \n"
 364             "distribution.\n"
 365             "\n"
 366             "The optional argument @var{state} must be of the type produced\n"
 367             "by @code{seed->random-state}. It defaults to the value of the\n"
 368             "variable @var{*random-state*}. This object is used to maintain\n"
 369             "the state of the pseudo-random-number generator and is altered\n"
 370             "as a side effect of the random operation.")
 371 #define FUNC_NAME s_scm_random
 372 {
 373   if (SCM_UNBNDP (state))
 374     state = SCM_VARIABLE_REF (scm_var_random_state);
 375   SCM_VALIDATE_RSTATE (2,state);
 376   if (SCM_INUMP (n))
 377     {
 378       unsigned long m = SCM_INUM (n);
 379       SCM_ASSERT_RANGE (1,n,m > 0);
 380       return SCM_MAKINUM (scm_c_random (SCM_RSTATE (state), m));
 381     }
 382   SCM_VALIDATE_NIM (1,n);
 383   if (SCM_REALP (n))
 384     return scm_make_real (SCM_REAL_VALUE (n)
 385                           * scm_c_uniform01 (SCM_RSTATE (state)));
 386   SCM_VALIDATE_SMOB (1, n, big);
 387   return scm_c_random_bignum (SCM_RSTATE (state), n);
 388 }
 389 #undef FUNC_NAME
 390
 391 SCM_DEFINE (scm_copy_random_state, "copy-random-state", 0, 1, 0,
 392             (SCM state),
 393             "Return a copy of the random state @var{state}.")
 394 #define FUNC_NAME s_scm_copy_random_state
 395 {
 396   if (SCM_UNBNDP (state))
 397     state = SCM_VARIABLE_REF (scm_var_random_state);
 398   SCM_VALIDATE_RSTATE (1,state);
 399   return make_rstate (scm_the_rng.copy_rstate (SCM_RSTATE (state)));
 400 }
 401 #undef FUNC_NAME
 402
 403 SCM_DEFINE (scm_seed_to_random_state, "seed->random-state", 1, 0, 0,
 404             (SCM seed),
 405             "Return a new random state using @var{seed}.")
 406 #define FUNC_NAME s_scm_seed_to_random_state
 407 {
 408   if (SCM_NUMBERP (seed))
 409     seed = scm_number_to_string (seed, SCM_UNDEFINED);
 410   SCM_VALIDATE_STRING (1,seed);
 411   return make_rstate (scm_c_make_rstate (SCM_STRING_CHARS (seed),
 412                                          SCM_STRING_LENGTH (seed)));
 413 }
 414 #undef FUNC_NAME
 415
 416 SCM_DEFINE (scm_random_uniform, "random:uniform", 0, 1, 0,
 417             (SCM state),
 418             "Return a uniformly distributed inexact real random number in\n"
 419             "[0,1).")
 420 #define FUNC_NAME s_scm_random_uniform
 421 {
 422   if (SCM_UNBNDP (state))
 423     state = SCM_VARIABLE_REF (scm_var_random_state);
 424   SCM_VALIDATE_RSTATE (1,state);
 425   return scm_make_real (scm_c_uniform01 (SCM_RSTATE (state)));
 426 }
 427 #undef FUNC_NAME
 428
 429 SCM_DEFINE (scm_random_normal, "random:normal", 0, 1, 0,
 430             (SCM state),
 431             "Return an inexact real in a normal distribution.  The\n"
 432             "distribution used has mean 0 and standard deviation 1.  For a\n"
 433             "normal distribution with mean m and standard deviation d use\n"
 434             "@code{(+ m (* d (random:normal)))}.")
 435 #define FUNC_NAME s_scm_random_normal
 436 {
 437   if (SCM_UNBNDP (state))
 438     state = SCM_VARIABLE_REF (scm_var_random_state);
 439   SCM_VALIDATE_RSTATE (1,state);
 440   return scm_make_real (scm_c_normal01 (SCM_RSTATE (state)));
 441 }
 442 #undef FUNC_NAME
 443
 444 #ifdef HAVE_ARRAYS
 445
 446 static void
 447 vector_scale (SCM v, double c)
 448 {
 449   int n = SCM_INUM (scm_uniform_vector_length (v));
 450   if (SCM_VECTORP (v))
 451     while (--n >= 0)
 452       SCM_REAL_VALUE (SCM_VELTS (v)[n]) *= c;
 453   else
 454     while (--n >= 0)
 455       ((double *) SCM_VELTS (v))[n] *= c;
 456 }
 457
 458 static double
 459 vector_sum_squares (SCM v)
 460 {
 461   double x, sum = 0.0;
 462   int n = SCM_INUM (scm_uniform_vector_length (v));
 463   if (SCM_VECTORP (v))
 464     while (--n >= 0)
 465       {
 466         x = SCM_REAL_VALUE (SCM_VELTS (v)[n]);
 467         sum += x * x;
 468       }
 469   else
 470     while (--n >= 0)
 471       {
 472         x = ((double *) SCM_VELTS (v))[n];
 473         sum += x * x;
 474       }
 475   return sum;
 476 }
 477
 478 /* For the uniform distribution on the solid sphere, note that in
 479  * this distribution the length r of the vector has cumulative
 480  * distribution r^n; i.e., u=r^n is uniform [0,1], so r can be
 481  * generated as r=u^(1/n).
 482  */
 483 SCM_DEFINE (scm_random_solid_sphere_x, "random:solid-sphere!", 1, 1, 0,
 484             (SCM v, SCM state),
 485             "Fills vect with inexact real random numbers\n"
 486             "the sum of whose squares is less than 1.0.\n"
 487             "Thinking of vect as coordinates in space of \n"
 488             "dimension n = (vector-length vect), the coordinates \n"
 489             "are uniformly distributed within the unit n-shere.\n"
 490             "The sum of the squares of the numbers is returned.")
 491 #define FUNC_NAME s_scm_random_solid_sphere_x
 492 {
 493   SCM_VALIDATE_VECTOR_OR_DVECTOR (1,v);
 494   if (SCM_UNBNDP (state))
 495     state = SCM_VARIABLE_REF (scm_var_random_state);
 496   SCM_VALIDATE_RSTATE (2,state);
 497   scm_random_normal_vector_x (v, state);
 498   vector_scale (v,
 499                 pow (scm_c_uniform01 (SCM_RSTATE (state)),
 500                      1.0 / SCM_INUM (scm_uniform_vector_length (v)))
 501                 / sqrt (vector_sum_squares (v)));
 502   return SCM_UNSPECIFIED;
 503 }
 504 #undef FUNC_NAME
 505
 506 SCM_DEFINE (scm_random_hollow_sphere_x, "random:hollow-sphere!", 1, 1, 0,
 507             (SCM v, SCM state),
 508             "Fills vect with inexact real random numbers\n"
 509             "the sum of whose squares is equal to 1.0.\n"
 510             "Thinking of vect as coordinates in space of \n"
 511             "dimension n = (vector-length vect), the coordinates\n"
 512             "are uniformly distributed over the surface of the \n"
 513             "unit n-shere.")
 514 #define FUNC_NAME s_scm_random_hollow_sphere_x
 515 {
 516   SCM_VALIDATE_VECTOR_OR_DVECTOR (1,v);
 517   if (SCM_UNBNDP (state))
 518     state = SCM_VARIABLE_REF (scm_var_random_state);
 519   SCM_VALIDATE_RSTATE (2,state);
 520   scm_random_normal_vector_x (v, state);
 521   vector_scale (v, 1 / sqrt (vector_sum_squares (v)));
 522   return SCM_UNSPECIFIED;
 523 }
 524 #undef FUNC_NAME
 525
 526
 527 SCM_DEFINE (scm_random_normal_vector_x, "random:normal-vector!", 1, 1, 0,
 528             (SCM v, SCM state),
 529             "Fills vect with inexact real random numbers that are\n"
 530             "independent and standard normally distributed\n"
 531             "(i.e., with mean 0 and variance 1).")
 532 #define FUNC_NAME s_scm_random_normal_vector_x
 533 {
 534   int n;
 535   SCM_VALIDATE_VECTOR_OR_DVECTOR (1,v);
 536   if (SCM_UNBNDP (state))
 537     state = SCM_VARIABLE_REF (scm_var_random_state);
 538   SCM_VALIDATE_RSTATE (2,state);
 539   n = SCM_INUM (scm_uniform_vector_length (v));
 540   if (SCM_VECTORP (v))
 541     while (--n >= 0)
 542       SCM_VELTS (v)[n] = scm_make_real (scm_c_normal01 (SCM_RSTATE (state)));
 543   else
 544     while (--n >= 0)
 545       ((double *) SCM_VELTS (v))[n] = scm_c_normal01 (SCM_RSTATE (state));
 546   return SCM_UNSPECIFIED;
 547 }
 548 #undef FUNC_NAME
 549
 550 #endif /* HAVE_ARRAYS */
 551
 552 SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0,
 553             (SCM state),
 554             "Return an inexact real in an exponential distribution with mean\n"
 555             "1.  For an exponential distribution with mean u use (* u\n"
 556             "(random:exp)).")
 557 #define FUNC_NAME s_scm_random_exp
 558 {
 559   if (SCM_UNBNDP (state))
 560     state = SCM_VARIABLE_REF (scm_var_random_state);
 561   SCM_VALIDATE_RSTATE (1,state);
 562   return scm_make_real (scm_c_exp1 (SCM_RSTATE (state)));
 563 }
 564 #undef FUNC_NAME
 565
 566 void
 567 scm_init_random ()
 568 {
 569   int i, m;
 570   /* plug in default RNG */
 571   scm_rng rng =
 572   {
 573     sizeof (scm_i_rstate),
 574     (unsigned long (*)()) scm_i_uniform32,
 575     (void (*)())          scm_i_init_rstate,
 576     (scm_rstate *(*)())    scm_i_copy_rstate
 577   };
 578   scm_the_rng = rng;
 579
 580   scm_tc16_rstate = scm_make_smob_type ("random-state", 0);
 581   scm_set_smob_free (scm_tc16_rstate, rstate_free);
 582
 583   for (m = 1; m <= 0x100; m <<= 1)
 584     for (i = m >> 1; i < m; ++i)
 585       scm_masktab[i] = m - 1;
 586
 587 #ifndef SCM_MAGIC_SNARFER
 588 #include "libguile/random.x"
 589 #endif
 590
 591   /* Check that the assumptions about bits per bignum digit are correct. */
 592 #if SIZEOF_INT == 4
 593   m = 16;
 594 #else
 595   m = 32;
 596 #endif
 597   if (m != SCM_BITSPERDIG)
 598     {
 599       fprintf (stderr, "Internal inconsistency: Confused about bignum digit size in random.c\n");
 600       exit (1);
 601     }
 602
 603   scm_add_feature ("random");
 604 }
 605
 606 /*
 607   Local Variables:
 608   c-file-style: "gnu"
 609   End:
 610 */