今ではもう利用されなくなった古典的な乱数生成アルゴリズムに、二乗中抜き法(middle-square method)というのがあります。平方採中法とも訳されます。
かのノイマンが1946年ごろに発明したとされる、世界初の疑似乱数生成アルゴリズムです。

使われなくなったのは欠点があるかららしいのですが、この記事では欠点についてTAOCP(The Art of Computer Programming)とその演習問題を通して見ていきます。

※ この記事のお題は、The Art of Computer Programming Volume 2 Seminumerical Algorithmsの最初の方で取り上げられている話題です。

二乗中抜き法

アルゴリズムとしてはとてもシンプルです。

適当な初期値を定め、最初の値とする。値に相応の桁数も定めておく。
次の値は、今の値を二乗して、中央から同じだけの桁数を取ったものとする（桁数が足りない場合は上位桁に0を詰める）

実装も簡単です。乱数生成器をRubyでざっと書いてみました。

def msm_enum(x0, digits = x0.to_s(10).bytesize)
  Enumerator.new do |yielder|
    x = x0
    loop {
      yielder << x
      # 二乗して文字列化する。桁数がちょうど倍になるように（必要なら）頭にゼロ埋めする
      square = sprintf("%0#{digits * 2}d", x ** 2)
      x = square[(digits/2) ... (digits/2 + digits)].to_i
    }
  end
end

実行例です。

>> msm_enum(1763).take(20)
=> [1763, 1081, 1685, 8392, 4256, 1135, 2882, 3059, 3574, 7734, 8147, 3736, 9576, 6997, 9580, 7764, 2796, 8176, 8469, 7239]

何がまずいのだろうか？

実行例はそれらしい乱数列が生成されており、そこまで悪くないようにも見えます。

別の初期値で動かしてみます。

>> msm_enum(1010).take(20)
=> [1010, 201, 404, 1632, 6634, 99, 98, 96, 92, 84, 70, 49, 24, 5, 0, 0, 0, 0, 0, 0]

途中から値が0だけになってしまっています。このアルゴリズムは、前の値から二乗によって次の値を生成します。ひとたび0になると、0の二乗は0となり、それに中抜き操作をしても0というわけです。これは良くない性質です。

さて、このようにゼロが続く数列を生み出してしまう場合はどれぐらいあるのか調べます。

require 'set'

def detect_zero(enum)
  set = Set.new
  while (value = enum.next) != 0
    return false if set.include?(value)
    set << value
  end
  true
end

detect_zeroは、「周期が見つかるか、ゼロが見つかる」まで乱数を生成してみるというメソッドです。ゼロが見つかったかどうかで真偽値を返します。

これを使って2桁以下の数すべての挙動を見てみます。

>> (0 .. 99).map {|x| detect_zero(msm_enum(x, 2)) }.count(true)
=> 62
>> (0 .. 99).map {|x| detect_zero(msm_enum(x, 2)) }.count(false)
=> 38

半数以上の数がゼロになってしまうことが分かります。

ゼロになる場合の解析例

ゼロになる問題について、TAOCPには少し異なる切り口からの演習問題も載っていました。

基数$b$で$2n$桁の数を使う二乗中抜き法には、数列内に上位$n$桁が0の数を含む場合、その後の数は次々と小さくなり、最終的に0を繰り返すという欠点があることを証明せよ

こんな数$X_i$が出てきたときのことを考えよ、と言われました。

X_i = \overbrace{0\ 0\ \cdots\ 0}^{n桁}\ \overbrace{b_n\ b_{n-1}\ \cdots\ b_0}^{n桁}

※ ここでは下$n$桁の各桁を、添字を付けて$b_i$と表記しています。

これを二乗して中抜きすることを考え、次の数$X_{i+1}$がどうなるか考察します。
出発点は、$X_i < b^n$です。この関係は、$X_i$が$n$桁の基数$b$の数であることから得られます。

二乗します。$X_i$は正数なので大小関係は変わりません。

X_i^2 < (b^n)^2 = b^{2n}

$X_i^2$を中抜きします。

元が$2n$桁なので2乗すると$4n$桁になっています
中央の$2n$桁を取り出すので、下$n$桁を切り落とします
- 上$n$桁は0になっていることに注意してください（つまり忘れて構いません）

切り落とすには$b^n$で割って床関数を付ければ良く、これが次の数$X_{i+1}$にあたるので

X_{i+1} = \lfloor X_i^2/b^n \rfloor \le X_i^2/b^n

です。最後の変形は一般に$\lfloor a \rfloor \le a$であることを使っています。

さて、$X_i < b^n$だったことを利用して、$X_i^2/b^n$を評価すると

X_i^2/b^n < X_i b^n / b^n = X_i

となるので、まとめると

X_{i+1} \le X_i^2/b^n < X_i

したがって

X_{i+1} < X_i

つまり、乱数列の値がどんどん小さくなっていくというわけです。少なくともこういう場合がある以上は、二乗中抜き法はあまり良いアルゴリズムとは言えなさそうですね……。

まとめ

二乗中抜き法について紹介しました
アルゴリズムについて簡単に解析し、特にゼロの発生パターンについて見てみることで、アルゴリズムの欠点を調べました
- このやり方はTAOCPの練習問題ほとんどそのままですが、ある程度手抜いています