hdu 4767 Bell

来源：互联网发布：初中英语单词大全软件编辑：程序博客网时间：2024/05/17 04:27

题目：hdu 4767 Bell

思路：两个公式 B(n+p)%p=(B(n)+B(n+1))%p B(p^m+n)%p=(m*B(n)+B(n+1))%p

首先明确的是，对于取模数，质因数分解，直接套用模板了。

#include <cstring>#include <iostream>#include <cmath>#include <algorithm>#include <cstdlib>#include <cstdio>#include <map>using namespace std;#define Times 10typedef long long LL;map<LL,int>mp;LL random(LL n){    return LL ((double)rand()/RAND_MAX*n+0.5);}//避免出现了相乘取余前造成的溢出LL multi(LL a,LL b,LL mod){    LL ans=0;    while(b)    {        if(b&1)        {            b--;            ans=(ans+a)%mod;        }        else        {            b/=2;            a=(a+a)%mod;        }    }    return ans;}LL Pow(LL a,LL b,LL mod){    LL ans=1;    while(b)    {        if(b&1)        {            b--;            ans=multi(ans,a,mod);        }        else        {            b/=2;            a=multi(a,a,mod);        }    }    return ans;}// 二次探测找到一个不是x=1或者x=p-1的解就一定不是素数bool witness(LL a,LL n){    LL d=n-1;    while(!(d&1))//幂是偶数时有探测的必要        d>>=1;    LL t=Pow(a,d,n);    while(d!=n-1 && t!=1 && t!=n-1)//如果存在一个非1或非p-1的解，就不是素数    {        t=multi(t,t,n);        d<<=1;    }    return t==n-1 || d&1 ;}bool miller_rabin(LL n){    if(n==2)        return true;    if(n<2 || !(n&1))        return false;    for(int i=1;i<=Times;i++)    {        LL a=random(n-2)+1;        if(!witness(a,n))            return false;    }    return true;}LL gcd(LL a,LL b){    if(b==0)        return a;    return gcd(b,a%b);}LL pollard_rho(LL n,int c){    LL x,y,d,i=1,k=2;    x=random(n-2)+1;    y=x;    while(1)    {        i++;        x=(multi(x,x,n)+c)%n;        d=gcd(y-x,n);        if(1<d && d<n)            return d;        if(y==x)            return n;        if(i==k)        {            y=x;            k<<=1;        }    }}void find(LL n,LL c){    if(n==1)        return ;    if(miller_rabin(n))    {        mp[n]++;        return ;    }    LL p=n;    while(p>=n)        p=pollard_rho(p,c--);    find(p,c);    find(n/p,c);}void Print(LL x,int num){    while(num--)        printf("%I64d ",x);}int main(){    long long p=95041567;    mp.clear();    find(p,20131001);    map<LL,int>::iterator it=mp.begin();    for(;it!=mp.end();it++)        Print(it->first,it->second);    printf("\n");    return 0;}

取模数可以质因数分解为 31*37*41*43*47，肯定是分成五个，然后用中国剩余定理合并的。

开始做的时候，用的是第二个公式，结果发现，递归层数太多，会超时。

// TLE 代码#include <cstring>#include <iostream>#include <cmath>#include <algorithm>#include <cstdio>using namespace std;#pragma comment(linker, "/STACK:102400000,102400000")typedef __int64 LL;LL w[5]={31,37,41,43,47};LL a[5]={0,0,0,0,0};LL MOD = 95041567 ;LL Bell[50][5];LL C[50][50][5];LL exgcd(LL a,LL b,LL &x,LL &y){    if(b==0)    {        x=1;        y=0;        return a;    }    LL ans=exgcd(b,a%b,x,y);    LL t=x;    x=y;    y=t-a/b*y;    return ans;}LL CRT(int r){    LL M=MOD;    LL i,d,x0,y0,ans=0;    for(int i=0;i<r;i++)    {        d=M/w[i];        exgcd(d,w[i],x0,y0);        ans=(ans+d*x0*a[i])%M;    }    while(ans<=0)        ans+=M;    return ans;}void init(){    memset(C,0,sizeof(C));    for(int k=0;k<5;k++) // C[i][j][k] 表示C(i,j)%w[k]    {        C[0][0][k]=1;        for(int i=0;i<50;i++)        {            C[i][0][k]=C[i][i][k]=1;            for(int j=1;j<i;j++)                C[i][j][k]=(C[i-1][j][k]+C[i-1][j-1][k])%w[k];        }    }    // 预处理出前50项分别取模的大小    for(int k=0;k<5;k++)    {        Bell[0][k]=1;        for(int i=1;i<50;i++) //复杂度10^3左右        {            Bell[i][k]=0;            for(int j=0;j<i;j++)                Bell[i][k]=(Bell[i][k]+Bell[j][k]*C[i-1][j][k]%w[k])%w[k];        }    }}LL tmp_Bell(LL n,int p) //返回 Bell[n][p]的值{    if(n<50)        return Bell[n][p];    LL cnt=1;    int mi=0;    while(cnt*w[p]<=n)    {        cnt*=w[p];        mi++;    }    n-=cnt;    return (mi*tmp_Bell(n,p)%w[p]+tmp_Bell(n+1,p))%w[p];}LL BellNumber(LL n){    for(int i=0;i<5;i++)    {        a[i]=tmp_Bell(n,i);    }    return CRT(5);}int main(){    init();    LL n;    int t;    scanf("%d",&t);    while(t--)    {        scanf("%I64d",&n);        printf("%I64d\n",BellNumber(n)%MOD);    }    return 0;}

上面的代码，很容易就卡死了，所以还是得看看第一个公式。

对于第一个公式，我们可以构造一个大小为当前p*p的矩阵。

构造矩阵后，很快就可以得到分别对各种p的取模值，然后用中国剩余定理合并。

详见代码。

#include <cstring>#include <iostream>#include <cmath>#include <algorithm>#include <cstdio>using namespace std;#pragma comment(linker, "/STACK:102400000,102400000")typedef __int64 LL;LL w[5]={31,37,41,43,47};LL a[5]={0,0,0,0,0};LL MOD = 95041567 ;LL Bell[50][5];LL C[50][50][5];struct Matrix{    LL m[50][50];}E,D;Matrix Multi(Matrix A,Matrix B,int m,int n,int r,int p){    Matrix ans;    for(int i=1;i<=m;i++)        for(int j=1;j<=n;j++)        {            ans.m[i][j]=0;            for(int k=1;k<=r;k++)                ans.m[i][j]=(ans.m[i][j]+A.m[i][k]*B.m[k][j])%w[p];        }    return ans;}Matrix Pow(Matrix A,LL k,int n,int p){    Matrix ans=E;    while(k)    {        if(k&1)        {            k--;            ans=Multi(ans,A,n,n,n,p);        }        else        {            k/=2;            A=Multi(A,A,n,n,n,p);        }    }    return ans;}LL exgcd(LL a,LL b,LL &x,LL &y){    if(b==0)    {        x=1;        y=0;        return a;    }    LL ans=exgcd(b,a%b,x,y);    LL t=x;    x=y;    y=t-a/b*y;    return ans;}LL CRT(int r){    LL M=MOD;    LL i,d,x0,y0,ans=0;    for(int i=0;i<r;i++)    {        d=M/w[i];        exgcd(d,w[i],x0,y0);        ans=(ans+d*x0*a[i])%M;    }    while(ans<=0)        ans+=M;    return ans;}void init(){    for(int i=1;i<50;i++)        for(int j=1;j<50;j++)            E.m[i][j]=(i==j);    memset(C,0,sizeof(C));    for(int k=0;k<5;k++) // C[i][j][k] 表示C(i,j)%w[k]    {        C[0][0][k]=1;        for(int i=0;i<50;i++)        {            C[i][0][k]=C[i][i][k]=1;            for(int j=1;j<i;j++)                C[i][j][k]=(C[i-1][j][k]+C[i-1][j-1][k])%w[k];        }    }    // 预处理出前50项分别取模的大小    for(int k=0;k<5;k++)    {        Bell[0][k]=1;        for(int i=1;i<50;i++) //复杂度10^3左右        {            Bell[i][k]=0;            for(int j=0;j<i;j++)                Bell[i][k]=(Bell[i][k]+Bell[j][k]*C[i-1][j][k]%w[k])%w[k];        }    }}LL tmp_Bell(LL n,int p) //返回 Bell[n][p]的值{    if(n<50)        return Bell[n][p];    Matrix tmp;    memset(tmp.m,0,sizeof(tmp.m));    for(int i=1;i<w[p];i++)        tmp.m[i][i]=1,tmp.m[i][i+1]=1;    tmp.m[w[p]][w[p]]=1;    tmp.m[w[p]][1]=1;    tmp.m[w[p]][2]=1;    LL cnt=n/w[p];    n-=w[p]*cnt;    Matrix ans;    for(int i=1;i<=w[p];i++)        ans.m[i][1]=Bell[i][p];    ans=Multi(Pow(tmp,cnt,w[p],p),ans,w[p],1,w[p],p);    return ans.m[n][1];}LL BellNumber(LL n){    for(int i=0;i<5;i++)    {        a[i]=tmp_Bell(n,i);    }    return CRT(5);}int main(){    init();    LL n;    int t;    scanf("%d",&t);    while(t--)    {        scanf("%I64d",&n);        printf("%I64d\n",BellNumber(n)%MOD);    }    return 0;}