exLucas
/* exLucas : 求 C(n,m)%P (P不保证是质数&&P<=1e6) 的问题 算法 设 P= p1^a1 * p2^a2 * ... pr^ar <-唯一分解定理 求出 { C(n,m) % p1^a1 C(n,m) % p2^a2 ... C(n,m) % pr^ar } 最后用中国剩余定理合并 分解成要求 C(n,m) % p^k(p为质数) n! so = ———————————— % p^k m!(n-m)! 由于式子是在模 p^k 意义下,所以分母要算乘法逆元. 我们无法求得m!,(n−m)!的逆元, why? m!可能含有 p 因子. 同余方程 a*x ≡ 1(%P)(即乘法逆元)有解 的充要条件为 <-裴蜀定理:gcd(a,b)=1 so我们换一个式子: n! ————— p^x —————————————— * p ^ (x-y-z) % p^k m! (n-m)! ——— * ———————— p^y p^z x 表示 n! 中包含多少个p因子,y,z同理. 现在问题转化为问题等价于求 n! ————— % p^k p^x 我们对n!进行变形: n!=1*2*3*...*n =(P*2P*3P...)(1*2*...) 左边是1∼n中所有 <=n的 P 的倍数,右边反之. 因为1∼n中有 ⌊Pn⌋ 个 P 的倍数, so n! = p^(⌊n/P⌋) * (1*2*3*...) * (1*2*3*...) ( n ) = p^(⌊n/P⌋) * (⌊n/P⌋)! *( ∏ i ) ( i=1,(i%p!=0) ) 显然后面这个 % P 是有循环节的. ( p^k )^⌊n/(p^k)⌋ ( n ) = p^(⌊n/P⌋) * (⌊n/P⌋)! * ( ∏ i ) ( ∏ i ) ( i=1,(i%p!=0) ) ( i=p^k*⌊n/(p^k)⌋,(i%p!=0) ) 其中 ( p^k )^⌊n/(p^k)⌋ ( ∏ i ) ( i=1,(i%p!=0) ) 是循环节1∼P中所有无 P 因子数的乘积. ( n ) ( ∏ i ) ( i=p^k*⌊n/(p^k)⌋,(i%p!=0) ) 是余项的乘积. 发现前面这个 P^⌊n/P⌋ 是要除掉的. 然而(⌊n/P⌋)!里可能还包含 P . 所以,我们定义函数: n! f(n) = ——————— P^x so ( p^k )^⌊n/(p^k)⌋ ( n ) f(n) = f(⌊n/P⌋) * (⌊n/P⌋) * ( ∏ i ) ( ∏ i ) ( i=1,(i%p!=0) ) ( i=p^k*⌊n/(p^k)⌋,(i%p!=0) ) 这样就好了.时间复杂度是O(logPn). 边界f(0)=1 看看原式 n! ————— p^x —————————————— * p ^ (x-y-z) % p^k m! (n-m)! ——— * ———————— p^y p^z f(n) ======= ——————————— * p ^ (x-y-z) % p^k f(m)*f(n-m) 由于f(m)一定与P^k互质,所以我们可以直接用exgcd求逆元. How to solve P ^ (x-y-z)? easy! For example : we want to get f(n) = n! / (p^x) 中的x 设g(n)=x(above) 再看看阶乘的式子 ( n ) n! = p^(⌊n/P⌋) * (⌊n/P⌋)! *( ∏ i ) ( i=1,(i%p!=0) ) 显然 we want p^(⌊n/P⌋) 而且(⌊nP⌋)!可能还有P因子。 所以我们可以得到以下递推式: g(n)=⌊n/P⌋+g(⌊n/P⌋) 时间复杂度是O(log(n/p)) 边界g(n)=0(n<P) 所以答案就是 n!/(p^x) = ———————————————————————————— * P^{g(n)-g(m)-g(n-m)} % (p^k) m!/(p^y) * (n-m)!/(p^z) 最后用中国剩余定理合并答案即可得到C(n,m) */ #include<bits/stdc++.h> #define Bessie moo~ #define int long long using namespace std; int read(){ int a=0,fl=1; char c; c=getchar(); while(c<'0'||c>'9')fl= c=='-' ? -1 : 1, c=getchar(); while('0'<=c&&c<='9')a=(a<<3)+(a<<1)+(c^48),c=getchar(); return a*fl; } void exgcd(int a,int b,int &x,int &y){ if(!b){ x=1,y=0; return ; } exgcd(b,a%b,x,y); int t=x; x=y,y=t-a/b*y; } int INV(int a,int MOD){ int x,y; exgcd(a,MOD,x,y); return (x+MOD)%MOD; } int qpow(int a,int b,int MOD){ int base=1; a%=MOD; while(b){ if(b&1)base=(base*a)%MOD; b>>=1; a=(a*a)%MOD; } return base; } int F(int n,int MOD,int PK){ if(n==0)return 1; int rou=1,//循环节 rem=1;//余项 for(int i=1;i<=PK;i++){ if(i%MOD)rou=rou*i%PK; } rou=qpow(rou,n/PK,PK); for(int i=PK*(n/PK);i<=n;i++){ if(i%MOD)rem=rem*(i%PK)%PK; } return F(n/MOD,MOD,PK)*rou%PK*rem%PK; } int G(int n,int MOD){ if(n<MOD)return 0; return G(n/MOD,MOD)+(n/MOD); } int C_PK(int n,int m,int MOD,int PK){ int fz=F(n,MOD,PK),fm1=INV(F(m,MOD,PK),PK),fm2=INV(F(n-m,MOD,PK),PK); int mi=qpow(MOD,G(n,MOD)-G(m,MOD)-G(n-m,MOD),PK); return fz*fm1%PK*fm2%PK*mi%PK; } int A[1005],B[1005]; //x≡B(MOD A) int exLucas(int n,int m,int MOD){ int ljc=MOD,tot=0; for(int tmp=2;tmp*tmp<=MOD;tmp++){ if(!(ljc%tmp)){ int PK=1;//p^k while(!(ljc%tmp)){ PK*=tmp; ljc/=tmp; } A[++tot]=PK; B[tot]=C_PK(n,m,tmp,PK); } } if(ljc!=1){ A[++tot]=ljc; B[tot]=C_PK(n,m,ljc,ljc); } int ans=0; for(int i=1;i<=tot;i++){//CRT int M=MOD/A[i],T=INV(M,A[i]); ans=(ans+B[i]*M%MOD*T%MOD)%MOD; } return ans; } signed main(){ int n=read(),m=read(),MOD=read(); printf("%lld",exLucas(n,m,MOD)); return 0; }
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 单线程的Redis速度为什么快?
· 展开说说关于C#中ORM框架的用法!
· Pantheons:用 TypeScript 打造主流大模型对话的一站式集成库
· SQL Server 2025 AI相关能力初探
· 为什么 退出登录 或 修改密码 无法使 token 失效