【笔记】 exlucas

扩展 lucas

求 \(\binom{n}{m} \bmod p\)，其中 \(1 \le m \le n \le 10^{18},2 \le p \le 10^6\)。

首先我们可以对 \(p\) 质因数分解 \(p=\prod p_i^{a_i}\)，对于每一个质因数我们求出答案，然后用 CRT 合并。

单独考虑每一个模数 \(p^k\)：

我们要求解 \(\binom{n}{m} \bmod p^k\)，其中 \(1 \le m \le n \le 10^{18},2 \le p \le 10^6\) 且 \(p\) 为质数。

障碍

找找是什么阻挡了我们求解的脚步？原因是 \((n-m)!m!\) 有可能找不到逆元。怎么解决呢？

天真的我曾经以为，当 \(n\) 足够大的时候，\(\binom{n}{m} \bmod p^k=0\)，原因是我觉得 \(n\) 会有很多 \(p\) 的因子，于是乎 \(\bmod p\) 完之后就成了 \(0\)。

且不论这个想法是否正确，至少它指明了一个方向：求出 \(\binom{n}{m}\) 有多少个 \(p\) 的因子。

\[\binom{n}{m} = \dfrac{n!}{(n-m)!m!}=\dfrac{\frac{n!}{p^x}}{\frac{(n-m)!}{p^y}\times \frac{m!}{p^z}} p^{x-y-z} \]

这样子，前面这部分就不含 \(p\) 的因子，就可以求出逆元了。

于是我们现在的问题分如下两部分：

1. 求解 \(n!\) 有多少个 \(p\) 的因子；
2. 求解 \(n!\) 去除 \(p\) 的因子后的值 \(\bmod p^k\)。

Task1

不难把 \(n!\) 拆解成如下部分：

\[n!\equiv p^{\lfloor \frac{n}{p} \rfloor}\prod_{i=1}^{\lfloor \frac{n}{p}\rfloor}i\times \left(\prod_{i=1\land i \not\equiv0\pmod p}^{p^k}i\right)^{\lfloor \frac{n}{p^k} \rfloor} \left( \prod_{i=p^k+1\land i\not\equiv0 \pmod p}^{n}i \right) \pmod{p^k} \]

后面两个带括号的部分都不带 \(p\) 因子，前面 \(p^{\lfloor \frac{n}{p} \rfloor}\) 贡献了 \(\lfloor \frac{n}{p} \rfloor\)，\(\prod_{i=1}^{\lfloor \frac{n}{p}\rfloor}i\) 是一个阶乘，可以递归下去。

不难发现可以 \(O(\log_p n)\) 计算 \(p\) 的指数。

Task2

按照 Task1 的式子计算。

不难发现，预处理一下 \(\prod_{i=1\land i \not\equiv0\pmod p}^{p^k}i\) 这个前缀积之后就又可以递归计算了，复杂度 \(O(\log_p n \log n)\)，多出来一个 \(\log n\) 是因为快速幂，如果需要多次计算组合数的话，可以预处理一下次幂，实现 \(O(1)\) 快速幂。

---

计算 \(\dfrac{\frac{n!}{p^x}}{\frac{(n-m)!}{p^y}\times \frac{m!}{p^z}} p^{x-y-z}\) 时要用 exgcd 求逆元，或者算 \(\varphi\) 求快速幂。

最后用 CRT 把各个质数次方的答案合并即可。

Code

#include <bits/stdc++.h>

#define debug(...) fprintf(stderr ,__VA_ARGS__)
#define __FILE(x)\
	freopen(#x".in" ,"r" ,stdin);\
	freopen(#x".out" ,"w" ,stdout)
#define LL long long

const int MX = 1e6 + 23;

LL read(){
	char k = getchar(); LL x = 0;
	while(k < '0' || k > '9') k = getchar();
	while(k >= '0' && k <= '9') x = x * 10 + k - '0' ,k = getchar();
	return x;
}

LL qpow(LL a ,LL b ,LL p){
	LL ans = 1;
	while(b){if(b & 1) ans = ans * a % p;
		a = a * a % p ,b >>= 1;
	}return ans;
}

void exgcd(LL a ,LL b ,LL &x ,LL &y){
	if(!b) return x = 1 ,y = 0 ,void();
	exgcd(b ,a % b ,y ,x); y -= a / b * x;
}

LL inv(LL x ,LL p){
	LL ret ,tmp;
	exgcd(x ,p ,ret ,tmp);
	return (ret + p) % p;
}

LL rfac[MX];
LL other(LL n ,LL p ,LL pk){
	if(n == 0) return 1LL;
	LL ans = other(n / p ,p ,pk);
	ans = ans * qpow(rfac[pk] ,n / pk ,pk) % pk;
	ans = ans * rfac[n % pk] % pk;
	return ans;
}

LL calcexp(LL n ,LL p){
	if(n < p) return 0LL;
	return n / p + calcexp(n / p ,p);
}

void prework(LL p ,LL pk){
	rfac[0] = 1;
	for(int i = 1 ; i <= pk ; ++i){
		if(i % p == 0) rfac[i] = rfac[i - 1];
		else rfac[i] = rfac[i - 1] * i % pk;
	}
}

LL a[MX] ,xp[MX] ,cnt;
LL b[MX];

LL exlucas(LL n ,LL m ,LL p){
	int f = p;
	for(int i = 2 ; i * i <= f ; ++i){
		if(f % i) continue;
		a[++cnt] = i;
		xp[cnt] = 1;
		while(f % i == 0){
			f /= i;
			xp[cnt] *= i;
		}
	}
	if(f != 1) a[++cnt] = f ,xp[cnt] = f;
	for(int u = 1 ; u <= cnt ; ++u){
		LL P = a[u] ,PK = xp[u];
		prework(P ,PK);
		LL ex = calcexp(n ,P) - calcexp(m ,P) - calcexp(n - m ,P);
	   	LL A = other(n ,P ,PK) * inv(other(m ,P ,PK) ,PK) % PK * inv(other(n - m ,P ,PK) ,PK) % PK;
		b[u] = qpow(P ,ex ,PK) * A % PK;
	}

	LL ans = 0;
	for(int u = 1 ; u <= cnt ; ++u){
		ans = (ans + (p / xp[u]) * inv(p / xp[u] ,xp[u]) % p * b[u]) % p;
	}
	return ans;
}

int main(){
	LL n = read() ,m = read() ,p = read();
	printf("%lld\n" ,exlucas(n ,m ,p));
	return 0;
}