Clone mt19937 state from 624 consecutive 32-bit output words

mt19937-clone.lisp

(defpackage #:snippets/mt19937-clone
  (:documentation
   "Clone an mt19937 random number generator (the generator used by
SBCL's random-state machinery) from 624 consecutive words.")
  (:use #:cl)
  ;; We build an SBCL-specific random state.
  (:import-from #:sb-kernel)
  (:export
   #:clone-random-state))

(in-package #:snippets/mt19937-clone)

(shadowing-import '(sb-kernel::random-state-state
                    sb-kernel::init-random-state))

;; ;; Debugging utils

(defmacro dsetq (place value-form)
  `(progn
     (format t "Setting ~S to ~S (form ~S)~%"
             ',place ,value-form ',value-form)
     (setq ,place ,value-form)))

(defun terprint (x)
  (print x)
  (terpri))

(defun hex (x)
  (format nil "~X" x))

(defun next (random-state)
  (random #x100000000 random-state))

(defun next-624 (random-state)
  (loop repeat 624
        collect (next random-state)))

(defun new-rand (x)
  (let ((s (make-random-state)))
    (init-random-state x (random-state-state s))
    s))

;; The bee's knees, ported from my Go cryptopals impl.

(defconstant m-7-00 #b00000000000000000000000001111111)
(defconstant m-7-07 #b00000000000000000011111110000000)
(defconstant m-7-14 #b00000000000111111100000000000000)
(defconstant m-7-21 #b00001111111000000000000000000000)
(defconstant m-4-28 #b11110000000000000000000000000000)
(defconstant m-10-11 #b00000000000111111111100000000000)
(defconstant m-10-00 #b00000000000000000000001111111111)
(defconstant m-15-00 #b00000000000000000111111111111111)
(defconstant m-15-15 #b00111111111111111000000000000000)
(defconstant m-2-30 #b11000000000000000000000000000000)

(defun ~ (x)
  (- #xFFFFFFFF x))

(defun untemper (y &aux (m 0))
  (declare (type (unsigned-byte 32) y m))
  ;; step 4 - y ^= y >> 18
  (setq y (logxor y (ash y -18)))
  ;; step 3 - y ^= (y << 15) & 0xEFC60000
  ;; ABCDEFGH IJKLMNOP QRSTUVWX YZabcdef
  ;; PQRSTUVW XYZabcde f0000000 00000000
  ;; yyxxxxxx xxxxxxxx xRSTUVWX YZabcdef
  (setq m y)
  (setq m (logand m m-15-00))
  (setq m (logior m (logxor (logand y m-15-15)
                            (logand (ash m 15) #xEFC60000))))
  (setq m (logior m (logxor (logand y m-2-30)
                            (logand (ash m 15)
                                    #xEFC60000
                                    (~ m-15-15)))))
  (setq y m)
  ;; step 2 - y ^= (y << 7) & 0x9D2C5680
  (setq m y)
  (setq m (logand m m-7-00))
  (setq m (logior m (logxor (logand y m-7-07)
                            (logand (ash m 7)
                                    #x9D2C5680))))
  (setq m (logior m (logxor (logand y m-7-14)
                            (logand (ash m 7)
                                    #x9D2C5680
                                    (~ m-7-07)))))
  (setq m (logior m
                  (logxor (logand y m-7-21)
                          (logand (ash m 7)
                                  #x9D2C5680
                                  (~ m-7-07)
                                  (~ m-7-14)))))
  (setq m (logior m
                  (logxor (logand y m-4-28)
                          (logand (ash m 7)
                                  #x9D2C5680
                                  (~ m-7-07)
                                  (~ m-7-14)
                                  (~ m-7-21)))))
  (setq y m)
  ;; step 1 - y ^= y >> 11
  ;; ABCDEFGH IJKLMNOP QRSTUVWX YZabcdef
  ;; 00000000 000ABCDE FGHIJKLM NOPQRSTU
  ;; ABCDEFGH IJKxxxxx xxxxxxyy yyyyyyyy
  (setq m (logxor y (ash y -11)))
  (setq m (logand m (~ m-10-00)))
  (setq m (logior m (logxor (logand y m-10-00)
                            (ash (logand m m-10-11) -11))))
  (setq y m)
  y)

(defun clone-random-state (vals)
  (assert (= (length vals) 624))
  (let ((r2 (new-rand 123)))
    (loop for v in vals
          for i upfrom 0
          do (setf (aref (random-state-state r2)
                         (+ i 3))
                   (untemper v)))
    r2))

;; Tests

(defun test-clone-random-state ()
  (let* ((vals (loop repeat 624
                     collect (next *random-state*)))
         (r2 (clone-random-state vals)))
    (dotimes (i 100)
      (assert (= (next *random-state*)
                 (next r2))))
    r2))

(defun test-untemper ()
  (let ((s (make-random-state)))
    (init-random-state 123 (random-state-state s))
    (let* ((y (next s))
           (actual (aref (random-state-state s) 3)))
      (assert (= actual (untemper y))))))