ab.com コミュニティ https://www.activebasic.com/forum/ |
|
[AB4] 乱数発生アルゴリズム「Mersenne Twister」 https://www.activebasic.com/forum/viewtopic.php?t=517 |
ページ 1 / 1 |
作成者: | Sinryow [ 2005年12月02日(金) 03:31 ] |
記事の件名: | [AB4] 乱数発生アルゴリズム「Mersenne Twister」 |
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に移植してみました。 ※関数リファレンスは省略しました。使用法ドキュメントをお読み下さい。 コードを表示する [ここをクリックすると内容が表示されます] 関数定義部 コード: ' ############################################################ ' ## ## ' ## メルセンヌツイスターによる乱数発生ライブラリ AB版 ## ' ## Sinryow http://www.sinryow.net/ ## ' ## ## ' ############################################################ ' 【注記】 ' 「/* ~ */」 の形のコメントの部分は,原プログラム ' http://www.math.sci.hiroshima-u.ac.jp/~ ... t19937ar.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/~ ... T/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/~ ... T/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 = _ (__mt_mt Xor ((__mt_mt[i-1] Xor (__mt_mt[i-1] >> 30)) * 1664525 As DWord)) _ + init_key[j] + j /* non linear */ '__mt_mt = __mt_mt 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 = _ (__mt_mt Xor ((__mt_mt[i-1] Xor (__mt_mt[i-1] >> 30)) * 1566083941 As DWord)) _ - i /* non linear */ '__mt_mt = __mt_mt 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 ※補足 AB4用ですが,「/* ~ */」型のコメントを「' ~ 」にして,Dim文中で変数の初期化を行わなければ,AB2・3でも使えるかもしれません。 |
作成者: | 河川屋 [ 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 と組まない限り、桁オーバーでアウト、と見るのが普通なんだろうけど、何で桁オーバーせず動くの? ヘルプ見ても、何故桁オーバーしないか、まるでわからん。 |
作成者: | Sinryow [ 2005年12月13日(火) 08:12 ] |
記事の件名: | |
引用: 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」は配列の要素数 |
作成者: | MT [ 2010年6月06日(日) 20:16 ] |
記事の件名: | MT厨 |
MTはええねMT |
ページ 1 / 1 | 全ての表示時間は UTC+09:00 です |
Powered by phpBB® Forum Software © phpBB Limited https://www.phpbb.com/ |