[Delphi][小ネタ] メルセンヌ・ツイスタ

はじめに

９日、１１日に引き続き乱数の話です。

今回はメルセンヌ・ツイスタについてです。

メルセンヌ・ツイスタとは

Wikipedia によると

メルセンヌ・ツイスタ (Mersenne twister、通称MT) は擬似乱数列生成器の1つである。
既存の疑似乱数列生成手法にある多くの欠点がなく、高品質の疑似乱数列を高速に生成できる。考案者らによる実装が修正BSDライセンスで公開されている。

とのことで、周期がメッチャ長い（$2^{19937} - 1$）、生成される乱数の質が高い、などの特徴があります。
デメリットは計算資源を結構使う（メモリ・CPU）事です。

実装

Delphi での実装はこんな感じです。
長いので折りたたみます。

メルセンヌ・ツイスタ

unit PK.Math.MersenneTwister;

{$OVERFLOWCHECKS OFF}

interface

uses
  System.SysUtils;

type
  TMersenneTwister = record
  private const
    N = 624;
    M = 397;
    MATRIX_A =   $9908B0df;
    UPPER_MASK = $80000000;
    LOWER_MASK = $7fffffff;
  private class var
    FMatrix: array[0.. N - 1] of UInt32;
    FIndex: Integer;
  private
    class procedure Twist; static;
  public
    class procedure Init; static;
    class procedure SetSeed(const ASeed: UInt32); static;
    class function Execute: UInt32; static;
  end;

implementation

{ TMersenneTwister }

class function TMersenneTwister.Execute: UInt32;
begin
  if FIndex >= N then
    Twist;

  var Y := FMatrix[FIndex];
  Inc(FIndex);

  // Tempering
  Y := Y xor (Y shr 11);
  Y := Y xor ((Y shl 7) and $9d2c5680);
  Y := Y xor ((Y shl 15) and $efc60000);
  Y := Y xor (Y shr 18);

  Result := Y;
end;

class procedure TMersenneTwister.Init;
begin
  RandSeed := 5489; // Default Seed (mt19937)
  SetSeed(RandSeed);
end;

class procedure TMersenneTwister.SetSeed(const ASeed: UInt32);
begin
  FMatrix[0] := ASeed;

  for var i := 1 to N - 1 do
    FMatrix[i] := 1812433253 * (FMatrix[i - 1] xor (FMatrix[i - 1] shr 30)) + i;

  FIndex := N;
end;

class procedure TMersenneTwister.Twist;
begin
  for var i := 0 to N - M - 1 do
  begin
    var Y := (FMatrix[i] and UPPER_MASK) or (FMatrix[i + 1] and LOWER_MASK);
    FMatrix[i] := FMatrix[i + M] xor (Y shr 1) xor (MATRIX_A * (Y and 1));
  end;

  for var i := N - M to N - 2 do
  begin
    var Y := (FMatrix[i] and UPPER_MASK) or (FMatrix[i + 1] and LOWER_MASK);
    FMatrix[i] := FMatrix[i + (M - N)] xor (Y shr 1) xor (MATRIX_A * (Y and 1));
  end;

  var Y := (FMatrix[N - 1] and UPPER_MASK) or (FMatrix[0] and LOWER_MASK);
  FMatrix[N - 1] := FMatrix[M - 1] xor (Y shr 1) xor (MATRIX_A * (Y and 1));

  FIndex := 0;
end;

end.

サンプル

前回、前々回と同じサンプルを動かしてみました。

結果はこうなりました。

0 975
1 1020
2 1027
3 993
4 969
5 968
6 953
7 990
8 1077
9 1029

特に代わり映えはしないですね。

まとめ

一般的なプログラムでメルセンヌ・ツイスタ程の強度を持つ乱数が必要な事はほとんど無いでしょう。
計算資源も多く使うため、使い所が難しいかもしれません。
また、メルセンヌ・ツイスタも暗号用には使えません。

反面、科学技術計算などでは周期の長さが役に立つかも知れません（モンテカルロ法を試す時とか）

乱数はゲーム業界では必須の機能で、機能に応じてアルゴリズムを使い分けたりしますが、SIer が手がける分野のプログラムではあまり使用しないと思います（宝くじぐらい？）。
そのため、紹介記事を書いてみました。