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 (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)"
34 ',(coerce (random-state-state state) 'list)))))
36 ;;; The state is stored in a (simple-array (unsigned-byte 32) (627))
37 ;;; wrapped in a random-state structure:
39 ;;; 0-1: Constant matrix A. [0, #x9908b0df]
43 ;;; Generate and initialize a new random-state array. Index is
44 ;;; initialized to 1 and the states to 32bit integers excluding zero.
46 ;;; Seed - A 32bit number, not zero.
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)
59 (declare (type (mod 625) k))
60 (setf (aref state (+ 3 k))
61 (logand (* 69069 (aref state (+ 3 (1- k)))) #xffffffff)))
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"))
70 (defun make-random-state (&optional state)
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
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))
81 (make-array 627 :element-type '(unsigned-byte 32))))
82 (/show0 "made NEW-STATE, about to DOTIMES")
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")
91 (copy-random-state *random-state*))
93 (/show0 "RANDOM-STATE-P clause")
94 (copy-random-state state))
97 (%make-random-state :state (init-random-state
98 (logand (get-universal-time)
103 ;;; This function generates a 32bit integer between 0 and #xffffffff
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)
114 (defun random-mt19937-update (state)
115 (declare (type (simple-array (unsigned-byte 32) (627)) state)
116 (optimize (speed 3) (safety 0)))
118 (declare (type (unsigned-byte 32) y))
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)))
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)))))
141 (defun random-chunk (state)
142 (declare (type random-state state))
143 (let* ((state (random-state-state state))
145 (declare (type (mod 628) k))
146 (when (= k mt19937-n)
147 (random-mt19937-update state)
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)))
158 ;;; Using inline VOP support, only available on the x86 so far.
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.
164 (defun random-chunk (state)
165 (declare (type random-state state))
166 (sb!vm::random-mt19937 (random-state-state state)))
168 ;;; Handle the single or double float case of RANDOM. We generate a
169 ;;; float between 0.0 and 1.0 by clobbering the significand of 1.0
170 ;;; with random bits, then subtracting 1.0. This hides the fact that
171 ;;; we have a hidden bit.
172 #!-sb-fluid (declaim (inline %random-single-float %random-double-float))
173 (declaim (ftype (function ((single-float (0f0)) random-state)
175 %random-single-float))
176 (defun %random-single-float (arg state)
177 (declare (type (single-float (0f0)) arg)
178 (type random-state state))
180 (- (make-single-float
181 (dpb (ash (random-chunk state)
182 (- sb!vm:single-float-digits n-random-chunk-bits))
183 sb!vm:single-float-significand-byte
184 (single-float-bits 1.0)))
186 (declaim (ftype (function ((double-float (0d0)) random-state)
188 %random-double-float))
192 (defun %random-double-float (arg state)
193 (declare (type (double-float (0d0)) arg)
194 (type random-state state))
195 (* (float (random-chunk state) 1d0) (/ 1d0 (expt 2 32))))
199 (defun %random-double-float (arg state)
200 (declare (type (double-float (0d0)) arg)
201 (type random-state state))
203 (- (sb!impl::make-double-float
204 (dpb (ash (random-chunk state)
205 (- sb!vm:double-float-digits n-random-chunk-bits 32))
206 sb!vm:double-float-significand-byte
207 (sb!impl::double-float-high-bits 1d0))
208 (random-chunk state))
211 ;;; using a faster inline VOP
213 (defun %random-double-float (arg state)
214 (declare (type (double-float (0d0)) arg)
215 (type random-state state))
216 (let ((state-vector (random-state-state state)))
218 (- (sb!impl::make-double-float
219 (dpb (ash (sb!vm::random-mt19937 state-vector)
220 (- sb!vm:double-float-digits n-random-chunk-bits
222 sb!vm:double-float-significand-byte
223 (sb!impl::double-float-high-bits 1d0))
224 (sb!vm::random-mt19937 state-vector))
230 ;;; a bitmask M wide enough that (= (LOGAND INCLUSIVE-LIMIT M) INCLUSIVE-LIMIT)
231 (declaim (inline %inclusive-random-integer-mask))
232 (defun %inclusive-random-integer-mask (inclusive-limit)
233 (1- (ash 1 (integer-length inclusive-limit))))
235 ;;; Transform a uniform sample from an at-least-as-large range into a
236 ;;; random sample in [0,INCLUSIVE-LIMIT] range by throwing away
239 ;;; Up through sbcl-1.0.4, the (RANDOM INTEGER) worked by taking (MOD
240 ;;; RAW-MERSENNE-OUTPUT INTEGER). That introduces enough bias that it
241 ;;; is unsuitable for some calculations (like the Metropolis Monte
242 ;;; Carlo simulations that I, WHN, worked on in grad school): in the
243 ;;; sbcl-1.0.4, the expectation value of a sample was (in the worst
244 ;;; part of the range) more than 1.0001 times the ideal expectation
245 ;;; value. Perhaps that was even ANSI-conformant, because ANSI says only
246 ;;; An approximately uniform choice distribution is used. If
247 ;;; LIMIT is an integer, each of the possible results occurs
248 ;;; with (approximate) probability 1/LIMIT.
249 ;;; But I (WHN again) said "yuck," so these days we try to get a
250 ;;; truly uniform distribution by translating RAW-MERSENNE-OUTPUT to
251 ;;; our output using the second method recommended by
252 ;;; <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/efaq.html>:
253 ;;; In the rare case that the rounding-off error becomes a problem,
254 ;;; obtain minimum integer n such that N<=2^n, generate integer
255 ;;; random numbers, take the most significant n bits, and discard
256 ;;; those more than or equal to N.
257 ;;; (The first method recommended there differs from the sbcl-1.0.4
258 ;;; method: the recommended method gives slightly different weights
259 ;;; to different integers, but distributes the overweighted integers
260 ;;; evenly across the range of outputs so that any skew would be
261 ;;; of order (/ MOST-POSITIVE-MACHINE-WORD), rather than the (/ 10000)
262 ;;; or so achieved by sbcl-1.0.4. That skew would probably be
263 ;;; acceptable --- it seems to be exactly the kind of deviation that
264 ;;; might have been anticipated in the ANSI CL standard. However, that
265 ;;; recommended calculation involves multiplication and division of
266 ;;; two-machine-word bignums, which is hard for us to do efficiently
267 ;;; here. Without special low-level hacking to support such short
268 ;;; bignums without consing, the accept-reject solution is not only
269 ;;; more exact, but also likely more efficient.)
270 (defmacro %inclusive-random-integer-accept-reject (raw-mersenne-output-expr
272 (with-unique-names (raw-mersenne-output inclusive-limit-once)
274 with ,inclusive-limit-once = ,inclusive-limit
275 for ,raw-mersenne-output = ,raw-mersenne-output-expr
276 until (<= ,raw-mersenne-output ,inclusive-limit-once)
277 finally (return ,raw-mersenne-output))))
279 ;;; an UNSIGNED-BYTE of N-CHUNKS chunks sampled from the Mersenne twister
280 (declaim (inline %random-chunks))
281 (defun %random-chunks (n-chunks state)
282 ;; KLUDGE: This algorithm will cons O(N^2) words when constructing
283 ;; an N-word result. To do better should be fundamentally easy, we
284 ;; just need to do some low-level hack preallocating the bignum and
285 ;; writing its words one by one.
287 ;; Note: The old sbcl-1.0.4 code did its analogous operation using a
288 ;; mysterious RANDOM-INTEGER-OVERLAP parameter, "the amount that we
289 ;; overlap chunks by when building a large integer to make up for
290 ;; the loss of randomness in the low bits." No such thing is called
292 ;; <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/efaq.html>:
293 ;; they say there that it's OK just to concatenate words of twister
294 ;; output with no overlap. Thus, crossing our fingers and hoping that
295 ;; the previous RNG author didn't have some subtle reason to need
296 ;; RANDOM-INTEGER-OVERLAP that we know not of, we just concatenate
298 (loop repeat n-chunks
299 for result = 0 then (logior (ash result n-random-chunk-bits)
300 (random-chunk state))
301 finally (return result)))
303 ;;; an UNSIGNED-BYTE of N-BITS bits sampled from the Mersenne twister
304 (declaim (inline %random-bits))
305 (defun %random-bits (n-bits state)
306 (multiple-value-bind (n-full-chunks n-extra-bits)
307 (floor n-bits n-random-chunk-bits)
308 (let ((full-chunks (%random-chunks n-full-chunks state)))
309 (if (zerop n-extra-bits)
312 (ash (logand (random-chunk state)
313 (1- (ash 1 n-extra-bits)))
314 (* n-full-chunks n-random-chunk-bits)))))))
316 ;;; the guts of (RANDOM (1+ INCLUSIVE-LIMIT))
317 (defun %inclusive-random-integer (inclusive-limit state)
318 (declare (optimize speed (space 1))) ; to ensure DEFTRANSFORM is enabled
319 (%inclusive-random-integer inclusive-limit state))
321 ;;; the guts of RANDOM for the known-to-return-FIXNUM case
323 ;;; We use INCLUSIVE-LIMIT instead of the (exclusive) LIMIT of RANDOM,
324 ;;; because we want it to be a FIXNUM even for the possibly-common
325 ;;; case of (RANDOM (1+ MOST-POSITIVE-FIXNUM)). (That case is what
326 ;;; one might use for generating random hash values, e.g.)
327 ;;; It also turns out to be just what's needed for INTEGER-LENGTH.
328 (declaim (maybe-inline %inclusive-random-fixnum))
329 (defun %inclusive-random-fixnum (inclusive-limit state)
330 (declare (type (and fixnum unsigned-byte) inclusive-limit))
331 (aver (<= inclusive-limit most-positive-random-chunk))
332 (let (;; If this calculation needs to be optimized further, a good
333 ;; start might be a DEFTRANSFORM which picks off the case of
334 ;; constant LIMIT and precomputes the MASK at compile time.
335 (mask (%inclusive-random-integer-mask inclusive-limit)))
336 (%inclusive-random-integer-accept-reject (logand (random-chunk state) mask)
339 ;;;; outer, dynamically-typed interface
341 (defun random (arg &optional (state *random-state*))
342 (declare (inline %random-single-float
344 #!+long-float %random-long-float))
346 ((and (fixnump arg) (plusp arg))
348 ;; The choice to inline this very common case of
349 ;; %INCLUSIVE-RANDOM-FIXNUM and not the less-common call
350 ;; below was just a guess (by WHN 2007-03-27), not based on
351 ;; benchmarking or anything.
352 (declare (inline %inclusive-random-fixnum))
353 (%inclusive-random-fixnum (1- arg) state)))
354 ((and (typep arg 'single-float) (> arg 0.0f0))
355 (%random-single-float arg state))
356 ((and (typep arg 'double-float) (> arg 0.0d0))
357 (%random-double-float arg state))
359 ((and (typep arg 'long-float) (> arg 0.0l0))
360 (%random-long-float arg state))
361 ((and (integerp arg) (plusp arg))
362 (if (= arg (1+ most-positive-fixnum))
363 (%inclusive-random-fixnum most-positive-fixnum state)
364 (%inclusive-random-integer (1- arg) state)))
366 (error 'simple-type-error
367 :expected-type '(or (integer 1) (float (0))) :datum arg
368 :format-control "~@<Argument is neither a positive integer nor a ~
369 positive float: ~2I~_~S~:>"
370 :format-arguments (list arg)))))