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 ;;; The state is stored in a (simple-array (unsigned-byte 32) (627))
  25 ;;; wrapped in a random-state structure:
  26 ;;;
  27 ;;;  0-1:   Constant matrix A. [0, #x9908b0df]
  28 ;;;  2:     Index k.
  29 ;;;  3-626: State.
  30
  31 ;;; Generate and initialize a new random-state array. Index is
  32 ;;; initialized to 1 and the states to 32bit integers excluding zero.
  33 ;;;
  34 ;;; Seed - A 32bit number, not zero.
  35 ;;;
  36 ;;; Apparently uses the generator Line 25 of Table 1 in
  37 ;;; [KNUTH 1981, The Art of Computer Programming, Vol. 2 (2nd Ed.), pp102]
  38 (defun init-random-state (&optional (seed 4357) state)
  39   (declare (type (integer 1 #xffffffff) seed))
  40   (let ((state (or state (make-array 627 :element-type '(unsigned-byte 32)))))
  41     (declare (type (simple-array (unsigned-byte 32) (627)) state))
  42     (setf (aref state 1) #x9908b0df)
  43     (setf (aref state 2) 1)
  44     (setf (aref state 3) seed)
  45     (do ((k 1 (1+ k)))
  46         ((>= k 624))
  47       (declare (type (mod 625) k))
  48       (setf (aref state (+ 3 k))
  49             (logand (* 69069 (aref state (+ 3 (1- k)))) #xffffffff)))
  50     state))
  51
  52 (defvar *random-state*)
  53 (defun !random-cold-init ()
  54   (/show0 "entering !RANDOM-COLD-INIT")
  55   (setf *random-state* (%make-random-state))
  56   (/show0 "returning from !RANDOM-COLD-INIT"))
  57
  58 (defun make-random-state (&optional state)
  59   #!+sb-doc
  60   "Make a random state object. If STATE is not supplied, return a copy
  61   of the default random state. If STATE is a random state, then return a
  62   copy of it. If STATE is T then return a random state generated from
  63   the universal time."
  64   (/show0 "entering !RANDOM-COLD-INIT")
  65   (flet ((copy-random-state (state)
  66            (/show0 "entering COPY-RANDOM-STATE")
  67            (let ((state (random-state-state state))
  68                  (new-state
  69                   (make-array 627 :element-type '(unsigned-byte 32))))
  70              (/show0 "made NEW-STATE, about to DOTIMES")
  71              (dotimes (i 627)
  72                (setf (aref new-state i) (aref state i)))
  73              (/show0 "falling through to %MAKE-RANDOM-STATE")
  74              (%make-random-state :state new-state))))
  75     (/show0 "at head of ETYPECASE in MAKE-RANDOM-STATE")
  76     (etypecase state
  77       (null
  78        (/show0 "NULL case")
  79        (copy-random-state *random-state*))
  80       (random-state
  81        (/show0 "RANDOM-STATE-P clause")
  82        (copy-random-state state))
  83       ((member t)
  84        (/show0 "T clause")
  85        (%make-random-state :state (init-random-state
  86                                    (logand (get-universal-time)
  87                                            #xffffffff)))))))
  88 \f
  89 ;;;; random entries
  90
  91 ;;; This function generates a 32bit integer between 0 and #xffffffff
  92 ;;; inclusive.
  93 #!-sb-fluid (declaim (inline random-chunk))
  94 ;;; portable implementation
  95 (defconstant mt19937-n 624)
  96 (defconstant mt19937-m 397)
  97 (defconstant mt19937-upper-mask #x80000000)
  98 (defconstant mt19937-lower-mask #x7fffffff)
  99 (defconstant mt19937-b #x9D2C5680)
 100 (defconstant mt19937-c #xEFC60000)
 101 #!-x86
 102 (defun random-mt19937-update (state)
 103   (declare (type (simple-array (unsigned-byte 32) (627)) state)
 104            (optimize (speed 3) (safety 0)))
 105   (let ((y 0))
 106     (declare (type (unsigned-byte 32) y))
 107     (do ((kk 3 (1+ kk)))
 108         ((>= kk (+ 3 (- mt19937-n mt19937-m))))
 109       (declare (type (mod 628) kk))
 110       (setf y (logior (logand (aref state kk) mt19937-upper-mask)
 111                       (logand (aref state (1+ kk)) mt19937-lower-mask)))
 112       (setf (aref state kk) (logxor (aref state (+ kk mt19937-m))
 113                                     (ash y -1) (aref state (logand y 1)))))
 114     (do ((kk (+ (- mt19937-n mt19937-m) 3) (1+ kk)))
 115         ((>= kk (+ (1- mt19937-n) 3)))
 116       (declare (type (mod 628) kk))
 117       (setf y (logior (logand (aref state kk) mt19937-upper-mask)
 118                       (logand (aref state (1+ kk)) mt19937-lower-mask)))
 119       (setf (aref state kk) (logxor (aref state (+ kk (- mt19937-m mt19937-n)))
 120                                     (ash y -1) (aref state (logand y 1)))))
 121     (setf y (logior (logand (aref state (+ 3 (1- mt19937-n)))
 122                             mt19937-upper-mask)
 123                     (logand (aref state 3) mt19937-lower-mask)))
 124     (setf (aref state (+ 3 (1- mt19937-n)))
 125           (logxor (aref state (+ 3 (1- mt19937-m)))
 126                   (ash y -1) (aref state (logand y 1)))))
 127   (values))
 128 #!-x86
 129 (defun random-chunk (state)
 130   (declare (type random-state state))
 131   (let* ((state (random-state-state state))
 132          (k (aref state 2)))
 133     (declare (type (mod 628) k))
 134     (when (= k mt19937-n)
 135       (random-mt19937-update state)
 136       (setf k 0))
 137     (setf (aref state 2) (1+ k))
 138     (let ((y (aref state (+ 3 k))))
 139       (declare (type (unsigned-byte 32) y))
 140       (setf y (logxor y (ash y -11)))
 141       (setf y (logxor y (ash (logand y (ash mt19937-b -7)) 7)))
 142       (setf y (logxor y (ash (logand y (ash mt19937-c -15)) 15)))
 143       (setf y (logxor y (ash y -18)))
 144       y)))
 145
 146 ;;; Using inline VOP support, only available on the x86 so far.
 147 ;;;
 148 ;;; FIXME: It would be nice to have some benchmark numbers on this.
 149 ;;; My inclination is to get rid of the nonportable implementation
 150 ;;; unless the performance difference is just enormous.
 151 #!+x86
 152 (defun random-chunk (state)
 153   (declare (type random-state state))
 154   (sb!vm::random-mt19937 (random-state-state state)))
 155 \f
 156 ;;; Handle the single or double float case of RANDOM. We generate a
 157 ;;; float between 0.0 and 1.0 by clobbering the significand of 1.0
 158 ;;; with random bits, then subtracting 1.0. This hides the fact that
 159 ;;; we have a hidden bit.
 160 #!-sb-fluid (declaim (inline %random-single-float %random-double-float))
 161 (declaim (ftype (function ((single-float (0f0)) random-state)
 162                           (single-float 0f0))
 163                 %random-single-float))
 164 (defun %random-single-float (arg state)
 165   (declare (type (single-float (0f0)) arg)
 166            (type random-state state))
 167   (* arg
 168      (- (make-single-float
 169          (dpb (ash (random-chunk state)
 170                    (- sb!vm:single-float-digits random-chunk-length))
 171               sb!vm:single-float-significand-byte
 172               (single-float-bits 1.0)))
 173         1.0)))
 174 (declaim (ftype (function ((double-float (0d0)) random-state)
 175                           (double-float 0d0))
 176                 %random-double-float))
 177
 178 ;;; 32-bit version
 179 #!+nil
 180 (defun %random-double-float (arg state)
 181   (declare (type (double-float (0d0)) arg)
 182            (type random-state state))
 183   (* (float (random-chunk state) 1d0) (/ 1d0 (expt 2 32))))
 184
 185 ;;; 53-bit version
 186 #!-x86
 187 (defun %random-double-float (arg state)
 188   (declare (type (double-float (0d0)) arg)
 189            (type random-state state))
 190   (* arg
 191      (- (sb!impl::make-double-float
 192          (dpb (ash (random-chunk state)
 193                    (- sb!vm:double-float-digits random-chunk-length
 194                       sb!vm:n-word-bits))
 195               sb!vm:double-float-significand-byte
 196               (sb!impl::double-float-high-bits 1d0))
 197          (random-chunk state))
 198         1d0)))
 199
 200 ;;; using a faster inline VOP
 201 #!+x86
 202 (defun %random-double-float (arg state)
 203   (declare (type (double-float (0d0)) arg)
 204            (type random-state state))
 205   (let ((state-vector (random-state-state state)))
 206     (* arg
 207        (- (sb!impl::make-double-float
 208            (dpb (ash (sb!vm::random-mt19937 state-vector)
 209                      (- sb!vm:double-float-digits random-chunk-length
 210                         sb!vm:n-word-bits))
 211                 sb!vm:double-float-significand-byte
 212                 (sb!impl::double-float-high-bits 1d0))
 213            (sb!vm::random-mt19937 state-vector))
 214           1d0))))
 215
 216 #!+long-float
 217 (declaim #!-sb-fluid (inline %random-long-float))
 218 #!+long-float
 219 (declaim (ftype (function ((long-float (0l0)) random-state) (long-float 0l0))
 220                 %random-long-float))
 221
 222 ;;; using a faster inline VOP
 223 #!+(and long-float x86)
 224 (defun %random-long-float (arg state)
 225   (declare (type (long-float (0l0)) arg)
 226            (type random-state state))
 227   (let ((state-vector (random-state-state state)))
 228     (* arg
 229        (- (sb!impl::make-long-float
 230            (sb!impl::long-float-exp-bits 1l0)
 231            (logior (sb!vm::random-mt19937 state-vector)
 232                    sb!vm:long-float-hidden-bit)
 233            (sb!vm::random-mt19937 state-vector))
 234           1l0))))
 235
 236 #!+(and long-float sparc)
 237 (defun %random-long-float (arg state)
 238   (declare (type (long-float (0l0)) arg)
 239            (type random-state state))
 240   (* arg
 241      (- (sb!impl::make-long-float
 242          (sb!impl::long-float-exp-bits 1l0)     ; X needs more work
 243          (random-chunk state) (random-chunk state) (random-chunk state))
 244         1l0)))
 245 \f
 246 ;;;; random integers
 247
 248 (defun %random-integer (arg state)
 249   (declare (type (integer 1) arg) (type random-state state))
 250   (let ((shift (- random-chunk-length random-integer-overlap)))
 251     (do ((bits (random-chunk state)
 252                (logxor (ash bits shift) (random-chunk state)))
 253          (count (+ (integer-length arg)
 254                    (- random-integer-extra-bits shift))
 255                 (- count shift)))
 256         ((minusp count)
 257          (rem bits arg))
 258       (declare (fixnum count)))))
 259
 260 (defun random (arg &optional (state *random-state*))
 261   (declare (inline %random-single-float %random-double-float
 262                    #!+long-float %random-long-float))
 263   (cond
 264     ((and (fixnump arg) (<= arg random-fixnum-max) (> arg 0))
 265      (rem (random-chunk state) arg))
 266     ((and (typep arg 'single-float) (> arg 0.0f0))
 267      (%random-single-float arg state))
 268     ((and (typep arg 'double-float) (> arg 0.0d0))
 269      (%random-double-float arg state))
 270     #!+long-float
 271     ((and (typep arg 'long-float) (> arg 0.0l0))
 272      (%random-long-float arg state))
 273     ((and (integerp arg) (> arg 0))
 274      (%random-integer arg state))
 275     (t
 276      (error 'simple-type-error
 277             :expected-type '(or (integer 1) (float (0))) :datum arg
 278             :format-control "~@<Argument is neither a positive integer nor a ~
 279                              positive float: ~2I~_~S~:>"
 280             :format-arguments (list arg)))))