扩展卢卡斯学习笔记

前言

不知道为什么对这种名字前面带个“扩展”的算法一直抱有一种畏惧心理。。。

扩展卢卡斯的作用就是在不保证\(p\)是质数，求\(C_n^m\% p\)。

但由于复杂度还是对\(p\)有一定依赖，因此还是做不到任意模数。

质因数分解+\(CRT\)合并

考虑我们给\(p\)做一个质因数分解，拆成若干\(p_i^{k_i}\)乘积的形式。

由于这些\(p_i^{k_i}\)两两互质，以它们计算出的答案可以直接\(CRT\)合并。

那么现在问题就被转化成了如何求解\(C_n^m\%p_i^{k_i}\)。

模\(p_i^{k_i}\)意义下的阶乘

首先我们考虑\(C_n^m=\frac{n!}{m!(n-m)!}\)，因此实际上只要分别求出\(n!,m!,(n-m)!\)的\(x\times p^y\)的形式，那么\(C_n^m\)也就可以表示成\(x\times p^y\)的形式，就很容易求出\(C_n^m\%p^k\)的值了。

要求\(p^y\)是很简单的，众所周知这就是\(\sum_{i}\lfloor\frac n{p^i}\rfloor\)。

核心还是在于求\(x\)，即\(n!\)除去\(p\)之后剩余的部分。

显然​\(n\)以内所有的\(k\times p\)都是含\(p\)的特殊项，把它们拎出来全都除以\(p\)便得到了\(\lfloor\frac np\rfloor!\)，那么直接递归处理就好了。

对于剩余部分，由于模\(p^k\)意义下的值存在一个\(p^k\)的周期，因此我们暴力计算\([1,p^k)\)中不含\(p\)的一般项乘积，计算\(\lfloor\frac n{p^k}\rfloor\)次只需给它做个快速幂，最后剩余一段不完整的直接暴力计算\([1,n\%p^k)\)中不含\(p\)的一般项乘积即可。

代码：\(O(\sum p_i^{k_i}logn)\)

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define LL long long
using namespace std;
LL n,m;int p;
namespace ExLucas//扩展卢卡斯
{
	I int QP(RI x,LL y,CI X) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
	I void exgcd(CI x,CI y,int& a,int& b) {y?(exgcd(y,x%y,b,a),b-=x/y*a):(a=1,b=0);}
	I int Inv(CI x,CI X) {RI a,b;return exgcd(x,X,a,b),(a%X+X)%X;}//非质数模数求逆元
	I int BF(CI n,CI p,CI X) {RI t=1;for(RI i=1;i<=n;++i) i%p&&(t=1LL*t*i%X);return t;}//暴力求不含p的一般项乘积
	I int Fac(Con LL& n,CI p,CI X)//递归求阶乘（除去p）
	{
		return n?(1LL*QP(BF(X-1,p,X),n/X,X)*BF(n%X,p,X)%X*Fac(n/p,p,X)%X):1;//利用周期性计算一般项，递归计算特殊项
	}
	I LL Get(LL n,CI p) {LL t=0;W(n) t+=(n/=p);return t;}//求出p的指数
	I int Calc(Con LL& n,Con LL& m,CI p,CI k)//计算C(n,m)%(p^k)
	{
		LL w=Get(n,p)-Get(m,p)-Get(n-m,p);if(w>=k) return 0;RI X=QP(p,k-w,1e9);//计算组合数中p的指数，剩余部分模数为p^(k-w)
		return 1LL*Fac(n,p,X)*Inv(1LL*Fac(m,p,X)*Fac(n-m,p,X)%X,X)%X*QP(p,w,1e9);//计算组合数剩余部分的值，再乘上p^w
	}
	I void CRT(int& r1,int& p1,CI r2,CI p2)//中国剩余定理合并
	{
		RI k=1LL*((r2-r1)%p2+p2)*Inv(p1,p2)%p2;r1=k*p1+r1,p1*=p2;
	}
	I int C(Con LL& n,Con LL& m,RI p)
	{
		RI i,x,k,t=0,X=1;for(i=2;i*i<=p;++i) if(!(p%i))//枚举质因子p
			{x=1,k=0;W(!(p%i)) p/=i,x*=i,++k;CRT(t,X,Calc(n,m,i,k),x);}//求解p^k与原有答案合并
		return p^1&&(CRT(t,X,Calc(n,m,p,1),p),0),t;//可能残余一个大质因子
	}
}
int main()
{
	return scanf("%lld%lld%d",&n,&m,&p),printf("%d\n",ExLucas::C(n,m,p)),0;
}