# HG changeset patch # User Steve Losh # Date 1486389677 0 # Node ID 03240af4df9bace2f710f6fa860b15266cfd30cf # Parent 486f9b2d6055816d1b4c31a15bd7193186f1bb58 Move PCG to its own lib diff -r 486f9b2d6055 -r 03240af4df9b src/random-numbers.lisp --- a/src/random-numbers.lisp Tue Jan 31 18:54:39 2017 +0000 +++ b/src/random-numbers.lisp Mon Feb 06 14:01:17 2017 +0000 @@ -154,147 +154,5 @@ (first (random-weighted-list weights 1))) -;; from cl-utilities - -;; If we're using the SB-ROTATE-BYTE extension, we should inline our -;; call and let SBCL handle optimization from there. -#+sbcl (declaim (inline rotate-byte)) - -(defun rotate-byte (count bytespec integer) - "Rotates a field of bits within INTEGER; specifically, returns an -integer that contains the bits of INTEGER rotated COUNT times -leftwards within the byte specified by BYTESPEC, and elsewhere -contains the bits of INTEGER. See http://www.cliki.net/ROTATE-BYTE" - (declare (optimize (speed 3) (safety 0) (space 0) (debug 1))) - #-sbcl - (let ((size (byte-size bytespec))) - (when (= size 0) - (return-from rotate-byte integer)) - (let ((count (mod count size))) - (flet ((rotate-byte-from-0 (count size integer) - (let ((bytespec (byte size 0))) - (if (> count 0) - (logior (ldb bytespec (ash integer count)) - (ldb bytespec (ash integer (- count size)))) - (logior (ldb bytespec (ash integer count)) - (ldb bytespec (ash integer (+ count size)))))))) - (dpb (rotate-byte-from-0 count size (ldb bytespec integer)) - bytespec - integer)))) - ;; On SBCL, we use the SB-ROTATE-BYTE extension. - #+sbcl (sb-rotate-byte:rotate-byte count bytespec integer)) - - -;;;; PCG -(defstruct (pcg (:constructor make-pcg%)) - (state 0 :type (unsigned-byte 64)) - (increment 0 :type (unsigned-byte 64))) - -(declaim (inline permute-xor-shift permute-rotate advance-state)) - -(defun permute-xor-shift (data) - (declare (optimize (speed 3) (debug 0) (safety 1)) - (type (unsigned-byte 37) data)) - (-<> data - (ash <> -18) - (logxor data <>) - (ldb (byte 32 0) <>))) - -(defun permute-rotate (data selector) - (declare (optimize (speed 3) (debug 0) (safety 1)) - (type (unsigned-byte 32) data) - (type (unsigned-byte 5) selector)) - (sb-rotate-byte:rotate-byte selector - (byte 32 0) - data)) - -(defun advance-state (pcg) - (declare (optimize (speed 3) (debug 0) (safety 1)) - (type pcg pcg)) - (setf (pcg-state pcg) - (mod (+ (* (pcg-state pcg) 6364136223846793005) - (pcg-increment pcg)) - (expt 2 64))) - nil) - - -; uint64_t oldstate = rng->state; -; rng->state = oldstate * 6364136223846793005ULL + rng->inc; -; bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb -; >>27 -; bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb........................... -; assign to uint32 fuckin lol -; .....bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb........................... - -; uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u; -; uint32_t rot = oldstate >> 59u; -; return (xorshifted >> rot) | (xorshifted << ((-rot) & 31)); -(declaim (ftype (function (pcg) (unsigned-byte 32)) pcg)) - -(defun-inlineable pcg-random (pcg) - (declare (optimize (speed 3) (debug 0) (safety 1)) - (type pcg pcg)) - (let* ((state (pcg-state pcg)) - (data (ldb (byte 37 (- 64 37)) state)) - (selector (ldb (byte 5 (- 64 5)) state))) - (advance-state pcg) - (-<> data - (permute-xor-shift <>) - (permute-rotate <> selector)))) - -(defun make-pcg (seed &optional (stream-id 0)) - (let* ((increment (logior 1 (ash stream-id 1))) - (pcg (make-pcg% :state 0 :increment increment))) - (pcg-random pcg) - (incf (pcg-state pcg) seed) - (modf (pcg-state pcg) 64) - (pcg-random pcg) - pcg)) - -(defun pcg-random-bounded (pcg bound) - (declare - (optimize (speed 3) (debug 0) (safety 1)) - (type pcg pcg) - (type (and (unsigned-byte 32) - (integer 1)) - bound) - (inline pcg-random)) - (loop - :with threshold = (mod (expt 2 32) bound) - :for n = (pcg-random pcg) - :when (>= n threshold) - :do (return (values (mod n bound))))) - -(declaim (ftype (function (pcg (integer 1 32)) - (unsigned-byte 32)) pcg-random-bits)) - -(defun-inline pcg-random-bits (pcg count) - (declare (optimize (speed 3) (debug 0) (safety 1)) - (inline pcg-random)) - (ldb (byte count 0) (pcg-random pcg))) - -(defun pcg-random-float (pcg) - (declare (optimize (speed 3) (debug 0) (safety 1)) - (type pcg pcg)) - ; https://en.wikipedia.org/wiki/Single-precision_floating-point_format - ; Singles have 24 bits of precision - (/ (pcg-random-bits pcg 24) - (coerce (expt 2 24) 'single-float))) - -; https://en.wikipedia.org/wiki/Double-precision_floating-point_format -; Doubles have 53 bits of precision -(defun pcg-random-double (pcg) - (declare (optimize (speed 3) (debug 0) (safety 1)) - (type pcg pcg)) - (/ (logior (ash (pcg-random-bits pcg 26) 27) - (pcg-random-bits pcg 27)) - (coerce (expt 2 53) 'double-float))) - -(defun pcg-random-range (pcg min max) - (+ min (pcg-random-bounded pcg (- max min)))) - -(defun pcg-random-range-inclusive (pcg min max) - (+ min (pcg-random-bounded pcg (1+ (- max min))))) - ;;;; Scratch ; (spectral)