From: William Harold Newman <wnewman@us...>  20070327 16:21:05

Update of /cvsroot/sbcl/sbcl/src/code In directory sc8prcvs8.sourceforge.net:/tmp/cvsserv11449/src/code Modified Files: random.lisp targetrandom.lisp Log Message: 1.0.4.6: rewrote RANDOM for INTEGER case to eliminate (small) bias Alas, it doesn't seem practical to test this properly. (Difficulty of building a nice test suite seems to be a pervasive problem for randomized code.) I could write a test case narrowly targeted to detect the old bug, but I don't think I could make even that single test simultaneously fast enough and reliable enough that it would be a good candidate for runtests.sh. And aiming at a more comprehensive set of tests just makes the makingitfastenough problem that much harder. Index: random.lisp =================================================================== RCS file: /cvsroot/sbcl/sbcl/src/code/random.lisp,v retrieving revision 1.5 retrieving revision 1.6 diff u d r1.5 r1.6  random.lisp 14 Jul 2005 16:30:38 0000 1.5 +++ random.lisp 27 Mar 2007 16:20:43 0000 1.6 @@ 9,21 +9,15 @@ (inpackage "SB!KERNEL") ;;; the size of the chunks returned by RANDOMCHUNK (def!constant randomchunklength 32)  ;;; the amount that we overlap chunks by when building a large integer ;;; to make up for the loss of randomness in the low bits (def!constant randomintegeroverlap 3)  ;;; extra bits of randomness that we generate before taking the value MOD the ;;; limit, to avoid loss of randomness near the limit (def!constant randomintegerextrabits 10)  ;;; the largest fixnum we can compute from one chunk of bits (def!constant randomfixnummax  (1 (ash 1 ( randomchunklength randomintegerextrabits)))) +;;; the size of the chunks returned from the fundamental random number +;;; generator +(def!constant nrandomchunkbits 32) +(def!constant mostpositiverandomchunk + (1 (ash 1 nrandomchunkbits))) +;;; our implementation of the RANDOMSTATE type specified by ANSI CL (sb!xc:defstruct (randomstate (:constructor %makerandomstate)  (:copier nil)) ; since shallow copy is wrong + ;; Shallow copy would be wrong: we must + ;; effectively COPYSEQ the STATE slot. + (:copier nil)) (state (initrandomstate) :type (simplearray (unsignedbyte 32) (627)))) Index: targetrandom.lisp =================================================================== RCS file: /cvsroot/sbcl/sbcl/src/code/targetrandom.lisp,v retrieving revision 1.14 retrieving revision 1.15 diff u d r1.14 r1.15  targetrandom.lisp 14 Jul 2005 16:30:40 0000 1.14 +++ targetrandom.lisp 27 Mar 2007 16:20:46 0000 1.15 @@ 69,7 +69,7 @@ (defun makerandomstate (&optional state) #!+sbdoc  "Make a random state object. If STATE is not supplied, return a copy + "Make a RANDOMSTATE 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,13 +164,6 @@ (defun randomchunk (state) (declare (type randomstate state)) (sb!vm::randommt19937 (randomstatestate state)))  #!sbfluid (declaim (inline bigrandomchunk)) (defun bigrandomchunk (state)  (declare (type randomstate state))  (logand (1 (expt 2 64))  (logior (ash (randomchunk state) 32)  (randomchunk 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 @@ 186,7 +179,7 @@ (* arg ( (makesinglefloat (dpb (ash (randomchunk state)  ( sb!vm:singlefloatdigits randomchunklength)) + ( sb!vm:singlefloatdigits nrandomchunkbits)) sb!vm:singlefloatsignificandbyte (singlefloatbits 1.0))) 1.0))) @@ 209,7 +202,7 @@ (* arg ( (sb!impl::makedoublefloat (dpb (ash (randomchunk state)  ( sb!vm:doublefloatdigits randomchunklength 32)) + ( sb!vm:doublefloatdigits nrandomchunkbits 32)) sb!vm:doublefloatsignificandbyte (sb!impl::doublefloathighbits 1d0)) (randomchunk state)) @@ 224,7 +217,7 @@ (* arg ( (sb!impl::makedoublefloat (dpb (ash (sb!vm::randommt19937 statevector)  ( sb!vm:doublefloatdigits randomchunklength + ( sb!vm:doublefloatdigits nrandomchunkbits sb!vm:nwordbits)) sb!vm:doublefloatsignificandbyte (sb!impl::doublefloathighbits 1d0)) @@ 234,24 +227,130 @@ ;;;; random integers (defun %randominteger (arg state)  (declare (type (integer 1) arg) (type randomstate state))  (let ((shift ( randomchunklength randomintegeroverlap)))  (do ((bits (randomchunk state)  (logxor (ash bits shift) (randomchunk state)))  (count (+ (integerlength arg)  ( randomintegerextrabits shift))  ( count shift)))  ((minusp count)  (rem bits arg))  (declare (fixnum count))))) +;;; a bitmask M wide enough that (= (LOGAND INCLUSIVELIMIT M) INCLUSIVELIMIT) +(declaim (inline %inclusiverandomintegermask)) +(defun %inclusiverandomintegermask (inclusivelimit) + (1 (ash 1 (integerlength inclusivelimit)))) + +;;; Transform a uniform sample from an atleastaslarge range into a +;;; random sample in [0,INCLUSIVELIMIT] range by throwing away +;;; toobig samples. +;;; +;;; Up through sbcl1.0.4, the (RANDOM INTEGER) worked by taking (MOD +;;; RAWMERSENNEOUTPUT 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 +;;; sbcl1.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 ANSIconformant, 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 RAWMERSENNEOUTPUT to +;;; our output using the second method recommended by +;;; <http://www.math.sci.hiroshimau.ac.jp/~mmat/MT/efaq.html>;: +;;; In the rare case that the roundingoff 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 sbcl1.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 (/ MOSTPOSITIVEMACHINEWORD), rather than the (/ 10000) +;;; or so achieved by sbcl1.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 +;;; twomachineword bignums, which is hard for us to do efficiently +;;; here. Without special lowlevel hacking to support such short +;;; bignums without consing, the acceptreject solution is not only +;;; more exact, but also likely more efficient.) +(defmacro %inclusiverandomintegeracceptreject (rawmersenneoutputexpr + inclusivelimit) + (withuniquenames (rawmersenneoutput inclusivelimitonce) + `(loop + with ,inclusivelimitonce = ,inclusivelimit + for ,rawmersenneoutput = ,rawmersenneoutputexpr + until (<= ,rawmersenneoutput ,inclusivelimitonce) + finally (return ,rawmersenneoutput)))) + +;;; an UNSIGNEDBYTE of NCHUNKS chunks sampled from the Mersenne twister +(declaim (inline %randomchunks)) +(defun %randomchunks (nchunks state) + ;; KLUDGE: This algorithm will cons O(N^2) words when constructing + ;; an Nword result. To do better should be fundamentally easy, we + ;; just need to do some lowlevel hack preallocating the bignum and + ;; writing its words one by one. + ;; + ;; Note: The old sbcl1.0.4 code did its analogous operation using a + ;; mysterious RANDOMINTEGEROVERLAP 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.hiroshimau.ac.jp/~mmat/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 + ;; RANDOMINTEGEROVERLAP that we know not of, we just concatenate + ;; chunks. + (loop repeat nchunks + for result = 0 then (logior (ash result nrandomchunkbits) + (randomchunk state)) + finally (return result))) + +;;; an UNSIGNEDBYTE of NBITS bits sampled from the Mersenne twister +(declaim (inline %randombits)) +(defun %randombits (nbits state) + (multiplevaluebind (nfullchunks nextrabits) + (floor nbits nrandomchunkbits) + (let ((fullchunks (%randomchunks nfullchunks state))) + (if (zerop nextrabits) + fullchunks + (logior fullchunks + (ash (logand (randomchunk state) + (1 (ash 1 nextrabits))) + (* nfullchunks nrandomchunkbits))))))) + +;;; the guts of (RANDOM (1+ INCLUSIVELIMIT)) +(defun %inclusiverandominteger (inclusivelimit state) + (declare (optimize speed (space 1))) ; to ensure DEFTRANSFORM is enabled + (%inclusiverandominteger inclusivelimit state)) + +;;; the guts of RANDOM for the knowntoreturnFIXNUM case +;;; +;;; We use INCLUSIVELIMIT instead of the (exclusive) LIMIT of RANDOM, +;;; because we want it to be a FIXNUM even for the possiblycommon +;;; case of (RANDOM (1+ MOSTPOSITIVEFIXNUM)). (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 INTEGERLENGTH. +(declaim (maybeinline %inclusiverandomfixnum)) +(defun %inclusiverandomfixnum (inclusivelimit state) + (declare (type (and fixnum unsignedbyte) inclusivelimit)) + (aver (<= inclusivelimit mostpositiverandomchunk)) + (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 (%inclusiverandomintegermask inclusivelimit))) + (%inclusiverandomintegeracceptreject (logand (randomchunk state) mask) + inclusivelimit))) + +;;;; outer, dynamicallytyped interface (defun random (arg &optional (state *randomstate*))  (declare (inline %randomsinglefloat %randomdoublefloat + (declare (inline %randomsinglefloat + %randomdoublefloat #!+longfloat %randomlongfloat)) (cond  ((and (fixnump arg) (<= arg randomfixnummax) (> arg 0))  (rem (randomchunk state) arg)) + ((and (fixnump arg) (plusp arg)) + (locally + ;; The choice to inline this very common case of + ;; %INCLUSIVERANDOMFIXNUM and not the lesscommon call + ;; below was just a guess (by WHN 20070327), not based on + ;; benchmarking or anything. + (declare (inline %inclusiverandomfixnum)) + (%inclusiverandomfixnum (1 arg) state))) ((and (typep arg 'singlefloat) (> arg 0.0f0)) (%randomsinglefloat arg state)) ((and (typep arg 'doublefloat) (> arg 0.0d0)) @@ 259,8 +358,10 @@ #!+longfloat ((and (typep arg 'longfloat) (> arg 0.0l0)) (%randomlongfloat arg state))  ((and (integerp arg) (> arg 0))  (%randominteger arg state)) + ((and (integerp arg) (plusp arg)) + (if (= arg (1+ mostpositivefixnum)) + (%inclusiverandomfixnum mostpositivefixnum state) + (%inclusiverandominteger (1 arg) state))) (t (error 'simpletypeerror :expectedtype '(or (integer 1) (float (0))) :datum arg 