Bell(hdu4767+矩陣+中國剩余定理+bell數+Stirling數+歐幾裡德)，hdu4767stirling

Bell

Time Limit:3000MS     Memory Limit:32768KB     64bit IO Format:%I64d & %I64u Submit Status Practice HDU 4767  

Description

What? MMM is learning Combinatorics!? 
Looks like she is playing with the bell sequence now: 
bell[n] = number of ways to partition of the set {1, 2, 3, ..., n} 

e.g. n = 3: 
{1, 2, 3} 
{1} {2 3} 
{1, 2} {3} 
{1, 3} {2} 
{1} {2} {3} 
so, bell[3] = 5 

MMM wonders how to compute the bell sequence efficiently.   

Input

T -- number of cases 
for each case: 
　　n (1 <= n < 2^31)   

Output

for each case: 
　　bell[n] % 95041567   

Sample Input

6 1 2 3 4 5 6  

Sample Output

1 2 5 15 52 203    

首先，我很高興終於把這題給解出來了！！！自己有幾個難點一直沒想通，今天無聊的時候一想。然後思路莫名奇妙的順了，一切都變的理所當然。。。然後，我趁熱打鐵，把思路一理，就被我A了！

題意：就是求n個數的非空集合的劃分方法數；

例如n = 3

{1, 2, 3}

{1}  {2, 3}

{1, 2} {3}

{1, 3} {2}

{1} {2} {3}

所以Bell(3) = 5

給你一個n，求Bell(n) % 95041567的值

 

這道題涉及的知識比較多，，，本人覺得有一定的難度系數，沒理解的，建議可以先從入門題刷起；

思路：

　　首先必須了解一個概念貝爾數（來自維基百科）& Stirling數（不懂的點鏈接）；

　　其次，如果你不懂Stirling數，那麼自己動手推一下，其實這個還是簡單的，杭電有一題簡單的題，可以練習下：

　　http://acm.hdu.edu.cn/showproblem.php?pid=2512 相信很多同學都做過了，只是當時不知道是Stirling數，這裡是第二類Stirling數；然後貝爾數呢，就是第二類Stirling數的和；

通過上面這個公式，我們可以利用矩陣計算大於n的貝爾數，方法是利用中國剩余定理，結合歐幾裡德和同余方程；首先我們輸入n(假設，n>p，p是95041567的質因子，eg:n=50)，對n<50的Bell數進行預處理，構造一個p*p的矩陣；

利用公式求出，他們的余數存入數組at[i];質因子m[5]={31,37,41,43,47};利用中國剩余定理；求得余數！

 

轉載請注明出處：尋找&星空の孩子   

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

 

  1 #include<stdio.h>
  2 #include<string.h>
  3 #define LL long long
  4 #define mod 95041567
  5 #define MM 50
  6 
  7 int m[5]={31,37,41,43,47};
  8 int Sti[MM][MM],bell[MM][MM];
  9 int at[5];//各項余數
 10 
 11 void stirling2()
 12 {
 13     memset(Sti,0,sizeof(Sti));
 14     Sti[0][0]=1;
 15     for(int i=1;i<=MM;i++)
 16     {
 17         for(int j=1;j<=i;j++)
 18         {
 19             Sti[i][j]=(int)(Sti[i-1][j-1]+((LL)j*(LL)Sti[i-1][j])%mod)%mod;
 20         }
 21     }
 22 }
 23 void init()
 24 {
 25     stirling2();
 26     for(int i=0;i<5;i++)
 27     {
 28         bell[i][0]=1;
 29         for(int j=1;j<m[i];j++)
 30         {
 31             bell[i][j]=0;
 32             for(int k=1;k<=j;k++)
 33             {
 34                 bell[i][j]=(bell[i][j]+Sti[j][k])%m[i];
 35             }
 36 //            printf("%d\t%d\n",j,bell[i][j]);
 37         }
 38     }
 39 }
 40 struct Matrix
 41 {
 42     int mat[MM][MM];
 43 };
 44 
 45 Matrix multiply(Matrix a,Matrix b,int M)
 46 {
 47     Matrix c;
 48     memset(c.mat,0,sizeof(c.mat));
 49     for(int i=0;i<M;i++)
 50     {
 51         for(int j=0;j<M;j++)
 52         {
 53             if(a.mat[i][j]==0)continue;
 54             for(int k=0;k<M;k++)
 55             {
 56                 if(b.mat[j][k]==0)continue;
 57                 c.mat[i][k]=(c.mat[i][k]+a.mat[i][j]*b.mat[j][k])%M;
 58             }
 59         }
 60     }
 61     return c;
 62 }
 63 Matrix quickmod(Matrix a,int n,int M)
 64 {
 65     Matrix res;
 66 
 67     for(int i=0;i<M;i++)
 68     {
 69         for(int j=0;j<M;j++)
 70             res.mat[i][j]=(i==j);
 71     }
 72 
 73     while(n)
 74     {
 75         if(n&1) res=multiply(res,a,M);
 76         n>>=1;
 77         a=multiply(a,a,M);
 78     }
 79     return res;
 80 }
 81 
 82 int work(int n,int M,int k)
 83 {
 84     Matrix per;//基礎矩陣；
 85     memset(per.mat,0,sizeof(per.mat));
 86     per.mat[0][M-1]=1;
 87 
 88     for(int i=1;i<M;i++)
 89     {
 90             per.mat[i][i]=per.mat[i][i-1]=1;
 91     }
 92 
 93     Matrix tmp=quickmod(per,n/(M-1),M);
 94 
 95     int ans[MM];
 96     for(int i=0;i<M;i++)
 97     {
 98         ans[i]=0;
 99         for(int j=0;j<M;j++)
100         {
101             ans[i]=(ans[i]+bell[k][j]*tmp.mat[i][j])%M;
102         }
103     }
104     return ans[n%(M-1)];
105 }
106 void exgcd(int a,int b,int& d,int& x,int &y)
107 {
108     if(!b){d=a;x=1;y=0;}
109     else
110     {
111         exgcd(b,a%b,d,y,x);
112         y-=x*(a/b);
113     }
114 }
115 int China(int r)
116 {
117     int Mc=1;
118     int Mi,d,x,y,as=0;
119     for(int i=0;i<r;i++)
120     {
121         Mc*=m[i];
122     }
123     for(int i=0;i<r;i++)
124     {
125         Mi=Mc/m[i];
126         exgcd(Mi,m[i],d,x,y);
127         as=(as+Mi*x*at[i])%Mc;
128     }
129     if(as<0) as+=Mc;
130     return as;
131 }
132 int main()
133 {
134     int T,n;
135     scanf("%d",&T);
136 
137     init();
138     while(T--)
139     {
140         scanf("%d",&n);
141         for(int i=0;i<5;i++)
142         {
143             at[i]=work(n,m[i],i);
144         }
145         int sol=China(5);
146         printf("%d\n",sol);
147     }
148 
149     return 0;
150 }

 

 

時間過的好快，居然要期末考了，，，在學校的日子也不多了，aking2015......為了我們共同的夢！fighting！