C言語による長周期のM系列乱数(Linear Feedback Shift Register)

M系列の乱数（最長線形帰還シフトレジスタ）を32ビットずつ出力するプログラムです。単純なM系 列 乱数でも，周期が長ければ，また，原始多項式の項の数が多ければ，Mersenne Twister のような高品質な乱数になり，Mersenne Twister のような複雑な理論が必要ないのではないかという疑問から作りました。独自のプログラムで見つけた，GF(2)上 の原始多項式を用いたもので，周期はMersenne Twisterと同様にメルセンヌ素数になります。Tempering 部分は整数乗算を用いた全単射にしていますが，最新のCPUには乗算専用の回路があるため若干遅くなる程度で，乱数の質は良くなりま す。

M系列乱数のコード(リアルタイム版)    mmm23209.c
GF(2)上の5項の原始多項式 によるM系列乱数のコードで，漸化式は an=an-23209+an-14389+an-8923+an-5491 です。周期はメルセンヌ素数 223209-1 で ，全周期では[23209/32]=725次元超立方体に均等分布します。漸化式の項の数を多くすることにより，ビットのかき混ぜを十分に行うと共に， Tempering 部分を整数乗算を用いた全単射にして良質の乱数を出力するようにしました。ワーキングメモリは Mersenne Twister の同周期版の2倍の726×2ワードです。コードには同周期で項の数が3～9の原始多項式をコメントアウトして埋め込んでありますから，速度を求めるか， 乱数の質を 求めるかに応じ て 簡単に変更できます。乱数の生成速度は， 原始多項式の項の数を少なくすれば Mersenne Twister より速くなります。また，零超過初期状態（初期状態に0が多いと、しばらくの間出力にも0が多くなる）からの回復は， 原始多項式の項の数 を 多くすれば，Mersenne Twister より早くなります。
乱数を1個ずつ生成しますので，応答の遅延はありません。

最速版    mmm19937.c
3項の原始多項式によるコードで，漸化式は an=an-19937+an-7083 です。 周期はメルセンヌ素数 219937-1 で，全周期では[19937/32]=623次元超立方体に均等分布します。ワーキングメモリは624×2ワードです。コードには他の周期の3項原始多項 式を コメントアウトして埋め込んであり，それに 簡単に変更できます。

超長周期版    mmm216091.c
5項の原始多項式によるコードで，漸化式は an=an-216091+an-133673+an-82613+an-51054 です。 周期はメルセンヌ素数 2216091-1≒1065050 で，全周期では[ 216091/32]=6752次元超立方体に均等分布します。ワーキングメモリは6753×2ワードです。SIMD-oriented Fast Mersenne Twister (SFMT) の最長周期のものと同周期で，世界最長周期の乱数と思われます。

3項M系列乱数のコード（省メモリ版）    ml_lfsr19937.c
3項の原始多項式によるコードで，漸化式は an=an-19937+an-7083 です。周期はメルセンヌ素数 219937-1 で ，全周期では[19937/32]=623次元超立方体に均等分布します。ワーキングメモリは624ワードで， 一気に624個の乱数を生成し，1つずつ出力しますが，mmm19937.c より遅くなります。コードには他の周期の3項原始多項式を コメントアウトして埋め込んであり，それに 簡単に変更できます。

ストリーム暗号用コード    mmm23209ang.c
mmm23209.c を改造して，ストリーム暗号用のコードを作りました。使用メモリは大きくなりましたが，シンプルな改造ですので，生成速度は元より若干遅くなった 程度です。

既約多項式探索用コード    hantei.c
既約多項式を探すために書いたコードです。

超立方体への均等分布の証明    distribution.pdf
M系列乱数の周期がメルセンヌ素数 2p-1 であるとき，32ビットずつ出力する乱数が，全周期では[p/32]次元超立方体に均等分布することの証明です。ここで，[x]はxを越えない最大の整数 を表します。

Tempering 部分について
以下のコードは，Tempering の効果を確かめるためのコードで，自然数の列の一部を Tempering だけによって，力ずくで乱数にするものです。このコードで発生させた乱数でも，NIST Special Publication 800-22rev1a の sts-2.1.1. の乱数テストをパスしました。速度とのバランスを考え，Tempering は3行にしていますが，行を増やせばさらに良質な乱数となります。
なお，0xF4B68C57UL などの係数は奇数ですから 232 と互いに素であり，X={n|nは整数，0≦n≦232-1} とするとき，Ｔempering はXからXへの全単射となっています。

#include <stdio.h>
int main(void)
{
    unsigned long i,y,buf[1],hajime=0xB3C64E6DUL;
    FILE *file = fopen("temp","wb");
    for ( i = hajime ; i<hajime+0x2000000; i++ ) {
        y=(i^(i>>16))*0xF4B68C57UL;
        y=(y^(y>>16))*0xC2A59E8DUL;
        y=(y^(y>>16))*0xB15D5363UL;
        buf[0]=y;
        fwrite(buf,sizeof(buf),1,file);
    }
    fclose(file);
    return 0;
}

実数の一様乱数
x をdouble型の変数として， 0≦x＜1で32ビット精度の乱数を出力するには通常は
    x = ransuu32()*(1.0/4294967296.0);
としますが，型変換を用いて
    unsigned long long a=ransuu32(),c;
    c = a<<20 ^ 0x3FF0000000000000ULL;
    x = (*(double*)&c)-1.0;
とすると若干速くなります。ここで，232=4294967296 です。

0≦x＜1で53ビット精度の実数の乱数を出力するには
    unsigned long long a=ransuu32(),b=ransuu32();
    x = (a<<21 ^ b)*(1.0/9007199254740992.0);
とすればよい。ここで，253=9007199254740992 で，(a<<21 ^ b) は2進法で53ビットの乱数を表します。

精度を52ビットに下げてよいなら
    unsigned long long a=ransuu32(),b=ransuu32(),c;
    c = (a<<20 ^ b) ^ 0x3FF0000000000000ULL;
    x = (*(double*)&c)-1.0;
とすれば高速です。

M系列乱数とGFSRの違い
M系列乱数とGFSRを同一視した記述を見かけますが，この2つは異なるものです。 GFSRは，M系列乱数を同時に複数走らせるもので，32ビットであれば32本のM系列乱数を1つに束ねたものです。「M系列乱数はビットが統計的に十分 にランダムでない」との批判も見られますが，これは，GFSRをM系列乱数と同一 視したための間違いであると思われます。

既約多項式
2n-1をメルセンヌ素数として， GF(2)上で xn+xk+1 が既約多項式となるものを示します。 2n-1が素数であるため，原始多項式にもなります。
xn+xk+1 が既約多項式であれば，xn+xn-k+1 も既約多項式なりますが，nが奇数であるため，k と n-k の一方は奇数，他方は偶数となります。したがって，kが奇数のものだけを示しました。

n     k
2     1
3     1
5     3
7     1     3
13    なし
17    3    5    11
19    なし
31    3    7    13    25
61     なし
89     51
107    なし
127     1     7     15     63    97
521     353   363   473    489
607     105   147   273
1279    861   1063
2203    なし
2281    715   915   1029
3217    67    2641
4253    なし
4423    271   369   649   1393   1419    2325   4053
9689    471   4187  7245  7853   9605
9941    なし
11213   なし
19937   881   7083  10095
21701   なし
23209   6619  9739  21679
44497   8575    23463
86243   なし
110503   53719   85273
132049   41469   52549   77595   98137   125049
216091   なし

ボックス・ミュラー法 (Box–Muller transform)    seiki.c
ボックス・ミュラー法は，[0，1) の実数の一様乱数から，正規分布（ガウス分布）に従う擬似乱数を生成するのに応用される方法です。　解説　逆関数法とボックス・ミュラー法
独自の高速な三角関数と対数関数のルーチンを用いることにより，math.hの組み込み関数を利用するときと比べて2倍ほど高速に正 規分布の乱数を発生できます。Box–Muller法の変形であるPolar法も組み込んでありますから，比較してみて下さい。
乱数は平均が0，標準偏差が1の正規分布に従 いますので，平均m，標準偏差sの正 規分布の乱数を得るには，s*ransuu_seiki()+m として下さい。
なお，カイ２乗分布，指数分布，ポアソン分布，二項分布も組み込んであります。

ご意見を，メールにて下記のアドレスにお寄せください。
MAIL.PNG - 671BYTES
(メールアドレスは，スパムメール対策のために画像になっています。)
最初 のページ に戻る