public
Description: A collection of functions arranged by programming language
Homepage: http://www.hackinghat.com/
Clone URL: git://github.com/hackinghat/common-library.git
common-library / clojure / mersenne_twister.clj
100644 105 lines (95 sloc) 4.201 kb
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
;;;; A Clojure implementation of the Mersenne Twister algorithm
;;;;
;;;; Copyright (c) 2009 Steve Knight <stknig@gmail.com>
;;;;
;;;; Copyright (C) 1997, 1999 Makoto Matsumoto and Takuji Nishimura.
;;;; matumoto@math.keio.ac.jp
;;;;
;;;; Permission is hereby granted, free of charge, to any person obtaining
;;;; a copy of this software and associated documentation files (the
;;;; "Software"), to deal in the Software without restriction, including
;;;; without limitation the rights to use, copy, modify, merge, publish,
;;;; distribute, sublicense, and/or sell copies of the Software, and to
;;;; permit persons to whom the Software is furnished to do so, subject to
;;;; the following conditions:
;;;;
;;;; The above copyright notice and this permission notice shall be
;;;; included in all copies or substantial portions of the Software.
;;;;
;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
;;;; MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
;;;; LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
;;;; OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
;;;; WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
;;;;
                                    
(ns com.hackinghat.common-library.mersenne-twister)
 
(def n 624)
(def m 397)
(def matrix-a 0x9908b0df)
(def upper-mask 0x80000000)
(def lower-mask 0x7fffffff)
(def tempering-mask-b 0x9d2c5680)
(def tempering-mask-c 0xefc60000)
(defn tempering-shift-u [y] (bit-shift-right y 11))
(defn tempering-shift-s [y] (bit-shift-left y 7))
(defn tempering-shift-t [y] (bit-shift-left y 15))
(defn tempering-shift-l [y] (bit-shift-right y 18))
 
(def mt (ref []))
(def mti (ref (inc n)))
 
(defn reset []
  (dosync
    (ref-set mt [])
    (ref-set mti (inc n)))
  nil)
    
(defn sgenrand [seed]
  "Seed the original array"
  (dosync
    (ref-set mt [])
    (loop [i (int 0) s seed]
      (if (< i n)
        (let [mti (bit-and s 0xffff0000)
              s-1 (inc (* 69069 s))]
           (ref-set mt (into @mt [(bit-or mti (bit-shift-right (bit-and s-1 0xffff0000) 16))]))
           (recur (inc i) (inc (* 69069 s-1))))))
    (ref-set mti n)))
 
(defn -genrand-inner [v1 v2 v3]
  "Factored from the reference implementation"
  (let [y (bit-or (bit-and v1 upper-mask)
                  (bit-and v2 lower-mask))]
    (bit-xor (bit-xor v3 (bit-shift-right y 1))
             (if (= 0 (bit-and y 0x1))
                0
                matrix-a))))
                
(defn -genrand-nth []
  "Returns the next number from the seuquence (factored from the reference implementation)"
  (dosync
    (if
      (>= @mti n)
        (do
          (if (= @mti (inc n))
            (sgenrand 4357))
          ;; Build a new vector of integers from the old one (and itself!)
          (ref-set mt (loop [k (int 0) v []]
                          (if (< k n)
                             (let [v-1 (cond
                                         (< k (- n m)) (conj v (-genrand-inner (nth @mt k) (nth @mt (inc k)) (nth @mt (+ k m))))
                                         (< k (dec n)) (conj v (-genrand-inner (nth @mt k) (nth @mt (inc k)) (nth v (+ k (- m n)))))
                                         (= k (dec n)) (conj v (-genrand-inner (nth @mt (dec n)) (nth v 0) (nth v (dec m)))))]
                             (recur (inc k) v-1))
                            v)))
           (ref-set mti 0)))
      ;; Now return the element and advance the counter
      (let [retval (int (nth @mt @mti))]
        (ref-set mti (inc @mti))
        retval)))
   
(defn genrand []
  "Generate the next number in the sequence, if this number is the nth in the sequence then
a new state table will be generated in this call"
  (let [y (-genrand-nth)
        y1 (bit-xor y (tempering-shift-u y))
        y2 (bit-xor y1 (bit-and (tempering-shift-s y1) tempering-mask-b))
        y3 (bit-xor y2 (bit-and (tempering-shift-t y2) tempering-mask-c))
        y4 (bit-xor y3 (tempering-shift-l y3))]
     y4))