正逆双方向に計算できる SFMT
正逆双方向に計算できる SFMT のクラスを作ってみた（※ただし周期が \(2^{19937}-1\) のものに限定した）。
逆算部分の理論はメルセンヌ・ツイスタの tempering の逆関数や TinyMT の更新関数の逆関数とほぼ変わらないので、詳細はそちらを参照。
速度を度外視して汎用性重視で作ったので高速化は各自で。
というか高速化プログラミング苦手だし、そもそも C# 使ってる時点で遅いし、みたいな。
それとおそらく SIMD-oriented な実装になってないので、その意味でも Fast じゃない。
結果が SFMT と同じになるだけ……
UInt128.cs
SFMT の更新の様子が分かりやすくなるように、 128ビットの符号なし整数型を使用する。
using System; using System.Collections.Generic; using System.Text; // C# namespace SFMTTest { internal struct UInt128 : IEquatable<UInt128> { #region constructor public UInt128(uint word3,uint word2,uint word1,uint word0) { this._w3 = word3; this._w2 = word2; this._w1 = word1; this._w0 = word0; } #endregion #region static field private const int WORD_SIZE = 32; #endregion #region static property public static UInt128 Zero { get { return new UInt128(0u,0u,0u,0u); } } public static UInt128 One { get { return new UInt128(0u,0u,0u,1u); } } #endregion #region field property // 第３ワード（最上位ワード） private readonly uint _w3; public uint Word3 { get { return _w3; } } // 第２ワード private readonly uint _w2; public uint Word2 { get { return _w2; } } // 第１ワード private readonly uint _w1; public uint Word1 { get { return _w1; } } // 第０ワード（最下位ワード） private readonly uint _w0; public uint Word0 { get { return _w0; } } #endregion #region property // 上位倍長ワード（上位64bit） public ulong DWordH { get { return ( (ulong)_w3 << 32 ) | _w2; } } // 下位倍長ワード（下位64bit） public ulong DWordL { get { return ( (ulong)_w1 << 32 ) | _w0; } } #endregion #region public bool IsZero { get { return _w0 == 0x0u && _w1 == 0x0u && _w2 == 0x0u && _w3 == 0x0u; } } #endregion #region indexer public uint this[int n] { get { switch( n & 3 ) { case 0: return _w0; case 1: return _w1; case 2: return _w2; case 3: return _w3; default: throw new IndexOutOfRangeException(); } } } #endregion #region property public string TextBin { get { return Convert.ToString(_w3,2).PadLeft(WORD_SIZE,'0') + Convert.ToString(_w2,2).PadLeft(WORD_SIZE,'0') + Convert.ToString(_w1,2).PadLeft(WORD_SIZE,'0') + Convert.ToString(_w0,2).PadLeft(WORD_SIZE,'0'); } } public string TextHex { get { return _w3.ToString("X08") + _w2.ToString("X08") + _w1.ToString("X08") + _w0.ToString("X08"); } } #endregion #region operator public static explicit operator UInt128(uint x) { return new UInt128(0u,0u,0u,x); } public static explicit operator UInt128(ulong x) { return new UInt128(0u,0u,(uint)( x>>32 ),(uint)x); } public static UInt128 operator&(UInt128 x,UInt128 y) { return new UInt128( x._w3 & y._w3, x._w2 & y._w2, x._w1 & y._w1, x._w0 & y._w0 ); } public static UInt128 operator^(UInt128 x,UInt128 y) { return new UInt128( x._w3 ^ y._w3, x._w2 ^ y._w2, x._w1 ^ y._w1, x._w0 ^ y._w0 ); } //public static UInt128 operator|(UInt128 x,UInt128 y) { // return new UInt128( // x._w3 | y._w3, // x._w2 | y._w2, // x._w1 | y._w1, // x._w0 | y._w0 // ); //} public static UInt128 operator<<(UInt128 x,int n) { n &= 0x7f; uint w3 = x._w3; uint w2 = x._w2; uint w1 = x._w1; uint w0 = x._w0; for( ;n>=WORD_SIZE;n-=WORD_SIZE ) { w3 = w2; w2 = w1; w1 = w0; w0 = 0x0u; } if( n > 0 ) { w3 = w3 << n | w2 >> ( WORD_SIZE-n ); w2 = w2 << n | w1 >> ( WORD_SIZE-n ); w1 = w1 << n | w0 >> ( WORD_SIZE-n ); w0 = w0 << n; } return new UInt128(w3,w2,w1,w0); } public static UInt128 operator>>(UInt128 x,int n) { n &= 0x7f; uint w3 = x._w3; uint w2 = x._w2; uint w1 = x._w1; uint w0 = x._w0; for( ;n>=WORD_SIZE;n-=WORD_SIZE ) { w0 = w1; w1 = w2; w2 = w3; w3 = 0x0u; } if( n > 0 ) { w0 = w0 >> n | w1 << ( WORD_SIZE-n ); w1 = w1 >> n | w2 << ( WORD_SIZE-n ); w2 = w2 >> n | w3 << ( WORD_SIZE-n ); w3 = w3 >> n; } return new UInt128(w3,w2,w1,w0); } #endregion #region override method public override string ToString() { return TextHex; } #endregion #region implement method public bool Equals(UInt128 other) { return this.Word0 == other.Word0 && this.Word1 == other.Word1 && this.Word2 == other.Word2 && this.Word3 == other.Word3; } #endregion } }
InvertibleSFMT19937.cs
分かりやすさを重視して上記の UInt128 構造体を使っている。ので遅い。
実際に使うなら uint[4] で実装したほうがいいと思う。
ポケモンサンムーンの乱数計算にしか使わないのなら、（32ビットの乱数生成をしないので） ulong[2] のほうがなおいい。タブンネ。
using System; using System.Collections.Generic; using System.Text; // C# namespace SFMTTest { public class InvertibleSFMT19937 { #region static constructor static InvertibleSFMT19937() { const uint MSK1 = 0xdfffffefu; // SFMT_MSK1 const uint MSK2 = 0xddfecb7fu; // SFMT_MSK2 const uint MSK3 = 0xbffaffffu; // SFMT_MSK3 const uint MSK4 = 0xbffffff6u; // SFMT_MSK4 const uint PARITY1 = 0x00000001u; // SFMT_PARITY1 const uint PARITY2 = 0x00000000u; // SFMT_PARITY2 const uint PARITY3 = 0x00000000u; // SFMT_PARITY3 const uint PARITY4 = 0x13c9e684u; // SFMT_PARITY4 const uint sr1Mask = 0xffffffffu >> SR1_BIT; const uint sl1Mask = 0xffffffffu << SL1_BIT; // initialize MSK_B MSK_B = new UInt128( // MSK_B = MSK4 & sr1Mask, // ( SFMT_MSK4 & ( 0xffffffffu >> SFMT_SR1 ) ) << 96 | MSK3 & sr1Mask, // ( SFMT_MSK3 & ( 0xffffffffu >> SFMT_SR1 ) ) << 64 | MSK2 & sr1Mask, // ( SFMT_MSK2 & ( 0xffffffffu >> SFMT_SR1 ) ) << 32 | MSK1 & sr1Mask // ( SFMT_MSK1 & ( 0xffffffffu >> SFMT_SR1 ) ) ); // initialize MSK_D MSK_D = new UInt128( // MSK_D = sl1Mask, // ( 0xffffffffu << SFMT_SL1 ) << 96 | sl1Mask, // ( 0xffffffffu << SFMT_SL1 ) << 64 | sl1Mask, // ( 0xffffffffu << SFMT_SL1 ) << 32 | sl1Mask // ( 0xffffffffu << SFMT_SL1 ) ); // initialize PARITY PARITY = new UInt128( // PARITY = PARITY4, // SFMT_PARITY4 << 96 | PARITY3, // SFMT_PARITY3 << 64 | PARITY2, // SFMT_PARITY2 << 32 | PARITY1 // SFMT_PARITY1 ); } #endregion #region constructor public InvertibleSFMT19937(uint seed) { InitGenRand(seed); } #endregion #region static field private const int MEXP = 19937; // SFMT_MEXP private const int N = MEXP / 128 + 1; // SFMT_N //private const int N32 = N * 4; // SFMT_N32 //private const int N64 = N * 2; // SFMT_N64 private const int POS1 = 122; // SFMT_POS1 private const int SL1_BIT = 18; // SFMT_SL1 // shift-left bit size private const int SL2_BYTE = 1; // SFMT_SL2 // shift-left byte size private const int SR1_BIT = 11; // SFMT_SR1 // shift-right bit size private const int SR2_BYTE = 1; // SFMT_SR2 // shift-right byte size private const int SL2_BIT = SL2_BYTE * 8; // // shift-left bit size private const int SR2_BIT = SR2_BYTE * 8; // // shift-right bit size private static readonly UInt128 MSK_B; // // initialize in static constructor private static readonly UInt128 MSK_D; // // initialize in static constructor private static readonly UInt128 PARITY; // // initialize in static constructor #endregion #region field private UInt128[] _state = new UInt128[N]; // 元のソースの sfmt_t.state に相当 private int _index = 0; // 元のソースの idx に相当 ※ idx = ( _index % N ) * 4 + _pos private int _pos = 0; // 同上 #endregion #region static method /// <summary>x を N で割った剰余（余り）を求める。</summary> /// <param name="x">剰余を求める値</param> /// <returns>x を N で割った剰余。 0 以上 N 未満の値を返す。</returns> private static int ModN(int x) { return x>=0 ? x % N : N - 1 - ( -1-x ) % N; } // 元のソースの do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c, w128_t *d) に相当 private static UInt128 DoRecursion(UInt128 a,UInt128 b,UInt128 c,UInt128 d) { UInt128 x = a << SL2_BIT; UInt128 y = c >> SR2_BIT; UInt128 r = a ^ x ^ ( ( b >> SR1_BIT ) & MSK_B ) ^ y ^ ( ( d << SL1_BIT ) & MSK_D ); return r; } // DoRecursion() の逆関数 private static UInt128 DoRecursionInv(UInt128 r,UInt128 b,UInt128 c,UInt128 d) { UInt128 y = c >> SR2_BIT; UInt128 a = r ^ ( ( b >> SR1_BIT ) & MSK_B ) ^ y ^ ( ( d << SL1_BIT ) & MSK_D ); for( int sh=SL2_BIT;sh<128;sh<<=1 ) { a ^= ( a << sh ); } return a; } #endregion #region method / init // 元のソースの sfmt_init_gen_rand(sfmt_t * sfmt, uint32_t seed) に相当 public void InitGenRand(uint seed) { uint[] arr = new uint[5]; arr[0] = seed; for( uint j=0;j<N;j++ ) { for( uint i=1;i<=4;i++ ) { arr[i] = 1812433253u * ( arr[i-1] ^ ( arr[i-1]>>30 ) ) + i + j * 4; } _state[j] = new UInt128(arr[3],arr[2],arr[1],arr[0]); arr[0] = arr[4]; } PeriodCertification(); _index = 0; _pos = 0; // テーブル全体を更新 for( int i=0;i<N;i++ ) { UpdateNext(); } return; } // 元のソースの period_certification(sfmt_t * sfmt) に相当 private void PeriodCertification() { // 立っているビット数の偶奇を数える UInt128 inner = _state[0] & PARITY; for( int i=64;i>0;i>>=1 ) { inner ^= inner >> i; } // 奇数個ならＯＫ if( ( inner.Word0 & 1u ) == 1u ) { return; } // 偶数個だとＮＧ → 修正 UInt128 work = UInt128.One; for( int i=0;i<128;i++ ) { if( !( work & PARITY ).IsZero ) { _state[0] ^= work; return; } work <<= 1; } return; } #endregion #region method / update // 元のソースの sfmt_gen_rand_all(sfmt_t * sfmt) 中の for ループの中身に相当 private void UpdateNext() { _state[ModN(_index)] = DoRecursion( _state[ModN(_index)], _state[ModN(_index+POS1)], _state[ModN(_index+N-2)], _state[ModN(_index+N-1)] ); _index++; _pos = 0; return; } // UpdateNext() の逆関数 private void UpdatePrev() { _index--; _pos = 3; _state[ModN(_index)] = DoRecursionInv( _state[ModN(_index)], _state[ModN(_index+POS1)], _state[ModN(_index+N-2)], _state[ModN(_index+N-1)] ); return; } #endregion #region method / get rand // 元のソースの sfmt_genrand_uint32(sfmt_t * sfmt) に相当 public uint NextUInt32() { uint r = _state[ModN(_index)][_pos++]; if( _pos >= 4 ) { UpdateNext(); } return r; } // 元のソースの sfmt_genrand_uint64(sfmt_t * sfmt) に相当 public ulong NextUInt64() { ulong next; switch( _pos & 3 ) { case 0: next = _state[ModN(_index)].DWordL; _pos = 2; return next; case 2: next = _state[ModN(_index)].DWordH; UpdateNext(); return next; default: // 元のソースの「 assert(sfmt->idx % 2 == 0); 」の部分に相当 throw new Exception("_posが奇数です"); } } // NextUInt32() の逆関数 public uint PrevUInt32() { _pos--; if( _pos < 0 ) { UpdatePrev(); } uint r = _state[ModN(_index)][_pos]; return r; } // NextUInt64() の逆関数 public ulong PrevUInt64() { ulong prev; switch( _pos & 3 ) { case 0: UpdatePrev(); _pos = 2; prev = _state[ModN(_index)].DWordH; return prev; case 2: _pos = 0; prev = _state[ModN(_index)].DWordL; return prev; default: // 元のソースの「 assert(sfmt->idx % 2 == 0); 」の部分に相当 throw new Exception("_posが奇数です"); } } #endregion } }
