Lucas定理初探

1.1 问题引入

已知(p)是一质数,求(dbinom{n}{m}pmod{p}).

关于组合数,它和排列数都是组合数学中的重要概念.这里会张贴有关这两个数的部分内容.

由于Lucas定理比较简单粗暴,这里直接给出表达式:

[oxed{dbinom{n}{m}equivdbinom{nmod p}{mmod p}dbinom{frac{n}{p}}{frac{m}{p}}pmod{p}}$$. 或者: $$oxed{dbinom{n}{m}equiv dbinom{n_1}{m_1}dbinom{n_2}{m_2}dotsdbinom{n_k}{m_k} pmod{p}}]

其中(n=(n_1n_2dots n_k)_p),(m=(m_1m_2dots m_k)_p),即(n_i)和(m_i)都是(n),(m)的(p)进制表示.

上面的公式其实是下面公式的递归版.仿照二进制,(n ext{ div } p = (n_1n_2n_3dots n_{k-1}n_k)_p ext{ div } p = (0n_1n_2dots n_{k-1})_p),相当于整个序列右移1位.而(n pmod{p}=n_k),就是取出序列的最后一位.我们来证明一下这个定理.

我们把(n)的序列最后一项单独抠出来,(n)可以表示成(n=lfloorfrac{n}{p} floor p + nmod{p}=(n_1n_2dots 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 eq 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)^{lfloorfrac{n}{p} floor p + npmod{p}})
(equiv (1+x)^{lfloorfrac{n}{p} floor}(1+x)^{npmod{p}})
(equiv (1+x^p)^{lfloorfrac{n}{p} floor}(1+x)^{npmod{p}})
(equiv sum_{i=0}^{lfloorfrac{n}{p} floor}dbinom{lfloorfrac{n}{p} floor}{i}(x^{p})^{i} cdot sum_{i=0}^{npmod{p}}dbinom{npmod{p}}{i}x^i)
(equiv G_{lfloorfrac{n}{p} floor}(x^p)cdot G_{npmod{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_{lfloorfrac{n}{p} floor}(x^p)cdot G_{npmod{p}}(x))中包含(x^k)的一项可以写成:
(a_kx^k=dbinom{lfloorfrac{n}{p} floor}{u}(x^{p})^{u}dbinom{npmod{p}}{v}x^v=dbinom{lfloorfrac{n}{p} floor}{u}dbinom{npmod{p}}{v}x^{pu+v} qquad (※))
则必然有(pu+v=k).
同样的,我们取(u=lfloorfrac{k}{p} floor),(v=kpmod{p})
带入((※))式:(a_kx^k=dbinom{lfloorfrac{n}{p} floor}{lfloorfrac{k}{p} floor}dbinom{npmod{p}}{kpmod{p}}x^k).
我们就得到了一个新的函数$$H_n(x)=G_{lfloorfrac{n}{p} floor}(x^p)cdot G_{npmod{p}}(x)=sum_{i=1}^{{n}dbinom{lfloorfrac{n}{p} floor}{lfloorfrac{i}{p} floor}dbinom{npmod{p}}{ipmod{p}}x}i$$.
由上面的推论,在模(p)意义下,(H_n(x) equiv G_n(x))恒成立,所以对于每一项(x^m),其系数均有:(dbinom{n}{m}equiv dbinom{lfloorfrac{n}{p} floor}{lfloorfrac{m}{p} floor}dbinom{npmod{p}}{mpmod{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}})，则原来的问题就变成了：$$egin{cases}xequiv dbinom{n}{m} pmod{p_1^{c_1}}xequiv dbinom{n}{m} pmod{p_2^{c_2}}\cdotsxequiv dbinom{n}{m} pmod{p_k^{c_k}}end{cases}$$

我们的问题就变成了求解这个同余方程组。利用中国剩余定理，我们知道最后解的形式一定是(xequiv 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_icdot 2p_icdot 3p_i cdots sp_i) prodlimits_{1leq k leq x , p_i mid k}k)
(=p_i^s(s!)prodlimits_{1leq k leq x , p_i mid k}k)

由于模的循环性质，有：(prodlimits_{1leq k leq p_i^{c_i}}kequiv prodlimits_{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}} floor)个。剩下的(x-lfloor frac{x}{p_i^{c_i}} floor)就直接暴力计算就可以了。这有一点点分块的味道。

对于每一个循环节，我们只需要暴力计算其中一个(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)的个数( ext{d}_p(x))。我们可以根据上面的过程，得到一个递归公式：

[ ext{d}_p(x)= ext{d}_p(lfloor frac{x}{p} floor)+ lfloor frac{x}{p} floor ag{1} ]

展开这个递归式可以得到( ext{d}_p(x))的计算公式：

[ ext{d}_p(x)=sum_{p^k leq x} lfloor frac{x}{p^k} floor ag{1'} ]