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

Contents of /src/code/rand.lisp

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.5