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,
8 ;;;; This software is part of the SBCL system. See the README file for
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.
17 (in-package "SB!KERNEL")
21 (def!method make-load-form ((random-state random-state) &optional environment)
22 (make-load-form-saving-slots random-state :environment environment))
24 ;;; The state is stored in a (simple-array (unsigned-byte 32) (627))
25 ;;; wrapped in a random-state structure:
27 ;;; 0-1: Constant matrix A. [0, #x9908b0df]
31 ;;; Generate and initialize a new random-state array. Index is
32 ;;; initialized to 1 and the states to 32bit integers excluding zero.
34 ;;; Seed - A 32bit number, not zero.
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)
47 (declare (type (mod 625) k))
48 (setf (aref state (+ 3 k))
49 (logand (* 69069 (aref state (+ 3 (1- k)))) #xffffffff)))
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"))
58 (defun make-random-state (&optional state)
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
64 (flet ((copy-random-state (state)
65 (let ((state (random-state-state state))
67 (make-array 627 :element-type '(unsigned-byte 32))))
69 (setf (aref new-state i) (aref state i)))
70 (%make-random-state :state new-state))))
71 (cond ((not state) (copy-random-state *random-state*))
72 ((random-state-p state) (copy-random-state state))
74 (%make-random-state :state (init-random-state
75 (logand (get-universal-time)
77 ;; FIXME: should be TYPE-ERROR?
78 (t (error "Argument is not a RANDOM-STATE, T or NIL: ~S" state)))))
82 ;;; This function generates a 32bit integer between 0 and #xffffffff
84 #!-sb-fluid (declaim (inline random-chunk))
85 ;;; portable implementation
86 (defconstant mt19937-n 624)
87 (defconstant mt19937-m 397)
88 (defconstant mt19937-upper-mask #x80000000)
89 (defconstant mt19937-lower-mask #x7fffffff)
90 (defconstant mt19937-b #x9D2C5680)
91 (defconstant mt19937-c #xEFC60000)
93 (defun random-mt19937-update (state)
94 (declare (type (simple-array (unsigned-byte 32) (627)) state)
95 (optimize (speed 3) (safety 0)))
97 (declare (type (unsigned-byte 32) y))
99 ((>= kk (+ 3 (- mt19937-n mt19937-m))))
100 (declare (type (mod 628) kk))
101 (setf y (logior (logand (aref state kk) mt19937-upper-mask)
102 (logand (aref state (1+ kk)) mt19937-lower-mask)))
103 (setf (aref state kk) (logxor (aref state (+ kk mt19937-m))
104 (ash y -1) (aref state (logand y 1)))))
105 (do ((kk (+ (- mt19937-n mt19937-m) 3) (1+ kk)))
106 ((>= kk (+ (1- mt19937-n) 3)))
107 (declare (type (mod 628) kk))
108 (setf y (logior (logand (aref state kk) mt19937-upper-mask)
109 (logand (aref state (1+ kk)) mt19937-lower-mask)))
110 (setf (aref state kk) (logxor (aref state (+ kk (- mt19937-m mt19937-n)))
111 (ash y -1) (aref state (logand y 1)))))
112 (setf y (logior (logand (aref state (+ 3 (1- mt19937-n)))
114 (logand (aref state 3) mt19937-lower-mask)))
115 (setf (aref state (+ 3 (1- mt19937-n)))
116 (logxor (aref state (+ 3 (1- mt19937-m)))
117 (ash y -1) (aref state (logand y 1)))))
120 (defun random-chunk (state)
121 (declare (type random-state state))
122 (let* ((state (random-state-state state))
124 (declare (type (mod 628) k))
125 (when (= k mt19937-n)
126 (random-mt19937-update state)
128 (setf (aref state 2) (1+ k))
129 (let ((y (aref state (+ 3 k))))
130 (declare (type (unsigned-byte 32) y))
131 (setf y (logxor y (ash y -11)))
132 (setf y (logxor y (ash (logand y (ash mt19937-b -7)) 7)))
133 (setf y (logxor y (ash (logand y (ash mt19937-c -15)) 15)))
134 (setf y (logxor y (ash y -18)))
137 ;;; Using inline VOP support, only available on the x86 so far.
139 ;;; FIXME: It would be nice to have some benchmark numbers on this.
140 ;;; My inclination is to get rid of the nonportable implementation
141 ;;; unless the performance difference is just enormous.
143 (defun random-chunk (state)
144 (declare (type random-state state))
145 (sb!vm::random-mt19937 (random-state-state state)))
147 ;;; Handle the single or double float case of RANDOM. We generate a
148 ;;; float between 0.0 and 1.0 by clobbering the significand of 1.0
149 ;;; with random bits, then subtracting 1.0. This hides the fact that
150 ;;; we have a hidden bit.
151 #!-sb-fluid (declaim (inline %random-single-float %random-double-float))
152 (declaim (ftype (function ((single-float (0f0)) random-state)
154 %random-single-float))
155 (defun %random-single-float (arg state)
156 (declare (type (single-float (0f0)) arg)
157 (type random-state state))
159 (- (make-single-float
160 (dpb (ash (random-chunk state)
161 (- sb!vm:single-float-digits random-chunk-length))
162 sb!vm:single-float-significand-byte
163 (single-float-bits 1.0)))
165 (declaim (ftype (function ((double-float (0d0)) random-state)
167 %random-double-float))
171 (defun %random-double-float (arg state)
172 (declare (type (double-float (0d0)) arg)
173 (type random-state state))
174 (* (float (random-chunk state) 1d0) (/ 1d0 (expt 2 32))))
178 (defun %random-double-float (arg state)
179 (declare (type (double-float (0d0)) arg)
180 (type random-state state))
182 (- (sb!impl::make-double-float
183 (dpb (ash (random-chunk state)
184 (- sb!vm:double-float-digits random-chunk-length
186 sb!vm:double-float-significand-byte
187 (sb!impl::double-float-high-bits 1d0))
188 (random-chunk state))
191 ;;; using a faster inline VOP
193 (defun %random-double-float (arg state)
194 (declare (type (double-float (0d0)) arg)
195 (type random-state state))
196 (let ((state-vector (random-state-state state)))
198 (- (sb!impl::make-double-float
199 (dpb (ash (sb!vm::random-mt19937 state-vector)
200 (- sb!vm:double-float-digits random-chunk-length
202 sb!vm:double-float-significand-byte
203 (sb!impl::double-float-high-bits 1d0))
204 (sb!vm::random-mt19937 state-vector))
208 (declaim #!-sb-fluid (inline %random-long-float))
210 (declaim (ftype (function ((long-float (0l0)) random-state) (long-float 0l0))
213 ;;; using a faster inline VOP
214 #!+(and long-float x86)
215 (defun %random-long-float (arg state)
216 (declare (type (long-float (0l0)) arg)
217 (type random-state state))
218 (let ((state-vector (random-state-state state)))
220 (- (sb!impl::make-long-float
221 (sb!impl::long-float-exp-bits 1l0)
222 (logior (sb!vm::random-mt19937 state-vector)
223 sb!vm:long-float-hidden-bit)
224 (sb!vm::random-mt19937 state-vector))
227 #!+(and long-float sparc)
228 (defun %random-long-float (arg state)
229 (declare (type (long-float (0l0)) arg)
230 (type random-state state))
232 (- (sb!impl::make-long-float
233 (sb!impl::long-float-exp-bits 1l0) ; X needs more work
234 (random-chunk state) (random-chunk state) (random-chunk state))
239 (defun %random-integer (arg state)
240 (declare (type (integer 1) arg) (type random-state state))
241 (let ((shift (- random-chunk-length random-integer-overlap)))
242 (do ((bits (random-chunk state)
243 (logxor (ash bits shift) (random-chunk state)))
244 (count (+ (integer-length arg)
245 (- random-integer-extra-bits shift))
249 (declare (fixnum count)))))
251 (defun random (arg &optional (state *random-state*))
253 "Generate a uniformly distributed pseudo-random number between zero
254 and Arg. State, if supplied, is the random state to use."
255 (declare (inline %random-single-float %random-double-float
256 #!+long-float %long-float))
258 ((and (fixnump arg) (<= arg random-fixnum-max) #!+high-security (> arg 0))
259 (rem (random-chunk state) arg))
260 ((and (typep arg 'single-float) #!+high-security (> arg 0.0S0))
261 (%random-single-float arg state))
262 ((and (typep arg 'double-float) #!+high-security (> arg 0.0D0))
263 (%random-double-float arg state))
265 ((and (typep arg 'long-float) #!+high-security (> arg 0.0L0))
266 (%random-long-float arg state))
267 ((and (integerp arg) #!+high-security (> arg 0))
268 (%random-integer arg state))
270 (error 'simple-type-error
271 :expected-type '(or (integer 1) (float (0))) :datum arg
272 :format-control "Argument is not a positive integer or a positive float: ~S"
273 :format-arguments (list arg)))))