在求组合数的时候， 我们可能遇到模数是非质数的情况， 这时正常的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)$

1
long long Mul(int n, int id)
2
{
3
  if (n == 0) return 1;
4
  long long ans = 1;
5
  if (n / pk[id])
6
  {
7
    for (int i = 2; i <= pk[id]; i++)
8
      if (i % exP[id])
9
        ans = ans * i % pk[id];
10
    // ans = pow_mod(TT[id][pk[id]], n / pk[id], pk[id]);
11
  }
12
  // ans = ans * TT[id][n % pk[id]] % pk[id];
13
  for (int i = 2; i <= n % pk[id]; i++)
14
    if (i % exP[id])
15
      ans = ans * i % pk[id];
16
  return ans * Mul(n / exP[id], id) % pk[id];
17
}
18
int exlucas(int n, int m, int id)
19
{
20
  if (n < m || n < 0 || m < 0) return 0;
21
  if (m == n || m == 0) return 1;
22
  long long a = Mul(n, id), b = Mul(m, id), c = Mul(n - m, id);
23
  int t = 0;
24
  for (int i = n; i; i /= exP[id]) t += i / exP[id];
25
  for (int i = m; i; i /= exP[id]) t -= i / exP[id];
26
  for (int i = n - m; i; i /= exP[id]) t -= i / exP[id];
27
  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];
28
}
29
long long CRT(int *a, int *b, int n)
30
{
31
    long long N = 1, Ni, now, ans = 0;
32
    for (int i = 1; i <= n; i++) N *= a[i];
33
    for (int i = 1; i <= n; i++)
34
    {
35
        Ni = N / a[i];
36
        now = pow_mod(Ni, phip[i] - 1, a[i]);
37
        ans = (ans + (b[i] * now % N) * Ni % N) % N;
38
    }
39
    return ans;
40
}
41
long long Calc(int n, int m)
42
{
43
  if (n < 0 || m < 0 || n < m) return 0;
44
  for (int i = 1; i <= t; i++)
45
  {
46
    b[i] = exlucas(n, m, i);
47
  }
48
  return CRT(pk, b, t);
49
}