Hdu4767 Bell (贝尔数 中国剩余定理 构造矩阵)

2013年长春网络赛的题。

比赛时很快就找到了bell数的性质，看到数据范围也马上想到要用矩阵做，只可惜当时没接触过中国剩余定理……

以下题解转自：http://fudq.blog.163.com/blog/static/19135023820139114642927

对于贝尔数，如果p是任意质数，则有B[p+n]=B
+B(n+1)%p

还有可以发现的是95041567=31*37*41*43*47。

利用以上信息我们便可解此题。

如果能够得到B(n)%31，B(n)%37，B(n)%41，B(n)%43，B(n)%47的结果，我们便可以用中国剩余定理求出B(n)%95041567的结果

现在问题转换成求B(n)%w[i]，由于w[i]是质数，便可利用上面Bell数的性质。

举个例子，如果w[i]=5，可以由B0 B1 B2 B3 B4得到B5 B6 B7 B8，又由B4 B5 B6 B7 B8得到B9 B10 B11 B12，依次类推。

这时可以利用矩阵二次幂快速求解，构造一个如下所示的01矩阵：

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

便可以由B0 B1 B2 B3 B4得到B4 B5 B6 B7 B8。

对于一个n，需要将构造矩阵乘上n/(w[i]-1)次，第n%(w[i]-1)列的结果便是B(n)%w[i]。

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

const int N=55;

int stir2[5]

;  //第二类斯特林数
int bell[5]
;   //贝尔数
int data[5];
int mod[] = {31, 37, 41, 43, 47};
int modd;   //矩阵计算时的模值

class Matrix
{
public:
	int a

;
	int n;   //矩阵大小

	void init (int x)
	{
		memset(a,0,sizeof(a));
		if (x) for (int i=0;i<N;i++)
			a[i][i]=1;
    }

	Matrix operator * (Matrix b)
	{
		Matrix p;
		p.n=b.n;   //
		p.init(0);
		for (int i=0;i<n;i++) for (int j=0;j<n;j++)
			for (int k=0;k<n;k++)
				(p.a[i][j]+=a[i][k]*b.a[k][j]%modd)%=modd;
		return p;
	}

	Matrix power (int t)
	{
		Matrix ans,p=*this;
		ans.n=p.n;  //
		ans.init(1);
		while (t)
		{
			if (t&1) ans=ans*p;
			p=p*p;
			t>>=1;
		}
		return ans;
	}
}a;

void Init () //预处理50以下的贝尔数
{
    int i,j,k;
    for (i=0;i<5;i++)
    {
        stir2[i][0][0]=1;
        bell[i][0]=bell[i][1]=1;
    }

    for (i=1;i<=50;i++)
    {
        for (j=0;j<5;j++)
            stir2[j][i][1]=stir2[j][i][i]=1;
        for (j=1;j<=i;j++) for (k=0;k<5;k++)
            stir2[k][i][j]=(stir2[k][i-1][j-1]+j*stir2[k][i-1][j])%mod[k];
    }
    for (i=1;i<=50;i++)
        for (k=0;k<5;k++)
        {
            bell[k][i]=0;
            for (j=0;j<=i;j++)
            bell[k][i]=(bell[k][i]+stir2[k][i][j])%mod[k];
        }
}

int Extended_Euclid (int a,int b,int &x,int &y)
{//扩展欧几里得算法,求ax+by=gcd(a,b)的一组解(x,y),d=gcd(a,b)
	int d;
	if (b==0)
	{
		x=1;y=0;
		return a;
	}
	d=Extended_Euclid(b,a%b,y,x);
	y-=a/b*x;
	return d;
}

int Chinese_Remainder (int a[],int w[],int len)
{//中国剩余定理  a[]存放余数  w[]存放两两互质的数
	int x,y,m;
	int lcm=1,i,ans=0;
	for (i=0;i<len;i++)
		lcm*=w[i];
	for (i=0;i<len;i++)
	{
		m=lcm/w[i];
		Extended_Euclid(w[i],m,x,y);
		ans=(ans+y*m*a[i])%lcm;
	}
	return (lcm+ans%lcm)%lcm;
}

int Cal (int n,int t)	//构造矩阵a,计算bell
%w[t]
{
	int ans=0,i;
	a.init(0);
	a.n=mod[t];
	modd=mod[t];
	for (i=1;i<mod[t];i++)
		a.a[i][i]=a.a[i-1][i]=1;
	a.a[ mod[t]-1 ][0]=1;
	int p1=n/(mod[t]-1);
	int p2=n%(mod[t]-1);

	a=a.power(p1);
	for (i=0;i<mod[t];i++)
		ans=(ans+bell[t][i]*a.a[i][p2])%mod[t];
	return ans;
}

int main ()
{
#ifdef ONLINE_JUDGE
#else
	freopen("read.txt","r",stdin);
#endif
    int T,n,i;
    scanf("%d",&T);
    Init ();
    while (T--)
    {
        scanf("%d",&n);
        if (n<50)
        {
            for (i=0;i<5;i++)
                data[i]=bell[i]
;
            printf("%d\n",Chinese_Remainder(data,mod,5));
            continue;
        }
        for (i=0;i<5;i++)
            data[i]=Cal(n,i);
        printf("%d\n",Chinese_Remainder(data,mod,5));
    }
	return 0;
}