randomアルゴリズムの改良

x022235 commented 12 years ago

動機

下記テスト結果を見ると 0.2.2.235 の random 関数の質はあまりよくない

make-random-state が時間 (秒単位) を seed にしているため、連続して make-random-state を呼び出した場合に同じ乱数列が生成される

(let ((s1 (make-random-state t))
    (s2 (make-random-state t)))
(dotimes (_ 10)
  (format t "~S~%" (list (random 100 s1) (random 100 s2)))))
(72 72)
(44 44)
(92 92)
(54 54)
(89 89)
(72 72)
(22 22)
(93 93)
(19 19)
(72 72)

ゴール

質のよい乱数を生成できるようにする
性能が劣化しないようにする
実現案
random 関数は [0, limit) の乱数を生成するという仕様であるため、乱数を double 値で生成してから limit で乗算するという実装になる
- そのため、IEEE 754 形式の倍精度浮動小数点数を直接生成できる dSFMT を利用する
- MT や SFMT だと整数から浮動小数点数への変換コストが発生する
seed の初期値を CryptGenRandom function で生成するようにする
性能評価

性能が劣化していないことを確認した。

(defun bench (n)
  (let ((s (get-internal-real-time)))
    (dotimes (_ n)
      (random 1))
    (- (get-internal-real-time)s)))
(compile 'bench)
(gc)
(bench 40000000)
;; 修正前
4773
;; 修正後
4711

互換性留意事項

内部ベクタのサイズが変わったので、S notation の表記も変わる
そのため、従来の xyzzy で生成した make-random-state の S notation を新しい xyzzy で eval できない
(defvar *my-random-state* #.(make-random-state)) というコードをコンパイルするとバイトコードに S notation で書かれるので、それは新しい xyzzy でロードできない
read-time-eval しないとバイトコードに埋め込まれることはないし、そのようなコードの存在は確認できなかったためこの非互換性は仕様とする
テスト結果

なお、以下のテストは make-random-state にシードとして整数を指定できるように改造した上でテストした。

最下位ビットの問題

#|
/* TEST1 */
main(void)
 {
    int i,j,k;
    for (i=0;i<4;i++) for (j=0;j<4;j++) {
       printf("%4d ",i+j*4 );
       srandom(i+j*4);
       for (k=0;k<62;k++) printf("%ld",random()&1);
       printf("\n");
   }
   return 0;
}
|#
(defun test1 ()
  (dotimes (i 4)
    (dotimes (j 4)
      (let* ((s (+ (* j 4) i))
             (r (make-random-state s)))
        (format t "~4D " s)
        (dotimes (k 62)
          (princ (logand 1 (random #x7FFFFFFF r))))
        (terpri)))))

実行結果

   0 01010101001101101110011000111000111110000010011101001000010100
   4 01010101001101101110011000111000111110000010011101001000010100
   8 01010101001101101110011000111000111110000010011101001000010100
  12 01010101001101101110011000111000111110000010011101001000010100
   1 00000011001110100110111011000111011001111001000110100011011000
   5 00000011001110100110111011000111011001111001000110100011011000
   9 00000011001110100110111011000111011001111001000110100011011000
  13 00000011001110100110111011000111011001111001000110100011011000
   2 01010101001101101110011000111000111110000010011101001000010100
   6 01010101001101101110011000111000111110000010011101001000010100
  10 01010101001101101110011000111000111110000010011101001000010100
  14 01010101001101101110011000111000111110000010011101001000010100
   3 00000011001110100110111011000111011001111001000110100011011000
   7 00000011001110100110111011000111011001111001000110100011011000
  11 00000011001110100110111011000111011001111001000110100011011000
  15 00000011001110100110111011000111011001111001000110100011011000

上位ビットの問題

#|
/* TEST2 */
using System;
class Program {
    static int SkippedRand(int seed,int NumOfSkip) {
        Random r=new Random(seed);
        for (int i=0;i<NumOfSkip;i++) r.Next();
        return r.Next();
    }
    static void Main(string[] args) {
        for (int i=1;i<=100;i++) {
            Console.Write("{0,2} 番目:",i);
            for (int s=0;s<64;s++)
                Console.Write("{0}",SkippedRand(s,i-1)>>30);
            Console.WriteLine();
        }
    }
}
|#
(defun skipped-rand (s n)
  (let ((r (make-random-state s)))
    (dotimes (_ n)
      (random #x7FFFFFFF r))
    (random #x7FFFFFFF r)))

(defun test2 ()
  (dotimes (i 100)
    (format t "~2D 番目: " i)
    (dotimes (s 64)
      (princ (ash (skipped-rand s i) -30)))
    (terpri)))

実行結果

 0 番目: 0110100100100100101101101101101001001001001011011011011010010010
 1 番目: 0111000111000110001110001110011100011100011000111000111001110001
 2 番目: 0011110000111100001111000011110000111100001111000111100001111000
 3 番目: 1001011010010110100101101001001011010010110100101101001011010010
 4 番目: 1001100110011001100110010011001100110011001100110011001100110011
 5 番目: 0011111000001111000001111100001111100000111100000111110000111110
 6 番目: 1000011110000111100001111000011100001111000011110000111100001111
 7 番目: 0101010101010101010101010101010010101010101010101010101010101010
 8 番目: 1000110011100110001100111001100011000110011100110001100111001100
 9 番目: 1111100000000000000000000000000000000000000000000000000001111111
10 番目: 0100100110110110110010010011011011001001001001101101100100100110
11 番目: 0110011001100110011001100110011001100110011001100110011001100110
12 番目: 1101101101101101100100100100100100100110110110110110110110010010
13 番目: 1010010101010101010010101010101011010101010101011010101010101001
14 番目: 0101101010010101001010100101010010101001010100101010010101001010
15 番目: 1001010100101010010101011010101101010110101010010101001010100101
16 番目: 1010010110110100101101001001011010010110110100101101001001011010
17 番目: 0111111000000111111000000111111000000111111000000111111000000111
18 番目: 0110011001100110011001100110011001100110110011001100110011001100
19 番目: 1111111111110000000000000000000000000001111111111111111111111111
20 番目: 1111111000000001111111000000001111111000000001111111000000001111
21 番目: 1110000111100001111100001111000011111000011110000111110000111100
22 番目: 0000000000000000000000000000000111111111111111111111111111111111
23 番目: 0011100011100011000111000111000111000111001110001110001110001110
24 番目: 0000001111111111000000000111111111000000000111111111000000000011
25 番目: 0000000011111111111111111111110000000000000000000000111111111111
26 番目: 1000110001100011000110011100111001110011100111001110011000110001
27 番目: 0110001110011100011000111001110001100011100111000110001110011100
28 番目: 0000000000000000000000000000000000000000000000000000000000000000
29 番目: 0101011010101011010101001010101101010101101010100101010110101010
30 番目: 0010010110110110110100100100101101101101101001001001001011011011
31 番目: 1100000111110000011111000001111110000011111000001111100000111110
32 番目: 1100100110110010011011001001101101100100110110010011011001001101
33 番目: 1000111100011110000111000011100001111000111100011110001111000011
34 番目: 0100101101101101101101101101101101001001001001001001001001001011
35 番目: 1010100101011010110101001010110101101010010101101010010100101011
36 番目: 1100100110110011011001101100110110011011001001100100110010011001
37 番目: 0110011011001001101100110110010011011001001100100110110010011001
38 番目: 0000000000011111111111000000000011111111111000000000011111111111
39 番目: 1010010101101010110101001010110101001010110101001010110101001010
40 番目: 0011100110011001100110011001100011001100110011001100110011100110
41 番目: 1010101010010101010101010101010101010101010101010101010101010010
42 番目: 0001110001110001110001110001110001110001110011100011100011100011
43 番目: 0001111000001111000001111000001111000001111000001111000001111000
44 番目: 1100000000000000000011111111111111111110000000000000000001111111
45 番目: 0011000110011100110001100111001100011001110011000110011100111001
46 番目: 0001111111111100000000000011111111111000000000000111111111110000
47 番目: 0000111100011110001111000011100001110000111100011110001111000011
48 番目: 0000111100001111000011110000111100001111000011110000111100001111
49 番目: 1010010110100101101101001011010010110110100101101001011011010010
50 番目: 1100111001110001100011100111001110001100011000111001110011100011
51 番目: 0001110001110001110001110001110001110001110011100011100011100011
52 番目: 1111111000000000111111111000000000111111111100000000011111111100
53 番目: 1111111100000000011111111100000000011111111100000000111111111000
54 番目: 1100110011100110001100110001100111001100111001100011001100011001
55 番目: 1011010101011010101011010101011010101011010101011010101011010101
56 番目: 0001111000011110000111100001111000011111000011110000111100001111
57 番目: 0110011001100110110011001100110110011001100110110011001100110110
58 番目: 0000000001111111111111110000000000000001111111111111110000000000
59 番目: 0010010010010010010010010010010010010010010010010010010010010010
60 番目: 1000111000110001110001110011100011100111000111001110001110011100
61 番目: 0011100011100011100011100011100011100011100011100011100011100011
62 番目: 0101010110101010010101011010101001010101101010100101010110101010
63 番目: 0110001100111001110011100110001100011000110001100111001110011100
64 番目: 0110011001100110011001100110011001100110011001100110010011001100
65 番目: 0001110011100011100011100111000111000111000110001110001110001100
66 番目: 0011111100000011111100000001111110000001111110000001111110000001
67 番目: 0101010010101001010101101010100101010110101010010101011010101001
68 番目: 0101010101010010101010101010101010101010101010101010101010101010
69 番目: 1101101001011011010010010110110100100101101101001011011010010010
70 番目: 1101010101010101010110101010101010101010010101010101010101011010
71 番目: 1011010101010101010101010101010101010101010101010101010101010101
72 番目: 0110011100110001100011001110011100110001100011001110011000110001
73 番目: 1001011010010110100101101001010010110100101101001011010010110100
74 番目: 0111001110001110001110001110011100011100011100011100111000111000
75 番目: 0111100000111110000011111000001111100001111100000111110000011111
76 番目: 1000111000111000110001110001110001110001110001110011100011100011
77 番目: 0000111111100000000111111100000000111111110000000111111110000000
78 番目: 1000000000011111111111000000000001111111111100000000000111111111
79 番目: 0101001011010110100101001011010110100101001010010110101101001010
80 番目: 0001111000111100011110001111000111100011100001110000111000011100
81 番目: 1100000000111111110000000011111110000000011111111000000001111111
82 番目: 1100011110001110001110000111000111000011100011100001110001110000
83 番目: 0110110110110110100100100100100100100100100100100100100100100100
84 番目: 0100101101001010010110100101101001011010010100101101001011010010
85 番目: 1010101001010101010101010101010101010101010100101010101010101010
86 番目: 1101101101001001001001001001001001001001011011011011011011011011
87 番目: 0101010101010101010101010101010101010101010101010101010101010101
88 番目: 1010110100101001011010110100101001011010110100101001010010110101
89 番目: 0101010101010101010101010101101010101010101010101010101010101001
90 番目: 0010010010010010010010010010010010010010011011011011011011011011
91 番目: 1011001100100110011001001100110110011001101100110010011001100100
92 番目: 1000110001100011000110001100011000110001110011100111001110011100
93 番目: 0011011011001001001101101100100100110110110010010011011011001001
94 番目: 1111000011111000011111000001111000001111100001111100001111100000
95 番目: 0110110011001100110011001100110011001100110011001100110011001100
96 番目: 0100110110010011011001001101100100110110010011011001001101100100
97 番目: 0011100110001100011001110011100111001100011000110011100111001100
98 番目: 1001100100110010011001001100100110010011001101100110110011011001
99 番目: 1111111000000111111100000011111110000000111111000000011111100000

参考

x022235 commented 12 years ago

最下位ビット

   0 01001111010101100010111100000000101011101100101111000010110111
   4 11100010101010111000101110101000101011010011100110001001111110
   8 10000101010000011100001011111000100011101110001011011111000000
  12 11100111001000001111000001011010111011100110001001110101010000
   1 01001010100001001010101010100110110011010001101000000001110011
   5 11101100001110101010010100101101011001001011010100010011110101
   9 10101011001010100010000011010011000110001101110000100110000111
  13 10101100011101111001010101011101010010011100110011000010011011
   2 01001111011011011110000111011111101011000111101100101011101110
   6 10101011010010011001000100111001101110000010111110001001111001
  10 00010010000001000101010011101110001011000101001110001001100111
  14 10011011010000101110100110100000001110001001101011011110010110
   3 10010000011100101100001101011111001100100101110101100010011001
   7 10011111001100111011001100111110010001001001001101011101000011
  11 01010100011000000000010010010010001111100101110110100111100010
  15 00001001110110001101110010100001110100100111100010101101001011

上位ビット

 0 番目: 1001101010100111001010011100111010001101000011100001110010010101
 1 番目: 1100101000000111000011110110011001000010111000010000110101111110
 2 番目: 1101110111001010111110011111100100101011001101111000010100000000
 3 番目: 1111010000011111010011100001000001011111000001110101010111111000
 4 番目: 1010101010100110101011011011001110001011111010010100000010010110
 5 番目: 1010100100000100010110110010010101001110101010111011001001111100
 6 番目: 1001110110111100000010100011010110100110001000100101101000001001
 7 番目: 1101110010111001011011111000001110001001011001000101110111101011
 8 番目: 0001100100000100011000000010101011011111000000011100010101111110
 9 番目: 1000101000100101100111101000001011010111110101111001110101111111
10 番目: 1001011100001011011010011001110010100111010110110111111011000001
11 番目: 0000000011000100000010101001001001111010010000110011000101011001
12 番目: 0000111100001110111100001101011100010000100100011100011101000110
13 番目: 0010001010110101100011000111011101010100010011000111100010101111
14 番目: 1010010000100100011111100000110011000100111100000010101001010101
15 番目: 0100101000010101111000001011111101000101110100001001010000011101
16 番目: 1000000010011110000110010111111011011001111010111111101011100001
17 番目: 0100110010001010111111010010110110111010010110000101010111011001
18 番目: 0100001000010100100101100001110110010101101111000100100100110010
19 番目: 0111111110010011101110101101110011111011111000011101010010001100
20 番目: 1010100110110011111000000001000001011101100011100100100000111010
21 番目: 1010011100111110010110010110100100111111100010011100110011110010
22 番目: 1110011111101000001111110011111000101101101101000000100101001000
23 番目: 0110110010111101111110011001100101100110101011101011000001111110
24 番目: 1001101011011011111110111001001110111001101001000001100010000111
25 番目: 0000000001010111101001000100101010111010011011000011101000111001
26 番目: 1110101010111101010010101111001000100110100000101110001101101011
27 番目: 1010110111111001100110111000110010011111010101011000111000010101
28 番目: 0001011011100011000010011000111010100111000101110110101110100011
29 番目: 0110001101010000001000001110101111010110011101100001010110100111
30 番目: 0111101101111100010001010100100011110011100010001010011110001100
31 番目: 1101101100010001100111010110000111101010001110111000001010101010
32 番目: 0010001001101010111001111010000011000101010011000100110110100101
33 番目: 0011100000111101000100010000100000111000101101011010010011010101
34 番目: 1101000011000000011111101111011000110110001011110010100000000010
35 番目: 1001011111111111001011011001110010001100100011101000010111010010
36 番目: 1010011101101010000001101010010010010111110010000010001001001001
37 番目: 0100111001111100000000001110000100000100100000011101001100110101
38 番目: 1000111001100010001100001100110110000111000100100000111110101101
39 番目: 1111100101000010101100100011111011011111011101110110011110011001
40 番目: 1110010100111111001010111010001011100010011111010111010011011111
41 番目: 0100011001011101001101100000101001010000110101111010111001101111
42 番目: 1101101111010000110111110110011000010011010100000111010110000011
43 番目: 1101000111101100011111111110100110111101110011010100111110100101
44 番目: 0000101001000001101110001011111011101001100000010111010001101011
45 番目: 1000100001101110001111011001110100010000111100100000001010010101
46 番目: 1100101100100100011101101100110100000010100110000111010100111011
47 番目: 1100011101000011001101010111000000011100000100010100001011101000
48 番目: 0101001101000011000111011110101100101110100000011011101000000110
49 番目: 1101010100111110100000100100011101001011100010010100110001000011
50 番目: 1111001001001001001101101111110010000111101110111001100110011101
51 番目: 1111011011000100000100001011001000110100111010101011011000001101
52 番目: 0000100101110110100011101010011001110011100011011101000110111000
53 番目: 1100010010010000011000111110101010100100000001110001110011100000
54 番目: 1000011110111110101110001101000110010000111100011001110011001111
55 番目: 1011110010001000110110101001111101010111011011011010111101001110
56 番目: 1000110101100000110000111010110001010011001011010000111110000001
57 番目: 0011100010010110111101111111101111100110010101001010111100010100
58 番目: 0101000011101101101100100011000011000110110000000010000110001001
59 番目: 0100111000100101100100100110100110000000111010010111111001110000
60 番目: 0001101001101110000010110010101101011101101001010011100110101101
61 番目: 0100100000010010011111100010010011100101001001010001100111110011
62 番目: 1001101001110001110110111111110011111011000100001110101110010100
63 番目: 1100111101101101011110110011011001111100100000001110111100111110
64 番目: 0111110001110111010101001111001001110010000001000111000010111110
65 番目: 0011011111010101011000101111100101100100110111100001110010000100
66 番目: 1101011010011100101100011011010100101111111100111101111000000000
67 番目: 0010011011100011101111110101011101011001101000001100001111101111
68 番目: 0110001010011011011011101001100000100001011111110110011010000000
69 番目: 0010011001011111100011001110001110110100111011010101101100011001
70 番目: 1000111101011111000001100110011000010101001110111101010001010001
71 番目: 0010001001000101111001110101110110101001100101110111101011001100
72 番目: 1111101101011000111101111100100111111101110110010110000101110000
73 番目: 0111011100010011111000010110110100101001000110000101100001101010
74 番目: 1111110001101010101010001101111010110011010101001110000001110010
75 番目: 0001011011101010101001100100010000101010110101111001001001101001
76 番目: 1000010110000100001101001110100010110100011010100110110000111011
77 番目: 1000010000000110011001011001100111111110010010101111110011011100
78 番目: 1110110010110000011101110100000001000001000000111001010111101100
79 番目: 0110001000001011011010110110010010101100111101110110010000010001
80 番目: 0110001001100000110010010011000000010000000111010110000111000110
81 番目: 1111110110101101000001111011100110011011011001000111100010110111
82 番目: 0110111010101110101010011111111110000010110000101011010001011101
83 番目: 1001010001001111101001110111110111110011111010101110001110001101
84 番目: 1010110110000100011001010000000101101001100111001011000111111100
85 番目: 0011001010100101110011101110100010111100100100110111111110001110
86 番目: 0100011000101001110101010100110100010111001110101000110101011000
87 番目: 1000110101100000001001000101111100001110010010100000111110011001
88 番目: 1010000001000000010101000111100010000100100110000010001010100010
89 番目: 1001011011001110110110101101011100000101010110101000100111011101
90 番目: 1001100110100010010010101111011111100101101011000111101001000001
91 番目: 0011011001111011001111111110100010001110001111000011100101110001
92 番目: 0101111011010000011001111011111001110100010111000111001111111101
93 番目: 0001000110110111001001001000110000100110000001101100000111111000
94 番目: 0010000111110100000111000100111010011010111010011100001100110010
95 番目: 1101001100111010100100100000011111010101110101001100110001011101
96 番目: 0010110001010111110011001011100010100100101011001110111101110011
97 番目: 1000000000001000110001101111011001110101101100000101111000000111
98 番目: 0010111100010010010101011111010000100000100100111000100101101111
99 番目: 0110011001110010010000101101010111100000000111101010110101111001

x022235 commented 12 years ago

とりあえず実装したが取り込むにはもう少し考慮が必要。

互換性の考慮。(defvar *state* #.(make-random-state t)) をコンパイルするとバイトコードに古い random-state が S notation で書かれているが、その場合どうするか?
- ランダムステートの初期化オブジェクトが不正です というエラーになる
- make-random-state なんと誰も使っていないだろうから気にしない?
*features* に :mt19937ar とか追加する?
random の引数が少数の場合は genrand_real2 を呼ぶようにしたほうが良い
リファレンスに暗号論的擬似乱数生成器ではないことを書いておく
というか http://msdn.microsoft.com/ja-jp/library/bb982398.aspx を使えばよかったのかｗ

southly commented 12 years ago

ソースを取り込むのであればSFMTの方がいろいろ改善されていて良いんじゃないでしょうか。

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index-jp.html

それともMTの方が枯れていて安心という判断でしょうか。

x022235 commented 12 years ago

大した理由はないのですが、「このプログラムは2007年1月31日に一般に公開されました。それで、いくつかバグがあるかも知れません。」という記述が気になったのと、あと参考にした ruby が mt だったからですかね。

SFMT も 2007 年からもう更新されていないから問題なさそうですね。

southly commented 12 years ago

なるほど。

自分はMTの内部ベクトルの大きさが若干気になってまして(特に*scratch*でのプリント時)、 SFMTなら現状の内部ベクトルと同程度の大きさを選べるかなと思ったのでした。

x022235 commented 12 years ago

現状と同程度ということは MEXP=2281 あたりでビルドするということでしょうか。

sizeof(long)   * 55               = 220       ; 0.2.2.235
sizeof(w128_t) * (2281 / 128 + 1) = 288       ; SFMT

生成する擬似乱数の質を良くしようというのがこの issue の背景なんですが、 scratch でのプリント時のためだけに周期を短くするのは本末転倒だと思います。

make-random-state を使う人はあまりいないでしょうし、SBCL でも make-random-state を評価すると長々とプリントされるので特に気にする必要性はないと思います。

southly commented 12 years ago

アルゴリズムにしか着目していませんでしたが、擬似乱数のの質と言った場合周期もあわせて考えるの普通でしたね。

というわけで納得です。

xyzzy-022 / xyzzy