Commit 831284f0 authored by Pierre Neidhardt's avatar Pierre Neidhardt
Browse files

Version 1.1.

parents
Introduction
============
MT19937 is a portable Mersenne Twister random number generator. It
is mainly a modification of CMUCL's random number generator with all
the CMUCL-specific parts taken out.
It is faster than the JMT Mersenne Twister implementation, but
significantly slower than native random number generators provided by
major Common Lisp implementations. For light use this shouldn't be a
problem, since it is still very fast.
It should be very stable, since it's based on stable CMUCL code. It
has been tested on CMUCL, SBCL, LispWorks, Allegro CL, GCL, CLISP, and
Corman Lisp.
What's the point?
=================
Why, you might ask, would you want to use a portable, slower random
number generator? The answer is consistancy. The portable version of
this code was originally created for Maxima, a computer algebra
system. They wanted the results of the random number generator to be
portable across several implementations and platforms, so that if you
used a certain seed on CMUCL the numbers generated would be the same
as you would get with the same seed on GCL or CLISP. This was more
important than achieving the maximum possible speed. You may have
similar problems.
Usage
=====
MT19937 is a plug-in replacement for the Common Lisp random-number
generator. The MT19937 package exports all the Common Lisp symbols
related to random number generation, so you just need to load MT19937
and call its functions instead of the built-in ones. For example:
;; Load MT19937
(asdf:oos 'asdf:load-op :mt19937)
;; Make random numbers with your implementation's random
;; number generator.
(random 1234567)
(random 42.56)
(random 3.1415d0)
;; Make random numbers using MT19937
(mt19937:random 1234567)
(mt19937:random 42.56)
(mt19937:random 3.1415d0)
;; MT19937 has its own random state
(eq *random-state*
mt19937:*random-state*) => nil
License
=======
MT19937 is in the public domain.
\ No newline at end of file
;; MT19937 benchmarking versus the RNG provided by the implementation.
;; As you can see, the standard RANDOM is much faster.
(defun test-mt19937-float ()
(loop repeat 250000
do (mt19937:random 1234.5678))
(loop repeat 250000
do (mt19937:random 1234.5678d0)))
(defun test-mt19937-integer ()
(loop repeat 1000000
do (mt19937:random 12345678)))
(defun test-cl-float ()
(loop repeat 250000
do (random 1234.5678))
(loop repeat 250000
do (random 1234.5678d0)))
(defun test-cl-integer ()
(loop repeat 1000000
do (random 12345678)))
(defun benchmark (&optional (stream *standard-output*))
"Benchmark MT19937:random versus cl:random and print the
results to STREAM (defaults to *standard-output*)"
(let ((*trace-output* stream))
(format stream "~&~%MT1997 integer test~%")
(time (test-mt19937-integer))
(format stream "~&~%Common Lisp integer test~%")
(time (test-cl-integer))
(format stream "~&~%~%MT1997 floating-point test~%")
(time (test-mt19937-float))
(format stream "~&~%Common Lisp floating-point test~%")
(time (test-cl-float))))
;; Here's an example of usage of MT19937. It's pretty simple.
;; Load MT19937
(asdf:oos 'asdf:load-op :mt19937)
;; Make random numbers with your implementation's random
;; number generator.
(random 1234567)
(random 42.56)
(random 3.1415d0)
;; Make random numbers using MT19937
(mt19937:random 1234567)
(mt19937:random 42.56)
(mt19937:random 3.1415d0)
;; MT19937 has its own random state
(eq *random-state*
mt19937:*random-state*)
\ No newline at end of file
;; -*- Lisp -*-
(defpackage #:mt19937-system
(:use #:common-lisp #:asdf))
(in-package #:mt19937-system)
(defsystem mt19937
:description "Portable MT19937 Mersenne Twister random number generator"
:author "Douglas T. Crosher and Raymond Toy"
:licence "Public domain"
:version "1.1"
:components ((:file "mt19937")))
\ No newline at end of file
;;; Mersenne Twister MT19937, adapted from CMUCL rand-mt19937.lisp -r1.11 (2003/03/06)
;;; CMUCL version by Douglas T. Crosher and Raymond Toy based
;;; on public domain code from Carnegie Mellon University.
;;; Modified for Maxima by Robert Dodier.
;;; (1) Construct floating point numbers using portable operations.
;;; (2) Construct large integers using all bits of each chunk.
(defpackage "MT19937"
(:use :common-lisp)
(:shadow #:random-state
#:random-state-p
#:random
#:*random-state*
#:make-random-state)
(:export #:random-state
#:random-state-p
#:random
#:*random-state*
#:make-random-state
#:%random-single-float
#:%random-double-float
#:random-chunk
#:init-random-state))
;;; Begin MT19937 implementation.
;;; **********************************************************************
;;;
;;; Support for the Mersenne Twister, MT19937, random number generator
;;; due to Matsumoto and Nishimura. This implementation has been
;;; placed in the public domain with permission from M. Matsumoto.
;;;
;;; Makoto Matsumoto and T. Nishimura, "Mersenne twister: A
;;; 623-dimensionally equidistributed uniform pseudorandom number
;;; generator.", ACM Transactions on Modeling and Computer Simulation,
;;; 1997, to appear.
(in-package "MT19937")
(defconstant mt19937-n 624)
(defconstant mt19937-m 397)
(defconstant mt19937-upper-mask #x80000000)
(defconstant mt19937-lower-mask #x7fffffff)
(defconstant mt19937-b #x9D2C5680)
(defconstant mt19937-c #xEFC60000)
;;;
;;;; Random state hackery:
;;; The state is stored in a (simple-array (unsigned-byte 32) (627))
;;; wrapped in a random-state structure:
;;;
;;; 0-1: Constant matrix A. [0, #x9908b0df]
;;; 2: Index k.
;;; 3-626: State.
;; GENERATE-SEED
;;
;; Generate a random seed that can be used for seeding the generator.
;; The current time is used as the seed.
(defun generate-seed ()
(logand (get-universal-time) #xffffffff))
;; New initializer proposed by Takuji Nishimura and Makota Matsumoto.
;; (See http://www.math.keio.ac.jp/~matumoto/MT2002/emt19937ar.html)
;;
;; This corrects a deficiency in the original initializer wherein the
;; MSB of the seed was not well represented in the state.
;;
;; The initialization routine is described below. Let s be the seed,
;; mt[] be the state vector. Then the algorithm is
;;
;; mt[0] = s & 0xffffffffUL
;;
;; for (k = 1; k < N; k++) {
;; mt[k] = 1812433253 * (mt[k-1] ^ (mt[k-1] >> 30)) + k
;; mt[k] &= 0xffffffffUL
;; }
;;
;; The multiplier is from Knuth TAOCP Vol2, 3rd Ed., p. 106.
;;
(defun int-init-random-state (&optional (seed 5489) state)
(declare (type (integer 1 #xffffffff) seed))
(let ((state (or state (make-array 627 :element-type '(unsigned-byte 32)))))
(declare (type (simple-array (unsigned-byte 32) (627)) state))
(setf (aref state 0) 0)
(setf (aref state 1) #x9908b0df)
(setf (aref state 2) mt19937-n)
(setf (aref state 3) seed)
(do ((k 1 (1+ k)))
((>= k 624))
(declare (type (mod 625) k))
(let ((prev (aref state (+ 3 (1- k)))))
(setf (aref state (+ 3 k))
(logand (+ (* 1812433253 (logxor prev (ash prev -30)))
k)
#xffffffff))))
state))
;; Initialize from an array.
;;
;; Here is the algorithm, in C. init_genrand is the initalizer above,
;; init_key is the seed vector of length key_length.
;;
;; init_genrand(19650218UL);
;; i=1; j=0;
;; k = (N>key_length ? N : key_length);
;; for (; k; k--) {
;; mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
;; + init_key[j] + j; /* non linear */
;; mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
;; i++; j++;
;; if (i>=N) {
;; mt[0] = mt[N-1]; i=1;
;; }
;; if (j>=key_length) {
;; j=0;
;; }
;; }
;; for (k=N-1; k; k--) {
;; mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
;; - i; /* non linear */
;; mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
;; i++;
;; if (i>=N) { mt[0] = mt[N-1]; i=1; }
;; }
;;
;; mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
;;
(defun vec-init-random-state (key &optional state)
(declare (type (array (unsigned-byte 32) (*)) key))
(let ((key-len (length key))
(state (init-random-state 19650218 state))
(i 1)
(j 0))
(loop for k from (max key-len mt19937-n) above 0 do
(let ((prev (aref state (+ 3 (1- i)))))
(setf (aref state (+ 3 i))
(ldb (byte 32 0)
(+ (aref key j) j
(logxor (aref state (+ 3 i))
(ldb (byte 32 0)
(* 1664525
(logxor prev (ash prev -30))))))))
(incf i)
(incf j)
(when (>= i mt19937-n)
(setf (aref state 3)
(aref state (+ 3 (- mt19937-n 1))))
(setf i 1))
(when (>= j key-len)
(setf j 0))))
(loop for k from (1- mt19937-n) above 0 do
(let ((prev (aref state (+ 3 (1- i)))))
(setf (aref state (+ 3 i))
(ldb (byte 32 0)
(- (logxor (aref state (+ 3 i))
(* 1566083941
(logxor prev (ash prev -30))))
i)))
(incf i)
(when (>= i mt19937-n)
(setf (aref state 3)
(aref state (+ 3 (- mt19937-n 1))))
(setf i 1))))
(setf (aref state 3) #x80000000)
state))
;;
(defun init-random-state (&optional (seed 5489) state)
"Generate an random state vector from the given SEED. The seed can be
either an integer or a vector of (unsigned-byte 32)"
(declare (type (or null integer
(array (unsigned-byte 32) (*)))
seed))
(etypecase seed
(integer
(int-init-random-state (ldb (byte 32 0) seed) state))
((array (unsigned-byte 32) (*))
(vec-init-random-state seed state))))
(defstruct (random-state
(:constructor make-random-object))
(state (init-random-state) :type (simple-array (unsigned-byte 32) (627))))
(defvar *random-state* (make-random-object))
(defun make-random-state (&optional state)
"Make a random state object. If STATE is not supplied, return a copy
of the default random state. If STATE is a random state, then return a
copy of STATE. If STATE is T then return a random state generated from
the universal time. To make a random state from an integer seed, try
``(make-random-object :state (init-random-state <seed>))''."
(flet ((copy-random-state (state)
(let ((state (random-state-state state))
(new-state
(make-array 627 :element-type '(unsigned-byte 32))))
(dotimes (i 627)
(setf (aref new-state i) (aref state i)))
(make-random-object :state new-state))))
(cond ((not state) (copy-random-state *random-state*))
((random-state-p state) (copy-random-state state))
((eq state t)
(make-random-object :state (init-random-state (generate-seed))))
(t (error "Argument is not a RANDOM-STATE, T or NIL: ~S" state)))))
;;;; Random entries:
;;; Size of the chunks returned by random-chunk.
;;;
(defconstant random-chunk-length 32)
;;; random-chunk -- Internal
;;;
;;; This function generaters a 32bit integer between 0 and #xffffffff
;;; inclusive.
;;;
(declaim (inline random-chunk))
;;;
;;; Portable implementation.
(defun random-mt19937-update (state)
(declare (type (simple-array (unsigned-byte 32) (627)) state)
(optimize (speed 3) (safety 0)))
(let ((y 0))
(declare (type (unsigned-byte 32) y))
(do ((kk 3 (1+ kk)))
((>= kk (+ 3 (- mt19937-n mt19937-m))))
(declare (type (mod 628) kk))
(setf y (logior (logand (aref state kk) mt19937-upper-mask)
(logand (aref state (1+ kk)) mt19937-lower-mask)))
(setf (aref state kk) (logxor (aref state (+ kk mt19937-m))
(ash y -1) (aref state (logand y 1)))))
(do ((kk (+ (- mt19937-n mt19937-m) 3) (1+ kk)))
((>= kk (+ (1- mt19937-n) 3)))
(declare (type (mod 628) kk))
(setf y (logior (logand (aref state kk) mt19937-upper-mask)
(logand (aref state (1+ kk)) mt19937-lower-mask)))
(setf (aref state kk) (logxor (aref state (+ kk (- mt19937-m mt19937-n)))
(ash y -1) (aref state (logand y 1)))))
(setf y (logior (logand (aref state (+ 3 (1- mt19937-n)))
mt19937-upper-mask)
(logand (aref state 3) mt19937-lower-mask)))
(setf (aref state (+ 3 (1- mt19937-n)))
(logxor (aref state (+ 3 (1- mt19937-m)))
(ash y -1) (aref state (logand y 1)))))
(values))
;;;
(defun random-chunk (state)
(declare (type random-state state)
(optimize (speed 3) (safety 0)))
(let* ((state (random-state-state state))
(k (aref state 2)))
(declare (type (mod 628) k))
(when (= k mt19937-n)
(random-mt19937-update state)
(setf k 0))
(setf (aref state 2) (1+ k))
(let ((y (aref state (+ 3 k))))
(declare (type (unsigned-byte 32) y))
(setf y (logxor y (ash y -11)))
(setf y (logxor y (ash (logand y (ash mt19937-b -7)) 7)))
(setf y (logxor y (ash (logand y (ash mt19937-c -15)) 15)))
(setf y (logxor y (ash y -18)))
y)))
;;; %RANDOM-SINGLE-FLOAT, %RANDOM-DOUBLE-FLOAT -- Interface
;;;
(declaim (inline %random-single-float %random-double-float))
(declaim (ftype (function ((single-float (0f0)) random-state)
(single-float 0f0))
%random-single-float))
;;;
;;;
(declaim (ftype (function ((double-float (0d0)) random-state)
(double-float 0d0))
%random-double-float))
;;;
;;;
(defun %random-single-float (arg state)
"Handle the single or double float case of RANDOM. We generate a float
in [0f0, 1f0) by clobbering the mantissa of 1f0 with random bits (23 bits);
this yields a number in [1f0, 2f0). Then 1f0 is subtracted."
(let*
((random-mantissa-bits (%random-integer (expt 2 23) state))
(random-unit-float (- (scale-float (float (+ (expt 2 23) random-mantissa-bits) 1f0) -23) 1f0)))
(* arg random-unit-float)))
(defun %random-double-float (arg state)
"Handle the single or double float case of RANDOM. We generate a float
in [0d0, 1d0) by clobbering the mantissa of 1d0 with random bits (52 bits);
this yields a number in [1d0, 2d0). Then 1d0 is subtracted."
(let*
((random-mantissa-bits (%random-integer (expt 2 52) state))
(random-unit-double (- (scale-float (float (+ (expt 2 52) random-mantissa-bits) 1d0) -52) 1d0)))
(* arg random-unit-double)))
;;;; Random integers:
;;; %RANDOM-INTEGER -- Internal
;;;
(defun %random-integer (arg state)
"Generates an integer greater than or equal to zero and less than Arg.
Successive chunks are concatenated without overlap to construct integers
larger than a single chunk. The return value has this property:
If two integers are generated from the same state with Arg equal to 2^m and 2^n,
respectively, then bit k is the same in both integers for 0 <= k < min(m,n).
Each call to %RANDOM-INTEGER consumes at least one chunk; bits left over
from previous chunks are not re-used."
(declare (type (integer 1) arg) (type random-state state))
(do*
((nchunks (ceiling (integer-length (1- arg)) random-chunk-length) (1- nchunks))
(new-bits 0 (random-chunk state))
(bits 0 (logior bits (ash new-bits shift)))
(shift 0 (+ shift random-chunk-length)))
((= 0 nchunks)
(rem bits arg))))
(defun random (arg &optional (state *random-state*))
"Generates a uniformly distributed pseudo-random number greater than or equal to zero
and less than Arg. State, if supplied, is the random state to use."
(declare (inline %random-single-float %random-double-float))
(cond
((and (typep arg 'single-float) (> arg 0.0F0))
(%random-single-float arg state))
((and (typep arg 'double-float) (> arg 0.0D0))
(%random-double-float arg state))
((and (integerp arg) (> arg 0))
(%random-integer arg state))
(t
(error 'simple-type-error
:expected-type '(or (integer 1) (float (0))) :datum arg
:format-control "Argument is not a positive integer or a positive float: ~S"
:format-arguments (list arg)))))
;; This should speed things up some on calls to random with constant
;; upper bound arguments, like (random 1234), by doing some error
;; checking at compile time and transforming it into an inline call
;; to one of the type-specific random number generation routines.
;; More speed should be possible, but this is better than nothing.
(define-compiler-macro random (&whole form arg &optional (state '*random-state*)
&environment env)
(cond
((typep arg 'single-float env)
(if (and (constantp arg env) (> arg 0.0F0))
`(%random-single-float ,arg ,state)))
((typep arg 'double-float env)
(if (and (constantp arg env) (> arg 0.0D0))
`(%random-double-float ,arg ,state)))
((typep arg 'integer env)
(if (and (constantp arg env) (> arg 0))
`(%random-integer ,arg ,state)))
(t form)))
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment