bzoj 3501 PA2008 Cliquers Strike Back——贝尔数

题目：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 进制位（从第二位，即 p¹ 开始），在该位上枚举从1到d[ i ]，做一次让角标 + pⁱ 的操作；

　　这样做完，自己的 b 数组存的就是 Bell[ d[m-1]*p^m-1+d[m-2]*p^m-2+...+0 ] ~ Bell[ d[m-1]*p^m-1+d[m-2]*p^m-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
",pw(m%mod,ans,mod));
  return 0;
}