# HG changeset patch # User Steve Losh # Date 1486769485 0 # Node ID c6a3ab353556d8a36d17019b89d0dbedec17b0b3 Initial commit diff -r 000000000000 -r c6a3ab353556 .ffignore --- /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 diff -r 000000000000 -r c6a3ab353556 .hgignore --- /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 diff -r 000000000000 -r c6a3ab353556 .lispwords --- /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*) diff -r 000000000000 -r c6a3ab353556 Makefile --- /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 diff -r 000000000000 -r c6a3ab353556 README.markdown --- /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:** +* **Mercurial:** +* **Git:** diff -r 000000000000 -r c6a3ab353556 cl-pcg.asd --- /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 " + + :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"))))) diff -r 000000000000 -r c6a3ab353556 docs/01-installation.markdown --- /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 diff -r 000000000000 -r c6a3ab353556 docs/02-usage.markdown --- /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] + diff -r 000000000000 -r c6a3ab353556 docs/03-reference.markdown --- /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` + diff -r 000000000000 -r c6a3ab353556 docs/04-changelog.markdown --- /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. diff -r 000000000000 -r c6a3ab353556 docs/api.lisp --- /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") + diff -r 000000000000 -r c6a3ab353556 docs/footer.markdown --- /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 @@ +Made with Lisp and love by [Steve Losh][] in Reykjavík, Iceland. + +[Steve Losh]: http://stevelosh.com/ + + diff -r 000000000000 -r c6a3ab353556 docs/index.markdown --- /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:** +* **Mercurial:** +* **Git:** diff -r 000000000000 -r c6a3ab353556 docs/title --- /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 diff -r 000000000000 -r c6a3ab353556 package.lisp --- /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 + )) diff -r 000000000000 -r c6a3ab353556 src/pcg.fasl Binary file src/pcg.fasl has changed diff -r 000000000000 -r c6a3ab353556 src/pcg.lisp --- /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. + + diff -r 000000000000 -r c6a3ab353556 vendor/make-quickutils.lisp --- /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") diff -r 000000000000 -r c6a3ab353556 vendor/quickutils-package.lisp --- /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) + diff -r 000000000000 -r c6a3ab353556 vendor/quickutils.lisp --- /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 ;;;;