其他分享
首页 > 其他分享> > 2022牛客暑期多校第三场 I. Ice Drinking

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

作者:互联网

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