扩展卢卡斯学习笔记

前言

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

扩展卢卡斯的作用就是在不保证(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 imes p^y)的形式，那么(C_n^m)也就可以表示成(x imes p^y)的形式，就很容易求出(C_n^m\%p^k)的值了。

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

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

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

对于剩余部分，由于模(p^k)意义下的值存在一个(p^k)的周期，因此我们暴力计算([1,p^k))中不含(p)的一般项乘积，计算(lfloorfrac n{p^k} floor)次只需给它做个快速幂，最后剩余一段不完整的直接暴力计算([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
",ExLucas::C(n,m,p)),0;
}