# HG changeset patch
# User Steve Losh <steve@stevelosh.com>
# Date 1485888879 0
# Node ID 486f9b2d6055816d1b4c31a15bd7193186f1bb58
# Parent  02235ffe8b98669b349c6292b9fd230c0baa1c93
Implement PCGs

diff -r 02235ffe8b98 -r 486f9b2d6055 src/names.lisp
--- a/src/names.lisp	Thu Jan 26 23:53:51 2017 +0000
+++ b/src/names.lisp	Tue Jan 31 18:54:39 2017 +0000
@@ -85,3 +85,4 @@
   (Lu cia nia)
   (Da niela)
   (Ga brie la))
+
diff -r 02235ffe8b98 -r 486f9b2d6055 src/random-numbers.lisp
--- a/src/random-numbers.lisp	Thu Jan 26 23:53:51 2017 +0000
+++ b/src/random-numbers.lisp	Tue Jan 31 18:54:39 2017 +0000
@@ -154,5 +154,147 @@
   (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)