バイナリ行列演算ライブラリの開発と Xorshift 疑似乱数生成器のパラメータ計算

Development of binary matrix calculation library and calculation of Xorshift RNG parameter

Abstract

  バイナリ行列※1は，0 と 1 の 2 値しかとらない行列で，xorshift 疑似乱数生成器の位数計算などに用いられる．しかしながら，高速にバイナリ行列を計算する専用のライブラリを確認できなかったため開発した．char 型の配列を利用した計算と比較して，1/8 のメモリで済み，bit 演算により バイナリ行列乗算を 24~ 600 倍程度高速に計算できると推測される※2．ライブラリとしての利便性を考慮して，演算子のオーバーロードを行うとともに，サンプルプログラムとして xorshift 疑似乱数生成器の位数を計算した．

※1. Wikipedia では，"logical matrix, binary matrix, relation matrix, Boolean matrix, or (0,1) matrix" の名称で紹介されている[8]．bits martix の名称は一般的でないと思われる．
※2. char 型の配列を利用したバイナリ行列演算を実装していないため，直接の比較はしていないが，利用した手法 [3] から，32x32 bits で約 24 倍，64x64 bits で約 60 倍，1024x1024 bits で約 120 倍，4096x4096 bits で約 590 倍，高速に計算できると推測される．

Introduction

  バイナリ行列演算は，情報理論の分野で利用されている．疑似乱数生成器では，位数と呼ばれる疑似乱数生成器がリセットされてから始めの状態に戻るまでの周期を計算するために，バイナリ行列演算を利用している．参考文献 [1] [2] より，xorshift random number generator (xorshift RNG) の位数を計算するためには，多くのバイナリ行列乗算を必要とすることが分かる．しかしながら，バイナリ行列演算専用のライブラリを確認できなかったため，高速なバイナリ行列乗算が可能なバイナリ行列演算ライブラリの開発を行う．

Previous research

■ 高速なバイナリ行列乗算

  高速なバイナリ行列乗算を実装するにあたり，参考文献 [3] の手法を利用した．(一部計算が一致しなかったため，実装を修正している)．提案された手法では，行列を 8x8 のブロックに分割して計算を行うことで，高速にバイナリ行列乗算を可能としている．ブロックは，8x8 であるため，実際に確保されるメモリサイズは，64 の倍数 bits となる．バイナリ行列の行サイズ・列サイズが 8 の倍数とならない場合は，8 の倍数となる最小サイズのブロックを確保して計算され，不要な bit は破棄される．このため，計算自体は任意のサイズで可能である．

■ 転置

  転置を実装するにあたり，参考文献 [4] [5] の分割統治法による手法を用いた．この手法では 8x8 のバイナリ行列に対して，数回の bit 演算で転置を行うことができる．

■ Xorshift RNG の位数計算

  疑似乱数生成器である Xorshift が George Marsaglia 氏 により提案され [1]，参考文献 [2] において，その数理がより平易に解説されている． これらより，位数が $2^{32}-1$ となる条件を調べ上げるには，$2^{32}-1$ でのみ位数 $T$ が単位バイナリ行列 $I$ に戻り，それ以外の $2^{32}-1$ の約数全てにおいて，$T\neq I$ を満たす必要がある．本投稿では，具体的に，[1] の page.2 下部に示されたパラメータ表の再現を行う．これは，[2] の「理解度確認テスト：優良パラメータの計算」に相当する．

■ 素数の計算

約数を計算するにあたり，まずは素数をエラトステネスの篩によって計算した．実装は [6] を参考にした．

■ 約数の計算

実装は [7] を参考にした．

Structure of programs

  データ構造は，uint64 配列 1 つを 8x8 のバイナリ行列としており，8x8 部分が 行優先，配列部分が 列優先として実装されている．

Implementation

  位数計算のサンプルプログラムを示す．このプログラムは，[1] の page.2 下部に示されたパラメータ表を再現する．全ての実装は，GitHub のリポジトリに示す．ここで利用している sstd は，私用のライブラリで，GitHub において開発を行っており，https://sstd-lib.blogspot.jp/ (/bitmatrix.html) に簡単なドキュメントを用意している．(サポートはない)

#define use_sstd_measureTime
#include <sstd/sstd.hpp>

bool isFullPeriod(uint N, uint a, uint b, uint c){

    //std::vector<uint64> divs = sstd::divisor(4294967296-1); sstd::printn(divs);
    // Below line returns the same result of this line, but took a little 
    // time to run everytime and takes a little heavy memory (about 2 GByte).
    std::vector<uint64> divs = {1, 3, 5, 15, 17, 51, 85, 255, 257, 771, 1285, 3855, 4369, 13107, 21845, 65535, 65537, 
    196611, 327685, 983055, 1114129, 3342387, 5570645, 16711935, 16843009, 50529027, 84215045, 252645135, 
    286331153, 858993459, 1431655765, 4294967295};

    sstd::bmat I  = sstd::eye(N, N);
    sstd::bmat La = sstd::LxShiftMat(N, a);
    sstd::bmat Rb = sstd::RxShiftMat(N, b);
    sstd::bmat Lc = sstd::LxShiftMat(N, c);
    sstd::bmat T = (I + La)*(I + Rb)*(I + Lc);
//  sstd::printn(T);
    
    for(uint i=0; i<divs.size()-1; i++){
        sstd::bmat Tp = T^(divs[i]); // XORSHIFT
        if(Tp==I){ return false; }
    }
    sstd::bmat Tp = T^(divs[divs.size()-1]); // XORSHIFT
    if(Tp==I){ return true;
    }  else  { return false; }
}
int main(){
    printf("This code calculates the table of \"George Marsaglia [2003], Xorshift RNGs, page.2\".\n");
    
    printf("■ measureTime_start---------------\n\n"); time_m timem; sstd::measureTime_start(timem);

    uint N = 32;
    uint num=0;
    printf("  a,  b,  c\n");
    for(uint a=0; a<N; a++){
        for(uint b=0; b<N; b++){
            for(uint c=0; c<N; c++){
                if(a<c && isFullPeriod(N, a, b, c)){
                    printf("|%2u, %2u, %2u", a, b, c);
                    num++;
                    if(num%9==0){ printf("|\n"); }
                }
            }
        }
    }
    printf("\n");

    printf("■ measureTime_stop----------------\n"); sstd::measureTime_stop(timem);
    return 0;
}

Results/Benchmarks

実行結果

$ ./exe
This code calculates the table of "George Marsaglia [2003], Xorshift RNGs, page.2".
■ measureTime_start---------------

  a,  b,  c
| 1,  3, 10| 1,  5, 16| 1,  5, 19| 1,  9, 29| 1, 11,  6| 1, 11, 16| 1, 19,  3| 1, 21, 20| 1, 27, 27|
| 2,  5, 15| 2,  5, 21| 2,  7,  7| 2,  7,  9| 2,  7, 25| 2,  9, 15| 2, 15, 17| 2, 15, 25| 2, 21,  9|
| 3,  1, 14| 3,  3, 26| 3,  3, 28| 3,  3, 29| 3,  5, 20| 3,  5, 22| 3,  5, 25| 3,  7, 29| 3, 13,  7|
| 3, 23, 25| 3, 25, 24| 3, 27, 11| 4,  3, 17| 4,  3, 27| 4,  5, 15| 5,  3, 21| 5,  7, 22| 5,  9,  7|
| 5,  9, 28| 5,  9, 31| 5, 13,  6| 5, 15, 17| 5, 17, 13| 5, 21, 12| 5, 27,  8| 5, 27, 21| 5, 27, 25|
| 5, 27, 28| 6,  1, 11| 6,  3, 17| 6, 17,  9| 6, 21,  7| 6, 21, 13| 7,  1,  9| 7,  1, 18| 7,  1, 25|
| 7, 13, 25| 7, 17, 21| 7, 25, 12| 7, 25, 20| 8,  7, 23| 8,  9, 23| 9,  5, 14| 9,  5, 25| 9, 11, 19|
| 9, 21, 16|10,  9, 21|10,  9, 25|11,  7, 12|11,  7, 16|11, 17, 13|11, 21, 13|12,  9, 23|13,  3, 17|
|13,  3, 27|13,  5, 19|13, 17, 15|14,  1, 15|14, 13, 15|15,  1, 29|17, 15, 20|17, 15, 23|17, 15, 26|

■ measureTime_stop----------------
--------------------------------
 Execution time:     9. 465 sec
--------------------------------

実行環境:
・Microsoft Windows 8.1 x64
・Cygwin, GCC 5.4.0
・Intel Core i5-5200U CPU @ 2.20GHz
・DDR3 SDRAM 4.0GB

Conclusion

  速度比較を自ら行っていないものの，成果物自体は実用的な計算速度を示すことが分かった．また，実装としても，演算子をオーバーロードすることで，試行錯誤が必要な部分を，非常に短いコード量で記述できることを示した．

Future work

  本ライブラリを利用して，[1] の残りの内容を再現する．

Acknowledgments.

  本記事は，Marsaglia, George 氏の論文 [1] を紹介した びりある 氏の記事 [2] に触発され執筆した．彼らと，その他参考文献に掲げた記事の著者に感謝申し上げる．

References

  [1] Marsaglia, George. 2003. Xorshift RNGs. Journal of Statistical Software
  [2] Google Chromeが採用した、擬似乱数生成アルゴリズム「xorshift」の数理
  [3] 高速なビット行列乗算
  [4] 8x8のビット行列を転置する - エンジニアのソフトウェア的愛情
  [5] ビット演算＞ビットマップ 転置 | 津田の開発な日記
  [6] Sim‏@Sim0000 ツイート 2017年5月9日 - 1:34, https://ideone.com/LsbblP
  [7] M.Hiroi's Home Page, Puzzle DE Programming, 約数 (divisor)
  [8] Logical matrix - Wikipedia - 2017年12月2日閲覧