/[cmucl]/src/code/rand-mt19937.lisp
ViewVC logotype

Contents of /src/code/rand-mt19937.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.23 - (show annotations)
Tue Apr 20 17:57:45 2010 UTC (3 years, 11 months ago) by rtoy
Branch: MAIN
CVS Tags: sparc-tramp-assem-base, release-20b-pre1, release-20b-pre2, sparc-tramp-assem-2010-07-19, GIT-CONVERSION, cross-sol-x86-merged, RELEASE_20b, cross-sol-x86-base, snapshot-2010-12, snapshot-2010-11, snapshot-2011-09, snapshot-2011-06, snapshot-2011-07, snapshot-2011-04, snapshot-2011-02, snapshot-2011-03, snapshot-2011-01, snapshot-2010-05, snapshot-2010-07, snapshot-2010-06, snapshot-2010-08, cross-sol-x86-2010-12-20, cross-sparc-branch-base, HEAD
Branch point for: cross-sparc-branch, RELEASE-20B-BRANCH, sparc-tramp-assem-branch, cross-sol-x86-branch
Changes since 1.22: +3 -3 lines
Change uses of _"foo" to (intl:gettext "foo").  This is because slime
may get confused with source locations if the reader macros are
installed.
1 ;;; -*- Mode: Lisp; Package: Kernel -*-
2 ;;;
3 ;;; **********************************************************************
4 ;;; This code was written by Douglas T. Crosher and Raymond Toy based
5 ;;; on public domain code from Carnegie Mellon University and has been
6 ;;; placed in the Public domain, and is provided 'as is'.
7 ;;;
8 (ext:file-comment
9 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/rand-mt19937.lisp,v 1.23 2010/04/20 17:57:45 rtoy Rel $")
10 ;;;
11 ;;; **********************************************************************
12 ;;;
13 ;;; Support for the Mersenne Twister, MT19937, random number generator
14 ;;; due to Matsumoto and Nishimura. This implementation has been
15 ;;; placed in the public domain with permission from M. Matsumoto.
16 ;;;
17 ;;; Makoto Matsumoto and T. Nishimura, "Mersenne twister: A
18 ;;; 623-dimensionally equidistributed uniform pseudorandom number
19 ;;; generator.", ACM Transactions on Modeling and Computer Simulation,
20 ;;; 1997, to appear.
21
22 (in-package "LISP")
23 (intl:textdomain "cmucl")
24
25 (export '(random-state random-state-p random *random-state*
26 make-random-state))
27
28 (in-package "KERNEL")
29 (export '(%random-single-float %random-double-float random-chunk init-random-state))
30
31 (sys:register-lisp-feature :random-mt19937)
32
33
34 ;;;; Random state hackery:
35
36 ;;; The state is stored in a (simple-array (unsigned-byte 32) (627))
37 ;;; wrapped in a random-state structure:
38 ;;;
39 ;;; 0-1: Constant matrix A. [0, #x9908b0df]
40 ;;; 2: Index k.
41 ;;; 3-626: State.
42
43 ;; GENERATE-SEED
44 ;;
45 ;; Generate a random seed that can be used for seeding the generator.
46 ;; If /dev/urandom is available, it is used to generate random data as
47 ;; the seed. Otherwise, the current time is used as the seed.
48 ;;
49 ;; The /dev/urandom device exists on Linux, FreeBSD, Solaris 8 and
50 ;; later. (You need to have patch 112438-01 for Solaris 8 to get that
51 ;; device, though). It returns pseudorandom data with entropy that the
52 ;; kernel collects from the environment. Unlike the related
53 ;; /dev/random, this device does not block when the entropy pool has
54 ;; been depleted.
55 (defun generate-seed (&optional (nwords 1))
56 ;; On some systems (as reported by Ole Rohne on cmucl-imp),
57 ;; /dev/urandom isn't what we think it is, so if it doesn't work,
58 ;; silently generate the seed from the current time.
59 (or (ignore-errors
60 (let ((words (make-array nwords :element-type '(unsigned-byte 32))))
61 (with-open-file (rand "/dev/urandom"
62 :direction :input
63 :element-type '(unsigned-byte 32))
64 (read-sequence words rand))
65 (if (= nwords 1)
66 (aref words 0)
67 words)))
68 (logand (get-universal-time) #xffffffff)))
69
70 ;; New initializer proposed by Takuji Nishimura and Makota Matsumoto.
71 ;; (See http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html)
72 ;;
73 ;; This corrects a deficiency in the original initializer wherein the
74 ;; MSB of the seed was not well represented in the state.
75 ;;
76 ;; The initialization routine is described below. Let s be the seed,
77 ;; mt[] be the state vector. Then the algorithm is
78 ;;
79 ;; mt[0] = s & 0xffffffffUL
80 ;;
81 ;; for (k = 1; k < N; k++) {
82 ;; mt[k] = 1812433253 * (mt[k-1] ^ (mt[k-1] >> 30)) + k
83 ;; mt[k] &= 0xffffffffUL
84 ;; }
85 ;;
86 ;; The multiplier is from Knuth TAOCP Vol2, 3rd Ed., p. 106.
87 ;;
88
89 (defun int-init-random-state (&optional (seed 5489) state)
90 (declare (type (integer 0 #xffffffff) seed))
91 (let ((state (or state (make-array 627 :element-type '(unsigned-byte 32)))))
92 (declare (type (simple-array (unsigned-byte 32) (627)) state))
93 (setf (aref state 0) 0)
94 (setf (aref state 1) #x9908b0df)
95 (setf (aref state 2) mt19937-n)
96 (setf (aref state 3) seed)
97 (do ((k 1 (1+ k)))
98 ((>= k 624))
99 (declare (type (mod 625) k))
100 (let ((prev (aref state (+ 3 (1- k)))))
101 (setf (aref state (+ 3 k))
102 (logand (+ (* 1812433253 (logxor prev (ash prev -30)))
103 k)
104 #xffffffff))))
105 state))
106
107 ;; Initialize from an array.
108 ;;
109 ;; Here is the algorithm, in C. init_genrand is the initalizer above,
110 ;; init_key is the seed vector of length key_length.
111 ;;
112 ;; init_genrand(19650218UL);
113 ;; i=1; j=0;
114 ;; k = (N>key_length ? N : key_length);
115 ;; for (; k; k--) {
116 ;; mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
117 ;; + init_key[j] + j; /* non linear */
118 ;; mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
119 ;; i++; j++;
120 ;; if (i>=N) {
121 ;; mt[0] = mt[N-1]; i=1;
122 ;; }
123 ;; if (j>=key_length) {
124 ;; j=0;
125 ;; }
126 ;; }
127 ;; for (k=N-1; k; k--) {
128 ;; mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
129 ;; - i; /* non linear */
130 ;; mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
131 ;; i++;
132 ;; if (i>=N) { mt[0] = mt[N-1]; i=1; }
133 ;; }
134 ;;
135 ;; mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
136 ;;
137
138 (defun vec-init-random-state (key &optional state)
139 (declare (type (array (unsigned-byte 32) (*)) key))
140 (let ((key-len (length key))
141 (state (init-random-state 19650218 state))
142 (i 1)
143 (j 0))
144 (loop for k from (max key-len mt19937-n) above 0 do
145 (let ((prev (aref state (+ 3 (1- i)))))
146 (setf (aref state (+ 3 i))
147 (ldb (byte 32 0)
148 (+ (aref key j) j
149 (logxor (aref state (+ 3 i))
150 (ldb (byte 32 0)
151 (* 1664525
152 (logxor prev (ash prev -30))))))))
153 (incf i)
154 (incf j)
155 (when (>= i mt19937-n)
156 (setf (aref state 3)
157 (aref state (+ 3 (- mt19937-n 1))))
158 (setf i 1))
159 (when (>= j key-len)
160 (setf j 0))))
161
162 (loop for k from (1- mt19937-n) above 0 do
163 (let ((prev (aref state (+ 3 (1- i)))))
164 (setf (aref state (+ 3 i))
165 (ldb (byte 32 0)
166 (- (logxor (aref state (+ 3 i))
167 (* 1566083941
168 (logxor prev (ash prev -30))))
169 i)))
170 (incf i)
171 (when (>= i mt19937-n)
172 (setf (aref state 3)
173 (aref state (+ 3 (- mt19937-n 1))))
174 (setf i 1))))
175 (setf (aref state 3) #x80000000)
176 state))
177
178 ;;
179 (defun init-random-state (&optional (seed 5489) state)
180 "Generate an random state vector from the given SEED. The seed can be
181 either an integer or a vector of (unsigned-byte 32)"
182 (declare (type (or null integer
183 (array (unsigned-byte 32) (*)))
184 seed))
185 (etypecase seed
186 (integer
187 (int-init-random-state (ldb (byte 32 0) seed) state))
188 ((array (unsigned-byte 32) (*))
189 (vec-init-random-state seed state))))
190
191 (defstruct (random-state
192 (:constructor make-random-object)
193 (:make-load-form-fun :just-dump-it-normally))
194 (state (init-random-state) :type (simple-array (unsigned-byte 32) (627))))
195
196 (defvar *random-state* (make-random-object))
197
198 (defun make-random-state (&optional state)
199 "Make a random state object. If STATE is not supplied, return a copy
200 of the default random state. If STATE is a random state, then return a
201 copy of it. If STATE is T then return a random state generated from
202 the universal time or /dev/urandom if available."
203 (flet ((copy-random-state (state)
204 (let ((state (random-state-state state))
205 (new-state
206 (make-array 627 :element-type '(unsigned-byte 32))))
207 (dotimes (i 627)
208 (setf (aref new-state i) (aref state i)))
209 (make-random-object :state new-state))))
210 (cond ((not state) (copy-random-state *random-state*))
211 ((random-state-p state) (copy-random-state state))
212 ((eq state t)
213 (make-random-object :state (init-random-state (generate-seed 627))))
214 (t (error (intl:gettext "Argument is not a RANDOM-STATE, T or NIL: ~S") state)))))
215
216 (defun rand-mt19937-initializer ()
217 (init-random-state (generate-seed)
218 (random-state-state *random-state*)))
219
220 (pushnew 'rand-mt19937-initializer ext:*after-save-initializations*)
221
222
223 ;;;; Random entries:
224
225 ;;; Size of the chunks returned by random-chunk.
226 ;;;
227 (defconstant random-chunk-length 32)
228
229 ;;; random-chunk -- Internal
230 ;;;
231 ;;; This function generaters a 32bit integer between 0 and #xffffffff
232 ;;; inclusive.
233 ;;;
234 (declaim (inline random-chunk))
235 ;;;
236 ;;; Portable implementation.
237 (defconstant mt19937-n 624)
238 (defconstant mt19937-m 397)
239 (defconstant mt19937-upper-mask #x80000000)
240 (defconstant mt19937-lower-mask #x7fffffff)
241 (defconstant mt19937-b #x9D2C5680)
242 (defconstant mt19937-c #xEFC60000)
243 ;;;
244 #-x86
245 (defun random-mt19937-update (state)
246 (declare (type (simple-array (unsigned-byte 32) (627)) state)
247 (optimize (speed 3) (safety 0)))
248 (let ((y 0))
249 (declare (type (unsigned-byte 32) y))
250 (do ((kk 3 (1+ kk)))
251 ((>= kk (+ 3 (- mt19937-n mt19937-m))))
252 (declare (type (mod 628) kk))
253 (setf y (logior (logand (aref state kk) mt19937-upper-mask)
254 (logand (aref state (1+ kk)) mt19937-lower-mask)))
255 (setf (aref state kk) (logxor (aref state (+ kk mt19937-m))
256 (ash y -1) (aref state (logand y 1)))))
257 (do ((kk (+ (- mt19937-n mt19937-m) 3) (1+ kk)))
258 ((>= kk (+ (1- mt19937-n) 3)))
259 (declare (type (mod 628) kk))
260 (setf y (logior (logand (aref state kk) mt19937-upper-mask)
261 (logand (aref state (1+ kk)) mt19937-lower-mask)))
262 (setf (aref state kk) (logxor (aref state (+ kk (- mt19937-m mt19937-n)))
263 (ash y -1) (aref state (logand y 1)))))
264 (setf y (logior (logand (aref state (+ 3 (1- mt19937-n)))
265 mt19937-upper-mask)
266 (logand (aref state 3) mt19937-lower-mask)))
267 (setf (aref state (+ 3 (1- mt19937-n)))
268 (logxor (aref state (+ 3 (1- mt19937-m)))
269 (ash y -1) (aref state (logand y 1)))))
270 (values))
271 ;;;
272 #-x86
273 (defun random-chunk (state)
274 (declare (type random-state state)
275 (optimize (speed 3) (safety 0)))
276 (let* ((state (random-state-state state))
277 (k (aref state 2)))
278 (declare (type (mod 628) k))
279 (when (= k mt19937-n)
280 (random-mt19937-update state)
281 (setf k 0))
282 (setf (aref state 2) (1+ k))
283 (let ((y (aref state (+ 3 k))))
284 (declare (type (unsigned-byte 32) y))
285 (setf y (logxor y (ash y -11)))
286 (setf y (logxor y (ash (logand y (ash mt19937-b -7)) 7)))
287 (setf y (logxor y (ash (logand y (ash mt19937-c -15)) 15)))
288 (setf y (logxor y (ash y -18)))
289 y)))
290
291 ;;; Using inline VOP support, only available on the x86 so far.
292 #+x86
293 (defun random-chunk (state)
294 (declare (type random-state state))
295 (vm::random-mt19937 (random-state-state state)))
296
297
298
299 ;;; %RANDOM-SINGLE-FLOAT, %RANDOM-DOUBLE-FLOAT -- Interface
300 ;;;
301 ;;; Handle the single or double float case of RANDOM. We generate a float
302 ;;; between 0.0 and 1.0 by clobbering the significand of 1.0 with random bits,
303 ;;; then subtracting 1.0. This hides the fact that we have a hidden bit.
304 ;;;
305 (declaim (inline %random-single-float %random-double-float))
306 (declaim (ftype (function ((single-float (0f0)) random-state)
307 (single-float 0f0))
308 %random-single-float))
309 ;;;
310 (defun %random-single-float (arg state)
311 (declare (type (single-float (0f0)) arg)
312 (type random-state state))
313 (* arg
314 (- (make-single-float
315 (dpb (ash (random-chunk state)
316 (- vm:single-float-digits random-chunk-length))
317 vm:single-float-significand-byte
318 (single-float-bits 1.0)))
319 1.0)))
320 ;;;
321 (declaim (ftype (function ((double-float (0d0)) random-state)
322 (double-float 0d0))
323 %random-double-float))
324 ;;;
325 ;;; 32bit version
326 ;;;
327 #+nil
328 (defun %random-double-float (arg state)
329 (declare (type (double-float (0d0)) arg)
330 (type random-state state))
331 (* (float (random-chunk state) 1d0) (/ 1d0 (expt 2 32))))
332 ;;;
333 ;;; 53bit version.
334 ;;;
335 #-x86
336 (defun %random-double-float (arg state)
337 (declare (type (double-float (0d0)) arg)
338 (type random-state state))
339 (* arg
340 (- (lisp::make-double-float
341 (dpb (ash (random-chunk state)
342 (- vm:double-float-digits random-chunk-length
343 vm:word-bits))
344 vm:double-float-significand-byte
345 (lisp::double-float-high-bits 1d0))
346 (random-chunk state))
347 1d0)))
348
349 ;;; Using a faster inline VOP.
350 #+x86
351 (defun %random-double-float (arg state)
352 (declare (type (double-float (0d0)) arg)
353 (type random-state state))
354 (let ((state-vector (random-state-state state)))
355 (* arg
356 (- (lisp::make-double-float
357 (dpb (ash (vm::random-mt19937 state-vector)
358 (- vm:double-float-digits random-chunk-length
359 vm:word-bits))
360 vm:double-float-significand-byte
361 (lisp::double-float-high-bits 1d0))
362 (vm::random-mt19937 state-vector))
363 1d0))))
364
365 #+long-float
366 (declaim (inline %random-long-float))
367 #+long-float
368 (declaim (ftype (function ((long-float (0l0)) random-state) (long-float 0l0))
369 %random-long-float))
370
371 ;;; Using a faster inline VOP.
372 #+(and long-float x86)
373 (defun %random-long-float (arg state)
374 (declare (type (long-float (0l0)) arg)
375 (type random-state state))
376 (let ((state-vector (random-state-state state)))
377 (* arg
378 (- (lisp::make-long-float
379 (lisp::long-float-exp-bits 1l0)
380 (logior (vm::random-mt19937 state-vector) vm:long-float-hidden-bit)
381 (vm::random-mt19937 state-vector))
382 1l0))))
383
384 #+(and long-float sparc)
385 (defun %random-long-float (arg state)
386 (declare (type (long-float (0l0)) arg)
387 (type random-state state))
388 (* arg
389 (- (lisp::make-long-float
390 (lisp::long-float-exp-bits 1l0) ; X needs more work
391 (random-chunk state) (random-chunk state) (random-chunk state))
392 1l0)))
393 #+double-double
394 (defun %random-double-double-float (arg state)
395 (declare (type (double-double-float (0w0)) arg)
396 (type random-state state))
397 ;; Generate a 31-bit integer, scale it and sum them up
398 (let* ((r 0w0)
399 (scale (scale-float 1d0 -31))
400 (mult scale))
401 (declare (double-float mult)
402 (type double-double-float r)
403 (optimize (speed 3) (inhibit-warnings 3)))
404 (dotimes (k 4)
405 (setf r (+ r (* mult (ldb (byte 31 0) (random-chunk state)))))
406 (setf mult (* mult scale)))
407 (* arg r)))
408
409
410 ;;;; Random integers:
411
412 ;;; Amount we overlap chunks by when building a large integer to make up for
413 ;;; the loss of randomness in the low bits.
414 ;;;
415 (defconstant random-integer-overlap 3)
416
417 ;;; Extra bits of randomness that we generate before taking the value MOD the
418 ;;; limit, to avoid loss of randomness near the limit.
419 ;;;
420 (defconstant random-integer-extra-bits 10)
421
422 ;;; Largest fixnum we can compute from one chunk of bits.
423 ;;;
424 (defconstant random-fixnum-max
425 (1- (ash 1 (- random-chunk-length random-integer-extra-bits))))
426
427
428 ;;; %RANDOM-INTEGER -- Internal
429 ;;;
430 (defun %random-integer (arg state)
431 (declare (type (integer 1) arg) (type random-state state))
432 (let ((shift (- random-chunk-length random-integer-overlap)))
433 (do ((bits (random-chunk state)
434 (logxor (ash bits shift) (random-chunk state)))
435 (count (+ (integer-length arg)
436 (- random-integer-extra-bits shift))
437 (- count shift)))
438 ((minusp count)
439 (rem bits arg))
440 (declare (fixnum count)))))
441
442 (defun random (arg &optional (state *random-state*))
443 "Generate a uniformly distributed pseudo-random number between zero
444 and Arg. State, if supplied, is the random state to use."
445 (declare (inline %random-single-float %random-double-float
446 #+long-float %long-float))
447 (cond
448 ((typep arg '(integer 1 #x100000000))
449 ;; Let the compiler deftransform take care of this case.
450 (random arg state))
451 ((and (typep arg 'single-float) (> arg 0.0F0))
452 (%random-single-float arg state))
453 ((and (typep arg 'double-float) (> arg 0.0D0))
454 (%random-double-float arg state))
455 #+long-float
456 ((and (typep arg 'long-float) (> arg 0.0L0))
457 (%random-long-float arg state))
458 #+double-double
459 ((and (typep arg 'double-double-float) (> arg 0.0w0))
460 (%random-double-double-float arg state))
461 ((and (integerp arg) (> arg 0))
462 (%random-integer arg state))
463 (t
464 (error 'simple-type-error
465 :expected-type '(or (integer 1) (float (0.0))) :datum arg
466 :format-control (intl:gettext "Argument is not a positive integer or a positive float: ~S")
467 :format-arguments (list arg)))))

  ViewVC Help
Powered by ViewVC 1.1.5