zoukankan      html  css  js  c++  java
  • OI中组合数的若干求法与CRT

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

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


    说明:本篇文章采用(inom{a}{b})而不是(C_{a}^b),以(p)指代模数,(fac_i)指代(i!)(inv_i)指代(i)(mod p)下的逆元,(invf_i)指代(i!)(mod p)下的逆元。

    一般性的组合数求法

    计算式:

    $$inom{m}{n}=frac{m!}{n! imes (m-n)!}$$

    一、 杨辉三角法

    [inom{m}{n}=inom{m-1}{n}+inom{m-1}{n-1} ]

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

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

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

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

    二、 预处理阶乘和阶乘逆元法

    根据计算式得到

    [inom{m}{n}=fac_m imes invf_n imes invf_{m-n} ]

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

    [invf_{i} equiv invf_{i+1} imes (i+1) pmod p ]

    证明:

    [fac_i imes invf_i equiv 1 pmod p ]

    [fac_{i-1} imes i imes invf_i equiv 1 pmod p ]

    只需要算一个就可以了

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

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

    麻烦一些的组合数求法

    三、 Lucas定理

    [inom{m}{n}=inom{m mod p}{n mod p} imes inom{m/p}{n/p} ]

    证明(有兴趣的可以看一下):

    (n=ap+r_1,m=bp+r_2)

    ((1+x)^nmod p)这样一个式子,我们对它进行变形(本质上对应将数进行p进制分解)

    [(1+x)^n ]

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

    [=((1+x)^p)^a imes(1+x)^{r_1} ]

    引理:

    ((a+b)^x=sum_{i=0}^xinom{x}{i}a^ib^{x-i})

    ②当(p)为质数时,任取(x in [1,p-1])(x in N),满足(inom{p}{x}equiv0 pmod p)

    则上式

    [≡(1+x^p)^a imes(1+x)^{r_1}pmod p ]

    [=sum_{i=0}^a inom{a}{i}x^{ip} imes sum_{j=0}^{r_1} inom{r_1}{j}x^j pmod p ]

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

    对左边这样展开的这一项(x^{bp+r_2})来说,它的系数为(inom{ap+r_1}{bq+r_2})

    而右边当且仅当(i=b,j=r_2)时,可以构成这个项,而它的系数为(inom{a}{b} imesinom{r_1}{r_2})

    即$$inom{m}{n}=inom{m mod p}{n mod p} imes inom{m/p}{n/p}$$则得证

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

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

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

    四、 lucas+CRT合并

    先说说前置的( t{CRT})合并吧

    ( t{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_1equiv 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}}$$

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

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

    中国剩余定理

    对同余方程组

    [left{egin{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} ight. ]

    (P_0=prodlimits_{i=1}^np_i)(P_i=frac{P_0}{p_i})

    则方程组有通解$$x equiv sumlimits_{i=1}^nP_i imes inv_{P_i} imes a_1 pmod {P_0}$$

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

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

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

    求$$inom{m}{n} mod p$$

    (p)唯一分解(p=prod p_i^{c_i})

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

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

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

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

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

    五、 exLucas

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

    首先这个东西跟( t{lucas})没关系,直接考虑组合数的计算式

    [frac{m!}{(m-n)!n!} mod p ]

    (p)唯一分解(p=prod p_i^{c_i})

    一样,考虑分解以后再合并

    [frac{m!}{(m-n)!n!} mod p_i^{c_i} ]

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

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

    [f(n)=f(lfloorfrac{n}{p_i} floor)+lfloorfrac{n}{p_i} floor ]

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

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

    因为(aequiv a+p_i^{c_i}),所以我们把(m!)阶乘分成(lfloorfrac{m}{p_i^{c_i}} floor)段,然后暴力算一整段的贡献并用快速幂乘起来。

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

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

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

    发现这个东西的复杂度中有一堆并不重要的(log),比如( t{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
    ",exlucas(n,m,mod));
        return 0;
    }
    

    六、 一些灵活一点的应用

    • 给定一个长为(n(nle 3 imes10^5))的有重序列,求这个序列在所有排列中的排名,对(p(ple 2 imes 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} imes sum_{d<a_1}c_d ]

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

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

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

    • 给出(n(nle 2 imes 10^5))个质数,求这些质数乘积的约数之积模(10^9+7)

    设乘积(N=prodlimits_{i=1}^kp_i^{c_i}),则约数个数(d=prodlimits_{i=1}^kc_i+1),那么答案为

    [prod_{i=1}^kp_i^{frac{d}{c_i+1}sumlimits_{j=0}^{c_i}j} mod p ]

    [=prod_{i=1}^kp_i^{frac{d imes c_i}{2}}mod p ]

    用扩展欧拉定理求一下指数,即(frac{d imes c_i}{2} mod varphi(p))

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

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

    其实就是求多重集的组合数,设有多重集({x_1cdot c_1,x_2cdot c_2,dots,x_kcdot c_k}),(x)是元素的值,而(c)是这个值的个数,设(C=sumlimits_{i=1}^k c_i)

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

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

    [frac{(n+k-1)!}{n! imes (k-1)!}=inom{n+k-1}{k-1} ]

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

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

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

    套用容斥原理了,即为

    [sum_{i=1}^{k}inom{n+k-c_i-2}{k-1}-sum_{i=1}^{k}inom{n+k-c_i-c_j-3}{k-1}+dots+(-1)^{k+1}cdotsum_{i=1}^{k}inom{n-C-1}{k-1} ]

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

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

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

    [frac{m imes (m-1) imes dots imes (m-n+1)}{n!} ]

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

    最后一点发现要用龟速乘,不过发现也可以先用( t{lucas})定理就不需要龟速乘了。

    题目地址

    PER-Permutation

    没找到。

    CF451E Devu and Flowers

  • 相关阅读:
    doges
    Unity Fps示例
    使用Unity的2D功能开发弹球游戏
    Unity UGUI 原理篇(二):Canvas Scaler 縮放核心
    UGUI 深度優化提升手遊效能
    关于Unity中的UGUI优化,你可能遇到这些问题
    git branch --set-upstream 本地关联远程分支
    git rm 与 git reset
    Git笔记之初识vi编辑器
    [内容分享]粗略判断Shader每条代码的成本
  • 原文地址:https://www.cnblogs.com/butterflydew/p/9849823.html
Copyright © 2011-2022 走看看