Hammersley QRNG

README.mess

@@ -1,5 +1,5 @@
 ## About random-state
-This is a collection of pseudo random number generators (PRNGs). While Common Lisp does provide a ``RANDOM`` function, it does not allow the user to pass an explicit seed, nor to portably exchange the random state between implementations. This can be a headache in cases like games, where a controlled seeding process is very useful.
+This is a collection of pseudo random number generators (PRNGs) and quasi-random number generators (QRNGs). While Common Lisp does provide a ``RANDOM`` function, it does not allow the user to pass an explicit seed, nor to portably exchange the random state between implementations. This can be a headache in cases like games, where a controlled seeding process is very useful.
 
 ## How To
 For both curiosity and convenience, this library offers multiple algorithms to generate random numbers, as well as a bunch of generally useful methods to produce desired ranges.

documentation.lisp

@@ -275,6 +275,17 @@ See HISTOGRAM"))
 
 ;; * RNGs
 (docs:define-docs
+  (type hammersley
+    "The Hammersley quasi-random number sequence.
+
+This is a multivariate generator with default dimensionality of 3.
+To change the dimensionality, pass a :LEAP initarg that is an array of
+the desired dimensionality. Each element in the LEAP array should be
+an integer greater or equal to 1, and can be used to advance the
+sequence more quickly for each dimension.
+
+See https://en.wikipedia.org/wiki/Low-discrepancy_sequence")
+
   (type linear-congruence
     "A very simple random number generator based on linear congruence.
 

hammersley.lisp

@@ -0,0 +1,23 @@
+(in-package #:org.shirakumo.random-state)
+
+(define-generator hammersley (make-list (length (hammersley-leap generator)) :initial-element 'single-float) (hash-generator)
+  ((leap (make-array 3 :element-type '(unsigned-byte 32) :initial-element 1) :type (simple-array (unsigned-byte 32) (*))))
+  (:copy
+   (make-hammersley :%seed (hammersley-%seed generator)
+                    :index (hammersley-index generator)
+                    :leap (make-array (length (hammersley-leap generator))
+                                      :element-type '(unsigned-byte 32)
+                                      :initial-contents (hammersley-leap generator))))
+  (:hash
+   (flet ((dim (base leap seed)
+            (let ((seed2 (+ seed (* index leap)))
+                  (base-inv (float (/ base) 0f0))
+                  (r 0.0))
+              (loop while (/= 0 seed2)
+                    do (incf r (* (mod seed2 base) base-inv))
+                       (setf base-inv (/ base-inv base))
+                       (setf seed2 (truncate seed2 base))
+                    finally (return r)))))
+     (let ((result (make-array (length leap) :element-type 'single-float)))
+       (dotimes (i (length result) result)
+         (setf (aref result i) (dim (prime (1+ i)) (aref leap i) (ldb (byte 8 i) seed))))))))

primes.lisp

@@ -0,0 +1,34 @@
+(in-package #:org.shirakumo.random-state)
+
+(let ((pis-1 (make-array 0 :adjustable T :fill-pointer T :element-type '(signed-byte 64)))
+      (pis (make-array 0 :adjustable T :fill-pointer T :element-type '(unsigned-byte 64)))
+      (primes (make-array 1 :adjustable T :fill-pointer T :element-type '(unsigned-byte 64) :initial-contents '(2))))
+  (defun compute-ruiz-pis-part1 (j)
+    (if (< j (length pis-1))
+        (aref pis-1 j)
+        (let* ((s-max (floor (sqrt j)))
+               (sigma1 (loop for s from 1 to s-max
+                             sum (- (floor (/ (1- j) s)) (floor (/ j s)))))
+               (result (floor (* (/ 2 j) (1+ sigma1)))))
+          (vector-push-extend result pis-1)
+          result)))
+
+  (defun compute-ruiz-pi (k)
+    (if (< k (length pis))
+        (aref pis k)
+        (let ((result (cond
+                        ((= k 1) 0)
+                        ((= k 2) (1+ (compute-ruiz-pis-part1 2)))
+                        (T (+ 1 (compute-ruiz-pis-part1 k) (compute-ruiz-pi (1- k)))))))
+          (vector-push-extend result pis)
+          result)))
+
+  (defun prime (n)
+    (declare (type (unsigned-byte 32) n))
+    (if (< n (length primes))
+        (aref primes n)
+        (let* ((n (1+ n))
+               (result (1+ (loop for k from 1 to (* 2 (1+ (floor (* n (log n)))))
+                                 summing (- 1 (floor (/ (compute-ruiz-pi k) n)))))))
+          (vector-push-extend result primes)
+          result))))

random-state.asd

@@ -11,6 +11,7 @@
                (:file "toolkit")
                (:file "generator")
                (:file "protocol")
+               (:file "primes")
                (:file "linear-congruence")
                (:file "mersenne-twister")
                (:file "middle-square")
@@ -23,6 +24,7 @@
                (:file "murmurhash")
                (:file "cityhash")
                (:file "xorshift")
+               (:file "hammersley")
                (:file "implementation")
                (:file "documentation"))
   :depends-on (:documentation-utils)