; mersenne twister
(define (logand a b)
(if (or (zero? a) (zero? b)) 0
(+ (* (logand (floor (/ a 2)) (floor (/ b 2))) 2)
(if (or (even? a) (even? b)) 0 1))))
(define (logior x y)
(cond ((= x y) x)
((zero? x) y)
((zero? y) x)
(else
(+ (* (logior (quotient x 2) (quotient y 2)) 2)
(if (and (even? x) (even? y)) 0 1)))))
(define (logxor a b)
(cond ((zero? a) b)
((zero? b) a)
(else
(+ (* (logxor (floor (/ a 2)) (floor (/ b 2))) 2)
(if (even? a)
(if (even? b) 0 1)
(if (even? b) 1 0))))))
(define (ash int cnt)
(if (negative? cnt)
(let ((n (expt 2 (- cnt))))
(if (negative? int)
(+ -1 (quotient (+ 1 int) n))
(quotient int n)))
(* (expt 2 cnt) int)))
(define n 624)
(define m 397)
(define upper-mask #x80000000)
(define lower-mask #x7FFFFFFF)
(define mt (make-vector n 0))
(define mti (+ n 1))
(define mag01 (vector #x0 #x9908B0DF))
(define (sgenrand seed)
(vector-set! mt 0 (logand seed #xFFFFFFFF))
(do ((mti 1 (+ mti 1))) ((= mti n))
(vector-set! mt mti
(logand (* 69069 (vector-ref mt (- mti 1))) #xFFFFFFFF)))
(set! mti n))
(define (genrand)
(when (<= n mti)
(when (= mti (+ n 1)) (sgenrand 4357))
(do ((kk 0 (+ kk 1))) ((= kk (- n m)))
(let ((y (logior
(logand (vector-ref mt kk) upper-mask)
(logand (vector-ref mt (+ kk 1)) lower-mask))))
(vector-set! mt kk
(logxor
(logxor (vector-ref mt (+ kk m)) (ash y -1))
(vector-ref mag01 (logand y #x1))))))
(do ((kk (- n m) (+ kk 1))) ((= kk (- n 1)))
(let ((y (logior
(logand (vector-ref mt kk) upper-mask)
(logand (vector-ref mt (+ kk 1)) lower-mask))))
(vector-set! mt kk
(logxor
(logxor (vector-ref mt (+ kk m (- n))) (ash y -1))
(vector-ref mag01 (logand y #x1))))))
(let ((y (logior
(logand (vector-ref mt (- n 1)) upper-mask)
(logand (vector-ref mt 0) lower-mask))))
(vector-set! mt (- n 1)
(logxor
(logxor (vector-ref mt (- m 1)) (ash y -1))
(vector-ref mag01 (logand y #x1)))))
(set! mti 0))
(let* ((y (vector-ref mt mti))
(y (logxor y (ash y -11)))
(y (logxor y (logand (ash y 7) #x9D2C5680)))
(y (logxor y (logand (ash y 15) #xEFC60000)))
(y (logxor y (ash y -18))))
(set! mti (+ mti 1))
y))
(define (test-rand)
(sgenrand 4357)
(do ((j 0 (+ j 1))) ((= j 1000))
(display (genrand))
(if (= (modulo j 8) 7) (newline) (display " "))))
(test-rand)