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 ;;;;   Vol. 8, No. 1, January pp.3-30 (1998) DOI:10.1145/272991.272995
   7 ;;;; http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
   8
   9 ;;;; This software is part of the SBCL system. See the README file for
  10 ;;;; more information.
  11 ;;;;
  12 ;;;; This software is derived from the CMU CL system, which was
  13 ;;;; written at Carnegie Mellon University and released into the
  14 ;;;; public domain. The software is in the public domain and is
  15 ;;;; provided with absolutely no warranty. See the COPYING and CREDITS
  16 ;;;; files for more information.
  17
  18 (in-package "SB!KERNEL")
  19 \f
  20 ;;;; Constants
  21 (defconstant mt19937-n 624)
  22 (defconstant mt19937-m 397)
  23 (defconstant mt19937-upper-mask #x80000000)
  24 (defconstant mt19937-lower-mask #x7FFFFFFF)
  25 (defconstant mt19937-a #x9908B0DF)
  26 (defconstant mt19937-b #x9D2C5680)
  27 (defconstant mt19937-c #xEFC60000)
  28
  29 ;;;; RANDOM-STATEs
  30
  31 ;;; The state is stored in a (simple-array (unsigned-byte 32) (627))
  32 ;;; wrapped in a random-state structure:
  33 ;;;
  34 ;;;  0-1:   Constant matrix A. [0, #x9908b0df]
  35 ;;;  2:     Index k.
  36 ;;;  3-626: State.
  37
  38 (deftype random-state-state () `(simple-array (unsigned-byte 32) (,(+ 3 mt19937-n))))
  39
  40 (def!method make-load-form ((random-state random-state) &optional environment)
  41   (make-load-form-saving-slots random-state :environment environment))
  42
  43 (def!method print-object ((state random-state) stream)
  44   (if (and *print-readably* (not *read-eval*))
  45       (error 'print-not-readable :object state)
  46       (format stream "#S(~S ~S #.~S)"
  47               'random-state
  48               ':state
  49               `(make-array (,(+ 3 mt19937-n))
  50                 :element-type
  51                 '(unsigned-byte 32)
  52                 :initial-contents
  53                 ',(coerce (random-state-state state) 'list)))))
  54
  55 ;;; Generate and initialize a new random-state array. Index is
  56 ;;; initialized to 1 and the states to 32bit integers excluding zero.
  57 ;;;
  58 ;;; Seed - A 32bit number.
  59 ;;;
  60 ;;; See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
  61 ;;; In the previous versions, MSBs of the seed affect only MSBs of the array.
  62 (defun init-random-state (&optional (seed 5489) state)
  63   (declare (type (unsigned-byte 32) seed))
  64   (let ((state (or state (make-array 627 :element-type '(unsigned-byte 32)))))
  65     (check-type state random-state-state)
  66     (setf (aref state 0) 0)
  67     (setf (aref state 1) mt19937-a)
  68     (setf (aref state 2) mt19937-n)
  69     (loop for i below mt19937-n
  70       for p from 3
  71       for s = seed then
  72       (logand #xFFFFFFFF
  73               (+ (* 1812433253
  74                     (logxor s (ash s -30)))
  75                  i))
  76       do (setf (aref state p) s))
  77     state))
  78
  79 (defvar *random-state*)
  80 (defun !random-cold-init ()
  81   (/show0 "entering !RANDOM-COLD-INIT")
  82   (setf *random-state* (%make-random-state))
  83   (/show0 "returning from !RANDOM-COLD-INIT"))
  84
  85 (defun make-random-state (&optional state)
  86   #!+sb-doc
  87   "Make a random state object. The optional STATE argument specifies a seed
  88   for deterministic pseudo-random number generation.
  89
  90   As per the Common Lisp standard,
  91   - If STATE is NIL or not supplied or is NIL, return a copy of the default
  92   *RANDOM-STATE*.
  93   - If STATE is a random state, return a copy of it.
  94   - If STATE is T, return a randomly initialized state (using operating-system
  95   provided randomness source where available, otherwise a poor substitute
  96   based on internal time and pid)
  97
  98   As an SBCL extension (starting with SBCL 1.0.33), we also support receiving
  99   as a seed an object of the following types:
 100   - (SIMPLE-ARRAY (UNSIGNED-BYTE 8) (*))
 101   - UNSIGNED-BYTE
 102   While we support arguments of any size and will mix the provided bits into
 103   the random state, it is probably overkill to provide more than 256 bits worth
 104   of actual information.
 105
 106   This particular SBCL version also accepts an argument of the following type:
 107   - (SIMPLE-ARRAY (UNSIGNED-BYTE 32) (*))
 108
 109   This particular SBCL version uses the popular MT19937 PRNG algorithm, and its
 110   internal state only effectively contains about 19937 bits of information.
 111   http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
 112 "
 113   (/show0 "entering MAKE-RANDOM-STATE")
 114   (etypecase state
 115     ;; Easy standard cases
 116     (null
 117      (/show0 "copying *RANDOM-STATE*")
 118      (%make-random-state :state (copy-seq (random-state-state *random-state*))))
 119     (random-state
 120      (/show0 "copying the provided RANDOM-STATE")
 121      (%make-random-state :state (copy-seq (random-state-state state))))
 122     ;; Standard case, less easy: try to randomly initialize a state.
 123     ((eql t)
 124      (/show0 "getting randomness from the operating system")
 125      (make-random-state
 126       (or
 127        ;; On unices, we try to read from /dev/urandom and pass the results
 128        ;; to our (simple-array (unsigned-byte 32) (*)) processor below.
 129        ;; More than 256 bits would provide a false sense of security.
 130        ;; If you need more bits than that, you probably also need
 131        ;; a better algorithm too.
 132        #!-win32
 133        (ignore-errors
 134          (with-open-file (r "/dev/urandom" :element-type '(unsigned-byte 32)
 135                             :direction :input :if-does-not-exist :error)
 136            (let ((a (make-array '(8) :element-type '(unsigned-byte 32))))
 137              (assert (= 8 (read-sequence a r)))
 138              a)))
 139        ;; When /dev/urandom is not available, we make do with time and pid
 140        ;; Thread ID and/or address of a CONS cell would be even better, but...
 141        (progn
 142          (/show0 "No /dev/urandom, using randomness from time and pid")
 143          (+ (get-internal-real-time)
 144             #!+unix (ash (sb!unix:unix-getpid) 32))))))
 145     ;; For convenience to users, we accept (simple-array (unsigned-byte 8) (*))
 146     ;; We just convert it to (simple-array (unsigned-byte 32) (*)) in a
 147     ;; completely straightforward way.
 148     ;; TODO: probably similarly accept other word sizes.
 149     ((simple-array (unsigned-byte 8) (*))
 150      (/show0 "getting random seed from byte vector (converting to 32-bit-word vector)")
 151      (let* ((l (length state))
 152             (m (ceiling l 4))
 153             (r (if (>= l 2496) 0 (mod l 4)))
 154             (y (make-array (list m) :element-type '(unsigned-byte 32))))
 155              (loop for i from 0 below (- m (if (zerop r) 0 1))
 156                for j = (* i 4) do
 157                (setf (aref y i)
 158                      (+ (aref state j)
 159                         (ash (aref state (+ j 1)) 8)
 160                         (ash (aref state (+ j 2)) 16)
 161                         (ash (aref state (+ j 3)) 24))))
 162              (unless (zerop r) ;; The last word may require special treatment.
 163                (let* ((p (1- m)) (q (* 4 p)))
 164                  (setf (aref y p)
 165                      (+ (aref state q)
 166                         (if (< 1 r) (ash (aref state (+ q 1)) 8) 0)
 167                         (if (= 3 r) (ash (aref state (+ q 2)) 16) 0)))))
 168              (make-random-state y)))
 169     ;; Also for convenience, we accept non-negative integers as seeds.
 170     ;; Small ones get passed to init-random-state, as before.
 171     ((unsigned-byte 32)
 172      (/show0 "getting random seed from 32-bit word")
 173      (%make-random-state :state (init-random-state state)))
 174     ;; Larger ones ones get trivially chopped into an array of (unsigned-byte 32)
 175     ((unsigned-byte)
 176      (/show0 "getting random seed from bignum (converting to 32-bit-word vector)")
 177      (loop with l = (ceiling (integer-length state) 32)
 178        with s = (make-array (list l) :element-type '(unsigned-byte 32))
 179        for i below l
 180        for p from 0 by 32
 181        do (setf (aref s i) (ldb (byte 32 p) state))
 182        finally (return (make-random-state s))))
 183     ;; Last but not least, when provided an array of 32-bit words, we truncate
 184     ;; it to 19968 bits and mix these into an initial state. We reuse the same
 185     ;; method as the authors of the original algorithm. See
 186     ;; http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c
 187     ;; NB: their mt[i] is our (aref s (+ 3 i))
 188     ((simple-array (unsigned-byte 32) (*))
 189      (/show0 "getting random seed from 32-bit-word vector")
 190      (let ((s (init-random-state 19650218))
 191            (i 1) (j 0) (l (length state)))
 192        (loop for k downfrom (max mt19937-n l) above 0 do
 193          (setf (aref s (+ i 3))
 194                (logand #xFFFFFFFF
 195                        (+ (logxor (aref s (+ i 3))
 196                                   (* 1664525
 197                                      (logxor (aref s (+ i 2))
 198                                              (ash (aref s (+ i 2)) -30))))
 199                           (aref state j) j))) ;; non-linear
 200          (incf i) (when (>= i mt19937-n) (setf (aref s 3) (aref s (+ 2 mt19937-n)) i 1))
 201          (incf j) (when (>= j l) (setf j 0)))
 202        (loop for k downfrom (1- mt19937-n) above 0 do
 203          (setf (aref s (+ i 3))
 204                (logand #xFFFFFFFF
 205                        (- (logxor (aref s (+ i 3))
 206                                   (* 1566083941
 207                                      (logxor (aref s (+ i 2))
 208                                              (ash (aref s (+ i 2)) -30))))
 209                           i))) ;; non-linear
 210          (incf i) (when (>= i mt19937-n) (setf (aref s 3) (aref s (+ 2 mt19937-n)) i 1)))
 211        (setf (aref s 3) #x80000000) ;; MSB is 1; assuring non-zero initial array
 212        (%make-random-state :state s)))))
 213 \f
 214 ;;;; random entries
 215
 216 ;;; This function generates a 32bit integer between 0 and #xffffffff
 217 ;;; inclusive.
 218 #!-sb-fluid (declaim (inline random-chunk))
 219 ;;; portable implementation
 220 #!-x86
 221 (defun random-mt19937-update (state)
 222   (declare (type random-state-state state)
 223            (optimize (speed 3) (safety 0)))
 224   (let ((y 0))
 225     (declare (type (unsigned-byte 32) y))
 226     (do ((kk 3 (1+ kk)))
 227         ((>= kk (+ 3 (- mt19937-n mt19937-m))))
 228       (declare (type (mod 628) kk))
 229       (setf y (logior (logand (aref state kk) mt19937-upper-mask)
 230                       (logand (aref state (1+ kk)) mt19937-lower-mask)))
 231       (setf (aref state kk) (logxor (aref state (+ kk mt19937-m))
 232                                     (ash y -1) (aref state (logand y 1)))))
 233     (do ((kk (+ (- mt19937-n mt19937-m) 3) (1+ kk)))
 234         ((>= kk (+ (1- mt19937-n) 3)))
 235       (declare (type (mod 628) kk))
 236       (setf y (logior (logand (aref state kk) mt19937-upper-mask)
 237                       (logand (aref state (1+ kk)) mt19937-lower-mask)))
 238       (setf (aref state kk) (logxor (aref state (+ kk (- mt19937-m mt19937-n)))
 239                                     (ash y -1) (aref state (logand y 1)))))
 240     (setf y (logior (logand (aref state (+ 3 (1- mt19937-n)))
 241                             mt19937-upper-mask)
 242                     (logand (aref state 3) mt19937-lower-mask)))
 243     (setf (aref state (+ 3 (1- mt19937-n)))
 244           (logxor (aref state (+ 3 (1- mt19937-m)))
 245                   (ash y -1) (aref state (logand y 1)))))
 246   (values))
 247 #!-x86
 248 (defun random-chunk (state)
 249   (declare (type random-state state))
 250   (let* ((state (random-state-state state))
 251          (k (aref state 2)))
 252     (declare (type (mod 628) k))
 253     (when (= k mt19937-n)
 254       (random-mt19937-update state)
 255       (setf k 0))
 256     (setf (aref state 2) (1+ k))
 257     (let ((y (aref state (+ 3 k))))
 258       (declare (type (unsigned-byte 32) y))
 259       (setf y (logxor y (ash y -11)))
 260       (setf y (logxor y (ash (logand y (ash mt19937-b -7)) 7)))
 261       (setf y (logxor y (ash (logand y (ash mt19937-c -15)) 15)))
 262       (setf y (logxor y (ash y -18)))
 263       y)))
 264
 265 ;;; Using inline VOP support, only available on the x86 so far.
 266 ;;;
 267 ;;; FIXME: It would be nice to have some benchmark numbers on this.
 268 ;;; My inclination is to get rid of the nonportable implementation
 269 ;;; unless the performance difference is just enormous.
 270 #!+x86
 271 (defun random-chunk (state)
 272   (declare (type random-state state))
 273   (sb!vm::random-mt19937 (random-state-state state)))
 274
 275 #!-sb-fluid (declaim (inline big-random-chunk))
 276 (defun big-random-chunk (state)
 277   (declare (type random-state state))
 278   (logand (1- (expt 2 64))
 279           (logior (ash (random-chunk state) 32)
 280                   (random-chunk state))))
 281 \f
 282 ;;; Handle the single or double float case of RANDOM. We generate a
 283 ;;; float between 0.0 and 1.0 by clobbering the significand of 1.0
 284 ;;; with random bits, then subtracting 1.0. This hides the fact that
 285 ;;; we have a hidden bit.
 286 #!-sb-fluid (declaim (inline %random-single-float %random-double-float))
 287 (declaim (ftype (function ((single-float (0f0)) random-state)
 288                           (single-float 0f0))
 289                 %random-single-float))
 290 (defun %random-single-float (arg state)
 291   (declare (type (single-float (0f0)) arg)
 292            (type random-state state))
 293   (* arg
 294      (- (make-single-float
 295          (dpb (ash (random-chunk state)
 296                    (- sb!vm:single-float-digits random-chunk-length))
 297               sb!vm:single-float-significand-byte
 298               (single-float-bits 1.0)))
 299         1.0)))
 300 (declaim (ftype (function ((double-float (0d0)) random-state)
 301                           (double-float 0d0))
 302                 %random-double-float))
 303
 304 ;;; 32-bit version
 305 #!+nil
 306 (defun %random-double-float (arg state)
 307   (declare (type (double-float (0d0)) arg)
 308            (type random-state state))
 309   (* (float (random-chunk state) 1d0) (/ 1d0 (expt 2 32))))
 310
 311 ;;; 53-bit version
 312 #!-x86
 313 (defun %random-double-float (arg state)
 314   (declare (type (double-float (0d0)) arg)
 315            (type random-state state))
 316   (* arg
 317      (- (sb!impl::make-double-float
 318          (dpb (ash (random-chunk state)
 319                    (- sb!vm:double-float-digits random-chunk-length 32))
 320               sb!vm:double-float-significand-byte
 321               (sb!impl::double-float-high-bits 1d0))
 322          (random-chunk state))
 323         1d0)))
 324
 325 ;;; using a faster inline VOP
 326 #!+x86
 327 (defun %random-double-float (arg state)
 328   (declare (type (double-float (0d0)) arg)
 329            (type random-state state))
 330   (let ((state-vector (random-state-state state)))
 331     (* arg
 332        (- (sb!impl::make-double-float
 333            (dpb (ash (sb!vm::random-mt19937 state-vector)
 334                      (- sb!vm:double-float-digits random-chunk-length
 335                         sb!vm:n-word-bits))
 336                 sb!vm:double-float-significand-byte
 337                 (sb!impl::double-float-high-bits 1d0))
 338            (sb!vm::random-mt19937 state-vector))
 339           1d0))))
 340
 341 \f
 342 ;;;; random integers
 343
 344 (defun %random-integer (arg state)
 345   (declare (type (integer 1) arg) (type random-state state))
 346   (let ((shift (- random-chunk-length random-integer-overlap)))
 347     (do ((bits (random-chunk state)
 348                (logxor (ash bits shift) (random-chunk state)))
 349          (count (+ (integer-length arg)
 350                    (- random-integer-extra-bits shift))
 351                 (- count shift)))
 352         ((minusp count)
 353          (rem bits arg))
 354       (declare (fixnum count)))))
 355
 356 (defun random (arg &optional (state *random-state*))
 357   (declare (inline %random-single-float %random-double-float
 358                    #!+long-float %random-long-float))
 359   (cond
 360     ((and (fixnump arg) (<= arg random-fixnum-max) (> arg 0))
 361      (rem (random-chunk state) arg))
 362     ((and (typep arg 'single-float) (> arg 0.0f0))
 363      (%random-single-float arg state))
 364     ((and (typep arg 'double-float) (> arg 0.0d0))
 365      (%random-double-float arg state))
 366     #!+long-float
 367     ((and (typep arg 'long-float) (> arg 0.0l0))
 368      (%random-long-float arg state))
 369     ((and (integerp arg) (> arg 0))
 370      (%random-integer arg state))
 371     (t
 372      (error 'simple-type-error
 373             :expected-type '(or (integer 1) (float (0))) :datum arg
 374             :format-control "~@<Argument is neither a positive integer nor a ~
 375                              positive float: ~2I~_~S~:>"
 376             :format-arguments (list arg)))))