牛客多校第三场 I 题题解

Ice Drinking

https://ac.nowcoder.com/acm/contest/33188/I

题意：给定 $n$ ，记随机变量 $X$ 为一个长度为 $n$ 的排列中满足 $a_{i} = i$ 的个数，求所有排列中 $E (X^{k})$ 。 $1 \leq n\leq 1\times 10^{18}$ ， $0 \leq k \leq n+5000$ 。

解法：由错排定义可得 $\displaystyle n!=\sum_{i=0}^n {n \choose i}D_{n-i}$ 。

\begin{aligned} E(X^k)=&\sum_{i=0}^n i^k \Pr[x=i]\\ =&\sum_{i=0}^n i^k {n \choose i} D_{n-i}\\ =&\sum_{i=0}^n D_{n-i}{n \choose i}\sum_{j=0}^k j!{i \choose j}{k \brace j}\\ =&\sum_{i=0}^n \sum_{j=0}^k {n \choose i}{i \choose j}j!{k \brace j}D_{n-i}\\ =&\sum_{i=0}^n \sum_{j=0}^k {n \choose j}{n-j \choose i-j}j!{k \brace j}D_{n-i}\\ =&\sum_{j=0}^k {n \choose j} {k \brace j}j!\left(\sum_{i=j}^n{n-j \choose i-j} D_{n-i}\right)\\ =&\sum_{j=0}^k {n \choose j} {k \brace j}j!\left(\sum_{i=0}^{n-j}{n-j \choose i} D_{n-i-j}\right)\\ =&\sum_{j=0}^k{n \choose j}{k \brace j}j!(n-j)!\\ =&\sum_{j=0}^k\dfrac{n!}{j!(n-j)!}{k \brace j}j!(n-j)!\\ =&n!\sum_{j=0}^n {k \brace j} \end{aligned}

因而对于 $n \geq k$ ，答案为贝尔数 $\displaystyle \varpi_n=\sum_{i=0}^k {k \brace i}$ ；对于 $n < k$ ，需要手动计算出最后几项斯特灵数，然后进行相减。

求贝尔数对 $p$ 取模，可以参考本文https://www.cnblogs.com/lfri/p/11546313.html：利用 $\varpi_{n+p}≡\varpi_n + \varpi_{(n+1)} \pmod p$ ， $p$ 是质数，其时间复杂度为 $\mathcal O(p^2 \log n)$ 。

对于斯特灵数的快速计算，对于本题 $k-n \leq 5\times 10^3$ ，可以采用二阶欧拉数 $\left\langle \left\langle \begin{matrix}n\\k\end{matrix} \right\rangle \right\rangle$ ：

{x \brace x-n}=\sum_{k=0}^n \left\langle \left\langle \begin{matrix}n\\k\end{matrix} \right\rangle \right\rangle {x+n-1-k \choose 2n}

其中二阶欧拉数递推式为 $\left\langle \left\langle \begin{matrix}n\\k\end{matrix} \right\rangle \right\rangle=(k+1)\left\langle \left\langle \begin{matrix}n-1\\k\end{matrix} \right\rangle \right\rangle+(2n-1-k)\left\langle \left\langle \begin{matrix}n-1\\k-1\end{matrix} \right\rangle \right\rangle$ 。该时间复杂度为 $\mathcal O((k-n)^2)$ 。

此外还有一个做法，来源于 hos.lyric：

{n \brace k} \bmod p= \begin{cases} \displaystyle {x \choose y} \bmod p,{\rm if}\ p|k,y \leftarrow \dfrac{k}{p},x \leftarrow \left \lfloor \dfrac{n-y}{p-1}\right \rfloor\\ \displaystyle {x \choose y}{n-1-(p-1)x-y \brace k} \bmod p, {\rm if}\ p \nmid k,y \leftarrow \left \lfloor \dfrac{k}{p}\right \rfloor,x \leftarrow \left \lfloor \dfrac{n-y-1}{p-1}\right \rfloor \end{cases}

具体证明可以参看http://matwbn.icm.edu.pl/ksiazki/aa/aa94/aa9413.pdf。

hos 对此的解释为：

I just guessed the pattern by looking the table for $p = 7$ . We can prove that Stirling1 mod p has a nice pattern because $x(x-1)...(x-(p-1)) \equiv x^p-x \pmod p$ , and s***irling2 is its inverse matrix, my intuition says it is correct.

如果提前预处理好 $p^{2}$ 的斯特灵数表和组合数表，利用 Lucas 定理可以 $\mathcal O(\log n)$ 的求出。

由于本题的模数 $862118861=857\times 997\times 1009$ ，因而可以通过中国剩余定理，对于每个质因子单独考虑，然后合并。最终复杂度为 $\mathcal O(p^2\log n+(k-n)^2)$ 。

#include <bits/stdc++.h>
using namespace std;
const int P = 862118861;
const int mod[3] = {857, 997, 1009};
class Bell
{
    const static int N = 1100;
    vector<vector<int>> Stirling;
    vector<int> bell;
    int mod;

public:
    Bell(int mod)
    {
        bell.resize(N + 5);
        Stirling.resize(N + 5);
        for (auto &i : Stirling)
            i.resize(N + 5);
        this->mod = mod;
        Stirling[0][0] = 1;
        for (int i = 0; i <= N; i++)
            for (int j = 1; j <= i; j++)
                Stirling[i][j] = (Stirling[i - 1][j - 1] + 1LL * j * Stirling[i - 1][j]) % mod;
        bell[0] = 1;
        for (int j = 1; j <= N; j++)
        {
            bell[j] = 0;
            for (int k = 1; k <= j; k++)
                bell[j] = (bell[j] + Stirling[j][k]) % mod;
        }
    }
    int query(long long n)
    {
        if(n < mod)
            return bell[n];
        vector<int> B(N + 5, 0), a(N + 5, 0), d(N + 5, 0);
        for (int i = 0; i <= mod;i++)
            B[i] = bell[i];
        int cnt = 0;
        while (n)
        {
            a[cnt++] = n % mod;
            n /= mod;
        }
        for (int i = 1; i < cnt; i++)
            for (int j = 1; j <= a[i]; j++)
            {
                for (int k = 0; k < mod; k++)
                    d[k] = (B[k + 1] + i * B[k] % mod) % mod;
                d[mod] = (d[0] + d[1]) % mod;
                for (int k = 0; k <= mod; k++)
                    B[k] = d[k];
            }
        return d[a[0]];
    }
};
class Binom
{
    int mod;
    const static int N = 1100;
    vector<vector<int>> C;

public:
    void set(int mod)
    {
        C.resize(N + 5);
        for (auto &i : C)
            i.resize(N + 5);
        this->mod = mod;
        for (int i = 0; i <= N;i++)
            C[i][i] = C[i][0] = 1;
        for (int i = 1; i <= N;i++)
            for (int j = 1; j < i;j++)
                C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % mod;
    }
    long long query(long long n, long long m)
    {
        if(n < mod && m < mod)
            return C[n][m];
        else
            return C[n % mod][m % mod] * query(n / mod, m / mod) % mod;
    }
};
class Euler
{
    const static int N = 5000;
    vector<vector<int>> euler;
    int mod;
    Binom b;

public:
    Euler(int mod)
    {
        euler.resize(N + 5);
        for (auto &i : euler)
            i.resize(N + 5);
        b.set(mod);
        this->mod = mod;
        euler[0][0] = 1;
        for (int i = 1; i <= N;i++)
        {
            euler[i][0] = 1;
            euler[i][i] = 0;
        }
        for (int i = 2; i <= N;i++)
            for (int j = 1; j < i;j++)
                euler[i][j] = (euler[i - 1][j] * (j + 1) + (2 * i - 1 - j) * euler[i - 1][j - 1]) % mod;
    }
    long long query_stirling(long long x, int n)
    {
        long long ans = 0;
        for (int k = 0; k <= n;k++)
            ans = (ans + euler[n][k] * b.query(x + n - 1 - k, 2 * n) % mod) % mod;
        return ans;
    }
};
int main()
{
    long long n, m;
    scanf("%lld%lld", &n, &m);
    long long ans[3] = {0};
    for (int i = 0; i < 3;i++)
    {
        Bell solve(mod[i]);
        ans[i] = solve.query(m);
        // printf("SOLVE:%d %lld\n", mod[i], ans[i]);
    }
    if (n < m)
        for (int i = 0; i < 3;i++)
        {
            Euler solve(mod[i]);
            for (int j = 0; j < m - n;j++)
                ans[i] = (ans[i] - solve.query_stirling(m, j) + mod[i]) % mod[i];
        }
    long long res = (463753553LL * ans[0] + 376150155LL * ans[1] + 22215154LL * ans[2]) % P;
    printf("%lld", res);
    return 0;
}