题目:https://www.lydsy.com/JudgeOnline/problem.php?id=3501
用贝尔三角形 p^2 地预处理 p 以内的贝尔数。可以模(mod-1)(它是每个分解下的质因子的倍数,所以不影响分开算的时候)。
用公式:\( Bell[n+p^{m}]=m*Bell[n]+Bell[n+1] (mod p) \) \( Bell[n+p]=Bell[n]+Bell[n+1] (mod p) \) 把 n 看成 p 进制,O( p^2 * log m ) 地算。
大概就是从低位走到高位。一开始自己的 b 数组是 Bell[ 0 ] ~ Bell[ p ] ;枚举每一个 p 进制位(从第二位,即 p1 开始),在该位上枚举从1到d[ i ],做一次让角标 + pi 的操作;
这样做完,自己的 b 数组存的就是 Bell[ d[m-1]*pm-1+d[m-2]*pm-2+...+0 ] ~ Bell[ d[m-1]*pm-1+d[m-2]*pm-2+...+p ] 的值。只要输出 b[ d[0] ] 就行了。
借鉴Claris的模板。http://www.cnblogs.com/clrs97/p/4714467.html
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #define ll long long using namespace std; const int mod=999999599,M=mod-1,N=7283; ll n,m; int p[5]={2,13,5281,7283},ans,f[N+5],s[2][N+5]; void upd(int &x,int md){x>=md?x-=md:0;} int pw(int x,int k,int md) {int ret=1;while(k){if(k&1)ret=(ll)ret*x%md;x=(ll)x*x%md;k>>=1;}return ret;} int calc(ll n,int p) { if(n<=N)return f[n]%p; int b[N+5],c[N+5],d[65],lm=0; for(int i=0;i<=p;i++)b[i]=f[i]%p; while(n)d[lm++]=n%p,n/=p; for(int i=1;i<lm;i++) for(int j=1;j<=d[i];j++) { for(int k=0;k<p;k++)c[k]=(i*b[k]+b[k+1])%p; c[p]=c[0]+c[1];upd(c[p],p); for(int k=0;k<=p;k++)b[k]=c[k]; } return b[d[0]]; } int main() { int i,j;bool fx; f[0]=s[1][0]=1; for(i=1,fx=0;i<=N;i++,fx=!fx)//i=1,len=2(0~i) for(f[i]=s[fx][0]=s[!fx][i-1],j=1;j<=i;j++) s[fx][j]=s[!fx][j-1]+s[fx][j-1],upd(s[fx][j],M);//%M?its lcm so ok scanf("%lld%lld",&n,&m); for(i=0;i<4;i++) ans=(ans+(ll)(M/p[i])*pw(M/p[i],p[i]-2,p[i])%M*calc(n,p[i]))%M; printf("%d\n",pw(m%mod,ans,mod)); return 0; }