MTの実装

偶然ですが、Movable Type(略してMT)を導入したいという衝動に駆られる前は、Mersenne Twister(略してMT)のC#による実装をしていました。ちなみにMersenne Twisterというのは、松本眞さんと西村拓士さんによって開発された非常に高性能な擬似乱数ジェネレータのことです。

一通り実装し終えたので、機能的にオリジナルと等価であることを確かめた後、32ビット符号無し整数で乱数を与えるgenrand_int32()関数を使い、Pentium4 1.70GHzマシンでパフォーマンステストをしてみました。結果がこちら:

• オリジナルC言語版 genrand_int32() : 1.718746e-008 seconds per call
• 今回のC#版 genrand_int32() : 2.440930e-008 seconds per call

大体3/4程度のパフォーマンスですが、まぁマネージドコードにしては頑張ってる方でしょう。

次に、unsigned intを返すgenrand_int32()を4294967296.0で割って、[ 0.0, 1.0 )に均等に分布する実数乱数を得る関数、genrand_real2()でも試してみました。

• オリジナルC言語版 genrand_real2() : 3.309488e-008 seconds per call
• 今回のC#版 genrand_real2() : 7.074817e-008 seconds per call

......何、これ。doubleにキャストして4294967296.0で割っただけなのに、差が付きすぎじゃあありませんか? ちなみに、C言語版では最適化によってgenrand_real2()のインライン展開が行われていましたが、それをしなくても3.590222e-008 seconds per callの速さです。

こんな時は逆アセンブルに限ります。ildasmで覗いてみると、この関数の中では

  IL_0000:  ldarg.0
  IL_0001:  call       instance uint32 MersenneTwister::genrand_int32()
  IL_0006:  conv.r.un
  IL_0007:  conv.r8
  IL_0008:  ldc.r8     2.3283064365386963e-010
  IL_0011:  mul
  IL_0012:  ret

という処理が行われていることが判りました。どうも.NET Framework 2.0のSDKにはMSIL(CIL)のリファレンスがついていないようで、それのためだけにバージョン1.1のSDKを引っ張ってくるのも面倒なのでECMA Specs - Part 3 (CIL)を見てこれにコメントを振ると、

  IL_0000:  ldarg.0  // this参照をスタックにプッシュ
  IL_0001:  call       instance uint32 MersenneTwister::genrand_int32()
  // this.genrand_int32()の呼び出し
  IL_0006:  conv.r.un  // スタックの先頭の
  // genrand_int32()の戻り値を、uint32からfloatに変換
  IL_0007:  conv.r8  // スタックの先頭のfloat値をfloat64値に変換
  IL_0008:  ldc.r8     2.3283064365386963e-010  // 定数をスタックにプッシュ
  IL_0011:  mul  // スタックの先頭の二つを掛け合わせて
  IL_0012:  ret  // 制御を返す

ということのようです。ぱっと見、conv.r.unとconv.r8はだぶっているようにも見えるのですが。ちなみに、CLR 2.0からの新機能のような気がしますが、Visual Studio 2005ではこのコードの機械語版も見られます(素晴らしい!)。尤も、C++のプログラムを逆アセンブルするほど読みやすくはないのですが。で、そのアセンブリ語訳が:

00000015  mov         ecx,edi 
00000017  call        dword ptr ds:[00A78CBCh] 
0000001d  mov         esi,eax 
0000001f  xor         eax,eax 
00000021  mov         dword ptr [esp],esi 
00000024  mov         dword ptr [esp+4],eax 
00000028  fild        qword ptr [esp] 
0000002b  fstp        qword ptr [esp] 
0000002e  fld         qword ptr [esp] 
00000031  fmul        dword ptr ds:[00DB1934h] 
00000037  add         esp,8 
0000003a  pop         esi  
0000003b  pop         edi  
0000003c  ret

アセンブリ言語レベルの最適化に詳しいわけではないですが、0x2b～0x2eが無駄な手順を踏んでいることくらいすぐわかります。これはfst qword ptr [esp]で一発です。その前に、わざわざespの指すアドレスに値をストアしたところで、使われていません。それに、0x1d～0x28はmov dword ptr [esp], eax、fild dword ptr [esp]とやれば済む話です。

これは困ったものです。C#のコンパイラには任せておけないとMSILを自分で書き換えるのでは、折角生産性を高めるためにC#を使っている意味がありません。

...が、これには抜け道が存在しました。uint32を直接float64に変換するのは無理なようですが、一旦int32を介せばスムーズに変換できるようなのです。即ち、今までは

public double genrand_real2() {
  return this.genrand_int32() * ( 1.0 / 4294967296.0 );
}

としていたところを

public double genrand_real2() {
  return (int)( this.genrand_int32() >> 1 ) * ( 1.0 / 2147483648.0 );
}

で書き換えるのです。精度は1ビット落ちるものの、これで速度が手に入れば安いものです。これを再びコンパイルし、逆アセンブルしたところ、MSILの段階でconv.r.unが消え、機械語レベルでも

00000013  mov         ecx,edi 
00000015  call        dword ptr ds:[00A78CBCh] 
0000001b  mov         esi,eax 
0000001d  shr         esi,1 
0000001f  mov         dword ptr [esp],esi 
00000022  fild        dword ptr [esp] 
00000025  fmul        dword ptr ds:[00DB1924h] 
0000002b  pop         ecx  
0000002c  pop         esi  
0000002d  pop         edi  
0000002e  ret

と、望ましいコードが得られました。勿論、そのパフォーマンスは

• 修正C#版 genrand_real2() : 3.681116e-008 seconds per call

と、オリジナルに迫る速度を弾き出してくれました。まぁ、Mersenne TwisterにはMMXやSSEを使って更に最適化されたバージョンもあるのですが、それらのことは考えないようにしましょう。これだけ出れば個人的には満足です。

今回の教訓 : 符号無し整数から浮動小数点数への直接の変換は行わないこと。