| history | TOP |
2005/12/28:作成
| overview | TOP |
最初に。このページの内容に独自性はありません。完全無欠にパクリです。
翻訳もそのまま、本当にそのままやってるんでコメントなどもそのままです。おまけに思いっきりDelphi版の実装も見つけてしまったりして...
しかもそのソースを覗き見して
{$OVERFLOWCHECKS OFF}
{$RANGECHECKS ON}
をソースに追加してます。やっぱりオーバーフローするよなぁというのもコレを見たことで胸をなでおろしてます。
| コード | TOP |
(*
===========================================================================
mt19937ar.c
mt19937ar.cの使い方 （2002/1, 2004/2改訂）
(略)
===========================================================================
*)
unit mt19937ar;
{$OVERFLOWCHECKS OFF}
{$RANGECHECKS ON}
interface
function GenRandInt32(): Cardinal;
function GenRandInt31(): Integer;
function GenRandReal1(): Double;
function GenRandReal2(): Double;
function GenRandReal3(): Double;
function GenRandRes53(): Double;
implementation
// Period parameters
const
MERSENNE_TWISTER_N = 624;
MERSENNE_TWISTER_M = 397;
MERSENNE_TWISTER_MATRIX_A = $9908b0df; // constant vector a
MERSENNE_TWISTER_UPPER_MASK = $80000000; // most significant w-r bits
MERSENNE_TWISTER_LOWER_MASK = $7fffffff; // least significant r bits
F_MASK = $FFFFFFFF;
var
// the array for the state vector (unsigned long)
Mt: array[0..MERSENNE_TWISTER_N-1] of Cardinal;
// mti==N+1 means mt[N] is not initialized
Mti: Cardinal = MERSENNE_TWISTER_N + 1;
//-----------------------------------------------------------------------------
// initializes mt[N] with a seed (unsigned long)
procedure InitGenRand(const Seed: Cardinal);
begin
Mt[0] := Seed and F_MASK;
Mti := 1;
while (Mti < MERSENNE_TWISTER_N) do
begin
Mt[Mti] := (1812433253 * (Mt[Mti-1] xor (Mt[Mti-1] shr 30)) + Mti);
// See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
// In the previous versions, MSBs of the seed affect
// only MSBs of the array mt[].
// 2002/01/09 modified by Makoto Matsumoto
Mt[Mti] := Mt[Mti] and F_MASK;
// for >32 bit machines
Inc(Mti);
end;
end;
//-----------------------------------------------------------------------------
function Max(const A, B: Integer): Integer;
begin
if (A >= B) then Result := A else Result := B;
end;
//-----------------------------------------------------------------------------
// initialize by an array with array-length
// init_key is the array for initializing keys
// key_length is its length
// slight change for C++, 2004/2/26
procedure InitByArray(const Seeds: array of Cardinal);
var
i, j, k: Cardinal;
SeedsLen: Cardinal;
begin
InitGenRand(19650218);
i := 1;
j := 0;
SeedsLen := Length(Seeds);
k := Max(SeedsLen, MERSENNE_TWISTER_N);
while (k > 0) do
begin
Mt[i] := (Mt[i] xor ((Mt[i-1] xor (Mt[i-1] shr 30)) * 1664525))
+ Seeds[j] + j; // non linear
Mt[i] := Mt[i] and F_MASK; // for WORDSIZE > 32 machines
Inc(i);
Inc(j);
if (i >= MERSENNE_TWISTER_N) then
begin
Mt[0] := Mt[MERSENNE_TWISTER_N-1];
i := 1;
end;
if (j >= SeedsLen) then
begin
j := 0;
end;
Dec(k);
end;
for k := MERSENNE_TWISTER_N - 1 downto 1 do
begin
Mt[i] := (Mt[i] xor ((Mt[i-1] xor (Mt[i-1] shr 30)) * 1566083941))
- i; // non linear
Mt[i] := Mt[i] and F_MASK; // for WORDSIZE > 32 machines
Inc(i);
if (i >= MERSENNE_TWISTER_N) then
begin
Mt[0] := Mt[MERSENNE_TWISTER_N-1];
i := 1;
end;
end;
// MSB is 1; assuring non-zero initial array
Mt[0] := $80000000;
end;
//-----------------------------------------------------------------------------
// generates a random number on [0,0xffffffff]-interval
function GenRandInt32(): Cardinal;
const
// mag01[x] = x * MATRIX_A for x=0,1
MAG01: array[0..1] of Cardinal = (0, MERSENNE_TWISTER_MATRIX_A);
var
kk: Integer;
begin
// generate N words at one time
if (Mti >= MERSENNE_TWISTER_N) then
begin
// if init_genrand() has not been called,
if (Mti = MERSENNE_TWISTER_N + 1) then
begin
// a default initial seed is used
InitGenRand(5489);
end;
kk := 0;
while (kk < MERSENNE_TWISTER_N - MERSENNE_TWISTER_M) do
begin
Result := (Mt[kk] and MERSENNE_TWISTER_UPPER_MASK)
or (Mt[kk+1] and MERSENNE_TWISTER_LOWER_MASK);
Mt[kk] := Mt[kk + MERSENNE_TWISTER_M]
xor (Result shr 1)
xor MAG01[Result and 1];
Inc(kk);
end;
while (kk < MERSENNE_TWISTER_N - 1) do
begin
Result := (Mt[kk] and MERSENNE_TWISTER_UPPER_MASK)
or (Mt[kk+1] and MERSENNE_TWISTER_LOWER_MASK);
Mt[kk] := Mt[kk + (MERSENNE_TWISTER_M - MERSENNE_TWISTER_N)]
xor (Result shr 1)
xor MAG01[Result and 1];
Inc(kk);
end;
Result := (Mt[MERSENNE_TWISTER_N-1] and MERSENNE_TWISTER_UPPER_MASK)
or (Mt[0] and MERSENNE_TWISTER_LOWER_MASK);
Mt[MERSENNE_TWISTER_N-1] := Mt[MERSENNE_TWISTER_M-1]
xor (Result shr 1)
xor MAG01[Result and 1];
Mti := 0;
end;
Result := Mt[Mti];
Inc(Mti);
// Tempering
Result := Result xor (Result shr 11);
Result := Result xor ((Result shl 7) and $9d2c5680);
Result := Result xor ((Result shl 15) and $efc60000);
Result := Result xor (Result shr 18);
end;
//-----------------------------------------------------------------------------
// generates a random number on [0,0x7fffffff]-interval
function GenRandInt31(): Integer;
begin
Result := (GenRandInt32() shr 1);
end;
//-----------------------------------------------------------------------------
// generates a random number on [0,1]-real-interval
function GenRandReal1(): Double;
begin
// divided by 2^32-1
Result := (GenRandInt32() * (1.0 / 4294967295.0));
end;
//-----------------------------------------------------------------------------
// generates a random number on [0,1)-real-interval
function GenRandReal2(): Double;
begin
// divided by 2^32
Result := (GenRandInt32() * (1.0 / 4294967296.0));
end;
//-----------------------------------------------------------------------------
// generates a random number on (0,1)-real-interval
function GenRandReal3(): Double;
begin
// divided by 2^32
Result := (GenRandInt32() + 0.5) * (1.0 / 4294967296.0);
end;
//-----------------------------------------------------------------------------
// generates a random number on [0,1) with 53-bit resolution
function GenRandRes53(): Double;
var
a, b: Cardinal;
begin
a := GenRandInt32() shr 5;
b := GenRandInt32() shr 6;
Result := (a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
end;
//-----------------------------------------------------------------------------
const
DefaultSeeds: array[0..3] of Cardinal = ($123, $234, $345, $456);
initialization
//InitGenRand(19650218);
InitByArray(DefaultSeeds);
end.
| 使い方 | TOP |
こんな感じで。これは本家が乱数出力サンプル(mt19937ar.out)を試験用に出しているんですが、それと同じ内容を出力するためのコードです。
uses
mt19937ar;
//-----------------------------------------------------------------------------
// These real versions are due to Isaku Wada, 2002/01/09 added
procedure TForm1.Button1Click(Sender: TObject);
var
i: Integer;
c: Cardinal;
r: double;
FileName: String;
FileStream: TFileStream;
Line: String;
begin
FileName := ChangeFileExt(ParamStr(0), '.txt');
FileStream := TFileStream.Create(
FileName,
fmCreate or fmShareDenyWrite);
for i := 0 to 999 do
begin
c := GenRandInt32();
Line := Format('%10u ', [c]);
FileStream.Write(Line[1], Length(Line));
if (i mod 5 = 4) then
begin
Line := #13#10;
FileStream.Write(Line[1], Length(Line));
end;
end;
for i := 0 to 999 do
begin
r := GenRandReal2();
Line := Format('%10.8f ', [r]);
FileStream.Write(Line[1], Length(Line));
if (i mod 5 = 4) then
begin
Line := #13#10;
FileStream.Write(Line[1], Length(Line));
end;
end;
FileStream.Free;
Application.MessageBox(PChar(FileName + 'に出力完了'), PChar('情報'),
MB_OK or MB_ICONINFORMATION);
end;
| 実行結果・ダウンロード | TOP |
1067595299 955945823 477289528 4107218783 4228976476 (990個の整数乱数を略) 2643151863 3896204135 2416995901 1397735321 3460025646 0.76275443 0.99000644 0.98670464 0.10143112 0.27933125 (990個の実数乱数を略) 0.33163422 0.85739792 0.59908488 0.74566046 0.72157152
ちゃんとあってるようです。
20051227MersenneTwisterTest.zip(167,509bytes)
| EOF | TOP |