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

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#elsefreopen("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;}