c6a3ab353556

Initial commit
[view raw] [browse files]
author Steve Losh <steve@stevelosh.com>
date Fri, 10 Feb 2017 23:31:25 +0000
parents
children 7781bd17fcdc
branches/tags (none)
files .ffignore .hgignore .lispwords Makefile README.markdown cl-pcg.asd docs/01-installation.markdown docs/02-usage.markdown docs/03-reference.markdown docs/04-changelog.markdown docs/api.lisp docs/footer.markdown docs/index.markdown docs/title package.lisp src/pcg.fasl src/pcg.lisp vendor/make-quickutils.lisp vendor/quickutils-package.lisp vendor/quickutils.lisp

Changes

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.ffignore	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,1 @@
+docs/build
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.hgignore	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,4 @@
+syntax: glob
+
+scratch.lisp
+docs/build
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.lispwords	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,1 @@
+(1 bind bind*)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Makefile	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,27 @@
+.PHONY: vendor pubdocs
+
+sourcefiles = $(shell ffind --full-path --literal .lisp)
+docfiles = $(shell ls docs/*.markdown)
+apidocs = $(shell ls docs/*reference*.markdown)
+
+# Vendor ----------------------------------------------------------------------
+vendor/quickutils.lisp: vendor/make-quickutils.lisp
+	cd vendor && sbcl --noinform --load make-quickutils.lisp  --eval '(quit)'
+
+vendor: vendor/quickutils.lisp
+
+
+# Documentation ---------------------------------------------------------------
+$(apidocs): $(sourcefiles)
+	sbcl --noinform --load docs/api.lisp  --eval '(quit)'
+
+docs/build/index.html: $(docfiles) $(apidocs) docs/title
+	cd docs && ~/.virtualenvs/d/bin/d
+
+docs: docs/build/index.html
+
+pubdocs: docs
+	hg -R ~/src/sjl.bitbucket.org pull -u
+	rsync --delete -a ./docs/build/ ~/src/sjl.bitbucket.org/cl-pcg
+	hg -R ~/src/sjl.bitbucket.org commit -Am 'cl-pcg: Update site.'
+	hg -R ~/src/sjl.bitbucket.org push
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/README.markdown	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,9 @@
+`cl-pcg` is a [permuted congruential generator][pcg] implementation in pure
+Common Lisp.
+
+[pcg]: http://www.pcg-random.org/
+
+* **License:** MIT/X11
+* **Documentation:** <https://sjl.bitbucket.io/cl-pcg/>
+* **Mercurial:** <https://bitbucket.org/sjl/cl-pcg/>
+* **Git:** <https://github.com/sjl/cl-pcg/>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cl-pcg.asd	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,17 @@
+(asdf:defsystem :cl-pcg
+  :description "A bare-bones Permuted Congruential Generator implementation in pure Common Lisp."
+
+  :author "Steve Losh <steve@stevelosh.com>"
+
+  :license "MIT/X11"
+  :version "1.0.0"
+
+  :depends-on (#+sbcl :sb-rotate-byte)
+
+  :serial t
+  :components ((:module "vendor" :serial t
+                :components ((:file "quickutils-package")
+                             (:file "quickutils")))
+               (:file "package")
+               (:module "src" :serial t
+                :components ((:file "pcg")))))
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/01-installation.markdown	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,8 @@
+Installation
+============
+
+`cl-pcg` is compatible with Quicklisp, but not *in* Quicklisp (yet?).  You can
+clone the repository into your [Quicklisp local-projects][local] directory for
+now.
+
+[local]: https://www.quicklisp.org/beta/faq.html#local-project
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/02-usage.markdown	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,9 @@
+Usage
+=====
+
+cl-pcg is...
+
+[pcg]: http://www.pcg-random.org/
+
+[TOC]
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/03-reference.markdown	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,14 @@
+# API Reference
+
+The following is a list of all user-facing parts of `cl-pcg`.
+
+If there are backwards-incompatible changes to anything listed here, they will
+be noted in the changelog and the author will feel bad.
+
+Anything not listed here is subject to change at any time with no warning, so
+don't touch it.
+
+[TOC]
+
+## Package `PCG`
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/04-changelog.markdown	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,11 @@
+Changelog
+=========
+
+Here's the list of changes in each released version.
+
+[TOC]
+
+v1.0.0
+------
+
+Initial release.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/api.lisp	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,20 @@
+(ql:quickload "cl-d-api")
+
+(defparameter *header*
+  "The following is a list of all user-facing parts of `cl-pcg`.
+
+If there are backwards-incompatible changes to anything listed here, they will
+be noted in the changelog and the author will feel bad.
+
+Anything not listed here is subject to change at any time with no warning, so
+don't touch it.
+
+")
+
+(d-api:generate-documentation
+  :cl-pcg
+  #p"docs/03-reference.markdown"
+  (list "PCG")
+  *header*
+  :title "API Reference")
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/footer.markdown	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,14 @@
+<i>Made with Lisp and love by [Steve Losh][] in Reykjavík, Iceland.</i>
+
+[Steve Losh]: http://stevelosh.com/
+
+<script>
+  (function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
+  (i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
+  m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
+  })(window,document,'script','//www.google-analytics.com/analytics.js','ga');
+
+  ga('create', 'UA-15328874-3', 'auto');
+  ga('send', 'pageview');
+
+</script>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/index.markdown	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,9 @@
+`cl-pcg` is a [permuted congruential generator][pcg] implementation in pure
+Common Lisp.
+
+[pcg]: http://www.pcg-random.org/
+
+* **License:** MIT/X11
+* **Documentation:** <https://sjl.bitbucket.io/cl-pcg/>
+* **Mercurial:** <https://bitbucket.org/sjl/cl-pcg/>
+* **Git:** <https://github.com/sjl/cl-pcg/>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/docs/title	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,1 @@
+cl-pcg
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/package.lisp	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,4 @@
+(defpackage :pcg
+  (:use :cl :pcg.quickutils)
+  (:export
+    ))
Binary file src/pcg.fasl has changed
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/pcg.lisp	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,341 @@
+(in-package :pcg)
+
+;;;; Constants ----------------------------------------------------------------
+(defconstant +multiplier+ 6364136223846793005)
+(defconstant +modulus+ (expt 2 64))
+(defconstant +limit+ (1- (expt 2 64)))
+
+
+;;;; Utils ---------------------------------------------------------------------
+
+(defmacro defun-inline (name &body body)
+  "Like `defun`, but declaims `name` to be `inline`."
+  `(progn
+     (declaim (inline ,name))
+     (defun ,name ,@body)
+     ',name))
+
+(defmacro check-types (&rest place-type-pairs)
+  `(progn ,@(loop :for (place type . nil) :on place-type-pairs :by #'cddr
+                  :collect `(check-type ,place ,type))))
+
+
+;; Embedded from cl-utilities, with the bugfix of putting the inline declamation
+;; FIRST so it'll actually work.
+(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)))
+  (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)))))
+
+(defun-inline rotate-byte (count bytespec integer)
+  #+sbcl
+  (sb-rotate-byte:rotate-byte count bytespec integer)
+  #-sbcl
+  (rotate-byte% count bytespec integer))
+
+
+(defun-inline % (n)
+  (mod n +modulus+))
+
+(defun-inline *+ (a b increment)
+  (+ (* a b) increment))
+
+
+(deftype u64 ()
+  '(unsigned-byte 64))
+
+
+;;;; Data ---------------------------------------------------------------------
+(defstruct (pcg (:constructor actually-make-pcg))
+  (state 0 :type u64)
+  (increment 0 :type u64))
+
+(defun-inline compute-increment (stream-id)
+  (logior 1 (ash stream-id 1)))
+
+
+;;;; Permutations -------------------------------------------------------------
+(defun-inline permute-xor-shift (data)
+  (declare (optimize speed)
+           (type (unsigned-byte 37) data))
+  ;; The reference implemtation does this:
+  ;;
+  ;;   uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
+  ;;
+  ;; Which is a bunch of shit packed into one line because C programmers don't
+  ;; like typing or something.
+  ;;
+  ;; oldstate starts as a 64-bit value, which looks roughly like this:
+  ;;
+  ;;   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
+  ;;       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
+  ;;       the output value.
+  ;; * x - 27 not-very-random garbage low-order bits that we throw away.
+  ;;
+  ;; What this permutation does is xor the top half of the good bits into the
+  ;; bottom half of the good bits to mix them up a bit:
+  ;;
+  ;;   oldstate:       SSSSSbbb bbbbbbbb bbbbbbbb bbbbbbbb bbbbbxxx xxxxxxxx xxxxxxxx xxxxxxxx
+  ;;   oldstate >> 18: 00000000 00000000 00SSSSSb bbbbbbbb bbbbbbbb bbbbbbbb bbbbbbxx xxxxxxxx ...lost...
+  ;;   result:         SSSSSbbb bbbbbbbb bbBBBBBB BBBBBBBB BBBBBxxx xxxxxxxx xxxxxxxx xxxxxxxx
+  ;;
+  ;; Then it shifts by 27 to drop the garbage bits:
+  ;;
+  ;;   SSSSS bbbbbbbb bbbbbBBB BBBBBBBB BBBBBBBB
+  ;;
+  ;; And finally it assigns the resulting 37-bit value to a uint32_t which
+  ;; clears out the 5 remaining high-order selector bits.
+  (ldb (byte 32 0)
+       (logxor data (ash data -18))))
+
+(defun-inline permute-rotate (data selector)
+  (declare (optimize speed)
+           (type (unsigned-byte 32) data)
+           (type (unsigned-byte 5) selector))
+  (rotate-byte selector (byte 32 0) data))
+
+
+;;;; State Advancement --------------------------------------------------------
+(defun-inline advance-state (pcg)
+  (declare (optimize speed)
+           (type pcg pcg))
+  (setf (pcg-state pcg)
+        (% (*+ (pcg-state pcg) +multiplier+ (pcg-increment pcg))))
+  nil)
+
+
+;;;; Low-Level API ------------------------------------------------------------
+(declaim
+  (ftype (function (pcg) (unsigned-byte 32)) pcg-random%)
+
+  ;; 2^32 streams ought to be enough for anybody
+  (ftype (function (u64 (unsigned-byte 32))
+                   pcg)
+         make-pcg%)
+
+  (ftype (function (pcg (integer 1 (#.(expt 2 32))))
+                   (unsigned-byte 32))
+         pcg-random-bounded%)
+
+  (ftype (function (pcg (integer 1 32))
+                   (unsigned-byte 32))
+         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%))
+
+
+(defun-inline pcg-random% (pcg)
+  "Return a random `(unsigned-byte 32)`.
+
+  As a side effect, the state of `pcg` will be advanced.
+
+  This is a low-level function that assumes you are passing in the correct types.
+
+  "
+  (declare (optimize speed))
+  (let* ((state (pcg-state pcg))
+         (data (ldb (byte 37 (- 64 37)) state))
+         (selector (ldb (byte 5 (- 64 5)) state)))
+    (advance-state pcg)
+    (permute-rotate (permute-xor-shift data)
+                    selector)))
+
+
+(defun-inline make-pcg% (seed stream-id)
+  "Create and return a new `pcg` for the given `seed` and `stream-id`.
+
+  This is a low-level function that assumes you are passing in the correct types.
+
+  "
+  (declare (optimize speed))
+  (let* ((increment (compute-increment stream-id))
+         (pcg (actually-make-pcg :state 0 :increment increment)))
+    (pcg-random% pcg)
+    (setf (pcg-state pcg)
+          (% (+ (pcg-state pcg) seed)))
+    (pcg-random% pcg)
+    pcg))
+
+
+(defun-inline 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.
+
+  This is a low-level function that assumes you are passing in the correct types.
+
+  "
+  (declare (optimize speed))
+  (loop :with threshold = (mod (expt 2 32) bound)
+        :for n = (pcg-random% pcg)
+        :when (>= n threshold)
+        :do (return (values (mod n bound)))))
+
+(defun-inline pcg-random-bits% (pcg count)
+  "Return a random `(unsigned-byte COUNT)`.
+
+  As a side effect, the state of `pcg` will be advanced.
+
+  `count` must be between `1` and `32` (though `32` would be identical to just
+  calling `pcg-random%`).
+
+  This is a low-level function that assumes you are passing in the correct types.
+
+  "
+  (declare (optimize speed))
+  (ldb (byte count 0) (pcg-random% pcg)))
+
+(defun-inline 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.
+
+  This is a low-level function that assumes you are passing in the correct types.
+
+  "
+  (declare (optimize speed))
+  (/ (pcg-random% pcg)
+     (coerce (expt 2 32) 'single-float)))
+
+
+(defun-inline pcg-advance% (pcg steps)
+  (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
+  ;;
+  ;; The original paper uses single-letter names like G, C, i, f, h because
+  ;; math, but we're all literate here so let's use words instead.
+  (loop
+    :with increment = (pcg-increment pcg)
+    :with multiplier-accumulator :of-type u64 = 1 ; G
+    :with increment-accumulator :of-type u64 = 0 ; C
+    :for i :of-type u64 = steps :then (ash i -1) ; i
+    :for inc :of-type u64 = increment :then (% (* inc (1+ mul))) ; f
+    :for mul :of-type u64 = +multiplier+ :then (% (expt mul 2)) ; h
+    :while (plusp i)
+    :when (oddp i)
+    :do (setf multiplier-accumulator (% (* multiplier-accumulator mul))
+              increment-accumulator (% (*+ increment-accumulator mul inc)))
+    :finally
+    (setf (pcg-state pcg)
+          (% (*+ (pcg-state pcg) multiplier-accumulator increment-accumulator)))))
+
+(defun-inline pcg-rewind% (pcg steps)
+  (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))
+
+
+(defparameter *global-pcg* (make-pcg))
+
+(deftype pcg-designator ()
+  '(or (eql t) pcg))
+
+(defun-inline resolve-pcg (pcg-designator)
+  (if (eq t pcg-designator)
+    *global-pcg*
+    pcg-designator))
+
+
+(defun pcg-random (pcg)
+  (check-types pcg pcg-designator)
+  (pcg-random% (resolve-pcg pcg)))
+
+(defun pcg-random-float (pcg)
+  (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)
+  (check-types pcg pcg-designator
+               bound (and (unsigned-byte 32)
+                          (integer 1)))
+  (pcg-random-bounded% (resolve-pcg pcg) bound))
+
+(defun pcg-random-range (pcg min max)
+  (check-types pcg pcg-designator
+               min (unsigned-byte 32)
+               max (unsigned-byte 32))
+  (assert (< min max) (min max))
+  (+ min (pcg-random-bounded% (resolve-pcg pcg) (- max min))))
+
+(defun pcg-random-range-inclusive (pcg min max)
+  (check-types pcg pcg-designator
+               min (unsigned-byte 32)
+               max (unsigned-byte 32))
+  (assert (<= min max) (min max))
+  (pcg-random-range (resolve-pcg pcg) min (1+ max)))
+
+(defun pcg-advance (pcg steps)
+  (check-types pcg pcg-designator
+              steps u64)
+  (pcg-advance% (resolve-pcg pcg) steps))
+
+(defun pcg-rewind (pcg steps)
+  (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*)))
+
+; (losh:gnuplot
+;   (data 10000)
+;   :x #'identity
+;   :y (lambda (i) (/ 1.0 10000))
+;   :line-width 1
+;   :smooth :cumulative)
+
+; (defun slow-advance (pcg n)
+;   (dotimes (_ n) (pcg-random pcg))
+;   nil)
+
+; (defun test (n)
+;   (let ((a (make-pcg))
+;         (b (make-pcg)))
+;     (pcg-advance% b n)
+;     (assert (= (pcg-random a) (pcg-random b)))))
+
+;;;; TODO ----------------------------------------------------------------------
+;;; https://experilous.com/1/blog/post/perfect-fast-random-floating-point-numbers
+;;; as an alternative float generation scheme.
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/vendor/make-quickutils.lisp	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,9 @@
+(ql:quickload 'quickutil)
+
+(qtlc:save-utils-as
+  "quickutils.lisp"
+  :utilities '(
+
+
+               )
+  :package "PCG.QUICKUTILS")
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/vendor/quickutils-package.lisp	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,12 @@
+(eval-when (:compile-toplevel :load-toplevel :execute)
+  (unless (find-package "PCG.QUICKUTILS")
+    (defpackage "PCG.QUICKUTILS"
+      (:documentation "Package that contains Quickutil utility functions.")
+      (:use :cl))))
+
+(in-package "PCG.QUICKUTILS")
+
+;; need to define this here so sbcl will shut the hell up about it being
+;; undefined when compiling quickutils.lisp.  computers are trash.
+(defparameter *utilities* nil)
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/vendor/quickutils.lisp	Fri Feb 10 23:31:25 2017 +0000
@@ -0,0 +1,17 @@
+;;;; This file was automatically generated by Quickutil.
+;;;; See http://quickutil.org for details.
+
+;;;; To regenerate:
+;;;; (qtlc:save-utils-as "quickutils.lisp" :ensure-package T :package "PCG.QUICKUTILS")
+
+(eval-when (:compile-toplevel :load-toplevel :execute)
+  (unless (find-package "PCG.QUICKUTILS")
+    (defpackage "PCG.QUICKUTILS"
+      (:documentation "Package that contains Quickutil utility functions.")
+      (:use #:cl))))
+
+(in-package "PCG.QUICKUTILS")
+
+NIL
+
+;;;; END OF quickutils.lisp ;;;;