使用递归分块法求解 洛谷U687409 阶乘消去p再模p的a次方

题目链接：https://www.luogu.com.cn/problem/U687409

解题思路：

求解 \((n!)_p \bmod p^a\) 的高效核心算法是递归分块法（常用于扩展卢卡斯定理 ExLucas）。

当 \(p^a\) 较小时（如 \(p^a \le 10^7\)），该算法可以在 \(O(p^a + \log_p n)\) 的时间复杂度内完美解决。

---

核心数学原理

利用递归思想，我们可以将 \((n!)_p\) 拆分为不含 \(p\) 的因数乘积和含有 \(p\) 的因数乘积两部分：

\[(n!)_p = \left( \prod_{1 \le k \le n, p \nmid k} k \right) \times \left( \lfloor n/p \rfloor ! \right)_p \]

令 \(f(n) = \prod_{1 \le k \le n, p \nmid k} k\)，则原式可展开为完全递归式：

\[(n!)_p = f(n) \times f(\lfloor n/p \rfloor) \times f(\lfloor n/p^2 \rfloor) \dots \]

对于 \(f(n) \bmod p^a\)，我们可以将 \(n\) 按照 \(p^a\) 大小进行分块：

• 设 \(n = q \cdot p^a + r\)，其中 \(r = n \bmod p^a\)。
• 每一个长度为 \(p^a\) 的完整块中，除去 \(p\) 的倍数后，其余数字在模 \(p^a\) 意义下构成的乘积是周期性循环的。
• 根据威尔逊定理的推广（高斯泛化），一个完整块的乘积 \(f(p^a) \equiv -1 \pmod{p^a}\)（唯独当 \(p=2, a \ge 3\) 时结果为 \(1\)）。

因此，分块加速公式为：

\[f(n) \equiv (-1)^q \times f(r) \pmod{p^a} \]

---

算法执行步骤

逐层迭代计算，直至 \(n=0\)：

1. 预处理：线性递推预处理出 \(0\) 到 \(p^a-1\) 范围内所有 \(f(r) \bmod p^a\) 的值，复杂度 \(O(p^a)\)。
2. 循环/递归：

• 取当前 \(r = n \bmod p^a\) 和 \(q = \lfloor n/p^a \rfloor\)。
  • 累乘当前层的贡献：\(f(n) \equiv (-1)^q \times f(r) \pmod{p^a}\)。
  • 令 \(n = \lfloor n/p \rfloor\)，进入下一轮循环，直到 \(n = 0\)。

3. 输出：所有层的贡献相乘即为最终答案。

示例程序

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int maxn = 1e7 + 5;

ll n, p, a, m = 1, f[maxn];

void init() {
    f[0] = 1;
    for (int i = 1; i < m; i++) {
        f[i] = f[i-1];
        if (i % p)
            (f[i] *= i) %= m;
    }
}

ll get_f(ll n) {
    if (n < m)
        return f[n];
    ll res = f[n % m], q = n / m;
    if (q % 2 && !(p == 2 && a >= 3))
        res = -res;
    return res;
}

int main() {
    cin >> n >> p >> a;
    for (int i = 0; i < a; i++) m *= p;
    init();
    ll ans = 1;
    for (ll i = n; i; i /= p)
        (ans *= get_f(i)) %= m;
    ans = (ans + m) % m;
    cout << ans << endl;
    return 0;
}