src/code/target-random.lisp

   1 ;;;; This implementation of RANDOM is based on the Mersenne Twister random
   2 ;;;; number generator "MT19937" due to Matsumoto and Nishimura. See:
   3 ;;;;   Makoto Matsumoto and T. Nishimura, "Mersenne twister: A
   4 ;;;;   623-dimensionally equidistributed uniform pseudorandom number
   5 ;;;;   generator.", ACM Transactions on Modeling and Computer Simulation,
   6 ;;;;   1997, to appear.
   7
   8 ;;;; This software is part of the SBCL system. See the README file for
   9 ;;;; more information.
  10 ;;;;
  11 ;;;; This software is derived from the CMU CL system, which was
  12 ;;;; written at Carnegie Mellon University and released into the
  13 ;;;; public domain. The software is in the public domain and is
  14 ;;;; provided with absolutely no warranty. See the COPYING and CREDITS
  15 ;;;; files for more information.
  16
  17 (in-package "SB!KERNEL")
  18 \f
  19 ;;;; RANDOM-STATEs
  20
  21 (def!method make-load-form ((random-state random-state) &optional environment)
  22   (make-load-form-saving-slots random-state :environment environment))
  23
  24 (def!method print-object ((state random-state) stream)
  25   (if (and *print-readably* (not *read-eval*))
  26       (error 'print-not-readable :object state)
  27       (format stream "#S(~S ~S #.~S)"
  28               'random-state
  29               ':state
  30               `(make-array 627
  31                 :element-type
  32                 '(unsigned-byte 32)
  33                 :initial-contents
  34                 ',(coerce (random-state-state state) 'list)))))
  35
  36 ;;; The state is stored in a (simple-array (unsigned-byte 32) (627))
  37 ;;; wrapped in a random-state structure:
  38 ;;;
  39 ;;;  0-1:   Constant matrix A. [0, #x9908b0df]
  40 ;;;  2:     Index k.
  41 ;;;  3-626: State.
  42
  43 ;;; Generate and initialize a new random-state array. Index is
  44 ;;; initialized to 1 and the states to 32bit integers excluding zero.
  45 ;;;
  46 ;;; Seed - A 32bit number, not zero.
  47 ;;;
  48 ;;; Apparently uses the generator Line 25 of Table 1 in
  49 ;;; [KNUTH 1981, The Art of Computer Programming, Vol. 2 (2nd Ed.), pp102]
  50 (defun init-random-state (&optional (seed 4357) state)
  51   (declare (type (integer 1 #xffffffff) seed))
  52   (let ((state (or state (make-array 627 :element-type '(unsigned-byte 32)))))
  53     (declare (type (simple-array (unsigned-byte 32) (627)) state))
  54     (setf (aref state 1) #x9908b0df)
  55     (setf (aref state 2) 1)
  56     (setf (aref state 3) seed)
  57     (do ((k 1 (1+ k)))
  58         ((>= k 624))
  59       (declare (type (mod 625) k))
  60       (setf (aref state (+ 3 k))
  61             (logand (* 69069 (aref state (+ 3 (1- k)))) #xffffffff)))
  62     state))
  63
  64 (defvar *random-state*)
  65 (defun !random-cold-init ()
  66   (/show0 "entering !RANDOM-COLD-INIT")
  67   (setf *random-state* (%make-random-state))
  68   (/show0 "returning from !RANDOM-COLD-INIT"))
  69
  70 (defun make-random-state (&optional state)
  71   #!+sb-doc
  72   "Make a RANDOM-STATE object. If STATE is not supplied, return a copy
  73   of the default random state. If STATE is a random state, then return a
  74   copy of it. If STATE is T then return a random state generated from
  75   the universal time."
  76   (/show0 "entering MAKE-RANDOM-STATE")
  77   (flet ((copy-random-state (state)
  78            (/show0 "entering COPY-RANDOM-STATE")
  79            (let ((state (random-state-state state))
  80                  (new-state
  81                   (make-array 627 :element-type '(unsigned-byte 32))))
  82              (/show0 "made NEW-STATE, about to DOTIMES")
  83              (dotimes (i 627)
  84                (setf (aref new-state i) (aref state i)))
  85              (/show0 "falling through to %MAKE-RANDOM-STATE")
  86              (%make-random-state :state new-state))))
  87     (/show0 "at head of ETYPECASE in MAKE-RANDOM-STATE")
  88     (etypecase state
  89       (null
  90        (/show0 "NULL case")
  91        (copy-random-state *random-state*))
  92       (random-state
  93        (/show0 "RANDOM-STATE-P clause")
  94        (copy-random-state state))
  95       ((member t)
  96        (/show0 "T clause")
  97        (%make-random-state :state (init-random-state
  98                                    (logand (get-universal-time)
  99                                            #xffffffff)))))))
 100 \f
 101 ;;;; random entries
 102
 103 ;;; This function generates a 32bit integer between 0 and #xffffffff
 104 ;;; inclusive.
 105 (defconstant n-random-mt19937-bits 32)
 106 (declaim (inline %random32))
 107 ;;; portable implementation
 108 (defconstant mt19937-n 624)
 109 (defconstant mt19937-m 397)
 110 (defconstant mt19937-upper-mask #x80000000)
 111 (defconstant mt19937-lower-mask #x7fffffff)
 112 (defconstant mt19937-b #x9D2C5680)
 113 (defconstant mt19937-c #xEFC60000)
 114 #!-x86
 115 (defun random-mt19937-update (state)
 116   (declare (type (simple-array (unsigned-byte 32) (627)) state)
 117            (optimize (speed 3) (safety 0)))
 118   (let ((y 0))
 119     (declare (type (unsigned-byte 32) y))
 120     (do ((kk 3 (1+ kk)))
 121         ((>= kk (+ 3 (- mt19937-n mt19937-m))))
 122       (declare (type (mod 628) kk))
 123       (setf y (logior (logand (aref state kk) mt19937-upper-mask)
 124                       (logand (aref state (1+ kk)) mt19937-lower-mask)))
 125       (setf (aref state kk) (logxor (aref state (+ kk mt19937-m))
 126                                     (ash y -1) (aref state (logand y 1)))))
 127     (do ((kk (+ (- mt19937-n mt19937-m) 3) (1+ kk)))
 128         ((>= kk (+ (1- mt19937-n) 3)))
 129       (declare (type (mod 628) kk))
 130       (setf y (logior (logand (aref state kk) mt19937-upper-mask)
 131                       (logand (aref state (1+ kk)) mt19937-lower-mask)))
 132       (setf (aref state kk) (logxor (aref state (+ kk (- mt19937-m mt19937-n)))
 133                                     (ash y -1) (aref state (logand y 1)))))
 134     (setf y (logior (logand (aref state (+ 3 (1- mt19937-n)))
 135                             mt19937-upper-mask)
 136                     (logand (aref state 3) mt19937-lower-mask)))
 137     (setf (aref state (+ 3 (1- mt19937-n)))
 138           (logxor (aref state (+ 3 (1- mt19937-m)))
 139                   (ash y -1) (aref state (logand y 1)))))
 140   (values))
 141 #!-x86
 142 (defun %random32 (state)
 143   (declare (type random-state state))
 144   (let* ((state (random-state-state state))
 145          (k (aref state 2)))
 146     (declare (type (mod 628) k))
 147     (when (= k mt19937-n)
 148       (random-mt19937-update state)
 149       (setf k 0))
 150     (setf (aref state 2) (1+ k))
 151     (let ((y (aref state (+ 3 k))))
 152       (declare (type (unsigned-byte 32) y))
 153       (setf y (logxor y (ash y -11)))
 154       (setf y (logxor y (ash (logand y (ash mt19937-b -7)) 7)))
 155       (setf y (logxor y (ash (logand y (ash mt19937-c -15)) 15)))
 156       (setf y (logxor y (ash y -18)))
 157       y)))
 158
 159 ;;; Using inline VOP support, only available on the x86 so far.
 160 ;;;
 161 ;;; FIXME: It would be nice to have some benchmark numbers on this.
 162 ;;; My inclination is to get rid of the nonportable implementation
 163 ;;; unless the performance difference is just enormous.
 164 #!+x86
 165 (defun %random32 (state)
 166   (declare (type random-state state))
 167   (sb!vm::random-mt19937 (random-state-state state)))
 168
 169 (declaim (inline %random-word))
 170 (defun %random-word (state)
 171   ;; KLUDGE: This #.(ECASE ...) is not the most flexible and elegant
 172   ;; construct one could imagine. It is intended as a quick fix to stand
 173   ;; in for The Right Thing which can't be coded so quickly.
 174   ;;
 175   ;; The Right Thing: The Mersenne Twister has been generalized to 64
 176   ;; bits, and seems likely to be generalized to any future common CPU
 177   ;; word width as well. Thus, it should be straightforward to
 178   ;; implement this as "Return one sample from the MT variant
 179   ;; corresponding to SB!VM:N-WORD-BITS."
 180   ;;
 181   ;; Meanwhile: Mock that up by pasting together as many samples from
 182   ;; the 32-bit Mersenne Twister as necessary.
 183   #.(ecase sb!vm:n-word-bits
 184       (32 '(%random32 state))
 185       (64 '(logior
 186             (%random32 state)
 187             (ash (%random32 state) 32)))))
 188 \f
 189 ;;; Handle the single or double float case of RANDOM. We generate a
 190 ;;; float between 0.0 and 1.0 by clobbering the significand of 1.0
 191 ;;; with random bits, then subtracting 1.0. This hides the fact that
 192 ;;; we have a hidden bit.
 193 #!-sb-fluid (declaim (inline %random-single-float %random-double-float))
 194 (declaim (ftype (function ((single-float (0f0)) random-state)
 195                           (single-float 0f0))
 196                 %random-single-float))
 197 (defun %random-single-float (arg state)
 198   (declare (type (single-float (0f0)) arg)
 199            (type random-state state))
 200   ;; KLUDGE: The hardwired-to-32 hackery here could be replaced by a
 201   ;; call to %RANDOM-BITS if there were sufficiently smart
 202   ;; DEFTRANSFORMs for it, but in sbcl-1.0.4.epsilon it looks as
 203   ;; though it would be a performance disaster.
 204   (aver (<= sb!vm:single-float-digits 32))
 205   (* arg
 206      (- (make-single-float
 207          (dpb (ash
 208                (%random32 state)
 209                (- sb!vm:single-float-digits 32))
 210               sb!vm:single-float-significand-byte
 211               (single-float-bits 1.0)))
 212         1.0)))
 213 (declaim (ftype (function ((double-float (0d0)) random-state)
 214                           (double-float 0d0))
 215                 %random-double-float))
 216
 217 ;;; 53-bit version
 218 #!-x86
 219 (defun %random-double-float (arg state)
 220   (declare (type (double-float (0d0)) arg)
 221            (type random-state state))
 222   ;; KLUDGE: As in %RANDOM-SIMNGLE-FLOAT, as of sbcl-1.0.4.epsilon
 223   ;; calling %RANDOM-BITS doesn't look reasonable, so we bang bits.
 224   (aver (<= sb!vm:single-float-digits 32))
 225   (* arg
 226      (- (sb!impl::make-double-float
 227          (dpb (ash (%random32 state)
 228                    (- sb!vm:double-float-digits 64))
 229               sb!vm:double-float-significand-byte
 230               (sb!impl::double-float-high-bits 1d0))
 231          (%random32 state))
 232         1d0)))
 233
 234 ;;; using a faster inline VOP
 235 #!+x86
 236 (defun %random-double-float (arg state)
 237   (declare (type (double-float (0d0)) arg)
 238            (type random-state state))
 239   (let ((state-vector (random-state-state state)))
 240     (* arg
 241        (- (sb!impl::make-double-float
 242            (dpb (ash (sb!vm::random-mt19937 state-vector)
 243                      (- sb!vm:double-float-digits
 244                         n-random-mt19937-bits
 245                         sb!vm:n-word-bits))
 246                 sb!vm:double-float-significand-byte
 247                 (sb!impl::double-float-high-bits 1d0))
 248            (sb!vm::random-mt19937 state-vector))
 249           1d0))))
 250 \f
 251 ;;;; random integers
 252
 253 ;;; a bitmask M wide enough that (= (LOGAND INCLUSIVE-LIMIT M) INCLUSIVE-LIMIT)
 254 (declaim (inline %inclusive-random-integer-mask))
 255 (defun %inclusive-random-integer-mask (inclusive-limit)
 256   (1- (ash 1 (integer-length inclusive-limit))))
 257
 258 ;;; Transform a uniform sample from an at-least-as-large range into a
 259 ;;; random sample in [0,INCLUSIVE-LIMIT] range by throwing away
 260 ;;; too-big samples.
 261 ;;;
 262 ;;; Up through sbcl-1.0.4, the (RANDOM INTEGER) worked by taking (MOD
 263 ;;; RAW-MERSENNE-OUTPUT INTEGER). That introduces enough bias that it
 264 ;;; is unsuitable for some calculations (like the Metropolis Monte
 265 ;;; Carlo simulations that I, WHN, worked on in grad school): in the
 266 ;;; sbcl-1.0.4, the expectation value of a sample was (in the worst
 267 ;;; part of the range) more than 1.0001 times the ideal expectation
 268 ;;; value. Perhaps that was even ANSI-conformant, because ANSI says only
 269 ;;;   An approximately uniform choice distribution is used. If
 270 ;;;   LIMIT is an integer, each of the possible results occurs
 271 ;;;   with (approximate) probability 1/LIMIT.
 272 ;;; But I (WHN again) said "yuck," so these days we try to get a
 273 ;;; truly uniform distribution by translating RAW-MERSENNE-OUTPUT to
 274 ;;; our output using the second method recommended by
 275 ;;; <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/efaq.html>:
 276 ;;;   In the rare case that the rounding-off error becomes a problem,
 277 ;;;   obtain minimum integer n such that N<=2^n, generate integer
 278 ;;;   random numbers, take the most significant n bits, and discard
 279 ;;;   those more than or equal to N.
 280 ;;; (The first method recommended there differs from the sbcl-1.0.4
 281 ;;; method: the recommended method gives slightly different weights
 282 ;;; to different integers, but distributes the overweighted integers
 283 ;;; evenly across the range of outputs so that any skew would be
 284 ;;; of order (/ MOST-POSITIVE-MACHINE-WORD), rather than the (/ 10000)
 285 ;;; or so achieved by sbcl-1.0.4. That skew would probably be
 286 ;;; acceptable --- it seems to be exactly the kind of deviation that
 287 ;;; might have been anticipated in the ANSI CL standard. However, that
 288 ;;; recommended calculation involves multiplication and division of
 289 ;;; two-machine-word bignums, which is hard for us to do efficiently
 290 ;;; here. Without special low-level hacking to support such short
 291 ;;; bignums without consing, the accept-reject solution is not only
 292 ;;; more exact, but also likely more efficient.)
 293 (defmacro %inclusive-random-integer-accept-reject (raw-mersenne-output-expr
 294                                                    inclusive-limit)
 295   (with-unique-names (raw-mersenne-output inclusive-limit-once)
 296   `(loop
 297     with ,inclusive-limit-once = ,inclusive-limit
 298     for ,raw-mersenne-output = ,raw-mersenne-output-expr
 299     until (<= ,raw-mersenne-output ,inclusive-limit-once)
 300     finally (return ,raw-mersenne-output))))
 301
 302 ;;; an UNSIGNED-BYTE of N-WORDS words sampled from the Mersenne twister
 303 (declaim (inline %random-words))
 304 (defun %random-words (n-words state)
 305   ;; KLUDGE: This algorithm will cons O(N^2) words when constructing
 306   ;; an N-word result. To do better should be fundamentally easy, we
 307   ;; just need to do some low-level hack preallocating the bignum and
 308   ;; writing its words one by one.
 309   ;;
 310   ;; Note: The old sbcl-1.0.4 code did its analogous operation using a
 311   ;; mysterious RANDOM-INTEGER-OVERLAP parameter, "the amount that we
 312   ;; overlap chunks by when building a large integer to make up for
 313   ;; the loss of randomness in the low bits." No such thing is called
 314   ;; for on
 315   ;;   <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/efaq.html>:
 316   ;; they say there that it's OK just to concatenate words of twister
 317   ;; output with no overlap. Thus, crossing our fingers and hoping that
 318   ;; the previous RNG author didn't have some subtle reason to need
 319   ;; RANDOM-INTEGER-OVERLAP that we know not of, we just concatenate
 320   ;; chunks.
 321   (loop repeat n-words
 322         for result = 0 then (logior (ash result sb!vm:n-word-bits)
 323                                     (%random-word state))
 324         finally (return result)))
 325
 326 ;;; an UNSIGNED-BYTE of N-BITS bits sampled from the Mersenne twister
 327 (declaim (inline %random-bits))
 328 (defun %random-bits (n-bits state)
 329   (multiple-value-bind (n-full-words n-extra-bits)
 330       (floor n-bits sb!vm:n-word-bits)
 331     (let ((full-chunks (%random-words n-full-words state)))
 332       (if (zerop n-extra-bits)
 333           full-chunks
 334           (logior full-chunks
 335                   (ash (logand (%random-word state)
 336                                (1- (ash 1 n-extra-bits)))
 337                        (* n-full-words
 338                           sb!vm:n-word-bits)))))))
 339
 340 ;;; the guts of (RANDOM (1+ INCLUSIVE-LIMIT))
 341 (defun %inclusive-random-integer (inclusive-limit state)
 342   (declare (optimize speed (space 1))) ; to ensure DEFTRANSFORM is enabled
 343   (%inclusive-random-integer inclusive-limit state))
 344
 345 ;;; the guts of RANDOM for the known-to-return-FIXNUM case
 346 ;;;
 347 ;;; We use INCLUSIVE-LIMIT instead of the (exclusive) LIMIT of RANDOM,
 348 ;;; because we want it to be a FIXNUM even for the possibly-common
 349 ;;; case of (RANDOM (1+ MOST-POSITIVE-FIXNUM)). (That case is what
 350 ;;; one might use for generating random hash values, e.g.)
 351 ;;; It also turns out to be just what's needed for INTEGER-LENGTH.
 352 (declaim (maybe-inline %inclusive-random-fixnum))
 353 (defun %inclusive-random-fixnum (inclusive-limit state)
 354   (declare (type (and fixnum unsigned-byte) inclusive-limit))
 355   (let (;; If this calculation needs to be optimized further, a good
 356         ;; start might be a DEFTRANSFORM which picks off the case of
 357         ;; constant LIMIT and precomputes the MASK at compile time.
 358         (mask (%inclusive-random-integer-mask inclusive-limit)))
 359     (%inclusive-random-integer-accept-reject (logand (%random-word state) mask)
 360                                              inclusive-limit)))
 361 \f
 362 ;;;; outer, dynamically-typed interface
 363
 364 (defun random (arg &optional (state *random-state*))
 365   (declare (inline %random-single-float
 366                    %random-double-float
 367                    #!+long-float %random-long-float))
 368   (cond
 369     ((and (fixnump arg) (plusp arg))
 370      (locally
 371          ;; The choice to inline this very common case of
 372          ;; %INCLUSIVE-RANDOM-FIXNUM and not the less-common call
 373          ;; below was just a guess (by WHN 2007-03-27), not based on
 374          ;; benchmarking or anything.
 375          (declare (inline %inclusive-random-fixnum))
 376        (%inclusive-random-fixnum (1- arg) state)))
 377     ((and (typep arg 'single-float) (> arg 0.0f0))
 378      (%random-single-float arg state))
 379     ((and (typep arg 'double-float) (> arg 0.0d0))
 380      (%random-double-float arg state))
 381     #!+long-float
 382     ((and (typep arg 'long-float) (> arg 0.0l0))
 383      (%random-long-float arg state))
 384     ((and (integerp arg) (plusp arg))
 385      (if (= arg (1+ most-positive-fixnum))
 386          (%inclusive-random-fixnum most-positive-fixnum state)
 387          (%inclusive-random-integer (1- arg) state)))
 388     (t
 389      (error 'simple-type-error
 390             :expected-type '(or (integer 1) (float (0))) :datum arg
 391             :format-control "~@<Argument is neither a positive integer nor a ~
 392                              positive float: ~2I~_~S~:>"
 393             :format-arguments (list arg)))))