bzoj 3501 PA2008 Cliquers Strike Back —— 贝尔数

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