Modified from rev 1.11 of CMUCL rand-mt19937.lisp source code. This loads and executes successfully in Clisp, GCL, and CMUCL. Updated 2005/01/12: construct multi-chunk integers by simple concatenation (no overlap, no thrownaway bits), & rework defpackage stuff. Updated 2005/01/10: reimplemented floating point stuff.
;;; 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)))))
;;; begin Maxima-specific stuff
(in-package "MAXIMA")
(defmfun $set_random_state (x)
"Copy the argument, and assign the copy to MT19937::*RANDOM-STATE*.
Returns '$done."
(setq mt19937::*random-state* (mt19937::make-random-state x))
'$done)
(defmfun $make_random_state (x)
"Returns a new random state object. If argument is an integer or array,
use argument to initialize random state. Otherwise punt to MT19937::MAKE-RANDOM-STATE."
(cond
((or (integerp x) (arrayp x))
(mt19937::make-random-object :state (mt19937::init-random-state x)))
(t
(mt19937::make-random-state x))))
(defmfun $random (x)
"Returns the next number from this generator.
Punt to MT19937::RANDOM."
(mt19937::random x))
;;; end Maxima-specific stuff
Last edited on Friday, January 14, 2005 9:05:24 am.
Page content is licensed under the terms of the GNU General Public License unless otherwise noted. Copyright for each edit is held by the person who made it.
  
 |
|