hdu 4767 Bell

最新推荐文章于 2019-09-18 23:25:00 发布

shiyuankongbu

最新推荐文章于 2019-09-18 23:25:00 发布

阅读量1.7k

点赞数

分类专栏：数学

本文链接：https://blog.csdn.net/shiyuankongbu/article/details/12218859

版权

数学同时被 2 个专栏收录

31 篇文章 0 订阅

订阅专栏

稍微有价值的题

11 篇文章 0 订阅

订阅专栏

题目：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 10
typedef 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;
}