Lucas定理初探
1.1 问题引入
已知\(p\)是一质数,求\(\dbinom{n}{m}\pmod{p}\).
关于组合数,它和排列数都是组合数学中的重要概念.这里会张贴有关这两个数的部分内容.
由于Lucas定理比较简单粗暴,这里直接给出表达式:
其中\(n=(n_1n_2\dots n_k)_p\),\(m=(m_1m_2\dots m_k)_p\),即\(n_i\)和\(m_i\)都是\(n\),\(m\)的\(p\)进制表示.
上面的公式其实是下面公式的递归版.仿照二进制,\(n \text{ div } p = (n_1n_2n_3\dots n_{k-1}n_k)_p \text{ div } p = (0n_1n_2\dots n_{k-1})_p\),相当于整个序列右移1位.而\(n \pmod{p}=n_k\),就是取出序列的最后一位.我们来证明一下这个定理.
我们把\(n\)的序列最后一项单独抠出来,\(n\)可以表示成\(n=\lfloor\frac{n}{p}\rfloor p + n\mod{p}=(n_1n_2\dots n_{k-1}0)_p+n_k\).
为了方便证明,我们把组合数套到二项式定理里面去.以下同余式均在模p意义下进行.
我们构造一个函数\(G_n(x)=(1+x)^n = \sum_{i=0}^n \dbinom{n}{i}x^i1^{n-i} = \sum_{i=0}^n \dbinom{n}{i}x^i\)
似乎这个函数和所谓的生成函数有关.
\(G_p(x)\equiv (1+x)^p \equiv \sum_{i=0}^p \dbinom{p}{i}x^i \equiv 1+ \sum_{i=1}^{p-1}\dbinom{p}{i}x^i + x^p\)
由于当\(i \neq 0 , p\)时,有\(\dbinom{p}{i} \equiv \dfrac{p!}{i!(p-i)!} \equiv p \dfrac{(p-1)!}{i!(p-i)!} \equiv 0\)
此时原式中间的那一大项式子可以直接消去,有:\((1+x)^p \equiv 1+x^p\).
由费马小定理,\((1+x)^p \equiv 1+x \pmod{p}\).
再结合上面的结论推到一下\(G_n(x)\)的表达式:
\(G_n(x)\equiv (1+x)^n \equiv (1+x)^{\lfloor\frac{n}{p}\rfloor p + n\pmod{p}}\)
\(\equiv (1+x)^{\lfloor\frac{n}{p}\rfloor}(1+x)^{n\pmod{p}}\)
\(\equiv (1+x^p)^{\lfloor\frac{n}{p}\rfloor}(1+x)^{n\pmod{p}}\)
\(\equiv \sum_{i=0}^{\lfloor\frac{n}{p}\rfloor}\dbinom{\lfloor\frac{n}{p}\rfloor}{i}(x^{p})^{i} \cdot \sum_{i=0}^{n\pmod{p}}\dbinom{n\pmod{p}}{i}x^i\)
\(\equiv G_{\lfloor\frac{n}{p}\rfloor}(x^p)\cdot G_{n\pmod{p}}(x)\)
其中函数\(G_s(x)\)中浓缩了一串数列:\(\{\dbinom{s}{0},\dbinom{s}{1},\dots,\dbinom{s}{s}\}\),每一项\(\dbinom{s}{i}\)是函数中\(x^i\)的系数.
因此,函数\(G_n(x)\)中\(x^k\)的系数是\(\dbinom{n}{k}\),假设复合函数\(G_{\lfloor\frac{n}{p}\rfloor}(x^p)\cdot G_{n\pmod{p}}(x)\)中包含\(x^k\)的一项可以写成:
\(a_kx^k=\dbinom{\lfloor\frac{n}{p}\rfloor}{u}(x^{p})^{u}\dbinom{n\pmod{p}}{v}x^v=\dbinom{\lfloor\frac{n}{p}\rfloor}{u}\dbinom{n\pmod{p}}{v}x^{pu+v} \qquad (※)\)
则必然有\(pu+v=k\).
同样的,我们取\(u=\lfloor\frac{k}{p}\rfloor\),\(v=k\pmod{p}\)
带入\((※)\)式:\(a_kx^k=\dbinom{\lfloor\frac{n}{p}\rfloor}{\lfloor\frac{k}{p}\rfloor}\dbinom{n\pmod{p}}{k\pmod{p}}x^k\).
我们就得到了一个新的函数$$H_n(x)=G_{\lfloor\frac{n}{p}\rfloor}(x^p)\cdot G_{n\pmod{p}}(x)=\sum_{i=1}{n}\dbinom{\lfloor\frac{n}{p}\rfloor}{\lfloor\frac{i}{p}\rfloor}\dbinom{n\pmod{p}}{i\pmod{p}}xi$$.
由上面的推论,在模\(p\)意义下,\(H_n(x) \equiv G_n(x)\)恒成立,所以对于每一项\(x^m\),其系数均有:\(\dbinom{n}{m}\equiv \dbinom{\lfloor\frac{n}{p}\rfloor}{\lfloor\frac{m}{p}\rfloor}\dbinom{n\pmod{p}}{m\pmod{p}}\)
2.1 拓展
对于\(p\)不是质数的情况，我们可以对p进行质因数分解：\(p=\prod_{i=1}^{k}p_i^{c_i}\)
那么原问题就变成了求解\(x=\dbinom{n}{m}\mod{(\prod_{i=1}^{k}p_i^{c_i})}\)
注意到每个\(p_i^{c_i}\)都必定是两两互素的数，我们先求\(\dbinom{n}{m}\pmod{p_i^{c_i}}\)，则原来的问题就变成了：$$\begin{cases}x\equiv \dbinom{n}{m} \pmod{p_1^{c_1}}\x\equiv \dbinom{n}{m} \pmod{p_2^{c_2}}\\cdots\x\equiv \dbinom{n}{m} \pmod{p_k^{c_k}}\end{cases}$$
我们的问题就变成了求解这个同余方程组。利用中国剩余定理，我们知道最后解的形式一定是\(x\equiv x_0 \pmod{(\prod_{i=1}^{k}p_i^{c_i})}\)。用中国剩余定理求出方程组的特解\(x_0\)就可以了。
如果求\(\dbinom{n}{m}\mod p_i^{c_i}\)，可以直接用阶乘计算：
\(\dbinom{n}{m}\mod{p_i^{c_i}}=\dfrac{n!\mod p_i^{c_i}}{(m!\mod{p_i^{c_i}})((n-m)!\mod{p_i^{c_i}})}\)
问题就是如何快速计算\(x!\mod{p_i^{c_i}}\)。
我们把阶乘的每一项展开来：\(x!=\prod_{k=1}^x{k}\)
试着每次从阶乘项里面提出\(p_i\)的倍数，进行计算：
\(x!=(p_i\cdot 2p_i\cdot 3p_i \cdots sp_i) \prod\limits_{1\leq k \leq x , p_i\nmid k}k\)
\(=p_i^s(s!)\prod\limits_{1\leq k \leq x , p_i\nmid k}k\)
由于模的循环性质，有：\(\prod\limits_{1\leq k \leq p_i^{c_i}}k\equiv \prod\limits_{p_i^{c_i} \leq k \leq 2p_i^{c_i}}k \equiv \cdots \pmod{p_i^{c_i}}\)
这样的循环节一共有\(\lfloor \frac{x}{p_i^{c_i}} \rfloor\)个。剩下的\(x-\lfloor \frac{x}{p_i^{c_i}} \rfloor\)就直接暴力计算就可以了。这有一点点分块的味道。
对于每一个循环节，我们只需要暴力计算其中一个\(p_i^{c_i}\)的部分，剩余的部分的答案是一样的。
典型地，我们采用网上流行的\(19!\mod 3^2\)来演示一下。
\(19!\equiv 1*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19\)
\(\equiv (1*2*4*5*7*8)*(10*11*13*14*16*17)*(19)*(3*6*9*12*15*18)\)
\(\equiv (1*2*4*5*7*8)*(1*2*4*5*7*8)*(1)*3^6(1*2*3*4*5*6)\)
这里出现了两个循环节和一个剩余的数。快速幂+暴力求出前面部分的值，后面\(3^6\)可以指数相减，阶乘可以递归处理。
每算一个\(n!\mod p_i^{c_i}\)时间复杂度是\(O(p_i^{c_i}\log n(\log_{p_i}n-c_i))\)。
这里就不加处理地粘贴了洛谷模板题。
#include<bits/stdc++.h>
#define RP(i,a,b) for(register int i=a;i<=b;++i)
#define DRP(i,a,b) for(register int i=a;i>=b;--i)
#define fre(z) freopen(z".in","r",stdin),freopen(z".out","w",stdout)
#define mem(a,b) memset(a,b,sizeof(a))
#define prique priority_queue
#define rg register
#include<cstdlib>
#include<ctime>
#ifdef WIN32
#define OT "%I64d"
#else
#define OT "%lld"
#endif
using namespace std;
typedef long long ll;
typedef double db;
template<class T> T qr(T type)
{
char ch=getchar();T ret=0;T q=1,res=1;
while(!isdigit(ch))
q=(ch=='-'?-1:q),ch=getchar();
while(isdigit(ch))
ret=ret*10+ch-'0',ch=getchar();
if(ch=='.')
{
ch=getchar();
while(isdigit(ch))
ret=ret*10+ch-'0',res=res*10,ch=getchar();
}
return q==-1?-ret/res:ret/res;
}
ll n,m,P;
inline ll exgcd(ll a,ll b,ll &x,ll &y)
{
if(b==0)
{
x=1,y=0;
return a;
}
ll d=exgcd(b,a%b,y,x);
y-=(a/b)*x;
return d;
}
inline void put(ll x)
{
if(x<0)
x=-x,putchar('-');
if(x>=10)
put(x/10);
putchar(x%10+'0');
}
inline ll qp(ll x,ll p,ll mod)
{
ll ans=1;
while(p)
{
if(p&1)
ans=(ans*x)%mod;
x=(x*x)%mod,p>>=1;
}
return ans;
}
ll MF(ll n,ll p,ll pc)//compute n! by pc norm
{
if(n==0)
return 1;
ll res=1;
for(rg ll i=2;i<=pc;++i)
if(i%p)//don't compute the p factor for it has been computed out of the function.
res=(res*i)%pc;//compute the single recurring period
res=qp(res,n/pc,pc);//compute the all recurring period
for(rg ll i=2;i<=n%pc;++i)//compute the res num. out of recper.
if(i%p)
res=(res*i)%pc;
return res*MF(n/p,p,pc)%pc;
}
inline ll inver(ll n,ll mod)
{
ll x,y;
exgcd(n,mod,x,y);
return (x+=mod)>mod?x-mod:x;
}
inline ll C(ll n,ll m,ll p,ll pc)
{
ll fc_n=MF(n,p,pc),fc_m=MF(m,p,pc),fc_del=MF(n-m,p,pc);//n!,m!,(n-m)! by pc norm
ll d=0;
for(rg ll i=n;i>0;i/=p)
d+=i/p;//d(n!)
for(rg ll i=m;i>0;i/=p)
d-=i/p;//d(m!)
for(rg ll i=n-m;i>0;i/=p)
d-=i/p;//d((n-m)!)
//k=d(n!)-d(m!)-d((n-m)!)
//d(n) is the fac p's num in n
return fc_n*qp(p,d,pc)%pc*inver(fc_m,pc)%pc*inver(fc_del,pc)%pc;
}
inline ll CRT(ll b,ll mod)
{
return b*inver(P/mod,mod)%P*(P/mod)%P;
}//add ans:aiMiti
/*
suppose P=\prod(p_i^{c_i})
the origin equation x=C(n,m)(mod P) is equivalance to:
{x=C(n,m)(mod p_i^{c_i})} for each factor is perpendicular
then we can solve the subans in turn and combine them.
ans=\sum_{i=1}^{n}{a_i M_i t_i}
M_i = P/(p_i^{c_i})
t_i is the inversion of M_i by p_i^{c_i} norm
*/
inline ll exlucas(ll n,ll m)//main function
{
ll res=0,cur=P,pc;
ll Sqr=sqrt(P)+5;
for(rg ll i=2;i<=Sqr;++i) //limit the searching range:sqrt(n) technique
{
if(cur%i==0)// i|cur
{
pc=1;
while(cur%i==0)
pc*=i,cur/=i;//p_i=i,p_i^{c_i}=pc
res=(res+CRT(C(n,m,i,pc),pc))%P;
}
}
if(cur>1)//cur=sqrt(P)
res=(res+CRT(C(n,m,cur,cur),cur))%P;
return res;
}
int main()
{
//fre("P4720");
n=qr(1ll),m=qr(1ll),P=qr(1ll);
put(exlucas(n,m));
return 0;
}
2.2 证明中出现的一些重要推论
阶乘\(n!\)中含有某个素因子\(p\)的个数\(\text{d}_p(x)\)。我们可以根据上面的过程，得到一个递归公式：
展开这个递归式可以得到\(\text{d}_p(x)\)的计算公式：
