-;;; 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-WORDS words sampled from the Mersenne twister
-(declaim (inline %random-words))
-(defun %random-words (n-words 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-words
- for result = 0 then (logior (ash result sb!vm:n-word-bits)
- (%random-word 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-words n-extra-bits)
- (floor n-bits sb!vm:n-word-bits)
- (let ((full-chunks (%random-words n-full-words state)))
- (if (zerop n-extra-bits)
- full-chunks
- (logior full-chunks
- (ash (logand (%random-word state)
- (1- (ash 1 n-extra-bits)))
- (* n-full-words
- sb!vm:n-word-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))
- (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-word state) mask)
- inclusive-limit)))
-\f
-;;;; outer, dynamically-typed interface