exLucas 算法学习笔记

一、问题

给定 \(n,m,p\)，求 \(C_n^m\bmod p\)，\(n,m\le 10^{18},p\le 10^6\)。
\(p\) 现在不一定是质数了，该怎么办？

二、解法

首先，数论题一个常见的做法：如果模数不一定是质数，那就把模数拆成若干个质数的积，然后分别求解，最后用中国剩余定理求解。

1. 拆 \(p\)

我们可以把 \(p\) 拆成 \(a_1^{b_1}\times a_2^{b_2}\times…\times a_k^{b_k}\)，其中所有 \(a_i\) 都是质数。
然后对于每一个 \(a_i^{b_i}\)，分别求出 \(C_n^m \bmod a_i^{b_i}\)，最后再用中国剩余定理合并。

2. 求 \(C_n^m \bmod a^{b}\)

这等同于求 \(\frac{n!}{m!\times(n-m)!}\bmod a^b\)。
但是我们不一定能求出 \(m!\times (n-m)!\) 的逆元，因为 \(m!\times (n-m)!\) 不一定与 \(a\) 互质。
所以我们要手动让它们互质。
原式可以表示为 \(\frac{\frac{n!}{a^x}}{\frac{m!}{a^y}\times\frac{(n-m)!}{a^z}}\times a^{x-y-z}\bmod a^b,x=v_a(n!),y=v_a(m!),z=v_a((n-m)!)\)。
然后我们可以用 \(\operatorname{F}(n,a)\) 表示 \(\frac{n!}{a^{v_a(n!)}}\)，用 \(\operatorname{G}(n,a)\) 表示 \(v_a(n!)\)，那么原式就是 \(\frac{\operatorname{F}(n,a)}{\operatorname{F}(m,a)\times \operatorname{F}(n-m,a)}\times a^{\operatorname{G}(n,a)-\operatorname{G}(m,a)-\operatorname{G}(n-m,a)}\bmod a^b\)。

3. 求 \(\operatorname{G}(n,a)\)

这可以用到 Legendre 公式，也就是 \(a\) 是质数，\(v_a(n!)=\sum\limits_{i=1}\limits^{\infty}\lfloor\frac{n!}{a^i}\rfloor\)。
所以直接按照这个求即可，时间复杂度 \(O(\log n)\)。

4. 求 \(\operatorname{F}(n,a)\)

由于这个数可能很大，我们求的其实是 \(\operatorname{F}(n,a)\bmod A\)，其中 \(A\) 就是 \(a^b\)。
好像就是 \(n!\times \operatorname{inv}({a^{\operatorname{G}(n,a)}})\bmod A\)。但是这样我们要预处理出 \(1!\to n!\)，时间复杂度 \(O(n)\)，超时。
考虑用另一种方法计算。我们注意到 \(n!=1\times 2\times …\times (A-1)\times A\times (A+1)\times …\times (2A-1)\times 2A\times (2A+1)\times…\times (\lfloor\frac{n}{A}\rfloor\times A-1)\times(\lfloor\frac{n}{A}\rfloor\times A)\times (\lfloor\frac{n}{A}\rfloor\times A+1)\times…\times n\)。
由于 \(\operatorname{F}(n,a)\) 中 \(a\) 的倍数全部被去掉（注意 \(a\) 是个质数），所以我们假设 \(\operatorname{H}(s,a)=\frac{s}{a^{v_a(s)}}\)。可以发现如果 \(\gcd(a,t)=1,t\mid s\)，那么 \(\operatorname{H}(s,a)=t\times \operatorname{H}(\frac{s}{t},a)\)。
所以
\(\begin{align}\operatorname{F}(n,a)&\equiv\operatorname{H}(1\times 2\times …\times (A-1)\times 1\times (A+1)\times …\times (2A-1)\times 2\times (2A+1)\times…\times (\lfloor\frac{n}{A}\rfloor\times A-1)\times\lfloor\frac{n}{A}\rfloor\times (\lfloor\frac{n}{A}\rfloor\times A+1))\times …\times n,a)\pmod A\\ &\equiv\operatorname{H}((1\times 2\times …\times (A-1))\times 1\times (1\times 2\times …\times (A-1))\times 2\times (1\times…\times (A-1))\times\lfloor\frac{n}{A}\rfloor\times (1\times …\times (n\bmod A)),a)\pmod A\\ &\equiv\operatorname{H}(((A-1)!)^{\lfloor\frac{n}{A}\rfloor}\times (n\bmod A)!\times …\times \lfloor\frac{n}{A}\rfloor!,a)\pmod A\\ &\equiv\operatorname{H}(((A-1)!)^{\lfloor\frac{n}{A}\rfloor}\times (n\bmod A)!\times …\times F(\lfloor\frac{n}{A}\rfloor,a),a)\pmod A\\ &\equiv((A-1)!)^{\lfloor\frac{n}{A}\rfloor}\times (n\bmod A)!\times …\times F(\lfloor\frac{n}{A}\rfloor,a)\pmod A\\ \end{align}\)
递归调用 \(O(\log n)\) 层，提前预处理出 \([1\to (A-1)]!\bmod A\) 每一层 \(O(\log n)\)（瓶颈在于快速幂），所以总时间复杂度 \(O(\log^2 n)\)。

5. 总时间复杂度

假设 \(c\) 是 \(p\) 的不同质因子个数，那么由于 \(2\times 3\times 5\times 7\times 11\times 13\times 17=1021020>10^6\)，所以 \(c\le 6\)。

1. 质因数分解，时间复杂度是 \(O(\sqrt p)\)。
2. 对于 \(i\in[1,c]\) 预处理 \([1\to (a_i^{b_i}-1)]!\bmod a_i^{b_i}\)，总时间复杂度是 \(O(\sum\limits_{i=1}\limits^c{a_i^{b_i}})\)。
3. 求 \(\operatorname{F},\operatorname{G}\)，因为至多调用计算函数 \(c\) 次，所以这个部分的上界是 \(O(c\log^2 n)\)。
4. 用中国剩余定理合并，一共合并 \(c\) 次，一次时间复杂度 \(O(\log p)\)，总时间复杂度 \(O(c\log p)\)。

所以总时间复杂度是 \(O(\sqrt p+c\log^2 n+c\log p+\sum\limits_{i=1}\limits^c{a_i^{b_i}})\le O(p)\)，并且常数很小。

三、代码

变量名与上文不完全对应。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll t,n,m,p,c=0,s[1000010],a[1010],b[1010];
inline void exgcd(ll &x,ll &y,ll a,ll b){
	if(!b){x=1;y=0;return;}
	exgcd(y,x,b,a%b);y-=a/b*x;
}
inline ll ksm(ll a,ll b,ll p){
	ll r=1;
	while(b){if(b&1)r=r*a%p;a=a*a%p;b>>=1;}
	return r;
}
inline ll inv(ll n,ll p){
	ll x,y;
	exgcd(x,y,n,p);
	return (x%p+p)%p;
}
inline ll G(ll n,ll p){
	ll r=0;
	while(n)n/=p,r+=n;
	return r;
}
inline ll F(ll n,ll p,ll pk){
	if(!n)return 1;
	return F(n/p,p,pk)*ksm(s[pk-1],n/pk,pk)%pk*s[n%pk]%pk;
}
inline ll qwq(ll n,ll m,ll p,ll pk){
	s[0]=1;
	for(ll i=1;i<pk;i++)
		if(i%p)s[i]=s[i-1]*i%pk;
		else s[i]=s[i-1];
	ll nowf=F(n,p,pk)*inv(F(m,p,pk)*F(n-m,p,pk)%pk,pk)%pk;
	ll nowg=ksm(p,G(n,p)-G(m,p)-G(n-m,p),pk);
	return nowf*nowg%pk;
}
inline ll exlucas(ll x,ll y,ll p){
	ll r=0,P=p;
	for(ll i=2;i*i<=p;i++)
		if(p%i==0){
			a[++c]=1;
			while(p%i==0)p/=i,a[c]*=i;
			b[c]=qwq(x,y,i,a[c]);
		}
	if(p>1)a[++c]=p,b[c]=qwq(x,y,p,p);
	for(ll i=1,x,y;i<=c;i++)r=(r+inv(P/a[i],a[i])*P/a[i]%P*b[i]%P)%P;
	return r;
}
int main(){
	cin>>n>>m>>p;
	cout<<exlucas(n,m,p)<<'\n';
	return 0;
}