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

痛定思痛。 2022-09-17 06:25 182阅读 0赞

2013年长春网络赛的题。

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

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

对于贝尔数，如果p是任意质数，则有B\[p+n\]=B\[n\]+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][N][N];  //第二类斯特林数
    int bell[5][N];   //贝尔数
    int data[5];
    int mod[] = {31, 37, 41, 43, 47};
    int modd;   //矩阵计算时的模值
    
    class Matrix
    {
    public:
    	int a[N][N];
    	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[n]%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][n];
                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;
    }