天泣記

線形合同法 (Linear Congruential Generator, LCG) は擬似乱数のアルゴリズムだが、下位ビットがあまり乱数になってないとか、良くないので使わない方がいいといわれている

なぜ良くないのか、ということを真面目に考えたことがなかったのだが、ふと考えてみた

線形合同法は以下によって生成される数列である

x{0} = seed
x{n+1} = (a * x{n} + c) % m

wikipedia の線形合同法の英語版によれば、 m は modulus, a は multiplier, c は increment という名前である

さて、x{n+1} とか x{n} と書くのは長いので、 y = x{n+1}, x = x{n} とおくことにする

y = (a * x + c) % m

ここで、よくある実装では m は 2^32 など、2のべき乗が選ばれる (wikipedia の線形合同法の英語版にはそういう例がたくさん載っている)

y = (a * x + c) % (2^n)

この場合、以下が成り立つ

y % 2 = (a * (x % 2) + c) % (2^n)

つまり、y の最下位ビット y % 2 は x の直前の状態 x の最下位ビット x % 2 だけに依存して決まる

ということは、状態が2種類しかないので、最下位ビットだけを取り出した数列は (定常状態としては) 0 -> 1 -> 0 -> 1 -> ... と 0 と 1 が交互に現れるか、 0 -> 0 -> ... と 0 がずっと続くか、 1 -> 1 -> ... と 1 がずっと続くか、という 3種類に限られる

つまり、最下位ビットが、最大周期2 の繰返しになるということで、これはよろしくない

2^n が悪いなら、3^n にすればいいのか、というと、

y = (a * x + c) % (3^n)
y % 3 = (a * (x % 3) + c) % (3^n)

となるので、3で割った余り、つまり、3進数における最下位桁が、直前の状態の最下位桁だけに依存して決まる、となる

そうすると、状態は3種類になるので、最大周期3 の繰返しになるので、これもよろしくない

では、(2^n)*(3^n) とか素因数をふたつ以上混ぜればいいのかというと、

y = (a * x + c) % (2^n)(3^n)
y % 2 = (a * (x % 2) + c) % (2^n)(3^n)
y % 3 = (a * (x % 3) + c) % (2^n)(3^n)

2進数として解釈したときの最下位桁が最大周期2の繰返しになり、かつ、 3進数として解釈したときの最下位桁が最大周期3の繰返しになる

つまりぜんぜん解決になっていない

一般に、y = (a * x + c) % m のとき、m を割りきれる d があると、 y % d = (a * (x % d) + c) % d が成り立つ

これは、Coq/SSReflect では以下のように証明できる (d %| m は d が m を割り切れる、という意味で、あと、%% は剰余で C の % に対応する)

From mathcomp Require Import all_ssreflect.

Goal forall a c m d x y, d %| m ->
  y = (a * x + c) %% m ->
  y %% d = (a * (x %% d) + c) %% d.
Proof.
  move=> a c m d x y Hdvdn ->.
  rewrite modn_dvdm; last by [].
  apply /eqP.
  by rewrite eqn_modDr modnMmr.
Qed.

ということで、modulus が 2^32 なら、2進数の最下位桁は周期が最大2になるし、下位2桁は周期が最大4になるし... 下位i桁は周期が最大 2^i になるし... となる

ということで、線形合同法の modulus は、 modulus を割り切る小さな数がないもの、つまり素数を選ぶことが望ましい

Wikipedia の線形合同法の日本語版の項には「Stephen K. Park と Keith W. Miller が、彼らのサーベイ中で「最低基準」として示したもの」として、

x{i+1} = (48271 * x{i}) % (2**31 - 1)

が示されており、この 2**31 - 1 = 2147483647 は素数である

wikipedia の線形合同法の英語版には、いろいろな線形合同法の実装で使われているパラメータがのっているが、 modulus はだいたいが 2のべき乗で、それ以外は 2**31 - 1 = 2147483647 が 3個、 2**3 * 7**5 = 134456 が 1個載っている

ところで、線形合同法における乱数の内部状態は直前に生成した乱数そのものなわけであるが、その状態は、32ビット整数とか、Cの整数型として普通に利用可能なものに保存されるだろう

また、乱数は、modulus 未満の整数である

とすると、(modulus が 2**n でない場合) 状態として表現可能な数の中には、乱数列に含まれない値が存在する。そういう値を種に設定したりすると何が起こるのだろうか、と考えると、まぁ、m は m を割りきるから、y = (a * (x % m) + c) % m であり、modulus 以上の値を指定することは、 mod m の値を指定することに対応するのだろう

あと気になるのは、modulus 未満の値を種にしても、周期が modulus でない場合はありうるという点である

ということで、

x{i+1} = (48271 * x{i}) % (2**31 - 1)

について、状態のなかにどういう列が入っているか (力ずくで) 調べてみた

% cat tst.c
#include <stdio.h>
#include <stdint.h>

uint32_t bitary[1 << (32-5)];
#define SETBIT(i) (bitary[(i) >> 5] |= 1UL << ((i) & 0x1f))
#define TESTBIT(i) ((bitary[(i) >> 5] & (1UL << ((i) & 0x1f))) != 0)

#define A 48271UL
#define C 0
#define MODULUS ((1UL << 31) - 1)

int main(int argc, char *argv)
{
  uint64_t i;
  for (i = 0; i < MODULUS; i++) {
    if (!TESTBIT(i)) {
      uint64_t j = i;
      uint64_t n = 0;
      printf("start=%lu ", i);
      fflush(stdout);
      while (!TESTBIT(j)) {
        SETBIT(j);
        n++;
        j = (j * A + C) % MODULUS;
      }
      printf("length=%lu end=%lu\n", n, j);
      fflush(stdout);
    }
  }
}
% gcc -O3 tst.c
% ./a.out
start=0 length=1 end=0
start=1 length=2147483646 end=1

つまり、状態が 0 というのが不動点で、これはずっと 0 で変わらないので乱数にならず、それ以外のすべての状態がひとつの乱数列をなしていることがわかる (まぁ、c = 0 なら 0 が不動点になるのは仕方ない)

また、Wikipedia に random0 として載っている

#define A 8121UL
#define C 28411UL
#define MODULUS 134456UL

というのを試すと、

start=0 length=134456 end=0

となったこれは 0..134455 の範囲がひとつの周期的な乱数列になっている

もちろん、134456 = 2**3 * 7**5 なので、この乱数は 2進数や7進数での下位桁が短い周期になっていてよろしくはない

では、modulus が素数であり、かつ、c が 0 でない例があるかと思って wikipedia の線形合同法の英語版を探すと、RtlUniform from Native API というのがあげられていた

#define A 2147483629UL
#define C 2147483587UL
#define MODULUS ((1UL << 31) - 1)

これを実行すると、以下のようになった

start=0 length=715827882 end=0
start=1 length=715827882 end=1
start=5 length=715827882 end=5
start=1243280003 length=1 end=1243280003

つまり、0..(MODULUS-1) という中に、4つの周期数列が含まれていて、 3個は周期 715827882 で、1個は周期 1 というわけである

まぁ、好んで使うものではない気がする

では、c = 0 だと種として 0 を選ぶと乱数でなくなってしまうので c を 0 以外に選び、下位桁の周期性を避けるために m を素数として選び、さらに周期が m になる、というのはあり得るだろうか、と考えたが、これはうまくいかないようだ

wikipediaには、 c != 0 の場合で周期を m になる場合について、 m が square-free integer の場合は a = 1 にならざるを得ない、と書いてある。

square-free integer というのは素因数分解したときに、同じ素数は高々1つしか現れないというやつで、素数は square-free integer である。

ということで、a = 1 にならざるを得ないということで、y = (x + c) % m になるわけだが、これは乱数としてはまったくよろしくない

天泣記

2020-12-30 (Wed)