Bell - HDU 4767 贝尔数

来源:互联网 发布:如何用ps制作淘宝店招 编辑:程序博客网 时间:2024/03/28 20:11

Bell

Time Limit: 6000/3000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 600    Accepted Submission(s): 252


Problem 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
6123456
 

Sample Output
1251552203
 


思路:参考http://fudq.blog.163.com/blog/static/19135023820139114642927。在这里也给出我的一份AC代码供大家参考。

AC代码如下:

#include<cstdio>#include<cstring>using namespace std;typedef long long ll;struct node{    ll f[50][50];};int T,t,n,MOD=95041567;ll num[55][55],bell[55],mod;ll prime[6]={0,31,37,41,43,47};ll M[10];node first,ret;node mul(node a,node b){    int i,j,k;    node m;    for(i=0;i<mod;i++)       for(j=0;j<mod;j++)       {           m.f[i][j]=0;           for(k=0;k<mod;k++)              m.f[i][j]=(m.f[i][j]+a.f[i][k]*b.f[k][j])%mod;       }    return m;}void init(){    int i,j,k;    num[1][1]=1;    for(i=2;i<=50;i++)       for(j=1;j<=i;j++)          num[i][j]=(num[i-1][j]*j+num[i-1][j-1])%MOD;    bell[0]=1;    for(i=1;i<=50;i++)       for(j=1;j<=i;j++)          bell[i]=(bell[i]+num[i][j])%MOD;}void solve(int pos){    int i,j,k;    mod=prime[pos];    memset(first.f,0,sizeof(first.f));    memset(ret.f,0,sizeof(ret.f));    for(i=0;i<mod;i++)       ret.f[0][i]=bell[i]%mod;    first.f[mod-1][0]=1;    i=0;    for(j=1;j<mod;j++)    {        first.f[i][j]=first.f[i+1][j]=1;        i++;    }    k=n/(mod-1);    while(k)    {        if(k&1)          ret=mul(ret,first);        first=mul(first,first);        k/=2;    }    M[pos]=ret.f[0][n%(mod-1)];}void ex_gcd(ll a,ll b,ll &d,ll &x,ll &y){    if(!b)    {        x=1;        y=0;        d=a;    }    else    {        ex_gcd(b,a%b,d,y,x);        y-=x*(a/b);    }}ll ex_crt(ll *m,ll *r,int n){    ll M=m[1],R=r[1],x,y,d;    int i;    for(i=2;i<=n;i++)    {        ex_gcd(M,m[i],d,x,y);        if((r[i]-R)%d)          return -1;        x=(r[i]-R)/d*x%(m[i]/d);        R+=x*M;        M=M/d*m[i];        R%=M;    }    return R>0 ? R : R+M;}int main(){    int i,j,k;    init();    ll ans;    scanf("%d",&T);    for(t=1;t<=T;t++)    {        scanf("%d",&n);        for(i=1;i<=5;i++)           solve(i);        ans=ex_crt(prime,M,5);        ans%=MOD;        printf("%I64d\n",ans);    }}


0 0
原创粉丝点击