2013长春网赛1009 hdu 4767 Bell（矩阵快速幂+中国剩余定理）

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=4767

题意：求集合{1, 2, 3, ..., n}有多少种划分情况bell[n]，最后结果bell[n] mod 95041567.

分析：首先了解三个概念：贝尔数   第二类斯特灵数   中国剩余定理

贝尔数是指基数为n的集合的划分方法的数目。

B_0=1,\quad B_1=1,\quad B_2=2,\quad B_3=5,\quad B_4=15,\quad B_5=52,\quad B_6=203,\quad\dots

贝尔数适合递推公式：

B_{n+1}=\sum_{k=0}^{n}{{n \choose k}B_k}.

每个贝尔数都是"第二类Stirling数"的和

B_n=\sum_{k=1}^n S(n,k).

贝尔数满足两个公式：（p为质数）

            1) B[n+p] = B[n] + B[n+1] （mod p） ;

            2) B[p^m+n] = m*B[n] + B[n+1] （mod p） .

将95041567质因数分解发现，95041567 = 31*37*41*43*47

所以B[n]%95041567可以分解为 B[n]%p(p=31,41,43,47)，

我们可以先求出B[n] mod p[i]的值a[i],这样问题就转化为 X=a[i](mod p[i])，

很明显这是几个一次同余方程，最后用中国剩余定理合并就可以了。

那要怎么求B[n]%p[i]呢？

利用上面的公式（1），我们发现这是 一个递推式，所以可以利用矩阵法来求解。

我们可以构造一个大小为当前p*p的矩阵。

这样我们就可以求出任意的B[n]了

我们可以先用贝尔数递推公式B_{n+1}=\sum_{k=0}^{n}{{n \choose k}B_k}.求出前50个贝尔数，因为p[i]<50，所以对于大于p[i]的贝尔数，由上面的矩阵法可以求得。

比如：| 1 1 0 0 .... 0 0 |      |   B[1]  |      | B[1+p] |

        | 0 1 1 0 .... 0 0 |      |   B[2]  |      | B[2+p] |

        | 0 0 1 1 .... 0 0 |      |   B[3]  |      | B[3+p] |

        | ....  .... .... ....  |  *  |   .....   |  =  |    .....   |

        | 0 0 0 0 .... 1 1 |      | B[p-1] |      | B[2p-1]|

        | 1 1 0 0 .... 0 1 |      |   B[p]  |      | B[2p]   |

若n=i+p,则只需求一次A*C=D，然后输出D[n-p]即D[i]就行了，

比如p[0]=31，如要求B[32]=B[1+31]，只需求一次A*C=D,然后输出D[1]，求B[51]则输出D[20]。

那么

若n=i+p^m,这只需求A^m*C=D，然后输出D[i]即可

到此算法结束

总结一下：

先用贝尔数递推公式求出前50个贝尔数，然后用矩阵快速幂分别求出bell[n]%p[i]赋给a[i]，然后用中国剩余定理合并结果就可以求出bell[n]%95041567了。

AC代码如下：

  1 #include<stdio.h>
  2 #include<string.h>
  3 #define LL __int64
  4 int w[5]={31,37,41,43,47};
  5 LL c[50][50][5],bell[50][5];
  6 LL a[5];
  7 struct Matrix{
  8     LL m[50][50];
  9 };
 10 void init()
 11 {
 12     int i,j,k;
 13     for(k=0;k<5;k++)
 14     {
 15         c[0][0][k]=1;
 16         for(i=1;i<50;i++)     //c[i][j][k] 表示c(i,j)%w[k]
 17         {
 18             c[i][0][k]=c[i][i][k]=1;
 19             for(j=1;j<i;j++)
 20                 c[i][j][k]=(c[i-1][j-1][k]+c[i-1][j][k])%w[k];
 21         }
 22     }
 23     // 预处理出前50项分别取模的大小 
 24     for(k=0;k<5;k++)
 25     {
 26         bell[0][k]=1;
 27         for(i=1;i<50;i++)
 28         {
 29             bell[i][k]=0;
 30             for(j=0;j<i;j++)
 31             {
 32                 bell[i][k]=(bell[i][k]+c[i-1][j][k]*bell[j][k])%w[k];
 33             }
 34         }
 35     }
 36 }
 37 LL exgcd(LL a,LL b,LL &x,LL &y)    //扩展欧几里得
 38 {
 39     if(b==0)
 40     {
 41         x=1; y=0;
 42         return a;
 43     }
 44     LL d=exgcd(b,a%b,x,y);
 45     LL t=x;
 46     x=y;
 47     y=t-a/b*y;
 48     return d;
 49 }
 50 LL CRT(int k)     //中国剩余定理
 51 {
 52     LL i,N=1,ans=0;
 53     LL t,x,y,d;
 54     for(i=0;i<k;i++)
 55         N*=w[i];
 56     for(i=0;i<k;i++)
 57     {
 58         t=N/w[i];
 59         d=exgcd(t,w[i],x,y);
 60         ans=(ans+x*t*a[i])%N;
 61     }
 62     return (ans+N)%N;
 63 }
 64 Matrix mul(Matrix x,Matrix y,int n,int mod)      //矩阵乘法
 65 {
 66     int i,j,k;
 67     Matrix tmp;
 68     memset(tmp.m,0,sizeof(tmp.m));
 69     for(i=1;i<=n;i++)
 70         for(j=1;j<=n;j++)
 71             for(k=1;k<=n;k++)
 72             {
 73                 tmp.m[i][j]+=(x.m[i][k]*y.m[k][j])%mod;
 74                 tmp.m[i][j]%=mod;
 75             }
 76     return tmp;
 77 }
 78 Matrix quickpow(Matrix res,Matrix A,int k,int n,int mod)   //矩阵快速幂模
 79 {
 80     int i;
 81     memset(res.m,0,sizeof(res.m));
 82     for(i=1;i<=n;i++)
 83         res.m[i][i]=1;
 84     while(k)
 85     {
 86         if(k&1)
 87             res=mul(res,A,n,mod);
 88         A=mul(A,A,n,mod);
 89         k>>=1;
 90     }
 91     return res;
 92 }
 93 LL BellNumber(int n,int p)    //求bell[n][p]即bell[n]%w[p]
 94 {
 95     int i;
 96     if(n<50)
 97         return bell[n][p];
 98     Matrix A,origin,ans;
 99     memset(A.m,0,sizeof(A.m));
100     memset(origin.m,0,sizeof(origin.m));
101     for(i=1;i<w[p];i++)
102         A.m[i][i]=A.m[i][i+1]=1;
103     A.m[w[p]][1]=1;
104     A.m[w[p]][2]=1;
105     A.m[w[p]][w[p]]=1;
106     for(i=1;i<=w[p];i++)
107         origin.m[i][1]=bell[i][p];
108     LL cnt=n/w[p];
109     n=n-w[p]*cnt;
110     Matrix res;
111     res=quickpow(res,A,cnt,w[p],w[p]);
112     ans=mul(res,origin,w[p],w[p]);
113     return ans.m[n][1];
114 }
115 void solve(int n)
116 {
117     int i;
118     for(i=0;i<5;i++)
119         a[i]=BellNumber(n,i);    //bell[n]mod各个质数的值
120     LL ans=CRT(5);
121     printf("%I64d\n",ans);
122 }
123 int main()
124 {
125     int t,n;
126     init();
127     scanf("%d",&t);
128     while(t--)
129     {
130         scanf("%d",&n);
131         solve(n);
132     }
133     return 0;
134 }
View Code

 

posted on 2013-10-05 00:05  jumpingfrog0  阅读(559)  评论(0)  编辑  收藏  举报

导航