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

Contents of /src/code/rand.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.10 - (hide annotations)
Wed Jun 11 18:32:14 1997 UTC (16 years, 10 months ago) by dtc
Branch: MAIN
Changes since 1.9: +12 -9 lines
New-random tuning; with latest propagate-float-type code can inline to
32 bits on x86 and 31 on other ports.
1 ram 1.5 ;;; -*- Mode: Lisp; Package: Kernel -*-
2 ram 1.1 ;;;
3     ;;; **********************************************************************
4 ram 1.3 ;;; 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 dtc 1.10 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/rand.lisp,v 1.10 1997/06/11 18:32:14 dtc Exp $")
9 ram 1.3 ;;;
10 ram 1.1 ;;; **********************************************************************
11     ;;;
12 ram 1.5 ;;; Functions to random number functions for Spice Lisp
13 ram 1.1 ;;;
14 ram 1.5 ;;; Originally written by David Adam. Python tuning, better large integer
15     ;;; randomness and efficient IEEE float support by Rob MacLachlan.
16 ram 1.1 ;;;
17 dtc 1.8 #-new-random
18 ram 1.5 (in-package "LISP")
19 dtc 1.8 #-new-random
20 ram 1.1 (export '(random-state random-state-p random *random-state*
21     make-random-state))
22    
23 dtc 1.8 #-new-random
24 ram 1.5 (in-package "KERNEL")
25 dtc 1.8 #-new-random
26 ram 1.5 (export '(%random-single-float %random-double-float random-chunk
27     random-fixnum-max))
28    
29    
30     ;;;; Random state hackery:
31    
32 dtc 1.8 #-new-random
33     (progn
34 ram 1.1 (defconstant random-const-a 8373)
35     (defconstant random-const-c 101010101)
36     (defconstant random-max 54)
37 ram 1.5
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 ram 1.1 (defvar rand-seed 0)
46 wlott 1.4 (defstruct (random-state
47     (:constructor make-random-object)
48     (:make-load-form-fun :just-dump-it-normally))
49 ram 1.5 (j 24 :type index)
50     (k 0 :type index)
51 ram 1.1 (seed (make-array (1+ random-max) :initial-contents
52     (do ((list-rands () (cons (rand1) list-rands))
53     (i 0 (1+ i)))
54 ram 1.5 ((> i random-max) list-rands)
55     (declare (fixnum i))))
56 ram 1.1 :type simple-vector))
57    
58    
59     ;;; Generates a random number from rand-seed.
60     (defun rand1 ()
61 ram 1.5 (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 wlott 1.2
66    
67     (defvar *random-state* (make-random-object))
68    
69 ram 1.1
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 ram 1.5 (t (error "Argument is not a RANDOM-STATE, T or NIL: ~S" state))))
91 ram 1.1
92 ram 1.5
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 ram 1.1 (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 ram 1.5 (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 ram 1.6 :format-control "Argument is not a positive real number: ~S"
207 ram 1.5 :format-arguments (list arg)))))
208 dtc 1.8
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 dtc 1.10 (export '(%random-single-float %random-double-float))
243 dtc 1.8
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 dtc 1.10 ;;; Random integers are created by concatenating a bunch of chunks
483     ;;; together. Chunks are generated without consing by truncating a
484     ;;; random float to an integer. With the latest propagate-float-type
485     ;;; code the compiler can inline truncate (signed-byte 32) allowing 31
486     ;;; bits, and (unsigned-byte 32) 32 bits on the x86. When not using the
487     ;;; propagate-float-type feature the best size that can be inlined is
488     ;;; 29 bits. Use of sizes greater than 29 bits causes consing in
489     ;;; argument passing to generic functions, so 29 bits is a good
490     ;;; default.
491 dtc 1.8 ;;;
492     (defconstant random-chunk-size 29)
493    
494     ;;; %RANDOM-INTEGER -- Internal
495     ;;;
496     ;;; Make a random integer. We extract the most significant bits from
497     ;;; the generator and concatenate enough of them together to make up
498     ;;; the desired integer.
499     (defun %random-integer (arg state)
500     (declare (type (integer 1) arg)
501     (type random-state state))
502     (flet ((random-chunk (state)
503 dtc 1.10 (values (truncate (* (float (expt 2 random-chunk-size) 1d0)
504 dtc 1.9 (new-random-float state))))))
505 dtc 1.8 (let ((shift random-chunk-size))
506     (do ((bits (random-chunk state)
507     (logxor (ash bits shift) (random-chunk state)))
508     (count (+ (integer-length arg)
509     (- 0 shift))
510     (- count shift)))
511     ((minusp count)
512     (rem bits arg))
513     (declare (fixnum count))))))
514    
515     (defun random (arg &optional (state *random-state*))
516     "Generate a uniformly distributed pseudo-random number between zero
517     and Arg. State, if supplied, is the random state to use."
518     (declare (inline %random-single-float %random-double-float))
519     (cond
520     ((<= arg 0)
521     (error 'simple-type-error
522     :expected-type '(or (integer 1) (float (0)))
523     :datum arg
524     :format-control "Argument is not a positive integer or a positive float: ~S"
525     :format-arguments (list arg)))
526     ((typep arg 'fixnum)
527     (values (truncate (%random-double-float (coerce arg 'double-float)
528     state))))
529     ((typep arg 'single-float)
530     (%random-single-float arg state))
531     ((typep arg 'double-float)
532     (%random-double-float arg state))
533     ((integerp arg)
534     (%random-integer arg state))
535     (t
536     (error 'simple-type-error
537     :expected-type '(or (integer 1) (float (0)))
538     :datum arg
539     :format-control "Argument is not a positive integer or a positive float: ~S"
540     :format-arguments (list arg)))))
541    
542     ) ; end progn

  ViewVC Help
Powered by ViewVC 1.1.5