数论-卢卡斯定理（lucas)与拓展卢卡斯定理 (exlucas)

卢卡斯定理（lucas)

[用途]

求解 C n m % p C_{n}^{m}\% p Cnm%p,其中m,n较大，p较小且为素数

[结论]

C n m ≡ C n / p m / p C n % p m % p ( m o d p ) C_{n}^{m}\equiv C_{n/p}^{m/p}C_{n\%p}^{m\%p}(mod \quad p) Cnm≡Cn/pm/pCn%pm%p(modp)

[证明]

假设： { n = s p + q m = t p + r \left\{\begin{array}{cc} n=sp+q\\ m=tp+r \end{array}\right. {n=sp+qm=tp+r
即证明 C n m ≡ C s t C q r ( m o d p ) C_{n}^{m}\equiv C_{s}^{t}C_{q}^{r}(mod \quad p) Cnm≡CstCqr(modp)
使用构造法：
( 1 + x ) n ≡ ( 1 + x ) s p + q ≡ [ ( 1 + x ) p ] s ( 1 + x ) q ( m o d p ) (1+x)^{n}\equiv(1+x)^{sp+q}\equiv[(1+x)^{p}]^{s}(1+x)^{q}(mod \quad p) (1+x)n≡(1+x)sp+q≡[(1+x)p]s(1+x)q(modp)
根据二项式定理
( 1 + x ) p = ∑ i = 0 p C p i x i (1+x)^{p}=\sum\limits_{i=0}^p{C_{p}^{i}x^{i}} (1+x)p=i=0∑pCpixi
除了第一项 C p 0 x 0 （即 1 ） C_{p}^{0}x^{0}（即1） Cp0x0（即1）与最后一项 C p p x p （即 x p ) C_{p}^{p}x^{p}（即x^{p}) Cppxp（即xp)不能整除p，中间项均能整除p,所以在 ( m o d p ) (mod \quad p) (modp)意义下 ( 1 + x ) p = 1 + x p (1+x)^{p}=1+x^{p} (1+x)p=1+xp
则 ( 1 + x ) n ≡ [ ( 1 + x ) p ] s ( 1 + x ) q ≡ ( 1 + x p ) s ( 1 + x ) q ( m o d p ) (1+x)^{n}\equiv[(1+x)^{p}]^{s}(1+x)^{q}\equiv(1+x^{p})^{s}(1+x)^{q}(mod \quad p) (1+x)n≡[(1+x)p]s(1+x)q≡(1+xp)s(1+x)q(modp)
即 ∑ k = 0 n C n k x k ≡ ∑ i = 0 s C s i x p i ∑ j = 0 q C q j x j ( m o d p ) \sum\limits_{k=0}^n{C_{n}^{k}x^{k}}\equiv\sum\limits_{i=0}^s{C_{s}^{i}x^{pi}}\sum\limits_{j=0}^q{C_{q}^{j}x^{j}}(mod \quad p) k=0∑nCnkxk≡i=0∑sCsixpij=0∑qCqjxj(modp)
我们在等式两边找 x m x^m xm项的系数
左:含 x m x^m xm项显然是 C n m C_{n}^{m} Cnm
右:含 x m x^m xm项只能由 x p i x j = x p i + j x^{pi}x^{j}=x^{pi+j} xpixj=xpi+j得到，即 m = p i + j m=pi+j m=pi+j，显然 i = t , j = r i=t,j=r i=t,j=r，此时对应系数恰好是 C s t C q r C_{s}^{t}C_{q}^{r} CstCqr
所以 C n m ≡ C s t C q r ( m o d p ) C_{n}^{m}\equiv C_{s}^{t}C_{q}^{r}(mod \quad p) Cnm≡CstCqr(modp)成立
即 C n m ≡ C n / p m / p C n % p m % p ( m o d p ) C_{n}^{m}\equiv C_{n/p}^{m/p}C_{n\%p}^{m\%p}(mod \quad p) Cnm≡Cn/pm/pCn%pm%p(modp)

[代码]

#include<iostream>
using namespace std;
typedef long long LL;
const LL N=1e5+2;
LL a[N];
void init(LL p)
{a[1]=1;for(int i=2;i<=p;++i)a[i]=a[i-1]*i%p;
}
void exgcd(LL a,LL b,LL &x,LL &y)
{if(!b){x=1;y=0;return;}exgcd(b,a%b,y,x);y-=a/b*x;
}
LL ksm(LL x,LL n,LL mod)
{LL ans=1;while(n){if(n&1)ans=ans*x%mod;n>>=1;x=x*x%mod;}return ans;
}
LL C(LL n,LL m,LL p)
{if(n==m||m==0)return 1;if(n<m)return 0;if(m*2>n)m=n-m;                        /*C(n,m)=c(n,n-m)*/return a[n]*ksm(a[m]*a[n-m],p-2,p)%p; /*求(a[m]*a[n-m])在(mod p)意义下的乘法逆元*//*拓展欧几里得与费马小定理均可*/ /*LL x,y;exgcd(a[m]*a[n-m],p,x,y);return (a[n]*x%p+p)%p;*/
}
LL lucas(LL n,LL m,LL p)
{if(!m)return 1;return lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}
int main()
{ios::sync_with_stdio(false);LL T,n,m,p;cin>>T;while(T--){cin>>n>>m>>p;init(p);cout<<lucas(n+m,m,p)<<endl;}return 0;
}

拓展卢卡斯定理（exlucas)

[用途]

求解 C n m % P C_{n}^{m}\% P Cnm%P,其中m,n较大，P较小且不一定为素数

[解析]

转化为中国剩余定理(crt）

根据算术基本定理

P = M 1 M 2 . . . M k = ∏ i = 0 k M i P=M_{1}M_{2}...M_{k}=\prod\limits_{i=0}^{k}M_{i} P=M1M2...Mk=i=0∏kMi
其中 M i = p i t i , M i M j ( i ! = j ) M_{i}=p_{i}^{t_{i}},M_{i}M_{j}(i!=j) Mi=piti,MiMj(i!=j)两两互素
设 X = C n m X=C_{n}^{m} X=Cnm
则 A i = C n m % M i A_{i}=C_{n}^{m}\%M_{i} Ai=Cnm%Mi
所以转化为中国剩余定理
X ≡ { A 1 ( m o d M 1 ) A 2 ( m o d M 2 ) . . . A k ( m o d M k ) X\ \equiv\left\{\begin{array}{cc} A_{1} \quad (mod \quad M_{1})\\ A_{2} \quad (mod \quad M_{2})\\ ...\\ A_{k} \quad (mod \quad M_{k}) \end{array}\right. X ≡⎩⎪⎪⎨⎪⎪⎧A1(modM1)A2(modM2)...Ak(modMk)
M i M_{i} Mi根据算术基本定理循环分解即可
重点的是 A i A_{i} Ai
A i = C n m % M i = n ! m ! ( n − m ) ! % M i \begin{aligned}A_{i} &=C_{n}^{m}\%M_{i}\\ &=\frac{n!}{m!(n-m)!}\%M_{i}\\ \end{aligned} Ai=Cnm%Mi=m!(n−m)!n!%Mi
可惜的是 m ! ( n − m ) ! 在模 M i ( M i = p i t i ) m!(n-m)!在模M_{i}(M_i=p_{i}^{t_{i}}) m!(n−m)!在模Mi(Mi=piti)意义下的乘法逆元不一定存在
我们可以先把阶层中含 p i p_{i} pi的项约去，这样就可以求逆元了，然后再乘回来即可
所以
A i = n ! m ! ( n − m ) ! % M i = n ! p i m ! p i ( n − m ) ! p i p i k % M i = n ! p i i n v ( m ! p i , M i ) i n v ( ( n − m ) ! p i , M i ) p i k % M i = ( n ! p i % M i ) i n v ( m ! p i % M i , M i ) i n v ( ( n − m ) ! p i % M i , M i ) p i k % M i \begin{aligned}A_{i} &=\frac{n!}{m!(n-m)!}\%M_{i}\\ &=\frac{\frac{n!}{p_i}}{\frac{m!}{p_i}\frac{(n-m)!}{p_i}}p_i^k\%M_{i}\\ &=\frac{n!}{p_i}inv(\frac{m!}{p_i},M_i)inv(\frac{(n-m)!}{p_i},M_i)p_i^k\%M_{i}\\ &=(\frac{n!}{p_i}\%M_{i})inv(\frac{m!}{p_i}\%M_{i},M_i)inv(\frac{(n-m)!}{p_i}\%M_{i},M_i)p_i^k\%M_{i} \end{aligned} Ai=m!(n−m)!n!%Mi=pim!pi(n−m)!pin!pik%Mi=pin!inv(pim!,Mi)inv(pi(n−m)!,Mi)pik%Mi=(pin!%Mi)inv(pim!%Mi,Mi)inv(pi(n−m)!%Mi,Mi)pik%Mi

令
n n = n ! p i % M i nn=\frac{n!}{p_i}\%M_{i} nn=pin!%Mi
m m = m ! p i % M i mm=\frac{m!}{p_i}\%M_{i} mm=pim!%Mi
n m = ( n − m ) ! p i % M i nm=\frac{(n-m)!}{p_i}\%M_{i} nm=pi(n−m)!%Mi
则 A i = n n ∗ i n v ( m m , M i ) ∗ i n v ( n m , M i ) ∗ p i k % M i A_i=nn*inv(mm,M_i)*inv(nm,M_i)*p_i^k\%M_i Ai=nn∗inv(mm,Mi)∗inv(nm,Mi)∗pik%Mi
现在的关键又转换为求 n n , m m , n m nn,mm,nm nn,mm,nm了
就nn来看
当n=19时, p i = 3 , t i = 2 p_i=3,t_i=2 pi=3,ti=2
n ! = 1 ∗ . . . ∗ 19 = 3 6 ∗ 6 ! ∗ [ ( 1 ∗ 2 ∗ 4 ∗ 5 ∗ 7 ∗ 8 ) ∗ ( 10 ∗ 11 ∗ 13 ∗ 14 ∗ 16 ∗ 17 ) ∗ 19 ) ] n!=1*...*19=3^6*6!*[(1*2*4*5*7*8)*(10*11*13*14*16*17)*19)] n!=1∗...∗19=36∗6!∗[(1∗2∗4∗5∗7∗8)∗(10∗11∗13∗14∗16∗17)∗19)]
n ! p i % M i = [ ( 1 ∗ 2 ∗ 4 ∗ 5 ∗ 7 ∗ 8 ) ∗ ( 10 ∗ 11 ∗ 13 ∗ 14 ∗ 16 ∗ 17 ) ∗ 19 ) ] ∗ 6 ！ % M i \frac{n!}{p_i}\%M_{i}=[(1*2*4*5*7*8)*(10*11*13*14*16*17)*19)]*6！\%M_i pin!%Mi=[(1∗2∗4∗5∗7∗8)∗(10∗11∗13∗14∗16∗17)∗19)]∗6！%Mi
而中括号这一堆是有循环节的 1 ∗ 2 ∗ 4 ∗ 5 ∗ 7 ∗ 8 ≡ 10 ∗ 11 ∗ 13 ∗ 14 ∗ 16 ∗ 17 ( m o d M i ) 1*2*4*5*7*8\equiv10*11*13*14*16*17(mod\quad M_i) 1∗2∗4∗5∗7∗8≡10∗11∗13∗14∗16∗17(modMi)
循环节为 n / M i n/M_i n/Mi，循环一组后然后快速幂就行了
剩下的个数就是 n % M i n\%M_i n%Mi这一类直接循环相乘
而 6 ! 6! 6!恰好是 ( n / p i ) ! (n/p_i)! (n/pi)!,直接递归即可
即 n ! p i % M i = [ . . . ] ∗ ( n / p i ) ! p i % M i \frac{n!}{p_i}\%M_{i}=[...]*\frac{(n/p_i)!}{p_i}\%M_{i} pin!%Mi=[...]∗pi(n/pi)!%Mi
[ . . . ] [...] [...]就是上面的循环节那块，分类进行求解,over！

[代码]

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL N=1e5+9;
LL A[N],M[N];
LL ksm(LL x,LL n,LL mod)
{LL ans=1;while(n){if(n&1)ans=ans*x%mod;n>>=1,x=x*x%mod;}return ans;
}
void exgcd(LL a,LL b,LL &x,LL &y)
{if(!b)x=1,y=0;else exgcd(b,a%b,y,x),y-=a/b*x;
}
LL inv(LL a,LL p)
{LL x,y;exgcd(a,p,x,y);return (x+p)%p?x:x+p;
}
LL get(LL n,LL pi,LL p)                                 /*求(与pi互素后的n!)%M[i]*/
{if(!n)return 1;LL ans=1;if(n/p){                                          /*判断有无循环节 */ for(LL i=2;i<=p;++i)if(i%pi)ans=ans*i%p;ans=ksm(ans,n/p,p);}for(LL i=2;i<=n%p;++i)if(i%pi)ans=ans*i%p;            /*循环节剩余部分*/ return ans*get(n/pi,pi,p)%p;
}
LL exlucas(LL n,LL m,LL pi,LL p)                        /*求A[i]*/
{LL nn=get(n,pi,p);                                    /*求(与pi互素后的n)%M[i]*/ LL mm=get(m,pi,p);                                    /*求(m!与pi互素后的m!)%M[i]*/ LL nm=get(n-m,pi,p);                               /*求(与pi互素后的(n-m)!)%M[i]*/ LL k=0;                                              /*含质因数pi的数量*/ for(LL i=n;i;i/=pi)k+=i/pi;for(LL i=m;i;i/=pi)k-=i/pi;for(LL i=n-m;i;i/=pi)k-=i/pi;return nn*inv(mm,p)*inv(nm,p)*ksm(pi,k,p)%p;
}
LL crt(LL len,LL Lcm)
{LL ans=0;for(LL i=1;i<=len;++i){LL Mi=Lcm/M[i];ans=((ans+A[i]*inv(Mi,M[i])*Mi)%Lcm+Lcm)%Lcm;}return ans;
}
int main()
{ios::sync_with_stdio(false);LL n,m,P,num;while(cin>>n>>m>>P){if(n<m){cout<<0<<endl;continue;}num=0;memset(A,0,sizeof(A));memset(M,0,sizeof(M));for(LL x=P,i=2;i<=P;++i)if(x%i==0){M[++num]=1;while(x%i==0){M[num]*=i;x/=i;}A[num]=exlucas(n,m,i,M[num])%P;} cout<<crt(num,P)<<endl;}return 0;
}