03240af4df9b

Move PCG to its own lib
[view raw] [browse files]
author Steve Losh <steve@stevelosh.com>
date Mon, 06 Feb 2017 14:01:17 +0000
parents 486f9b2d6055
children 6b1e5154d5de
branches/tags (none)
files src/random-numbers.lisp

Changes

--- 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)