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.
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)
for each case:
n (1 <= n < 2^31)
Output
for each case:
bell[n] % 95041567
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!