1 ! Copyright (C) 2005, 2008 Doug Coleman.
2 ! See http://factorcode.org/license.txt for BSD license.
3 ! mersenne twister based on
4 ! http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c
5 USING: arrays kernel math namespaces sequences system init
6 accessors math.ranges random circular math.bitwise
8 IN: random.mersenne-twister
12 TUPLE: mersenne-twister seq i ;
16 : mt-a HEX: 9908b0df ; inline
18 : calculate-y ( n seq -- y )
20 [ [ 1+ ] [ nth ] bi* 31 bits ] 2bi bitor ; inline
22 : (mt-generate) ( n seq -- next-mt )
25 [ 2/ ] [ odd? mt-a 0 ? ] bi bitxor
27 [ mt-m + ] [ nth ] bi*
30 : mt-generate ( mt -- )
33 [ (mt-generate) ] [ set-nth ] 2bi
37 : init-mt-formula ( i seq -- f(seq[i]) )
38 dupd nth dup -30 shift bitxor 1812433253 * + 1+ 32 bits ;
40 : init-mt-rest ( seq -- )
42 [ init-mt-formula ] [ >r 1+ r> set-nth ] 2bi
45 : init-mt-seq ( seed -- seq )
46 32 bits mt-n 0 <array> <circular>
47 [ set-first ] [ init-mt-rest ] [ ] tri ;
49 : mt-temper ( y -- yt )
51 dup 7 shift HEX: 9d2c5680 bitand bitxor
52 dup 15 shift HEX: efc60000 bitand bitxor
53 dup -18 shift bitxor ; inline
55 : next-index ( mt -- i )
56 dup i>> dup mt-n < [ nip ] [ drop mt-generate 0 ] if ;
60 : <mersenne-twister> ( seed -- obj )
61 init-mt-seq 0 mersenne-twister boa
64 M: mersenne-twister seed-random ( mt seed -- )
65 init-mt-seq >>seq drop ;
67 M: mersenne-twister random-32* ( mt -- r )
69 [ seq>> nth mt-temper ]
70 [ [ 1+ ] change-i drop ] tri ;
75 [ 32 random-bits ] with-system-random
76 <mersenne-twister> random-generator set-global
77 ] "bootstrap.random" add-init-hook