Diff of Modified rand-mt19937.lisp source code wrt to original. x86-dependent stuff and /dev/urandom have fallen by the wayside. Declarations modified to make Clisp and GCL happy.
--- rand-mt19937.lisp-orig 2005-01-07 21:24:05.000000000 -0700 +++ rand-mt19937.lisp 2005-01-14 09:55:22.000000000 -0700 @@ -1,10 +1,26 @@ -;;; -*- Mode: Lisp; Package: Kernel -*- -;;; -;;; ********************************************************************** -;;; This code was written by Douglas T. Crosher and Raymond Toy based -;;; on public domain code from Carnegie Mellon University and has been -;;; placed in the Public domain, and is provided 'as is'. -;;; -(ext:file-comment - "$Header: /project/cmucl/cvsroot/src/code/rand-mt19937.lisp,v 1.11 2003/03/06 09:31:06 emarsden Exp $") -;;; +;;; 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. @@ -22,6 +38 @@ -(in-package "LISP") -(export '(random-state random-state-p random *random-state* - make-random-state)) - -(in-package "KERNEL") -(export '(%random-single-float %random-double-float random-chunk init-random-state)) +(in-package "MT19937") @@ -29 +40,7 @@ - +(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) +;;; @@ -42,21 +59,4 @@ -;; If /dev/urandom is available, it is used to generate random data as -;; the seed. Otherwise, the current time is used as the seed. -;; -;; The /dev/urandom device exists on Linux, FreeBSD, Solaris 8 and -;; later. (You need to have patch 112438-01 for Solaris 8 to get that -;; device, though). It returns pseudorandom data with entropy that the -;; kernel collects from the environment. Unlike the related -;; /dev/random, this device does not block when the entropy pool has -;; been depleted. -(defun generate-seed (&optional (nwords 1)) - (cond ((probe-file "/dev/urandom") - (let ((words (make-array nwords :element-type '(unsigned-byte 32)))) - (with-open-file (rand "/dev/urandom" - :direction :input - :element-type '(unsigned-byte 32)) - (read-sequence words rand)) - (if (= nwords 1) - (aref words 0) - words))) - (t - (logand (get-universal-time) #xffffffff)))) +;; The current time is used as the seed. + +(defun generate-seed () + (logand (get-universal-time) #xffffffff)) @@ -186,2 +186 @@ - (:constructor make-random-object) - (:make-load-form-fun :just-dump-it-normally)) + (:constructor make-random-object)) @@ -195,2 +194,3 @@ - copy of it. If STATE is T then return a random state generated from - the universal time or /dev/urandom if available." + 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>))''." @@ -210,7 +209,0 @@ -(defun rand-mt19937-initializer () - (init-random-state (generate-seed) - (random-state-state *random-state*))) - -(pushnew 'rand-mt19937-initializer ext:*after-save-initializations*) - - @@ -231,8 +223,0 @@ -(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) -;;; -#-x86 @@ -266 +250,0 @@ -#-x86 @@ -285,8 +268,0 @@ -;;; Using inline VOP support, only available on the x86 so far. -#+x86 -(defun random-chunk (state) - (declare (type random-state state)) - (vm::random-mt19937 (random-state-state state))) - - - @@ -295,4 +270,0 @@ -;;; Handle the single or double float case of RANDOM. We generate a float -;;; between 0.0 and 1.0 by clobbering the significand of 1.0 with random bits, -;;; then subtracting 1.0. This hides the fact that we have a hidden bit. -;;; @@ -304,10 +275,0 @@ -(defun %random-single-float (arg state) - (declare (type (single-float (0f0)) arg) - (type random-state state)) - (* arg - (- (make-single-float - (dpb (ash (random-chunk state) - (- vm:single-float-digits random-chunk-length)) - vm:single-float-significand-byte - (single-float-bits 1.0))) - 1.0))) @@ -319,7 +280,0 @@ -;;; 32bit version -;;; -#+nil -(defun %random-double-float (arg state) - (declare (type (double-float (0d0)) arg) - (type random-state state)) - (* (float (random-chunk state) 1d0) (/ 1d0 (expt 2 32)))) @@ -327,15 +282,8 @@ -;;; 53bit version. -;;; -#-x86 -(defun %random-double-float (arg state) - (declare (type (double-float (0d0)) arg) - (type random-state state)) - (* arg - (- (lisp::make-double-float - (dpb (ash (random-chunk state) - (- vm:double-float-digits random-chunk-length - vm:word-bits)) - vm:double-float-significand-byte - (lisp::double-float-high-bits 1d0)) - (random-chunk state)) - 1d0))) +(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))) @@ -343,2 +290,0 @@ -;;; Using a faster inline VOP. -#+x86 @@ -346,41 +292,7 @@ - (declare (type (double-float (0d0)) arg) - (type random-state state)) - (let ((state-vector (random-state-state state))) - (* arg - (- (lisp::make-double-float - (dpb (ash (vm::random-mt19937 state-vector) - (- vm:double-float-digits random-chunk-length - vm:word-bits)) - vm:double-float-significand-byte - (lisp::double-float-high-bits 1d0)) - (vm::random-mt19937 state-vector)) - 1d0)))) - -#+long-float -(declaim (inline %random-long-float)) -#+long-float -(declaim (ftype (function ((long-float (0l0)) random-state) (long-float 0l0)) - %random-long-float)) - -;;; Using a faster inline VOP. -#+(and long-float x86) -(defun %random-long-float (arg state) - (declare (type (long-float (0l0)) arg) - (type random-state state)) - (let ((state-vector (random-state-state state))) - (* arg - (- (lisp::make-long-float - (lisp::long-float-exp-bits 1l0) - (logior (vm::random-mt19937 state-vector) vm:long-float-hidden-bit) - (vm::random-mt19937 state-vector)) - 1l0)))) - -#+(and long-float sparc) -(defun %random-long-float (arg state) - (declare (type (long-float (0l0)) arg) - (type random-state state)) - (* arg - (- (lisp::make-long-float - (lisp::long-float-exp-bits 1l0) ; X needs more work - (random-chunk state) (random-chunk state) (random-chunk state)) - 1l0))) + "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))) @@ -388 +299,0 @@ - @@ -391,16 +301,0 @@ -;;; Amount we overlap chunks by when building a large integer to make up for -;;; the loss of randomness in the low bits. -;;; -(defconstant random-integer-overlap 3) - -;;; Extra bits of randomness that we generate before taking the value MOD the -;;; limit, to avoid loss of randomness near the limit. -;;; -(defconstant random-integer-extra-bits 10) - -;;; Largest fixnum we can compute from one chunk of bits. -;;; -(defconstant random-fixnum-max - (1- (ash 1 (- random-chunk-length random-integer-extra-bits)))) - - @@ -409,0 +305,7 @@ + "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." @@ -411,9 +313,7 @@ - (let ((shift (- random-chunk-length random-integer-overlap))) - (do ((bits (random-chunk state) - (logxor (ash bits shift) (random-chunk state))) - (count (+ (integer-length arg) - (- random-integer-extra-bits shift)) - (- count shift))) - ((minusp count) - (rem bits arg)) - (declare (fixnum count))))) + (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)))) @@ -422,4 +322,3 @@ - "Generate a uniformly distributed pseudo-random number between zero - and Arg. State, if supplied, is the random state to use." - (declare (inline %random-single-float %random-double-float - #+long-float %long-float)) + "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)) @@ -427,2 +325,0 @@ - ((and (fixnump arg) (<= arg random-fixnum-max) (> arg 0)) - (rem (random-chunk state) arg)) @@ -433,3 +329,0 @@ - #+long-float - ((and (typep arg 'long-float) (> arg 0.0L0)) - (%random-long-float arg state)) @@ -442,0 +337,27 @@ + +;;; 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:08:02 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.
  
 |
|