卢卡斯定理:
( n m ) ≡ ( n % p n % p ) ( ⌊ n p ⌋ ⌊ m p ⌋ ) ( m o d p ) ) \boxed{\dbinom n m\equiv \dbinom{n\%p}{n\%p} \dbinom{\lfloor \dfrac n p\rfloor }{\lfloor \dfrac m p\rfloor }(\mod p))} (mn)≡(n%pn%p)(⌊pm⌋⌊pn⌋)(modp))
证明:二项式定理+费马小定理.
{ 1 + x p ≡ 1 + x ( 1 + x ) p ≡ 1 + x → ( 1 + x ) p ≡ 1 + x p ( m o d p ) → ( 1 + x ) n ≡ ( 1 + x ) ⌊ n p ⌋ p ( 1 + x p ) n % p ≡ ( 1 + x p ) ⌊ n p ⌋ ( 1 + x n % p ) ( m o d p ) \begin{cases} 1+x^p\equiv 1+x\\(1+x)^p \equiv 1+x\end{cases} \rightarrow (1+x)^p \equiv 1+x^p(\mod p)\rightarrow (1+x)^n\equiv(1+x)^{\lfloor \dfrac n p \rfloor p}(1 +x^p)^{n\%p}\equiv (1+x^p)^{\lfloor \dfrac n p \rfloor }(1+x^{n\%p})(\mod p) {1+xp≡1+x(1+x)p≡1+x→(1+x)p≡1+xp(modp)→(1+x)n≡(1+x)⌊pn⌋p(1+xp)n%p≡(1+xp)⌊pn⌋(1+xn%p)(modp).(费马小定理)
然后我们带入二项式定理,考究系数的关系.
为了方便,以后的下去整忽略不写.
∑ i = 0 n C n i x i ≡ ∑ i = 0 n / p C n / p i C n i x i p ∗ ∑ j = 0 n % p C n % p j x j ( m o d p ) \sum_{i=0}^n C_n^i x^i\equiv \sum_{i=0}^{ n/p} C_{n/p}^i C_n^i x^{ip} *\sum_{j=0}^{n\% p} C_{n\% p}^j x^j(\mod p) ∑i=0nCnixi≡∑i=0n/pCn/piCnixip∗∑j=0n%pCn%pjxj(modp)
∀ i ∈ [ 0 , n ] ∩ N , i = s p + t ( p , t ∈ N ) , C n i x i ≡ C n / p s x s p C n % p t x t \forall i\in [0,n]\cap \N,i=sp+t(p,t\in \N),C_n^ix^i\equiv C_{n/p}^s x^{sp}C_{n\%p}^tx^t ∀i∈[0,n]∩N,i=sp+t(p,t∈N),Cnixi≡Cn/psxspCn%ptxt
消去 x x x,则可得到卢卡斯定理.
实现:递归处理.复杂度: O ( log p n ) O(\log_p n) O(logpn).(不计预处理阶乘的时间)
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=1e5+10;
ll jc[N],jc_inv[N],ans;
ll power_mod(ll a,ll b,ll c)
{
ll ans=1%c;a%=c;
while(b>0)
{
if(b&1)ans=ans*a%c;
a=a*a%c;b=b>>1;
}
return ans;
}
int t,p,n,m;
void dfs(int a,int b)
{
if(!b)return;
dfs(a/p,b/p);
ans=ans*(jc[a%p]*jc_inv[a%p-b%p]*jc_inv[b%p]%p)%p;
}
int main()
{
scanf("%d",&t);
jc[0]=jc_inv[0]=1;
while(t--)
{
scanf("%d%d%d",&n,&m,&p);
ans=1;
for(int i=1;i<=p;i++)
{
jc[i]=jc[i-1]*i%p;
jc_inv[i]=power_mod(jc[i],p-2,p);
}
dfs(n+m,m);
printf("%lld\n",ans);
}
return 0;
}
拓展卢卡斯:
已知 n , m , p n,m,p n,m,p,求 ( n m ) m o d p \boxed{\dbinom n m \mod p} (mn)modp(其中 p p p不为质数)
推导:设 p = ∏ p i k i p=\prod p_i^{k_i} p=∏piki,则我们只要求解所有 ( n m ) m o d p i k i \dbinom n m \mod p_i^{k_i} (mn)modpiki,最后用中国剩余定理合并即可.
此时我们就求 n ! m ! ( n − m ) ! m o d p k \dfrac {n!}{m!(n-m)!}\mod p^k m!(n−m)!n!modpk(为了方便这里的p是质数,与上面的意义不同).
由于不一定有逆元,所以我们可以先把 p p p的倍数约去,即 n ! p x m ! p y ( n − m ) ! p z ∗ p x − y − z m o d p k \dfrac {\dfrac {n!}{p^x}}{\dfrac{m!}{p^y} \dfrac {(n-m)!} {p^z} }*p^{x-y-z}\mod p^k pym!pz(n−m)!pxn!∗px−y−zmodpk.
定义 F ( x ) F(x) F(x)表示 x ! x! x!除去 p p p的因子后 m o d p k \mod p^k modpk的结果.则有:
F ( x ) = p n / p F ( n / p ) ∏ i = 1 n / p k ∗ p k i ( i m o d p ≠ 0 ) ∏ i = n / p k ∗ p k + 1 n i ( i m o d p ≠ 0 ) F(x)=p^{n/p}F(n/p) \prod_{i=1}^{n/p^k*p^k} i(i\mod p\ne 0)\prod_{i=n/p^k*p^k+1}^n i(i\mod p\ne 0) F(x)=pn/pF(n/p)∏i=1n/pk∗pki(imodp=0)∏i=n/pk∗pk+1ni(imodp=0).
PS:前面部分表示含 p p p的因子的情况.
中间部分有周期性,我们可以减小枚举范围.
而且 p p p要消去,所以可以忽略.
所以 F ( x ) = F ( n / p ) ∗ ( ∏ i = 1 p k i ( i m o d p ≠ 0 ) ) n / p k ∗ ∏ i = 1 n − n / p k ∗ p k i ( i m o d p ≠ 0 ) F(x)=F(n/p)*(\prod_{i=1}^{p^k} i(i\mod p\ne 0))^{n/p^k}*\prod_{i=1}^{n-n/p^k*p^k} i(i\mod p\ne 0) F(x)=F(n/p)∗(∏i=1pki(imodp=0))n/pk∗∏i=1n−n/pk∗pki(imodp=0)
递归处理即可.递归层数 log p n \log_p n logpn,单层扫描 p k p^k pk.
总复杂度: O ( log p n p k ) O(\log_p^n p^k) O(logpnpk).
之后我们定义 G ( x ) G(x) G(x)表示 x x x中含有多少个 p p p个因子,复杂度: O ( log p n ) O(\log_p n) O(logpn).
单次求解 C n m m o d p i k i C_n^m \mod p_i^{k_i} Cnmmodpiki的复杂度: O ( p k log p n ) O(p^k \log_p^n) O(pklogpn).
合并复杂度: O ( log p k T ) , T O(\log p^k T),T O(logpkT),T为质因子个数.
总复杂度可以认为是: O ( p log p ) O(p\log p) O(plogp).
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int N=66;
typedef long long ll;
ll power(ll a,ll b,ll mod) {
ll c=1;
while(b&&c) {
if(b&1) c=c*a%mod;
b /= 2; a=a*a%mod;
}
return c;
}
void exgcd(ll a,ll b,ll &x,ll &y) {
if(!a) {x=0; y=1; return ;}
exgcd(b%a,a,y,x); x-=b/a*y;
}
ll inv(ll a,ll p) {
ll x,y; exgcd(a,p,x,y);
return (x%p+p)%p;
}
ll F(ll n,ll p,ll pk) {
if(!n) return 1;
ll s=1;
for(int i=1;i<pk;i++)
if(i%p) s=s*i%pk;
s=power(s,n/pk,pk);
for(int i=n%pk; i;i--)
if(i%p) s=s*i%pk;
return F(n/p,p,pk)*s%pk;
}
ll G(ll n,ll p) {
ll s=0;
while(n) s+=(n/=p);
return s;
}
ll C(ll n,ll m,ll p,ll pk) {
ll a=F(n,p,pk),b=inv(F(m,p,pk),pk),c=inv(F(n-m,p,pk),pk),d=power(p,G(n,p)-G(m,p)-G(n-m,p),pk);
return a*b%pk*c%pk*d%pk;
}
ll A[N],B[N];
ll lucas(ll n,ll m,ll p) {
int cnt=0;ll x=p;
for(int i=2;i*i<=x;i++) if(x%i==0) {
ll y=1;
while(x%i==0) y*=i,x/=i;
A[++cnt]=y; B[cnt]=C(n,m,i,y);
}
if(x>1) A[++cnt]=x,B[cnt]=C(n,m,x,x);
ll ans=0;
for(int i=1;i<=cnt;i++)
(ans += B[i]*(p/A[i])%p*inv(p/A[i],A[i])%p) %= p;
return ans;
}
int main() {
ll n,m,p; scanf("%lld%lld%lld",&n,&m,&p);
printf("%lld\n",lucas(n,m,p)); return 0;
}