nowcoder50596 Combination （含组合数取模讲解）

题目链接

点击跳转

从阶乘说起

当$p=10007$的时候，如何计算$n!\mathrm{m}\mathrm{o}\mathrm{d}p$？

显然，当$n\ge p$时，$n!=0$，否则我可以线性求出$n!$

注意到，$\left(\genfrac{}{}{0ex}{}{n}{m}\right)=\frac{n!}{m!(n-m)!}$，如何计算$\frac{n!}{m!(n-m)!}$呢？

如果我算出$n!\mathrm{m}\mathrm{o}\mathrm{d}p$，再算出$m!\mathrm{m}\mathrm{o}\mathrm{d}p$，再算出$(n-m)\mathrm{m}\mathrm{o}\mathrm{d}p$，最后求一下分母的逆元，看起来貌似是很常用的流程，但是当$n\ge p$的时候，算出来答案是$0$。

这是为啥呢，原因很简单，$x$在模$p$意义下逆元存在当且仅当$(x,p)=1$，而$n!$是$p$的倍数，所以这里求逆元的算法失效了。

那怎么办呢？

我们可以计算出分子、分母中$p$的次数，然后用分子的减去分母的，就得到最后答案中$p$的次数，因为组合数都是整数，所以这个次数不能为负。如果这个次数为正，那么$\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}p=0$，否则，答案显然不为$0$，我们有别的方法来计算这个东西。

只要把分子、分母中出去$p$之外的乘积在模$p$意义下的值，以及$p$的幂次都求出来，我们就具备了求解组合数在模$p$意义下的值的所有信息。

现在要解决的问题是，如何将$n!$表示成$x\times p^q$，其中$(x,p)=1$

因为$p$是素数，所以$1,2,\dots ,p-1$中都不包含$p$，$p$中包含$p$的一次幂。接下来是$p+1,p+2,\dots ,2p-1$这一段也不包含$p$，而且这一段的乘积在模$p$意义下和第一段相等，接下来一个数字是$2p$…以此类推，可以发现$n!$被分成$\lfloor \frac{n}{p}\rfloor$个整段，以及最后一段$(n\mathrm{m}\mathrm{o}\mathrm{d}p)$个数字的连乘，前面每段都是$(p-1)!$和一个$p$的倍数。

也就是
$$n!=(p-1){!}^{\lfloor \frac{n}{p}\rfloor }\times (n\%p)!\times p^{\lfloor \frac{n}{p}\rfloor }\times (\lfloor \frac{n}{p}\rfloor )!$$

其中$\%$表示求余数

卢卡斯定理

我已经得到了上面所描述的重要等式，将其带入$\left(\genfrac{}{}{0ex}{}{n}{m}\right)=\frac{n!}{m!(n-m)!}$并化简，得到：

$$\left(\genfrac{}{}{0ex}{}{n}{m}\right)=p{!}^{\lfloor n/p\rfloor -\lfloor m/p\rfloor -\lfloor (n-m)/p\rfloor }\times \frac{(n\%p)!}{(m\%p)!((n-m)\%p)!}\times \frac{\lfloor n/p\rfloor !}{\lfloor m/p\rfloor !\lfloor (n-m)/p\rfloor !}$$

如果$\lfloor n/p\rfloor -\lfloor m/p\rfloor -\lfloor (n-m)/p\rfloor =0$，那么$\left(\genfrac{}{}{0ex}{}{n}{m}\right)\%p=0$，否则$\left(\genfrac{}{}{0ex}{}{n}{m}\right)=\left(\genfrac{}{}{0ex}{}{\lfloor n/p\rfloor }{\lfloor m/p\rfloor }\right)\left(\genfrac{}{}{0ex}{}{n\%p}{m\%p}\right)$

$\lfloor n/p\rfloor -\lfloor m/p\rfloor -\lfloor (n-m)/p\rfloor >0$，仅当$n\%p<m\%p$，如果我规定“当$m>n$时$\left(\genfrac{}{}{0ex}{}{n}{m}\right)=0$”，那么两种情况就可以统一写成

$$\left(\genfrac{}{}{0ex}{}{n}{m}\right)=\left(\genfrac{}{}{0ex}{}{\lfloor n/p\rfloor }{\lfloor m/p\rfloor }\right)\left(\genfrac{}{}{0ex}{}{n\%p}{m\%p}\right)$$

其中$\%$表示取余数

这就是卢卡斯定理了

💡ひらめき：

从上述讨论，可以看出只有$m$的$p$进制表示的每一位都不大于$n$的$p$进制表示的对应位的时候，$\left(\genfrac{}{}{0ex}{}{n}{m}\right)$才不等于$0$

那么判断$\left(\genfrac{}{}{0ex}{}{n}{m}\right)$为奇数的充要条件也就成了$n\&m=m$，其中$\&$是按位与

对于更一般的情况

当$p$不是素数时，就存在一个问题，上述等式虽然仍成立，但是我们在递归的终点却没法用逆元去计算$\left(\genfrac{}{}{0ex}{}{n}{m}\right)$，那咋办呢？

这样，我把$p$分解成若干个素数的幂，做完了之后中国剩余定理合并一下就可以了

现在的问题变成，如何计算$\left(\genfrac{}{}{0ex}{}{n}{m}\right)\%t^c$？，其中$t$是素数

利用这个等式：$n!=(p-1){!}^{\lfloor \frac{n}{p}\rfloor }\times (n\%p)!\times p^{\lfloor \frac{n}{p}\rfloor }\times (\lfloor \frac{n}{p}\rfloor )!$，便可以递归计算出$n!,m!,(n-m)!$中$t$的幂次，以及和$t$互质的部分在模$t^c$意义下的值，最后合并一下就好了。

代码

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
#define iinf 0x3f3f3f3f
#define linf (1ll<<60)
#define eps 1e-8
#define maxn 1000010
#define maxe 1000010
#define cl(x) memset(x,0,sizeof(x))
#define rep(_,__) for(_=1;_<=(__);_++)
#define em(x) emplace(x)
#define emb(x) emplace_back(x)
#define emf(x) emplace_front(x)
#define fi first
#define se second
#define de(x) cerr<<#x<<" = "<<x<<endl
using namespace std;
using namespace __gnu_pbds;
typedef long long ll;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
ll read(ll x=0)
{
    ll c, f(1);
    for(c=getchar();!isdigit(c);c=getchar())if(c=='-')f=-f;
    for(;isdigit(c);c=getchar())x=x*10+c-0x30;
    return f*x;
}
struct EasyMath
{
    ll prime[maxn], phi[maxn], mu[maxn];
    bool mark[maxn];
    ll fastpow(ll a, ll b, ll c)
    {
        ll t(a%c), ans(1ll);
        for(;b;b>>=1,t=t*t%c)if(b&1)ans=ans*t%c;
        return ans;
    }
    void exgcd(ll a, ll b, ll &x, ll &y)
    {
        if(!b){x=1,y=0;return;}
        ll xx, yy;
        exgcd(b,a%b,xx,yy);
        x=yy, y=xx-a/b*yy;
    }
    ll inv(ll x, ll p)  //p是素数
    {return fastpow(x%p,p-2,p);}
    ll inv2(ll a, ll p)
    {
        ll x, y;
        exgcd(a,p,x,y);
        return (x+p)%p;
    }
    void shai(ll N)
    {
        ll i, j;
        for(i=2;i<=N;i++)mark[i]=false;
        *prime=0;
        phi[1]=mu[1]=1;
        for(i=2;i<=N;i++)
        {
            if(!mark[i])prime[++*prime]=i, mu[i]=-1, phi[i]=i-1;
            for(j=1;j<=*prime and i*prime[j]<=N;j++)
            {
                mark[i*prime[j]]=true;
                if(i%prime[j]==0)
                {
                    phi[i*prime[j]]=phi[i]*prime[j];
                    break;
                }
                mu[i*prime[j]]=-mu[i];
                phi[i*prime[j]]=phi[i]*(prime[j]-1);
            }
        }
    }
    ll CRT(vector<ll> a, vector<ll> m) //要求模数两两互质
    {
        ll M=1, ans=0, n=a.size(), i;
        for(i=0;i<n;i++)M*=m[i];
        for(i=0;i<n;i++)(ans+=a[i]*(M/m[i])%M*inv2(M/m[i],m[i]))%=M;
        return ans;
    }
}em;
struct CombinatorialNumber_mod
{
    ll p[20], q[20], fact[20][maxn], tot, t[20];
    void init(ll P)
    {
        tot=0;
        ll i, j;
        for(i=2;i*i<=P;i++)
            if(P%i==0)
            {
                p[++tot]=i;
                q[tot]=0;
                while(P%i==0)q[tot]++, P/=i;
            }
        if(P>1)p[++tot]=P, q[tot]=1;
        rep(i,tot)
        {
            fact[i][0]=1;
            t[i]=1;
            rep(j,q[i])t[i]*=p[i];
            rep(j,maxn-1)
            {
                if(j%p[i]==0)fact[i][j]=fact[i][j-1];
                else fact[i][j]=fact[i][j-1]*j%t[i];
            }
        }
    }
    pll fact_mod(ll n, ll id)
    {
        if(n<p[id])return pll(0,fact[id][n]);
        pll ans = pll( n/p[id], em.fastpow(fact[id][t[id]],n/t[id],t[id]) );
        (ans.second*=fact[id][n%t[id]])%t[id];
        auto nex = fact_mod(n/p[id],id);
        ans.first+=nex.first;
        (ans.second*=nex.second)%=t[id];
        return ans;
    }
    ll exlucas(ll n, ll m, ll id)
    {
        ll cnt;
        auto a=fact_mod(n,id), b=fact_mod(m,id), c=fact_mod(n-m,id);
        a.first-=b.first+c.first;
        (a.second*=em.inv2(b.second*c.second%t[id],t[id]))%t[id];
        return em.fastpow(p[id],a.first,t[id])*a.second%t[id];
    }
    ll calc(ll n, ll m)
    {
        if(m>n or m<0 or n<0)return 0;
        vector<ll> a, v(t+1,t+tot+1);
        ll i;
        rep(i,tot)a.emb(exlucas(n,m,i));
        return em.CRT(a,v);
    }
}Cmod;
#define mod 10007ll
int main()
{
    ll n, m, T=read();
    Cmod.init(mod);
    while(T--)
    {
        n=read(), m=read();
        ll ans=Cmod.calc(n,m);
        printf("%lld\n",ans);
    }
    return 0;
}