Mersenne twister 随机数算法实现 in Scheme

这个实现基本上是从 Wiki 上的 Python 版翻译过来的,大量使用了赋值。


;; 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))))

我使用了自己定义的 for 宏,以及 range 函数来实现 Python 风格的 for 循环,下面是相关的定义:

(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)))))

 使用了 R6RS 特有的一些函数及语法,使用时不要忘记在头部加上 (import (rnrs),如果还依赖别的库请查阅 R6RS 文档。

posted @ 2017-01-06 09:28  fmcdr  阅读(1286)  评论(0编辑  收藏  举报