【学习笔记】扩展Lucas定理

0.概述

用于解决组合数取模的问题。请先学习中国剩余定理，有点作用的。

1.分解

现在我们要求 $\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}p$ 的结果，$p$ 是个小家伙，$n,m$ 长的很魁梧。

现在我们找到 $p$ 的质因数分解 $$p=\prod _i{p}_i^{t_i}$$

我们如果能求出 $$\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{t_i}$$

再用中国剩余定理合成回去即可了！所以下文就只讨论对 $p^k(p\in \mathbb{P})$ 取模的情况。

2.硬算

用 $$\left(\genfrac{}{}{0ex}{}{n}{m}\right)=\frac{n!}{m!(n-m)!}$$

我们试着单独考虑 $p$ 的次数，就可以求逆元了。用 $f(n)=\frac{n}{p^x}\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$ 来表示，其中 $x$ 是满足 $p^x|n$ 最大的 $x$（即去除 $n$ 的质因数分解的 $p$ 那一项）。

我们将 $n!$ 分成两部分，一部分是 $p$ 的倍数，另一部分不是。即 $$n!=\prod _{i=1}^{\lfloor \frac{n}{p}\rfloor }(ip)\cdot \prod _{\neg (p|i)}^{i\le n}i$$

后面那一部分大有搞头。毕竟你要对 $p^k$ 取模，做乘法之前也可以先行取模。我们根据 “除以 $p^k$ 的商” 进行分类，就可以看出
$$\prod _{\neg (p|i)}^{i\le n}i\equiv {\left(\prod _{\neg (p|i)}^{i\le p^k}i\right)}^{\lfloor \frac{n}{p^k}\rfloor }\times \prod _{\neg (p|i)}^{i\le (n\mathrm{m}\mathrm{o}\mathrm{d}{p}^k)}i(\mathrm{m}\mathrm{o}\mathrm{d}{p}^k)$$

然后前面那一坨把 $p$ 提出来，得到
$$\prod _{i=1}^{\lfloor \frac{n}{p}\rfloor }(ip)=p^{\lfloor \frac{n}{p}\rfloor }\prod _{i=1}^{\lfloor \frac{n}{p}\rfloor }i=p^{\lfloor \frac{n}{p}\rfloor }\times \lfloor \frac{n}{p}\rfloor !$$

注意到 $\lfloor \frac{n}{p}\rfloor !$ 中可能仍然有 $p$ 的因数，于是要递归处理。总结一下就是
$$f(n)\equiv f\left(\lfloor \frac{n}{p}\rfloor \right)\times {\left(\prod _{\neg (p|i)}^{i\le p^k}i\right)}^{\lfloor \frac{n}{p^k}\rfloor }\times \prod _{\neg (p|i)}^{i\le (n\mathrm{m}\mathrm{o}\mathrm{d}{p}^k)}i(\mathrm{m}\mathrm{o}\mathrm{d}{p}^k)$$

复杂度呢？如果直接按照这个式子递归，很明显 $T(n)=T(\lfloor \frac{n}{p}\rfloor )+\mathcal{O}(p^k\log n)=\mathcal{O}(p^k{\log }^{2}n)$ 。如果你再细心一些，你发现 ${\prod }_{\neg (p|i)}^{i\le p^k}i$ 都是以幂的形式乘在一起。包括 $n\mathrm{m}\mathrm{o}\mathrm{d}{p}^k$ 也唯一决定了最后一项的值。所以你考虑 预处理 $x\le p_k$ 的前缀积，而后可以 $\mathcal{O}(1)$ 算后两项的值。所以是 $\mathcal{O}(\log n)$ 单次的。

别指望它是万能的，如果 $p=2^{32}$ ，那么复杂度还是 $\mathcal{O}(n)$ 的了。中国剩余定理合并的时间是不用考虑的，很小。

把 $p$ 的因数提了出来，$p$ 的指数是几呢？众所周知，
$$\underset{x\in {\mathbb{N}}^{+},p^x|n}{\max }x=\sum _{i=1}^{+\infty }\lfloor \frac{n}{p^i}\rfloor$$

然后就完了。可以提一下 $k=1$ 的特殊情况，$$\left(\genfrac{}{}{0ex}{}{n}{m}\right)\equiv \left(\genfrac{}{}{0ex}{}{n/p}{m/p}\right)\left(\genfrac{}{}{0ex}{}{n\mathrm{m}\mathrm{o}\mathrm{d}p}{m\mathrm{m}\mathrm{o}\mathrm{d}p}\right)(\mathrm{m}\mathrm{o}\mathrm{d}p),p\in \mathbb{P}$$

除法是下取整，$\mathbb{P}$ 是质数集合。

代码

这个代码是洛谷板题的，使用暴力递归法。

inline int qkpow(int_ base,int_ q,int Mod){
	int ans = 1; base %= Mod;
	for(; q; q>>=1,base=1ll*base*base%Mod)
		if(q&1) ans = 1ll*ans*base%Mod;
	return ans;
}

int exgcd(int a,int b,int_ &x,int_ &y){
	if(b == 0){
		x = 1, y = 0; return a;
	}
	int d = exgcd(b,a%b,y,x);
	y -= (a/b)*x; return d;
}

int getInv(int x,int p){
	int_ res, useless;
	if(exgcd(x%p,p,res,useless) != 1)
		return -1; // impossible
	return (res%p+p)%p;
}

int CRT(int n,int m[],int c[]){
	int M = 1, ans = 0;
	for(int i=1; i<=n; ++i)
		M *= m[i];
	for(int i=1; i<=n; ++i){
		int tmp = M/m[i];
		tmp *= getInv(tmp,m[i]);
		ans = (ans+1ll*tmp*c[i])%M;
	}
	return ans;
}

int function(int_ n,int p,int pk){
	if(n <= 1) return 1; int res = 1;
	if(n >= pk){
		for(int i=1; i<=pk; ++i)
			if(i%p != 0)
				res = 1ll*res*i%pk;
		res = qkpow(res,n/pk,pk);
	}
	for(int i=1; i<=n%pk; ++i)
		if(i%p != 0)
			res = 1ll*res*i%pk;
	return 1ll*res*function(n/p,p,pk)%pk;
}

int_ calc(int_ n,int p){
	if(n < p) return 0; return n/p+calc(n/p,p);
}
int exLucas(int_ n,int_ m,int p,int pk){
	int_ pt = calc(n,p)-calc(m,p)-calc(n-m,p);
	int_ jb = function(n,p,pk);
	jb = 1ll*jb*getInv(function(m,p,pk),pk)%pk;
	jb = 1ll*jb*getInv(function(n-m,p,pk),pk)%pk;
	return 1ll*jb*qkpow(p,pt,pk)%pk;
}

int mods[50], yu[50];
int Lucas(int_ n,int_ m,int p){
	int cnt = 0;
	for(int i=2,j; 1ll*i*i<=p; ++i){
		if(p%i != 0) continue;
		for(j=1; p%i==0; j*=i) p /= i;
		mods[++ cnt] = j;
		yu[cnt] = exLucas(n,m,i,j);
	}
	if(p != 1){
		mods[++ cnt] = p;
		yu[cnt] = exLucas(n,m,p,p);
	}
	return CRT(cnt,mods,yu);
}

预处理版适用于模数固定但需要多次计算组合数的情况。复杂度 $\mathcal{O}(k{p}^k\log p)$ 预处理 $+\mathcal{O}(\log n)$ 每次询问。

template < int p, int pk >
struct LUCAS{
	int jc[pk+1]; // factorial
	int invjc[pk+1]; // inv of jc
	LUCAS(){
		int phi = pk/p*(p-1);
		jc[0] = invjc[0] = 1;
		jc[1] = invjc[1] = 1;
		for(int i=2; i<=pk; ++i){
			jc[i] = jc[i-1]; // copy
			invjc[i] = invjc[i-1];
			if(i%p == 0) continue;
			jc[i] = 1ll*jc[i]*i%pk;
			invjc[i] = 1ll*invjc[i]
				*qkpow(i,phi-1,pk)%pk;
		}
	}
	int Lucas(int_ n,int_ m){
		if(n < m || m < 0) return 0;
		int z = 0; int_ ans = 1;
		for(int_ i=n; i; i/=p){
			ans = ans*jc[i%pk]%pk;
			z += i/pk; // how many [1,pk]
		}
		for(int_ i=m; i; i/=p){
			ans = ans*invjc[i%pk]%pk;
			z -= i/pk; // cashier
		}
		for(int_ i=n-m; i; i/=p){
			ans = ans*invjc[i%pk]%pk;
			z -= i/pk; // counter
		}
		if(z > 0) // something good
			ans = ans*qkpow(jc[pk],z,pk)%pk;
		else if(z < 0) // can't be worse
			ans = ans*qkpow(invjc[pk],-z,pk)%pk;
		z = 0; // 现在求 C(n,m) 中 p 的指数
		for(int_ i=p; i<=n; i*=p)
			z += (n/i)-(m/i)-((n-m)/i);
		ans = ans*qkpow(p,z,pk)%pk;
		return ans; // 外接 CRT
	}
};