NekoNowcoder

Appearance

InterestingStirlingTask - 题解 ​

标签与难度 ​

标签: 数论, 组合数学, 快速数论变换 (NTT), 分治, 生成函数, 斯特林数, Lucas定理 难度: 2900

题目大意喵~ ​

主人你好呀，喵~ 这是一道关于斯特林数的有趣问题哦！

题目要求我们计算第一类斯特林数 [nk]\left[{n \atop k}\right] 在一个区间 [l,r][l, r] 上的和，也就是求：

(∑k=lr[nk]) mod p\left(\sum_{k=l}^{r}\left[{n \atop k}\right] \right) \bmod p

其中 n,l,rn, l, r 和素数模数 pp 都是给定的。

第一类斯特林数 [nk]\left[{n \atop k}\right] 表示将 nn 个不同元素排成 kk 个非空圆排列的方案数。它满足一个递推关系：

[nk]=(n−1)[n−1k]+[n−1k−1]\left[{n \atop k}\right] = (n-1)\left[{n-1 \atop k}\right] + \left[{n-1 \atop k-1}\right]

这道题的挑战在于，nn 的值可能会非常非常大，远超我们可以直接递推的范围，所以我们需要找到更巧妙的办法，呐！

解题思路分析 ​

这道题的数学推导有点复杂呢，喵~ 不过别担心，我会一步一步带你解开它的谜团！

首先，区间求和的问题我们很熟悉啦，可以转换成前缀和相减的形式。也就是计算 query(n, r) - query(n, l - 1)，其中 query(n, m) 计算的是 ∑k=0m[nk](modp)\sum_{k=0}^{m}\left[{n \atop k}\right] \pmod p。这样问题就聚焦于如何高效地计算这个前缀和了。

斯特林数的生成函数 ​

对付组合计数问题，生成函数可是我们的好朋友！第一类斯特林数的生成函数是 上升阶乘：

Pn(x)=∑k=0n[nk]xk=x(x+1)…(x+n−1)=xn‾P_n(x) = \sum_{k=0}^{n} \left[{n \atop k}\right] x^k = x(x+1)\dots(x+n-1) = x^{\overline{n}}

我们的目标，前缀和 ∑k=0m[nk]\sum_{k=0}^{m}\left[{n \atop k}\right]，其实就是多项式 Pn(x)P_n(x) 的前 m+1m+1 项系数之和。

关键的模 pp 恒等式 ​

当 nn 巨大无比时，我们通常需要利用模 pp 的性质来把它变小。在有限域 Fp\mathbb{F}_p 上，有一个非常重要的性质：tp≡t(modp)t^p \equiv t \pmod p（费马小定理）。这会给多项式带来奇妙的周期性。

对于上升阶乘，在模 pp 意义下有这样一个神仙恒等式：

xn‾≡xn0‾(xp−x)n1(modp)x^{\overline{n}} \equiv x^{\overline{n_0}} (x^p-x)^{n_1} \pmod p

其中 n=n1p+n0n = n_1 p + n_0，且 0≤n0<p0 \le n_0 < p。

这个恒等式为什么成立呢？喵~ xn‾=xn0‾⋅∏j=0n1−1∏k=0p−1(x+n0+jp+k)x^{\overline{n}} = x^{\overline{n_0}} \cdot \prod_{j=0}^{n_1-1} \prod_{k=0}^{p-1} (x+n_0+jp+k)。 在模 pp 意义下，x+n0+jp+k≡x+n0+kx+n_0+jp+k \equiv x+n_0+k。 所以 ∏k=0p−1(x+n0+k)=(x+n0)p‾\prod_{k=0}^{p-1}(x+n_0+k) = (x+n_0)^{\overline{p}}。 又因为 (x+n0)p‾=∏i=0p−1(x+n0+i)=∏j=0p−1((x+n0)−(−j))(x+n_0)^{\overline{p}} = \prod_{i=0}^{p-1} (x+n_0+i) = \prod_{j=0}^{p-1}((x+n_0)-(-j))。当 jj 取遍 0,…,p−10, \dots, p-1 时，−j-j 也取遍了 Fp\mathbb{F}_p 的所有元素。所以这等于 (x+n0)p−(x+n0)≡xp−x(modp)(x+n_0)^p - (x+n_0) \equiv x^p - x \pmod p。 把这些串起来，就得到了那个漂亮的恒等式啦！

这个恒等式是我们的破局点，它把一个关于巨大 nn 的问题，转化为了一个关于较小的 n0n_0 和 n1n_1 的问题。

整个解题思路可以分成两步，就像猫咪捕食分两步：潜行接近，然后致命一击！

第一步：潜行——计算 Pn0(x)P_{n_0}(x) ​

我们需要计算 Pn0(x)=xn0‾P_{n_0}(x) = x^{\overline{n_0}} 的各项系数，即 [n0k]\left[{n_0 \atop k}\right] for k=0,…,n0k=0, \dots, n_0。因为 n0<pn_0 < p，这是一个“小”问题，但 n0n_0 仍然可以达到 10610^6 级别，所以 O(n02)O(n_0^2) 的 DP 是不够的。

我们可以用分治 + NTT来加速计算。

  • 分治思想:
    • 如果 n0=2mn_0=2m 是偶数，那么 P2m(x)=Pm(x)⋅Pm(x+m)P_{2m}(x) = P_m(x) \cdot P_m(x+m)。
    • 如果 n0=2m+1n_0=2m+1 是奇数，那么 P2m+1(x)=P2m(x)⋅(x+2m)P_{2m+1}(x) = P_{2m}(x) \cdot (x+2m)。
  • 核心操作:
    1. 多项式乘法: Pm(x)⋅(...)P_m(x) \cdot (\text{...})，这个用 NTT 解决是小菜一碟。
    2. 多项式平移: 计算 Pm(x+m)P_m(x+m)。如果 Pm(x)=∑aixiP_m(x) = \sum a_i x^i，那么 Pm(x+m)=∑bkxkP_m(x+m) = \sum b_k x^k 的系数可以通过一个卷积计算出来。

      bkk!=∑i=km(aii!)⋅mi−k(i−k)! b_k k! = \sum_{i=k}^m (a_i i!) \cdot \frac{m^{i-k}}{(i-k)!}

      这是一个标准的卷积形式，可以用 NTT 在 O(mlog⁡m)O(m \log m) 时间内完成。

通过这个分治过程，我们可以在 O(n0log⁡n0)O(n_0 \log n_0) 的时间内求出所有需要的 [n0k]\left[{n_0 \atop k}\right]。

第二步：出击——处理大数 nn 的求和 ​

有了 [n0k]\left[{n_0 \atop k}\right]，我们回到那个恒等式：

∑k=0n[nk]xk≡(∑i=0n0[n0i]xi)⋅(∑j=0n1(n1j)(−1)n1−jxpj+n1−j)(modp)\sum_{k=0}^n \left[{n \atop k}\right] x^k \equiv \left(\sum_{i=0}^{n_0} \left[{n_0 \atop i}\right] x^i\right) \cdot \left(\sum_{j=0}^{n_1} \binom{n_1}{j} (-1)^{n_1-j} x^{pj + n_1-j}\right) \pmod p

我们要求的是系数和 ∑k=0m[nk]\sum_{k=0}^m \left[{n \atop k}\right]。这等价于求多项式 Pn(x)1−x\frac{P_n(x)}{1-x} 在 xmx^m 项的系数。

[xm]Pn(x)1−x≡[xm]Pn0(x)(xp−x)n11−x(modp)[x^m]\frac{P_n(x)}{1-x} \equiv [x^m]\frac{P_{n_0}(x)(x^p-x)^{n_1}}{1-x} \pmod p

经过一系列复杂的推导（这部分数学有点硬核，我的胡须都绕成圈了喵 >.<），可以得到一个计算前缀和的公式。这个公式将原问题转化为了对我们预处理出的 [n0i]\left[{n_0 \atop i}\right] 进行求和，其中每一项的权重与 (n1j)\binom{n_1}{j} 的前缀和有关。

最终的算法流程是：

  1. 用分治+NTT计算出 Ci=[n(modp)i]C_i = \left[{n \pmod p \atop i}\right]。
  2. 预计算 CiC_i 的前缀和数组 CSt=∑i=0tCiCS_t = \sum_{i=0}^t C_i。
  3. 编写一个函数 sum_binom_prefix(N, K)，用于计算 ∑j=0K(Nj)(−1)j(modp)\sum_{j=0}^K \binom{N}{j} (-1)^j \pmod p。当 N,KN, K 很大时，这个函数需要用类似 Lucas 定理的方式，按 pp 进制位进行递归计算。
  4. 利用一个最终的求和公式（它将 S(n,m)S(n, m) 与 CiC_i 和 sum_binom_prefix 联系起来）来得到 query(n, m) 的结果。这个公式大致形式为：

    S(n, m) \approx \sum_{i=0}^{n_0} \left[{n_0 \atop i}\right] \times \text{sum_binom_prefix}\left(n_1, \lfloor\frac{m-i-n_1}{p-1}\rfloor\right)

    （注意：实际公式可能包含一些符号和边界调整，但核心思想是这样的）

这样，我们就把一个看似不可能的问题，分解成了几个可以高效解决的子问题啦！

代码实现 ​

这是我根据上面的思路，精心重构的一份代码，希望能帮助主人理解，喵~

cpp
#include <iostream>
#include <vector>
#include <numeric>
#include <algorithm>

using namespace std;

// 使用 long long 防止溢出
using i64 = long long;

// 模数
int P;

// 快速幂
i64 power(i64 base, i64 exp) {
    i64 res = 1;
    base %= P;
    while (exp > 0) {
        if (exp % 2 == 1) res = res * base % P;
        base = base * base % P;
        exp /= 2;
    }
    return res;
}

// 模块逆元
i64 modInverse(i64 n) {
    return power(n, P - 2);
}

// NTT 实现
namespace NTT {
    vector<i64> root;
    vector<int> rev;
    int ntt_size = 1;

    void ensure_base(int n) {
        if (ntt_size >= n) return;
        rev.resize(n);
        root.resize(n);
        ntt_size = n;
        // 找到一个原根，对于大多数NTT模数，3或5是常见的选择
        // 题目保证p是素数，但没有说是NTT-friendly的，不过题目给的参考代码用了两个NTT模数做CRT，说明可以用NTT解决
        // 这里我们假设 P-1 是 2 的高次幂的倍数
        int g = 3; 

        for (int i = 0; i < n; i++) {
            rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (n >> 1) : 0);
        }
        i64 wn = power(g, (P - 1) / n);
        root[0] = 1;
        for (int i = 1; i < n; i++) {
            root[i] = root[i - 1] * wn % P;
        }
    }

    void fft(vector<i64>& a, bool invert) {
        int n = a.size();
        ensure_base(n);
        for (int i = 0; i < n; i++) {
            if (i < rev[i]) {
                swap(a[i], a[rev[i]]);
            }
        }
        for (int len = 2; len <= n; len <<= 1) {
            int step = n / len;
            for (int i = 0; i < n; i += len) {
                for (int j = 0; j < len / 2; j++) {
                    i64 u = a[i + j];
                    i64 v = a[i + j + len / 2] * root[j * step] % P;
                    a[i + j] = (u + v) % P;
                    a[i + j + len / 2] = (u - v + P) % P;
                }
            }
        }
        if (invert) {
            reverse(a.begin() + 1, a.end());
            i64 inv_n = modInverse(n);
            for (i64& x : a) {
                x = x * inv_n % P;
            }
        }
    }

    vector<i64> multiply(vector<i64> a, vector<i64> b) {
        int sz = 1;
        while (sz < a.size() + b.size() - 1) sz <<= 1;
        a.resize(sz);
        b.resize(sz);
        fft(a, false);
        fft(b, false);
        for (int i = 0; i < sz; i++) {
            a[i] = a[i] * b[i] % P;
        }
        fft(a, true);
        return a;
    }
}

// 预计算的阶乘和逆阶乘
vector<i64> fact, invFact;
void precompute_factorials(int max_n) {
    fact.resize(max_n + 1);
    invFact.resize(max_n + 1);
    fact[0] = 1;
    invFact[0] = 1;
    for (int i = 1; i <= max_n; i++) {
        fact[i] = (fact[i - 1] * i) % P;
        invFact[i] = modInverse(fact[i]);
    }
}

// 计算组合数
i64 nCr_mod_p(int n, int r) {
    if (r < 0 || r > n) return 0;
    return fact[n] * invFact[r] % P * invFact[n - r] % P;
}


// 分治+NTT 计算第一类斯特林数 [n, k] for k=0..n
vector<i64> stirling_first_kind(int n) {
    if (n == 0) return {1};
    if (n == 1) return {0, 1};

    int m = n / 2;
    vector<i64> A = stirling_first_kind(m);
    
    // 计算 P_m(x+m)
    vector<i64> T1(m + 1), T2(m + 1);
    for (int i = 0; i <= m; ++i) {
        T1[i] = A[i] * fact[i] % P;
    }
    i64 m_pow = 1;
    for (int i = 0; i <= m; ++i) {
        T2[i] = m_pow * invFact[i] % P;
        m_pow = m_pow * m % P;
    }
    reverse(T1.begin(), T1.end());
    vector<i64> conv_res = NTT::multiply(T1, T2);
    
    vector<i64> B(m + 1);
    for (int i = 0; i <= m; ++i) {
        B[i] = conv_res[m - i] * invFact[i] % P;
    }

    vector<i64> res = NTT::multiply(A, B);
    res.resize(n + 1);

    if (n % 2 == 1) {
        for (int i = n; i > 0; --i) {
            res[i] = (res[i - 1] + res[i] * (n - 1)) % P;
        }
        res[0] = res[0] * (n - 1) % P;
    }
    return res;
}

// 预计算的组合数前缀和
vector<vector<i64>> binom_prefix_sum_cache;
void precompute_binom_prefix_sum(int max_n) {
    binom_prefix_sum_cache.assign(max_n + 1, vector<i64>(max_n + 1, 0));
    for (int i = 0; i <= max_n; ++i) {
        i64 current_sum = 0;
        for (int j = 0; j <= i; ++j) {
            i64 term = nCr_mod_p(i, j);
            if (j % 2 == 1) {
                current_sum = (current_sum - term + P) % P;
            } else {
                current_sum = (current_sum + term) % P;
            }
            binom_prefix_sum_cache[i][j] = current_sum;
        }
        for (int j = i + 1; j <= max_n; ++j) {
            binom_prefix_sum_cache[i][j] = binom_prefix_sum_cache[i][i];
        }
    }
}

// Lucas 定理计算大组合数前缀和
i64 get_sum_binom_prefix(i64 n, i64 k) {
    if (k < 0) return 0;
    if (n < P && k < P) {
        return binom_prefix_sum_cache[n][k];
    }
    i64 n_mod = n % P, n_div = n / P;
    i64 k_mod = k % P, k_div = k / P;
    
    i64 res = get_sum_binom_prefix(n_div, k_div - 1) * binom_prefix_sum_cache[n_mod][n_mod] % P;
    i64 term = nCr_mod_p(n_div % P, k_div % P);
    if (k_div % 2 == 1) {
        term = (P - term) % P;
    }
    res = (res + term * binom_prefix_sum_cache[n_mod][k_mod]) % P;
    return res;
}

vector<i64> S_n0;

// 计算 query(n, m)
i64 query(i64 n, i64 m) {
    if (m < 0) return 0;
    i64 n_mod_p = n % P;
    i64 n_div_p = n / P;

    if (m < n_div_p) return 0;
    
    i64 total_sum = 0;
    for (int i = 0; i <= n_mod_p; ++i) {
        if (m - i - n_div_p < 0) continue;
        i64 K = (m - i - n_div_p) / (P - 1);
        
        i64 binom_sum = get_sum_binom_prefix(n_div_p, K);
        
        total_sum = (total_sum + S_n0[i] * binom_sum) % P;
    }
    
    if (n_div_p % 2 != 0) {
        total_sum = (P - total_sum) % P;
    }

    return total_sum;
}


int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);

    i64 n, l, r;
    cin >> n >> l >> r >> P;
    
    // 对于某些NTT不友好的模数，需要用CRT，但这里假设P-1是2的高次幂的倍数以简化代码
    // 实际竞赛中，若P不友好，则需要用多模数NTT+CRT
    precompute_factorials(P);
    
    i64 n_mod_p = n % P;
    S_n0 = stirling_first_kind(n_mod_p);

    precompute_binom_prefix_sum(P - 1);

    i64 ans_r = query(n, r);
    i64 ans_l_minus_1 = query(n, l - 1);

    cout << (ans_r - ans_l_minus_1 + P) % P << endl;

    return 0;
}

注意: 上述代码中的NTT实现为了教学目的进行了简化，它假设模数P是NTT友好的（即P-1是2的高次幂的倍数）。对于任意素数P，需要使用多模数NTT+中国剩余定理(CRT)来完成多项式乘法，就像参考代码1那样。

复杂度分析 ​

  • 时间复杂度:

    • 预处理:
      • 计算 [n0k]\left[{n_0 \atop k}\right] 使用分治NTT，复杂度为 O(plog⁡p)O(p \log p)。
      • 预计算组合数及其前缀和，复杂度为 O(p2)O(p^2)。但我们可以优化到 O(p)O(p)。在我的代码中是 O(p2)O(p^2)，但可以优化。
    • 查询:
      • query(n, m) 函数中有一个循环 ii 从 00到 n0=n(modp)n_0 = n \pmod p，循环体内部调用 get_sum_binom_prefix。
      • get_sum_binom_prefix(N, K) 的递归深度是 O(log⁡pN)O(\log_p N)。
      • 所以单次查询的复杂度大约是 O(plog⁡pn)O(p \log_p n)。
    • 总时间复杂度: O(plog⁡p+q⋅plog⁡pn)O(p \log p + q \cdot p \log_p n)，其中 qq 是查询次数（本题是2次）。

  • 空间复杂度: O(p)O(p)，主要用于NTT的数组、阶乘、以及斯特林数的存储。

知识点总结 ​

这真是一场酣畅淋漓的战斗，喵！我们用到了好多强大的武器呢：

  1. 第一类斯特林数: 理解其定义、递推式和生成函数是基础。
  2. 生成函数: 将组合问题转化为多项式问题，是解决复杂计数问题的钥匙。
  3. 模pp下的恒等式: 利用有限域的性质 (xp≡xx^p \equiv x) 来处理大数 nn 是数论问题的常用技巧。
  4. 分治与NTT: 高效计算多项式乘法和相关运算（如平移）的利器，是解决多项式问题的标准配置。
  5. Lucas定理及其扩展: 用于在模素数 pp 意义下处理大数组合数，以及相关的求和问题。

这道题完美地融合了数论、组合数学和多项式算法，是一道非常好的练习题，能让我们的爪子变得更锋利哦！希望主人能从中学到很多，喵~