X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Ftarget-random.lisp;h=f0b740fd3afb4e501ec1058ce011d286d11ff042;hb=95591ed483dbb8c0846c129953acac1554f28809;hp=3dc3325bc33522565e071e6e6f97fd92c1c0960a;hpb=0a7604d54581d2c846838c26ce6a7993629586fa;p=sbcl.git diff --git a/src/code/target-random.lisp b/src/code/target-random.lisp index 3dc3325..f0b740f 100644 --- a/src/code/target-random.lisp +++ b/src/code/target-random.lisp @@ -69,7 +69,7 @@ (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." @@ -164,6 +164,13 @@ (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)))) ;;; 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 @@ -179,7 +186,7 @@ (* arg (- (make-single-float (dpb (ash (random-chunk state) - (- sb!vm:single-float-digits n-random-chunk-bits)) + (- sb!vm:single-float-digits random-chunk-length)) sb!vm:single-float-significand-byte (single-float-bits 1.0))) 1.0))) @@ -202,7 +209,7 @@ (* arg (- (sb!impl::make-double-float (dpb (ash (random-chunk state) - (- sb!vm:double-float-digits n-random-chunk-bits 32)) + (- sb!vm:double-float-digits random-chunk-length 32)) sb!vm:double-float-significand-byte (sb!impl::double-float-high-bits 1d0)) (random-chunk state)) @@ -217,7 +224,7 @@ (* arg (- (sb!impl::make-double-float (dpb (ash (sb!vm::random-mt19937 state-vector) - (- sb!vm:double-float-digits n-random-chunk-bits + (- sb!vm:double-float-digits random-chunk-length sb!vm:n-word-bits)) sb!vm:double-float-significand-byte (sb!impl::double-float-high-bits 1d0)) @@ -227,130 +234,24 @@ ;;;; random integers -;;; 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 -;;; : -;;; 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 - ;; : - ;; 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))) - -;;;; outer, dynamically-typed interface +(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))))) (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) (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 (fixnump arg) (<= arg random-fixnum-max) (> arg 0)) + (rem (random-chunk state) arg)) ((and (typep arg 'single-float) (> arg 0.0f0)) (%random-single-float arg state)) ((and (typep arg 'double-float) (> arg 0.0d0)) @@ -358,10 +259,8 @@ #!+long-float ((and (typep arg 'long-float) (> arg 0.0l0)) (%random-long-float 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))) + ((and (integerp arg) (> arg 0)) + (%random-integer arg state)) (t (error 'simple-type-error :expected-type '(or (integer 1) (float (0))) :datum arg