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")
24 (def!method make-load-form ((random-state random-state) &optional environment)
25 (make-load-form-saving-slots random-state :environment environment))
27 ;;; The state is stored in a (simple-array (unsigned-byte 32) (627))
28 ;;; wrapped in a random-state structure:
30 ;;; 0-1: Constant matrix A. [0, #x9908b0df]
34 ;;; Generate and initialize a new random-state array. Index is
35 ;;; initialized to 1 and the states to 32bit integers excluding zero.
37 ;;; Seed - A 32bit number, not zero.
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)
50 (declare (type (mod 625) k))
51 (setf (aref state (+ 3 k))
52 (logand (* 69069 (aref state (+ 3 (1- k)))) #xffffffff)))
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"))
61 (defun make-random-state (&optional state)
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
67 (flet ((copy-random-state (state)
68 (let ((state (random-state-state state))
70 (make-array 627 :element-type '(unsigned-byte 32))))
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))
77 (%make-random-state :state (init-random-state
78 (logand (get-universal-time)
80 ;; FIXME: should be TYPE-ERROR?
81 (t (error "Argument is not a RANDOM-STATE, T or NIL: ~S" state)))))
85 ;;; This function generates a 32bit integer between 0 and #xffffffff
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)
96 (defun random-mt19937-update (state)
97 (declare (type (simple-array (unsigned-byte 32) (627)) state)
98 (optimize (speed 3) (safety 0)))
100 (declare (type (unsigned-byte 32) y))
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)))
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)))))
123 (defun random-chunk (state)
124 (declare (type random-state state))
125 (let* ((state (random-state-state state))
127 (declare (type (mod 628) k))
128 (when (= k mt19937-n)
129 (random-mt19937-update state)
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)))
140 ;;; Using inline VOP support, only available on the x86 so far.
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.
146 (defun random-chunk (state)
147 (declare (type random-state state))
148 (sb!vm::random-mt19937 (random-state-state state)))
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)
157 %random-single-float))
158 (defun %random-single-float (arg state)
159 (declare (type (single-float (0f0)) arg)
160 (type random-state state))
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)))
168 (declaim (ftype (function ((double-float (0d0)) random-state)
170 %random-double-float))
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))))
181 (defun %random-double-float (arg state)
182 (declare (type (double-float (0d0)) arg)
183 (type random-state state))
185 (- (sb!impl::make-double-float
186 (dpb (ash (random-chunk state)
187 (- sb!vm:double-float-digits random-chunk-length
189 sb!vm:double-float-significand-byte
190 (sb!impl::double-float-high-bits 1d0))
191 (random-chunk state))
194 ;;; using a faster inline VOP
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)))
201 (- (sb!impl::make-double-float
202 (dpb (ash (sb!vm::random-mt19937 state-vector)
203 (- sb!vm:double-float-digits random-chunk-length
205 sb!vm:double-float-significand-byte
206 (sb!impl::double-float-high-bits 1d0))
207 (sb!vm::random-mt19937 state-vector))
211 (declaim #!-sb-fluid (inline %random-long-float))
213 (declaim (ftype (function ((long-float (0l0)) random-state) (long-float 0l0))
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)))
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))
230 #!+(and long-float sparc)
231 (defun %random-long-float (arg state)
232 (declare (type (long-float (0l0)) arg)
233 (type random-state state))
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))
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))
252 (declare (fixnum count)))))
254 (defun random (arg &optional (state *random-state*))
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))
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))
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))
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)))))