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 (defconstant n-random-mt19937-bits 32)
106 (declaim (inline %random32))
107 ;;; portable implementation
108 (defconstant mt19937-n 624)
109 (defconstant mt19937-m 397)
110 (defconstant mt19937-upper-mask #x80000000)
111 (defconstant mt19937-lower-mask #x7fffffff)
112 (defconstant mt19937-b #x9D2C5680)
113 (defconstant mt19937-c #xEFC60000)
115 (defun random-mt19937-update (state)
116 (declare (type (simple-array (unsigned-byte 32) (627)) state)
117 (optimize (speed 3) (safety 0)))
119 (declare (type (unsigned-byte 32) y))
121 ((>= kk (+ 3 (- mt19937-n mt19937-m))))
122 (declare (type (mod 628) kk))
123 (setf y (logior (logand (aref state kk) mt19937-upper-mask)
124 (logand (aref state (1+ kk)) mt19937-lower-mask)))
125 (setf (aref state kk) (logxor (aref state (+ kk mt19937-m))
126 (ash y -1) (aref state (logand y 1)))))
127 (do ((kk (+ (- mt19937-n mt19937-m) 3) (1+ kk)))
128 ((>= kk (+ (1- mt19937-n) 3)))
129 (declare (type (mod 628) kk))
130 (setf y (logior (logand (aref state kk) mt19937-upper-mask)
131 (logand (aref state (1+ kk)) mt19937-lower-mask)))
132 (setf (aref state kk) (logxor (aref state (+ kk (- mt19937-m mt19937-n)))
133 (ash y -1) (aref state (logand y 1)))))
134 (setf y (logior (logand (aref state (+ 3 (1- mt19937-n)))
136 (logand (aref state 3) mt19937-lower-mask)))
137 (setf (aref state (+ 3 (1- mt19937-n)))
138 (logxor (aref state (+ 3 (1- mt19937-m)))
139 (ash y -1) (aref state (logand y 1)))))
142 (defun %random32 (state)
143 (declare (type random-state state))
144 (let* ((state (random-state-state state))
146 (declare (type (mod 628) k))
147 (when (= k mt19937-n)
148 (random-mt19937-update state)
150 (setf (aref state 2) (1+ k))
151 (let ((y (aref state (+ 3 k))))
152 (declare (type (unsigned-byte 32) y))
153 (setf y (logxor y (ash y -11)))
154 (setf y (logxor y (ash (logand y (ash mt19937-b -7)) 7)))
155 (setf y (logxor y (ash (logand y (ash mt19937-c -15)) 15)))
156 (setf y (logxor y (ash y -18)))
159 ;;; Using inline VOP support, only available on the x86 so far.
161 ;;; FIXME: It would be nice to have some benchmark numbers on this.
162 ;;; My inclination is to get rid of the nonportable implementation
163 ;;; unless the performance difference is just enormous.
165 (defun %random32 (state)
166 (declare (type random-state state))
167 (sb!vm::random-mt19937 (random-state-state state)))
169 (declaim (inline %random-word))
170 (defun %random-word (state)
171 ;; KLUDGE: This #.(ECASE ...) is not the most flexible and elegant
172 ;; construct one could imagine. It is intended as a quick fix to stand
173 ;; in for The Right Thing which can't be coded so quickly.
175 ;; The Right Thing: The Mersenne Twister has been generalized to 64
176 ;; bits, and seems likely to be generalized to any future common CPU
177 ;; word width as well. Thus, it should be straightforward to
178 ;; implement this as "Return one sample from the MT variant
179 ;; corresponding to SB!VM:N-WORD-BITS."
181 ;; Meanwhile: Mock that up by pasting together as many samples from
182 ;; the 32-bit Mersenne Twister as necessary.
183 #.(ecase sb!vm:n-word-bits
184 (32 '(%random32 state))
187 (ash (%random32 state) 32)))))
189 ;;; Handle the single or double float case of RANDOM. We generate a
190 ;;; float between 0.0 and 1.0 by clobbering the significand of 1.0
191 ;;; with random bits, then subtracting 1.0. This hides the fact that
192 ;;; we have a hidden bit.
193 #!-sb-fluid (declaim (inline %random-single-float %random-double-float))
194 (declaim (ftype (function ((single-float (0f0)) random-state)
196 %random-single-float))
197 (defun %random-single-float (arg state)
198 (declare (type (single-float (0f0)) arg)
199 (type random-state state))
200 ;; KLUDGE: The hardwired-to-32 hackery here could be replaced by a
201 ;; call to %RANDOM-BITS if there were sufficiently smart
202 ;; DEFTRANSFORMs for it, but in sbcl-1.0.4.epsilon it looks as
203 ;; though it would be a performance disaster.
204 (aver (<= sb!vm:single-float-digits 32))
206 (- (make-single-float
209 (- sb!vm:single-float-digits 32))
210 sb!vm:single-float-significand-byte
211 (single-float-bits 1.0)))
213 (declaim (ftype (function ((double-float (0d0)) random-state)
215 %random-double-float))
219 (defun %random-double-float (arg state)
220 (declare (type (double-float (0d0)) arg)
221 (type random-state state))
222 ;; KLUDGE: As in %RANDOM-SIMNGLE-FLOAT, as of sbcl-1.0.4.epsilon
223 ;; calling %RANDOM-BITS doesn't look reasonable, so we bang bits.
224 (aver (<= sb!vm:single-float-digits 32))
226 (- (sb!impl::make-double-float
227 (dpb (ash (%random32 state)
228 (- sb!vm:double-float-digits 64))
229 sb!vm:double-float-significand-byte
230 (sb!impl::double-float-high-bits 1d0))
234 ;;; using a faster inline VOP
236 (defun %random-double-float (arg state)
237 (declare (type (double-float (0d0)) arg)
238 (type random-state state))
239 (let ((state-vector (random-state-state state)))
241 (- (sb!impl::make-double-float
242 (dpb (ash (sb!vm::random-mt19937 state-vector)
243 (- sb!vm:double-float-digits
244 n-random-mt19937-bits
246 sb!vm:double-float-significand-byte
247 (sb!impl::double-float-high-bits 1d0))
248 (sb!vm::random-mt19937 state-vector))
253 ;;; a bitmask M wide enough that (= (LOGAND INCLUSIVE-LIMIT M) INCLUSIVE-LIMIT)
254 (declaim (inline %inclusive-random-integer-mask))
255 (defun %inclusive-random-integer-mask (inclusive-limit)
256 (1- (ash 1 (integer-length inclusive-limit))))
258 ;;; Transform a uniform sample from an at-least-as-large range into a
259 ;;; random sample in [0,INCLUSIVE-LIMIT] range by throwing away
262 ;;; Up through sbcl-1.0.4, the (RANDOM INTEGER) worked by taking (MOD
263 ;;; RAW-MERSENNE-OUTPUT INTEGER). That introduces enough bias that it
264 ;;; is unsuitable for some calculations (like the Metropolis Monte
265 ;;; Carlo simulations that I, WHN, worked on in grad school): in the
266 ;;; sbcl-1.0.4, the expectation value of a sample was (in the worst
267 ;;; part of the range) more than 1.0001 times the ideal expectation
268 ;;; value. Perhaps that was even ANSI-conformant, because ANSI says only
269 ;;; An approximately uniform choice distribution is used. If
270 ;;; LIMIT is an integer, each of the possible results occurs
271 ;;; with (approximate) probability 1/LIMIT.
272 ;;; But I (WHN again) said "yuck," so these days we try to get a
273 ;;; truly uniform distribution by translating RAW-MERSENNE-OUTPUT to
274 ;;; our output using the second method recommended by
275 ;;; <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/efaq.html>:
276 ;;; In the rare case that the rounding-off error becomes a problem,
277 ;;; obtain minimum integer n such that N<=2^n, generate integer
278 ;;; random numbers, take the most significant n bits, and discard
279 ;;; those more than or equal to N.
280 ;;; (The first method recommended there differs from the sbcl-1.0.4
281 ;;; method: the recommended method gives slightly different weights
282 ;;; to different integers, but distributes the overweighted integers
283 ;;; evenly across the range of outputs so that any skew would be
284 ;;; of order (/ MOST-POSITIVE-MACHINE-WORD), rather than the (/ 10000)
285 ;;; or so achieved by sbcl-1.0.4. That skew would probably be
286 ;;; acceptable --- it seems to be exactly the kind of deviation that
287 ;;; might have been anticipated in the ANSI CL standard. However, that
288 ;;; recommended calculation involves multiplication and division of
289 ;;; two-machine-word bignums, which is hard for us to do efficiently
290 ;;; here. Without special low-level hacking to support such short
291 ;;; bignums without consing, the accept-reject solution is not only
292 ;;; more exact, but also likely more efficient.)
293 (defmacro %inclusive-random-integer-accept-reject (raw-mersenne-output-expr
295 (with-unique-names (raw-mersenne-output inclusive-limit-once)
297 with ,inclusive-limit-once = ,inclusive-limit
298 for ,raw-mersenne-output = ,raw-mersenne-output-expr
299 until (<= ,raw-mersenne-output ,inclusive-limit-once)
300 finally (return ,raw-mersenne-output))))
302 ;;; an UNSIGNED-BYTE of N-WORDS words sampled from the Mersenne twister
303 (declaim (inline %random-words))
304 (defun %random-words (n-words state)
305 ;; KLUDGE: This algorithm will cons O(N^2) words when constructing
306 ;; an N-word result. To do better should be fundamentally easy, we
307 ;; just need to do some low-level hack preallocating the bignum and
308 ;; writing its words one by one.
310 ;; Note: The old sbcl-1.0.4 code did its analogous operation using a
311 ;; mysterious RANDOM-INTEGER-OVERLAP parameter, "the amount that we
312 ;; overlap chunks by when building a large integer to make up for
313 ;; the loss of randomness in the low bits." No such thing is called
315 ;; <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/efaq.html>:
316 ;; they say there that it's OK just to concatenate words of twister
317 ;; output with no overlap. Thus, crossing our fingers and hoping that
318 ;; the previous RNG author didn't have some subtle reason to need
319 ;; RANDOM-INTEGER-OVERLAP that we know not of, we just concatenate
322 for result = 0 then (logior (ash result sb!vm:n-word-bits)
323 (%random-word state))
324 finally (return result)))
326 ;;; an UNSIGNED-BYTE of N-BITS bits sampled from the Mersenne twister
327 (declaim (inline %random-bits))
328 (defun %random-bits (n-bits state)
329 (multiple-value-bind (n-full-words n-extra-bits)
330 (floor n-bits sb!vm:n-word-bits)
331 (let ((full-chunks (%random-words n-full-words state)))
332 (if (zerop n-extra-bits)
335 (ash (logand (%random-word state)
336 (1- (ash 1 n-extra-bits)))
338 sb!vm:n-word-bits)))))))
340 ;;; the guts of (RANDOM (1+ INCLUSIVE-LIMIT))
341 (defun %inclusive-random-integer (inclusive-limit state)
342 (declare (optimize speed (space 1))) ; to ensure DEFTRANSFORM is enabled
343 (%inclusive-random-integer inclusive-limit state))
345 ;;; the guts of RANDOM for the known-to-return-FIXNUM case
347 ;;; We use INCLUSIVE-LIMIT instead of the (exclusive) LIMIT of RANDOM,
348 ;;; because we want it to be a FIXNUM even for the possibly-common
349 ;;; case of (RANDOM (1+ MOST-POSITIVE-FIXNUM)). (That case is what
350 ;;; one might use for generating random hash values, e.g.)
351 ;;; It also turns out to be just what's needed for INTEGER-LENGTH.
352 (declaim (maybe-inline %inclusive-random-fixnum))
353 (defun %inclusive-random-fixnum (inclusive-limit state)
354 (declare (type (and fixnum unsigned-byte) inclusive-limit))
355 (let (;; If this calculation needs to be optimized further, a good
356 ;; start might be a DEFTRANSFORM which picks off the case of
357 ;; constant LIMIT and precomputes the MASK at compile time.
358 (mask (%inclusive-random-integer-mask inclusive-limit)))
359 (%inclusive-random-integer-accept-reject (logand (%random-word state) mask)
362 ;;;; outer, dynamically-typed interface
364 (defun random (arg &optional (state *random-state*))
365 (declare (inline %random-single-float
367 #!+long-float %random-long-float))
369 ((and (fixnump arg) (plusp arg))
371 ;; The choice to inline this very common case of
372 ;; %INCLUSIVE-RANDOM-FIXNUM and not the less-common call
373 ;; below was just a guess (by WHN 2007-03-27), not based on
374 ;; benchmarking or anything.
375 (declare (inline %inclusive-random-fixnum))
376 (%inclusive-random-fixnum (1- arg) state)))
377 ((and (typep arg 'single-float) (> arg 0.0f0))
378 (%random-single-float arg state))
379 ((and (typep arg 'double-float) (> arg 0.0d0))
380 (%random-double-float arg state))
382 ((and (typep arg 'long-float) (> arg 0.0l0))
383 (%random-long-float arg state))
384 ((and (integerp arg) (plusp arg))
385 (if (= arg (1+ most-positive-fixnum))
386 (%inclusive-random-fixnum most-positive-fixnum state)
387 (%inclusive-random-integer (1- arg) state)))
389 (error 'simple-type-error
390 :expected-type '(or (integer 1) (float (0))) :datum arg
391 :format-control "~@<Argument is neither a positive integer nor a ~
392 positive float: ~2I~_~S~:>"
393 :format-arguments (list arg)))))