ab.com コミュニティ

ActiveBasicを通したコミュニケーション
現在時刻 - 2017年9月22日(金) 11:47

All times are UTC+09:00




新しいトピックを投稿する  トピックへ返信する  [ 4 件の記事 ] 
作成者 メッセージ
投稿記事Posted: 2005年12月02日(金) 03:31 
オフライン

登録日時: 2005年5月31日(火) 09:34
記事: 141
住所: 北海道
Mersenne Twister(MTと略す)とは,松本眞(現広島大教授)・西村拓士(現山形大助手)氏による乱数発生アルゴリズムです。
詳細はこちらをどうぞ→Mersenne Twister Home Page

少なくとも,ABで使っている乱数よりは性能がいいそうです。
また,不必要に乗算を使用しないためかABの乱数よりも少しだけ速いです。

<処理時間比較>
MT法 ABのrand (時間:msec)
100,000回 71 76
500,000回 371 394
1,000,000回 759 878
5,000,000回 3855 4276

大学の授業に出てきたのですが,興味があったので原版のCプログラムをABに移植してみました。
※関数リファレンスは省略しました。使用法ドキュメントをお読み下さい。

[hide=コードを表示する]関数定義部
コード:
' ############################################################
' ##                                                        ##
' ##  メルセンヌツイスターによる乱数発生ライブラリ AB版     ##
' ##  Sinryow http://www.sinryow.net/                       ##
' ##                                                        ##
' ############################################################

' 【注記】
' 「/* ~ */」 の形のコメントの部分は,原プログラム
' http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c
' に書かれていたコメントそのものです。私が追加したコメントは全
' て「' ~」の形になっています。
'
' また英語の訳の正確性についての保証は致しかねます。もし間違い
' がございましたら,私のホームページ,あるいは私のメールアドレ
' ス mail (at-mark) sinryow.netにご連絡下さい。

/* 
   A C-program for MT19937, with initialization improved 2002/1/26.
   Coded by Takuji Nishimura and Makoto Matsumoto.

   Before using, initialize the state by using init_genrand(seed)  
   or init_by_array(init_key, key_length).

   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
   All rights reserved.						  

   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions
   are met:

     1. Redistributions of source code must retain the above copyright
        notice, this list of conditions and the following disclaimer.

     2. Redistributions in binary form must reproduce the above copyright
        notice, this list of conditions and the following disclaimer in the
        documentation and/or other materials provided with the distribution.

     3. The names of its contributors may not be used to endorse or promote 
        products derived from this software without specific prior written 
        permission.

   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


   Any feedback is very welcome.
   http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
   email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
*/

' 【訳】
' MT19937(メルセンヌツイスター)のCのプログラムです。初期化の
' 方法を2002年1月26日に改良しました。
'
' 使用する前に,init_genrand(seed) あるいは init_by_array(init_key, key_length).
' を用いて乱数の状態を初期化して下さい。
' 
'   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
'   All rights reserved.						  
'
' 以下の条件にあてはまる限り,ソースコードでもバイナリ形式でも,
' 改変があってもなくても再配布・利用が可能です。(訳注1)
' 
' <利用条件>
' 1. ソースコードの再配布においては,上記の著作権表示,この利
'	用条件および下記の免責条項を残さなければならない。
' 2. バイナリ形式での再配布においては,それと共に配布されるド
'	キュメント等に,上記の著作権表示,この利用条件および下記
'	の免責条項を複製して載せなければならない。
' 3. このプログラムの作成に関わった者は,このプログラムを利用
'	して作られた別のプログラムについては,事前に明記された許
'	可がない限り,そのプログラムを発展させた者の名前として用
'	いてはならない。
'
' <免責条項>
' 本ソフトウェアは著作権者によって,「現に提供されている状態の
' ままで」提供されるものとする。本ソフトウェアについては、明示
' あるいは黙示を問わず,商用品として通常そなえるべき品質をそな
' えているとの保証も,特定の目的に適合するとの保証を含め,何の
' 保証もないものとする。事由のいかんを問わず,損害発生の原因い
' かんを問わず,且つ,責任の根拠が契約であるか厳格責任であるか
' (過失その他の)不法行為であるかを問わず,著作権者も寄与者も,
' 仮にそのような損害が発生する可能性を知らされていたとしても,
' 本ソフトウェアの使用から発生した直接的な・間接的な・偶発的な・
' 特別な・懲罰的な・結果的な損害のいずれに対しても(代替品また
' はサービスの提供,使用機会,データまたは利益の損失の補償,ま
' たは,業務の中断に対する補償を含め)責任を一切負わないものと
' する。
'
' ご意見・ご感想はこちらへお願いします。(訳注2)
' http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
' メール: m-mat @ math.sci.hiroshima-u.ac.jp
'		 (スパム防止のためスペースを入れてあります)
'
' 【訳注1】
' 日本語訳の部分は法的には有効ではありませんので,実際にこのプ
' ログラムの再配布を行う場合は,「英語の部分を」再配布するもの
' に付加する必要があります。
'
' 【訳注2】
' AB版のプログラムについてのお問い合わせは,私(Sinryow)にお
' 願いします。松本氏(C言語版作者)には問い合わせないようお願
' いします。

Const __MT_N = 624
Const __MT_M = 397
Const __MT_MATRIX_A = &H9908b0df   /* constant vector a */
Const __MT_UPPER_MASK = &H80000000 /* most significant w-r bits */
Const __MT_LOWER_MASK = &H7fffffff /* least significant r bits */

Dim __mt_mt[ELM(__MT_N)] As DWord /* the array for the state vector  */
Dim __mt_mti = (__MT_N)+1 As Long /* mti==N+1 means mt[N] is not initialized */

/* initializes mt[N] with a seed */
Sub init_genrand(s As DWord)
	__mt_mt[0] = s
	'__mt_mt[0] = s And &Hffffffff
	
	For __mt_mti = 1 To (__MT_N)-1
		__mt_mt[(__mt_mti)] = _
			(1812433253 As DWord * (__mt_mt[(__mt_mti)-1] Xor (__mt_mt[(__mt_mti)-1] >> 30)) _
			+ (__mt_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_mt[__mt_mti] = __mt_mt[__mt_mti] And &Hffffffff
		/* for >32 bit machines */
	Next
End Sub

/* 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 */
Sub init_by_array(init_key As DWordPtr, key_length As Long)
	Dim i As Long, j As Long, k As Long
	init_genrand(19650218)
	i=1: j=0
	If (__MT_N)>key_length Then
		k = (__MT_N)
	Else
		k = key_length
	End If
	
	While k
		__mt_mt[i] = _
			(__mt_mt[i] Xor ((__mt_mt[i-1] Xor (__mt_mt[i-1] >> 30)) * 1664525 As DWord)) _
			+ init_key[j] + j /* non linear */
		'__mt_mt[i] = __mt_mt[i] And &Hffffffff; /* for WORDSIZE > 32 machines */
		i = i+1: j = j+1
		If i>=(__MT_N) Then
			__mt_mt[0] = __mt_mt[(__MT_N)-1]
			i=1
		End If
		If j>=key_length Then j=0

		k = k-1
	Wend
	
	For k=(__MT_N)-1 To 1 Step -1
		__mt_mt[i] = _
			(__mt_mt[i] Xor ((__mt_mt[i-1] Xor (__mt_mt[i-1] >> 30)) * 1566083941 As DWord)) _
			- i /* non linear */
		'__mt_mt[i] = __mt_mt[i] And &Hffffffff /* for WORDSIZE > 32 machines */
		i = i+1
		If i>=(__MT_N) Then
			__mt_mt[0] = __mt_mt[(__MT_N)-1]
			i=1
		End If
	Next
	
	__mt_mt[0] = &H80000000; /* MSB is 1; assuring non-zero initial array */ 
End Sub

/* generates a random number on [0,0xffffffff]-interval */

Dim __mt_mag01[ELM(2)]=[0, __MT_MATRIX_A] As DWord
/* mag01[x] = x * MATRIX_A  for x=0,1 */

Function genrand_int32() As DWord
	Dim y As DWord

	If __mt_mti >= (__MT_N) Then
		/* generate N words at one time */
		Dim kk As Long

		If __mt_mti = (__MT_N)+1 Then init_genrand(5489)
			/* if init_genrand() has not been called, */
			/* a default initial seed is used */

		For kk=0 To (__MT_N)-(__MT_M)-1
			y = (__mt_mt[kk] And __MT_UPPER_MASK) Or (__mt_mt[kk+1] And __MT_LOWER_MASK)
			__mt_mt[kk] = __mt_mt[kk+(__MT_M)] Xor (y >> 1) Xor __mt_mag01[y And &H1]
		Next
		
		While kk<(__MT_N)-1
			y = (__mt_mt[kk] And __MT_UPPER_MASK) Or (__mt_mt[kk+1] And __MT_LOWER_MASK)
			__mt_mt[kk] = __mt_mt[kk+(__MT_M)-(__MT_N)] Xor (y >> 1) Xor __mt_mag01[y And &H1]
			kk = kk+1
		Wend
		y = (__mt_mt[(__MT_N)-1] And __MT_UPPER_MASK) Or (__mt_mt[0] And __MT_LOWER_MASK)
		__mt_mt[(__MT_N)-1] = __mt_mt[(__MT_M)-1] Xor (y >> 1) Xor __mt_mag01[y And &H1]

		__mt_mti = 0
	End If

	y = __mt_mt[__mt_mti]
	__mt_mti = __mt_mti+1
	
	/* Tempering */
	y = y Xor (y >> 11)
	y = y Xor ((y << 7) And &H9d2c5680 As DWord)
	y = y Xor ((y << 15) And &Hefc60000 As DWord)
	y = y Xor ((y >> 18))

	genrand_int32 = y
End Function

/* generates a random number on [0,0x7fffffff]-interval */
Function genrand_int31() As Long
	genrand_int31 = genrand_int32()>>1
End Function

/* generates a random number on [0,1]-real-interval */
Function genrand_real1() As Double
	genrand_real1 = (genrand_int32() As Double)*(1.0 As Double/4294967295.0 As Double)
	/* divided by 2^32-1 */ 
End Function

/* generates a random number on [0,1)-real-interval */
Function genrand_real2() As Double
	genrand_real2 = (genrand_int32() As Double)*(1.0 As Double/4294967296.0 As Double)
	/* divided by 2^32 */
End Function

/* generates a random number on (0,1)-real-interval */
Function genrand_real3() As Double
	genrand_real3 = ((genrand_int32() As Double) + 0.5)*(1.0 As Double/4294967296.0 As Double)
	/* divided by 2^32 */
End Function

/* generates a random number on [0,1) with 53-bit resolution*/
Function genrand_res53() As Double
	Dim a = genrand_int32()>>5 As DWord, b = genrand_int32()>>6 As DWord
	genrand_res53 =
		((a As Double)*(67108864.0 As Double)+(b As Double))* _
		(1.0 As Double/9007199254740992.0 As Double)
End Function

/* These real versions are due to Isaku Wada, 2002/01/09 added */
使用例
コード:
#include "Mersenne.idx"

#console
Declare Function timeGetTime Lib "winmm.dll" () As Long
' ↓ ここからプログラムが実行されます

Dim i As Long
Dim init[ELM(4)]=[&H123, &H234, &H345, &H456] As DWord
Dim length=4 As Long

init_by_array(init, length)

' 以下,乱数発生テスト
/*
Print "1000 outputs of genrand_int32()"
For i=0 To 999
	Print genrand_int32()
Next
Print
Print "1000 outputs of genrand_real2()"
For i=0 To 999
	Print genrand_real2()
Next
*/

' 以下,速度測定用
Dim StartTime As Long
Dim EndTime As Long
Dim REPEATTIME As Long
Input "Repeat"; REPEATTIME

StartTime=timeGetTime()

For i=1 To REPEATTIME
	genrand_int32()
Next

EndTime = timeGetTime()
Print "Mersenne", EndTime-StartTime

StartTime=timeGetTime()

For i=1 To REPEATTIME
	rand()
Next

EndTime = timeGetTime()
Print "rand()",, EndTime-StartTime


Dim s As String
Input s
[/hide]
※補足
AB4用ですが,「/* ~ */」型のコメントを「' ~ 」にして,Dim文中で変数の初期化を行わなければ,AB2・3でも使えるかもしれません。

_________________
' ============================================================
' Sinryow Game Home Page - http://www.sinryow.net/
' Sinryow ActiveBasic Center - http://ab.sinryow.net/
' ============================================================


通報する
ページトップ
 記事の件名:
投稿記事Posted: 2005年12月13日(火) 02:15 
> 少なくとも,ABで使っている乱数よりは性能がいい
> ABのrand

まず、比較が間違っています。
Abにおける正式仕様はrandではなくてRnd。
randは隠し関数なので、いつ、こっそり変えられても文句を言ってはならないのです。
だから、Rndで見ないとなんない。

>>>>
ヘルプより。
乱数は 0 以上、1 未満の実数で表されます。

乱数関連のABソース再掲。(AB 4.12.01)
Function Rnd() As Double
Rnd=CDbl(rand()*rand() Mod 100000)/100000
End Function
Dim _System_RndNext As DWord
Function rand() As DWord
_System_RndNext=_System_RndNext*1103515245+12345
rand=((_CUDbl(_System_RndNext)/65536) As QWord) Mod 32768
End Function
>>>>>

以上より、ABのRndはハチャメチャもいいところ。
まず、ヘルプによれば、どう見ても、
・0~1の一様乱数を返し、その精度は倍精度相当である。
としか読めません。でも、
>Rnd=CDbl(rand()*rand() Mod 100000)/100000
であるから、「10進5桁の一様乱数」なんです。 ただしこれは許容するとします。
(ヘルプからこういう仕様であると判断することは不可能です。でも許容するものとします。)
で、ホントに「10進5桁の一様乱数」なのか?とんでもない。
「10進5桁の一様乱数」イコール0~99999の整数の一様乱数 ですよね?
で、rand()は、0~32768の整数の一様乱数 だから、
(rand()*rand() Mod 100000) の1の位を考えればはっきりしています。
1の位が偶数である確率:75%。 奇数である確率:25%。
こんなに偏った数値を返すのにランダムに返すとは何のこっちゃ!そりゃないよ。
Function Rnd() As Double
Rnd=(rand()+rand()\32768)/32768
End Function
としたほうがどれだけマシなことか。
(単精度の範囲で一様乱数をキープ。でも、結果は単精度という注釈必要。)


あと、
Function rand() As DWord  ですが、こちらは、一応は、混合合同法という
古典的ながらそれなりに効力のあるアルゴリズムなのでコメント無し。
ただし、
Function rand() As DWord
_System_RndNext=_System_RndNext*1103515245+12345
rand=((_CUDbl(_System_RndNext)/65536) As QWord) Mod 32768
End Function
というのは何か変で、
Function rand() As Int64
_System_RndNext=_System_RndNext*1103515245+12345
_System_RndNext=_System_RndNext Mod 2^32
rand=((_CUDbl(_System_RndNext)/65536) As QWord) Mod 32768
End Function
と組まない限り、桁オーバーでアウト、と見るのが普通なんだろうけど、何で桁オーバーせず動くの?
ヘルプ見ても、何故桁オーバーしないか、まるでわからん。


通報する
ページトップ
   
 記事の件名:
投稿記事Posted: 2005年12月13日(火) 08:12 
オフライン

登録日時: 2005年5月31日(火) 09:34
記事: 141
住所: 北海道
引用:
Mersenne Twisterについては、以前このフォーラムでSinryow氏によりAB版に移植されたものが公開されていましたが(http://www.discoversoft.net/forum/viewt ... 5%F3%A5%CC)、
その使用法について質問させていただきます。

MTの乱数を初期化する際、32ビットより長い初期シードを用いたい場合は、init_by_array()により任意長の配列を用いて初期化することができるとされていますが、具体的にはどのように記述するのでしょうか?

自分としては、128ビットの値を32ビット×4の配列にし、それを初期シードとしてMTに渡したいと思っているのですが、init_by_array()の引数をどのように記述すれば良いのかが良くわからないのです。

(ABには直接関係ない内容かとは思いましたが)よろしくお願いいたします。
以下のようになるかと思います。試していませんが。
コード:
Dim Seed[ELM(4)] As DWord
Seed[0] = &H83024923
Seed[1] = &H14325091
Seed[2] = &H01328943
Seed[3] = &H65129462
init_by_array(Seed, 4) ' 「4」は配列の要素数

_________________
' ============================================================
' Sinryow Game Home Page - http://www.sinryow.net/
' Sinryow ActiveBasic Center - http://ab.sinryow.net/
' ============================================================


通報する
ページトップ
 記事の件名: MT厨
投稿記事Posted: 2010年6月06日(日) 20:16 
MTはええねMT


通報する
ページトップ
   
期間内表示:  ソート  
新しいトピックを投稿する  トピックへ返信する  [ 4 件の記事 ] 

All times are UTC+09:00


オンラインデータ

このフォーラムを閲覧中のユーザー: なし & ゲスト[1人]


トピック投稿:  可
返信投稿:  可
記事編集: 不可
記事削除: 不可
ファイル添付: 不可

検索:
ページ移動:  
cron
Powered by phpBB® Forum Software © phpBB Limited
Japanese translation principally by KONISHI Yohsuke