题目：https://www.lydsy.com/JudgeOnline/problem.php?id=3501
用贝尔三角预处理贝尔数，拆模数并在 ( p ) 进制下使用公式，因为这样每次角标增加的是 ( p^{k} )；
循环使用数组非常优美！0~p 的角标背后是许多 ( p ) 的整数次幂，而角标那个数字是它的 ( p^0 ) 上的数，所以最后取 ( b[d[0]] )；
Claris 写得太好了！http://www.cnblogs.com/clrs97/p/4714467.html
代码如下：
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> using namespace std; typedef long long ll; int const xn=8005,mod=999999598; int p[4]={2,13,5281,7283},s[2][xn],f[xn],b[xn],c[xn],d[65]; ll n,m; ll pw(ll a,int b,int md) { ll ret=1; a=a%md; for(;b;b>>=1,a=(a*a)%md)if(b&1)ret=(ret*a)%md; return ret; } int upt(int x){while(x>=mod)x-=mod; while(x<0)x+=mod; return x;} void init() { f[0]=f[1]=s[0][0]=1; s[0][1]=2; for(int i=2,x=1,j;i<=p[3];i++,x^=1) for(j=1,f[i]=s[x][0]=s[x^1][i-1];j<=i;j++) s[x][j]=upt(s[x^1][j-1]+s[x][j-1]); } int cal(ll n,int p) { for(int i=0;i<=p;i++)b[i]=f[i]%p;//<= int m=0; while(n)d[m++]=n%p,n/=p; m--; for(int i=1;i<=m;i++)//i=1 -> d:0~m for(int j=1;j<=d[i];j++)//d[i]*p^k { for(int k=0;k<p;k++)c[k]=(i*b[k]+b[k+1])%p; c[p]=(c[0]+c[1])%p;//c[p]! for(int k=0;k<=p;k++)b[k]=c[k];//<=!!! } return b[d[0]];//id -> last digit } int main() { scanf("%lld%lld",&n,&m); init(); if(n<=p[3]){printf("%lld ",pw(m,f[n],mod+1)); return 0;}// ll ans=0; for(int i=0;i<=3;i++) { int r=cal(n,p[i]),w=mod/p[i]; ans=(ans+(ll)w*pw(w,p[i]-2,p[i])%mod*r)%mod;//w //%p[i] } printf("%lld ",pw(m,ans,mod+1)); return 0; }