From: Lutz Euler <leuler@us...>  20120501 17:00:24

The branch "master" has been updated in SBCL: via 436b2ab0276f547e8537b6c1fb52b11fa1f53975 (commit) from fa1f8141814d146ed69630dcd08a749058ef5119 (commit)  Log  commit 436b2ab0276f547e8537b6c1fb52b11fa1f53975 Author: Lutz Euler <lutz.euler@...> Date: Tue May 1 18:59:12 2012 +0200 Better equidistributed and faster/less consing integer RANDOM. Up to now the implementation of RANDOM with an integer argument just generated a few more random bits than the length of the argument and took this value MOD the argument. This led to a slightly uneven distribution of the possible values unless the argument was a power of two. Moreover, for bignums, the algorithm was quadratic both in time and space dependent on the number of bits of the argument. Instead generate random integers using an acceptreject loop and change the bignum implementation to an algorithm that is linear in time and space. I took some inspiration from WHN's attempt at an acceptreject loop implementation in commit 0a7604d54581d2c846838c26ce6a7993629586fa and following. Thanks to Christophe Rhodes for reviewing this patch! Some details: The implementation works correctly with both a random chunk size equal to the word size and equal to half the word size. This is currently necessary as a 32bit pseudo random generator is used both under 32 and under 64 bit word size. In the generic RANDOM case, fixnum and bignum limits are differentiated: With a fixnum limit an acceptreject loop on a masked random chunk is always used. Under 64 bit word size two random chunks are used only if the limit is so large that one doesn't suffice. This never conses. With a bignum limit four cases are differentiated to minimize consing. If just one random chunk is needed to supply a sufficient number of bits the implementation only conses if the result is indeed a bignum: * If the limit is a power of two, a chunk is generated and shifted to get the correct number of bits. * If the limit is not a power of two an acceptreject loop with shifting is used. If more than one random chunk is needed, a bignum is always consed even if it happens to normalize to a fixnum: * If the limit is a power of two a straightforward algorithm is used to fill a newly created bignum with random bits. * If the limit is not a power of two an acceptreject loop is used that detects rejection early by starting from the most significant bits, thus generating on the average only one random chunk more than needed to fill the result once. The test for power of two is nonconsing, too. In the case of a compiletime constant integer argument (of at most word size) a DEFTRANSFORM triggers, that, in the general case, compiles an acceptreject loop. For values of the limit where this sufficiently reduces the rejection probability the largest multiple of the limit fitting in one or two random chunks is used instead inside the loop. To bring the result in the correct range a division is then necessary (which the compiler converts into a multiplication). Powers of two are optimized by leaving out the rejection test. In those cases where a word has more bits than a random chunk, the generated expression uses two chunks only if necessary.  NEWS  13 +++ buildorder.lispexpr  2 + packagedatalist.lispexpr  5 + src/code/bignumrandom.lisp  183 ++++++++++++++++++++++++++++++++++++++++++ src/code/random.lisp  14 + src/code/targetrandom.lisp  65 ++++++++++ src/compiler/floattran.lisp  124 ++++++++++++++++ tests/random.pure.lisp  4 + 8 files changed, 318 insertions(+), 92 deletions() diff git a/NEWS b/NEWS index 86d8250..5b4ac62 100644  a/NEWS +++ b/NEWS @@ 1,5 +1,18 @@ ;;;; * coding: utf8; fillcolumn: 78 * changes relative to sbcl1.0.56: + * RANDOM enhancements and bug fixes: + ** bug fix: the range and distribution of random integers could be + catastrophically wrong when the compiler derived the type of its + argument as a disjoint set of small integers. + ** bug fix: the distribution of random integers is now completely + uniform even when the specified limit is not a power of two. + (Previously some values could be about 0.1 % more probable than + others in the worst case.) + ** RANDOM on large integer arguments is generally faster and conses + less than before; this is visible for fixnums above a length of + about 24 bits, but extremely so for bignums: the old implementation + used time and space quadratical in the size of the argument there, + the new one is linear. * enhancement: redesigned protocol for quitting SBCL. SBEXT:EXIT is the new main entry point, SBEXT:QUIT is deprecated. * enhancement: additions to the SBTHREAD API: RETURNFROMTHREAD, diff git a/buildorder.lispexpr b/buildorder.lispexpr index e2ec49c..c31d518 100644  a/buildorder.lispexpr +++ b/buildorder.lispexpr @@ 656,6 +656,8 @@ ("src/code/targetsap" :nothost) ; uses SAPINT type ("src/code/targetpackage" :nothost) ; needs "code/package" ("src/code/targetrandom" :nothost) ; needs "code/random" + ("src/code/bignumrandom" :nothost) ; needs "code/random" and + ; "code/bignum" ("src/code/targethashtable" :nothost) ; needs "code/hashtable" ("src/code/reader" :nothost) ; needs "code/readtable" ("src/code/targetstream" :nothost) ; needs WHITESPACEP from "code/reader" diff git a/packagedatalist.lispexpr b/packagedatalist.lispexpr index cbc0b4a..13ede30 100644  a/packagedatalist.lispexpr +++ b/packagedatalist.lispexpr @@ 180,7 +180,7 @@ of SBCL which maintained the CMUCLstyle split into two packages.)" "FLOATBIGNUMRATIO" "MAKESMALLBIGNUM" "MULTIPLYBIGNUMANDFIXNUM" "MULTIPLYBIGNUMS" "MULTIPLYFIXNUMS" "NEGATEBIGNUM"  "SUBTRACTBIGNUM" "SXHASHBIGNUM")) + "%RANDOMBIGNUM" "SUBTRACTBIGNUM" "SXHASHBIGNUM")) #s(sbcold:packagedata :name "SB!C" @@ 1865,6 +1865,7 @@ is a good idea, but see SBSYS re. blurring of boundaries." #!+longfloat "%RANDOMLONGFLOAT" "%RANDOMSINGLEFLOAT" "STATICCLASSOID" "%FUNCALLABLEINSTANCEINFO" "RANDOMCHUNK" "BIGRANDOMCHUNK" + "NRANDOMCHUNKBITS" "LAYOUTCLOSHASHLIMIT" "BUILTINCLASSOIDDIRECTSUPERCLASSES" "BUILTINCLASSOIDTRANSLATION" "RANDOMLAYOUTCLOSHASH" @@ 1873,7 +1874,7 @@ is a good idea, but see SBSYS re. blurring of boundaries." "%SETFUNCALLABLEINSTANCELAYOUT" "BASICSTRUCTURECLASSOID" "REGISTERLAYOUT"  "FUNCALLABLEINSTANCE" "RANDOMFIXNUMMAX" + "FUNCALLABLEINSTANCE" "MAKESTATICCLASSOID" "%MAKESYMBOL" "%FUNCALLABLEINSTANCEFUNCTION" "SYMBOLHASH" diff git a/src/code/bignumrandom.lisp b/src/code/bignumrandom.lisp new file mode 100644 index 0000000..821d4a7  /dev/null +++ b/src/code/bignumrandom.lisp @@ 0,0 +1,183 @@ +;;;; generation of random bignums +;;;; +;;;; The implementation assumes that the random chunk size is either +;;;; equal to the word size or equal to half the word size. + +;;;; This software is part of the SBCL system. See the README file for +;;;; more information. +;;;; +;;;; This software is derived from the CMU CL system, which was +;;;; written at Carnegie Mellon University and released into the +;;;; public domain. The software is in the public domain and is +;;;; provided with absolutely no warranty. See the COPYING and CREDITS +;;;; files for more information. + +(inpackage "SB!BIGNUM") + +;;; Return T if the least significant NBITS bits of BIGNUM are all +;;; zero, else NIL. If the integerlength of BIGNUM is less than NBITS, +;;; the result is NIL, too. +(declaim (inline bignumlowerbitszerop)) +(defun bignumlowerbitszerop (bignum nbits) + (declare (type bignumtype bignum) + (type bitindex nbits)) + (multiplevaluebind (nfulldigits nbitspartialdigit) + (floor nbits digitsize) + (declare (type bignumindex nfulldigits)) + (when (> (%bignumlength bignum) nfulldigits) + (dotimes (index nfulldigits) + (declare (type bignumindex index)) + (unless (zerop (%bignumref bignum index)) + (returnfrom bignumlowerbitszerop nil))) + (zerop (logand (1 (ash 1 nbitspartialdigit)) + (%bignumref bignum nfulldigits)))))) + +;;; Return a nonnegative integer of DIGITSIZE many pseudo random bits. +(declaim (inline randombignumdigit)) +(defun randombignumdigit (state) + (if (= nrandomchunkbits digitsize) + (randomchunk state) + (bigrandomchunk state))) + +;;; Return a nonnegative integer of NBITS many pseudo random bits. +;;; NBITS must be nonnegative and less than DIGITSIZE. +(declaim (inline randombignumpartialdigit)) +(defun randombignumpartialdigit (nbits state) + (declare (type (integer 0 #.(1 digitsize)) nbits) + (type randomstate state)) + (logand (1 (ash 1 nbits)) + (if (<= nbits nrandomchunkbits) + (randomchunk state) + (bigrandomchunk state)))) + +;;; Create a (nonnegative) bignum by concatenating RANDOMCHUNK and +;;; BITCOUNT many pseudo random bits, normalise and return it. +;;; RANDOMCHUNK must fit into a bignum digit. +(declaim (inline concatenaterandombignum)) +(defun concatenaterandombignum (randomchunk bitcount state) + (declare (type bignumelementtype randomchunk) + (type (integer 0 #.sb!xc:mostpositivefixnum) bitcount) + (type randomstate state)) + (let* ((ntotalbits (+ 1 nrandomchunkbits bitcount)) ; sign bit + (length (ceiling ntotalbits digitsize)) + (bignum (%allocatebignum length))) + (multiplevaluebind (nrandomdigits nrandombits) + (floor bitcount digitsize) + (declare (type bignumindex nrandomdigits)) + (dotimes (index nrandomdigits) + (setf (%bignumref bignum index) + (randombignumdigit state))) + (if (zerop nrandombits) + (setf (%bignumref bignum nrandomdigits) randomchunk) + (progn + (setf (%bignumref bignum nrandomdigits) + (%logior (randombignumpartialdigit nrandombits + state) + (%ashl randomchunk nrandombits))) + (let ((shift ( digitsize nrandombits))) + (when (< shift nrandomchunkbits) + (setf (%bignumref bignum (1+ nrandomdigits)) + (%digitlogicalshiftright randomchunk shift))))))) + (%normalizebignum bignum length))) + +;;; Create and return a (nonnegative) bignum of NBITS many pseudo +;;; random bits. The result is normalised, so may be a fixnum, too. +(declaim (inline makerandombignum)) +(defun makerandombignum (nbits state) + (declare (type (and fixnum (integer 0)) nbits) + (type randomstate state)) + (let* ((ntotalbits (1+ nbits)) ; sign bit + (length (ceiling ntotalbits digitsize)) + (bignum (%allocatebignum length))) + (declare (type bignumindex length)) + (multiplevaluebind (ndigits nbitspartialdigit) + (floor nbits digitsize) + (declare (type bignumindex ndigits)) + (dotimes (index ndigits) + (setf (%bignumref bignum index) + (randombignumdigit state))) + (unless (zerop nbitspartialdigit) + (setf (%bignumref bignum ndigits) + (randombignumpartialdigit nbitspartialdigit state)))) + (%normalizebignum bignum length))) + +;;; Create and return a pseudo random bignum less than ARG. The result +;;; is normalised, so may be a fixnum, too. We try to keep the number of +;;; times RANDOMCHUNK is called and the amount of storage consed to a +;;; minimum. +;;; Four cases are differentiated: +;;; * If ARG is a power of two and only one random chunk is needed to +;;; supply a sufficient number of bits, a chunk is generated and +;;; shifted to get the correct number of bits. This only conses if the +;;; result is indeed a bignum. This case can only occur if the size of +;;; the random chunks is equal to the word size. +;;; * If ARG is a power of two and multiple chunks are needed, we call +;;; MAKERANDOMBIGNUM. Here a bignum is always consed even if it +;;; happens to normalize to a fixnum, which can't be avoided. +;;; * If ARG is not a power of two but one random chunk suffices an +;;; acceptreject loop is used. Each time through the loop a chunk is +;;; generated and shifted to get the correct number of bits. This only +;;; conses if the final accepted result is indeed a bignum. This case +;;; too can only occur if the size of the random chunks is equal to the +;;; word size. +;;; * If ARG is not a power of two and multiple chunks are needed an +;;; acceptreject loop is used that detects rejection early by +;;; starting the generation with a random chunk aligned to the most +;;; significant bits of ARG. If the random value is larger than the +;;; corresponding chunk of ARG we don't need to generate the full +;;; amount of random bits but can retry immediately. If the random +;;; value is smaller than the ARG chunk we know already that the +;;; result will be accepted independently of what the remaining random +;;; bits will be, so we generate them and return. Only in the rare +;;; case that the random value and the ARG chunk are equal we need to +;;; generate and compare the complete random number and risk to reject +;;; it. +(defun %randombignum (arg state) + (declare (type (integer #.(1+ sb!xc:mostpositivefixnum)) arg) + (type randomstate state) + (inline bignumlowerbitszerop)) + (let ((nbits (bignumintegerlength arg))) + (declare (type (integer #.sb!vm:nfixnumbits) nbits)) + ;; Don't use (ZEROP (LOGAND ARG (1 ARG))) to test if ARG is a power + ;; of two as that would cons. + (cond ((bignumlowerbitszerop arg (1 nbits)) + ;; ARG is a power of two. We need one bit less than its + ;; INTEGERLENGTH. Not using (DECF NBITS) here allows the + ;; compiler to make optimal use of the type declaration for + ;; NBITS above. + (let ((nbits (1 nbits))) + (if (<= nbits nrandomchunkbits) + (%digitlogicalshiftright (randomchunk state) + ( nrandomchunkbits nbits)) + (makerandombignum nbits state)))) + ((<= nbits nrandomchunkbits) + (let ((shift ( nrandomchunkbits nbits)) + (arg (%bignumref arg 0))) + (loop + (let ((bits (%digitlogicalshiftright (randomchunk state) + shift))) + (when (< bits arg) + (return bits)))))) + (t + ;; ARG is not a power of two and we need more than one random + ;; chunk. + (let* ((shift ( nbits nrandomchunkbits)) + (argfirstchunk (ldb (byte nrandomchunkbits shift) + arg))) + (loop + (let ((randomchunk (randomchunk state))) + ;; If the random value is larger than the corresponding + ;; chunk from the most significant bits of ARG we can + ;; retry immediately; no need to generate the remaining + ;; random bits. + (unless (> randomchunk argfirstchunk) + ;; We need to generate the complete random number. + (let ((bits (concatenaterandombignum randomchunk + shift state))) + ;; While the second comparison below subsumes the + ;; first, the first is faster and will nearly + ;; always be true, so it's worth it to try it + ;; first. + (when (or (< randomchunk argfirstchunk) + (< bits arg)) + (return bits))))))))))) diff git a/src/code/random.lisp b/src/code/random.lisp index 373a3b2..cfcc7fc 100644  a/src/code/random.lisp +++ b/src/code/random.lisp @@ 10,19 +10,7 @@ (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)))) +(def!constant nrandomchunkbits 32) (sb!xc:defstruct (randomstate (:constructor %makerandomstate) (:copier nil)) ; since shallow copy is wrong diff git a/src/code/targetrandom.lisp b/src/code/targetrandom.lisp index 57c57a7..b4cb278 100644  a/src/code/targetrandom.lisp +++ b/src/code/targetrandom.lisp @@ 292,9 +292,8 @@ #!sbfluid (declaim (inline bigrandomchunk)) (defun bigrandomchunk (state) (declare (type randomstate state))  (logand (1 (expt 2 64))  (logior (ash (randomchunk state) 32)  (randomchunk state)))) + (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 @@ 310,7 +309,7 @@ (* arg ( (makesinglefloat (dpb (ash (randomchunk state)  ( sb!vm:singlefloatdigits randomchunklength)) + ( sb!vm:singlefloatdigits nrandomchunkbits)) sb!vm:singlefloatsignificandbyte (singlefloatbits 1.0))) 1.0))) @@ 333,7 +332,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)) @@ 348,7 +347,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)) @@ 356,26 +355,46 @@ 1d0)))) ;;;; random integers +;;;; random fixnums (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))))) +;;; Generate and return a pseudo random fixnum less than ARG. To achieve +;;; equidistribution an acceptreject loop is used. +;;; No extra effort is made to detect the case of ARG being a power of +;;; two where rejection is not possible, as the cost of checking for +;;; this case is the same as doing the rejection test. When ARG is +;;; larger than (expt 2 NRANDOMCHUNKBITS), which can only happen if +;;; the random chunk size is half the word size, two random chunks are +;;; used in each loop iteration, otherwise only one. Finally, the +;;; rejection probability could often be reduced by not masking the +;;; chunk but rejecting only values as least as large as the largest +;;; multiple of ARG that fits in a chunk (or two), but this is not done +;;; as the speed gains due to needing fewer loop iterations are by far +;;; outweighted by the cost of the two divisions required (one to find +;;; the multiplier and one to bring the result into the correct range). +#!sbfluid (declaim (inline %randomfixnum)) +(defun %randomfixnum (arg state) + (declare (type (integer 1 #.sb!xc:mostpositivefixnum) arg) + (type randomstate state)) + (if (= arg 1) + 0 + (let* ((nbits (integerlength (1 arg))) + (mask (1 (ash 1 nbits)))) + (macrolet ((acceptrejectloop (generator) + `(loop + (let ((bits (logand mask (,generator state)))) + (when (< bits arg) + (return bits)))))) + (assert (<= nbits (* 2 nrandomchunkbits))) + (if (<= nbits nrandomchunkbits) + (acceptrejectloop randomchunk) + (acceptrejectloop bigrandomchunk)))))) (defun random (arg &optional (state *randomstate*))  (declare (inline %randomsinglefloat %randomdoublefloat + (declare (inline %randomfixnum %randomsinglefloat %randomdoublefloat #!+longfloat %randomlongfloat)) (cond  ((and (fixnump arg) (<= arg randomfixnummax) (> arg 0))  (rem (randomchunk state) arg)) + ((and (fixnump arg) (> arg 0)) + (%randomfixnum arg state)) ((and (typep arg 'singlefloat) (> arg 0.0f0)) (%randomsinglefloat arg state)) ((and (typep arg 'doublefloat) (> arg 0.0d0)) @@ 383,8 +402,8 @@ #!+longfloat ((and (typep arg 'longfloat) (> arg 0.0l0)) (%randomlongfloat arg state))  ((and (integerp arg) (> arg 0))  (%randominteger arg state)) + ((and (bignump arg) (> arg 0)) + (%randombignum arg state)) (t (error 'simpletypeerror :expectedtype '(or (integer 1) (float (0))) :datum arg diff git a/src/compiler/floattran.lisp b/src/compiler/floattran.lisp index 8d2eaed..b5cd536 100644  a/src/compiler/floattran.lisp +++ b/src/compiler/floattran.lisp @@ 46,59 +46,79 @@ (frob %randomsinglefloat singlefloat) (frob %randomdoublefloat doublefloat)) ;;; 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. +;;; Return an expression to generate an integer of NBITS many random +;;; bits, using the minimal number of random chunks possible. +(defun generaterandomexprforpowerof2 (nbits state) + (declare (type (integer 1 #.sb!vm:nwordbits) nbits)) + (multiplevaluebind (nchunkbits chunkexpr) + (cond ((<= nbits nrandomchunkbits) + (values nrandomchunkbits `(randomchunk ,state))) + ((<= nbits (* 2 nrandomchunkbits)) + (values (* 2 nrandomchunkbits) `(bigrandomchunk ,state))) + (t + (error "Unexpectedly small NRANDOMCHUNKBITS"))) + (if (< nbits nchunkbits) + `(logand ,(1 (ash 1 nbits)) ,chunkexpr) + chunkexpr))) + +;;; This transform for compiletime constant wordsized integers +;;; generates an acceptreject loop to achieve equidistribution of the +;;; returned values. Several optimizations are done: If NUM is a power +;;; of two no loop is needed. If the random chunk size is half the word +;;; size only one chunk is used where sufficient. For values of NUM +;;; where it is possible and results in faster code, the rejection +;;; probability is reduced by accepting all values below the largest +;;; multiple of the limit that fits into one or two chunks and and doing +;;; a division to get the random value into the desired range. (deftransform random ((num &optional state)  ((integer 1 #.(expt 2 sb!vm::nwordbits)) &optional *))  ;; FIXME: I almost conditionalized this as #!+sbdoc. Find some way  ;; of automatically finding #!+sbdoc in proximity to DEFTRANSFORM  ;; to let me scan for places that I made this mistake and didn't  ;; catch myself.  "use inline (UNSIGNEDBYTE 32) operations"  (let ((type (lvartype num))  (limit (expt 2 sb!vm::nwordbits))  (randomchunk (ecase sb!vm::nwordbits  (32 'randomchunk)  (64 'sb!kernel::bigrandomchunk))))  (if (numerictypep type)  (let ((numhigh (numerictypehigh (lvartype num))))  (aver numhigh)  (cond ((constantlvarp num)  ;; Check the worst case sum absolute error for the  ;; random number expectations.  (let ((rem (rem limit numhigh)))  (unless (< (/ (* 2 rem ( numhigh rem))  numhigh limit)  (expt 2 ( sb!kernel::randomintegerextrabits)))  (giveupir1transform  "The random number expectations are inaccurate."))  (if (= numhigh limit)  `(,randomchunk (or state *randomstate*))  #!(or x86 x8664)  `(rem (,randomchunk (or state *randomstate*)) num)  #!+(or x86 x8664)  ;; Use multiplication, which is faster.  `(values (sb!bignum::%multiply  (,randomchunk (or state *randomstate*))  num)))))  ((> numhigh randomfixnummax)  (giveupir1transform  "The range is too large to ensure an accurate result."))  #!+(or x86 x8664)  ((< numhigh limit)  `(values (sb!bignum::%multiply  (,randomchunk (or state *randomstate*))  num)))  (t  `(rem (,randomchunk (or state *randomstate*)) num))))  ;; KLUDGE: a relatively conservative treatment, but better  ;; than a bug (reported by PFD sbcldevel towards the end of  ;; 200411.  (giveupir1transform  "Argument type is too complex to optimize for.")))) + ((constantarg (integer 1 #.(expt 2 sb!vm:nwordbits))) + &optional *) + * + :policy (and (> speed compilationspeed) + (> speed space))) + "optimize to inlined RANDOMCHUNK operations" + (let ((num (lvarvalue num))) + (if (= num 1) + 0 + (flet ((chunknbitsandexpr (nbits) + (cond ((<= nbits nrandomchunkbits) + (values nrandomchunkbits + '(randomchunk (or state *randomstate*)))) + ((<= nbits (* 2 nrandomchunkbits)) + (values (* 2 nrandomchunkbits) + '(bigrandomchunk (or state *randomstate*)))) + (t + (error "Unexpectedly small NRANDOMCHUNKBITS"))))) + (if (zerop (logand num (1 num))) + ;; NUM is a power of 2. + (let ((nbits (integerlength (1 num)))) + (multiplevaluebind (nchunkbits chunkexpr) + (chunknbitsandexpr nbits) + (if (< nbits nchunkbits) + `(logand ,(1 (ash 1 nbits)) ,chunkexpr) + chunkexpr))) + ;; Generate an acceptreject loop. + (let ((nbits (integerlength num))) + (multiplevaluebind (nchunkbits chunkexpr) + (chunknbitsandexpr nbits) + (if (or (> (* num 3) (expt 2 nchunkbits)) + (logbitp ( nbits 2) num)) + ;; Division can't help as the quotient is below 3, + ;; or is too costly as the rejection probability + ;; without it is already small (namely at most 1/4 + ;; with the given test, which is experimentally a + ;; reasonable threshold and cheap to test for). + `(loop + (let ((bits ,(generaterandomexprforpowerof2 + nbits '(or state *randomstate*)))) + (when (< bits num) + (return bits)))) + (let ((d (truncate (expt 2 nchunkbits) num))) + `(loop + (let ((bits ,chunkexpr)) + (when (< bits ,(* num d)) + (return (values (truncate bits ,d))))))))))))))) + ;;;; float accessors diff git a/tests/random.pure.lisp b/tests/random.pure.lisp index 101333e..829ff4d 100644  a/tests/random.pure.lisp +++ b/tests/random.pure.lisp @@ 114,8 +114,8 @@ nconc (list (1 (expt 2 i)) (expt 2 i) (1+ (expt 2 i))))  ,@(loop for i from (1 sbkernel::randomchunklength)  to (* sbkernel::randomchunklength 4) + ,@(loop for i from (1 sbkernel::nrandomchunkbits) + to (* sbkernel::nrandomchunkbits 4) collect (* 3 (expt 2 i))) ,@(loop for i from 2 to sbvm:nwordbits for n = (expt 16 i)  hooks/postreceive  SBCL 