[筆記] 擴展Lucas定理
[筆記] 擴展\(Lucas\)定理
\(Lucas\)定理:\(\binom{n}{m} \equiv \binom{n/P}{m/P} \binom{n \% P}{m \% P}\pmod{P}\)\((P\ is \ prime)\)
Theory
那麽如果\(p\)不是一個質數怎麽辦?
當我們需要計算\(C_n^m\mod p\)，其中\(p = p_1^{q_1}\times p_2^{q_2}\times ...\times p_k^{q_k}\)，我們可以求出：\(C_n^m\equiv a_i\pmod {p_i^{q_i}} (1\lt i \lt k)\)
然後對於方程組：
\(x\equiv a_i \pmod {p_i^{q_i}}(1\lt i\lt k)\)
我們可以求出滿足條件的最小的\(x\),記為\(x_0\)那麽我們有： \(C_n^m\equiv x_0\pmod p\)
但是，我們發現，\(p_i^{q_i}\)並不是一個素數，它是某個素數的某次方。
下面我們介紹如何計算\(C_n^m \mod p^t(t\ge2,p \ is \ prime)\)
我們知道，\(C_n^m=\frac {n!}{m!(n-m)!}\)，若我們可以計算出\(m!\mod p^t\)，我們就能計算出\((n-m)!\mod p^t\)以及\((n-m)!\mod p^t\)
我們不妨設\(x=n!\mod p^t,y=m!\mod p^t,z=(n-m)!\mod p^t,\)
那麽答案就是
\(x\cdot inv(y,p^t)\cdot inv(z,p^t)\)那麽下面問題就轉化成如何計算\(n!\mod p^t\).
例如\(p=3,t=2,n=19\)
\(n!=1\times2\times3\times4\times5\times6\times7\times8\times ...\times19 \\ \ \ \ =(1\times2\times4\times5\times7\times8\times...\times16\times17\times19)\times(3\times6\times9\times12\times15\times18)\\\ \ \ =(1\times2\times4\times5\times7\times8\times...\times16\times17\times19)\times3^6\times(1\times2\times3\times4\times5\times6)\)
部分恰好是\((n/p)!\)，於是遞歸即可。前半部分是以\(p^t\)為周期的\((1\times2\times4\times5\times7\times8)\equiv(10\times11\times13\times14\times16\times17)\pmod9\).下面是孤立的\(19\)，可以知道孤立出來的長度不超過\(p^t\),於是直接計算即可。對於最後剩下的\(3^6\)這些數我們只要計算出\(n!,m!,(n-m)!\)裏含有多少個\(p\)(不妨設\(x,y,z\))，那麽\(x?y?z\)就是\(C_n^m\)中\(p\)的個數，直接計算就行。
Code
// luogu-judger-enable-o2
#include <iostream>
#include <cstdio>
#include <cmath>
typedef long long ll;
const ll N = 1e6 + 5;
ll g_x, g_y, cnt;
ll d[N], r[N];
ll q_pow(ll a, ll b, ll p){
ll w = 1;
while(b){
if(b & 1)
w = (a * w) % p;
b >>= 1;//1
a = (a * a) % p;
}
return w % p;
}
ll exgcd(ll a, ll b){
if(b == 0){
g_x = 1;
g_y = 0;
return a;
}
ll gcd = exgcd(b, a % b);
ll t = g_x;
g_x = g_y;
g_y = t - a / b * g_y;
return gcd;
}/*exgcd求逆元 , ab + mt = 1, 前提: gcd(a, m) = 1;
用exgcd求逆元有個好處,不用讓模數為素數,只要模數和這個a互質就好*/
ll inv(ll a, ll p){
exgcd(a, p);
return (g_x + p) % p;
}
ll fac(ll n, ll pi, ll pk){
if(!n) return 1; //2 //遞歸邊界
ll res = 1;
for(register ll i = 2; i <= pk; ++i){
if(i % pi)//3 pk
res = (res * i) % pk;
}/*因為循環節長度最多為pk,所以只需要算一遍pk,選這些數字裏面不是pi倍數的數字
(因為是倍數的,我們已經處理掉了*/
res = q_pow(res, n / pk, pk);//有n/pk個循環
for(register ll i = 2; i <= n % pk; ++i){
if(i % pi)//3 pk
res = (res * i) % pk;//剩下的暴力做
}
return (res * fac(n / pi, pi, pk)) % pk; //遞歸繼續
}
void cal(ll n, ll m, ll pi, ll pk){
ll up = fac(n, pi, pk), d1 = fac(n - m, pi, pk), d2 = fac(m, pi, pk);
ll k = 0;
for(register ll i = n; i; i /= pi) k += i / pi;
for(register ll i = m; i; i /= pi) k -= i / pi;
for(register ll i = n - m; i; i /= pi) k -= i / pi;//這三行都是統計pi倍數的個數,後兩個要減去是因為組合數的計算不就有個除嘛
r[++cnt] = q_pow(pi, k, pk) % pk * up % pk * inv(d1, pk) % pk * inv(d2, pk) % pk;
d[cnt] = pk;//CRT,r表示余數,d表示除數
}
ll mul(ll a, ll b, ll p){
ll f = 1;
if(a < 0) f = -f, a = -a;
if(b < 0) f = -f, b = -b;
ll w = 0;
while(b){
if(b & 1)
w = (w + a) % p;
b >>= 1;
a = (a + a) % p;
}
return w * f;
}//龜速乘
ll exCRT(){
for(register ll i = 2; i <= cnt; ++i){
ll C = r[1] - r[i];
ll D = exgcd(d[i], d[1]);
if(C % D) return -1;
ll k1 = mul(g_y, C / D, d[1] / D * d[i]);
ll x0 = mul(-k1 , d[1], d[1] / D * d[i] ) + r[1];
d[1] = d[1] / D * d[i], r[1] = x0;
r[1] = (r[1] % d[1] + d[1]) % d[1];
}
return r[1];//很好的exCRT
}
ll exlucas(ll n, ll m, ll p){
ll lim = sqrt(p) + 1;
ll tmp = p;
ll pk;
for(register ll i = 2; i <= lim; ++i){
if(tmp % i == 0){
pk = 1;
while(tmp % i == 0){
pk *= i;
tmp /= i;//為了得出pk
}
cal(n, m, i, pk);
}
}//唯一分解
if(tmp > 1) cal(n, m, tmp, tmp);//4 唯一分解後可能會留下一個大素數
return exCRT() % p;//5 pk
}
int main(){
ll n, m, p;
scanf("%lld %lld %lld", &n, &m, &p);
printf("%lld\n", exlucas(n, m, p));
}Wrong
- 快速冪中的指數忘記右移
- \(fac\)遞歸邊界忽略掉了
- 循環節計算時,是取那些不是\(pi\)倍數的,而不是不是\(pk\)倍數的
- \(exCRT\)最後模的是\(p\),不是\(pk\)!!!
愛你喲
[筆記] 擴展Lucas定理