;; Mersenne twister algorithm from Wikipedia

;; returns a closure that returns a pseudo-random integer

;; for each call

;;

(define (make-MT19937 seed)

  ;; some bitwise procedure alias for short

  (define << bitwise-arithmetic-shift-left)

  (define >> bitwise-arithmetic-shift-right)

  (define xor bitwise-xor)

  (letrec ((mt (make-vector 624))

           (index 624)

           ;; reset index

           (twist

            (lambda ()

              (for i in (range 624)

                   (let ((y (bitwise-and

                             #xffffffff

                             (+ (bitwise-and (vector-ref mt i)

                                             #x80000000)

                                (bitwise-and (vector-ref mt (mod (+ i 1) 624))

                                             #x7fffffff)))))

                     (vector-set! mt i (xor (vector-ref mt (mod (+ i 397) 624))

                                            (>> y 1)))

                     (when (odd? y)

                           (vector-set!

                            mt i (xor (vector-ref mt i) #x9908b0df)))))

              (set! index 0)))

           ;; generates a number

           (extract_number

            (lambda ()

              (when (>= index 624)

                    (twist))

              (let ((y (vector-ref mt index)))

                (set! y (xor y (>> y 11)))

                (set! y (xor y (bitwise-and (<< y 7) 2636928640)))

                (set! y (xor y (bitwise-and (<< y 15) 4022730752)))

                (set! y (xor y (>> y 18)))

                (set! index (+ index 1))

                (bitwise-and #xffffffff y)))))

    (vector-set! mt 0 seed)

    ;; initialize the vector

    (for i in (range 1 624)

         (vector-set!

          mt

          i

          (bitwise-and (+ i

                          (* 1812433253

                             (bitwise-xor (vector-ref mt (- i 1))

                                          (>> (vector-ref mt (- i 1)) 30))))

                       #xffffffff)))

    ;; return a closure

    (lambda ()

      (extract_number))))

;; It may be better to set the seed as the system clock

;; but that involves different implementations

(define generator (make-MT19937 4294967296))  ;; the seed

(define (randint . arg)

  (if (null? arg)

      (generator)

      (mod (generator) (car arg))))

(define-syntax for

  (syntax-rules ()

    ;; loop in list

    ;; (for i in '(a b c) do something...)

    ((_ i in lst body ...)

     (let loop ((l lst))

       (unless (null? l)

               (let ((i (car l)))

                 body ...

                 (loop (cdr l))))))))

(define range

  (let ((make-range

         (lambda (first end step)

           (if (or (= step 0)

                   (> (abs (- (+ first step) end))

                      (abs (- first end))))

               (error 'range "wrong `step' leads to an infinite loop")

               (let iter ((cnt first) (result '()))

                 (cond ((or (and (> step 0) (>= cnt end))

                            (and (< step 0) (<= cnt end)))

                        (reverse result))

                       (else (iter (+ cnt step) (cons cnt result)))))))))

    (case-lambda

     ((a) (make-range 0 a 1))

     ((a b) (make-range a b 1))

     ((a b c) (make-range a b c)))))