exLucas定理 | Maserhe

数论组合

数学

发布日期: 2020-04-13

作者: Maserhe

文章字数: 1.5k

阅读时长: 7 分

模板
这好像不是个定理……
看起来是个算法，但是不知道为什么大家都叫它定理，所以我也跟着写了……
还有，这玩意儿好像跟Lucas定理没有半毛钱关系，不知道为什么叫exLucas……

前置知识

为了学会这个“定理”，你需要了解几样东西：

质因数分解
CRT（中国剩余定理）
逆元
其他基础的数论
组合数的定义（？？？）
好的视力，或者一个放大镜（因为有的又有分数又有上下标，很容易看不清）

这些知识如果不会，可以参考我的另一篇博客（5、6除外）

规定

下面的讲解有几个规定：

$p$为模数
除了$p$以外的所有含$p$的字母默认为质数
所有的数都是非负整数
$i$的范围默认为距离这个$i$最近（上面）的一次定义的$i$的范围（可能在式子中）
与第四条类似，①式默认为距离这个①式最近（上面）的一次定义的①式

分解p

首先，我们要先将$p$进行质因数分解
设$p=\prod \limits_{i=1}^{t}p_i^{\alpha_i}$
这样，我们就只需要计算$\binom{n}{m}modp_i^{\alpha_i}(*)$，最后再用一下CRT就好了

算(∗)式！

看到这个式子，有的同学就会说：看！$p_i$是质数！直接用Lucas定理算！
……首先，我说过，这玩意儿跟Lucas定理没有半毛钱关系；其次，它还有个次数$\alpha_i$啊！
首先，我们由组合数的定义可以知道：$\binom{n}{m}=\frac{n!}{m!(n-m)!}$
这时候，又有同学要说了：$\frac{1}{m!(n-m)!}$直接求逆元啊！
……求逆元的前提是互质啊！你看$m!(n-m)!$这个阶乘，它能跟$p_i$互质吗？
所以，我们得把这玩意儿中的$p_i$提出来……
假设$p^{k1}\mid\mid n!,p^{k2}\mid\mid m!,p^{k3}\mid\mid (n-m)!$
那么$\binom{n}{m}=\frac{\frac{n!}{p^{k1}}}{\frac{m!}{p^{k2}}\times\frac{(n-m)!}{p^{k3}}}\times p^{k1-k2-k3}$

阶乘中质数的个数

假设$p^\alpha\mid\mid n!$
首先考虑$n!%p^k$咋算
$n!%{p^k}={p^{n/p}}\times (n/p)!\times (\prod\limits_{i,gcd(i,p)=1}^{p^k}i)^{n/p^k}\times(\prod\limits_{i,gcd(i,p)=1}^{p%k}i)%p^k$
其中，$\prod\limits_{i,gcd(i,p)=1}^{p^k}i$是循环的，直接暴力枚举，然后快速幂
剩下的$\prod\limits_{i,gcd(i,p)=1}^{p%k}i$是多余的部分，也直接暴力枚举就完了
最后递归一下

long long fac(long long n,long long p,long long ppa)//计算n!%p^ppa
{
    if (n==0) return 1;
    long long cir=1/*循环节*/,rem=1/*余数*/;
    for (long long i=1;i<=ppa;i++) if(i%p) cir=cir*i%ppa;
    cir=POW(cir,n/ppa,ppa);
    for(long long i=ppa*(n/ppa);i<=n;i++) if(i%p) rem=rem*(i%ppa)%ppa;
    return fac(n/p,p,ppa)*cir%ppa*rem%ppa;
}

接着，我们需要算出$\alpha$
这个就非常简单了，我们假设$p^{\alpha’}\mid\mid \left[\frac{n}{p}\right]!$
那么$\alpha=\left[\frac{n}{p}\right]+\alpha’$

long long sum_fac(long long n,long long p)//n!中p的个数
{
    return n<p?0:sum_fac(n/p,p)+(n/p);
}

回到①式，$\binom{n}{m}=\frac{\frac{n!}{p^{k1}}}{\frac{m!}{p^{k2}}\times\frac{(n-m)!}{p^{k3}}}\times p^{k1-k2-k3}$
所以我们可以得出整个计算①式的程序

long long fac(long long n,long long p,long long ppa)//计算n!
{
    if (n==0) return 1;
    long long cir=1/*循环节*/,rem=1/*余数*/;
    for (long long i=1;i<=ppa;i++) if(i%p) cir=cir*i%ppa;
    cir=POW(cir,n/ppa,ppa);
    for(long long i=ppa*(n/ppa);i<=n;i++) if(i%p) rem=rem*(i%ppa)%ppa;
    return fac(n/p,p,ppa)*cir%ppa*rem%ppa;
}
long long sum_fac(long long n,long long p)//n!中p的个数
{
    return n<p?0:sum_fac(n/p,p)+(n/p);
}
long long C(long long n,long long m,long long p,long long ppa)//计算Cnm%pi^ai
{
    long long fz=fac(n,p,ppa),fm1=ny(fac(m,p,ppa),ppa),fm2=ny(fac(n-m,p,ppa),ppa);
    long long mi=POW(p,sum_fac(n,p)-sum_fac(m,p)-sum_fac(n-m,p),ppa);
    return fz*fm1%ppa*fm2%ppa*mi%ppa;
}

解决了①式，剩下的部分就迎刃而解了

#include<algorithm>
#include<iostream>
using namespace std;
long long n,m,p,cnt/*个数*/,pr[1010]/*质数*/,al[1010]/*指数*/;
void exgcd(long long a,long long b,long long &x,long long &y)//扩展欧几里得算法
{
    if (!b) return (void)(x=1,y=0);
    exgcd(b,a%b,x,y);
    long long tmp=x;x=y;y=tmp-a/b*y;
}
long long ny(long long a,long long p)//逆元
{
    long long x,y;
    exgcd(a,p,x,y);
    return (x+p)%p;
}
long long POW(long long a,long long b,long long p)//快速幂
{
    long long t=1;
    a%=p;
    for(;b;b>>=1){
        if(b&1) t=t*a%p;
        a=a*a%p;
    }
    return t;
}
long long fac(long long n,long long p,long long ppa)//计算n!
{
    if (n==0) return 1;
    long long cir=1/*循环节*/,rem=1/*余数*/;
    for (long long i=1;i<=ppa;i++) if(i%p) cir=cir*i%ppa;
    cir=POW(cir,n/ppa,ppa);
    for(long long i=ppa*(n/ppa);i<=n;i++) if(i%p) rem=rem*(i%ppa)%ppa;
    return fac(n/p,p,ppa)*cir%ppa*rem%ppa;
}
long long sum_fac(long long n,long long p)//n!中p的个数
{
    return n<p?0:sum_fac(n/p,p)+(n/p);
}
long long C(long long n,long long m,long long p,long long ppa)//计算Cnm%pi^ai
{
    long long fz=fac(n,p,ppa),fm1=ny(fac(m,p,ppa),ppa),fm2=ny(fac(n-m,p,ppa),ppa);
    long long mi=POW(p,sum_fac(n,p)-sum_fac(m,p)-sum_fac(n-m,p),ppa);
    return fz*fm1%ppa*fm2%ppa*mi%ppa;
}
void pfd(long long n,long long m)//分解p
{
    long long P=p;
    for(long long i=2;i*i<=p;i++){
        if(!(P%i)){
            long long ppa=1;
            while(!(P%i)) ppa*=i,P/=i;
            pr[++cnt]=ppa;
            al[cnt]=C(n,m,i,ppa);
        }
    }
    if(P!=1) pr[++cnt]=P,al[cnt]=C(n,m,P,P);
}
long long crt()//中国剩余定理
{
    long long ans=0;
    for(long long i=1;i<=cnt;i++){
        long long M=p/pr[i],T=ny(M,pr[i]);
        ans=(ans+al[i]*M%p*T%p)%p;
    }
    return ans;
}
long long exlucas(long long n,long long m)//扩展卢卡斯 
{
    pfd(n,m);
    return crt();
}
int main()
{
    cin>>n>>m>>p;
    cout<<exlucas(n,m);
}