组合数学的实质还是DP,但是从通式角度处理的话有利于FFT等的实现。
首先推荐$Candy?$的球划分问题集合:
http://www.cnblogs.com/candy99/p/6400735.html
以下部分节选自
http://blog.****.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}}=\sum\limits_{i=0}^{k}(-1)^{k-i}s(k,i)\cdot n^i$,$n^{\overline{k}}=\sum\limits_{i=0}^{k}s(k,i)\cdot n^i$
也就是$s(n,k)=[x^k]\prod\limits_{i=0}^{n-1}(x+i)$
于是一种显然的做法就是分治NTT,复杂度$O(n\log^2 n)$。
#include<cstdio>
#include<algorithm>
#define rep(i,l,r) for (int i=(l); i<=(r); i++)
using namespace std; const int N=,mod=;
int n,A,B,top,stk[],rev[N],fac[N],inv[N],S[N],tmp[][N]; int ksm(int a,int b){
int res=;
for (; b; a=1ll*a*a%mod,b>>=)
if (b & ) res=1ll*res*a%mod;
return res;
} void NTT(int a[],int n,int f){
for (int i=; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
for (int i=; i<n; i<<=){
int wn=ksm(,f ? (mod-)/(i<<) : (mod-)-(mod-)/(i<<));
for (int p=i<<,j=; j<n; j+=p){
int w=;
for (int k=; k<i; k++,w=1ll*w*wn%mod){
int x=a[j+k],y=1ll*w*a[i+j+k]%mod;
a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod;
}
}
}
if (f) return;
int inv=ksm(n,mod-);
for (int i=; i<n; i++) a[i]=1ll*a[i]*inv%mod;
} int solve(int l,int r,int a[]){
if (l==r){ a[]=l; a[]=; return ; }
int ls=stk[top--],rs=stk[top--],mid=(l+r)>>;
int l1=solve(l,mid,tmp[ls]),l2=solve(mid+,r,tmp[rs]),n=,L=;
for (; n<=l1+l2; n<<=) L++;
for (int i=; i<n; i++) rev[i]=(rev[i>>]>>)|((i&)<<(L-));
NTT(tmp[ls],n,); NTT(tmp[rs],n,);
for (int i=; i<n; i++) a[i]=1ll*tmp[ls][i]*tmp[rs][i]%mod;
NTT(a,n,); stk[++top]=ls; stk[++top]=rs;
for (int i=; i<n; i++) tmp[ls][i]=tmp[rs][i]=;
return l1+l2;
} int main(){
freopen("960g.in","r",stdin);
freopen("960g.out","w",stdout);
scanf("%d%d%d",&n,&A,&B);
if (n==){ if (A== && B==) puts(""); else puts(""); return ; }
rep(i,,) stk[i]=i; top=;
solve(,n-,S);
fac[]=; rep(i,,A+B) fac[i]=1ll*fac[i-]*i%mod;
inv[A+B]=ksm(fac[A+B],mod-);
for (int i=A+B-; ~i; i--) inv[i]=1ll*inv[i+]*(i+)%mod;
printf("%lld\n",1ll*S[A+B-]*fac[A+B-]%mod*inv[A-]%mod*inv[B-]%mod);
return ;
}
O(nlog^2n)
另一种做法是:令$F_n(x)=\prod\limits_{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$
$\begin{aligned} F_n(x+n)&=\sum_{i=0}^{n-1}a_i (x+n)^i\\ &=\sum_{i=0}^{n-1}a_i\sum_{j=0}^i{i\choose j}n^{i-j}x^j\\&=\sum_{i=0}^{n-1}(\sum_{j=i}^n {j\choose i}n^{j-i}a_j)x^i \end{aligned}$
括号内的部分拆开之后,可以分成差相等的两个部分,意味着可以翻转之后卷积。
那么每次递归前一半,后面一半卷积求解即可得到答案。
时间复杂度为一个log。(来自这里)
#include<cstdio>
#include<algorithm>
#define rep(i,l,r) for (int i=(l); i<=(r); i++)
using namespace std; const int N=,mod=;
int n,A,B,top,rev[N],fac[N],inv[N],pw[N],S[N],a[N],b[N]; int ksm(int a,int b){
int res=;
for (; b; a=1ll*a*a%mod,b>>=)
if (b & ) res=1ll*res*a%mod;
return res;
} void NTT(int a[],int n,int f){
for (int i=; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
for (int i=; i<n; i<<=){
int wn=ksm(,f ? (mod-)/(i<<) : (mod-)-(mod-)/(i<<));
for (int p=i<<,j=; j<n; j+=p){
int w=;
for (int k=; k<i; k++,w=1ll*w*wn%mod){
int x=a[j+k],y=1ll*w*a[i+j+k]%mod;
a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod;
}
}
}
if (f) return;
int inv=ksm(n,mod-);
for (int i=; i<n; i++) a[i]=1ll*a[i]*inv%mod;
} void solve(int l){
if (l<=){ S[l]=; return; }
if (l&){
solve(l-);
for (int i=l; i; i--) S[i]=(S[i-]+1ll*S[i]*(l-))%mod;
}else{
int m=l>>,n=,L=; solve(m);
for (; n<=l; n<<=) L++;
pw[]=; rep(i,,m) pw[i]=1ll*pw[i-]*m%mod;
rep(i,,m) a[i]=1ll*S[i]*fac[i]%mod;
rep(i,,m) b[i]=1ll*pw[m-i]*inv[m-i]%mod;
for (int i=; i<n; i++) rev[i]=(rev[i>>]>>)|((i&)<<(L-));
NTT(a,n,); NTT(b,n,);
for (int i=; i<n; i++) a[i]=1ll*a[i]*b[i]%mod;
NTT(a,n,);
rep(i,,m) a[i]=1ll*a[i+m]*inv[i]%mod,b[i]=S[i];
rep(i,m+,n) a[i]=b[i]=;
NTT(a,n,); NTT(b,n,);
for (int i=; i<n; i++) a[i]=1ll*a[i]*b[i]%mod;
NTT(a,n,);
rep(i,,l) S[i]=a[i];
for (int i=; i<n; i++) a[i]=b[i]=;
}
} int main(){
freopen("960g.in","r",stdin);
freopen("960g.out","w",stdout);
scanf("%d%d%d",&n,&A,&B); int m=max(A+B,n);
fac[]=; rep(i,,m) fac[i]=1ll*fac[i-]*i%mod;
inv[m]=ksm(fac[m],mod-);
for (int i=m-; ~i; i--) inv[i]=1ll*inv[i+]*(i+)%mod;
solve(n-); int ans=1ll*fac[A+B-]*inv[A-]%mod*inv[B-]%mod*S[A+B-]%mod;
printf("%d\n",ans);
return ;
}
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)=\sum\limits_{k=0}^{i} (-1)^k \binom{i}{k}(i-k)^n$
显然就有$$S(n,i)=\sum\limits_{k=0}^{i} (-1)^k\binom{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}^n\binom{n}{k}B_k$
Catalan数
$$C_n=\frac{1}{n+1}\binom{2n}{n}=\frac{(2n)!}{(n+1)!n!}=\binom{2n}{n}-\binom{2n}{n+1}=\frac{1}{n+1} \sum\limits_{i=0}^{n}\binom{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=1\quad\quad\sum_{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$$
自然数幂和