libguile/random.c

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