exLucas 算法学习笔记

一、问题

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

二、解法

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

1. 拆 $p$

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

2. 求 $C_n^{m}mod{a}^b$

这等同于求 $\frac{n!}{m!\times (n-m)!}mod{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}mod{a}^b,x=v_a(n!),y=v_a(m!),z=v_a((n-m)!)$。
然后我们可以用 $\mathrm{F}(n,a)$ 表示 $\frac{n!}{a^{v_a(n!)}}$，用 $\mathrm{G}(n,a)$ 表示 $v_a(n!)$，那么原式就是 $\frac{\mathrm{F}(n,a)}{\mathrm{F}(m,a)\times \mathrm{F}(n-m,a)}\times a^{\mathrm{G}(n,a)-\mathrm{G}(m,a)-\mathrm{G}(n-m,a)}mod{a}^b$。

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

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

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

由于这个数可能很大，我们求的其实是 $\mathrm{F}(n,a)modA$，其中 $A$ 就是 $a^b$。
好像就是 $n!\times \mathrm{inv}(a^{\mathrm{G}(n,a)})modA$。但是这样我们要预处理出 $1!\to n!$，时间复杂度 $O(n)$，超时。
考虑用另一种方法计算。我们注意到 $n!=1\times 2\times \dots \times (A-1)\times A\times (A+1)\times \dots \times (2A-1)\times 2A\times (2A+1)\times \dots \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 \dots \times n$。
由于 $\mathrm{F}(n,a)$ 中 $a$ 的倍数全部被去掉（注意 $a$ 是个质数），所以我们假设 $\mathrm{H}(s,a)=\frac{s}{a^{v_a(s)}}$。可以发现如果 $gcd(a,t)=1,t\mid s$，那么 $\mathrm{H}(s,a)=t\times \mathrm{H}(\frac{s}{t},a)$。
所以
$\begin{array}{}\text{(1)}& \mathrm{F}(n,a)& \equiv \mathrm{H}(1\times 2\times \dots \times (A-1)\times 1\times (A+1)\times \dots \times (2A-1)\times 2\times (2A+1)\times \dots \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 \dots \times n,a)(\mathrm{mod}A)\text{(2)}& & \equiv \mathrm{H}((1\times 2\times \dots \times (A-1))\times 1\times (1\times 2\times \dots \times (A-1))\times 2\times (1\times \dots \times (A-1))\times \lfloor \frac{n}{A}\rfloor \times (1\times \dots \times (nmodA)),a)(\mathrm{mod}A)\text{(3)}& & \equiv \mathrm{H}(((A-1)!{)}^{\lfloor \frac{n}{A}\rfloor }\times (nmodA)!\times \dots \times \lfloor \frac{n}{A}\rfloor !,a)(\mathrm{mod}A)\text{(4)}& & \equiv \mathrm{H}(((A-1)!{)}^{\lfloor \frac{n}{A}\rfloor }\times (nmodA)!\times \dots \times F(\lfloor \frac{n}{A}\rfloor ,a),a)(\mathrm{mod}A)\text{(5)}& & \equiv ((A-1)!{)}^{\lfloor \frac{n}{A}\rfloor }\times (nmodA)!\times \dots \times F(\lfloor \frac{n}{A}\rfloor ,a)(\mathrm{mod}A)\end{array}$
递归调用 $O(\log n)$ 层，提前预处理出 $[1\to (A-1)]!modA$ 每一层 $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)]!mod{a}_i^{b_i}$，总时间复杂度是 $O(\sum _{i=1}^c{a}_i^{b_i})$。
3. 求 $\mathrm{F},\mathrm{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 _{i=1}^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;
}