[ Mersenne Twister ] 梅森旋轉算法

梅森旋轉算法是一種可以快速產生高級偽隨機數的算法，修正了很多以前的隨機速算法的缺陷(Ex:線性同於算法)。

這是專門用在蒙地卡羅測試用的，不要拿它來做持久化隨機二叉數的亂數產生器，會太慢

最為廣泛使用Mersenne Twister的一種變體是MT19937，可以產生32位整數序列，從C++11開始，C++也可以使用這種算法。在Boost C++,Glib和NAG數值庫中，作為插件提供。但是一般在競賽時可能會有不支援的情況發生，所以將它做成模板當作參考。

MT19937_64則是可以產生64位整數序列

以下提供模板:

	#ifndef SUNMOON_MERSENNE_TWISTER
	#define SUNMOON_MERSENNE_TWISTER
	template<typename T,T w,int n,T m,T r,T a,T u,T d,T s,T b,T t,T c,T l,T f>
	class mersenne_twister{
	private:
	int index;
	const T lower_mask=(1ull<<r)-1ull;
	const T upper_mask=~lower_mask;
	T mt[n];
	inline void twist(){
	for(int i=0;i<n;++i){
	T y=(mt[i]&upper_mask)+(mt[(i+1)%n]&lower_mask);
	mt[i]=mt[(i+m)%n]^y>>1;
	if(y%2!=0)mt[i]=mt[i]^a;
	}
	index=0;
	}
	public:
	mersenne_twister(T seed):index(n){
	mt[0]=seed;
	for(int i=1;i<n;++i)mt[i]=f*(mt[i-1]^mt[i-1]>>(w-2))+i;
	}
	inline T operator()(){
	if(index>=n)twist();
	T y=mt[index++];
	y=y^y>>u&d;
	y=y^y<<s&b;
	y=y^y<<t&c;
	y=y^y>>l;
	return y;
	}
	};
	typedef mersenne_twister<unsigned,32,624,397,31,0x9908b0dfUL,11,
	0xffffffffUL,7,0x9d2c5680UL,15,0xefc60000UL,18,1812433253UL> MT19937;
	typedef mersenne_twister<unsigned long long,64,312,156,31,0xb5026f5aa96619e9ULL,29,0x5555555555555555ULL
	,17,0x71d67fffeda60000ULL,37,0xfff7eee000000000ULL,43,6364136223846793005ULL> MT19937_64;
	#endif

view raw mersenne_twister.cpp hosted with ❤ by GitHub

[ chinese remainder theorem ] 中國剩餘定理

我是看維基百科實現的，所以就直接貼上模板
注意:

• crt函數的m[i]是模數
• crt函數的a[i]是原本數模m[i]後的值
• 必須要保證m兩兩互質
• 容易溢位

模板:

	#ifndef CHINESE_REMAINDER_THEOREM
	#define CHINESE_REMAINDER_THEOREM
	template<typename T>
	inline T Euler(T n){
	T ans=n;
	for(T i=2;i*i<=n;++i){
	if(n%i==0){
	ans=ans/i*(i-1);
	while(n%i==0)n/=i;
	}
	}
	if(n>1)ans=ans/n*(n-1);
	return ans;
	}
	template<typename T>
	inline T pow_mod(T n,T k,T m){
	T ans=1;
	for(n=(n>=m?n%m:n);k;k>>=1){
	if(k&1)ans=ans*n%m;
	n=n*n%m;
	}
	return ans;
	}
	template<typename T>
	inline T crt(std::vector<T> &m,std::vector<T> &a){
	T M=1,tM,ans=0;
	for(int i=0;i<(int)m.size();++i)M*=m[i];
	for(int i=0;i<(int)a.size();++i){
	tM=M/m[i];
	ans=(ans+(a[i]*tM%M)*pow_mod(tM,Euler(m[i])-1,m[i])%M)%M;
	/*如果m[i]是質數，Euler(m[i])-1=m[i]-2，就不用算Euler了*/
	}
	return ans;
	}
	#endif

view raw Chinese remainder theorem.cpp hosted with ❤ by GitHub