[Maxima-commits] CVS: maxima/src rand-mt19937.lisp,NONE,1.1

[Robert D. <rob...@us...>]

Update of /cvsroot/maxima/maxima/src
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv11530/src

Added Files:
	rand-mt19937.lisp 
Log Message:
MT19937 random number generator, adapted from CMUCL sources,
and test script.

--- NEW FILE: rand-mt19937.lisp ---
;;; 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