From: William Harold Newman <wnewman@us...>  20070328 17:59:31

Update of /cvsroot/sbcl/sbcl/src/code In directory sc8prcvs8.sourceforge.net:/tmp/cvsserv30996/src/code Modified Files: random.lisp targetrandom.lisp Log Message: 1.0.4.8: "What's old is new again" or "what was he thinking?" At least in the wisdom of hindsight I can see that my new RANDOM is not ready for release, so I have reverted it. This version should be essentially 1.0.4.5 with a new version number. (It also has an extra, orphaned file in src/compiler/integertran.lisp. When I revisit the problem later, that should be cleaned up one way or another, meanwhile I don't feel like deleting it and then likely recreating it later.) Sorry about the confusion.:( Index: random.lisp =================================================================== RCS file: /cvsroot/sbcl/sbcl/src/code/random.lisp,v retrieving revision 1.7 retrieving revision 1.8 diff u d r1.7 r1.8  random.lisp 28 Mar 2007 14:32:12 0000 1.7 +++ random.lisp 28 Mar 2007 17:59:18 0000 1.8 @@ 9,9 +9,21 @@ (inpackage "SB!KERNEL") ;;; our implementation of the RANDOMSTATE type specified by ANSI CL +;;; 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)))) + (sb!xc:defstruct (randomstate (:constructor %makerandomstate)  ;; Shallow copy would be wrong: we must  ;; effectively COPYSEQ the STATE slot.  (:copier nil)) + (:copier nil)) ; since shallow copy is wrong (state (initrandomstate) :type (simplearray (unsignedbyte 32) (627)))) Index: targetrandom.lisp =================================================================== RCS file: /cvsroot/sbcl/sbcl/src/code/targetrandom.lisp,v retrieving revision 1.16 retrieving revision 1.17 diff u d r1.16 r1.17  targetrandom.lisp 28 Mar 2007 14:32:12 0000 1.16 +++ targetrandom.lisp 28 Mar 2007 17:59:18 0000 1.17 @@ 69,7 +69,7 @@ (defun makerandomstate (&optional state) #!+sbdoc  "Make a RANDOMSTATE 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." @@ 102,8 +102,7 @@ ;;; This function generates a 32bit integer between 0 and #xffffffff ;;; inclusive. (defconstant nrandommt19937bits 32) (declaim (inline %random32)) +#!sbfluid (declaim (inline randomchunk)) ;;; portable implementation (defconstant mt19937n 624) (defconstant mt19937m 397) @@ 139,7 +138,7 @@ (ash y 1) (aref state (logand y 1))))) (values)) #!x86 (defun %random32 (state) +(defun randomchunk (state) (declare (type randomstate state)) (let* ((state (randomstatestate state)) (k (aref state 2))) @@ 162,29 +161,16 @@ ;;; My inclination is to get rid of the nonportable implementation ;;; unless the performance difference is just enormous. #!+x86 (defun %random32 (state) +(defun randomchunk (state) (declare (type randomstate state)) (sb!vm::randommt19937 (randomstatestate state))) (declaim (inline %randomword)) (defun %randomword (state)  ;; KLUDGE: This #.(ECASE ...) is not the most flexible and elegant  ;; construct one could imagine. It is intended as a quick fix to stand  ;; in for The Right Thing which can't be coded so quickly.  ;;  ;; The Right Thing: The Mersenne Twister has been generalized to 64  ;; bits, and seems likely to be generalized to any future common CPU  ;; word width as well. Thus, it should be straightforward to  ;; implement this as "Return one sample from the MT variant  ;; corresponding to SB!VM:NWORDBITS."  ;;  ;; Meanwhile: Mock that up by pasting together as many samples from  ;; the 32bit Mersenne Twister as necessary.  #.(ecase sb!vm:nwordbits  (32 '(%random32 state))  (64 '(logior  (%random32 state)  (ash (%random32 state) 32))))) +#!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 @@ 197,16 +183,10 @@ (defun %randomsinglefloat (arg state) (declare (type (singlefloat (0f0)) arg) (type randomstate state))  ;; KLUDGE: The hardwiredto32 hackery here could be replaced by a  ;; call to %RANDOMBITS if there were sufficiently smart  ;; DEFTRANSFORMs for it, but in sbcl1.0.4.epsilon it looks as  ;; though it would be a performance disaster.  (aver (<= sb!vm:singlefloatdigits 32)) (* arg ( (makesinglefloat  (dpb (ash  (%random32 state)  ( sb!vm:singlefloatdigits 32)) + (dpb (ash (randomchunk state) + ( sb!vm:singlefloatdigits randomchunklength)) sb!vm:singlefloatsignificandbyte (singlefloatbits 1.0))) 1.0))) @@ 214,21 +194,25 @@ (doublefloat 0d0)) %randomdoublefloat)) +;;; 32bit version +#!+nil +(defun %randomdoublefloat (arg state) + (declare (type (doublefloat (0d0)) arg) + (type randomstate state)) + (* (float (randomchunk state) 1d0) (/ 1d0 (expt 2 32)))) + ;;; 53bit version #!x86 (defun %randomdoublefloat (arg state) (declare (type (doublefloat (0d0)) arg) (type randomstate state))  ;; KLUDGE: As in %RANDOMSIMNGLEFLOAT, as of sbcl1.0.4.epsilon  ;; calling %RANDOMBITS doesn't look reasonable, so we bang bits.  (aver (<= sb!vm:singlefloatdigits 32)) (* arg ( (sb!impl::makedoublefloat  (dpb (ash (%random32 state)  ( sb!vm:doublefloatdigits 64)) + (dpb (ash (randomchunk state) + ( sb!vm:doublefloatdigits randomchunklength 32)) sb!vm:doublefloatsignificandbyte (sb!impl::doublefloathighbits 1d0))  (%random32 state)) + (randomchunk state)) 1d0))) ;;; using a faster inline VOP @@ 240,140 +224,34 @@ (* arg ( (sb!impl::makedoublefloat (dpb (ash (sb!vm::randommt19937 statevector)  ( sb!vm:doublefloatdigits  nrandommt19937bits + ( sb!vm:doublefloatdigits randomchunklength sb!vm:nwordbits)) sb!vm:doublefloatsignificandbyte (sb!impl::doublefloathighbits 1d0)) (sb!vm::randommt19937 statevector)) 1d0)))) + ;;;; random integers ;;; 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 NWORDS words sampled from the Mersenne twister (declaim (inline %randomwords)) (defun %randomwords (nwords 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 nwords  for result = 0 then (logior (ash result sb!vm:nwordbits)  (%randomword state))  finally (return result)))  ;;; an UNSIGNEDBYTE of NBITS bits sampled from the Mersenne twister (declaim (inline %randombits)) (defun %randombits (nbits state)  (multiplevaluebind (nfullwords nextrabits)  (floor nbits sb!vm:nwordbits)  (let ((fullchunks (%randomwords nfullwords state)))  (if (zerop nextrabits)  fullchunks  (logior fullchunks  (ash (logand (%randomword state)  (1 (ash 1 nextrabits)))  (* nfullwords  sb!vm:nwordbits)))))))  ;;; 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))  (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 (%randomword state) mask)  inclusivelimit)))  ;;;; outer, dynamicallytyped interface +(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))))) (defun random (arg &optional (state *randomstate*))  (declare (inline %randomsinglefloat  %randomdoublefloat + (declare (inline %randomsinglefloat %randomdoublefloat #!+longfloat %randomlongfloat)) (cond  ((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 (fixnump arg) (<= arg randomfixnummax) (> arg 0)) + (rem (randomchunk state) arg)) ((and (typep arg 'singlefloat) (> arg 0.0f0)) (%randomsinglefloat arg state)) ((and (typep arg 'doublefloat) (> arg 0.0d0)) @@ 381,10 +259,8 @@ #!+longfloat ((and (typep arg 'longfloat) (> arg 0.0l0)) (%randomlongfloat arg state))  ((and (integerp arg) (plusp arg))  (if (= arg (1+ mostpositivefixnum))  (%inclusiverandomfixnum mostpositivefixnum state)  (%inclusiverandominteger (1 arg) state))) + ((and (integerp arg) (> arg 0)) + (%randominteger arg state)) (t (error 'simpletypeerror :expectedtype '(or (integer 1) (float (0))) :datum arg 