16番会議室「玉石混淆みんなで作るSample蔵」に寄せられたサンプル
"高速で高品位な乱数生成ユニット"
この発言に対し以下のコメントが寄せられています
#01201 ひの さん RE:高速で高品位な乱数生成ユニット
#01203 ひの さん 高速で高品位な乱数生成ユニット(D3版)
Mersenne Twister という新しく(といっても二年前)発明されたアルゴリ
ズムによる乱数生成ユニットです。
従来のアルゴリズムに比べて高速かつ高品位です。今後はこれが標準的な擬
似乱数生成アルゴリズムになっていくと思います。詳しくは、
http://www.math.keio.ac.jp/~nisimura/random/int/mt19937int.c
を参照してください。以下のソースはこのページにあるCのソースを元に移
植したものです(関数名は変更)。なお、私は Delphi 4 を持っていないので
動作チェックはFree Pascalで行い、オリジナル版と同じ出力になることを確
かめました。Delphi 3 以前は32ビット符号なし整数がサポートされていない
のでこのソースはこのままでは使えません(変更が必要)。
Free Pascalでは標準の乱数ユニットに比べて5倍高速でした(^^)v。
Delphi 3 では、標準の乱数ルーチンに比べて25%高速でした(^^)v
Delphi 4, Delphi 5 をお使いの方、ぜひベンチマークを取ってください。
なお、Delphi 4 で書かれたDLL版を、山本秀樹氏が
http://www.asahi-net.or.jp/~jz6h-ymmt/toolbox/mt.htm
で公開しています。
(*-----------------------------------------------------------------
Free Pascal unit for Mersenne Twister random number generator
Original source :
http://www.math.keio.ac.jp/~nisimura/random/int/mt19937int.c
Mersenne Twister Home Page :
http://www.math.keio.ac.jp/~matumoto/mt.html
Ported to Free Pascal unit by M. Hinoue (2000/11/27)
------------------------------------------------------------------*)
unit MTRand;
(* A C-program for MT19937: Integer version (1999/10/28) *)
(* genrand() generates one pseudorandom unsigned integer (32bit) *)
(* which is uniformly distributed among 0 to 2^32-1 for each *)
(* call. sgenrand(seed) sets initial values to the working area *)
(* of 624 words. Before genrand(), sgenrand(seed) must be *)
(* called once. (seed is any 32-bit integer.) *)
(* Coded by Takuji Nishimura, considering the suggestions by *)
(* Topher Cooper and Marc Rieffel in July-Aug. 1997. *)
(* This library is free software; you can redistribute it and/or *)
(* modify it under the terms of the GNU Library General Public *)
(* License as published by the Free Software Foundation; either *)
(* version 2 of the License, or (at your option) any later *)
(* version. *)
(* This library is distributed in the hope that it will be useful, *)
(* but WITHOUT ANY WARRANTY; without even the implied warranty of *)
(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *)
(* See the GNU Library General Public License for more details. *)
(* You should have received a copy of the GNU Library General *)
(* Public License along with this library; if not, write to the *)
(* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA *)
(* 02111-1307 USA *)
(* Copyright (C) 1997, 1999 Makoto Matsumoto and Takuji Nishimura. *)
(* When you use this, send an email to: matumoto@math.keio.ac.jp *)
(* with an appropriate reference to your work. *)
(* REFERENCE *)
(* M. Matsumoto and T. Nishimura, *)
(* "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform *)
(* Pseudo-Random Number Generator", *)
(* ACM Transactions on Modeling and Computer Simulation, *)
(* Vol. 8, No. 1, January 1998, pp 3--30. *)
(*------------*)
Interface
(*------------*)
(* initialize *)
procedure MTRandomSeed(Seed : Cardinal);
(* return real random number *)
function MTRandom : extended;
(* return integer random number [0..N-1] *)
function MTRandomInt(N : Longint) : Longint;
(*----------------*)
Implementation
(*----------------*)
const
(* Period parameters *)
N = 624;
M = 397;
MATRIX_A = $9908b0df; (* constant vector a *)
UPPER_MASK = $80000000; (* most significant w-r bits *)
LOWER_MASK = $7fffffff; (* least significant r bits *)
(* Tempering parameters *)
TEMPERING_MASK_B = $9d2c5680;
TEMPERING_MASK_C = $efc60000;
var
mti : integer;
mt : array[0..N-1] of Cardinal; (* the array for the state vector
*)
(*---------------------------------------
initialize the array with a seed
-----------------------------------------*)
procedure MTRandomSeed(seed : Cardinal);
var
i : integer;
begin
for i := 0 to N - 1 do begin
mt[i] := seed and $ffff0000;
seed := 69069 * seed + 1;
mt[i] := mt[i] or ((seed and $ffff0000) shr 16);
seed := 69069 * seed + 1;
end;
mti := N;
end;
(*----------------------------
return real random number
-----------------------------*)
function MTRandom : extended;
const
mag01 : array[0..1] of Cardinal = ($0, MATRIX_A);
(* mag01[x] = x * MATRIX_A for x=0,1 *)
var
y : Cardinal;
kk : integer;
begin
if mti >= N then begin (* generate N words at one time *)
for kk := 0 to N-M-1 do begin
y := (mt[kk] and UPPER_MASK) or (mt[kk+1] and LOWER_MASK);
mt[kk] := mt[kk+M] xor (y shr 1) xor mag01[y and $1];
end;
for kk := N - M to N-2 do begin
y := (mt[kk] and UPPER_MASK) or (mt[kk+1] and LOWER_MASK);
mt[kk] := mt[kk+(M-N)] xor (y shr 1) xor mag01[y and $1];
end;
y := (mt[N-1] and UPPER_MASK) or (mt[0] and LOWER_MASK);
mt[N-1] := mt[M-1] xor (y shr 1) xor mag01[y and $1];
mti := 0;
end;{if}
y := mt[mti];
Inc(mti);
y := y xor (y shr 11);
y := y xor ((y shl 7) and TEMPERING_MASK_B);
y := y xor ((y shl 15) and TEMPERING_MASK_C);
MTRandom := (y xor (y shr 18)) * 2.3283064365386962891e-10;
end;
(*---------------------------------------
return integer random number [0..N-1]
-----------------------------------------*)
function MTRandomInt(N : Longint) : Longint;
begin
MTRandomInt := Trunc(N * MTRandom);
end;
(*-------------------
initialize unit
---------------------*)
begin
MTRandomSeed(4357);
end.
2000/11/27(Mon) 10:18pm GFD03044 ひの
Original document by ひの 氏 ID:(GFD03044)
ここにあるドキュメントは NIFTY SERVEの Delphi Users' Forum の16番会議室「玉石混淆みんなで作るSample蔵」に投稿されたサンプルです。これらのサンプルはボーランド株式会社がサポートする公式のものではありません。また、必ずしも動作が検証されているものではありません。これらのサンプルを使用したことに起因するいかなる損害も投稿者、およびフォーラムスタッフはその責めを負いません。使用者のリスクの範疇でご使用下さい。
Copyright 1996-2002 Delphi Users' Forum
|