扩展卢卡斯定理

扩展卢卡斯定理

前置芝士

卢卡斯定理，中国剩余定理

作用

和Lucas定理一样，只是 (C_m^n\%p) 中的 (p) 不一定是质数

结论

令 (p=p_1^{k_1}*p_2^{k_2}*...*p_q^{k_q})

列出同余方程组

[egin{cases}ansequiv c_1~(mod~p_1^{k_1})\ansequiv c_2~(mod~p_2^{k_2})\~~~~~~~~~~~~dots\ansequiv c_q~(mod~p_q^{k_q})end{cases} ]

其中 (c_1...c_q) 是对于每一个 (C_n^m\%p_i^{k_i}) 求出的答案，然后根据中国剩余定理合并，可见 (p_i^{k_i}) 不是质数，对于如何求 (C_n^m\%p_i^{k_i}) 将在证明中给出

证明

1. 由于同于方程组在模 (p=p_1^{k_1}*p_2^{k_2}*...*p_q^{k_q}) 意义下有唯一解，可以证明上面做法的正确性

2. 由于 (p=p_1^{k_1}*p_2^{k_2}*...*p_q^{k_q}) 实质上是将 (p) 进行质因数分解，所以 (p_i^{k_i}) 满足两两互质，可以直接用中国剩余定理合并

3. 如何求 (C_n^m\%p_i^{k_i})

   根据 (C_n^m=dfrac{n!}{m!(n-m)!}) ，只要分别求出 (n!\%p_i^{k_i},m!\%p_i^{k_i},(n-m)!\%p_i^{k_i}) 就可以通过逆元求出 (C_n^m\%p_i^{k_i})

4. 如何求 (n!\%p_i^{k_i})

   我们以 (n=19,p_i=3,k_i=2) 为例

   (n!=1∗2∗3∗4∗5∗6∗7∗8∗9∗10∗11∗12∗13∗14∗15∗16∗17∗18∗19\=(1∗2∗4∗5∗7∗8∗10∗11∗13∗14∗16∗17∗19)∗36∗(1∗2∗3∗4∗5∗6))

   根据这个例子发现，求解n!可以分为3部分：第一部分是 (p_i)的幂的部分，也就是 (3^6)即 (p_i^{lfloorfrac n{p_i} floor})，可以直接求解；第二部分是一个新的阶乘，也就是6!即 (lfloorfrac n{p_i} floor!)，可以递归下去求解；第三部分是除前两部分之外剩下的数

5. 考虑第三部分如何求解

   发现第三部分在模 (p_i^{k_i}) 意义下是以 (p_i^{k_i}) 为周期的，即:

   ((1∗2∗4∗5∗7∗8)≡(10∗11∗13∗14∗16∗17)~(mod~p_i^{k_i})) ，所以只要求出 (p_i^{k_i}) 的长度即可；

   但是还剩下一个孤立的19，可以发现剩下孤立的数长度不会超过 (p_i^{k_i}) ，只需要暴力求解即可

6. 最后一个问题是对于求出的 (m!\%p_i^{k_i}) 和 ((n-m)!\%p_i^{k_i}) 有可能与 (p_i^{k_i}) 不互质，无法求逆元

   所以要将 ((n-m)!\%p_i^{k_i}) 和 (m!\%p_i^{k_i}) 中质因子 (p_i) 先全部除去，求出逆元后再全部乘回去

   计算n!中质因子p的个数x的公式为 (x=lfloorfrac np floor+lfloorfrac n{p^2} floor+lfloorfrac n{p^3} floor+...)
   递推式也可以写为 (f(n)=f(lfloorfrac np floor)+lfloorfrac np floor)

代码

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
#define LL long long

LL n,m,MOD,ans;

LL fast_pow(LL a,LL p,LL Mod)
{
    LL ans=1LL;
    for (;p;p>>=1,a=a*a%Mod)
        if (p&1)
            ans=ans*a%Mod;
    return ans;
}
void exgcd(LL a,LL b,LL &x,LL &y)
{
    if (!b) x=1LL,y=0LL;
    else exgcd(b,a%b,y,x),y-=a/b*x;
}
LL inv(LL A,LL Mod)
{
    if (!A) return 0LL;
    LL a=A,b=Mod,x=0LL,y=0LL;
    exgcd(a,b,x,y);
    x=((x%b)+b)%b;
    if (!x) x+=b;
    return x;
}
LL Mul(LL n,LL pi,LL pk)
{
    if (!n) return 1LL;
    LL ans=1LL;
    if (n/pk)
    {
        for (LL i=2;i<=pk;++i)
            if (i%pi) ans=ans*i%pk;
        ans=fast_pow(ans,n/pk,pk);
    }
    for (LL i=2;i<=n%pk;++i)
        if (i%pi) ans=ans*i%pk;
    return ans*Mul(n/pi,pi,pk)%pk;
}
LL C(LL n,LL m,LL Mod,LL pi,LL pk)
{
    if (m>n) return 0LL;
    LL a=Mul(n,pi,pk),b=Mul(m,pi,pk),c=Mul(n-m,pi,pk);
    LL k=0LL,ans;
    for (LL i=n;i;i/=pi) k+=i/pi;
    for (LL i=m;i;i/=pi) k-=i/pi;
    for (LL i=n-m;i;i/=pi) k-=i/pi;
    ans=a*inv(b,pk)%pk*inv(c,pk)%pk*fast_pow(pi,k,pk)%pk;
    return ans*(Mod/pk)%Mod*inv(Mod/pk,pk)%Mod;
}
int main()
{
    scanf("%I64d%I64d%I64d",&n,&m,&MOD);		//因为cf里面好像lld不太行，所以要用I64d
    for (LL x=MOD,i=2;i<=MOD;++i)
        if (x%i==0)
        {
            LL pk=1LL;
            while (x%i==0) pk*=i,x/=i;
            ans=(ans+C(n,m,MOD,i,pk))%MOD;
        }
    printf("%I64d
",ans);
}