[Bell数] HDU 4767 Bell & BZOJ 3501 PA2008 Cliquers Strike Back

ACdreamer大牛的详细介绍

Bell数 大概是这样的

• 有递推公式 可以分治FFT求
  
• 生成函数很优美 可以用多项式科技求
  
• 是第二类斯特林数的和 可以用Bell三角形预处理
  
• 神奇的同余性质 可以计算对小质数取模的值 也可以CRT合并 这里的p是不大于100的素数
  
  
• 模素数p的周期
  

这里就是求bell模一个数的值 可以发现 $95041567 = 31*37*41*43*47$
那么我们可以对每个小质数取模 然后CRT
小质数取模 我们可以用做 直接矩阵快速幂
这个网上很多就不多说了
但是我们也可以用做

inline int calc(int n,int p,int *B){
  int c[N],d[N],a[N],cnt=0;
  while (n) a[cnt++]=n%p,n/=p;
  for (int i=0;i<=p;i++) c[i]=B[i]%p;
  for (int i=1;i<cnt;i++)
    for (int j=1;j<=a[i];j++){
      for (int k=0;k<p;k++)
    d[k]=(c[k+1]+i*c[k])%p;
      d[p]=(d[0]+d[1])%p;
      for (int k=0;k<=p;k++)
    c[k]=d[k];
    }
  return c[a[0]];    
}

这个方法在ACdreamer的blog里有写 但没有解释的很详细 我在被rxd教导后的理解就是每次迭代 相当于把下标加上了一个$p^m$这样最后取下标为 n <script type="math/tex" id="MathJax-Element-3">n</script>的就好了

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

const int P[]={31,37,41,43,47};

namespace CRT{
  inline int Pow(int a,int b,int p){
    a%=p;
    int ret=1; for (;b;b>>=1,a=a*a%p) if (b&1) ret=ret*a%p; return ret;
  }
  const int N=15;
  int n,P;
  int a[N],p[N];
  inline int Solve(){
    P=1; for (int i=1;i<=n;i++) P*=p[i];
    ll ret=0;
    for (int i=1;i<=n;i++){
      int c=P/p[i];
      c=(ll)c*Pow(c,p[i]-2,p[i])%P;
      (ret+=(ll)a[i]*c%P)%=P;
    }
    return ret;
  }
  inline void clear(){ n=0; }
  inline void add(int _p,int _a) { a[++n]=_a; p[n]=_p; }
}

const int N=65;

int T[N],B[5][N];
inline void Pre(int n,int p,int *B){
  T[1]=1; B[0]=B[1]=1;
  for (int i=2;i<=n;i++){
    T[i]=B[i-1];
    for (int j=i-1;j;j--) (T[j]+=T[j+1])%=p;
    B[i]=T[1];
  }
}

inline int calc(int n,int p,int *B){
  int c[N],d[N],a[N],cnt=0;
  while (n) a[cnt++]=n%p,n/=p;
  for (int i=0;i<=p;i++) c[i]=B[i]%p;
  for (int i=1;i<cnt;i++)
    for (int j=1;j<=a[i];j++){
      for (int k=0;k<p;k++)
    d[k]=(c[k+1]+i*c[k])%p;
      d[p]=(d[0]+d[1])%p;
      for (int k=0;k<=p;k++)
    c[k]=d[k];
    }
  return c[a[0]];    
}

int main(){
  int Q,n; 
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  for (int i=0;i<5;i++) Pre(50,P[i],B[i]);
  scanf("%d",&Q);
  while (Q--){
    scanf("%d",&n);
    for (int i=0;i<5;i++)
      CRT::add(P[i],calc(n,P[i],B[i]));
    printf("%d\n",CRT::Solve());
    CRT::clear();
  }
  return 0;
}

UPD：BZOJ 3501
怎么感觉我的预处理被卡常数了啊

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

const int PP=999999599;
const int P[]={2,13,5281,7283};

namespace CRT{
  inline int Pow(int a,int b,int p){
    a%=p;
    int ret=1; for (;b;b>>=1,a=a*a%p) if (b&1) ret=ret*a%p; return ret;
  }
  const int N=15;
  int n,P;
  int a[N],p[N];
  inline int Solve(){
    P=1; for (int i=1;i<=n;i++) P*=p[i];
    ll ret=0;
    for (int i=1;i<=n;i++){
      int c=P/p[i];
      c=(ll)c*Pow(c,p[i]-2,p[i])%P;
      (ret+=(ll)a[i]*c%P)%=P;
    }
    return ret;
  }
  inline void clear(){ n=0; }
  inline void add(int _p,int _a) { a[++n]=_a; p[n]=_p; }
}

const int N=10005;

int B[N],s[2][N];

inline void Pre(int n,int PP){
  B[0]=B[1]=s[0][0]=1,s[0][1]=2;
  for (int i=2,x=1;i<N;i++,x^=1){
    B[i]=s[x][0]=s[x^1][i-1];
   for (int j=1;j<=i;j++)
    s[x][j]=(s[x^1][j-1]+s[x][j-1])%PP;
   }
}

int c[N],d[N],a[N],cnt=0;

inline ll calc(ll n,int p){
  cnt=0;
  while (n) a[cnt++]=n%p,n/=p;
  for (int i=0;i<=p;i++) c[i]=B[i]%p;
  for (int i=1;i<cnt;i++)
    for (int j=1;j<=a[i];j++){
      for (int k=0;k<p;k++)
    d[k]=(c[k+1]+i*c[k])%p;
      d[p]=(d[0]+d[1])%p;
      for (int k=0;k<=p;k++)
    c[k]=d[k];
    }
  return c[a[0]];
}

inline ll Pow(ll a,ll b){
  ll ret=1;
  for (;b;b>>=1,a=a*a%PP)
    if (b&1)
      ret=ret*a%PP;
  return ret;
}

int main(){
  ll n,X,m; 
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  scanf("%lld%lld",&n,&m); Pre(7284,PP-1);
  for (int i=0;i<4;i++) CRT::add(P[i],calc(n,P[i]));
  X=CRT::Solve();
  m%=PP;
  printf("%lld\n",Pow(m,X));
  return 0;
}