Lucas定理&&扩展Lucas

Lucas定理

求解$C_n^m\mathrm{m}\mathrm{o}\mathrm{d}p$ （p为素数）

蒟蒻不会证明，记着递推公式就好=_=
$Lucas(n,m)=C(n,m)\mathrm{m}\mathrm{o}\mathrm{d}p$$=C(n\mathrm{m}\mathrm{o}\mathrm{d}p$$,m\mathrm{m}\mathrm{o}\mathrm{d}p$$)\ast Lucas(n/p,m/p)\mathrm{m}\mathrm{o}\mathrm{d}p$
边界为$Lucas(i,0)=1$

int C(int n,int m)
{
    int ans=1;
    for(int i=1;i<=m;++i)
    ans*=(n-m+i)*inv[i]%p;
    return ans;
} 

int lucas(int n,int m)
{
    if(m==0) return 1;
    else return C(n%p,m%p)*lucas(n/p,m/p)%p;
}

---

扩展Lucas

求解$C_n^m\mathrm{m}\mathrm{o}\mathrm{d}p$ （p不一定为素数）

先导知识 中国剩余定理、扩欧、逆元、组合数、质因数分解

首先讲p分解质因数
$p=p_1^{k_1}+p_2^{k_2}+...+p_q^{k_q}$

于是可以列出同余方程组
$$\{\begin{array}{r}ans\equiv C_n^m\mathrm{m}\mathrm{o}\mathrm{d}{p}_1^{k_1}(\mathrm{m}\mathrm{o}\mathrm{d}{p}_1^{k_1})\\ ans\equiv C_n^m\mathrm{m}\mathrm{o}\mathrm{d}{p}_2^{k_2}(\mathrm{m}\mathrm{o}\mathrm{d}{p}_2^{k_2})\\ ...\\ ans\equiv C_n^m\mathrm{m}\mathrm{o}\mathrm{d}{p}_q^{kq}(\mathrm{m}\mathrm{o}\mathrm{d}{p}_q^{kq})\end{array}$$

显然ans就是我们所要的答案
因为$p_i^{k_i}$是分解质因数得到的，所以必定两两互质
只要用中国剩余定理合并就好

我们要讨论的是怎样求出$C_n^m\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{k_i}$
直接考虑从定义式入手
$C_n^m\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{k_i}=\frac{n!\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{k_i}}{m!\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{k_i}\ast (n-m)!\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{k_i}}$

直接看例子会比较容易理解,看这个例子$n=19,pi=3,ki=2$

先试着把$n!$中所有$p_i$的倍数分离出来
得$3\ast 6\ast 9\ast 12\ast 15\ast 18$
再把这里面的$p_i$全部提出来
得$3^6\ast (1\ast 2\ast 3\ast 4\ast 5\ast 6)=3^6\ast 6!$

形式化的描述这部分就是$p_i^{\lfloor n/p_i\rfloor }\ast \lfloor n/p_i\rfloor !$
$p_i$得幂次方部分可以直接快速幂，阶乘部分可以继续递归

处理完上面来看剩下一部分
$1\ast 2\ast 4\ast 5\ast 7\ast 8\ast 10\ast 11\ast 13\ast 14\ast 16\ast 17\ast 19$
仔细观察可以发现这一串在膜$p{i}^{ki}$意义下是有周期的
$(1\ast 2\ast 4\ast 5\ast 7\ast 8)\equiv (10\ast 11\ast 13\ast 14\ast 16\ast 17)(\mathrm{m}\mathrm{o}\mathrm{d}p{i}^{ki})$
每一个循环节长度为$p{i}^{ki}$
只要求出第一个循环节，然后在快速幂$\lfloor n/p{i}^{k_i}\rfloor$次方就可以了

最后发现还剩下一个19
也就是说我们算完上面两部分还有剩余$n-\lfloor n/p{i}^{k_i}\rfloor \ast p{i}^{k_i}$个数字
不过这些数字长度显然不超过$p{i}^{ki}$，直接暴力就好

最后再次总结一下分解的几部分

1.一个幂次方以及一个阶乘，$p_i^{\lfloor n/p_i\rfloor }\ast \lfloor n/p_i\rfloor !$

2.一个循环节的$\lfloor n/p_k\rfloor$次方

3.剩下的$n-\lfloor n/p{i}^{k_i}\rfloor \ast p{i}^{k_i}$个数字暴力乘

洛谷P4720 【模板】扩展卢卡斯

#include<iostream>
#include<cmath>
#include<algorithm>
#include<queue>
#include<cstring>
#include<cstdio>
using namespace std;
typedef long long lt;

lt read()
{
    lt f=1,x=0;
    char ss=getchar();
    while(ss<'0'||ss>'9'){if(ss=='-')f=-1;ss=getchar();}
    while(ss>='0'&&ss<='9'){x=x*10+ss-'0';ss=getchar();}
    return f*x;
}

const int maxn=50001;
lt a[maxn],b[maxn],cnt;

lt qpow(lt ai,lt k,lt mod)
{
    lt mul=1;
    while(k>0)
    {
        if(k&1)mul=(mul*ai)%mod;
        ai=(ai*ai)%mod;
        k>>=1;
    }
    return mul;
}

lt fac(lt n,lt pi,lt pk)
{
    if(!n) return 1;
    lt mul=1;
    
    for(lt i=2;i<=pk;++i)//分解阶乘第二部分，循环节
    if(i%pi)mul=(mul*i)%pk;
    mul=qpow(mul,n/pk,pk);
    
    for(lt i=2;i<=n%pk;++i)//分解阶乘第三部分，求剩余数字
    if(i%pi)mul=(mul*i)%pk;
    
    return mul*fac(n/pi,pi,pk)%pk;//分解阶乘第一部分的另一个阶乘递归
}

void exgcd(lt a,lt b,lt &x,lt &y)
{
    if(b==0){ x=1; y=0; return;}
    exgcd(b,a%b,x,y);
    lt tp=x;
    x=y; y=tp-a/b*y;
}

lt inv(lt a,lt b)
{
    lt x,y;
    exgcd(a,b,x,y);
    return (x%b+b)%b;
}

lt C(lt n,lt m,lt pi,lt pk)
{
    lt facn=fac(n,pi,pk);//分别求n,m,n-m膜pi^ki的阶乘
    lt facm=fac(m,pi,pk);
    lt facnm=fac(n-m,pi,pk);
    
    lt kk=0;
    for(lt i=n;i;i/=pi)kk+=i/pi;//上述分解阶乘第一部分的pi幂次方
    for(lt i=m;i;i/=pi)kk-=i/pi;
    for(lt i=n-m;i;i/=pi)kk-=i/pi;
    
    return facn*inv(facm,pk)%pk*inv(facnm,pk)%pk*qpow(pi,kk,pk)%pk;//注意求逆元
}

void div(lt n,lt m,lt x)
{
    for(lt i=2;i<=sqrt(x);++i)
    {
        if(x%i==0)
        {
            lt pi=i,ki=0;
            while(x%i==0)x/=i,ki++;
            b[++cnt]=qpow(pi,ki,1e7);
            a[cnt]=C(n,m,pi,b[cnt]);//C(n,m)%pi^ki
        }
    }
    if(x>1)b[++cnt]=x,a[cnt]=C(n,m,x,b[cnt]);
}

lt china()
{
    lt ans=0,M=1,x,y;
    for(int i=1;i<=cnt;++i) M*=b[i];
    for(int i=1;i<=cnt;++i)
    {
        lt tp=M/b[i];
        exgcd(tp,b[i],x,y);
        x=(x%b[i]+b[i])%b[i];
        ans=(ans+tp*x*a[i])%M;
    }
    return (ans+M)%M;
}

int main()
{
    lt n=read(),m=read(),p=read();
    div(n,m,p);
    printf("%lld",china());
    return 0;
}