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 #!-sb-fluid (declaim (inline random-chunk))
 106 ;;; portable implementation
 107 (defconstant mt19937-n 624)
 108 (defconstant mt19937-m 397)
 109 (defconstant mt19937-upper-mask #x80000000)
 110 (defconstant mt19937-lower-mask #x7fffffff)
 111 (defconstant mt19937-b #x9D2C5680)
 112 (defconstant mt19937-c #xEFC60000)
 113 #!-x86
 114 (defun random-mt19937-update (state)
 115   (declare (type (simple-array (unsigned-byte 32) (627)) state)
 116            (optimize (speed 3) (safety 0)))
 117   (let ((y 0))
 118     (declare (type (unsigned-byte 32) y))
 119     (do ((kk 3 (1+ kk)))
 120         ((>= kk (+ 3 (- mt19937-n mt19937-m))))
 121       (declare (type (mod 628) kk))
 122       (setf y (logior (logand (aref state kk) mt19937-upper-mask)
 123                       (logand (aref state (1+ kk)) mt19937-lower-mask)))
 124       (setf (aref state kk) (logxor (aref state (+ kk mt19937-m))
 125                                     (ash y -1) (aref state (logand y 1)))))
 126     (do ((kk (+ (- mt19937-n mt19937-m) 3) (1+ kk)))
 127         ((>= kk (+ (1- mt19937-n) 3)))
 128       (declare (type (mod 628) kk))
 129       (setf y (logior (logand (aref state kk) mt19937-upper-mask)
 130                       (logand (aref state (1+ kk)) mt19937-lower-mask)))
 131       (setf (aref state kk) (logxor (aref state (+ kk (- mt19937-m mt19937-n)))
 132                                     (ash y -1) (aref state (logand y 1)))))
 133     (setf y (logior (logand (aref state (+ 3 (1- mt19937-n)))
 134                             mt19937-upper-mask)
 135                     (logand (aref state 3) mt19937-lower-mask)))
 136     (setf (aref state (+ 3 (1- mt19937-n)))
 137           (logxor (aref state (+ 3 (1- mt19937-m)))
 138                   (ash y -1) (aref state (logand y 1)))))
 139   (values))
 140 #!-x86
 141 (defun random-chunk (state)
 142   (declare (type random-state state))
 143   (let* ((state (random-state-state state))
 144          (k (aref state 2)))
 145     (declare (type (mod 628) k))
 146     (when (= k mt19937-n)
 147       (random-mt19937-update state)
 148       (setf k 0))
 149     (setf (aref state 2) (1+ k))
 150     (let ((y (aref state (+ 3 k))))
 151       (declare (type (unsigned-byte 32) y))
 152       (setf y (logxor y (ash y -11)))
 153       (setf y (logxor y (ash (logand y (ash mt19937-b -7)) 7)))
 154       (setf y (logxor y (ash (logand y (ash mt19937-c -15)) 15)))
 155       (setf y (logxor y (ash y -18)))
 156       y)))
 157
 158 ;;; Using inline VOP support, only available on the x86 so far.
 159 ;;;
 160 ;;; FIXME: It would be nice to have some benchmark numbers on this.
 161 ;;; My inclination is to get rid of the nonportable implementation
 162 ;;; unless the performance difference is just enormous.
 163 #!+x86
 164 (defun random-chunk (state)
 165   (declare (type random-state state))
 166   (sb!vm::random-mt19937 (random-state-state state)))
 167 \f
 168 ;;; Handle the single or double float case of RANDOM. We generate a
 169 ;;; float between 0.0 and 1.0 by clobbering the significand of 1.0
 170 ;;; with random bits, then subtracting 1.0. This hides the fact that
 171 ;;; we have a hidden bit.
 172 #!-sb-fluid (declaim (inline %random-single-float %random-double-float))
 173 (declaim (ftype (function ((single-float (0f0)) random-state)
 174                           (single-float 0f0))
 175                 %random-single-float))
 176 (defun %random-single-float (arg state)
 177   (declare (type (single-float (0f0)) arg)
 178            (type random-state state))
 179   (* arg
 180      (- (make-single-float
 181          (dpb (ash (random-chunk state)
 182                    (- sb!vm:single-float-digits n-random-chunk-bits))
 183               sb!vm:single-float-significand-byte
 184               (single-float-bits 1.0)))
 185         1.0)))
 186 (declaim (ftype (function ((double-float (0d0)) random-state)
 187                           (double-float 0d0))
 188                 %random-double-float))
 189
 190 ;;; 32-bit version
 191 #!+nil
 192 (defun %random-double-float (arg state)
 193   (declare (type (double-float (0d0)) arg)
 194            (type random-state state))
 195   (* (float (random-chunk state) 1d0) (/ 1d0 (expt 2 32))))
 196
 197 ;;; 53-bit version
 198 #!-x86
 199 (defun %random-double-float (arg state)
 200   (declare (type (double-float (0d0)) arg)
 201            (type random-state state))
 202   (* arg
 203      (- (sb!impl::make-double-float
 204          (dpb (ash (random-chunk state)
 205                    (- sb!vm:double-float-digits n-random-chunk-bits 32))
 206               sb!vm:double-float-significand-byte
 207               (sb!impl::double-float-high-bits 1d0))
 208          (random-chunk state))
 209         1d0)))
 210
 211 ;;; using a faster inline VOP
 212 #!+x86
 213 (defun %random-double-float (arg state)
 214   (declare (type (double-float (0d0)) arg)
 215            (type random-state state))
 216   (let ((state-vector (random-state-state state)))
 217     (* arg
 218        (- (sb!impl::make-double-float
 219            (dpb (ash (sb!vm::random-mt19937 state-vector)
 220                      (- sb!vm:double-float-digits n-random-chunk-bits
 221                         sb!vm:n-word-bits))
 222                 sb!vm:double-float-significand-byte
 223                 (sb!impl::double-float-high-bits 1d0))
 224            (sb!vm::random-mt19937 state-vector))
 225           1d0))))
 226
 227 \f
 228 ;;;; random integers
 229
 230 ;;; a bitmask M wide enough that (= (LOGAND INCLUSIVE-LIMIT M) INCLUSIVE-LIMIT)
 231 (declaim (inline %inclusive-random-integer-mask))
 232 (defun %inclusive-random-integer-mask (inclusive-limit)
 233   (1- (ash 1 (integer-length inclusive-limit))))
 234
 235 ;;; Transform a uniform sample from an at-least-as-large range into a
 236 ;;; random sample in [0,INCLUSIVE-LIMIT] range by throwing away
 237 ;;; too-big samples.
 238 ;;;
 239 ;;; Up through sbcl-1.0.4, the (RANDOM INTEGER) worked by taking (MOD
 240 ;;; RAW-MERSENNE-OUTPUT INTEGER). That introduces enough bias that it
 241 ;;; is unsuitable for some calculations (like the Metropolis Monte
 242 ;;; Carlo simulations that I, WHN, worked on in grad school): in the
 243 ;;; sbcl-1.0.4, the expectation value of a sample was (in the worst
 244 ;;; part of the range) more than 1.0001 times the ideal expectation
 245 ;;; value. Perhaps that was even ANSI-conformant, because ANSI says only
 246 ;;;   An approximately uniform choice distribution is used. If
 247 ;;;   LIMIT is an integer, each of the possible results occurs
 248 ;;;   with (approximate) probability 1/LIMIT.
 249 ;;; But I (WHN again) said "yuck," so these days we try to get a
 250 ;;; truly uniform distribution by translating RAW-MERSENNE-OUTPUT to
 251 ;;; our output using the second method recommended by
 252 ;;; <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/efaq.html>:
 253 ;;;   In the rare case that the rounding-off error becomes a problem,
 254 ;;;   obtain minimum integer n such that N<=2^n, generate integer
 255 ;;;   random numbers, take the most significant n bits, and discard
 256 ;;;   those more than or equal to N.
 257 ;;; (The first method recommended there differs from the sbcl-1.0.4
 258 ;;; method: the recommended method gives slightly different weights
 259 ;;; to different integers, but distributes the overweighted integers
 260 ;;; evenly across the range of outputs so that any skew would be
 261 ;;; of order (/ MOST-POSITIVE-MACHINE-WORD), rather than the (/ 10000)
 262 ;;; or so achieved by sbcl-1.0.4. That skew would probably be
 263 ;;; acceptable --- it seems to be exactly the kind of deviation that
 264 ;;; might have been anticipated in the ANSI CL standard. However, that
 265 ;;; recommended calculation involves multiplication and division of
 266 ;;; two-machine-word bignums, which is hard for us to do efficiently
 267 ;;; here. Without special low-level hacking to support such short
 268 ;;; bignums without consing, the accept-reject solution is not only
 269 ;;; more exact, but also likely more efficient.)
 270 (defmacro %inclusive-random-integer-accept-reject (raw-mersenne-output-expr
 271                                                    inclusive-limit)
 272   (with-unique-names (raw-mersenne-output inclusive-limit-once)
 273   `(loop
 274     with ,inclusive-limit-once = ,inclusive-limit
 275     for ,raw-mersenne-output = ,raw-mersenne-output-expr
 276     until (<= ,raw-mersenne-output ,inclusive-limit-once)
 277     finally (return ,raw-mersenne-output))))
 278
 279 ;;; an UNSIGNED-BYTE of N-CHUNKS chunks sampled from the Mersenne twister
 280 (declaim (inline %random-chunks))
 281 (defun %random-chunks (n-chunks state)
 282   ;; KLUDGE: This algorithm will cons O(N^2) words when constructing
 283   ;; an N-word result. To do better should be fundamentally easy, we
 284   ;; just need to do some low-level hack preallocating the bignum and
 285   ;; writing its words one by one.
 286   ;;
 287   ;; Note: The old sbcl-1.0.4 code did its analogous operation using a
 288   ;; mysterious RANDOM-INTEGER-OVERLAP parameter, "the amount that we
 289   ;; overlap chunks by when building a large integer to make up for
 290   ;; the loss of randomness in the low bits." No such thing is called
 291   ;; for on
 292   ;;   <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/efaq.html>:
 293   ;; they say there that it's OK just to concatenate words of twister
 294   ;; output with no overlap. Thus, crossing our fingers and hoping that
 295   ;; the previous RNG author didn't have some subtle reason to need
 296   ;; RANDOM-INTEGER-OVERLAP that we know not of, we just concatenate
 297   ;; chunks.
 298   (loop repeat n-chunks
 299         for result = 0 then (logior (ash result n-random-chunk-bits)
 300                                     (random-chunk state))
 301         finally (return result)))
 302
 303 ;;; an UNSIGNED-BYTE of N-BITS bits sampled from the Mersenne twister
 304 (declaim (inline %random-bits))
 305 (defun %random-bits (n-bits state)
 306   (multiple-value-bind (n-full-chunks n-extra-bits)
 307       (floor n-bits n-random-chunk-bits)
 308     (let ((full-chunks (%random-chunks n-full-chunks state)))
 309       (if (zerop n-extra-bits)
 310           full-chunks
 311           (logior full-chunks
 312                   (ash (logand (random-chunk state)
 313                                (1- (ash 1 n-extra-bits)))
 314                        (* n-full-chunks n-random-chunk-bits)))))))
 315
 316 ;;; the guts of (RANDOM (1+ INCLUSIVE-LIMIT))
 317 (defun %inclusive-random-integer (inclusive-limit state)
 318   (declare (optimize speed (space 1))) ; to ensure DEFTRANSFORM is enabled
 319   (%inclusive-random-integer inclusive-limit state))
 320
 321 ;;; the guts of RANDOM for the known-to-return-FIXNUM case
 322 ;;;
 323 ;;; We use INCLUSIVE-LIMIT instead of the (exclusive) LIMIT of RANDOM,
 324 ;;; because we want it to be a FIXNUM even for the possibly-common
 325 ;;; case of (RANDOM (1+ MOST-POSITIVE-FIXNUM)). (That case is what
 326 ;;; one might use for generating random hash values, e.g.)
 327 ;;; It also turns out to be just what's needed for INTEGER-LENGTH.
 328 (declaim (maybe-inline %inclusive-random-fixnum))
 329 (defun %inclusive-random-fixnum (inclusive-limit state)
 330   (declare (type (and fixnum unsigned-byte) inclusive-limit))
 331   (aver (<= inclusive-limit most-positive-random-chunk))
 332   (let (;; If this calculation needs to be optimized further, a good
 333         ;; start might be a DEFTRANSFORM which picks off the case of
 334         ;; constant LIMIT and precomputes the MASK at compile time.
 335         (mask (%inclusive-random-integer-mask inclusive-limit)))
 336     (%inclusive-random-integer-accept-reject (logand (random-chunk state) mask)
 337                                              inclusive-limit)))
 338 \f
 339 ;;;; outer, dynamically-typed interface
 340
 341 (defun random (arg &optional (state *random-state*))
 342   (declare (inline %random-single-float
 343                    %random-double-float
 344                    #!+long-float %random-long-float))
 345   (cond
 346     ((and (fixnump arg) (plusp arg))
 347      (locally
 348          ;; The choice to inline this very common case of
 349          ;; %INCLUSIVE-RANDOM-FIXNUM and not the less-common call
 350          ;; below was just a guess (by WHN 2007-03-27), not based on
 351          ;; benchmarking or anything.
 352          (declare (inline %inclusive-random-fixnum))
 353        (%inclusive-random-fixnum (1- arg) state)))
 354     ((and (typep arg 'single-float) (> arg 0.0f0))
 355      (%random-single-float arg state))
 356     ((and (typep arg 'double-float) (> arg 0.0d0))
 357      (%random-double-float arg state))
 358     #!+long-float
 359     ((and (typep arg 'long-float) (> arg 0.0l0))
 360      (%random-long-float arg state))
 361     ((and (integerp arg) (plusp arg))
 362      (if (= arg (1+ most-positive-fixnum))
 363          (%inclusive-random-fixnum most-positive-fixnum state)
 364          (%inclusive-random-integer (1- arg) state)))
 365     (t
 366      (error 'simple-type-error
 367             :expected-type '(or (integer 1) (float (0))) :datum arg
 368             :format-control "~@<Argument is neither a positive integer nor a ~
 369                              positive float: ~2I~_~S~:>"
 370             :format-arguments (list arg)))))