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

Contents of /src/code/rand.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.6 - (hide annotations)
Thu Jun 24 12:57:50 1993 UTC (20 years, 9 months ago) by ram
Branch: MAIN
Changes since 1.5: +2 -2 lines
format-string -> format-control
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     ;;; If you want to use this code or any part of CMU Common Lisp, please contact
7     ;;; Scott Fahlman or slisp-group@cs.cmu.edu.
8     ;;;
9     (ext:file-comment
10 ram 1.6 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/rand.lisp,v 1.6 1993/06/24 12:57:50 ram Exp $")
11 ram 1.3 ;;;
12 ram 1.1 ;;; **********************************************************************
13     ;;;
14 ram 1.5 ;;; Functions to random number functions for Spice Lisp
15 ram 1.1 ;;;
16 ram 1.5 ;;; Originally written by David Adam. Python tuning, better large integer
17     ;;; randomness and efficient IEEE float support by Rob MacLachlan.
18 ram 1.1 ;;;
19 ram 1.5 (in-package "LISP")
20 ram 1.1 (export '(random-state random-state-p random *random-state*
21     make-random-state))
22    
23 ram 1.5 (in-package "KERNEL")
24     (export '(%random-single-float %random-double-float random-chunk
25     random-fixnum-max))
26    
27    
28     ;;;; Random state hackery:
29    
30 ram 1.1 (defconstant random-const-a 8373)
31     (defconstant random-const-c 101010101)
32     (defconstant random-max 54)
33 ram 1.5
34     ;;; Inclusive upper bound on the size of fixnum kept in the state (and returned
35     ;;; by random-chunk.) Must be even.
36     ;;;
37     (defconstant random-upper-bound (- most-positive-fixnum 3))
38     (defconstant random-chunk-length (integer-length random-upper-bound))
39     (deftype random-chunk () `(integer 0 ,random-upper-bound))
40    
41 ram 1.1 (defvar rand-seed 0)
42 wlott 1.4 (defstruct (random-state
43     (:constructor make-random-object)
44     (:make-load-form-fun :just-dump-it-normally))
45 ram 1.5 (j 24 :type index)
46     (k 0 :type index)
47 ram 1.1 (seed (make-array (1+ random-max) :initial-contents
48     (do ((list-rands () (cons (rand1) list-rands))
49     (i 0 (1+ i)))
50 ram 1.5 ((> i random-max) list-rands)
51     (declare (fixnum i))))
52 ram 1.1 :type simple-vector))
53    
54    
55     ;;; Generates a random number from rand-seed.
56     (defun rand1 ()
57 ram 1.5 (declare (optimize (inhibit-warnings 3)))
58     (setq rand-seed
59     (mod (+ (* rand-seed random-const-a) random-const-c)
60     (1+ random-upper-bound))))
61 wlott 1.2
62    
63     (defvar *random-state* (make-random-object))
64    
65 ram 1.1
66     (defun copy-state (cur-state)
67     (let ((state (make-random-object
68     :seed (make-array 55)
69     :j (random-state-j cur-state)
70     :k (random-state-k cur-state))))
71     (do ((i 0 (1+ i)))
72     ((= i 55) state)
73     (declare (fixnum i))
74     (setf (aref (random-state-seed state) i)
75     (aref (random-state-seed cur-state) i)))))
76    
77     (defun make-random-state (&optional state)
78     "Make a random state object. If State is not supplied, return a copy
79     of the default random state. If State is a random state, then return a
80     copy of it. If state is T then return a random state generated from
81     the universal time."
82     (cond ((not state) (copy-state *random-state*))
83     ((random-state-p state) (copy-state state))
84     ((eq state t) (setq rand-seed (get-universal-time))
85     (make-random-object))
86 ram 1.5 (t (error "Argument is not a RANDOM-STATE, T or NIL: ~S" state))))
87 ram 1.1
88 ram 1.5
89     ;;;; Random entries:
90    
91     (declaim (start-block random %random-single-float %random-double-float
92     random-chunk))
93    
94     ;;; random-chunk -- Internal
95     ;;;
96     ;;; This function generates fixnums between 0 and random-upper-bound,
97     ;;; inclusive. For the algorithm to work random-upper-bound must be an
98     ;;; even positive fixnum. State is the random state to use.
99     ;;;
100     (declaim (ftype (function (random-state) random-chunk) random-chunk))
101     (defun random-chunk (state)
102     (let* ((seed (random-state-seed state))
103     (j (random-state-j state))
104     (k (random-state-k state))
105     (a (- (- random-upper-bound
106     (the random-chunk
107     (svref seed
108     (setf (random-state-j state)
109     (if (= j 0) random-max (1- j))))))
110     (the random-chunk
111     (svref seed
112     (setf (random-state-k state)
113     (if (= k 0) random-max (1- k))))))))
114     (declare (fixnum a))
115     (setf (svref seed k)
116     (the random-chunk (if (minusp a) (- a) (- random-upper-bound a))))))
117    
118    
119     ;;; %RANDOM-SINGLE-FLOAT, %RANDOM-DOUBLE-FLOAT -- Interface
120     ;;;
121     ;;; Handle the single or double float case of RANDOM. We generate a float
122     ;;; between 0.0 and 1.0 by clobbering the significand of 1.0 with random bits,
123     ;;; then subtracting 1.0. This hides the fact that we have a hidden bit.
124     ;;;
125     (declaim (inline %random-single-float %random-double-float))
126     (defun %random-single-float (arg state)
127     (declare (type (single-float (0f0)) arg)
128     (type (or random-state null) state))
129     (* arg
130     (- (make-single-float
131     (dpb (ash (random-chunk (or state *random-state*))
132     (- vm:single-float-digits random-chunk-length))
133     vm:single-float-significand-byte
134     (single-float-bits 1.0)))
135     1.0)))
136     ;;;
137     (defun %random-double-float (arg state)
138     (declare (type (double-float (0d0)) arg)
139     (type (or random-state null) state))
140     (let ((state (or state *random-state*)))
141     (* arg
142     (- (make-double-float
143     (dpb (ash (random-chunk state)
144     (- vm:double-float-digits random-chunk-length
145     vm:word-bits))
146     vm:double-float-significand-byte
147     (double-float-high-bits 1d0))
148     (logxor (ash (random-chunk state)
149     (- vm:word-bits random-chunk-length))
150     (ash (random-chunk state)
151     (- random-chunk-length vm:word-bits))))
152     1d0))))
153    
154    
155     ;;;; Random integers:
156    
157     ;;; Amount we overlap chunks by when building a large integer to make up for
158     ;;; the loss of randomness in the low bits.
159     ;;;
160     (defconstant random-integer-overlap 3)
161    
162     ;;; Extra bits of randomness that we generate before taking the value MOD the
163     ;;; limit, to avoid loss of randomness near the limit.
164     ;;;
165     (defconstant random-integer-extra-bits 10)
166    
167     ;;; Largest fixnum we can compute from one chunk of bits.
168     ;;;
169     (defconstant random-fixnum-max
170     (1- (ash 1 (- random-chunk-length random-integer-extra-bits))))
171    
172    
173     ;;; %RANDOM-INTEGER -- Internal
174     ;;;
175     (defun %random-integer (arg state)
176     (declare (type (integer 1) arg) (type random-state state))
177     (let ((shift (- random-chunk-length random-integer-overlap)))
178     (do ((bits (random-chunk state)
179     (logxor (ash bits shift) (random-chunk state)))
180     (count (+ (integer-length arg)
181     (- random-integer-extra-bits shift))
182     (- count shift)))
183     ((minusp count)
184     (rem bits arg))
185     (declare (fixnum count)))))
186    
187 ram 1.1 (defun random (arg &optional (state *random-state*))
188     "Generate a uniformly distributed pseudo-random number between zero
189     and Arg. State, if supplied, is the random state to use."
190 ram 1.5 (declare (inline %random-single-float %random-double-float))
191     (cond
192     ((and (fixnump arg) (<= arg random-fixnum-max))
193     (rem (random-chunk state) arg))
194     ((typep arg 'single-float)
195     (%random-single-float arg state))
196     ((typep arg 'double-float)
197     (%random-double-float arg state))
198     ((integerp arg)
199     (%random-integer arg state))
200     (t
201     (error 'simple-type-error :expected-type '(real (0)) :datum arg
202 ram 1.6 :format-control "Argument is not a positive real number: ~S"
203 ram 1.5 :format-arguments (list arg)))))

  ViewVC Help
Powered by ViewVC 1.1.5