/[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 - (hide 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 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 rtoy 1.11.54.1 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/rand.lisp,v 1.11.54.1 2010/02/08 17:15:48 rtoy 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 rtoy 1.11.54.1 (intl:textdomain "cmucl")
20    
21 dtc 1.8 #-new-random
22 ram 1.1 (export '(random-state random-state-p random *random-state*
23     make-random-state))
24    
25 dtc 1.8 #-new-random
26 ram 1.5 (in-package "KERNEL")
27 dtc 1.8 #-new-random
28 ram 1.5 (export '(%random-single-float %random-double-float random-chunk
29     random-fixnum-max))
30    
31    
32     ;;;; Random state hackery:
33    
34 dtc 1.8 #-new-random
35     (progn
36 ram 1.1 (defconstant random-const-a 8373)
37     (defconstant random-const-c 101010101)
38     (defconstant random-max 54)
39 ram 1.5
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 ram 1.1 (defvar rand-seed 0)
48 wlott 1.4 (defstruct (random-state
49     (:constructor make-random-object)
50     (:make-load-form-fun :just-dump-it-normally))
51 ram 1.5 (j 24 :type index)
52     (k 0 :type index)
53 ram 1.1 (seed (make-array (1+ random-max) :initial-contents
54     (do ((list-rands () (cons (rand1) list-rands))
55     (i 0 (1+ i)))
56 ram 1.5 ((> i random-max) list-rands)
57     (declare (fixnum i))))
58 ram 1.1 :type simple-vector))
59    
60    
61     ;;; Generates a random number from rand-seed.
62     (defun rand1 ()
63 ram 1.5 (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 wlott 1.2
68    
69     (defvar *random-state* (make-random-object))
70    
71 ram 1.1
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 toy 1.11 (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 ram 1.1
98 ram 1.5
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 ram 1.1 (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 ram 1.5 (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 ram 1.6 :format-control "Argument is not a positive real number: ~S"
213 ram 1.5 :format-arguments (list arg)))))
214 dtc 1.8
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 dtc 1.10 (export '(%random-single-float %random-double-float))
249 dtc 1.8
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 toy 1.11 (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 dtc 1.8
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 dtc 1.10 ;;; 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 dtc 1.8 ;;;
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 dtc 1.10 (values (truncate (* (float (expt 2 random-chunk-size) 1d0)
512 dtc 1.9 (new-random-float state))))))
513 dtc 1.8 (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