src/code/bignum-random.lisp

   1 ;;;; generation of random bignums
   2 ;;;;
   3 ;;;; The implementation assumes that the random chunk size is either
   4 ;;;; equal to the word size or equal to half the word size.
   5
   6 ;;;; This software is part of the SBCL system. See the README file for
   7 ;;;; more information.
   8 ;;;;
   9 ;;;; This software is derived from the CMU CL system, which was
  10 ;;;; written at Carnegie Mellon University and released into the
  11 ;;;; public domain. The software is in the public domain and is
  12 ;;;; provided with absolutely no warranty. See the COPYING and CREDITS
  13 ;;;; files for more information.
  14
  15 (in-package "SB!BIGNUM")
  16
  17 ;;; Return T if the least significant N-BITS bits of BIGNUM are all
  18 ;;; zero, else NIL. If the integer-length of BIGNUM is less than N-BITS,
  19 ;;; the result is NIL, too.
  20 (declaim (inline bignum-lower-bits-zero-p))
  21 (defun bignum-lower-bits-zero-p (bignum n-bits)
  22   (declare (type bignum-type bignum)
  23            (type bit-index n-bits))
  24   (multiple-value-bind (n-full-digits n-bits-partial-digit)
  25       (floor n-bits digit-size)
  26     (declare (type bignum-index n-full-digits))
  27     (when (> (%bignum-length bignum) n-full-digits)
  28       (dotimes (index n-full-digits)
  29         (declare (type bignum-index index))
  30         (unless (zerop (%bignum-ref bignum index))
  31           (return-from bignum-lower-bits-zero-p nil)))
  32       (zerop (logand (1- (ash 1 n-bits-partial-digit))
  33                      (%bignum-ref bignum n-full-digits))))))
  34
  35 ;;; Return a nonnegative integer of DIGIT-SIZE many pseudo random bits.
  36 (declaim (inline random-bignum-digit))
  37 (defun random-bignum-digit (state)
  38   (if (= n-random-chunk-bits digit-size)
  39       (random-chunk state)
  40       (big-random-chunk state)))
  41
  42 ;;; Return a nonnegative integer of N-BITS many pseudo random bits.
  43 ;;; N-BITS must be nonnegative and less than DIGIT-SIZE.
  44 (declaim (inline random-bignum-partial-digit))
  45 (defun random-bignum-partial-digit (n-bits state)
  46   (declare (type (integer 0 #.(1- digit-size)) n-bits)
  47            (type random-state state))
  48   (logand (1- (ash 1 n-bits))
  49           (if (<= n-bits n-random-chunk-bits)
  50               (random-chunk state)
  51               (big-random-chunk state))))
  52
  53 ;;; Create a (nonnegative) bignum by concatenating RANDOM-CHUNK and
  54 ;;; BIT-COUNT many pseudo random bits, normalise and return it.
  55 ;;; RANDOM-CHUNK must fit into a bignum digit.
  56 (declaim (inline concatenate-random-bignum))
  57 (defun concatenate-random-bignum (random-chunk bit-count state)
  58   (declare (type bignum-element-type random-chunk)
  59            (type (integer 0 #.sb!xc:most-positive-fixnum) bit-count)
  60            (type random-state state))
  61   (let* ((n-total-bits (+ 1 n-random-chunk-bits bit-count)) ; sign bit
  62          (length (ceiling n-total-bits digit-size))
  63          (bignum (%allocate-bignum length)))
  64     (multiple-value-bind (n-random-digits n-random-bits)
  65         (floor bit-count digit-size)
  66       (declare (type bignum-index n-random-digits))
  67       (dotimes (index n-random-digits)
  68         (setf (%bignum-ref bignum index)
  69               (random-bignum-digit state)))
  70       (if (zerop n-random-bits)
  71           (setf (%bignum-ref bignum n-random-digits) random-chunk)
  72           (progn
  73             (setf (%bignum-ref bignum n-random-digits)
  74                   (%logior (random-bignum-partial-digit n-random-bits
  75                                                         state)
  76                            (%ashl random-chunk n-random-bits)))
  77             (let ((shift (- digit-size n-random-bits)))
  78               (when (< shift n-random-chunk-bits)
  79                 (setf (%bignum-ref bignum (1+ n-random-digits))
  80                       (%digit-logical-shift-right random-chunk shift)))))))
  81     (%normalize-bignum bignum length)))
  82
  83 ;;; Create and return a (nonnegative) bignum of N-BITS many pseudo
  84 ;;; random bits. The result is normalised, so may be a fixnum, too.
  85 (declaim (inline make-random-bignum))
  86 (defun make-random-bignum (n-bits state)
  87   (declare (type (and fixnum (integer 0)) n-bits)
  88            (type random-state state))
  89   (let* ((n-total-bits (1+ n-bits)) ; sign bit
  90          (length (ceiling n-total-bits digit-size))
  91          (bignum (%allocate-bignum length)))
  92     (declare (type bignum-index length))
  93     (multiple-value-bind (n-digits n-bits-partial-digit)
  94         (floor n-bits digit-size)
  95       (declare (type bignum-index n-digits))
  96       (dotimes (index n-digits)
  97         (setf (%bignum-ref bignum index)
  98               (random-bignum-digit state)))
  99       (unless (zerop n-bits-partial-digit)
 100         (setf (%bignum-ref bignum n-digits)
 101               (random-bignum-partial-digit n-bits-partial-digit state))))
 102     (%normalize-bignum bignum length)))
 103
 104 ;;; Create and return a pseudo random bignum less than ARG. The result
 105 ;;; is normalised, so may be a fixnum, too. We try to keep the number of
 106 ;;; times RANDOM-CHUNK is called and the amount of storage consed to a
 107 ;;; minimum.
 108 ;;; Four cases are differentiated:
 109 ;;; * If ARG is a power of two and only one random chunk is needed to
 110 ;;;   supply a sufficient number of bits, a chunk is generated and
 111 ;;;   shifted to get the correct number of bits. This only conses if the
 112 ;;;   result is indeed a bignum. This case can only occur if the size of
 113 ;;;   the random chunks is equal to the word size.
 114 ;;; * If ARG is a power of two and multiple chunks are needed, we call
 115 ;;;   MAKE-RANDOM-BIGNUM. Here a bignum is always consed even if it
 116 ;;;   happens to normalize to a fixnum, which can't be avoided.
 117 ;;; * If ARG is not a power of two but one random chunk suffices an
 118 ;;;   accept-reject loop is used. Each time through the loop a chunk is
 119 ;;;   generated and shifted to get the correct number of bits. This only
 120 ;;;   conses if the final accepted result is indeed a bignum. This case
 121 ;;;   too can only occur if the size of the random chunks is equal to the
 122 ;;;   word size.
 123 ;;; * If ARG is not a power of two and multiple chunks are needed an
 124 ;;;   accept-reject loop is used that detects rejection early by
 125 ;;;   starting the generation with a random chunk aligned to the most
 126 ;;;   significant bits of ARG. If the random value is larger than the
 127 ;;;   corresponding chunk of ARG we don't need to generate the full
 128 ;;;   amount of random bits but can retry immediately. If the random
 129 ;;;   value is smaller than the ARG chunk we know already that the
 130 ;;;   result will be accepted independently of what the remaining random
 131 ;;;   bits will be, so we generate them and return. Only in the rare
 132 ;;;   case that the random value and the ARG chunk are equal we need to
 133 ;;;   generate and compare the complete random number and risk to reject
 134 ;;;   it.
 135 (defun %random-bignum (arg state)
 136   (declare (type (integer #.(1+ sb!xc:most-positive-fixnum)) arg)
 137            (type random-state state)
 138            (inline bignum-lower-bits-zero-p))
 139   (let ((n-bits (bignum-integer-length arg)))
 140     (declare (type (integer #.sb!vm:n-fixnum-bits) n-bits))
 141     ;; Don't use (ZEROP (LOGAND ARG (1- ARG))) to test if ARG is a power
 142     ;; of two as that would cons.
 143     (cond ((bignum-lower-bits-zero-p arg (1- n-bits))
 144            ;; ARG is a power of two. We need one bit less than its
 145            ;; INTEGER-LENGTH. Not using (DECF N-BITS) here allows the
 146            ;; compiler to make optimal use of the type declaration for
 147            ;; N-BITS above.
 148            (let ((n-bits (1- n-bits)))
 149              (if (<= n-bits n-random-chunk-bits)
 150                  (%digit-logical-shift-right (random-chunk state)
 151                                              (- n-random-chunk-bits n-bits))
 152                  (make-random-bignum n-bits state))))
 153           ((<= n-bits n-random-chunk-bits)
 154            (let ((shift (- n-random-chunk-bits n-bits))
 155                  (arg (%bignum-ref arg 0)))
 156              (loop
 157                (let ((bits (%digit-logical-shift-right (random-chunk state)
 158                                                        shift)))
 159                  (when (< bits arg)
 160                    (return bits))))))
 161           (t
 162            ;; ARG is not a power of two and we need more than one random
 163            ;; chunk.
 164            (let* ((shift (- n-bits n-random-chunk-bits))
 165                   (arg-first-chunk (ldb (byte n-random-chunk-bits shift)
 166                                         arg)))
 167              (loop
 168                (let ((random-chunk (random-chunk state)))
 169                  ;; If the random value is larger than the corresponding
 170                  ;; chunk from the most significant bits of ARG we can
 171                  ;; retry immediately; no need to generate the remaining
 172                  ;; random bits.
 173                  (unless (> random-chunk arg-first-chunk)
 174                    ;; We need to generate the complete random number.
 175                    (let ((bits (concatenate-random-bignum random-chunk
 176                                                           shift state)))
 177                      ;; While the second comparison below subsumes the
 178                      ;; first, the first is faster and will nearly
 179                      ;; always be true, so it's worth it to try it
 180                      ;; first.
 181                      (when (or (< random-chunk arg-first-chunk)
 182                                (< bits arg))
 183                        (return bits)))))))))))