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