OI中组合数的若干求法与CRT

ButterflyDew · 2018-10-25 12:56:46 · 算法·理论

OI中组合数的若干求法与CRT

只是下决心整理一下子呢~

说明：本篇文章采用\binom{a}{b}而不是C_{a}^b，以p指代模数，fac_i指代i!，inv_i指代i在\mod p下的逆元，invf_i指代i!在\mod p下的逆元。

一般性的组合数求法

计算式：

\binom{m}{n}=\frac{m!}{n!\times (m-n)!}

一、 \color{red} \text{杨辉三角法}

\binom{m}{n}=\binom{m-1}{n}+\binom{m-1}{n-1}

证明可以从计算式或者组合意义入手，不过多赘述了。

一般情况下，这个适用性并不是很广，但打起来还是很方便的，也不容易错。

值得一提的是，在有关容斥的一些证明中，经常用到这个式子消去一些项，比如证明\sum_{d|n}\mu(d)=[n=1]中

复杂度O(n^2)\sim O(1)，对p无要求

二、 \color{red} \text{预处理阶乘和阶乘逆元法}

根据计算式得到

\binom{m}{n}=fac_m\times invf_n \times invf_{m-n}

我们事先预处理阶乘及阶乘逆元，关于预处理阶乘逆元，有

invf_{i} \equiv invf_{i+1}\times (i+1) \pmod p

证明：

fac_i \times invf_i \equiv 1 \pmod p

fac_{i-1}\times i \times invf_i \equiv 1 \pmod p

只需要算一个就可以了

复杂度O(n+\log p)\sim O(1)，p为质数

这里p为质数才能保证一定有解，若不为质数，可能求不出来，因为不存在逆元。

麻烦一些的组合数求法

三、 \color{red} \tt{Lucas}\text{定理}

\binom{m}{n}=\binom{m \bmod p}{n \bmod p} \times \binom{m/p}{n/p}

证明（有兴趣的可以看一下）：

设n=ap+r_1,m=bp+r_2

取(1+x)^n\bmod p这样一个式子，我们对它进行变形（本质上对应将数进行p进制分解）

(1+x)^n

=(1+x)^{ap+r_1}

=((1+x)^p)^a\times(1+x)^{r_1}

引理：

①(a+b)^x=\sum_{i=0}^x\binom{x}{i}a^ib^{x-i}

②当p为质数时，任取x \in [1,p-1]且x \in N,满足\binom{p}{x}\equiv0 \pmod p。

则上式

≡(1+x^p)^a\times(1+x)^{r_1}\pmod p

=\sum_{i=0}^a \binom{a}{i}x^{ip}\times \sum_{j=0}^{r_1} \binom{r_1}{j}x^j \pmod p

整理得(1+x)^n=\sum_{i=0}^n\binom{n}{i}x^i≡\sum_{i=0}^a \binom{a}{i}x^{ip}\times \sum_{j=0}^{r1} \binom{r_1}{j}x^j \pmod p

对左边这样展开的这一项x^{bp+r_2}来说，它的系数为\binom{ap+r_1}{bq+r_2}

而右边当且仅当i=b,j=r_2时，可以构成这个项，而它的系数为\binom{a}{b}\times\binom{r_1}{r_2}

即\binom{m}{n}=\binom{m \bmod p}{n \bmod p} \times \binom{m/p}{n/p}则得证

用法上，这个要配合阶乘和阶乘逆元完成，当m小于p时，我们就要按照计算式直接算了。因此，\tt{lucas}定理适用于m,n较大的情况

这个定理其实用处不是很大，但是实现和形式上比较简单，所以可以按实际情况考虑使用。它还有一个优点，有时候可以避免n,m过大而强制使用龟速乘。

复杂度O(p+\log p)\sim O(\log n),p为质数

四、 \color{red} \tt{lucas}\text{定理}+\tt{CRT}\text{合并}

先说说前置的\tt{CRT}合并吧

\tt{CRT}合并

考虑两个同余方程

x \equiv b_1 \pmod {p_1}

x \equiv b_2 \pmod {p_2}

设x=kp_1+b_1,t=gcd(p_1,p_2)

将x代入方程2

kp_1+b_1\equiv b_2 \pmod {p_2}

此时解k

kp_1-yp_2=b_2-b_1

由裴蜀定理等，若k有解k_0(即t|b_2-b_1)，则k \equiv k_0 \pmod {\frac{p_1p_2}{t}}

则x \equiv k_0p_1+b_1 \pmod \frac{p_1p_2}{t}

其实这个就是\tt{EXCRT}了。我们对可能不互质的p，从第一个方程开始进行合并，只要有一个无解，整个同余方程组就无解了。

按理说\tt{EXCRT}是可以代替\tt{CRT}的所有功能的，但是\tt{CRT}打起来相对较简单，而且在扩展\tt{lucas}中用\tt{CRT}就够了，所以还是了解一下比较好。

中国剩余定理

对同余方程组

\left\{\begin{aligned}x \equiv a_1 \pmod {p_1}\\ x \equiv a_2 \pmod {p_2} \\ \dots \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\x \equiv a_n \pmod {p_n} \end{aligned}\right.

设P_0=\prod\limits_{i=1}^np_i，P_i=\frac{P_0}{p_i}

则方程组有通解x \equiv \sum\limits_{i=1}^nP_i\times inv_{P_i}\times a_1 \pmod {P_0}

理解起来不难，对P_i满足了可以整除其他的p_i，不会对其他方程产生影响，乘上自己的逆元和a_i，就满足了自己的方程。

两个方法的复杂度是一样的，均为O(n\log p)，但实际上普通的\tt{CRT}用的更多一些，注意有时候避免乘起来爆\text{long long}要用龟速乘

把前置技能点上了以后，就可以说一下方法4了

求\binom{m}{n} \bmod p

对p唯一分解p=\prod p_i^{c_i}

显然我们对每一个p_i^{c_i}得到一个答案之后，是可以用\tt{CRT}合并的。

但是发现模数不是质数的情况可能无解，因此实际上这种方法的局限性很大。

当题目已经给定模数而不是输入模数时，可能会用的到，比如古代猪文这个题。

为什么又强调\tt{lucas}定理呢？大概是因为分解以后模数不会太大的考虑。

这个方法的复杂度其实不好具体计算，但跑起来大概像这个O(\max p_i)\sim O(\log p)

五、 \color{red} \tt{exLucas}

其实这个东西挺麻烦的，能不写应该不会去写这玩意儿。

首先这个东西跟\tt{lucas}没关系，直接考虑组合数的计算式

\frac{m!}{(m-n)!n!} \bmod p

对p唯一分解p=\prod p_i^{c_i}

像四一样，考虑分解以后再合并

\frac{m!}{(m-n)!n!} \bmod p_i^{c_i}

因为p_i可能是m!的约数之类的，所以不妨先把分子分母的p_i都给拿下来。

首先，我们考虑m!里面有多少个p_i,设f(n)代表n!中有多少个，考虑一个递推式

f(n)=f(\lfloor\frac{n}{p_i}\rfloor)+\lfloor\frac{n}{p_i}\rfloor

理解起来很简单，先把所有p_i的倍数拿出来，都只产生1的贡献以后，这些p_i的倍数除上一个p_i后产生的序列就构成了子问题。

然后我们求剩下来的除掉所有p_i以后的阶乘

因为a\equiv a+p_i^{c_i}，所以我们把m!阶乘分成\lfloor\frac{m}{p_i^{c_i}}\rfloor段，然后暴力算一整段的贡献并用快速幂乘起来。

为了方便，我们依然不计算p_i的倍数，并把它当做一个子问题递归计算。

这个的求阶乘复杂度算起来其实是一个等比数列，是O(n)的

然后剩下的计算对求出来的东西求逆元并和剩下的p_i相乘，就得到了在模p_i^{c_i}意义下的答案，把所有的求出来以后直接\tt{CRT}合并即可。

发现这个东西的复杂度中有一堆并不重要的\log，比如\tt{CRT}，逆元什么的。真正复杂度瓶颈在于求阶乘的那个O(n)，可以近似的认为复杂度是O(\max p_i^{c_i})的。

贴一下我写的代码吧，写的不好看轻喷

   #include <cstdio>
   #define ll long long
   ll mod;
   void exgcd(ll &x,ll &y,ll a,ll b)
   {
       if(!b){x=1,y=0;return;}
       exgcd(x,y,b,a%b);
       ll tmp=x;
       x=y;
       y=tmp-a/b*y;
   }
   ll inv(ll a,ll p)
   {
       ll x,y;
       exgcd(x,y,a,p);
       return x;
   }
   ll quickpow(ll d,ll k,ll p)
   {
       ll f=1;
       while(k) {if(k&1)f=f*d%p;d=d*d%p;k>>=1;}
       return f;
   }
   ll fac(ll n,ll p,ll tp)
   {
       if(!n) return 1;
       ll f=1,res=1;
       for(ll i=1;i<tp;i++)
       {
           if(i%p) (f*=i)%=tp;
           if(i==n%tp) res=f;
       }
       f=quickpow(f,n/tp,tp);
       return fac(n/p,p,tp)*f%tp*res%tp;
   }
   ll C(ll m,ll n,ll p,ll tp)
   {
       ll ct=0;
       for(ll i=m;i;i/=p) ct+=i/p;
       for(ll i=n;i;i/=p) ct-=i/p;
       for(ll i=m-n;i;i/=p) ct-=i/p;
       return fac(m,p,tp)*inv(fac(n,p,tp),tp)%tp*inv(fac(m-n,p,tp),tp)%tp*quickpow(p,ct,tp)%tp;
   }
   ll CRT(ll m,ll n,ll p,ll tp)
   {
       return C(m,n,p,tp)*(mod/tp)%mod*inv(mod/tp,tp)%mod;
   }
   ll exlucas(ll m,ll n,ll p)
   {
       ll ans=0;
       for(ll i=2;i*i<=p;i++)
       {
           if(p%i==0)
           {
               ll tp=1;
               while(p%i==0) tp*=i,p/=i;
               (ans+=CRT(m,n,i,tp))%=mod;
           }
       }
       if(p!=1) (ans+=CRT(m,n,p,p))%=mod;
       return (ans%mod+mod)%mod;
   }
   int main()
   {
       ll n,m;
       scanf("%lld%lld%lld",&n,&m,&mod);
       printf("%lld\n",exlucas(n,m,mod));
       return 0;
   }

六、 \color{red} \text{一些灵活一点的应用}

给定一个长为n(n\le 3\times10^5)的有重序列，求这个序列在所有排列中的排名，对p(p\le 2\times 10^9)取模

其实就是有重复元素的排列问题。

像康托展开那样，讨论每一位的贡献。

不妨设当前处理到第i位的子问题了，令a_1,a_2,\dots,a_m代表从第i位到第n位的排列，b_1,b_2,\dots,b_m代表i到第n位的第一个排列，c_i代表数值为i的数字的出现次数。

若a_1=b_1，则当前位无贡献，处理下一位的子问题。

若a_1>b_1，则计算所有小于a_1的数值排在第一位后，后面的全排列，即为

\frac{m!}{\prod c_i} \times \sum_{d<a_1}c_d

这里c_i代表从第i为到之后的每一位的出现次数，你把后面的c_d乘下去就是后面的排列了。

后面的求和用树状数组维护，前面的删一个维护一个就可以了。

发现p_i^{c_i}可能很大，而n比较小，所以我们可以直接预处理到n的阶乘，当然也要先把p_i拿出来。

给出n(n\le 2\times 10^5)个质数，求这些质数乘积的约数之积模10^9+7

设乘积N=\prod\limits_{i=1}^kp_i^{c_i}，则约数个数d=\prod\limits_{i=1}^kc_i+1，那么答案为

\prod_{i=1}^kp_i^{\frac{d}{c_i+1}\sum\limits_{j=0}^{c_i}j} \bmod p

=\prod_{i=1}^kp_i^{\frac{d\times c_i}{2}}\bmod p

用扩展欧拉定理求一下指数，即\frac{d\times c_i}{2} \bmod \varphi(p)

把\varphi(p)拆成2\times 500000003做然后合并。对那个2随便从d中间拿下来一个就可以了，不需要多麻烦，因为只有一个。

有n(n\le 20)个集合，第i个集合有c_i(0\le c_i \le 10^{12})个相同元素，不同集合的元素互不相同，从中取出m(0\le m\le 10^{14})个元素，问有多少种取法，对10^9+7取膜

其实就是求多重集的组合数，设有多重集\{x_1\cdot c_1,x_2\cdot c_2,\dots,x_k\cdot c_k\},x是元素的值，而c是这个值的个数,设C=\sum\limits_{i=1}^k c_i。

首先我们不考虑每种元素的上限时（也就是当每种元素都有无限个的时候）的方案数

从组合的意义出发，相当于把n个元素有序的划分成k组，每组可以为空的方案数，即为

\frac{(n+k-1)!}{n! \times (k-1)!}=\binom{n+k-1}{k-1}

考虑x_i取了至少c_i+1个的方案数，首先我们把这c_i+1个取出来，然后再随便取,为\binom{n+k-c_i-2}{k-1}

则至少某一种元素取多了的方案数为\sum_{i=1}^{k}\binom{n+k-c_i-2}{k-1}，注意这里的方案数是有算重的。

类比可以得到至少两种元素取多了的方案数等等

套用容斥原理了，即为

\sum_{i=1}^{k}\binom{n+k-c_i-2}{k-1}-\sum_{i=1}^{k}\binom{n+k-c_i-c_j-3}{k-1}+\dots+(-1)^{k+1}\cdot\sum_{i=1}^{k}\binom{n-C-1}{k-1}

上式代表的是至少某一种元素取多了的方案数，注意这一步的“至少”和之前的意义并不一样。

则总方案为无限制的多重集的总排列数减去至少某一个元素多取了个数的方案数。

在计算这个题的组合数时，发现m,p都很大，而n很小，我们可以直接把计算式给约分，就是在算

\frac{m\times (m-1)\times \dots \times (m-n+1)}{n!}

处理一下底下的逆元直接爆算就可以了。

最后一点发现要用龟速乘，不过发现也可以先用\tt{lucas}定理就不需要龟速乘了。

题目地址

PER-Permutation

没找到。

CF451E Devu and Flowers