Monte Carlo and Las Vegas mapping
author |
Steve Losh <steve@stevelosh.com> |
date |
Sat, 09 Jul 2016 16:29:35 +0000 |
parents |
277982735a9d
|
children |
6f72eefef02e
|
branches/tags |
(none) |
files |
src/random-numbers.lisp |
Changes
--- 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))))