;;;; Makoto Matsumoto and T. Nishimura, "Mersenne twister: A
;;;; 623-dimensionally equidistributed uniform pseudorandom number
;;;; generator.", ACM Transactions on Modeling and Computer Simulation,
-;;;; 1997, to appear.
+;;;; Vol. 8, No. 1, January pp.3-30 (1998) DOI:10.1145/272991.272995
+;;;; http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
;;;; This software is part of the SBCL system. See the README file for
;;;; more information.
(in-package "SB!KERNEL")
\f
+;;;; Constants
+(defconstant mt19937-n 624)
+(defconstant mt19937-m 397)
+(defconstant mt19937-upper-mask #x80000000)
+(defconstant mt19937-lower-mask #x7FFFFFFF)
+(defconstant mt19937-a #x9908B0DF)
+(defconstant mt19937-b #x9D2C5680)
+(defconstant mt19937-c #xEFC60000)
+
;;;; RANDOM-STATEs
+;;; The state is stored in a (simple-array (unsigned-byte 32) (627))
+;;; wrapped in a random-state structure:
+;;;
+;;; 0-1: Constant matrix A. [0, #x9908b0df]
+;;; 2: Index k.
+;;; 3-626: State.
+
+(deftype random-state-state () `(simple-array (unsigned-byte 32) (,(+ 3 mt19937-n))))
+
(def!method make-load-form ((random-state random-state) &optional environment)
(make-load-form-saving-slots random-state :environment environment))
(def!method print-object ((state random-state) stream)
(if (and *print-readably* (not *read-eval*))
- (error 'print-not-readable :object state)
+ (print-not-readable-error state stream)
(format stream "#S(~S ~S #.~S)"
'random-state
':state
- `(make-array 627
+ `(make-array ,(+ 3 mt19937-n)
:element-type
'(unsigned-byte 32)
:initial-contents
',(coerce (random-state-state state) 'list)))))
-;;; The state is stored in a (simple-array (unsigned-byte 32) (627))
-;;; wrapped in a random-state structure:
-;;;
-;;; 0-1: Constant matrix A. [0, #x9908b0df]
-;;; 2: Index k.
-;;; 3-626: State.
-
;;; Generate and initialize a new random-state array. Index is
;;; initialized to 1 and the states to 32bit integers excluding zero.
;;;
-;;; Seed - A 32bit number, not zero.
+;;; Seed - A 32bit number.
;;;
-;;; Apparently uses the generator Line 25 of Table 1 in
-;;; [KNUTH 1981, The Art of Computer Programming, Vol. 2 (2nd Ed.), pp102]
-(defun init-random-state (&optional (seed 4357) state)
- (declare (type (integer 1 #xffffffff) seed))
+;;; See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
+;;; In the previous versions, MSBs of the seed affect only MSBs of the array.
+(defun init-random-state (&optional (seed 5489) state)
+ (declare (type (unsigned-byte 32) seed))
(let ((state (or state (make-array 627 :element-type '(unsigned-byte 32)))))
- (declare (type (simple-array (unsigned-byte 32) (627)) state))
- (setf (aref state 1) #x9908b0df)
- (setf (aref state 2) 1)
- (setf (aref state 3) seed)
- (do ((k 1 (1+ k)))
- ((>= k 624))
- (declare (type (mod 625) k))
- (setf (aref state (+ 3 k))
- (logand (* 69069 (aref state (+ 3 (1- k)))) #xffffffff)))
+ (check-type state random-state-state)
+ (setf (aref state 0) 0)
+ (setf (aref state 1) mt19937-a)
+ (setf (aref state 2) mt19937-n)
+ (loop for i below mt19937-n
+ for p from 3
+ for s = seed then
+ (logand #xFFFFFFFF
+ (+ (* 1812433253
+ (logxor s (ash s -30)))
+ i))
+ do (setf (aref state p) s))
state))
(defvar *random-state*)
(defun make-random-state (&optional state)
#!+sb-doc
- "Make a RANDOM-STATE object. If STATE is not supplied, return a copy
- of the default random state. If STATE is a random state, then return a
- copy of it. If STATE is T then return a random state generated from
- the universal time."
+ "Make a random state object. The optional STATE argument specifies a seed
+for deterministic pseudo-random number generation.
+
+As per the Common Lisp standard,
+- If STATE is NIL or not supplied, return a copy of the default
+ *RANDOM-STATE*.
+- If STATE is a random state, return a copy of it.
+- If STATE is T, return a randomly initialized state (using operating-system
+ provided randomness where available, otherwise a poor substitute based on
+ internal time and PID).
+
+See SB-EXT:SEED-RANDOM-STATE for a SBCL extension to this functionality."
(/show0 "entering MAKE-RANDOM-STATE")
- (flet ((copy-random-state (state)
- (/show0 "entering COPY-RANDOM-STATE")
- (let ((state (random-state-state state))
- (new-state
- (make-array 627 :element-type '(unsigned-byte 32))))
- (/show0 "made NEW-STATE, about to DOTIMES")
- (dotimes (i 627)
- (setf (aref new-state i) (aref state i)))
- (/show0 "falling through to %MAKE-RANDOM-STATE")
- (%make-random-state :state new-state))))
- (/show0 "at head of ETYPECASE in MAKE-RANDOM-STATE")
- (etypecase state
- (null
- (/show0 "NULL case")
- (copy-random-state *random-state*))
- (random-state
- (/show0 "RANDOM-STATE-P clause")
- (copy-random-state state))
- ((member t)
- (/show0 "T clause")
- (%make-random-state :state (init-random-state
- (logand (get-universal-time)
- #xffffffff)))))))
+ (check-type state (or boolean random-state))
+ (seed-random-state state))
+
+(defun fallback-random-seed ()
+ ;; When /dev/urandom is not available, we make do with time and pid
+ ;; Thread ID and/or address of a CONS cell would be even better, but...
+ (/show0 "No /dev/urandom, using randomness from time and pid")
+ (+ (get-internal-real-time)
+ (ash (sb!unix:unix-getpid) 32)))
+
+#!-win32
+(defun os-random-seed ()
+ (or
+ ;; On unices, we try to read from /dev/urandom and pass the results
+ ;; to our (simple-array (unsigned-byte 32) (*)) processor below.
+ ;; More than 256 bits would provide a false sense of security.
+ ;; If you need more bits than that, you probably also need
+ ;; a better algorithm too.
+ (ignore-errors
+ (with-open-file (r "/dev/urandom" :element-type '(unsigned-byte 32)
+ :direction :input :if-does-not-exist :error)
+ (let ((a (make-array '(8) :element-type '(unsigned-byte 32))))
+ (assert (= 8 (read-sequence a r)))
+ a)))
+ (fallback-random-seed)))
+
+#!+win32
+(defun os-random-seed ()
+ (/show0 "Getting randomness from CryptGenRandom")
+ (or (sb!win32:crypt-gen-random 32)
+ (fallback-random-seed)))
+
+(defun seed-random-state (&optional state)
+ #!+sb-doc
+ "Make a random state object. The optional STATE argument specifies a seed
+for deterministic pseudo-random number generation.
+
+As per the Common Lisp standard for MAKE-RANDOM-STATE,
+- If STATE is NIL or not supplied, return a copy of the default
+ *RANDOM-STATE*.
+- If STATE is a random state, return a copy of it.
+- If STATE is T, return a randomly initialized state (using operating-system
+ provided randomness where available, otherwise a poor substitute based on
+ internal time and pid).
+
+As a supported SBCL extension, we also support receiving as a seed an object
+of the following types:
+- (SIMPLE-ARRAY (UNSIGNED-BYTE 8) (*))
+- UNSIGNED-BYTE
+While we support arguments of any size and will mix the provided bits into
+the random state, it is probably overkill to provide more than 256 bits worth
+of actual information.
+
+This particular SBCL version also accepts an argument of the following type:
+(SIMPLE-ARRAY (UNSIGNED-BYTE 32) (*))
+
+This particular SBCL version uses the popular MT19937 PRNG algorithm, and its
+internal state only effectively contains about 19937 bits of information.
+http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
+"
+ (etypecase state
+ ;; Easy standard cases
+ (null
+ (/show0 "copying *RANDOM-STATE*")
+ (%make-random-state :state (copy-seq (random-state-state *random-state*))))
+ (random-state
+ (/show0 "copying the provided RANDOM-STATE")
+ (%make-random-state :state (copy-seq (random-state-state state))))
+ ;; Standard case, less easy: try to randomly initialize a state.
+ ((eql t)
+ (/show0 "getting randomness from the operating system")
+ (seed-random-state (os-random-seed)))
+ ;; For convenience to users, we accept (simple-array (unsigned-byte 8) (*))
+ ;; We just convert it to (simple-array (unsigned-byte 32) (*)) in a
+ ;; completely straightforward way.
+ ;; TODO: probably similarly accept other word sizes.
+ ((simple-array (unsigned-byte 8) (*))
+ (/show0 "getting random seed from byte vector (converting to 32-bit-word vector)")
+ (let* ((l (length state))
+ (m (ceiling l 4))
+ (r (if (>= l 2496) 0 (mod l 4)))
+ (y (make-array (list m) :element-type '(unsigned-byte 32))))
+ (loop for i from 0 below (- m (if (zerop r) 0 1))
+ for j = (* i 4) do
+ (setf (aref y i)
+ (+ (aref state j)
+ (ash (aref state (+ j 1)) 8)
+ (ash (aref state (+ j 2)) 16)
+ (ash (aref state (+ j 3)) 24))))
+ (unless (zerop r) ;; The last word may require special treatment.
+ (let* ((p (1- m)) (q (* 4 p)))
+ (setf (aref y p)
+ (+ (aref state q)
+ (if (< 1 r) (ash (aref state (+ q 1)) 8) 0)
+ (if (= 3 r) (ash (aref state (+ q 2)) 16) 0)))))
+ (seed-random-state y)))
+ ;; Also for convenience, we accept non-negative integers as seeds.
+ ;; Small ones get passed to init-random-state, as before.
+ ((unsigned-byte 32)
+ (/show0 "getting random seed from 32-bit word")
+ (%make-random-state :state (init-random-state state)))
+ ;; Larger ones ones get trivially chopped into an array of (unsigned-byte 32)
+ ((unsigned-byte)
+ (/show0 "getting random seed from bignum (converting to 32-bit-word vector)")
+ (loop with l = (ceiling (integer-length state) 32)
+ with s = (make-array (list l) :element-type '(unsigned-byte 32))
+ for i below l
+ for p from 0 by 32
+ do (setf (aref s i) (ldb (byte 32 p) state))
+ finally (return (seed-random-state s))))
+ ;; Last but not least, when provided an array of 32-bit words, we truncate
+ ;; it to 19968 bits and mix these into an initial state. We reuse the same
+ ;; method as the authors of the original algorithm. See
+ ;; http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c
+ ;; NB: their mt[i] is our (aref s (+ 3 i))
+ ((simple-array (unsigned-byte 32) (*))
+ (/show0 "getting random seed from 32-bit-word vector")
+ (let ((s (init-random-state 19650218))
+ (i 1) (j 0) (l (length state)))
+ (loop for k downfrom (max mt19937-n l) above 0 do
+ (setf (aref s (+ i 3))
+ (logand #xFFFFFFFF
+ (+ (logxor (aref s (+ i 3))
+ (* 1664525
+ (logxor (aref s (+ i 2))
+ (ash (aref s (+ i 2)) -30))))
+ (aref state j) j))) ;; non-linear
+ (incf i) (when (>= i mt19937-n) (setf (aref s 3) (aref s (+ 2 mt19937-n)) i 1))
+ (incf j) (when (>= j l) (setf j 0)))
+ (loop for k downfrom (1- mt19937-n) above 0 do
+ (setf (aref s (+ i 3))
+ (logand #xFFFFFFFF
+ (- (logxor (aref s (+ i 3))
+ (* 1566083941
+ (logxor (aref s (+ i 2))
+ (ash (aref s (+ i 2)) -30))))
+ i))) ;; non-linear
+ (incf i) (when (>= i mt19937-n) (setf (aref s 3) (aref s (+ 2 mt19937-n)) i 1)))
+ (setf (aref s 3) #x80000000) ;; MSB is 1; assuring non-zero initial array
+ (%make-random-state :state s)))))
\f
;;;; random entries
;;; inclusive.
#!-sb-fluid (declaim (inline random-chunk))
;;; portable implementation
-(defconstant mt19937-n 624)
-(defconstant mt19937-m 397)
-(defconstant mt19937-upper-mask #x80000000)
-(defconstant mt19937-lower-mask #x7fffffff)
-(defconstant mt19937-b #x9D2C5680)
-(defconstant mt19937-c #xEFC60000)
#!-x86
(defun random-mt19937-update (state)
- (declare (type (simple-array (unsigned-byte 32) (627)) state)
+ (declare (type random-state-state state)
(optimize (speed 3) (safety 0)))
(let ((y 0))
(declare (type (unsigned-byte 32) y))
(defun random-chunk (state)
(declare (type random-state state))
(sb!vm::random-mt19937 (random-state-state state)))
+
+#!-sb-fluid (declaim (inline big-random-chunk))
+(defun big-random-chunk (state)
+ (declare (type random-state state))
+ (logior (ash (random-chunk state) 32)
+ (random-chunk state)))
\f
;;; Handle the single or double float case of RANDOM. We generate a
;;; float between 0.0 and 1.0 by clobbering the significand of 1.0
1d0))))
\f
-;;;; random integers
-
-;;; a bitmask M wide enough that (= (LOGAND INCLUSIVE-LIMIT M) INCLUSIVE-LIMIT)
-(declaim (inline %inclusive-random-integer-mask))
-(defun %inclusive-random-integer-mask (inclusive-limit)
- (1- (ash 1 (integer-length inclusive-limit))))
-
-;;; Transform a uniform sample from an at-least-as-large range into a
-;;; random sample in [0,INCLUSIVE-LIMIT] range by throwing away
-;;; too-big samples.
-;;;
-;;; Up through sbcl-1.0.4, the (RANDOM INTEGER) worked by taking (MOD
-;;; RAW-MERSENNE-OUTPUT INTEGER). That introduces enough bias that it
-;;; is unsuitable for some calculations (like the Metropolis Monte
-;;; Carlo simulations that I, WHN, worked on in grad school): in the
-;;; sbcl-1.0.4, the expectation value of a sample was (in the worst
-;;; part of the range) more than 1.0001 times the ideal expectation
-;;; value. Perhaps that was even ANSI-conformant, because ANSI says only
-;;; An approximately uniform choice distribution is used. If
-;;; LIMIT is an integer, each of the possible results occurs
-;;; with (approximate) probability 1/LIMIT.
-;;; But I (WHN again) said "yuck," so these days we try to get a
-;;; truly uniform distribution by translating RAW-MERSENNE-OUTPUT to
-;;; our output using the second method recommended by
-;;; <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/efaq.html>:
-;;; In the rare case that the rounding-off error becomes a problem,
-;;; obtain minimum integer n such that N<=2^n, generate integer
-;;; random numbers, take the most significant n bits, and discard
-;;; those more than or equal to N.
-;;; (The first method recommended there differs from the sbcl-1.0.4
-;;; method: the recommended method gives slightly different weights
-;;; to different integers, but distributes the overweighted integers
-;;; evenly across the range of outputs so that any skew would be
-;;; of order (/ MOST-POSITIVE-MACHINE-WORD), rather than the (/ 10000)
-;;; or so achieved by sbcl-1.0.4. That skew would probably be
-;;; acceptable --- it seems to be exactly the kind of deviation that
-;;; might have been anticipated in the ANSI CL standard. However, that
-;;; recommended calculation involves multiplication and division of
-;;; two-machine-word bignums, which is hard for us to do efficiently
-;;; here. Without special low-level hacking to support such short
-;;; bignums without consing, the accept-reject solution is not only
-;;; more exact, but also likely more efficient.)
-(defmacro %inclusive-random-integer-accept-reject (raw-mersenne-output-expr
- inclusive-limit)
- (with-unique-names (raw-mersenne-output inclusive-limit-once)
- `(loop
- with ,inclusive-limit-once = ,inclusive-limit
- for ,raw-mersenne-output = ,raw-mersenne-output-expr
- until (<= ,raw-mersenne-output ,inclusive-limit-once)
- finally (return ,raw-mersenne-output))))
+;;;; random fixnums
-;;; an UNSIGNED-BYTE of N-CHUNKS chunks sampled from the Mersenne twister
-(declaim (inline %random-chunks))
-(defun %random-chunks (n-chunks state)
- ;; KLUDGE: This algorithm will cons O(N^2) words when constructing
- ;; an N-word result. To do better should be fundamentally easy, we
- ;; just need to do some low-level hack preallocating the bignum and
- ;; writing its words one by one.
- ;;
- ;; Note: The old sbcl-1.0.4 code did its analogous operation using a
- ;; mysterious RANDOM-INTEGER-OVERLAP parameter, "the amount that we
- ;; overlap chunks by when building a large integer to make up for
- ;; the loss of randomness in the low bits." No such thing is called
- ;; for on
- ;; <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/efaq.html>:
- ;; they say there that it's OK just to concatenate words of twister
- ;; output with no overlap. Thus, crossing our fingers and hoping that
- ;; the previous RNG author didn't have some subtle reason to need
- ;; RANDOM-INTEGER-OVERLAP that we know not of, we just concatenate
- ;; chunks.
- (loop repeat n-chunks
- for result = 0 then (logior (ash result n-random-chunk-bits)
- (random-chunk state))
- finally (return result)))
-
-;;; an UNSIGNED-BYTE of N-BITS bits sampled from the Mersenne twister
-(declaim (inline %random-bits))
-(defun %random-bits (n-bits state)
- (multiple-value-bind (n-full-chunks n-extra-bits)
- (floor n-bits n-random-chunk-bits)
- (let ((full-chunks (%random-chunks n-full-chunks state)))
- (if (zerop n-extra-bits)
- full-chunks
- (logior full-chunks
- (ash (logand (random-chunk state)
- (1- (ash 1 n-extra-bits)))
- (* n-full-chunks n-random-chunk-bits)))))))
-
-;;; the guts of (RANDOM (1+ INCLUSIVE-LIMIT))
-(defun %inclusive-random-integer (inclusive-limit state)
- (declare (optimize speed (space 1))) ; to ensure DEFTRANSFORM is enabled
- (%inclusive-random-integer inclusive-limit state))
-
-;;; the guts of RANDOM for the known-to-return-FIXNUM case
-;;;
-;;; We use INCLUSIVE-LIMIT instead of the (exclusive) LIMIT of RANDOM,
-;;; because we want it to be a FIXNUM even for the possibly-common
-;;; case of (RANDOM (1+ MOST-POSITIVE-FIXNUM)). (That case is what
-;;; one might use for generating random hash values, e.g.)
-;;; It also turns out to be just what's needed for INTEGER-LENGTH.
-(declaim (maybe-inline %inclusive-random-fixnum))
-(defun %inclusive-random-fixnum (inclusive-limit state)
- (declare (type (and fixnum unsigned-byte) inclusive-limit))
- (aver (<= inclusive-limit most-positive-random-chunk))
- (let (;; If this calculation needs to be optimized further, a good
- ;; start might be a DEFTRANSFORM which picks off the case of
- ;; constant LIMIT and precomputes the MASK at compile time.
- (mask (%inclusive-random-integer-mask inclusive-limit)))
- (%inclusive-random-integer-accept-reject (logand (random-chunk state) mask)
- inclusive-limit)))
-\f
-;;;; outer, dynamically-typed interface
+;;; Generate and return a pseudo random fixnum less than ARG. To achieve
+;;; equidistribution an accept-reject loop is used.
+;;; No extra effort is made to detect the case of ARG being a power of
+;;; two where rejection is not possible, as the cost of checking for
+;;; this case is the same as doing the rejection test. When ARG is
+;;; larger than (expt 2 N-RANDOM-CHUNK-BITS), which can only happen if
+;;; the random chunk size is half the word size, two random chunks are
+;;; used in each loop iteration, otherwise only one. Finally, the
+;;; rejection probability could often be reduced by not masking the
+;;; chunk but rejecting only values as least as large as the largest
+;;; multiple of ARG that fits in a chunk (or two), but this is not done
+;;; as the speed gains due to needing fewer loop iterations are by far
+;;; outweighted by the cost of the two divisions required (one to find
+;;; the multiplier and one to bring the result into the correct range).
+#!-sb-fluid (declaim (inline %random-fixnum))
+(defun %random-fixnum (arg state)
+ (declare (type (integer 1 #.sb!xc:most-positive-fixnum) arg)
+ (type random-state state))
+ (if (= arg 1)
+ 0
+ (let* ((n-bits (integer-length (1- arg)))
+ (mask (1- (ash 1 n-bits))))
+ (macrolet ((accept-reject-loop (generator)
+ `(loop
+ (let ((bits (logand mask (,generator state))))
+ (when (< bits arg)
+ (return bits))))))
+ (assert (<= n-bits (* 2 n-random-chunk-bits)))
+ (if (<= n-bits n-random-chunk-bits)
+ (accept-reject-loop random-chunk)
+ (accept-reject-loop big-random-chunk))))))
(defun random (arg &optional (state *random-state*))
- (declare (inline %random-single-float
- %random-double-float
+ (declare (inline %random-fixnum %random-single-float %random-double-float
#!+long-float %random-long-float))
(cond
- ((and (fixnump arg) (plusp arg))
- (locally
- ;; The choice to inline this very common case of
- ;; %INCLUSIVE-RANDOM-FIXNUM and not the less-common call
- ;; below was just a guess (by WHN 2007-03-27), not based on
- ;; benchmarking or anything.
- (declare (inline %inclusive-random-fixnum))
- (%inclusive-random-fixnum (1- arg) state)))
+ ((and (fixnump arg) (> arg 0))
+ (%random-fixnum arg state))
((and (typep arg 'single-float) (> arg 0.0f0))
(%random-single-float arg state))
((and (typep arg 'double-float) (> arg 0.0d0))
#!+long-float
((and (typep arg 'long-float) (> arg 0.0l0))
(%random-long-float arg state))
- ((and (integerp arg) (plusp arg))
- (if (= arg (1+ most-positive-fixnum))
- (%inclusive-random-fixnum most-positive-fixnum state)
- (%inclusive-random-integer (1- arg) state)))
+ ((and (bignump arg) (> arg 0))
+ (%random-bignum arg state))
(t
(error 'simple-type-error
:expected-type '(or (integer 1) (float (0))) :datum arg
projects / sbcl.git / blobdiff