• Stirling数,Bell数,Catalan数,Bernoulli数


    组合数学的实质还是DP,但是从通式角度处理的话有利于FFT等的实现。

    首先推荐$Candy?$的球划分问题集合:

    http://www.cnblogs.com/candy99/p/6400735.html

    以下部分节选自

    http://blog.csdn.net/sr_19930829/article/details/40888349

    第一类Stirling数

      定理:第一类Stirling数$s(p,k)$计数的是把p个对象排成k个非空循环排列的方法数。

           证明:把上述定理叙述中的循环排列叫做圆圈。递推公式为:

           $s(p,p)=1 (p>=0)$    有p个人和p个圆圈,每个圆圈就只有一个人

           $s(p,0)=0 (p>=1)$    如果至少有1个人,那么任何的安排都至少包含一个圆圈$$s(p,k)=(p-1)*s(p-1,k)+s(p-1,k-1)$$       设人被标上$1,2,ldots,p$。将这p个人排成k个圆圈有两种情况。第一种排法是在一个圆圈里只有标号为$p$的人自己,排法有$s(p-1,k-1)$个。第二种排法中,$p$至少和另一个人在一个圆圈里。这些排法可以通过把$1,2,ldots,p-1$排成k个圆圈再把p放在$1,2,ldots,p-1$任何一人的左边得到,因此第二种类型的排法共有$(p-1)*s(p-1,k)$种排法。

           在证明中我们所做的就是把${1,2,ldots,p}$划分到k个非空且不可区分的盒子,然后将每个盒子中的元素排成一个循环排列。

       组合数通式:$$s(k,i)=frac{C_n^k k!}{sum_{i=0}^{k}(-1)^{i+k}n^i}$$

       第一类斯特林数快速求法:

        首先我们有:$n^{underline{k}}=sumlimits_{i=0}^{k}(-1)^{k-i}s(k,i)cdot n^i$,$n^{overline{k}}=sumlimits_{i=0}^{k}s(k,i)cdot n^i$

        也就是$s(n,k)=[x^k]prodlimits_{i=0}^{n-1}(x+i)$

        于是一种显然的做法就是分治NTT,复杂度$O(nlog^2 n)$。

     1 #include<cstdio>
     2 #include<algorithm>
     3 #define rep(i,l,r) for (int i=(l); i<=(r); i++)
     4 using namespace std;
     5 
     6 const int N=300010,mod=998244353;
     7 int n,A,B,top,stk[52],rev[N],fac[N],inv[N],S[N],tmp[52][N];
     8 
     9 int ksm(int a,int b){
    10     int res=1;
    11     for (; b; a=1ll*a*a%mod,b>>=1)
    12         if (b & 1) res=1ll*res*a%mod;
    13     return res;
    14 }
    15 
    16 void NTT(int a[],int n,int f){
    17     for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
    18     for (int i=1; i<n; i<<=1){
    19         int wn=ksm(3,f ? (mod-1)/(i<<1) : (mod-1)-(mod-1)/(i<<1));
    20         for (int p=i<<1,j=0; j<n; j+=p){
    21             int w=1;
    22             for (int k=0; k<i; k++,w=1ll*w*wn%mod){
    23                 int x=a[j+k],y=1ll*w*a[i+j+k]%mod;
    24                 a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod;
    25             }
    26         }
    27     }
    28     if (f) return;
    29     int inv=ksm(n,mod-2);
    30     for (int i=0; i<n; i++) a[i]=1ll*a[i]*inv%mod;
    31 }
    32 
    33 int solve(int l,int r,int a[]){
    34     if (l==r){ a[0]=l; a[1]=1; return 1; }
    35     int ls=stk[top--],rs=stk[top--],mid=(l+r)>>1;
    36     int l1=solve(l,mid,tmp[ls]),l2=solve(mid+1,r,tmp[rs]),n=1,L=0;
    37     for (; n<=l1+l2; n<<=1) L++;
    38     for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    39     NTT(tmp[ls],n,1); NTT(tmp[rs],n,1);
    40     for (int i=0; i<n; i++) a[i]=1ll*tmp[ls][i]*tmp[rs][i]%mod;
    41     NTT(a,n,0); stk[++top]=ls; stk[++top]=rs;
    42     for (int i=0; i<n; i++) tmp[ls][i]=tmp[rs][i]=0;
    43     return l1+l2;
    44 }
    45 
    46 int main(){
    47     freopen("960g.in","r",stdin);
    48     freopen("960g.out","w",stdout);
    49     scanf("%d%d%d",&n,&A,&B);
    50     if (n==1){ if (A==1 && B==1) puts("1"); else puts("0"); return 0; }
    51     rep(i,1,50) stk[i]=i; top=50;
    52     solve(0,n-2,S);
    53     fac[0]=1; rep(i,1,A+B) fac[i]=1ll*fac[i-1]*i%mod;
    54     inv[A+B]=ksm(fac[A+B],mod-2);
    55     for (int i=A+B-1; ~i; i--) inv[i]=1ll*inv[i+1]*(i+1)%mod;
    56     printf("%lld
    ",1ll*S[A+B-2]*fac[A+B-2]%mod*inv[A-1]%mod*inv[B-1]%mod);
    57     return 0;
    58 }
    O(nlog^2n)

        另一种做法是:令$F_n(x)=prodlimits_{i=0}^{n-1}(x+i)$,则有$F_{2n}(x)=F_n(x)F_{n}(x+n)$。
        $F_n(x)$我们递归求得其答案,现在考虑如何利用$F_n(x)$快速求出$F_n(x+n)$。
        在这里我们假设$F_n(x)=sum_{i=0}^{n-1} a_i x^i$
           $egin{aligned} F_n(x+n)&=sum_{i=0}^{n-1}a_i (x+n)^i\ &=sum_{i=0}^{n-1}a_isum_{j=0}^i{ichoose j}n^{i-j}x^j\&=sum_{i=0}^{n-1}(sum_{j=i}^n {jchoose i}n^{j-i}a_j)x^i end{aligned}$
        括号内的部分拆开之后,可以分成差相等的两个部分,意味着可以翻转之后卷积。
        那么每次递归前一半,后面一半卷积求解即可得到答案。
        时间复杂度为一个log。(来自这里

     1 #include<cstdio>
     2 #include<algorithm>
     3 #define rep(i,l,r) for (int i=(l); i<=(r); i++)
     4 using namespace std;
     5 
     6 const int N=300010,mod=998244353;
     7 int n,A,B,top,rev[N],fac[N],inv[N],pw[N],S[N],a[N],b[N];
     8 
     9 int ksm(int a,int b){
    10     int res=1;
    11     for (; b; a=1ll*a*a%mod,b>>=1)
    12         if (b & 1) res=1ll*res*a%mod;
    13     return res;
    14 }
    15 
    16 void NTT(int a[],int n,int f){
    17     for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
    18     for (int i=1; i<n; i<<=1){
    19         int wn=ksm(3,f ? (mod-1)/(i<<1) : (mod-1)-(mod-1)/(i<<1));
    20         for (int p=i<<1,j=0; j<n; j+=p){
    21             int w=1;
    22             for (int k=0; k<i; k++,w=1ll*w*wn%mod){
    23                 int x=a[j+k],y=1ll*w*a[i+j+k]%mod;
    24                 a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod;
    25             }
    26         }
    27     }
    28     if (f) return;
    29     int inv=ksm(n,mod-2);
    30     for (int i=0; i<n; i++) a[i]=1ll*a[i]*inv%mod;
    31 }
    32 
    33 void solve(int l){
    34     if (l<=1){ S[l]=1; return; }
    35     if (l&1){
    36         solve(l-1);
    37         for (int i=l; i; i--) S[i]=(S[i-1]+1ll*S[i]*(l-1))%mod;
    38     }else{
    39         int m=l>>1,n=1,L=0; solve(m);
    40         for (; n<=l; n<<=1) L++;
    41         pw[0]=1; rep(i,1,m) pw[i]=1ll*pw[i-1]*m%mod;
    42         rep(i,0,m) a[i]=1ll*S[i]*fac[i]%mod;
    43         rep(i,0,m) b[i]=1ll*pw[m-i]*inv[m-i]%mod;
    44         for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    45         NTT(a,n,1); NTT(b,n,1);
    46         for (int i=0; i<n; i++) a[i]=1ll*a[i]*b[i]%mod;
    47         NTT(a,n,0);
    48         rep(i,0,m) a[i]=1ll*a[i+m]*inv[i]%mod,b[i]=S[i];
    49         rep(i,m+1,n) a[i]=b[i]=0;
    50         NTT(a,n,1); NTT(b,n,1);
    51         for (int i=0; i<n; i++) a[i]=1ll*a[i]*b[i]%mod;
    52         NTT(a,n,0);
    53         rep(i,0,l) S[i]=a[i];
    54         for (int i=0; i<n; i++) a[i]=b[i]=0;
    55     }
    56 }
    57 
    58 int main(){
    59     freopen("960g.in","r",stdin);
    60     freopen("960g.out","w",stdout);
    61     scanf("%d%d%d",&n,&A,&B); int m=max(A+B,n);
    62     fac[0]=1; rep(i,1,m) fac[i]=1ll*fac[i-1]*i%mod;
    63     inv[m]=ksm(fac[m],mod-2);
    64     for (int i=m-1; ~i; i--) inv[i]=1ll*inv[i+1]*(i+1)%mod;
    65     solve(n-1); int ans=1ll*fac[A+B-2]*inv[A-1]%mod*inv[B-1]%mod*S[A+B-2]%mod;
    66     printf("%d
    ",ans);
    67     return 0;
    68 }
    O(nlogn)

    第二类Stirling数

       定理:第二类Stirling数$S(p,k)$计数的是把p元素集合划分到k个不可区分的盒子里且没有空盒子的划分个数。

            证明:元素在拿些盒子并不重要,唯一重要的是各个盒子里装的是什么,而不管哪个盒子装了什么。

            递推公式有:$$S(p,p)=1 (p>=0) $$ $$ S(p,0)=0 (p>=1) $$ $$S(p,k)=k*S(p-1,k)+S(p-1,k-1)$$考虑将前p个正整数,$1,2,ldots,p$的集合作为要被划分的集合,把${1,2,ldots,p}$分到k个非空且不可区分的盒子的划分有两种情况:

           (1)那些使得p自己单独在一个盒子的划分,存在有$S(p-1,k-1)$种划分个数

           (2)那些使得p不单独自己在一个盒子的划分,存在有 $k*S(p-1,k)$种划分个数

            考虑第二种情况,p不单独自己在一个盒子,也就是p和其他元素在一个集合里面,也就是说在没有放p之前,有p-1个元素已经分到了k个非空且不可区分的盒子里面(划分个数为$S(p-1,k)$,那么现在问题是把p放在哪个盒子里面呢,有k种选择,所以存在有$k*S(p-1,k)$。

      第二类斯特林数组合数通式求法:

        先考虑每个盒子都可以区分的方案数:$S'(n,i)=S(n,i)*i!$

        首先,原问题不允许有空盒子,通过容斥解除这个限制:$S'(n,i)=sumlimits_{k=0}^{i} (-1)^k inom{i}{k}(i-k)^n$

        显然就有$$S(n,i)=sumlimits_{k=0}^{i} (-1)^kinom{i}{k}(i-k)^n *frac{1}{i!}$$适合FFT的写法:$$S(n,m)=sum_{i=0}^{m}frac{(-1)^i}{i!}frac{(m-i)^n}{(m-i)!}$$

        应用:http://www.cnblogs.com/candy99/p/6648754.html

    Bell数

            定理:Bell数$B(p)$是将p元素集合分到非空且不可区分盒子的划分个数(没有说分到几个盒子里面)。$$B(p)=S(p,0)+S(p,1)+.....+S(p,k)$$

            所以要求Bell数就要先求出第二类Stiring数。

      $B_0=1$     $B_{n+1}=sum_{k=0}^ninom{n}{k}B_k$

    Catalan数

    $$C_n=frac{1}{n+1}inom{2n}{n}=frac{(2n)!}{(n+1)!n!}=inom{2n}{n}-inom{2n}{n+1}=frac{1}{n+1} sumlimits_{i=0}^{n}inom{n}{i}^2$$$$C_n=frac{4n-2}{n+1}C_{n-1}$$$$C_n=sum_{i=0}^{n-1}C_i C_{n-i-1}quad C_{n+1}=sum_{i=0}^{n}C_i C{n-i}$$

    Bernoulli数

    $$B_0=1quadquadsum_{k=0}^{n}C_{n+1}^k B_k =0$$$$B_n=-frac{1}{n+1}(C_{n+1}^0 B_0+C_{n+1}^1 B_1+ cdots +C_{n+1}^{n-1}B_{n-1})$$$$sum_{i=1}^{n} i^k=frac{1}{k+1} sum_{i=1}^{k+1} C_{k+1}^{i}B_{k+1-i} (n+1)^i$$

    自然数幂和

      http://blog.csdn.net/doyouseeman/article/details/50826293

      http://blog.csdn.net/a_crazy_czy/article/details/50926374

  • 相关阅读:
    (转)simple-framework(MaliSDK框架分析)
    (转)libhybris及EGL Platform-在Glibc生态中重用Android的驱动
    (原)在firefly_rk3288开发板上解决openGL在设置32位色深以后出现花屏的问题
    参考论坛:Mali kernel driver TX011-SW-99002-r5p1-00rel0 for firefly
    (转)android媒体--stagefright概述【一】
    (转)android系统架构及源码目录结构
    (原)linux下利用cmake来编译jthread开源库
    (原)U盘可见容量不能被识别的处理方法
    学会阅读Java字节码
    JVM
  • 原文地址:https://www.cnblogs.com/HocRiser/p/8420867.html
Copyright © 2020-2023  润新知