P4720 【模板】扩展卢卡斯

文章目录

R e s u l t Result Result
H y p e r l i n k Hyperlink Hyperlink
D e s c r i p t i o n Description Description
S o l u t i o n Solution Solution
C o d e Code Code

R e s u l t Result Result

H y p e r l i n k Hyperlink Hyperlink

https://www.luogu.com.cn/problem/P4720

D e s c r i p t i o n Description Description

求 C n m m o d p C_n^m\mod p Cnmmodp

数据范围： n , m ≤ 1 0 18 , p ≤ 1 0 6 n,m\leq 10^{18},p\leq 10^6 n,m≤1018,p≤106

S o l u t i o n Solution Solution

由于 p p p不是质数，我们考虑将其转化为若干质数的乘积，即用唯一分解定理

得 p = ∏ p i c i p=\prod p_i^{c_i} p=∏pici

若如果我们能求出 C n m m o d p i c i C_n^m\mod p_i^{c_i} Cnmmodpici的值，再用 C R T CRT CRT合并，就可以得到答案

有 C n m = n ! m ! ( n − m ) ! C_n^m=\frac{n!}{m!(n-m)!} Cnm=m!(n−m)!n!
要求 n ! m ! ( n − m ) ! m o d p i c i \frac{n!}{m!(n-m)!}\mod p_i^{c_i} m!(n−m)!n!modpici
考虑把其中与 p i p_i pi有关的提出来，设 g ( n ) g(n) g(n)表示 n ! n! n!分解质因数后，质数 p i p_i pi的次数
得到 n ! g ( n ) m ! g ( m ) ( n − m ) ! g ( n − m ) p i g ( n ) − g ( m ) − g ( n − m ) \frac{\frac{n!}{g(n)}}{\frac{m!}{g(m)}\frac{(n-m)!}{g(n-m)}}p_i^{g(n)-g(m)-g(n-m)} g(m)m!g(n−m)(n−m)!g(n)n!pig(n)−g(m)−g(n−m)
用 E x g c d Exgcd Exgcd处理下面那一块，问题转换成为求 n ! p g ( n ) m o d p i c i \frac {n!}{p^{g(n)}}\mod p_i^{c_i} pg(n)n!modpici，显然这是有规律的

例如求 22 ! m o d 3 2 22!\mod 3^2 22!mod32
22 ! = 1 × 2 × 3 × 4 × 5 × 6 × 7 × 8 × 9 × 10 × 11 × 12 × 13 × 14 × 15 × 16 × 17 × 18 × 19 × 20 × 21 × 22 22!=1\times 2\times 3\times 4\times 5\times 6\times 7\times 8\times 9\times 10\times 11\times 12\times 13\times 14\times 15\times 16\times 17\times 18\times 19\times 20\times 21\times 22 22!=1×2×3×4×5×6×7×8×9×10×11×12×13×14×15×16×17×18×19×20×21×22

考虑提取出 3 3 3的倍数
22 ! = ( 3 × 6 × 9 × 12 × 15 × 18 × 21 ) × 22!=(3\times 6\times 9\times 12\times 15\times 18\times 21)\times 22!=(3×6×9×12×15×18×21)×剩下的
22 ! = 3 7 7 ! × 22!=3^77!\times 22!=377!×剩下的
22 ! = 3 7 7 ! × 1 × 2 × 4 × 5 × 7 × 8 × 10 × 11 × 13 × 14 × 16 × 17 × 19 × 20 × 22 22!=3^77!\times 1\times 2\times 4\times 5\times 7\times 8\times 10\times 11\times 13\times 14\times 16\times 17\times 19\times 20\times 22 22!=377!×1×2×4×5×7×8×10×11×13×14×16×17×19×20×22
在取模的意义下，运用结合律
22 ! = 3 7 7 ! × ( 1 × 2 × 4 × 5 × 7 × 8 ) × ( 10 × 11 × 13 × 14 × 16 ) × 17 × 19 × 20 × 22 22!=3^77!\times (1\times 2\times 4\times 5\times 7\times 8)\times (10\times 11\times 13\times 14\times 16)\times 17\times 19\times 20\times 22 22!=377!×(1×2×4×5×7×8)×(10×11×13×14×16)×17×19×20×22
括号那两坨在取模意义下是相等的
即
22 ! ≡ 3 7 7 ! × ( 1 × 2 × 4 × 5 × 7 × 8 ) 2 × 19 × 20 × 22 ( m o d p i c i ) 22!\equiv3^77!\times (1\times 2\times 4\times 5\times 7\times 8)^2\times 19\times 20\times 22(\mod p_i^{c_i}) 22!≡377!×(1×2×4×5×7×8)2×19×20×22(modpici)

考虑将其推演到 n ! n! n!，令 p k = p i c i pk=p_i^{c_i} pk=pici
n ! ≡ p i ⌊ n p i ⌋ ⌊ n p i ⌋ ! × ( ∏ i = 1 , i 不是 p i 的倍数 p k i ) ⌊ n p k ⌋ × ( ∏ i = p k × ⌊ n p k ⌋ n i ) ( m o d p k ) n!\equiv p_i^{\lfloor \frac n{p_i}\rfloor}\lfloor \frac n{p_i}\rfloor!\times (\prod_{i=1,i不是p_i的倍数}^{pk} i)^{\lfloor \frac n{pk}\rfloor}\times(\prod _{i=pk\times \lfloor \frac n{pk}\rfloor}^n i)(\mod pk) n!≡pi⌊pin⌋⌊pin⌋!×(∏i=1,i不是pi的倍数pki)⌊pkn⌋×(∏i=pk×⌊pkn⌋ni)(modpk)

而我们要求的是 n ! g ( n ) m o d p k \frac{n!}{g(n)}\mod pk g(n)n!modpk，也就是还要除一个 g ( n ) g(n) g(n)
所以我们令 f ( n ) = n ! g ( n ) m o d p k f(n)=\frac {n!}{g(n)}\mod pk f(n)=g(n)n!modpk，带入到上面那个式子，就可以得到

f ( n ) = f ( ⌊ n p i ⌋ ) × ( ∏ i = 1 , i 不是 p i 的倍数 p k i ) ⌊ n p k ⌋ × ( ∏ i = p k × ⌊ n p k ⌋ n i ) ( m o d p k ) f(n)=f(\lfloor \frac {n}{p_i}\rfloor)\times (\prod_{i=1,i不是p_i的倍数}^{pk} i)^{\lfloor \frac n{pk}\rfloor}\times(\prod _{i=pk\times \lfloor \frac n{pk}\rfloor}^n i)(\mod pk) f(n)=f(⌊pin⌋)×(∏i=1,i不是pi的倍数pki)⌊pkn⌋×(∏i=pk×⌊pkn⌋ni)(modpk)

g ( n ) g(n) g(n)的递推式显而易见，为 g ( n ) = ⌊ n p i ⌋ + g ( ⌊ n p i ⌋ ) g(n)=\lfloor \frac n{p_i}\rfloor+g(\lfloor \frac n{p_i}\rfloor) g(n)=⌊pin⌋+g(⌊pin⌋)（除去那个因子，递归计算）
而我们每次递归都会除掉 p i p_i pi的若干次方，递归的方式也和上式一致，所以最终除去的部分就是 p i g ( n ) p_i^{g(n)} pig(n)符合题意

好了，现在我们已经可以求出 C n m ≡ f ( n ) f ( n ) ( f ( n − m ) p i g ( n ) − g ( m ) − ( g ( n − m ) ( m o d p i c i ) C_n^m\equiv\frac{f(n)}{f(n)(f(n-m)}p_i^{g(n)-g(m)-(g(n-m)}(\mod p_i^{c_i}) Cnm≡f(n)(f(n−m)f(n)pig(n)−g(m)−(g(n−m)(modpici)了
考虑用 C R T CRT CRT将每个质因子求出来的答案合并

先写出 C R T CRT CRT的一般形式
{ ( d + x ) ≡ b 1 ( m o d a 1 ) ( d + x ) ≡ b 2 ( m o d a 2 ) ( d + x ) ≡ b 3 ( m o d a 3 ) … … \left\{\begin{matrix} (d+x)\equiv b_1 & (mod\ a_1)\\ (d+x)\equiv b_2 & (mod\ a_2)\\ (d+x)\equiv b_3 & (mod\ a_3)\\ …… \end{matrix}\right. ⎩⎪⎪⎨⎪⎪⎧(d+x)≡b1(d+x)≡b2(d+x)≡b3……(mod a1)(mod a2)(mod a3)

它能求出上述同余方程组 x x x的最小正整数解
具体做法为
令
M = ∏ a M=\prod a M=∏a
M i = M a i M_i=\frac {M}{a_i} Mi=aiM
M i ‘ ≡ M i − 1 ( m o d a i ) M_i`\equiv M_i^{-1}(mod\ a_i) Mi‘≡Mi−1(mod ai)
则
x + d ≡ ∑ ( b i × M i × M i ‘ ) ( m o d M ) x+d\equiv \sum (b_i\times M_i\times M_i^`)(\mod M) x+d≡∑(bi×Mi×Mi‘)(modM)

套用到这题来，首先我们有
{ C n m ≡ b 1 ( m o d p 1 c 1 ) C n m ≡ b 2 ( m o d p 2 c 2 ) C n m ≡ b 3 ( m o d p 3 c 3 ) … … \left\{\begin{matrix} C_n^m\equiv b_1 & (mod\ p_1^{c_1})\\ C_n^m\equiv b_2 & (mod\ p_2^{c_2})\\ C_n^m\equiv b_3 & (mod\ p_3^{c_3})\\ …… \end{matrix}\right. ⎩⎪⎪⎨⎪⎪⎧Cnm≡b1Cnm≡b2Cnm≡b3……(mod p1c1)(mod p2c2)(mod p3c3)
其中， b i b_i bi就是我们上面费尽心思求的那坨东西， M = ∏ p i c i = p M=\prod p_i^{c_i}=p M=∏pici=p
套用 C R T CRT CRT，即可得到 C n m C_n^m Cnm

时间复杂度： O ( p l o g p n ) O(plog_pn) O(plogpn)

C o d e Code Code

#include<cctype>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define LL long long
using namespace std;LL n,m,p;
inline LL read()
{char c;LL d=1,f=0;while(c=getchar(),!isdigit(c)) if(c=='-') d=-1;f=(f<<3)+(f<<1)+c-48;while(c=getchar(),isdigit(c)) f=(f<<3)+(f<<1)+c-48;return d*f;
}
inline LL ksm(LL x,LL y,LL p)
{LL res=1;for(;y;y>>=1,(x*=x)%=p) if(y&1) (res*=x)%=p;return res;
}
inline void Exgcd(LL a,LL b,LL &x,LL &y)
{if(b==0) return (void)(x=1,y=0);Exgcd(b,a%b,x,y);LL tmp=x;x=y;y=tmp-a/b*y;return;
}
inline LL INV(LL a,LL p)
{LL x,y;Exgcd(a,p,x,y);return (x+p)%p;
}
inline LL Fac(LL n,LL p,LL pk)
{if(n<2) return 1;LL prod_xh=1,prod_last=1;for(LL i=1;i<=pk;i++) if(i%p) (prod_xh*=i)%=pk;prod_xh=ksm(prod_xh,n/pk,pk);for(LL i=pk*(n/pk);i<=n;i++) if(i%p) (prod_last*=i%pk)%=pk;return Fac(n/p,p,pk)*prod_xh%pk*prod_last%pk;
}
inline LL G(LL n,LL p)
{if(n<p) return 0;return G(n/p,p)+(n/p);
}
inline LL C(LL n,LL m,LL p,LL pk){return Fac(n,p,pk)*INV(Fac(m,p,pk),pk)%pk*INV(Fac(n-m,p,pk),pk)%pk*ksm(p,G(n,p)-G(m,p)-G(n-m,p),pk)%pk;}
LL A[65],B[65],tot;
inline LL Exlucas(LL n,LL m,LL p)
{LL x=p;for(LL i=2;i*i<=x;i++)if(x%i==0){LL tmp=1;do{x/=i;tmp*=i;}while(x%i==0);B[++tot]=C(n,m,i,tmp);A[tot]=tmp;}if(x!=1){B[++tot]=C(n,m,x,x);A[tot]=x;}LL res=0,M=p;for(register int i=1;i<=tot;i++){LL Mi=M/A[i],Mi_=INV(Mi,A[i]);(res+=B[i]*Mi%M*Mi_%M)%=M;}return res;
}
signed main()
{n=read();m=read();p=read();printf("%lld",Exlucas(n,m,p));
}