HDU4767 Bell 中国剩余定理 贝尔数 第二类斯特灵数

这道题要是出在现场赛 估计做出的人要少很多很多，要查阅很多资料才做的出来，至少对于我这个菜狗来说 必须要查资料

http://zh.wikipedia.org/wiki/%E6%96%AF%E7%89%B9%E7%81%B5%E6%95%B0

首先要了解第二类斯特灵数，还有贝尔数，

接下来是解题的关键：

若p为任意质数，则有Bell(n+p) = Bell(n) + Bell(n+1) (mod  p)；

我们来分解题目给的模的数 95041567=31*37*41*43*47，这五个数都是素数，所以可以用上述的式子来解题目，同时又因为 这五个数两两互素，所以可以利用中国剩余定理来解 x=a1(mod m1)这个同余方程组，求出 mod（m1*m2...*mi）d的值，这样就可以得到答案了

但是在使用中国剩余定理之前 我们要求出 Bell（n）%各个素数的 值，这里需要矩阵快速幂

举个例子 若 素数为5此时可以构造这样一个矩阵

0 1 0 0 0
0 1 1 0 0
0 0 1 1 0
0 0 0 1 1
1 0 0 0 1

再列举一个3的矩阵，你就会发现一定的规律，这样就可以求出 Bell(n)%各个素数的值了

#include<iostream>
#include<cstdio>
#include<list>
#include<algorithm>
#include<cstring>
#include<string>
#include<queue>
#include<stack>
#include<map>
#include<vector>
#include<cmath>
#include<memory.h>
#include<set>

#define ll long long
#define LL __int64
#define eps 1e-8

//const ll INF=9999999999999;

#define inf 0xfffffff

using namespace std;

//vector<pair<int,int> > G;
//typedef pair<int,int> P;
//vector<pair<int,int>> ::iterator iter;
//
//map<ll,int>mp;
//map<ll,int>::iterator p;
//
//vector<int>G[30012];

int MOD[5]={31,37,41,43,47};

int a[12][62][62],b[12][62];
int c[62];

int n;

void getstrling()//获取第二类斯特灵数，可以看维基百科
{
	for(int i=0;i<5;i++)
		a[i][0][0]=1;
	for(int i=1;i<=50;i++)
	{
		for(int j=0;j<5;j++)
			a[j][i][1]=1;
		for(int j=1;j<=i;j++)
		{
			for(int k=0;k<5;k++)
				a[k][i][j]=(a[k][i-1][j-1]+j*a[k][i-1][j])%MOD[k];
		}
		for(int j=0;j<5;j++)
		{
			b[j][i]=0;
			for(int k=1;k<=i;k++)
				b[j][i]=(b[j][i]+a[j][i][k])%MOD[j];
		}
	}
}

int Matrixquick(int id,int mod)
{
	int x[62],xx[62],y[62][62],yy[62][62];
	int len=MOD[id]+1;
	for(int i=1;i<=len+1;i++)
		for(int j=1;j<=len+1;j++)
			y[i][j]=0;
	for(int i=1;i<len;i++)
		y[i+1][i]=1;
	y[2][len]=y[3][len]=1;
	for(int i=1;i<=len;i++)
		x[i]=b[id][i];//构造矩阵的过程，这个矩阵有点规律的，
	if(MOD[id]>=n)
		return x[n];
	int tempn=n-1;
	while(tempn)//再写一个函数反而觉得麻烦 直接讲快速幂套在里面了
	{
		if(tempn&1)
		{
			for(int i=1;i<=len;i++)
			{
				xx[i]=0;
				for(int j=1;j<=len;j++)
					xx[i]+=x[j]*y[j][i];
			}
			for(int i=1;i<=len;i++)
				x[i]=xx[i]%mod;
		}
		for(int i=1;i<=len;i++)
		{
			for(int j=1;j<=len;j++)
			{
				yy[i][j]=0;
				for(int k=1;k<=len;k++)
					yy[i][j]+=y[i][k]*y[k][j];
				yy[i][j]%=mod;
			}
		}
		for(int i=1;i<=len;i++)
			for(int j=1;j<=len;j++)
				y[i][j]=yy[i][j];
		tempn>>=1;
	}
	return x[1];
}

int exgcd(int a,int &x,int b,int &y)
{
	if(b==0)
	{
		x=1;
		y=0;
		return a;
	}
	int r=exgcd(b,x,a%b,y);
	int tmp=x;
	x=y;
	y=tmp-a/b*y;
	return r;
}

int china()//中国剩余定理
{
	int M=1;
	int ans=0;
	for(int i=0;i<5;i++)
		M*=MOD[i];
	for(int i=0;i<5;i++)
	{
		int Mi=M/MOD[i];
		int x0,y0;
		int gcd=exgcd(Mi,x0,MOD[i],y0);
		ans=(ans+x0*Mi*c[i])%M;
	}
	return (ans+M)%M;
}

int main(void)
{
	getstrling();
	int t;
	scanf("%d",&t);
	while(t--)
	{
		scanf("%d",&n);
		for(int i=0;i<5;i++)
			c[i]=Matrixquick(i,MOD[i]);
		int ans=china();
		printf("%d\n",ans);
	}
}