1.1 问题引入
已知(p)是一质数,求(dbinom{n}{m}pmod{p}).
关于组合数,它和排列数都是组合数学中的重要概念.这里会张贴有关这两个数的部分内容.
由于Lucas定理比较简单粗暴,这里直接给出表达式:
其中(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}}xi$$.
由上面的推论,在模(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))的计算公式：