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

Contents of /src/code/rand.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.11.54.1 - (show annotations)
Mon Feb 8 17:15:48 2010 UTC (4 years, 2 months ago) by rtoy
Branch: intl-branch
CVS Tags: intl-branch-working-2010-02-19-1000, intl-branch-working-2010-02-11-1000, intl-branch-2010-03-18-1300
Changes since 1.11: +3 -1 lines
Add (intl:textdomain "cmucl") to the files to set the textdomain.
1 ;;; -*- Mode: Lisp; Package: Kernel -*-
2 ;;;
3 ;;; **********************************************************************
4 ;;; This code was written as part of the CMU Common Lisp project at
5 ;;; Carnegie Mellon University, and has been placed in the public domain.
6 ;;;
7 (ext:file-comment
8 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/rand.lisp,v 1.11.54.1 2010/02/08 17:15:48 rtoy Exp $")
9 ;;;
10 ;;; **********************************************************************
11 ;;;
12 ;;; Functions to random number functions for Spice Lisp
13 ;;;
14 ;;; Originally written by David Adam. Python tuning, better large integer
15 ;;; randomness and efficient IEEE float support by Rob MacLachlan.
16 ;;;
17 #-new-random
18 (in-package "LISP")
19 (intl:textdomain "cmucl")
20
21 #-new-random
22 (export '(random-state random-state-p random *random-state*
23 make-random-state))
24
25 #-new-random
26 (in-package "KERNEL")
27 #-new-random
28 (export '(%random-single-float %random-double-float random-chunk
29 random-fixnum-max))
30
31
32 ;;;; Random state hackery:
33
34 #-new-random
35 (progn
36 (defconstant random-const-a 8373)
37 (defconstant random-const-c 101010101)
38 (defconstant random-max 54)
39
40 ;;; Inclusive upper bound on the size of fixnum kept in the state (and returned
41 ;;; by random-chunk.) Must be even.
42 ;;;
43 (defconstant random-upper-bound (- most-positive-fixnum 3))
44 (defconstant random-chunk-length (integer-length random-upper-bound))
45 (deftype random-chunk () `(integer 0 ,random-upper-bound))
46
47 (defvar rand-seed 0)
48 (defstruct (random-state
49 (:constructor make-random-object)
50 (:make-load-form-fun :just-dump-it-normally))
51 (j 24 :type index)
52 (k 0 :type index)
53 (seed (make-array (1+ random-max) :initial-contents
54 (do ((list-rands () (cons (rand1) list-rands))
55 (i 0 (1+ i)))
56 ((> i random-max) list-rands)
57 (declare (fixnum i))))
58 :type simple-vector))
59
60
61 ;;; Generates a random number from rand-seed.
62 (defun rand1 ()
63 (declare (optimize (inhibit-warnings 3)))
64 (setq rand-seed
65 (mod (+ (* rand-seed random-const-a) random-const-c)
66 (1+ random-upper-bound))))
67
68
69 (defvar *random-state* (make-random-object))
70
71
72 (defun copy-state (cur-state)
73 (let ((state (make-random-object
74 :seed (make-array 55)
75 :j (random-state-j cur-state)
76 :k (random-state-k cur-state))))
77 (do ((i 0 (1+ i)))
78 ((= i 55) state)
79 (declare (fixnum i))
80 (setf (aref (random-state-seed state) i)
81 (aref (random-state-seed cur-state) i)))))
82
83 (defun make-random-state (&optional state)
84 "Make a random state object. If State is not supplied, return a copy
85 of the default random state. If State is a random state, then return a
86 copy of it. If state is T then return a random state generated from
87 the universal time."
88 (cond ((not state) (copy-state *random-state*))
89 ((random-state-p state) (copy-state state))
90 ((eq state t) (setq rand-seed (get-universal-time))
91 (make-random-object))
92 (t (error 'simple-type-error
93 :expected-type '(or null t (satisfies random-state-p))
94 :datum state
95 :format-control "Argument is not a RANDOM-STATE, T or NIL: ~S"
96 :format-arguments (list state)))))
97
98
99 ;;;; Random entries:
100
101 (declaim (start-block random %random-single-float %random-double-float
102 random-chunk))
103
104 ;;; random-chunk -- Internal
105 ;;;
106 ;;; This function generates fixnums between 0 and random-upper-bound,
107 ;;; inclusive. For the algorithm to work random-upper-bound must be an
108 ;;; even positive fixnum. State is the random state to use.
109 ;;;
110 (declaim (ftype (function (random-state) random-chunk) random-chunk))
111 (defun random-chunk (state)
112 (let* ((seed (random-state-seed state))
113 (j (random-state-j state))
114 (k (random-state-k state))
115 (a (- (- random-upper-bound
116 (the random-chunk
117 (svref seed
118 (setf (random-state-j state)
119 (if (= j 0) random-max (1- j))))))
120 (the random-chunk
121 (svref seed
122 (setf (random-state-k state)
123 (if (= k 0) random-max (1- k))))))))
124 (declare (fixnum a))
125 (setf (svref seed k)
126 (the random-chunk (if (minusp a) (- a) (- random-upper-bound a))))))
127
128
129 ;;; %RANDOM-SINGLE-FLOAT, %RANDOM-DOUBLE-FLOAT -- Interface
130 ;;;
131 ;;; Handle the single or double float case of RANDOM. We generate a float
132 ;;; between 0.0 and 1.0 by clobbering the significand of 1.0 with random bits,
133 ;;; then subtracting 1.0. This hides the fact that we have a hidden bit.
134 ;;;
135 (declaim (inline %random-single-float %random-double-float))
136 (defun %random-single-float (arg state)
137 (declare (type (single-float (0f0)) arg)
138 (type (or random-state null) state))
139 (* arg
140 (- (make-single-float
141 (dpb (ash (random-chunk (or state *random-state*))
142 (- vm:single-float-digits random-chunk-length))
143 vm:single-float-significand-byte
144 (single-float-bits 1.0)))
145 1.0)))
146 ;;;
147 (defun %random-double-float (arg state)
148 (declare (type (double-float (0d0)) arg)
149 (type (or random-state null) state))
150 (let ((state (or state *random-state*)))
151 (* arg
152 (- (make-double-float
153 (dpb (ash (random-chunk state)
154 (- vm:double-float-digits random-chunk-length
155 vm:word-bits))
156 vm:double-float-significand-byte
157 (double-float-high-bits 1d0))
158 (logxor (ash (random-chunk state)
159 (- vm:word-bits random-chunk-length))
160 (ash (random-chunk state)
161 (- random-chunk-length vm:word-bits))))
162 1d0))))
163
164
165 ;;;; Random integers:
166
167 ;;; Amount we overlap chunks by when building a large integer to make up for
168 ;;; the loss of randomness in the low bits.
169 ;;;
170 (defconstant random-integer-overlap 3)
171
172 ;;; Extra bits of randomness that we generate before taking the value MOD the
173 ;;; limit, to avoid loss of randomness near the limit.
174 ;;;
175 (defconstant random-integer-extra-bits 10)
176
177 ;;; Largest fixnum we can compute from one chunk of bits.
178 ;;;
179 (defconstant random-fixnum-max
180 (1- (ash 1 (- random-chunk-length random-integer-extra-bits))))
181
182
183 ;;; %RANDOM-INTEGER -- Internal
184 ;;;
185 (defun %random-integer (arg state)
186 (declare (type (integer 1) arg) (type random-state state))
187 (let ((shift (- random-chunk-length random-integer-overlap)))
188 (do ((bits (random-chunk state)
189 (logxor (ash bits shift) (random-chunk state)))
190 (count (+ (integer-length arg)
191 (- random-integer-extra-bits shift))
192 (- count shift)))
193 ((minusp count)
194 (rem bits arg))
195 (declare (fixnum count)))))
196
197 (defun random (arg &optional (state *random-state*))
198 "Generate a uniformly distributed pseudo-random number between zero
199 and Arg. State, if supplied, is the random state to use."
200 (declare (inline %random-single-float %random-double-float))
201 (cond
202 ((and (fixnump arg) (<= arg random-fixnum-max))
203 (rem (random-chunk state) arg))
204 ((typep arg 'single-float)
205 (%random-single-float arg state))
206 ((typep arg 'double-float)
207 (%random-double-float arg state))
208 ((integerp arg)
209 (%random-integer arg state))
210 (t
211 (error 'simple-type-error :expected-type '(real (0)) :datum arg
212 :format-control "Argument is not a positive real number: ~S"
213 :format-arguments (list arg)))))
214
215 ) ; end progn
216
217 ;;;
218 ;;; **********************************************************************
219 ;;;
220 ;;; Functions for generating random numbers for CMU Common Lisp.
221 ;;;
222 ;;; Written by Raymond Toy. This is a replacement of the original
223 ;;; rand.lisp in CMUCL. According to tests done by Peter VanEynde,
224 ;;; the original generator fails some tests of randomness in
225 ;;; Marsaglia's diehard test suite. This generator passes those
226 ;;; tests.
227 ;;;
228 ;;; Also, the original generator lacks documentation. This remedies
229 ;;; that, and includes references for the algorithms used. Hopefully
230 ;;; better algorithms have been chosen that are at least as good and
231 ;;; at least as fast as the original. In fact, the heart of the
232 ;;; generator is geared towards generating floating point numbers
233 ;;; which should be faster than the original. This should also give
234 ;;; faster integer numbers since longer floating-point random numbers
235 ;;; are generated. Tests have shown that this new generator is always
236 ;;; a 10-20% faster and upto 3 times faster for double-floats.
237 ;;;
238
239 #+new-random
240 (in-package "LISP")
241 #+new-random
242 (export '(random-state random-state-p random *random-state*
243 make-random-state))
244
245 #+new-random
246 (in-package "KERNEL")
247 #+new-random
248 (export '(%random-single-float %random-double-float))
249
250
251 #+new-random
252 (progn
253 ;;;; Random state hackery used to initialize the random state
254 (defvar rand-seed 54217137)
255 (declaim (type (integer 0 900000000) rand-seed))
256
257 ;;; This is the function RMARIN, in James' paper. The routine is
258 ;;; slightly extended to handle the longer floats that we are using
259 ;;; here.
260
261 (defun init-random-state ()
262 (flet ((init-aux (y0 y1 y2 z0)
263 (declare (type (mod 179) y0 y1 y2)
264 (type (mod 169) z0))
265 (let ((r (make-array 32 :element-type 'double-float)))
266 (dotimes (n 32
267 r)
268 (declare (fixnum n))
269 (let ((s 0d0)
270 (tt 0.5d0))
271 (declare (double-float s tt))
272 (dotimes (m #.(float-digits 1d0))
273 (declare (fixnum m))
274 (let ((new-y (mod (* y0 y1 y2) 179))
275 (new-z (mod (+ 1 (* 53 z0)) 169)))
276 (when (> (mod (* new-y new-z) 63) 32)
277 (incf s tt))
278 (setf tt (* tt 0.5d0))
279 (setf y0 y1)
280 (setf y1 y2)
281 (setf y2 new-y)
282 (setf z0 new-z)))
283 (setf (aref r n) s))))))
284
285 (multiple-value-bind (ij kl)
286 (floor rand-seed 30082)
287 (declare (type (integer 0 29918) ij)
288 (type (integer 0 30081) kl))
289 (multiple-value-bind (i j)
290 (floor ij 177)
291 (declare (type (mod 179) i j))
292 (incf i 2)
293 (incf j 2)
294 (multiple-value-bind (k L)
295 (floor kl 169)
296 (declare (fixnum k L))
297 (incf k 1)
298 (init-aux i j k l))))))
299
300 (defstruct (random-state
301 (:constructor make-random-object)
302 (:make-load-form-fun :just-dump-it-normally))
303 (index 0
304 :type (mod 32))
305 (weyl 0.0d0
306 :type double-float)
307 (cache (make-array 32
308 :element-type 'double-float
309 :initial-contents (init-random-state))
310 :type (simple-array double-float (32)))
311 (borrow 0.0d0
312 :type double-float))
313
314 (defvar *random-state* (make-random-object))
315
316 (defun copy-state (cur-state)
317 (let ((state (make-random-object
318 :index (random-state-index cur-state)
319 :weyl (random-state-weyl cur-state)
320 :borrow (random-state-borrow cur-state)
321 :cache (make-array 32 :element-type 'double-float))))
322 (dotimes (k 32
323 state)
324 (declare (fixnum k))
325 (setf (aref (random-state-cache state) k)
326 (aref (random-state-cache cur-state) k)))))
327
328 (defun make-random-state (&optional state)
329 "Make a random state object. If State is not supplied, return a copy
330 of the default random state. If State is a random state, then return a
331 copy of it. If state is T then return a random state generated from
332 the universal time."
333 (cond ((not state)
334 (copy-state *random-state*))
335 ((random-state-p state)
336 (copy-state state))
337 ((eq state t)
338 (setq rand-seed (mod (get-universal-time) 900000000))
339 (make-random-object))
340 (t (error 'simple-type-error
341 :expected-type '(or null t (satisfies random-state-p))
342 :datum state
343 :format-control "Argument is not a RANDOM-STATE, T or NIL: ~S"
344 :format-arguments (list state)))))
345
346
347 ;;;; Random entries
348
349 #+nil ; Not yet as block compiling tickles compiler trouble.
350 (declaim (start-block random %random-single-float %random-double-float
351 new-random-float))
352
353 ;;; random-float -- Internal
354 ;;;
355 ;;; Generates a double-float random number between 0 and 1.
356 ;;;
357 ;;; The technique used is a subtract-with-borrow generator combined
358 ;;; with a Weyl generator.
359 ;;;
360 ;;; The subtract-with-borrow generator is described in the Matlab News
361 ;;; and Notes, Fall 1995, based on suggestions from George
362 ;;; Marsaglia[1]. The generator works as follows.
363 ;;;
364 ;;; The i'th random number is given by
365 ;;;
366 ;;; z[i] = z[i + 20] - z[i + 5] - b
367 ;;;
368 ;;; The indices, i, i + 20, and i + 5 are all interpreted mod 32,
369 ;;; i.e., use just the 5 least significant bits. (Note that we really
370 ;;; only need the last 20 values of z, not the last 32. However,
371 ;;; masking out 5 bits should be faster than comparing against 20.)
372 ;;;
373 ;;; The quantity b is carried over from the previous iteration and is
374 ;;; either 0 or a small quantity. If z[i] is positive, b is set to
375 ;;; zero. If z[i] is negative, we make it positive by adding 1 and we
376 ;;; set b to be 2^-53, or double-float-epsilon for double precision
377 ;;; IEEE floats.
378 ;;;
379 ;;; This generator is claimed to have a period of almost 2^1430
380 ;;; values. However, it has a defect: all numbers are integer
381 ;;; multiples of 2^-53. Consequently many of the floating-point
382 ;;; numbers are not represented.
383 ;;;
384 ;;; The article describes one method for overcoming this difficulty,
385 ;;; but does not give the full details. Therefore, we use the
386 ;;; technique of combining two generators. The basic idea is to
387 ;;; combine the output from the subtract-with-borrow generator above
388 ;;; with the Weyl generator: x[i] = x[i-1] - y mod 1, where y is any
389 ;;; (fixed) irrational number. This technique is used in RANMAR [2]
390 ;;; and ACARRYPC [3].
391 ;;;
392 ;;; In RANMAR, y = 7654321 / 2^24. In ACARRYPC, y = 632436069 / 2^32
393 ;;; is used. We use y = 7097293079245107 / 2^53, where y was randomly
394 ;;; selected.
395 ;;;
396 ;;; Thus, the random number returned is (z[i] - x[i]) mod 1.
397 ;;;
398 ;;; While we could have used a Lisp version of RANMAR, we choose not
399 ;;; to because it only produces single-float results. The generator
400 ;;; given in Matlab News and Notes generates double-float results and
401 ;;; only requires about 32 double-floats of state compared to 97
402 ;;; single-floats for RANMAR. James also gives a subtract-with-borrow
403 ;;; generator, but it is also only for single-floats.
404 ;;;
405 ;;; Marsaglia's ACARRYPC (below), generates 32-bit integers.
406 ;;;
407 ;;; This generator uses Matlab's generator combined with the technique
408 ;;; given by Marsaglia for combining a Weyl generator with the
409 ;;; subtract-with-borrow generator. This idea also seems to be used
410 ;;; in RANMAR.
411 ;;;
412 ;;; References:
413 ;;;
414 ;;; [1] Cleve Moler, Matlab News and Notes, Fall 1995.
415 ;;;
416 ;;; [2] F. James, "A Review of Pseudorandom Number Generators",
417 ;;; Computer Physics Communications 60 (1990).
418 ;;;
419 ;;; [3] George Marsaglia, B. Narasimhan, and Arif Zaman, "A Random
420 ;;; Number Generator for PC's", Computer Physics Communications 60
421 ;;; (1990) 345-349.
422
423 ;;; According to Marsaglia, a suitable constant is any integer
424 ;;; relatively prime to 2^53. We randomly selected 7097293079245107.
425 ;;;
426 (defconstant weyl-constant #.(scale-float (float 7097293079245107 1d0) -53))
427
428 (declaim (ftype (function (random-state)
429 (double-float 0.0d0 (1.0d0)))
430 new-random-float))
431
432 (declaim (inline new-random-float))
433
434 (defun new-random-float (state)
435 (declare (type random-state state))
436 (let* ((z (random-state-cache state))
437 (index (random-state-index state))
438 (new-z (- (aref z (logand (+ index 20) #x1f))
439 (aref z (logand (+ index 5) #x1f))
440 (random-state-borrow state))))
441 (declare (double-float new-z)
442 (type (simple-array double-float (32)) z)
443 (type (mod 32) index))
444 ;; Do the subtract-with-borrow generator
445 (setf (random-state-borrow state)
446 (cond ((minusp new-z)
447 (incf new-z 1.0d0)
448 double-float-epsilon)
449 (t
450 0d0)))
451 (setf (aref z index) new-z)
452 (setf (random-state-index state) (logand (1+ index) #x1f))
453
454 ;; Weyl generator
455 (let ((w (- (random-state-weyl state) weyl-constant)))
456 (declare (double-float w))
457 (let ((w1 (if (minusp w) (1+ w) w)))
458 (declare (double-float w1))
459 (setf (random-state-weyl state) w1)
460
461 ;; Mix them together
462 (decf new-z w1)
463 (when (minusp new-z)
464 (incf new-z 1d0))
465 new-z))))
466
467
468 ;;; %RANDOM-SINGLE-FLOAT, %RANDOM-DOUBLE-FLOAT -- Interface
469 ;;;
470 ;;; Handle the single or double float case of RANDOM. These
471 ;;; functions are required because the compiler will optimize calls
472 ;;; to random to one of the following functions if possible.
473 ;;;
474
475 (declaim (inline %random-single-float %random-double-float))
476 (defun %random-double-float (arg state)
477 (declare (type (double-float (0d0)) arg)
478 (type random-state state))
479 (* arg (new-random-float state)))
480
481 (defun %random-single-float (arg state)
482 (declare (type (single-float (0f0)) arg)
483 (type random-state state))
484 (coerce (%random-double-float (coerce arg 'double-float) state)
485 'single-float))
486
487
488 ;;;; Random integers:
489
490 ;;; Random integers are created by concatenating a bunch of chunks
491 ;;; together. Chunks are generated without consing by truncating a
492 ;;; random float to an integer. With the latest propagate-float-type
493 ;;; code the compiler can inline truncate (signed-byte 32) allowing 31
494 ;;; bits, and (unsigned-byte 32) 32 bits on the x86. When not using the
495 ;;; propagate-float-type feature the best size that can be inlined is
496 ;;; 29 bits. Use of sizes greater than 29 bits causes consing in
497 ;;; argument passing to generic functions, so 29 bits is a good
498 ;;; default.
499 ;;;
500 (defconstant random-chunk-size 29)
501
502 ;;; %RANDOM-INTEGER -- Internal
503 ;;;
504 ;;; Make a random integer. We extract the most significant bits from
505 ;;; the generator and concatenate enough of them together to make up
506 ;;; the desired integer.
507 (defun %random-integer (arg state)
508 (declare (type (integer 1) arg)
509 (type random-state state))
510 (flet ((random-chunk (state)
511 (values (truncate (* (float (expt 2 random-chunk-size) 1d0)
512 (new-random-float state))))))
513 (let ((shift random-chunk-size))
514 (do ((bits (random-chunk state)
515 (logxor (ash bits shift) (random-chunk state)))
516 (count (+ (integer-length arg)
517 (- 0 shift))
518 (- count shift)))
519 ((minusp count)
520 (rem bits arg))
521 (declare (fixnum count))))))
522
523 (defun random (arg &optional (state *random-state*))
524 "Generate a uniformly distributed pseudo-random number between zero
525 and Arg. State, if supplied, is the random state to use."
526 (declare (inline %random-single-float %random-double-float))
527 (cond
528 ((<= arg 0)
529 (error 'simple-type-error
530 :expected-type '(or (integer 1) (float (0)))
531 :datum arg
532 :format-control "Argument is not a positive integer or a positive float: ~S"
533 :format-arguments (list arg)))
534 ((typep arg 'fixnum)
535 (values (truncate (%random-double-float (coerce arg 'double-float)
536 state))))
537 ((typep arg 'single-float)
538 (%random-single-float arg state))
539 ((typep arg 'double-float)
540 (%random-double-float arg state))
541 ((integerp arg)
542 (%random-integer arg state))
543 (t
544 (error 'simple-type-error
545 :expected-type '(or (integer 1) (float (0)))
546 :datum arg
547 :format-control "Argument is not a positive integer or a positive float: ~S"
548 :format-arguments (list arg)))))
549
550 ) ; end progn

  ViewVC Help
Powered by ViewVC 1.1.5