# HG changeset patch
# User Steve Losh <steve@stevelosh.com>
# Date 1468081775 0
# Node ID 3c908eea3940ee3d313e0425a975406682b31e24
# Parent  277982735a9d957e9d31f60818b92146e0d1bd91
Monte Carlo and Las Vegas mapping

diff -r 277982735a9d -r 3c908eea3940 src/random-numbers.lisp
--- a/src/random-numbers.lisp	Sat Jul 09 16:00:42 2016 +0000
+++ b/src/random-numbers.lisp	Sat Jul 09 16:29:35 2016 +0000
@@ -1,6 +1,7 @@
 (in-package #:sand.random-numbers)
 
 
+;;;; Types, etc
 (declaim (optimize (speed 1) (safety 1) (debug 3)))
 
 (deftype positive-fixnum () `(integer 1 ,most-positive-fixnum))
@@ -8,6 +9,7 @@
 (deftype nonnegative-fixnum () `(integer 1 ,most-positive-fixnum))
 (deftype nonpositive-fixnum () `(integer ,most-negative-fixnum -1))
 
+;;;; Utils
 (defun* +mod ((x nonnegative-fixnum)
               (y nonnegative-fixnum)
               (m positive-fixnum))
@@ -16,17 +18,19 @@
     (- x (- m y))))
 
 
+;;;; Random Number Generators
 (defun* make-linear-congruential-rng
     ((modulus positive-fixnum)
      (multiplier nonnegative-fixnum)
      (increment nonnegative-fixnum)
      (seed nonnegative-fixnum))
   (let ((val seed))
-    (lambda (incr)
-      (loop :repeat incr :do
-            (setf val (mod (+ (* multiplier val)
-                              increment)
-                           modulus))))))
+    (lambda (msg)
+      (ecase msg
+        (:next (setf val (mod (+ (* multiplier val)
+                                 increment)
+                              modulus)))
+        (:modulus modulus)))))
 
 (defun* make-linear-congruential-rng-fast%
     ((modulus positive-fixnum)
@@ -35,12 +39,24 @@
      (seed nonnegative-fixnum))
   (declare (optimize (speed 3) (safety 0) (debug 0)))
   (let ((val seed))
-    (lambda (incr)
-      (declare (positive-fixnum incr))
-      (loop :repeat incr :do
-            (setf val (mod (+ (the nonnegative-fixnum (* multiplier val))
-                              increment)
-                           modulus))))))
+    (lambda (msg)
+      (ecase msg
+        (:next (setf val (mod (+ (the nonnegative-fixnum (* multiplier val))
+                                 increment)
+                              modulus)))
+        (:modulus modulus)))))
+
+
+(declaim (inline rng-next rng-modulus))
+
+(defun* rng-next ((generator function))
+  (:returns positive-fixnum)
+  (funcall generator :next))
+
+(defun* rng-modulus ((generator function))
+  (:returns positive-fixnum)
+  (funcall generator :modulus))
+
 
 (define-compiler-macro make-linear-congruential-rng
     (&whole form
@@ -53,3 +69,50 @@
     form))
 
 
+(defparameter *generator* (make-linear-congruential-rng 601 15 4 354))
+
+
+(defun rand ()
+  (rng-next *generator*))
+
+(defun rand-float ()
+  (float (/ (rng-next *generator*)
+            (rng-modulus *generator*))))
+
+
+;;;; Mapping
+;;; The Monte Carlo method is bad because it's biased, but it's fast.
+;;;
+;;; Basically we take our generator that generates say 1-8, and map the range
+;;; ABC onto it:
+;;;
+;;;    1 2 3 4 5 6 7 8
+;;;    A B C A B C A B
+;;;
+;;; Notice that it's not uniform.
+(defun* monte-carlo ((width positive-fixnum))
+  (:returns positive-fixnum)
+  (mod (rng-next *generator*) width))
+
+
+;;; The Las Vegas method is a bit slower, but unbiased.  We group the random
+;;; numbers into contiguous buckets, with the last "partial bucket" being
+;;; excess.  If we hit that one we just loop and try again:
+;;;
+;;;    1 2 3 4 5 6 7 8
+;;;    A A B B C C retry
+(defun* las-vegas ((width positive-fixnum))
+  (:returns positive-fixnum)
+  (let* ((modulus (rng-modulus *generator*))
+         (bucket-width (truncate (/ modulus width))))
+    (iterate
+      (for bucket = (truncate (/ (rng-next *generator*)
+                                 bucket-width)))
+      (finding bucket :such-that (< bucket width)))))
+
+
+(defun rand-range-bad (min max)
+  (+ min (monte-carlo (- max min))))
+
+(defun rand-range (min max)
+  (+ min (las-vegas (- max min))))