(defun make-random-state (&optional state)
#!+sb-doc
- "Make a random state object. If STATE is not supplied, return a copy
+ "Make a RANDOM-STATE object. If STATE is not supplied, return a copy
of the default random state. If STATE is a random state, then return a
copy of it. If STATE is T then return a random state generated from
the universal time."
(defun random-chunk (state)
(declare (type random-state state))
(sb!vm::random-mt19937 (random-state-state state)))
-
-#!-sb-fluid (declaim (inline big-random-chunk))
-(defun big-random-chunk (state)
- (declare (type random-state state))
- (logand (1- (expt 2 64))
- (logior (ash (random-chunk state) 32)
- (random-chunk state))))
\f
;;; Handle the single or double float case of RANDOM. We generate a
;;; float between 0.0 and 1.0 by clobbering the significand of 1.0
(* arg
(- (make-single-float
(dpb (ash (random-chunk state)
- (- sb!vm:single-float-digits random-chunk-length))
+ (- sb!vm:single-float-digits n-random-chunk-bits))
sb!vm:single-float-significand-byte
(single-float-bits 1.0)))
1.0)))
(* arg
(- (sb!impl::make-double-float
(dpb (ash (random-chunk state)
- (- sb!vm:double-float-digits random-chunk-length 32))
+ (- sb!vm:double-float-digits n-random-chunk-bits 32))
sb!vm:double-float-significand-byte
(sb!impl::double-float-high-bits 1d0))
(random-chunk state))
(* arg
(- (sb!impl::make-double-float
(dpb (ash (sb!vm::random-mt19937 state-vector)
- (- sb!vm:double-float-digits random-chunk-length
+ (- sb!vm:double-float-digits n-random-chunk-bits
sb!vm:n-word-bits))
sb!vm:double-float-significand-byte
(sb!impl::double-float-high-bits 1d0))
\f
;;;; random integers
-(defun %random-integer (arg state)
- (declare (type (integer 1) arg) (type random-state state))
- (let ((shift (- random-chunk-length random-integer-overlap)))
- (do ((bits (random-chunk state)
- (logxor (ash bits shift) (random-chunk state)))
- (count (+ (integer-length arg)
- (- random-integer-extra-bits shift))
- (- count shift)))
- ((minusp count)
- (rem bits arg))
- (declare (fixnum count)))))
+;;; a bitmask M wide enough that (= (LOGAND INCLUSIVE-LIMIT M) INCLUSIVE-LIMIT)
+(declaim (inline %inclusive-random-integer-mask))
+(defun %inclusive-random-integer-mask (inclusive-limit)
+ (1- (ash 1 (integer-length inclusive-limit))))
+
+;;; Transform a uniform sample from an at-least-as-large range into a
+;;; random sample in [0,INCLUSIVE-LIMIT] range by throwing away
+;;; too-big samples.
+;;;
+;;; Up through sbcl-1.0.4, the (RANDOM INTEGER) worked by taking (MOD
+;;; RAW-MERSENNE-OUTPUT INTEGER). That introduces enough bias that it
+;;; is unsuitable for some calculations (like the Metropolis Monte
+;;; Carlo simulations that I, WHN, worked on in grad school): in the
+;;; sbcl-1.0.4, the expectation value of a sample was (in the worst
+;;; part of the range) more than 1.0001 times the ideal expectation
+;;; value. Perhaps that was even ANSI-conformant, because ANSI says only
+;;; An approximately uniform choice distribution is used. If
+;;; LIMIT is an integer, each of the possible results occurs
+;;; with (approximate) probability 1/LIMIT.
+;;; But I (WHN again) said "yuck," so these days we try to get a
+;;; truly uniform distribution by translating RAW-MERSENNE-OUTPUT to
+;;; our output using the second method recommended by
+;;; <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/efaq.html>:
+;;; In the rare case that the rounding-off error becomes a problem,
+;;; obtain minimum integer n such that N<=2^n, generate integer
+;;; random numbers, take the most significant n bits, and discard
+;;; those more than or equal to N.
+;;; (The first method recommended there differs from the sbcl-1.0.4
+;;; method: the recommended method gives slightly different weights
+;;; to different integers, but distributes the overweighted integers
+;;; evenly across the range of outputs so that any skew would be
+;;; of order (/ MOST-POSITIVE-MACHINE-WORD), rather than the (/ 10000)
+;;; or so achieved by sbcl-1.0.4. That skew would probably be
+;;; acceptable --- it seems to be exactly the kind of deviation that
+;;; might have been anticipated in the ANSI CL standard. However, that
+;;; recommended calculation involves multiplication and division of
+;;; two-machine-word bignums, which is hard for us to do efficiently
+;;; here. Without special low-level hacking to support such short
+;;; bignums without consing, the accept-reject solution is not only
+;;; more exact, but also likely more efficient.)
+(defmacro %inclusive-random-integer-accept-reject (raw-mersenne-output-expr
+ inclusive-limit)
+ (with-unique-names (raw-mersenne-output inclusive-limit-once)
+ `(loop
+ with ,inclusive-limit-once = ,inclusive-limit
+ for ,raw-mersenne-output = ,raw-mersenne-output-expr
+ until (<= ,raw-mersenne-output ,inclusive-limit-once)
+ finally (return ,raw-mersenne-output))))
+
+;;; an UNSIGNED-BYTE of N-CHUNKS chunks sampled from the Mersenne twister
+(declaim (inline %random-chunks))
+(defun %random-chunks (n-chunks state)
+ ;; KLUDGE: This algorithm will cons O(N^2) words when constructing
+ ;; an N-word result. To do better should be fundamentally easy, we
+ ;; just need to do some low-level hack preallocating the bignum and
+ ;; writing its words one by one.
+ ;;
+ ;; Note: The old sbcl-1.0.4 code did its analogous operation using a
+ ;; mysterious RANDOM-INTEGER-OVERLAP parameter, "the amount that we
+ ;; overlap chunks by when building a large integer to make up for
+ ;; the loss of randomness in the low bits." No such thing is called
+ ;; for on
+ ;; <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/efaq.html>:
+ ;; they say there that it's OK just to concatenate words of twister
+ ;; output with no overlap. Thus, crossing our fingers and hoping that
+ ;; the previous RNG author didn't have some subtle reason to need
+ ;; RANDOM-INTEGER-OVERLAP that we know not of, we just concatenate
+ ;; chunks.
+ (loop repeat n-chunks
+ for result = 0 then (logior (ash result n-random-chunk-bits)
+ (random-chunk state))
+ finally (return result)))
+
+;;; an UNSIGNED-BYTE of N-BITS bits sampled from the Mersenne twister
+(declaim (inline %random-bits))
+(defun %random-bits (n-bits state)
+ (multiple-value-bind (n-full-chunks n-extra-bits)
+ (floor n-bits n-random-chunk-bits)
+ (let ((full-chunks (%random-chunks n-full-chunks state)))
+ (if (zerop n-extra-bits)
+ full-chunks
+ (logior full-chunks
+ (ash (logand (random-chunk state)
+ (1- (ash 1 n-extra-bits)))
+ (* n-full-chunks n-random-chunk-bits)))))))
+
+;;; the guts of (RANDOM (1+ INCLUSIVE-LIMIT))
+(defun %inclusive-random-integer (inclusive-limit state)
+ (declare (optimize speed (space 1))) ; to ensure DEFTRANSFORM is enabled
+ (%inclusive-random-integer inclusive-limit state))
+
+;;; the guts of RANDOM for the known-to-return-FIXNUM case
+;;;
+;;; We use INCLUSIVE-LIMIT instead of the (exclusive) LIMIT of RANDOM,
+;;; because we want it to be a FIXNUM even for the possibly-common
+;;; case of (RANDOM (1+ MOST-POSITIVE-FIXNUM)). (That case is what
+;;; one might use for generating random hash values, e.g.)
+;;; It also turns out to be just what's needed for INTEGER-LENGTH.
+(declaim (maybe-inline %inclusive-random-fixnum))
+(defun %inclusive-random-fixnum (inclusive-limit state)
+ (declare (type (and fixnum unsigned-byte) inclusive-limit))
+ (aver (<= inclusive-limit most-positive-random-chunk))
+ (let (;; If this calculation needs to be optimized further, a good
+ ;; start might be a DEFTRANSFORM which picks off the case of
+ ;; constant LIMIT and precomputes the MASK at compile time.
+ (mask (%inclusive-random-integer-mask inclusive-limit)))
+ (%inclusive-random-integer-accept-reject (logand (random-chunk state) mask)
+ inclusive-limit)))
+\f
+;;;; outer, dynamically-typed interface
(defun random (arg &optional (state *random-state*))
- (declare (inline %random-single-float %random-double-float
+ (declare (inline %random-single-float
+ %random-double-float
#!+long-float %random-long-float))
(cond
- ((and (fixnump arg) (<= arg random-fixnum-max) (> arg 0))
- (rem (random-chunk state) arg))
+ ((and (fixnump arg) (plusp arg))
+ (locally
+ ;; The choice to inline this very common case of
+ ;; %INCLUSIVE-RANDOM-FIXNUM and not the less-common call
+ ;; below was just a guess (by WHN 2007-03-27), not based on
+ ;; benchmarking or anything.
+ (declare (inline %inclusive-random-fixnum))
+ (%inclusive-random-fixnum (1- arg) state)))
((and (typep arg 'single-float) (> arg 0.0f0))
(%random-single-float arg state))
((and (typep arg 'double-float) (> arg 0.0d0))
#!+long-float
((and (typep arg 'long-float) (> arg 0.0l0))
(%random-long-float arg state))
- ((and (integerp arg) (> arg 0))
- (%random-integer arg state))
+ ((and (integerp arg) (plusp arg))
+ (if (= arg (1+ most-positive-fixnum))
+ (%inclusive-random-fixnum most-positive-fixnum state)
+ (%inclusive-random-integer (1- arg) state)))
(t
(error 'simple-type-error
:expected-type '(or (integer 1) (float (0))) :datum arg
'(,fun num (or state *random-state*)))))
(frob %random-single-float single-float)
(frob %random-double-float double-float))
-
-;;; Mersenne Twister RNG
-;;;
-;;; FIXME: It's unpleasant to have RANDOM functionality scattered
-;;; through the code this way. It would be nice to move this into the
-;;; same file as the other RANDOM definitions.
-(deftransform random ((num &optional state)
- ((integer 1 #.(expt 2 sb!vm::n-word-bits)) &optional *))
- ;; FIXME: I almost conditionalized this as #!+sb-doc. Find some way
- ;; of automatically finding #!+sb-doc in proximity to DEFTRANSFORM
- ;; to let me scan for places that I made this mistake and didn't
- ;; catch myself.
- "use inline (UNSIGNED-BYTE 32) operations"
- (let ((type (lvar-type num))
- (limit (expt 2 sb!vm::n-word-bits))
- (random-chunk (ecase sb!vm::n-word-bits
- (32 'random-chunk)
- (64 'sb!kernel::big-random-chunk))))
- (if (numeric-type-p type)
- (let ((num-high (numeric-type-high (lvar-type num))))
- (aver num-high)
- (cond ((constant-lvar-p num)
- ;; Check the worst case sum absolute error for the
- ;; random number expectations.
- (let ((rem (rem limit num-high)))
- (unless (< (/ (* 2 rem (- num-high rem))
- num-high limit)
- (expt 2 (- sb!kernel::random-integer-extra-bits)))
- (give-up-ir1-transform
- "The random number expectations are inaccurate."))
- (if (= num-high limit)
- `(,random-chunk (or state *random-state*))
- #!-(or x86 x86-64)
- `(rem (,random-chunk (or state *random-state*)) num)
- #!+(or x86 x86-64)
- ;; Use multiplication, which is faster.
- `(values (sb!bignum::%multiply
- (,random-chunk (or state *random-state*))
- num)))))
- ((> num-high random-fixnum-max)
- (give-up-ir1-transform
- "The range is too large to ensure an accurate result."))
- #!+(or x86 x86-64)
- ((< num-high limit)
- `(values (sb!bignum::%multiply
- (,random-chunk (or state *random-state*))
- num)))
- (t
- `(rem (,random-chunk (or state *random-state*)) num))))
- ;; KLUDGE: a relatively conservative treatment, but better
- ;; than a bug (reported by PFD sbcl-devel towards the end of
- ;; 2004-11.
- '(rem (random-chunk (or state *random-state*)) num))))
\f
;;;; float accessors
projects / sbcl.git / commitdiff