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
 168 #!-sb-fluid (declaim (inline big-random-chunk))
 169 (defun big-random-chunk (state)
 170   (declare (type random-state state))
 171   (logand (1- (expt 2 64))
 172           (logior (ash (random-chunk state) 32)
 173                   (random-chunk state))))
 174 \f
 175 ;;; Handle the single or double float case of RANDOM. We generate a
 176 ;;; float between 0.0 and 1.0 by clobbering the significand of 1.0
 177 ;;; with random bits, then subtracting 1.0. This hides the fact that
 178 ;;; we have a hidden bit.
 179 #!-sb-fluid (declaim (inline %random-single-float %random-double-float))
 180 (declaim (ftype (function ((single-float (0f0)) random-state)
 181                           (single-float 0f0))
 182                 %random-single-float))
 183 (defun %random-single-float (arg state)
 184   (declare (type (single-float (0f0)) arg)
 185            (type random-state state))
 186   (* arg
 187      (- (make-single-float
 188          (dpb (ash (random-chunk state)
 189                    (- sb!vm:single-float-digits random-chunk-length))
 190               sb!vm:single-float-significand-byte
 191               (single-float-bits 1.0)))
 192         1.0)))
 193 (declaim (ftype (function ((double-float (0d0)) random-state)
 194                           (double-float 0d0))
 195                 %random-double-float))
 196
 197 ;;; 32-bit version
 198 #!+nil
 199 (defun %random-double-float (arg state)
 200   (declare (type (double-float (0d0)) arg)
 201            (type random-state state))
 202   (* (float (random-chunk state) 1d0) (/ 1d0 (expt 2 32))))
 203
 204 ;;; 53-bit version
 205 #!-x86
 206 (defun %random-double-float (arg state)
 207   (declare (type (double-float (0d0)) arg)
 208            (type random-state state))
 209   (* arg
 210      (- (sb!impl::make-double-float
 211          (dpb (ash (random-chunk state)
 212                    (- sb!vm:double-float-digits random-chunk-length 32))
 213               sb!vm:double-float-significand-byte
 214               (sb!impl::double-float-high-bits 1d0))
 215          (random-chunk state))
 216         1d0)))
 217
 218 ;;; using a faster inline VOP
 219 #!+x86
 220 (defun %random-double-float (arg state)
 221   (declare (type (double-float (0d0)) arg)
 222            (type random-state state))
 223   (let ((state-vector (random-state-state state)))
 224     (* arg
 225        (- (sb!impl::make-double-float
 226            (dpb (ash (sb!vm::random-mt19937 state-vector)
 227                      (- sb!vm:double-float-digits random-chunk-length
 228                         sb!vm:n-word-bits))
 229                 sb!vm:double-float-significand-byte
 230                 (sb!impl::double-float-high-bits 1d0))
 231            (sb!vm::random-mt19937 state-vector))
 232           1d0))))
 233
 234 \f
 235 ;;;; random integers
 236
 237 (defun %random-integer (arg state)
 238   (declare (type (integer 1) arg) (type random-state state))
 239   (let ((shift (- random-chunk-length random-integer-overlap)))
 240     (do ((bits (random-chunk state)
 241                (logxor (ash bits shift) (random-chunk state)))
 242          (count (+ (integer-length arg)
 243                    (- random-integer-extra-bits shift))
 244                 (- count shift)))
 245         ((minusp count)
 246          (rem bits arg))
 247       (declare (fixnum count)))))
 248
 249 (defun random (arg &optional (state *random-state*))
 250   (declare (inline %random-single-float %random-double-float
 251                    #!+long-float %random-long-float))
 252   (cond
 253     ((and (fixnump arg) (<= arg random-fixnum-max) (> arg 0))
 254      (rem (random-chunk state) arg))
 255     ((and (typep arg 'single-float) (> arg 0.0f0))
 256      (%random-single-float arg state))
 257     ((and (typep arg 'double-float) (> arg 0.0d0))
 258      (%random-double-float arg state))
 259     #!+long-float
 260     ((and (typep arg 'long-float) (> arg 0.0l0))
 261      (%random-long-float arg state))
 262     ((and (integerp arg) (> arg 0))
 263      (%random-integer arg state))
 264     (t
 265      (error 'simple-type-error
 266             :expected-type '(or (integer 1) (float (0))) :datum arg
 267             :format-control "~@<Argument is neither a positive integer nor a ~
 268                              positive float: ~2I~_~S~:>"
 269             :format-arguments (list arg)))))