Lucas定理&扩展Lucas学习笔记

Lucas定理

内容

若$p$为素数,则对于任意整数$1\leq m\leq n$,都有$C_n^m\equiv C_{n\%p}^{m\%p}*C_{n/p}^{m/p}\pmod p$

也就是说,把$n$和$m$表示成$p$进制数,每一位分别算组合数,再乘起来。

证明比较复杂,我不会

代码

inline int Lucas(int n,int m,int p) {
    if (!m) return 1;
    return C(n%p,m%p,p)*Lucas(n/p,m/p,p)%p;
}

C()是计算组合数的函数,大家应该都知道怎么写。

扩展Lucas

Lucas定理只有在模质数意义下才有效!不是质数就GG了!​

于是有了扩展Lucas,它可以解决模数不为质数的情况。

前方高能预警!

前置知识

  • 扩展欧几里得算法$\texttt{exgcd}$
  • 中国剩余定理$\texttt{CRT}$

算法

还是设模数为$p$。

分解质因数

非常简单,就是把$p$分解成$\prod p_i^{c_i}$的形式。

计算组合数

在这一步中,我们只考虑计算$C_n^m\mod p_i^{c_i}$的值。

直接上$\texttt{Lucas}$定理显然不可行,考虑直接拆式子:

$\large C_n^m=\frac{n!}{(n-m)!m!}$

问题在于怎么快速求出阶乘和阶乘的逆元。

阶乘的逆元用$\texttt{exgcd}$随便搞一下即可,问题进一步转换为求阶乘。

先将阶乘中所有$p$的倍数提出来,显然只有$\left\lfloor\frac{n}{p}\right\rfloor$个。每一项除去$p$,得到了$\left\lfloor\frac{n}{p}\right\rfloor!$,然后就可以递归求解了。

剩下的非$p$的倍数的项,显然存在长度小于$p_i^{c_i}$的循环节,可以丢到一起去算。

最后还剩下几项凑不成一个循环节,这一部分直接暴力算就好了。

还需要考虑$p$的次幂一共出现了多少次。设$f(n)$表示$n!$中有多少个$p$的因数,可以得到:

$\large f(n)=f(\left\lfloor\frac{n}{p}\right\rfloor)+\left\lfloor\frac{n}{p}\right\rfloor$

边界为$f(x)=0\;(x<p)$。

结合代码理解一下(当然只放了关键函数):

inline int fact(int n,int pi,int pk) {
    if (!n) return 1;
    int ans=1;
    for (re int i=2;i<=pk;++i)
        if (i%pi) ans=ans*i%pk;
    ans=FastPow(ans,n/pk,pk);
    for (re int i=2;i<=n%pk;++i)
        if (i%pi) ans=ans*i%pk;
    return ans*fact(n/pi,pi,pk)%pk;
}

inline int C(int n,int m,int pi,int pk) {
    int d1=fact(n,pi,pk),d2=fact(m,pi,pk),d3=fact(n-m,pi,pk);
    int k=0;
    for (re int i=n;i;i/=pi) k+=i/pi;
    for (re int i=m;i;i/=pi) k-=i/pi;
    for (re int i=n-m;i;i/=pi) k-=i/pi;
    return d1*getInv(d2,pk)%pk*getInv(d3,pk)%pk*FastPow(pi,k,pk)%pk;
}

还是看不懂?请结合下面的完整代码理解。

合并

我们已经求出了一堆的$x\mod p_i^{c_i}=a_i$这样的式子。

发现这是个线性同余方程组,于是用$\texttt{CRT}$合并,解出$x$就是答案。

代码

inline int exgcd(int a,int b,int& x,int& y) {
    if (!b) { x=1,y=0; return a; }
    int d=exgcd(b,a%b,x,y);
    int z=x; x=y; y=z-(a/b)*y;
    return d;
}

inline int FastPow(int a,int b,int m) {
    int ans=1;
    for (;b;b>>=1,a=a*a%m)
        if (b&1) ans=ans*a%m;
    return ans;
}

inline int fact(int n,int pi,int pk) {
    if (!n) return 1;
    int ans=1;
    for (re int i=2;i<=pk;++i)
        if (i%pi) ans=ans*i%pk;
    ans=FastPow(ans,n/pk,pk);
    for (re int i=2;i<=n%pk;++i)
        if (i%pi) ans=ans*i%pk;
    return ans*fact(n/pi,pi,pk)%pk;
}

inline int getInv(int n,int m) {
    int x,y; exgcd(n,m,x,y);
    return x>0?x:x+m;
}

inline int CRT(int b,int m) { return b*getInv(p/m,m)%p*(p/m)%p; }

inline int C(int n,int m,int pi,int pk) { //C(n,m)%pk
    int d1=fact(n,pi,pk),d2=fact(m,pi,pk),d3=fact(n-m,pi,pk);
    int k=0;
    for (re int i=n;i;i/=pi) k+=i/pi;
    for (re int i=m;i;i/=pi) k-=i/pi;
    for (re int i=n-m;i;i/=pi) k-=i/pi;
    return d1*getInv(d2,pk)%pk*getInv(d3,pk)%pk*FastPow(pi,k,pk)%pk;
}

inline int EX_Lucas(int n,int m) {
    int ans=0,tmp=p,pk;
    int end=sqrt(p)+5;
    for (re int i=2;i<=end;++i) {
        if (tmp%i==0) {
            pk=1; while (tmp%i==0) pk*=i,tmp/=i;
            ans=(ans+CRT(C(n,m,i,pk),pk))%p;
        }
    }
    if (tmp>1) ans=(ans+CRT(C(n,m,tmp,tmp),tmp))%p;
    return ans;
}

如果你看不懂,背下来就行了

最后修改:2019 年 11 月 12 日 03 : 11 PM

发表评论