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