Silence notes about being specialised EQ templates on x86oids
[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 ;;;;   Vol. 8, No. 1, January pp.3-30 (1998) DOI:10.1145/272991.272995
7 ;;;; http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
8
9 ;;;; This software is part of the SBCL system. See the README file for
10 ;;;; more information.
11 ;;;;
12 ;;;; This software is derived from the CMU CL system, which was
13 ;;;; written at Carnegie Mellon University and released into the
14 ;;;; public domain. The software is in the public domain and is
15 ;;;; provided with absolutely no warranty. See the COPYING and CREDITS
16 ;;;; files for more information.
17
18 (in-package "SB!KERNEL")
19 \f
20 ;;;; Constants
21 (defconstant mt19937-n 624)
22 (defconstant mt19937-m 397)
23 (defconstant mt19937-upper-mask #x80000000)
24 (defconstant mt19937-lower-mask #x7FFFFFFF)
25 (defconstant mt19937-a #x9908B0DF)
26 (defconstant mt19937-b #x9D2C5680)
27 (defconstant mt19937-c #xEFC60000)
28
29 ;;;; RANDOM-STATEs
30
31 ;;; The state is stored in a (simple-array (unsigned-byte 32) (627))
32 ;;; wrapped in a random-state structure:
33 ;;;
34 ;;;  0-1:   Constant matrix A. [0, #x9908b0df]
35 ;;;  2:     Index k.
36 ;;;  3-626: State.
37
38 (deftype random-state-state () `(simple-array (unsigned-byte 32) (,(+ 3 mt19937-n))))
39
40 (def!method make-load-form ((random-state random-state) &optional environment)
41   (make-load-form-saving-slots random-state :environment environment))
42
43 (def!method print-object ((state random-state) stream)
44   (if (and *print-readably* (not *read-eval*))
45       (print-not-readable-error state stream)
46       (format stream "#S(~S ~S #.~S)"
47               'random-state
48               ':state
49               `(make-array ,(+ 3 mt19937-n)
50                 :element-type
51                 '(unsigned-byte 32)
52                 :initial-contents
53                 ',(coerce (random-state-state state) 'list)))))
54
55 ;;; Generate and initialize a new random-state array. Index is
56 ;;; initialized to 1 and the states to 32bit integers excluding zero.
57 ;;;
58 ;;; Seed - A 32bit number.
59 ;;;
60 ;;; See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
61 ;;; In the previous versions, MSBs of the seed affect only MSBs of the array.
62 (defun init-random-state (&optional (seed 5489) state)
63   (declare (type (unsigned-byte 32) seed))
64   (let ((state (or state (make-array 627 :element-type '(unsigned-byte 32)))))
65     (check-type state random-state-state)
66     (setf (aref state 0) 0)
67     (setf (aref state 1) mt19937-a)
68     (setf (aref state 2) mt19937-n)
69     (loop for i below mt19937-n
70       for p from 3
71       for s = seed then
72       (logand #xFFFFFFFF
73               (+ (* 1812433253
74                     (logxor s (ash s -30)))
75                  i))
76       do (setf (aref state p) s))
77     state))
78
79 (defvar *random-state*)
80 (defun !random-cold-init ()
81   (/show0 "entering !RANDOM-COLD-INIT")
82   (setf *random-state* (%make-random-state))
83   (/show0 "returning from !RANDOM-COLD-INIT"))
84
85 (defun make-random-state (&optional state)
86   #!+sb-doc
87   "Make a random state object. The optional STATE argument specifies a seed
88 for deterministic pseudo-random number generation.
89
90 As per the Common Lisp standard,
91 - If STATE is NIL or not supplied, return a copy of the default
92   *RANDOM-STATE*.
93 - If STATE is a random state, return a copy of it.
94 - If STATE is T, return a randomly initialized state (using operating-system
95   provided randomness where available, otherwise a poor substitute based on
96   internal time and pid).
97
98 See SB-EXT:SEED-RANDOM-STATE for a SBCL extension to this functionality."
99   (/show0 "entering MAKE-RANDOM-STATE")
100   (check-type state (or boolean random-state))
101   (seed-random-state state))
102
103 (defun fallback-random-seed ()
104   ;; When /dev/urandom is not available, we make do with time and pid
105   ;; Thread ID and/or address of a CONS cell would be even better, but...
106   (/show0 "No /dev/urandom, using randomness from time and pid")
107   (+ (get-internal-real-time)
108      (ash (sb!unix:unix-getpid) 32)))
109
110 #!-win32
111 (defun os-random-seed ()
112   (or
113    ;; On unices, we try to read from /dev/urandom and pass the results
114    ;; to our (simple-array (unsigned-byte 32) (*)) processor below.
115    ;; More than 256 bits would provide a false sense of security.
116    ;; If you need more bits than that, you probably also need
117    ;; a better algorithm too.
118    (ignore-errors
119     (with-open-file (r "/dev/urandom" :element-type '(unsigned-byte 32)
120                                       :direction :input :if-does-not-exist :error)
121       (let ((a (make-array '(8) :element-type '(unsigned-byte 32))))
122         (assert (= 8 (read-sequence a r)))
123         a)))
124    (fallback-random-seed)))
125
126 #!+win32
127 (defun os-random-seed ()
128   (/show0 "Getting randomness from CryptGenRandom")
129   (or (sb!win32:crypt-gen-random 32)
130       (fallback-random-seed)))
131
132 (defun seed-random-state (&optional state)
133   #!+sb-doc
134   "Make a random state object. The optional STATE argument specifies a seed
135 for deterministic pseudo-random number generation.
136
137 As per the Common Lisp standard for MAKE-RANDOM-STATE,
138 - If STATE is NIL or not supplied, return a copy of the default
139   *RANDOM-STATE*.
140 - If STATE is a random state, return a copy of it.
141 - If STATE is T, return a randomly initialized state (using operating-system
142   provided randomness where available, otherwise a poor substitute based on
143   internal time and pid).
144
145 As a supported SBCL extension, we also support receiving as a seed an object
146 of the following types:
147 - (SIMPLE-ARRAY (UNSIGNED-BYTE 8) (*))
148 - UNSIGNED-BYTE
149 While we support arguments of any size and will mix the provided bits into
150 the random state, it is probably overkill to provide more than 256 bits worth
151 of actual information.
152
153 This particular SBCL version also accepts an argument of the following type:
154 (SIMPLE-ARRAY (UNSIGNED-BYTE 32) (*))
155
156 This particular SBCL version uses the popular MT19937 PRNG algorithm, and its
157 internal state only effectively contains about 19937 bits of information.
158 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
159 "
160   (etypecase state
161     ;; Easy standard cases
162     (null
163      (/show0 "copying *RANDOM-STATE*")
164      (%make-random-state :state (copy-seq (random-state-state *random-state*))))
165     (random-state
166      (/show0 "copying the provided RANDOM-STATE")
167      (%make-random-state :state (copy-seq (random-state-state state))))
168     ;; Standard case, less easy: try to randomly initialize a state.
169     ((eql t)
170      (/show0 "getting randomness from the operating system")
171      (seed-random-state (os-random-seed)))
172     ;; For convenience to users, we accept (simple-array (unsigned-byte 8) (*))
173     ;; We just convert it to (simple-array (unsigned-byte 32) (*)) in a
174     ;; completely straightforward way.
175     ;; TODO: probably similarly accept other word sizes.
176     ((simple-array (unsigned-byte 8) (*))
177      (/show0 "getting random seed from byte vector (converting to 32-bit-word vector)")
178      (let* ((l (length state))
179             (m (ceiling l 4))
180             (r (if (>= l 2496) 0 (mod l 4)))
181             (y (make-array (list m) :element-type '(unsigned-byte 32))))
182              (loop for i from 0 below (- m (if (zerop r) 0 1))
183                for j = (* i 4) do
184                (setf (aref y i)
185                      (+ (aref state j)
186                         (ash (aref state (+ j 1)) 8)
187                         (ash (aref state (+ j 2)) 16)
188                         (ash (aref state (+ j 3)) 24))))
189              (unless (zerop r) ;; The last word may require special treatment.
190                (let* ((p (1- m)) (q (* 4 p)))
191                  (setf (aref y p)
192                      (+ (aref state q)
193                         (if (< 1 r) (ash (aref state (+ q 1)) 8) 0)
194                         (if (= 3 r) (ash (aref state (+ q 2)) 16) 0)))))
195              (seed-random-state y)))
196     ;; Also for convenience, we accept non-negative integers as seeds.
197     ;; Small ones get passed to init-random-state, as before.
198     ((unsigned-byte 32)
199      (/show0 "getting random seed from 32-bit word")
200      (%make-random-state :state (init-random-state state)))
201     ;; Larger ones ones get trivially chopped into an array of (unsigned-byte 32)
202     ((unsigned-byte)
203      (/show0 "getting random seed from bignum (converting to 32-bit-word vector)")
204      (loop with l = (ceiling (integer-length state) 32)
205        with s = (make-array (list l) :element-type '(unsigned-byte 32))
206        for i below l
207        for p from 0 by 32
208        do (setf (aref s i) (ldb (byte 32 p) state))
209        finally (return (seed-random-state s))))
210     ;; Last but not least, when provided an array of 32-bit words, we truncate
211     ;; it to 19968 bits and mix these into an initial state. We reuse the same
212     ;; method as the authors of the original algorithm. See
213     ;; http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c
214     ;; NB: their mt[i] is our (aref s (+ 3 i))
215     ((simple-array (unsigned-byte 32) (*))
216      (/show0 "getting random seed from 32-bit-word vector")
217      (let ((s (init-random-state 19650218))
218            (i 1) (j 0) (l (length state)))
219        (loop for k downfrom (max mt19937-n l) above 0 do
220          (setf (aref s (+ i 3))
221                (logand #xFFFFFFFF
222                        (+ (logxor (aref s (+ i 3))
223                                   (* 1664525
224                                      (logxor (aref s (+ i 2))
225                                              (ash (aref s (+ i 2)) -30))))
226                           (aref state j) j))) ;; non-linear
227          (incf i) (when (>= i mt19937-n) (setf (aref s 3) (aref s (+ 2 mt19937-n)) i 1))
228          (incf j) (when (>= j l) (setf j 0)))
229        (loop for k downfrom (1- mt19937-n) above 0 do
230          (setf (aref s (+ i 3))
231                (logand #xFFFFFFFF
232                        (- (logxor (aref s (+ i 3))
233                                   (* 1566083941
234                                      (logxor (aref s (+ i 2))
235                                              (ash (aref s (+ i 2)) -30))))
236                           i))) ;; non-linear
237          (incf i) (when (>= i mt19937-n) (setf (aref s 3) (aref s (+ 2 mt19937-n)) i 1)))
238        (setf (aref s 3) #x80000000) ;; MSB is 1; assuring non-zero initial array
239        (%make-random-state :state s)))))
240 \f
241 ;;;; random entries
242
243 ;;; This function generates a 32bit integer between 0 and #xffffffff
244 ;;; inclusive.
245 #!-sb-fluid (declaim (inline random-chunk))
246 ;;; portable implementation
247 #!-x86
248 (defun random-mt19937-update (state)
249   (declare (type random-state-state state)
250            (optimize (speed 3) (safety 0)))
251   (let ((y 0))
252     (declare (type (unsigned-byte 32) y))
253     (do ((kk 3 (1+ kk)))
254         ((>= kk (+ 3 (- mt19937-n mt19937-m))))
255       (declare (type (mod 628) kk))
256       (setf y (logior (logand (aref state kk) mt19937-upper-mask)
257                       (logand (aref state (1+ kk)) mt19937-lower-mask)))
258       (setf (aref state kk) (logxor (aref state (+ kk mt19937-m))
259                                     (ash y -1) (aref state (logand y 1)))))
260     (do ((kk (+ (- mt19937-n mt19937-m) 3) (1+ kk)))
261         ((>= kk (+ (1- mt19937-n) 3)))
262       (declare (type (mod 628) kk))
263       (setf y (logior (logand (aref state kk) mt19937-upper-mask)
264                       (logand (aref state (1+ kk)) mt19937-lower-mask)))
265       (setf (aref state kk) (logxor (aref state (+ kk (- mt19937-m mt19937-n)))
266                                     (ash y -1) (aref state (logand y 1)))))
267     (setf y (logior (logand (aref state (+ 3 (1- mt19937-n)))
268                             mt19937-upper-mask)
269                     (logand (aref state 3) mt19937-lower-mask)))
270     (setf (aref state (+ 3 (1- mt19937-n)))
271           (logxor (aref state (+ 3 (1- mt19937-m)))
272                   (ash y -1) (aref state (logand y 1)))))
273   (values))
274 #!-x86
275 (defun random-chunk (state)
276   (declare (type random-state state))
277   (let* ((state (random-state-state state))
278          (k (aref state 2)))
279     (declare (type (mod 628) k))
280     (when (= k mt19937-n)
281       (random-mt19937-update state)
282       (setf k 0))
283     (setf (aref state 2) (1+ k))
284     (let ((y (aref state (+ 3 k))))
285       (declare (type (unsigned-byte 32) y))
286       (setf y (logxor y (ash y -11)))
287       (setf y (logxor y (ash (logand y (ash mt19937-b -7)) 7)))
288       (setf y (logxor y (ash (logand y (ash mt19937-c -15)) 15)))
289       (setf y (logxor y (ash y -18)))
290       y)))
291
292 ;;; Using inline VOP support, only available on the x86 so far.
293 ;;;
294 ;;; FIXME: It would be nice to have some benchmark numbers on this.
295 ;;; My inclination is to get rid of the nonportable implementation
296 ;;; unless the performance difference is just enormous.
297 #!+x86
298 (defun random-chunk (state)
299   (declare (type random-state state))
300   (sb!vm::random-mt19937 (random-state-state state)))
301
302 #!-sb-fluid (declaim (inline big-random-chunk))
303 (defun big-random-chunk (state)
304   (declare (type random-state state))
305   (logior (ash (random-chunk state) 32)
306           (random-chunk state)))
307 \f
308 ;;; Handle the single or double float case of RANDOM. We generate a
309 ;;; float between 0.0 and 1.0 by clobbering the significand of 1.0
310 ;;; with random bits, then subtracting 1.0. This hides the fact that
311 ;;; we have a hidden bit.
312 #!-sb-fluid (declaim (inline %random-single-float %random-double-float))
313 (declaim (ftype (function ((single-float (0f0)) random-state)
314                           (single-float 0f0))
315                 %random-single-float))
316 (defun %random-single-float (arg state)
317   (declare (type (single-float (0f0)) arg)
318            (type random-state state))
319   (* arg
320      (- (make-single-float
321          (dpb (ash (random-chunk state)
322                    (- sb!vm:single-float-digits n-random-chunk-bits))
323               sb!vm:single-float-significand-byte
324               (single-float-bits 1.0)))
325         1.0)))
326 (declaim (ftype (function ((double-float (0d0)) random-state)
327                           (double-float 0d0))
328                 %random-double-float))
329
330 ;;; 32-bit version
331 #!+nil
332 (defun %random-double-float (arg state)
333   (declare (type (double-float (0d0)) arg)
334            (type random-state state))
335   (* (float (random-chunk state) 1d0) (/ 1d0 (expt 2 32))))
336
337 ;;; 53-bit version
338 #!-x86
339 (defun %random-double-float (arg state)
340   (declare (type (double-float (0d0)) arg)
341            (type random-state state))
342   (* arg
343      (- (sb!impl::make-double-float
344          (dpb (ash (random-chunk state)
345                    (- sb!vm:double-float-digits n-random-chunk-bits 32))
346               sb!vm:double-float-significand-byte
347               (sb!impl::double-float-high-bits 1d0))
348          (random-chunk state))
349         1d0)))
350
351 ;;; using a faster inline VOP
352 #!+x86
353 (defun %random-double-float (arg state)
354   (declare (type (double-float (0d0)) arg)
355            (type random-state state))
356   (let ((state-vector (random-state-state state)))
357     (* arg
358        (- (sb!impl::make-double-float
359            (dpb (ash (sb!vm::random-mt19937 state-vector)
360                      (- sb!vm:double-float-digits n-random-chunk-bits
361                         sb!vm:n-word-bits))
362                 sb!vm:double-float-significand-byte
363                 (sb!impl::double-float-high-bits 1d0))
364            (sb!vm::random-mt19937 state-vector))
365           1d0))))
366
367 \f
368 ;;;; random fixnums
369
370 ;;; Generate and return a pseudo random fixnum less than ARG. To achieve
371 ;;; equidistribution an accept-reject loop is used.
372 ;;; No extra effort is made to detect the case of ARG being a power of
373 ;;; two where rejection is not possible, as the cost of checking for
374 ;;; this case is the same as doing the rejection test. When ARG is
375 ;;; larger than (expt 2 N-RANDOM-CHUNK-BITS), which can only happen if
376 ;;; the random chunk size is half the word size, two random chunks are
377 ;;; used in each loop iteration, otherwise only one. Finally, the
378 ;;; rejection probability could often be reduced by not masking the
379 ;;; chunk but rejecting only values as least as large as the largest
380 ;;; multiple of ARG that fits in a chunk (or two), but this is not done
381 ;;; as the speed gains due to needing fewer loop iterations are by far
382 ;;; outweighted by the cost of the two divisions required (one to find
383 ;;; the multiplier and one to bring the result into the correct range).
384 #!-sb-fluid (declaim (inline %random-fixnum))
385 (defun %random-fixnum (arg state)
386   (declare (type (integer 1 #.sb!xc:most-positive-fixnum) arg)
387            (type random-state state))
388   (if (= arg 1)
389       0
390       (let* ((n-bits (integer-length (1- arg)))
391              (mask (1- (ash 1 n-bits))))
392         (macrolet ((accept-reject-loop (generator)
393                      `(loop
394                         (let ((bits (logand mask (,generator state))))
395                           (when (< bits arg)
396                             (return bits))))))
397           (assert (<= n-bits (* 2 n-random-chunk-bits)))
398           (if (<= n-bits n-random-chunk-bits)
399               (accept-reject-loop random-chunk)
400               (accept-reject-loop big-random-chunk))))))
401
402 (defun random (arg &optional (state *random-state*))
403   (declare (inline %random-fixnum %random-single-float %random-double-float
404                    #!+long-float %random-long-float))
405   (cond
406     ((and (fixnump arg) (> arg 0))
407      (%random-fixnum arg state))
408     ((and (typep arg 'single-float) (> arg 0.0f0))
409      (%random-single-float arg state))
410     ((and (typep arg 'double-float) (> arg 0.0d0))
411      (%random-double-float arg state))
412     #!+long-float
413     ((and (typep arg 'long-float) (> arg 0.0l0))
414      (%random-long-float arg state))
415     ((and (bignump arg) (> arg 0))
416      (%random-bignum arg state))
417     (t
418      (error 'simple-type-error
419             :expected-type '(or (integer 1) (float (0))) :datum arg
420             :format-control "~@<Argument is neither a positive integer nor a ~
421                              positive float: ~2I~_~S~:>"
422             :format-arguments (list arg)))))