学好数理化，走遍天下都不怕

exLucas定理

发表于 更新于 分类于

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

前置知识

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

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

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

规定

下面的讲解有几个规定：

  1. $p$为模数
  2. 除了$p$以外的所有含$p$的字母默认为质数
  3. 所有的数都是非负整数
  4. $i$的范围默认为距离这个$i$最近（上面）的一次定义的$i$的范围（可能在式子中）
  5. 与第四条类似，①式默认为距离这个①式最近（上面）的一次定义的①式
分解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$是多余的部分，也直接暴力枚举就完了
最后递归一下

1
2
3
4
5
6
7
8
9
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’$
1
2
3
4
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}$
所以我们可以得出整个计算①式的程序

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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;
}

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
#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);
}