标签：p1 sum Lint 多校 Ice fac binom 第三场 mod

2022牛客暑期多校第三场 I. Ice Drinking

题意

按随机顺序摆放 \(1,2,3,...,n\)，设随机变量 \(x\) 为数字与位置（第几个）相等的个数，给定非负整数 \(k\)，求 \(x^k\) 的期望。

分析

设错排方案为 \(P(n)\)，根据组合意义，枚举正确的个数，剩下的全部错排的方案数之和就是全排列

\[n! = \sum_{i=0}^n\binom{n}{i}P(i) \]

二项式反演得到错排的容斥公式

\[P(n)=\sum_{i=0}^n(-1)^{n-i}\binom{n}{i}i! \]

那么恰好有 \(j\) 个位置正确的情况数为

\[\binom{n}{j}P(n-j) =\sum_{i=0}^{n-j}(-1)^{n-i-j}\binom{n}{j}\binom{n-j}{i}i! \]

于是期望为

\[\frac{1}{n!}\sum_{j=0}^n\sum_{i=0}^{n-j}(-1)^{n-i-j}\binom{n}{j}\binom{n-j}{i}i!j^k \]

稍微化简一下

\[\sum_{j=0}^n\sum_{i=0}^{n-j}(-1)^{n-i-j}\frac{1}{j!(n-j-i)!}j^k \]

众所周知，\(j^k\) 可以利用第二类斯特林数的组合意义，用第二类斯特林数的和式代替（组合含义为 \(j\) 个不同盒子，\(k\) 个不同小球的放法等价于枚举非空盒个数 \(t\)，在 \(j\) 个不同盒子里选择 \(t\) 个盒子，\(k\) 个不同小球划分到 \(t\) 个相同盒子里（盒子不允许为空）的方案乘上非空盒的全排列）

\[j^k=\sum_{t=0}^j \binom{j}{t}t!S(k,t) \]

我们不直接代入，继续反演

\[j!S(k,j)=\sum_{t=0}^j (-1)^{j-t}\binom{j}{t}t^k \]

把 \(j!\) 移到右边去，得到 \(S(k,j)\) 的表达式

\[S(k,j)=\sum_{t=0}^j (-1)^{j-t}\frac{1}{t!(j-t)!}t^k \]

再次观察期望的式子，发现形式很像，还差一点。

\[\sum_{j=0}^n\sum_{i=0}^{n-j}(-1)^{n-i-j}\frac{1}{j!(n-j-i)!}j^k \]

把 \(i,j\) 计算次序交换一下

\[\sum_{i=0}^n\sum_{j=0}^{n-i}(-1)^{n-i-j}\frac{1}{j!(n-j-i)!}j^k \]

发现了不得了的事情，里层求和符号的形式与 \(S(k,n-i)\) 一模一样。

于是期望就变成了

\[\sum_{i=0}^n S(k,n-i)=\sum_{i=0}^n S(k,i) \]

仅仅推出这个式子还不够，\(n\) 范围非常大，难以直接计算。

由 \(n,k\) 范围可以知道 \(k\) 最多比 \(n\) 大 \(5000\)，也就是如果知道一行的斯特林数之和，最多只要减掉 \(5000\) 个数。

一行斯特林数的和就是贝尔数，所以 \(n > k\) 就是第 \(k\) 个贝尔数，而贝尔数有一个很好的同余性质，可以快速计算模素数意义下的贝尔数，复杂度是 \(O(p^2 \log_p n)\) 的，而模数 \(862118861 = 857 * 997 * 1009\)，我们可以求出三个质因子为模数的答案，然后中国剩余定理解决。具体做法，参照这个链接 Bell数。

问题就变成如何逆推斯特林数的一行的后面几个，只需要 \(O((k-n)^2)\) 的递推就够了。这在《具体数学》这本书里有详细的讲解。

由于我写这道题时还不知道《具体数学》书里有相应记载，所以自行采用了另外一种方法进行递推（与书上方法不同，推出的结论看起来也不太相同，但是殊途同归，想看书上方法的请自行查阅这本书的“特殊的数”这一章节），以下是方法。

众所周知，斯特林数有一个递推式（组合意义也十分明确）

\[S(n,m)=S(n-1,m-1)+mS(n-1,m)\\ S(0,0)=1\\ S(n,0)=0(n>0) \]

我们尝试找一找规律，
有

\[S(n,n)=S(n-1,n-1)=...=S(0,0)=1 \]

还有

\[\begin{aligned} &S(n,n-1)\\ =&S(n-1,n-2)+(n-1)\\ =&S(n-2,n-3)+(n-2)+(n-1)\\ ...\\ =&S(1,0)+1+2+...+(n-1)\\ =&\binom{n}{2} \end{aligned} \]

还有

\[\begin{aligned} &S(n,n-2)\\ =&S(n-1,n-3)+(n-2)\binom{n-1}{2}\\ =&S(n-2,n-4)+(n-3)\binom{n-2}{2}+(n-2)\binom{n-1}{2}\\ ...\\ =&S(2,0)+1\cdot\binom{2}{2}+2\cdot\binom{3}{2}+...+(n-2)\binom{n-1}{2}\\ =&\sum_{i=1}^{n-2}i\binom{i+1}{2}\\ =&\sum_{i=1}^{n-2}\binom{i+1}{2}+\sum_{i=1}^{n-2}(i-1)\binom{i+1}{2}\\ =&\sum_{i=1}^{n-2}\binom{i+1}{2}+\frac{3!}{2!}\sum_{i=1}^{n-2}\binom{i+1}{3}\\ =&\binom{n}{3}+3\binom{n}{4} \end{aligned} \]

还有

\[\begin{aligned} &S(n,n-3)\\ =&S(n-1,n-4)+(n-3)\left(\binom{n-1}{3}+3\binom{n-1}{4}\right)\\ =&S(n-2,n-5)+(n-4)\left(\binom{n-2}{3}+3\binom{n-2}{4}\right)\\ ...\\ =&\sum_{i=1}^{n-3}i\left(\binom{i+2}{3}+3\binom{i+2}{4}\right)\\ =&\sum_{i=1}^{n-3}i\binom{i+2}{3}+3\sum_{i=1}^{n-3}i\binom{i+2}{4}\\ =&\sum_{i=1}^{n-3}(i-1)\binom{i+2}{3}+\sum_{i=1}^{n-3}\binom{i+2}{3}+3\sum_{i=1}^{n-3}(i-2)\binom{i+2}{4}+6\sum_{i=1}^{n-3}\binom{i+2}{4}\\ =&\frac{4!}{3!}\sum_{i=1}^{n-3}\binom{i+2}{4}+\sum_{i=1}^{n-3}\binom{i+2}{3}+3\cdot\frac{5!}{4!}\sum_{i=1}^{n-3}\binom{i+2}{5}+6\sum_{i=1}^{n-3}\binom{i+2}{4}\\ =&4\binom{n}{5}+\binom{n}{4}+15\binom{n}{6}+6\binom{n}{5}\\ =&\binom{n}{4}+10\binom{n}{5}+15\binom{n}{6} \end{aligned} \]

推了 \(3\) 个，发现这些都是若干组合数的和（先不管系数具体多少），记系数为 \(p_{i,j}\) （角标是为了后面推这个系数的递推式用的），具体地归纳一下，有

\[S(n,n-k)=\sum_{i=1}^kp_{k,i}\binom{n}{k+i} \]

怎么求系数 \(p_{i,j}\) 呢？利用数学归纳法，就能求递推式了。

假定如下式子对于任意 \(n\) 成立（\(k\) 固定）

\[S(n,n-k+1)=\sum_{i=1}^{k-1}p_{k-1,i}\binom{n}{k-1+i} \]

我们开始由上面的式子推 \(S(n,n-k)\)，方法和上面找规律的推法很像。

\[\begin{aligned} &S(n,n-k)\\ =&S(n-1,n-k-1)+(n-k)S(n-1,n-k)\\ =&S(n-2,n-k-2)+(n-k-1)S(n-2,n-k-1)+(n-k)S(n-1,n-k)\\ ...\\ =&\sum_{i=1}^{n-k}iS(i+k-1,i)\\ =&\sum_{i=1}^{n-k}\sum_{j=1}^{k-1}i\cdot p_{k-1,j}\binom{i+k-1}{k-1+j}\\ =&\sum_{j=1}^{k-1}\sum_{i=1}^{n-k}i\cdot p_{k-1,j}\binom{i+k-1}{k-1+j}\\ =&\sum_{j=1}^{k-1}\sum_{i=1}^{n-k}(i-j)\cdot p_{k-1,j}\binom{i+k-1}{k-1+j}+\sum_{j=1}^{k-1}\sum_{i=1}^{n-k}j\cdot p_{k-1,j}\binom{i+k-1}{k-1+j}\\ =&\sum_{j=1}^{k-1}(k+j)\cdot p_{k-1,j}\sum_{i=1}^{n-k}\binom{i+k-1}{k+j}+\sum_{j=1}^{k-1}j\cdot p_{k-1,j}\sum_{i=1}^{n-k}\binom{i+k-1}{k-1+j}\\ =&\sum_{j=1}^{k-1}(k+j)\cdot p_{k-1,j}\binom{n}{k+j+1}+\sum_{j=1}^{k-1}j\cdot p_{k-1,j}\binom{n}{k+j}\\ =&\sum_{j=2}^{k}(k+j-1)\cdot p_{k-1,j-1}\binom{n}{k+j}+\sum_{j=1}^{k-1}j\cdot p_{k-1,j}\binom{n}{k+j}\\ \end{aligned} \]

为了优雅（方便）起见，定义 \(p_{i,0}=p_{i,i+1}=0\)

\[\begin{aligned} &S(n,n-k)\\ =&\sum_{j=1}^{k}(k+j-1)\cdot p_{k-1,j-1}\binom{n}{k+j}+\sum_{j=1}^{k}j\cdot p_{k-1,j}\binom{n}{k+j}\\ =&\sum_{j=1}^{k}\left(j\cdot p_{k-1,j}+(k+j-1)\cdot p_{k-1,j-1}\right)\binom{n}{k+j}\\ \end{aligned} \]

和我们猜想的结论对比一下

\[S(n,n-k)=\sum_{i=1}^kp_{k,i}\binom{n}{k+i} \]

递推式就得到了（讲真，这个式子很优美）

\[p_{k,i}=i\cdot p_{k-1,i}+(k+i-1)\cdot p_{k-1,i-1} \]

上面的数归也就非常完美的证明的这个猜想。

因为 \(S(n,n-1)=\binom{n}{2}\)，所以 \(p_{1,1}=1\)，加之定义的 \(p_{i,0}=p_{i,i+1}=0\)，所有的递推边界条件就齐了。

用卢卡斯定理+中国剩余定理计算组合数取模，这道题就结束了（非常毒瘤）。

代码

#include <iostream>
#include <vector>
using namespace std;
typedef long long Lint;
constexpr Lint fpow(Lint a, Lint b, Lint mod) {
    a %= mod;
    Lint res = 1;
    for (; b; b >>= 1) {
        if (b & 1)
            res = res * a % mod;
        a = a * a % mod;
    }
    return res;
}
constexpr Lint mod = 857 * 997 * 1009;
constexpr Lint p1 = 857;
constexpr Lint p2 = 997;
constexpr Lint p3 = 1009;
constexpr Lint mod_p1 = mod / p1;
constexpr Lint mod_p2 = mod / p2;
constexpr Lint mod_p3 = mod / p3;
constexpr Lint inv_mod_p1 = fpow(mod_p1, p1 - 2, p1);
constexpr Lint inv_mod_p2 = fpow(mod_p2, p2 - 2, p2);
constexpr Lint inv_mod_p3 = fpow(mod_p3, p3 - 2, p3);
constexpr Lint c1 = mod_p1 * inv_mod_p1 % mod;
constexpr Lint c2 = mod_p2 * inv_mod_p2 % mod;
constexpr Lint c3 = mod_p3 * inv_mod_p3 % mod;
Lint CRT(Lint a1, Lint a2, Lint a3) {
    return (a1 * c1 % mod + a2 * c2 % mod + a3 * c3 % mod) % mod;
}
Lint fac_p1[p1], fac_p2[p2], fac_p3[p3];
Lint inv_fac_p1[p1], inv_fac_p2[p2], inv_fac_p3[p3];
Lint bell_p1[p1 + 1], bell_p2[p2 + 1], bell_p3[p3 + 1];
Lint temp[2][2000];
Lint C_p(Lint n, Lint k, Lint p) {
    if (n < k || n < 0 || k < 0)
        return 0;
    if (p == p1) {
        return fac_p1[n] * inv_fac_p1[k] % p * inv_fac_p1[n - k] % p;
    } else if (p == p2) {
        return fac_p2[n] * inv_fac_p2[k] % p * inv_fac_p2[n - k] % p;
    } else {
        return fac_p3[n] * inv_fac_p3[k] % p * inv_fac_p3[n - k] % p;
    }
}
Lint Lucas(Lint n, Lint k, Lint p) {
    if (n < k || n < 0 || k < 0)
        return 0;
    return n == 0 ? 1 : Lucas(n / p, k / p, p) * C_p(n % p, k % p, p) % p;
}
Lint C(Lint n, Lint k) {
    if (n < k || n < 0 || k < 0)
        return 0;
    return CRT(Lucas(n, k, p1), Lucas(n, k, p2), Lucas(n, k, p3));
}
void init() {
    fac_p1[0] = fac_p2[0] = fac_p3[0] = 1;
    for (int i = 1; i < p1; i++) {
        fac_p1[i] = fac_p1[i - 1] * i % p1;
    }
    for (int i = 1; i < p2; i++) {
        fac_p2[i] = fac_p2[i - 1] * i % p2;
    }
    for (int i = 1; i < p3; i++) {
        fac_p3[i] = fac_p3[i - 1] * i % p3;
    }
    inv_fac_p1[p1 - 1] = fpow(fac_p1[p1 - 1], p1 - 2, p1);
    inv_fac_p2[p2 - 1] = fpow(fac_p2[p2 - 1], p2 - 2, p2);
    inv_fac_p3[p3 - 1] = fpow(fac_p3[p3 - 1], p3 - 2, p3);
    for (int i = p1 - 2; i >= 0; i--) {
        inv_fac_p1[i] = inv_fac_p1[i + 1] * (i + 1) % mod;
    }
    for (int i = p2 - 2; i >= 0; i--) {
        inv_fac_p2[i] = inv_fac_p2[i + 1] * (i + 1) % mod;
    }
    for (int i = p3 - 2; i >= 0; i--) {
        inv_fac_p3[i] = inv_fac_p3[i + 1] * (i + 1) % mod;
    }

    bell_p1[0] = bell_p2[0] = bell_p3[0] = 1;
    temp[0][0] = 1;
    for (int i = 1; i <= p1; i++) {
        int u = i & 1;
        bell_p1[i] = temp[u][0] = temp[u ^ 1][i - 1];
        for (int j = 1; j <= i; j++) {
            temp[u][j] = (temp[u][j - 1] + temp[u ^ 1][j - 1]) % p1;
        }
    }
    temp[0][0] = 1;
    for (int i = 1; i <= p2; i++) {
        int u = i & 1;
        bell_p2[i] = temp[u][0] = temp[u ^ 1][i - 1];
        for (int j = 1; j <= i; j++) {
            temp[u][j] = (temp[u][j - 1] + temp[u ^ 1][j - 1]) % p2;
        }
    }
    temp[0][0] = 1;
    for (int i = 1; i <= p3; i++) {
        int u = i & 1;
        bell_p3[i] = temp[u][0] = temp[u ^ 1][i - 1];
        for (int j = 1; j <= i; j++) {
            temp[u][j] = (temp[u][j - 1] + temp[u ^ 1][j - 1]) % p3;
        }
    }
}
Lint get_Bell_p(Lint n, Lint p) {
    if (p == p1)
        return bell_p1[n];
    else if (p == p2)
        return bell_p2[n];
    else
        return bell_p3[n];
}

Lint Bell_p(Lint n, Lint p) {
    vector<Lint> a;
    for (; n; n /= p) {
        a.push_back(n % p);
    }
    for (int i = 0; i <= p; i++) {
        temp[0][i] = get_Bell_p(i, p);
    }
    for (int i = 1; i < a.size(); i++) {
        for (int j = 1; j <= a[i]; j++) {
            for (int k = 0; k < p; k++) {
                temp[1][k] = (i * temp[0][k] % p + temp[0][k + 1]) % p;
            }
            temp[1][p] = (temp[1][0] + temp[1][1]) % p;
            for (int k = 0; k <= p; k++)
                temp[0][k] = temp[1][k];
        }
    }
    return temp[0][a[0]];
}

Lint Bell(Lint n) {
    return CRT(Bell_p(n, p1), Bell_p(n, p2), Bell_p(n, p3));
}

Lint p[5005][5005];
int main() {
    init();
    Lint n, k;
    cin >> n >> k;
    if (n >= k) {
        cout << Bell(k) << '\n';
        return 0;
    }
    Lint sum = Bell(k);
    p[1][1] = 1;
    for (int i = 2; i <= k - n; i++) {
        for (int j = 1; j <= i; j++) {
            p[i][j] = (i + j - 1) * p[i - 1][j - 1] % mod + j * p[i - 1][j] % mod;
            if (p[i][j] >= mod)
                p[i][j] -= mod;
        }
    }
    Lint rm = 1;
    for (int i = 1; i < k - n; i++) {
        for (int j = 1; j <= i; j++) {
            rm += p[i][j] * C(k, i + j) % mod;
            if (rm >= mod)
                rm -= mod;
        }
    }
    sum = (sum - rm + mod) % mod;
    cout << sum << '\n';
    return 0;
}

标签：p1,sum,Lint,多校,Ice,fac,binom,第三场,mod
来源： https://www.cnblogs.com/Bamboo-Wind/p/16523152.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。