[学习笔记][详解]扩展卢卡斯

扩展卢卡斯

upd:2020.10.15 重写

一.什么是扩展Lucas

解决组合数对非质数取模。

二.算法流程

由于P可能为合数，先将P分解质因数。分解为\(P=p_1^{k_1}*p_2^{k_2}...p_n^{k_n}\)对于任意两个\(p_i^{k_i}\),是互质的，并且有几个方程。形如

\[C_n^m\equiv C_n^m\%p_i^{k_i}(mod\;p_i^{k_i}) \]

可用中国剩余定理（CRT）将各个方程合并（不同质因子之间肯定互质，模数乘起来就会变作P）。此处挂一个算法主体与CRT

ll CRT(ll r,ll pi){return r*Inv(mod/pi,pi)%mod*(mod/pi)%mod;}
ll Exlucas(ll a,ll b){
	ll tmp=mod,res=0;
	for(ll i=2;i*i<=mod;++i)
		if(tmp%i==0){
			ll pk=1;
			while(!(tmp%i))tmp/=i,pk*=i;
			(res+=CRT(C(a,b,i,pk),pk))%=mod;
		}
	if(tmp>1)(res+=CRT(C(a,b,tmp,tmp),tmp))%=mod;
	return res%mod;
}

(mod为模数)

三.阶乘的本质

\[n!=p^{⌊\frac{n}{p}⌋}*⌊\frac{n}{p}⌋!* ∏_{p∤i}^{n}i \]

其中 p 可取任何值。（即使大于n也可以，想想为什么）

很好理解，就只是选了一个数 p ，将 \(1-n\) 中的 p 的倍数全部提取出来，然后再将系数提出来，分成系数相乘和 p 的幂次，剩下不是p的倍数的原样处理。

举个例子，

\[\begin{align} &1*2*3*4*5*6*7*8*9\\ =&(3*6*9)*(1*2*4*5*7*8)\\ =&(1*3*2*3*3*3)*(1*2*4*5*7*8)\\ =&(1*2*3)*3^3*(1*2*4*5*7*8) \end{align} \]

这个柿子具有可以递归的性质，这很重要。至于有什么具体的用处请见下文。

四.求解组合数模单个质数乘方的值

先把组合数的式子展开。

\[\frac{n!}{(n-m)!m!}\%p_i^{k_i} \]

分母目前无法用逆元处理，因为并不与 \({p_i}^{k_i}\) 互质。为了达到互质的目的，将所有的 \(p_i\) 因子提取出来，其中 \(f(n,p_i)\) 为 \(n!\) 中因子 \(p_i\) 的个数，于是。

\[\frac{n!}{(n-m)!m!}\equiv\frac{\frac{n!}{p_i^{f(n,p_i)}}}{\frac{m!}{p_i^{f(m,p_i)}}*\frac{(n-m)!}{p_i^{f(n-m,p_i)}}}*p_i^{f(n,p_i)-f(m,p_i)-f(n-m,p_i)} \pmod{p_i^{k_i}} \]

已经将多余的因子影响去掉了，可以逆元。此时需要一个求解\(x!~\%~p^k\)的函数，并在其中不乘上 \(p\) 的倍数。

挂一个C的写法。（其中jc指求解\(x!~\%~p^k\)的函数）

ll C(ll a,ll b,ll pi,ll pk){
	ll k=f(a,pi)-f(b,pi)-f(a-b,pi);
	if(!KSM(pi,k,pk))return 0;
	return jc(a,pi,pk)*Inv(jc(b,pi,pk),pk)%pk*Inv(jc(a-b,pi,pk),pk)%pk*KSM(pi,k,pk)%pk;
}

至于 f 既可以用递归，也可以非递归。递归形式为

\[f(n,p)=f(\lfloor\frac{n}{p}\rfloor,p)+\lfloor\frac{n}{p}\rfloor \]

（可以参见三大点进行理解，非递归放在第五大点中）

四.写出阶乘函数（去掉 p 的影响）

\[n!=p^{⌊\frac{n}{p}⌋}*⌊\frac{n}{p}⌋!* ∏_{p∤i}^{}n_i \]

首先柿子很明显有递归性质。含 \(p\) 的第一项舍去（不计p的影响），第二项递归。考虑第三项。

因为是对于 \({p_i}^{k_i}\) 进行取模，可以对于要乘的一系列连续的数进行划分，长度为 \({p_i}^{k_i}\) 。用柿子来解释（第一个结合律展开一下消掉含有\(p^k\)的柿子）。

\[\begin{align} \because\prod\limits_{i=u*{p}^{k}+1,p∤i}^{(u+1)*{p}^{k}}i\equiv\prod\limits_{i=u*{p}^{k}+1,p∤i}^{(u+1)*{p}^{k}}i+p^k \pmod{p^k} \\\therefore\prod\limits_{i=u*{p}^{k}+1,p∤i}^{(u+1)*{p}^{k}}i\equiv\prod\limits_{i=(u+1)*{p}^{k}+1,p∤i}^{(u+2)*{p}^{k}}i \pmod{p^k} \end{align} \]

划分之后只取一段然后快速幂即可。注意处理最后凑不齐一段的数。

ll jc(ll a,ll pi,ll pk){
	if(!a)return 1;
	ll res=1;
	m_for(i,2,pk)if(i%pi)res=(res*i)%pk;
	res=KSM(res,a/pk,pk);
	m_for(i,2,a%pk)if(i%pi)res=(res*i)%pk;
	return res*jc(a/pi,pi,pk)%pk;
}

五.Code

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define m_for(i,a,b)for(int i=a;i<=b;++i)
const int N=50010;
ll mod;
ll n,k;
ll f(ll a,ll pi){return a?f(a/pi,pi)+(a/pi):0;}
/*
ll f(ll a,ll pi){
	ll res=0;
	while(a)res+=a/pi,a/=pi;
	return res;
}
*/
ll Exgcd(ll a,ll b,ll &x,ll &y){
	if(!b)x=1,y=0;
	else Exgcd(b,a%b,y,x),y-=a/b*x;
}
ll Inv(int a,int p){
	ll x,y;Exgcd(a,p,x,y);
	return (x%p+p)%p;
}
inline ll KSM(ll a,ll b,ll p){
	ll res=1;
	for(;b;b>>=1,(a*=a)%=p)if(b&1)(res*=a)%p;
	return res%p;
}
ll jc(ll a,ll pi,ll pk){
	if(!a)return 1;
	ll res=1;
	m_for(i,2,pk)if(i%pi)res=(res*i)%pk;
	res=KSM(res,a/pk,pk);
	m_for(i,2,a%pk)if(i%pi)res=(res*i)%pk;
	return res*jc(a/pi,pi,pk)%pk;
}
ll C(ll a,ll b,ll pi,ll pk){
	ll k=f(a,pi)-f(b,pi)-f(a-b,pi);
	if(!KSM(pi,k,pk))return 0;
	return jc(a,pi,pk)*Inv(jc(b,pi,pk),pk)%pk*Inv(jc(a-b,pi,pk),pk)%pk*KSM(pi,k,pk)%pk;
}
ll CRT(ll r,ll pi){return r*Inv(mod/pi,pi)%mod*(mod/pi)%mod;}
ll Exlucas(ll a,ll b){
	ll tmp=mod,res=0;
	for(ll i=2;i*i<=mod;++i)
		if(tmp%i==0){
			ll pk=1;
			while(!(tmp%i))tmp/=i,pk*=i;
			(res+=CRT(C(a,b,i,pk),pk))%=mod;
		}
	if(tmp>1)(res+=CRT(C(a,b,tmp,tmp),tmp))%=mod;
	return res%mod;
}
int main(){
	scanf("%lld%lld%lld",&n,&k,&mod);
	printf("%lld",Exlucas(n,k));
	return 0;
}