[BZOJ4872][Shoi2017]分手是祝愿（期望dp）

期望太强大了！

Source

https://www.lydsy.com/JudgeOnline/problem.php?id=4872

Solution

先考虑求出操作次数最少（也就是每个开关最多操作一次）的方案。
考虑到操作第 $i$ 个开关，会使状态改变的灯一定为 $i$ 或在 $i$ 的前面。
因此，可以从第 $n$ 个开关倒序到第 $1$ 个开关进行讨论，当走到第 $i$ 个开关时，如果第 $i$ 个灯是亮的，那么就对第 $i$ 个开关进行操作并改变编号为 $i$ 的约数的灯的状态。
设最少需要 $cnt$ 步才能使所有灯都灭掉。
当 $cnt\le k$ 时，期望操作次数为 $cnt$ 。
否则将问题抽象成： $n$ 个数，有 $cnt$ 个数为 $1$ ，其余为 $0$ 。每次会等概率随机地选择一个数进行取反，如果在某时刻 $1$ 的个数恰好为 $k$ 则不再随机，直接将这 $k$ 个 $1$ 取反。求将所有数变成 $0$ 的期望取反次数。
DP 。定义状态： $f[i]$ 表示当前有 $i$ 个 $1$ ，将 $1$ 的个数减少 $1$ （即变为 $i-1$ ）的期望步数。
根据期望的线性性质，可以从 $f$ 得出将所有数变为 $0$ 的期望步数：

$$k+\sum_{i=k+1}^{cnt}f[i]$$
$f$ 的边界为 $f[n]=1$ 。
而对于一个 $k<i<n$ ，
考虑第一次等概率随机选出的数是 $1$ 或 $0$ ：
选出的数是 $1$ 时，这个 $1$ 变成 $0$ 。 $n$ 个数中 $1$ 的数量减少 $1$ ，已经达到目的。概率为 $\frac in$ ，步数为 $1$ 。
选出的数是 $0$ 时，这个 $0$ 变成 $1$ 。 $n$ 个数中 $1$ 的数量增加 $1$ ，这时候则需要将 $1$ 的数量减少 $1$ 后再减少 $1$ 。概率为 $\frac{n-i}n$ ，步数为 $1+f[i+1]+f[i]$ 。
于是得出 $f[i]$ 与 $f[i+1]$ 的关系式：
$$f[i]=\frac in\times 1+\frac{n-i}n\times(1+f[i+1]+f[i])$$
这还不能达到递推的目的，因此移项，将 $f[i]$ 移到左边：
$$f[i]=1+\frac{n-i}i\times(1+f[i+1])$$
得到 $f$ 的递推式。
注意最后乘上 $n!$ 。

Code

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define For(i, a, b) for (i = a; i <= b; i++)
#define Rof(i, a, b) for (i = a; i >= b; i--)
using namespace std;
inline int read() {
    int res = 0; bool bo = 0; char c;
    while (((c = getchar()) < '0' || c > '9') && c != '-');
    if (c == '-') bo = 1; else res = c - 48;
    while ((c = getchar()) >= '0' && c <= '9')
        res = (res << 3) + (res << 1) + (c - 48);
    return bo ? ~res + 1 : res;
}
const int N = 1e5 + 5, PYZ = 1e5 + 3;
int n, k, a[N], cnt, np = 1, f[N];
void orzpyzdalao(int x) {
    int i, S = sqrt(x); For (i, 1, S) {
        if (x % i) continue; a[i] ^= 1; if (i * i != x) a[x / i] ^= 1;
    }
}
int qpow(int a, int b) {
    int res = 1; while (b) b & 1 ? res = 1ll * res * a % PYZ : 0,
        a = 1ll * a * a % PYZ, b >>= 1; return res;
}
int main() {
    int i; n = read(); k = read(); For (i, 1, n) a[i] = read();
    Rof (i, n, 1) if (a[i]) orzpyzdalao(i), cnt++;
    For (i, 1, n) np = 1ll * np * i % PYZ;
    if (cnt <= k) return printf("%d\n", 1ll * cnt * np % PYZ), 0;
    f[n] = 1; Rof (i, n - 1, k + 1)
        f[i] = (1 + 1ll * (n - i) * qpow(i, PYZ - 2) % PYZ
            * (1 + f[i + 1]) % PYZ) % PYZ;
    int sum = k; Rof (i, cnt, k + 1) sum = (sum + f[i]) % PYZ;
    cout << 1ll * sum * np % PYZ << endl; return 0;
}