# HG changeset patch
# User Steve Losh <steve@stevelosh.com>
# Date 1489537399 0
# Node ID 94eb90b8d9f36522d512070b85afacd2f1dffbe7
# Parent  d8537fcd89e54f4983c386553c2f3be7c57d3d11
Add some docstrings

diff -r d8537fcd89e5 -r 94eb90b8d9f3 src/pcg.lisp
--- a/src/pcg.lisp	Tue Mar 14 13:33:53 2017 +0000
+++ b/src/pcg.lisp	Wed Mar 15 00:23:19 2017 +0000
@@ -7,7 +7,6 @@
 
 
 ;;;; Utils ---------------------------------------------------------------------
-
 (defmacro defun-inline (name &body body)
   "Like `defun`, but declaims `name` to be `inline`."
   `(progn
@@ -52,6 +51,9 @@
 (deftype u64 ()
   '(unsigned-byte 64))
 
+(deftype u32 ()
+  '(unsigned-byte 32))
+
 
 ;;;; Data ---------------------------------------------------------------------
 (defstruct (pcg (:constructor actually-make-pcg))
@@ -66,7 +68,7 @@
 (defun-inline permute-xor-shift (data)
   (declare (optimize speed)
            (type (unsigned-byte 37) data))
-  ;; The reference implemtation does this:
+  ;; The reference implementation does this:
   ;;
   ;;   uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
   ;;
@@ -78,7 +80,7 @@
   ;;   SSSSSbbb bbbbbbbb bbbbbbbb bbbbbbbb bbbbbxxx xxxxxxxx xxxxxxxx xxxxxxxx
   ;;
   ;; * S - 5 "selector" bits for the permutation.  Confusingly the xorshift
-  ;;       permutation doesn't use these to select anything, but instead mixes
+  ;;       permutation doesn't use these to SELECT anything, but instead mixes
   ;;       them into the other random bits.  The next permutation uses them to
   ;;       decide how far to rotate.
   ;; * b - 32 decently random bits, which are going to be permuted and used for
@@ -120,26 +122,27 @@
 
 ;;;; Low-Level API ------------------------------------------------------------
 (declaim
-  (ftype (function (pcg) (unsigned-byte 32)) pcg-random%)
+  (ftype (function (pcg) u32) pcg-random%)
 
   ;; 2^32 streams ought to be enough for anybody
-  (ftype (function (u64 (unsigned-byte 32))
+  (ftype (function (u64 u32)
                    pcg)
          make-pcg%)
 
-  (ftype (function (pcg (integer 1 (#.(expt 2 32))))
-                   (unsigned-byte 32))
+  (ftype (function (pcg (integer 1 (#.(expt 2 32)))) u32)
          pcg-random-bounded%)
 
-  (ftype (function (pcg (integer 1 32))
-                   (unsigned-byte 32))
+  (ftype (function (pcg (signed-byte 32) (signed-byte 32))
+                   (signed-byte 32))
+         pcg-random-range%
+         pcg-random-range-inclusive%)
+
+  (ftype (function (pcg (integer 1 32)) u32)
          pcg-random-bits%)
 
   (ftype (function (pcg u64)) pcg-advance% pcg-rewind%)
 
-  (ftype (function (pcg) single-float) pcg-random-float%)
-
-  (ftype (function (pcg) double-float) pcg-random-double%))
+  (ftype (function (pcg) single-float) pcg-random-float%))
 
 
 (defun-inline pcg-random% (pcg)
@@ -189,6 +192,14 @@
         :when (>= n threshold)
         :do (return (values (mod n bound)))))
 
+(defun-inline pcg-random-range% (pcg min max)
+  (declare (optimize speed))
+  (+ min (pcg-random-bounded% pcg (- max min))))
+
+(defun-inline pcg-random-range-inclusive% (pcg min max)
+  (declare (optimize speed))
+  (+ min (pcg-random-bounded% pcg (1+ (- max min)))))
+
 (defun-inline pcg-random-bits% (pcg count)
   "Return a random `(unsigned-byte COUNT)`.
 
@@ -217,6 +228,13 @@
 
 
 (defun-inline pcg-advance% (pcg steps)
+  "Advance the state of `pcg` by `steps` steps.
+
+  This function returns `nil` and is only useful for its side effects.
+
+  This is a low-level function that assumes you are passing in the correct types.
+
+  "
   (declare (optimize speed))
   ;; See "Random Number Generation with Arbitrary Strides", Forrest B. Brown:
   ;; https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/anl-rn-arb-stride.pdf
@@ -239,16 +257,36 @@
           (% (*+ (pcg-state pcg) multiplier-accumulator increment-accumulator)))))
 
 (defun-inline pcg-rewind% (pcg steps)
+  "Rewind the state of `pcg` by `steps` steps.
+
+  This function returns `nil` and is only useful for its side effects.
+
+  This is a low-level function that assumes you are passing in the correct types.
+
+  "
   (declare (optimize speed))
   (when (plusp steps)
     (pcg-advance% pcg (- +limit+ (1- steps)))))
 
 
 ;;;; High-Level API -----------------------------------------------------------
-(defun make-pcg (&key (seed 0) (stream-id 0))
-  (check-types stream-id (unsigned-byte 32)
-               seed u64)
-  (make-pcg% seed stream-id))
+(defun resolve-seed (seed)
+  (if (null seed)
+    (random (expt 2 64) (make-random-state t))
+    seed))
+
+(defun make-pcg (&key (seed nil) (stream-id 0))
+  "Create and return a new PCG.
+
+  If `seed` is `nil`, a fresh random seed will be generated with the
+  implementation's `cl:random` function, as if by:
+
+    (random ... (make-random-state t))
+
+  "
+  (check-types stream-id u32
+               seed (or null u64))
+  (make-pcg% (resolve-seed seed) stream-id))
 
 
 (defparameter *global-pcg* (make-pcg))
@@ -263,68 +301,117 @@
 
 
 (defun pcg-random (pcg)
+  "Return a random `(unsigned-byte 32)`.
+
+  As a side effect, the state of `pcg` will be advanced.
+
+  "
   (check-types pcg pcg-designator)
   (pcg-random% (resolve-pcg pcg)))
 
 (defun pcg-random-float (pcg)
+  "Return a random `single-float` between `0.0` and `1.0`.
+
+  As a side effect, the state of `pcg` will be advanced.
+
+  "
   (check-types pcg pcg-designator)
   (pcg-random-float% (resolve-pcg pcg)))
 
-(defun pcg-random-double (pcg)
-  (check-types pcg pcg-designator)
-  (pcg-random-double% (resolve-pcg pcg)))
 
 (defun pcg-random-bounded (pcg bound)
+  "Return a random integer between `0` (inclusive) and `bound` (exclusive).
+
+  As a side effect, the state of `pcg` will be advanced.
+
+  "
   (check-types pcg pcg-designator
-               bound (and (unsigned-byte 32)
-                          (integer 1)))
+               bound (and u32 (integer 1)))
   (pcg-random-bounded% (resolve-pcg pcg) bound))
 
 (defun pcg-random-range (pcg min max)
+  "Return a random integer between `min` (inclusive) and `max` (exclusive).
+
+  As a side effect, the state of `pcg` will be advanced.
+
+  "
   (check-types pcg pcg-designator
-               min (unsigned-byte 32)
-               max (unsigned-byte 32))
+               min (signed-byte 32)
+               max (signed-byte 32))
   (assert (< min max) (min max))
-  (+ min (pcg-random-bounded% (resolve-pcg pcg) (- max min))))
+  (pcg-random-range% (resolve-pcg pcg) min max))
 
 (defun pcg-random-range-inclusive (pcg min max)
+  "Return a random integer between `min` (inclusive) and `max` (inclusive).
+
+  As a side effect, the state of `pcg` will be advanced.
+
+  "
   (check-types pcg pcg-designator
-               min (unsigned-byte 32)
-               max (unsigned-byte 32))
+               min (signed-byte 32)
+               max (signed-byte 32))
   (assert (<= min max) (min max))
-  (pcg-random-range (resolve-pcg pcg) min (1+ max)))
+  (pcg-random-range-inclusive% (resolve-pcg pcg) min max))
+
 
 (defun pcg-advance (pcg steps)
+  "Advance the state of `pcg` by `steps` steps.
+
+  This function returns `nil` and is only useful for its side effects.
+
+  "
   (check-types pcg pcg-designator
               steps u64)
   (pcg-advance% (resolve-pcg pcg) steps))
 
 (defun pcg-rewind (pcg steps)
+  "Rewind the state of `pcg` by `steps` steps.
+
+  This function returns `nil` and is only useful for its side effects.
+
+  "
   (check-types pcg pcg-designator
                steps u64)
   (pcg-rewind% (resolve-pcg pcg) steps))
 
 
 ;;;; Scratch ------------------------------------------------------------------
-; (defun data (n)
-;   (loop :repeat n :collect (pcg-random-single-float% *p*)))
+;; (defparameter *p* (make-pcg))
+
+;; (defun data (n)
+;;   (loop :repeat n :collect (pcg-random-float% *p*)))
 
-; (losh:gnuplot
-;   (data 10000)
-;   :x #'identity
-;   :y (lambda (i) (/ 1.0 10000))
-;   :line-width 1
-;   :smooth :cumulative)
+;; (losh:gnuplot
+;;   (data 10000)
+;;   :x #'identity
+;;   :y (lambda (i) (/ 1.0 10000))
+;;   :line-width 2
+;;   :smooth :cumulative)
+
+;; (defun slow-advance (pcg n)
+;;   (dotimes (_ n) (pcg-random pcg))
+;;   nil)
 
-; (defun slow-advance (pcg n)
-;   (dotimes (_ n) (pcg-random pcg))
-;   nil)
+;; (defun test (n)
+;;   (let ((a (make-pcg))
+;;         (b (make-pcg)))
+;;     (time (slow-advance a n))
+;;     (time (pcg-advance% b n))
+;;     (assert (= (pcg-random a) (pcg-random b)))))
 
-; (defun test (n)
-;   (let ((a (make-pcg))
-;         (b (make-pcg)))
-;     (pcg-advance% b n)
-;     (assert (= (pcg-random a) (pcg-random b)))))
+;; (losh:gnuplot
+;;   (sort (losh:hash-table-contents
+;;           (losh:proportions
+;;             (loop :repeat 1000000
+;;                   :collect (pcg-random-range-inclusive *p* -50 200))))
+;;         #'<
+;;         :key #'first)
+;;   :x #'first
+;;   :y #'second
+;;   :min-y 0.00
+;;   :max-y 0.01
+;;   :style :lines
+;;   :line-width 1)
 
 ;;;; TODO ----------------------------------------------------------------------
 ;;; https://experilous.com/1/blog/post/perfect-fast-random-floating-point-numbers