扩展Lucas定理求组合数 | NekoMio's Blog

在求组合数的时候， 我们可能遇到模数是非质数的情况， 这时正常的Lucas可能无法解决问题。 所以我们要用到扩展Lucas定理

我们令 $p = p_{1}^{k_{1} } + p_{2}^{k_{2} } + p_{3}^{k_{3} } + \cdots + p_{q}^{k_{q} }$
可得同余方程
$$ \left \lbrace { \begin{array}{c} ans \equiv c_1\pmod { {p_1}^{k_1} } \\ ans\equiv c_2\pmod { {p_2}^{k_2} } \\ ans\equiv c_2\pmod { {p_2}^{k_2} } \\ \ldots \\ ans\equiv c_q\pmod { {p_q}^{k_q} } \end{array} } \right.$$

然后我们考虑求 $C_{n}^{m} % p_{i}^{k_{i} }$ 的值
因为 $C_{n}^{m} = \frac{n!}{n!(n - m)!}$
我们只要求出$n! % p_{i}^{k_{i} }$, $m! % p_{i}^{k_{i} }$, $(n - m)! % p_{i}^{k_{i} }$

我们考虑如何求 $x! % p_{i}^{k_{i} }$
以 $x = 19, p_{i} = 2, k_{i}=2$ 为例
$x! = 1 \times 2 \times 3 \times 4 \times 5 \times 6 \times 7 \times 8 \times 9 \times 10 \times 11 \times 12 \times 13 \times 14 \times 15 \times 16 \times 17 \times 18 \times 19$
$= (1 \times 2 \times 4 \times 5 \times 7 \times 8 \times 10 \times 11 \times 13 \times 14 \times 16 \times 17 \times 19) \times (1 \times 2 \times 3 \times 4 \times 5 \times 6) \times 3^6$

求解$n!$可以分为3部分: 第一部分是$p_i$的幂的部分可以直接求解 ${p_i}^{\lfloor{n\over p_i}\rfloor}$
第二部分是一个新的阶乘${\lfloor{n\over p_i}\rfloor}!$
发现第三部分在模$p_{i}^{k_{i} }$意义下是以$p_{i}^{k_{i} }$为周期的然后就可以较轻松的求出了

最后一个问题是对于求出的$m! % p_{i}^{k_{i} }$， $(n - m)! % p_{i}^{k_{i} }$ 有可能与 $p_{i}^{k_{i} }$ 不互质。
我们需要将 $p_{i}^x$ 拆出来考虑就可以了

计算$n!$中质因子$p$的个数$x$的公式为$x=\lfloor{n\over p}\rfloor+\lfloor{n\over p^2}\rfloor+\lfloor{n\over p^3}\rfloor+\ldots$
递推式也可以写为$f(n)=f(\lfloor{n\over p}\rfloor)+\lfloor{n\over p}\rfloor$

pk 为 $p^k$, exP为$p_{i}$ phip 为 $\phi(p^k)$

long long Mul(int n, int id)
{
	if (n == 0) return 1;
	long long ans = 1;
	if (n / pk[id])
	{
		for (int i = 2; i <= pk[id]; i++)
			if (i % exP[id])
				ans = ans * i % pk[id];
		// ans = pow_mod(TT[id][pk[id]], n / pk[id], pk[id]);
	}
	// ans = ans * TT[id][n % pk[id]] % pk[id];
	for (int i = 2; i <= n % pk[id]; i++)
		if (i % exP[id])
			ans = ans * i % pk[id];
	return ans * Mul(n / exP[id], id) % pk[id];
}
int exlucas(int n, int m, int id)
{
	if (n < m || n < 0 || m < 0) return 0;
	if (m == n || m == 0) return 1;
	long long a = Mul(n, id), b = Mul(m, id), c = Mul(n - m, id);
	int t = 0;
	for (int i = n; i; i /= exP[id]) t += i / exP[id];
	for (int i = m; i; i /= exP[id]) t -= i / exP[id];
	for (int i = n - m; i; i /= exP[id]) t -= i / exP[id];
	return a * pow_mod(b, phip[id] - 1, pk[id]) % pk[id] * pow_mod(c, phip[id] - 1, pk[id]) % pk[id] * pow_mod(exP[id], t, pk[id]) % pk[id];
}
long long CRT(int *a, int *b, int n)
{
    long long N = 1, Ni, now, ans = 0;
    for (int i = 1; i <= n; i++) N *= a[i];
    for (int i = 1; i <= n; i++)
    {
        Ni = N / a[i];
        now = pow_mod(Ni, phip[i] - 1, a[i]);
        ans = (ans + (b[i] * now % N) * Ni % N) % N;    
    }
    return ans;
}
long long Calc(int n, int m)
{
	if (n < 0 || m < 0 || n < m) return 0;
	for (int i = 1; i <= t; i++)
	{
		b[i] = exlucas(n, m, i);
	}
	return CRT(pk, b, t);
}