Note: The MT 19937 rng is included in Maxima in versions 5.9.2 and later. The Maxima functions are make_random_state, set_random_state, and random. See the Maxima manual for more information.

Modified rand-mt19937.lisp source code here.

Diff wrt to original rand-mt19937.lisp here.

The modified rand-mt19937.lisp passes two tests. (1) The first 1,000,000 chunks produced by the modified code are the same as the sequence of chunks produced by the unmodified code. This is a CMUCL-only test. See Compare chunks below. (2) The first 10,000 integers and floats of various sizes have the same bits as one would expect from the corresponding sequence of chunks. Clisp and CMUCL pass this test, GCL does too if one ignores %random-single-float (since that yields a 64 bit number in GCL). In addition the first 10,000 chunks are the same as the chunks produced by CMUCL. Clisp and GCL pass this test also. See Compare non-chunks below.

Here's a short test script (to run in Maxima) to compare the CMUCL MT to JMT. This scripts uses jmt.lisp from http://www.lisp-p.org/jmt/.

load ("jmt.lisp");
load ("rand-mt19937.lisp");

:lisp (defmfun $mt19937_random (x) (mt19937::random x))

:lisp (defmfun $jmt_random (x) (mt-random x))

showtime: true$

l: makelist (random (123456), i, 1, 10000)$

l: makelist (mt19937_random (123456), i, 1, 10000)$

l: makelist (jmt_random (123456), i, 1, 10000)$

Compare chunks Here is a comparison of the sequence of chunks produced by CMUCL MT19937 to that produced by the modified code. This test assumes the package MAXIMA is defined. "cmucl-rand-mt19937.lisp" is the original rand-mt19937.lisp with minimal modifications to make it load from the Lisp prompt.

(compile-file "rand-mt19937.lisp")
(compile-file "cmucl-rand-mt19937.lisp")

(load "rand-mt19937.x86f")         ; substantially modified for Maxima
(load "cmucl-rand-mt19937.x86f")   ; minimal modifications, just enough to get it to load from cmd prompt

(defun set-mt19937 (seed)
  (let ((my-state (mt19937::make-random-object :state (mt19937::init-random-state seed))))
    (setq mt19937::*random-state* (mt19937::make-random-state my-state)))
  t)

(defun set-cmucl-mt19937 (seed)
  (let ((my-state (cmucl-mt19937::make-random-object :state (cmucl-mt19937::init-random-state seed))))
    (setq cmucl-mt19937::*random-state* (cmucl-mt19937::make-random-state my-state)))
  t)

(set-mt19937 9876)
(set-cmucl-mt19937 9876)

; Sequence of chunks should be the same.

(time
  (dotimes (i 1000000 t)
    (let
      ((x (mt19937::random-chunk mt19937::*random-state*))
       (y (cmucl-mt19937::random-chunk cmucl-mt19937::*random-state*)))
      (cond
        ((not (= x y)) (format t "i: ~S, chunks differ: ~S, ~S~%" i y x) (return nil))))))

Compare non-chunks Here is a comparison of the sequences of integers and floats. This test assumes the package MAXIMA is defined.

(load "rand-mt19937.lisp")

(defun set-mt19937 (seed)
  (let ((my-state (mt19937::make-random-object :state (mt19937::init-random-state seed))))
    (setq mt19937::*random-state* (mt19937::make-random-state my-state)))
  t)

; Bits for a sequence of m-bit integers should be same as low m bits of n-bit integers for m < n.

(setq nn 10000)

(format t "first generate ~S samples of varying numbers of bits ...~%" nn)

(set-mt19937 54321)
(setq a23 (make-array nn :element-type 'integer))
(dotimes (i nn)
  (setf (aref a23 i) (mt19937::%random-integer #X800000 mt19937::*random-state*)))

(set-mt19937 54321)
(setq asingle (make-array nn :element-type 'integer))
(dotimes (i nn)
  (let ((x (mt19937::%random-single-float 1f0 mt19937::*random-state*)))
    (setf (aref asingle i) (logand #X7FFFFF (integer-decode-float (1+ x))))))

(set-mt19937 54321)
(setq achunk (make-array nn :element-type 'integer))
(dotimes (i nn)
  (setf (aref achunk i) (mt19937::random-chunk mt19937::*random-state*)))

(set-mt19937 54321)
(setq a32 (make-array nn :element-type 'integer))
(dotimes (i nn)
  (setf (aref a32 i) (mt19937::%random-integer #X100000000 mt19937::*random-state*)))

(set-mt19937 54321)
(setq a52 (make-array nn :element-type 'integer))
(dotimes (i nn)
  (setf (aref a52 i) (mt19937::%random-integer #X10000000000000 mt19937::*random-state*)))

(set-mt19937 54321)
(setq adouble (make-array nn :element-type 'integer))
(dotimes (i nn)
  (let ((x (mt19937::%random-double-float 1d0 mt19937::*random-state*)))
    (setf (aref adouble i) (logand #XFFFFFFFFFFFFF (integer-decode-float (1+ x))))))

(set-mt19937 54321)
(setq a64 (make-array nn :element-type 'integer))
(dotimes (i nn)
  (setf (aref a64 i) (mt19937::%random-integer #X10000000000000000 mt19937::*random-state*)))

#+cmu   (setq filename "tmp-cmucl.out")
#+clisp (setq filename "tmp-clisp.out")
#+gcl   (setq filename "tmp-gcl.out")

(format t "... now dump arrays to ~S ...~%" filename)

(with-open-file (s filename :direction :output :if-exists :supersede)
  (dotimes (i nn)
    (cond
      ((evenp i)
        (format s "~8X  ~8X  ~8X  ~8X  ~16X  ~16X  ~16X~%" (aref a23 i) (aref asingle i) (aref achunk i) (aref a32 i) (aref a52 (ash i -1)) (aref adouble (ash i -1)) (aref a64 (ash i -1))))
      (t
        (format s "~8X  ~8X  ~8X  ~8X~%" (aref a23 i) (aref asingle i) (aref achunk i) (aref a32 i))))))

(format t "... now verify that shorter values are same as corresponding bits in longer ones.~%")

(dotimes (i nn)
  (cond
    ((evenp i)
     (cond
       ((not (= (aref a23 i) (rem (aref achunk i) (expt 2 23))))
        (format t "oops, a23[~S]~%" i))
       ((not (= (aref asingle i) (rem (aref achunk i) (expt 2 23))))
        (format t "oops, asingle[~S]~%" i))
       ((not (= (aref a32 i) (aref achunk i)))
        (format t "oops, a32[~S]~%" i))
       ((not (= (aref a52 (ash i -1)) (rem (aref a64 (ash i -1)) (expt 2 52))))
        (format t "oops, a52[~S]~%" (ash i -1)))
       ((not (= (aref adouble (ash i -1)) (rem (aref a64 (ash i -1)) (expt 2 52))))
        (format t "oops, adouble[~S]~%" (ash i -1)))))
    (t
      (let ((aa (aref a64 (ash i -1))) (bb (+ (ash (aref a32 i) 32) (aref a32 (1- i)))))
        (cond
          ((not (= (aref a23 i) (rem (aref achunk i) (expt 2 23))))
           (format t "oops, a23[~S]~%" i))
          ((not (= (aref asingle i) (rem (aref achunk i) (expt 2 23))))
           (format t "oops, asingle[~S]~%" i))
          ((not (= (aref a32 i) (aref achunk i)))
           (format t "oops, a32[~S]~%" i))
          ((not (= aa bb))
           (format t "oops, a64[~S]~%" (ash i -1))))))))