733 字
4 分钟
扩展Lucas定理求组合数
2018-04-20

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

我们令 p=p1k1+p2k2+p3k3++pqkqp = p_{1}^{k_{1} } + p_{2}^{k_{2} } + p_{3}^{k_{3} } + \cdots + p_{q}^{k_{q} }
可得同余方程
{ansc1(modp1k1)ansc2(modp2k2)ansc2(modp2k2)anscq(modpqkq)\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.

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

我们考虑如何求 x!%pikix! \% p_{i}^{k_{i} }
x=19,pi=2,ki=2x = 19, p_{i} = 2, k_{i}=2 为例
x!=1×2×3×4×5×6×7×8×9×10×11×12×13×14×15×16×17×18×19x! = 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×2×4×5×7×8×10×11×13×14×16×17×19)×(1×2×3×4×5×6)×36= (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!n!可以分为3部分: 第一部分是pip_i的幂的部分可以直接求解 pinpi{p_i}^{\lfloor{n\over p_i}\rfloor}
第二部分是一个新的阶乘npi!{\lfloor{n\over p_i}\rfloor}!
发现第三部分在模pikip_{i}^{k_{i} }意义下是以pikip_{i}^{k_{i} }为周期的然后就可以较轻松的求出了

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

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

pkpkp^k, exPpip_{i} phipϕ(pk)\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);
}
扩展Lucas定理求组合数
https://www.nekomio.com/posts/148/
作者
NekoMio
发布于
2018-04-20
许可协议
CC BY-NC-SA 4.0