卢卡斯定理学习笔记

编辑文章

概述

卢卡斯定理主要用于求解组合数取模问题。

公式

原命题见: https://zh.wikipedia.org/wiki/%E5%8D%A2%E5%8D%A1%E6%96%AF%E5%AE%9A%E7%90%86

原命题等价于:

$$\binom{m}{n}=\binom{\lfloor \dfrac{m}{p}\rfloor}{\lfloor \dfrac{n}{p}\rfloor}\times \binom{m\mod p}{n\mod p} \pmod p$$

其中 $p$ 是质数。

$\binom{m\mod p}{n\mod p}$ 中上下两项都 $< p$ ,当 $p\le 10^{7}$ 时,可以用 $O(n)$ 处理阶乘进行计算。剩下的可以不断递归进行计算。

inline int C(int n,int m)
{
    if (m>n) return 0;
    return fac[n]*(inv(fac[m])*inv(fac[n-m])%ha)%ha;
}

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

拓展卢卡斯定理

当 $p$ 很大或不是质数时,显然不能用上述方法算。

当 $p$ 不是质数,可以对其进行分解

$$p=\prod_{i=1}^n p_i^{k_i}$$

然后计算每一个 $\mod p_i^{k_i}$ 后的值,最后用中国剩余定理合并答案即可。

当 $p$ 很大时,需要对阶乘的处理进行优化。

假设当前的模数为 $p^k$ ,需要求 $n!$ 。

可以先将 $p$ 的倍数提出来,并记录个数为 $s$ 。

int s=0; //pi的个数
for (int i=n;i;i/=p) s+=i/p;

$p$ 的倍数 $\div p$ 后剩下的数的乘积即为 $\lfloor \dfrac{n}{p} \rfloor !$

然后可以对剩下的数按照 $\div p^k$ 的值进行分块,长度为 $p^k$ ,计算一个块内乘积 $\mod p$ 的值就可以得到所有块的值。

还剩下 $n\mod p^k$ 个数不在块内,这些数可以单独处理。

inline int fac(int x,int pi,int pk)
{
    if (!x) return 1;
    int ans=1;
    for (int i=2;i<=pk;i++) //处理整块
    {
        if (i%pi==0) continue;
        ans=ans*i%pk;
    }
    ans=qpow(ans,x/pk,pk); //所有块的乘积
    int len=x%pk;
    for (int i=2;i<=len;i++) //单独处理不在块内的元素
    {
        if (i%pi==0) continue;
        ans=ans*i%pk;
    }
    return ans*fac(x/pi,pi,pk)%pk; //乘上/p剩下的数
}

组合数的计算公式为

$$\binom{m}{n}=\dfrac{n!}{m!\times (n-m)!}$$

显然涉及到除法,那么之前统计的 $s$ 就有用了。对于 $n!$ 就加,$m!$ 和 $(n-m)!$ 就减去 $p$ 的个数就行了。

inline int C(int n,int m,int pi,int pk)
{
    int s=0; //pi的个数
    for (int i=n;i;i/=pi) s+=i/pi;
    for (int i=m;i;i/=pi) s-=i/pi;
    for (int i=n-m;i;i/=pi) s-=i/pi;
    return fac(n,pi,pk)*inv(fac(m,pi,pk),pk)%pk*inv(fac(n-m,pi,pk),pk)%pk*qpow(pi,s,pk)%pk;
}

主体部分就是对 $p$ 进行分解,然后对每个因数记录 $p_i$ 为 pi,$p_i^{k_i}$ 为 pk,再对组合数进行计算,最后中国剩余定理合并答案。

模板:洛谷4720 【模板】扩展卢卡斯

#include<bits/stdc++.h>
#define ll long long

using namespace std;

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

ll p,zs[1005],cnt;
bool ss[1005];

inline void get_prime()
{
    memset(ss,1,sizeof(ss));
    ss[0]=ss[1]=0;
    ll sqn=sqrt(p);
    for (ll i=2;i<=sqn;i++)
    {
        if (ss[i]) zs[++cnt]=i;
        for (ll j=1;j<=cnt && zs[j]*i<=sqn;j++)
        {
            ss[zs[j]*i]=0;
            if (i%zs[j]==0) break;
        }
    }
}

ll exgcd(ll l,ll r,ll &x,ll &y)
{
    if (r==0)
    {
        x=1; y=0;
        return l;
    }
    ll ans=exgcd(r,l%r,y,x);
    y-=l/r*x;
    return ans;
}

inline ll qpow(ll x,ll y,ll ha)
{
    ll ans=1;
    while (y)
    {
        if (y&1) ans=ans*x%ha;
        y>>=1;
        x=x*x%ha;
    }
    return ans%ha;
}

inline ll inv(ll now,ll ha)
{
    ll x=0,y=0;
    exgcd(now,ha,x,y);
    return (x%ha+ha)%ha;
}

inline ll fac(ll x,ll pi,ll pk)
{
    if (!x) return 1;
    ll ans=1;
    for (ll i=2;i<=pk;i++) //整块
    {
        if (i%pi==0) continue;
        ans=ans*i%pk;
    }
    ans=qpow(ans,x/pk,pk);
    ll len=x%pk;
    for (ll i=2;i<=len;i++)
    {
        if (i%pi==0) continue;
        ans=ans*i%pk;
    }
    return ans*fac(x/pi,pi,pk)%pk;
}

inline ll C(ll n,ll m,ll pi,ll pk)
{
    ll s=0; //pi的个数
    for (ll i=n;i;i/=pi) s+=i/pi;
    for (ll i=m;i;i/=pi) s-=i/pi;
    for (ll i=n-m;i;i/=pi) s-=i/pi;
    return fac(n,pi,pk)*inv(fac(m,pi,pk),pk)%pk*inv(fac(n-m,pi,pk),pk)%pk*qpow(pi,s,pk)%pk;
}

inline ll exlucas(ll n,ll m)
{
    ll ans=0,p2=p;
    for (ll i=1;i<=cnt;i++)
    {
        if (p2%zs[i]) continue;
        ll pk=1;
        while (p2%zs[i]==0) p2/=zs[i],pk*=zs[i];
        ll ai=C(n,m,zs[i],pk),mi=p/pk;
        ans=(ans+ai*mi%p*inv(mi,pk)%p)%p;
    }
    if (p2>1)
    {
        ll ai=C(n,m,p2,p2),mi=p/p2;
        ans=(ans+ai*mi%p*inv(mi,p2)%p)%p;
    }
    return ans;
}

signed main()
{
    ll n,m;
    n=read(); m=read(); p=read();
    get_prime();
    return !printf("%lld",exlucas(n,m));
}

新评论

称呼不能为空
邮箱格式不合法
网站格式不合法
内容不能为空

第一次提交的评论将在审核后显示。