0.8.15.7
[sbcl.git] / src / code / target-random.lisp
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,
6 ;;;;   1997, to appear.
7
8 ;;;; This software is part of the SBCL system. See the README file for
9 ;;;; more information.
10 ;;;;
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.
16
17 (in-package "SB!KERNEL")
18 \f
19 ;;;; RANDOM-STATEs
20
21 (def!method make-load-form ((random-state random-state) &optional environment) 
22   (make-load-form-saving-slots random-state :environment environment))
23
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)"
28               'random-state
29               ':state
30               `(make-array 627 
31                 :element-type 
32                 '(unsigned-byte 32)
33                 :initial-contents 
34                 ',(coerce (random-state-state state) 'list)))))
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 and initialize a new random-state array. Index is
44 ;;; initialized to 1 and the states to 32bit integers excluding zero.
45 ;;;
46 ;;; Seed - A 32bit number, not zero.
47 ;;;
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)
57     (do ((k 1 (1+ k)))
58         ((>= k 624))
59       (declare (type (mod 625) k))
60       (setf (aref state (+ 3 k))
61             (logand (* 69069 (aref state (+ 3 (1- k)))) #xffffffff)))
62     state))
63
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"))
69
70 (defun make-random-state (&optional state)
71   #!+sb-doc
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
75   the universal time."
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))
80                  (new-state
81                   (make-array 627 :element-type '(unsigned-byte 32))))
82              (/show0 "made NEW-STATE, about to DOTIMES")
83              (dotimes (i 627)
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")
88     (etypecase state
89       (null
90        (/show0 "NULL case")
91        (copy-random-state *random-state*))
92       (random-state
93        (/show0 "RANDOM-STATE-P clause")
94        (copy-random-state state))
95       ((member t)
96        (/show0 "T clause")
97        (%make-random-state :state (init-random-state
98                                    (logand (get-universal-time)
99                                            #xffffffff)))))))
100 \f
101 ;;;; random entries
102
103 ;;; This function generates a 32bit integer between 0 and #xffffffff
104 ;;; inclusive.
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)
113 #!-x86
114 (defun random-mt19937-update (state)
115   (declare (type (simple-array (unsigned-byte 32) (627)) state)
116            (optimize (speed 3) (safety 0)))
117   (let ((y 0))
118     (declare (type (unsigned-byte 32) y))
119     (do ((kk 3 (1+ kk)))
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)))
134                             mt19937-upper-mask)
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)))))
139   (values))
140 #!-x86
141 (defun random-chunk (state)
142   (declare (type random-state state))
143   (let* ((state (random-state-state state))
144          (k (aref state 2)))
145     (declare (type (mod 628) k))
146     (when (= k mt19937-n)
147       (random-mt19937-update state)
148       (setf k 0))
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)))
156       y)))
157
158 ;;; Using inline VOP support, only available on the x86 so far.
159 ;;;
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.
163 #!+x86
164 (defun random-chunk (state)
165   (declare (type random-state state))
166   (sb!vm::random-mt19937 (random-state-state state)))
167 \f
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)
174                           (single-float 0f0))
175                 %random-single-float))
176 (defun %random-single-float (arg state)
177   (declare (type (single-float (0f0)) arg)
178            (type random-state state))
179   (* arg
180      (- (make-single-float
181          (dpb (ash (random-chunk state)
182                    (- sb!vm:single-float-digits random-chunk-length))
183               sb!vm:single-float-significand-byte
184               (single-float-bits 1.0)))
185         1.0)))
186 (declaim (ftype (function ((double-float (0d0)) random-state)
187                           (double-float 0d0))
188                 %random-double-float))
189
190 ;;; 32-bit version
191 #!+nil
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))))
196
197 ;;; 53-bit version
198 #!-x86
199 (defun %random-double-float (arg state)
200   (declare (type (double-float (0d0)) arg)
201            (type random-state state))
202   (* arg
203      (- (sb!impl::make-double-float
204          (dpb (ash (random-chunk state)
205                    (- sb!vm:double-float-digits random-chunk-length
206                       sb!vm:n-word-bits))
207               sb!vm:double-float-significand-byte
208               (sb!impl::double-float-high-bits 1d0))
209          (random-chunk state))
210         1d0)))
211
212 ;;; using a faster inline VOP
213 #!+x86
214 (defun %random-double-float (arg state)
215   (declare (type (double-float (0d0)) arg)
216            (type random-state state))
217   (let ((state-vector (random-state-state state)))
218     (* arg
219        (- (sb!impl::make-double-float
220            (dpb (ash (sb!vm::random-mt19937 state-vector)
221                      (- sb!vm:double-float-digits random-chunk-length
222                         sb!vm:n-word-bits))
223                 sb!vm:double-float-significand-byte
224                 (sb!impl::double-float-high-bits 1d0))
225            (sb!vm::random-mt19937 state-vector))
226           1d0))))
227
228 \f
229 ;;;; random integers
230
231 (defun %random-integer (arg state)
232   (declare (type (integer 1) arg) (type random-state state))
233   (let ((shift (- random-chunk-length random-integer-overlap)))
234     (do ((bits (random-chunk state)
235                (logxor (ash bits shift) (random-chunk state)))
236          (count (+ (integer-length arg)
237                    (- random-integer-extra-bits shift))
238                 (- count shift)))
239         ((minusp count)
240          (rem bits arg))
241       (declare (fixnum count)))))
242
243 (defun random (arg &optional (state *random-state*))
244   (declare (inline %random-single-float %random-double-float
245                    #!+long-float %random-long-float))
246   (cond
247     ((and (fixnump arg) (<= arg random-fixnum-max) (> arg 0))
248      (rem (random-chunk state) arg))
249     ((and (typep arg 'single-float) (> arg 0.0f0))
250      (%random-single-float arg state))
251     ((and (typep arg 'double-float) (> arg 0.0d0))
252      (%random-double-float arg state))
253     #!+long-float
254     ((and (typep arg 'long-float) (> arg 0.0l0))
255      (%random-long-float arg state))
256     ((and (integerp arg) (> arg 0))
257      (%random-integer arg state))
258     (t
259      (error 'simple-type-error
260             :expected-type '(or (integer 1) (float (0))) :datum arg
261             :format-control "~@<Argument is neither a positive integer nor a ~
262                              positive float: ~2I~_~S~:>"
263             :format-arguments (list arg)))))