正逆双方向に計算できるメルセンヌ・ツイスタ

正逆双方向に計算できるメルセンヌ・ツイスタのクラスを作ってみた。
理論などの詳細は下記参照。

• メルセンヌ・ツイスタのtemperingの逆関数に関する考察 - Plus Le Blog
• メルセンヌ・ツイスタの逆算に関する考察 - Plus Le Blog

速度を度外視して汎用性重視で作ったので高速化は各自で。
というか高速化プログラミング苦手だし、そもそもC#使ってる時点で遅いし、みたいな。

using System;
using System.Collections.Generic;
using System.Text;

// C#
namespace Tmp {
	public class InvertibleMT {
		#region constructor
		/// <summary>
		/// <para>seed から MT を作る。</para>
		/// <para>初期化方法はポケモンGen5準拠。</para>
		/// </summary>
		/// <param name="seed">seed</param>
		public InvertibleMT(uint seed) {
			this.table[0] = seed;
			for( int i=1;i<N;i++ ) {
				this.table[i] = ( ( this.table[i-1]>>30 ) ^ this.table[i-1] ) * 0x6c078965u + (uint)i;
			}
			for( int i=0;i<N;i++ ) {
				this.Next();
			}
			this.index = 0;
		}
		/// <summary>
		/// <para>【デバッグ用】</para>
		/// <para>テーブルを直接指定して MT を作る。</para>
		/// </summary>
		/// <param name="data">乱数の初期テーブル。多い場合は最初の N 個を使う。不足分は 0 で補う。</param>
		public InvertibleMT(IEnumerable<uint> data) {
			if( data != null ) {
				int i = 0;
				foreach( uint x in data ) {
					if( i >= N ) {
						break;
					}
					this.table[i++] = x;
				}
			}
		}
		#endregion

		#region const field
		private const int N = 624;
		private const int M = 397;
		private const uint MATRIX_A   = 0x9908b0dfu;
		private const uint UPPER_MASK = 0x80000000u;
		private const uint LOWER_MASK = 0x7fffffffu;
		private const uint TEMPERING_MASK_B  = 0x9d2c5680u;
		private const uint TEMPERING_MASK_B2 = 0x94284000u;		// = B & (B<<7)
		private const uint TEMPERING_MASK_B3 = 0x14200000u;		// = B & (B<<7) & (B<<14)
		private const uint TEMPERING_MASK_B4 = 0x10000000u;		// = B & (B<<7) & (B<<14) & (B<<21)
		private const uint TEMPERING_MASK_C  = 0xefc60000u;
		#endregion
		#region static field
		private static readonly uint[] mag01 = new uint[2] { 0x0u,MATRIX_A };
		#endregion
		#region field
		private readonly uint[] table = new uint[N];
		#endregion
		#region field property
		private int index = 0;
		/// <summary>現在のインデックスの位置</summary>
		public int Index { get { return index; } }
		#endregion

		#region static method
		/// <summary>１つ次の乱数値を計算する。</summary>
		/// <param name="data">要素数 N の配列</param>
		/// <param name="index">計算の起点の位置</param>
		/// <returns>１つ次の乱数値</returns>
		public static uint CalcNext(uint[] data,int index) {
			if( data == null ) {
				throw new ArgumentNullException();
			} else if( data.Length != N ) {
				throw new ArgumentException();
			}

			uint x = ( data[ModN(index)] & UPPER_MASK ) | ( data[ModN(index+1)] & LOWER_MASK );
			uint next = data[ModN(index+M)] ^ ( x >> 1 ) ^ mag01[x&0x1];
			return next;
		}
		/// <summary>１つ前の乱数値を計算する。</summary>
		/// <param name="data">要素数 N の配列</param>
		/// <param name="index">計算の起点の位置</param>
		/// <returns>１つ前の乱数値</returns>
		public static uint CalcPrev(uint[] data,int index) {
			if( data == null ) {
				throw new ArgumentNullException();
			} else if( data.Length != N ) {
				throw new ArgumentException();
			}

			uint tail = data[ModN(index+N-1)] ^ data[ModN(index+M-1)];
			uint body = tail ^ mag01[tail>>31];
			uint head = data[ModN(index+N)] ^ data[ModN(index+M)];
			uint prev
				= ( ( head <<  1 ) & UPPER_MASK )
				^ ( ( body <<  1 ) & LOWER_MASK )
				^   ( tail >> 31 );
			return prev;
		}
		/// <summary>１つ前の乱数値を逆算可能か否か調べる。</summary>
		/// <param name="data">要素数 N の配列</param>
		/// <returns>１つ前の乱数値を逆算可能な場合、true</returns>
		public static bool IsInvertible(uint[] data) {
			if( data == null || data.Length != N ) {
				return false;
			}

			uint tail = data[N-1] ^ data[M-1];
			uint body = tail ^ mag01[tail>>31];
			uint head = tail;
			return ( ( body ^ ( data[0] >>  1 ) ) << 2 ) == 0u
				&& ( ( head ^ ( data[0] << 31 ) ) & UPPER_MASK ) == 0u;
		}
		/// <summary>メルセンヌ・ツイスタの tempering 部分。</summary>
		/// <param name="x">tempering する値</param>
		/// <returns>tempering した値</returns>
		public static uint Tempering(uint x) {
			x ^= ( x >> 11 );
			x ^= ( x <<  7 ) & TEMPERING_MASK_B;
			x ^= ( x << 15 ) & TEMPERING_MASK_C;
			x ^= ( x >> 18 );
			return x;
		}
		/// <summary>メルセンヌ・ツイスタの tempering 部分の逆関数。</summary>
		/// <param name="x">tempering された値</param>
		/// <returns>tempering する前の値</returns>
		public static uint TemperingInv(uint x) {
			x ^=    ( x >> 18 );
			x ^=    ( x << 15 ) & TEMPERING_MASK_C;
			x ^=  ( ( x <<  7 ) & TEMPERING_MASK_B )
				^ ( ( x << 14 ) & TEMPERING_MASK_B2 )
				^ ( ( x << 21 ) & TEMPERING_MASK_B3 )
				^ ( ( x << 28 ) & TEMPERING_MASK_B4 );
			x ^=    ( x >> 11 ) ^ ( x >> 22 );
			return x;
		}
		/// <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;
		}
		#endregion

		#region method
		/// <summary>次の乱数値を返す。</summary>
		/// <returns>次の乱数値</returns>
		public uint Next() {
			//// ※index がオーバーフローして負になった場合の処理
			////   でも例外投げるのも面倒なので意味のない値でもいいから返すことにする
			//if( index <= -N ) {
			//    throw new InvalidOperationException();
			//}

			uint result = table[ModN(index)];
			table[ModN(index)] = CalcNext(table,index);		// table を次の状態に更新
			index++;
			return Tempering(result);
		}
		/// <summary>前の乱数値を返す。</summary>
		/// <returns>前の乱数値</returns>
		public uint Prev() {
			//// ※MTによって生成された最初の N 個の値を MT[0], MT[1], ..., MT[N-1] とするとき
			////   MT[-N] 以前は不可逆なので値が求めらない（無意味な値が出てくる）
			////   でも例外投げるのも面倒なので意味のない値でもいいから返すことにする
			//if( index <= -N ) {
			//    throw new InvalidOperationException();
			//}

			index--;
			uint result = CalcPrev(table,index);
			table[ModN(index)] = result;		// table を前の状態に更新
			return Tempering(result);
		}
		#endregion

	}
}

----------------------------------------------------------------

以下は、上のを書く前の古い記事（一応残しておいた）。
速報的な位置づけのプログラムなので参考程度に。

メルセンヌ・ツイスタを逆算するプログラム

以下に示すのは、メルセンヌ・ツイスタによって生成された連続する \(N\left(=624\right)\) 個の値から、その直前の \(N\) 個の値を求めるプログラムである。
ちなみにでたらめな \(N\) 個の値を入力した場合の挙動は未定義。タブンネ。

using System;
using System.Collections.Generic;
using System.Text;

// c#
namespace Tmp {
	class Program {
		private const int N = 624;
		private const int M = 397;
		private const uint MATRIX_A   = 0x9908b0dfu;
		private const uint UPPER_MASK = 0x80000000u;
		private const uint LOWER_MASK = 0x7fffffffu;
		private static uint[] mag01 = { 0x0,MATRIX_A };

		static void Main(string[] args) {
			// 準備
			Random rand = new Random();
			int skip = rand.Next(0,1000);
			uint seed = (uint)rand.Next();
			MT mt = new MT(seed);

			// 初期消費
			for( int i=0;i<skip;i++ ) {
				mt.Next();
			}

			// 答え合わせ用のテーブル
			uint[] answer = new uint[N];
			for( int i=0;i<N;i++ ) {
				answer[i] = mt.Next();
			}

			// 問題用のテーブル
			uint[] table = new uint[N];
			for( int i=0;i<N;i++ ) {
				uint x = mt.Next();
				table[i] = TemperingInv(x);		// tempering 前の状態に戻しておく
			}

			// １つ前のテーブルを求める（つまり table から answer を求める）
			for( int i=N-1;i>=0;i-- ) {
				uint head = ( ( table[( i+N )%N] ^ table[( i+M )%N] ^ mag01[table[( i+1 )%N]&0x1] ) << 1 ) & UPPER_MASK;
				uint tail = table[( i+N-1 )%N] ^ table[( i+M-1 )%N];
				uint body = ( ( tail ^ mag01[tail>>31] ) << 1 ) & LOWER_MASK;
				table[i] = head ^ body ^ ( tail>>31 );
			}

			// 答え合わせ
			bool isWrong = false;
			for( int i=0;i<N;i++ ) {
				uint x = Tempering(table[i]);	// 求めたテーブルは tempering 前の状態なので tempering する
				uint y = answer[i];
				if( x != y ) {
					isWrong = true;
				}
			}
			if( isWrong ) {
				Console.WriteLine("ちがうよ！");
			} else {
				Console.WriteLine("おめでと！");
			}
			Console.WriteLine();
			Console.WriteLine("skip : {0}",skip);
			Console.WriteLine("seed : 0x{0:x08}",seed);
			Console.WriteLine();
			Console.WriteLine("finish");
			Console.ReadLine();
			return;
		}

		// メルセンヌ・ツイスタの tempering 部分
		private static uint Tempering(uint x) {
			x ^= ( x >> 11 );
			x ^= ( x <<  7 ) & 0x9d2c5680u;
			x ^= ( x << 15 ) & 0xefc60000u;
			x ^= ( x >> 18 );
			return x;
		}
		// メルセンヌ・ツイスタの tempering の逆関数
		private static uint TemperingInv(uint x) {
			x ^=    ( x >> 18 );
			x ^=    ( x << 15 ) & 0xefc60000u;
			x ^=  ( ( x <<  7 ) & 0x9d2c5680u )
				^ ( ( x << 14 ) & 0x94284000u )
				^ ( ( x << 21 ) & 0x14200000u )
				^ ( ( x << 28 ) & 0x10000000u );
			x ^=    ( x >> 11 ) ^ ( x >> 22 );
			return x;
		}

	}

	public class MT {
		public MT(uint seed) {
			this.table[0] = seed;
			for( tableIndex=1;tableIndex<N;tableIndex++ ) {
				this.table[tableIndex] = ( ( this.table[tableIndex-1]>>30 ) ^ this.table[tableIndex-1] ) * 0x6C078965U + (uint)tableIndex;
			}
		}

		private const int N = 624;
		private const int M = 397;
		private const uint MATRIX_A   = 0x9908b0df;
		private const uint UPPER_MASK = 0x80000000;
		private const uint LOWER_MASK = 0x7fffffff;
		private const uint TEMPERING_MASK_B = 0x9d2c5680;
		private const uint TEMPERING_MASK_C = 0xefc60000;
		private static uint[] mag01 = { 0x0,MATRIX_A };
		private uint[] table = new uint[N];
		private short tableIndex;

		private static uint TEMPERING_SHIFT_U(uint y) { return ( y >> 11 ); }
		private static uint TEMPERING_SHIFT_S(uint y) { return ( y <<  7 ); }
		private static uint TEMPERING_SHIFT_T(uint y) { return ( y << 15 ); }
		private static uint TEMPERING_SHIFT_L(uint y) { return ( y >> 18 ); }

		public uint Next() {
			uint y;
			if( tableIndex >= N ) {
				// テーブルの更新
				short i = 0;
				for( ;i<N-M;i++ ) {
					y = ( table[i] & UPPER_MASK ) | ( table[i+1] & LOWER_MASK );
					table[i] = table[i+M] ^ ( y>>1 ) ^ mag01[y&0x1];
				}
				for( ;i<N-1;i++ ) {
					y = ( table[i] & UPPER_MASK ) | ( table[i+1] & LOWER_MASK );
					table[i] = table[i+( M-N )] ^ ( y>>1 ) ^ mag01[y&0x1];
				}
				{
					y = ( table[N-1] & UPPER_MASK ) | ( table[0] & LOWER_MASK );
					table[N-1] = table[M-1] ^ ( y>>1 ) ^ mag01[y&0x1];
				}
				tableIndex = 0;
			}

			y = table[tableIndex++];
			y ^= TEMPERING_SHIFT_U(y);
			y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
			y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
			y ^= TEMPERING_SHIFT_L(y);
			return y;
		}

	}
}

同じ定数やメソッドがコピペして使われてるあたり最悪だけど、テスト用に作っただけだから気にしないで///
もう少しスマートに書けないものかなぁ。

関連リンク

• ライセンス - Plus Le Blog