[六省联考2017]分手是祝愿

观察数据范围可以发现，题目有一半的分是满足 (k = n) 的性质的，这肯定是有用的，我们可以从这里开始下手。

因为特殊条件是 (k = n)，我们可以猜测每个状态的操作最小次数应该是不大于 (n) 的，近一步地我们可以发现使得操作次数最小的操作方案。因为标号小的操作不会影响到编号大的，而编号大的必须要让它变为 (0)，因此我们从编号最大的开始，如果当前状态为 (1) 我们就在这个位置操作一次，这样显然是最优的，因为我们只能通过操作更大的位置使得这个位置变为 (0)，假设我们现在考虑的是最大的为 (1) 的位置，而如果我们操作了一个更大的为 (0) 的位置，那么这个位置和操作位置中间的 (1) 我们还需要操作成 (0)，相当于没有操作，因此这样是不优的，所以上面那个贪心是正确的。

由上面的那个贪心证明，可以发现操作方案是唯一的。并且我们可以发现这样一件事，我们的操作方案相当于是在指定的一些位置上进行操作，假设我们在这些指定位置上操作一次，那么最优操作次数肯定会 (-1)，相反，如果在不是指定的位置上操作一次，那么最优操作次数一定会 (+1)（因为操作方案唯一，在不是指定位置上操作最优操作次数一定会增加，而我们再在这个位置上操作一次就可以让序列变回来）。因此我们当前的最优操作次数的改变实际上于原序列无关，而至于当前的最优操作次数有关。于是我们可以得到一个 (dp)，令 (dp_i) 表示当前的最优操作次数为 (i) 次使得最优操作次数达到 (k) 还需要的期望次数。那么我们有转移 (dp_{i} = frac{i}{n}dp_{i - 1} + frac{n - i}{n}dp_{i + 1} + 1)，我们当然可以暴力上高斯消元来搞这个东西。但时间复杂度并不允许。

可以发现上面的 (dp) 难以转移是因为每次我们变量之间组成了一个方程组，必须使用高斯消元，回忆我们之前做过的期望 (dp)，可以发现都是向同一个方向之间转移，这样我们可以通过简单的移项来获得变量之间的关系。那么我们有没有办法让这个 (dp) 一样地只往一个方向转移呢？不难发现我们应该改写 (dp) 状态以达到这个目的，一个经典的 (trick) 就是我们可以考虑将状态改写成从当前状态到原转移方程另一个状态的期望值，即我们令 (dp_i) 为最优次数从 (i) 转移到 (i - 1) 的期望次数。于是我们就有转移 (dp_i = frac{i}{n} imes 1 + frac{n - i}{n}(dp_{i + 1} + dp_i + 1))，表示我们有 (i) 种选择花费一步转移到 (i - 1)，有 (n - i) 种选择花费 (1) 先转移到 (i + 1)，花费 (dp_{i + 1}) 转移回来再花费 (dp_i) 转移到 (i - 1)。于是通过简单的移项我们可以得到 (dp_i = frac{n - i}{n}dp_{i + 1} + frac{n}{i})，根据期望的线性性可以知道 (sumlimits_{i = k + 1} ^ {limit} dp_i) 即为最终答案，(limit) 为原序列的最优操作次数。

#include<bits/stdc++.h>
using namespace std;
#define N 100000 + 5
#define Mod 100003
#define rep(i, l, r) for(int i = l; i <= r; ++i)
#define dep(i, l, r) for(int i = r; i >= l; --i)
int n, k, fac, ans, limit, a[N], dp[N];
int read(){
    char c; int x = 0, f = 1;
    c = getchar();
    while(c > '9' || c < '0'){ if(c == '-') f = -1; c = getchar();}
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x * f;
}
int Inc(int a, int b){
    return (a += b) >= Mod ? a - Mod : a;
}
int Mul(int a, int b){
    return 1ll * a * b % Mod;
}
int Qpow(int a, int b){
    int ans = 1;
    while(b){
        if(b & 1) ans = Mul(ans, a);
        a = Mul(a, a), b >>= 1;
    }
    return ans;
}
int Div(int a, int b){
    return Mul(a, Qpow(b, Mod - 2));
}
int main(){
    n = read(), k = read(), fac = 1;
    rep(i, 1, n) a[i] = read(), fac = 1ll * fac * i % Mod;
    dep(i, 1, n){
        int cnt = 0;
        for(int j = i; j <= n; j += i) if(a[j]) ++cnt;
        if(cnt & 1) a[i] = 1, ++limit;
        else a[i] = 0;
    }
    if(limit <= k) printf("%d", Mul(limit, fac));
    else{
        dp[n] = 1, ans = k;
        dep(i, k + 1, n - 1) dp[i] = Inc(Mul(Div(n - i, i), dp[i + 1]), Div(n, i));
        rep(i, k + 1, limit) ans = Inc(ans, dp[i]);
        printf("%d", Mul(ans, fac));
    }
    return 0;
}

实际上这道题还有另一种方法，我们继续考虑上面的那个转移方程 (dp_{i} = frac{i}{n}dp_{i - 1} + frac{n - i}{n}dp_{i + 1} + 1) 可知 (dp_n = dp_{n - 1} + 1)，然后有 (dp_{n - 1} = frac{n - 1}{n}dp_{i - 2} + frac{1}{n}dp_{n} + 1)，将上面那个柿子带进来有，(dp_{n - 1} = frac{n - 1}{n}dp_{n - 2} + frac{1}{n}(dp_{n - 1} + 1) + 1)，移项有 (dp_{n - 1} = dp_{n - 2} + frac{n + 1}{n - 1})，因为我们得到的是 (f_x = f_{x - 1} + b) 的形式，因此将这个带到下一个柿子还会是这样的形式，因为移项后可以将系数刚好约掉，因此我们可以令 (dp_i = dp_{i - 1} + b_i)，带入最开始的转移方程有 (dp_i = frac{i}{n}(dp_i - b_i) + frac{n - i}{n}(dp_i + k_{i + 1}))，移项可得 (k_i = frac{n - i}{i}k_{i + 1})，于是我们先递推出 (k)，后递推出 (dp) 即可。

这里告诉我们在转移比较复杂的情况下我们可以从特殊的转移开始考虑，或许可以发现其中的规律。

实际上这两种方法是等价的，上面那种方法是建立在只用两种转移方向的前提下，而因为我们一个点向外转移条件概率和必为 (1)，因此第二种方法能刚好约掉是必然的。

#include<bits/stdc++.h>
using namespace std;
#define N 100000 + 5
#define Mod 100003
#define rep(i, l, r) for(int i = l; i <= r; ++i)
#define dep(i, l, r) for(int i = r; i >= l; --i)
int n, k, fac, limit, ans, a[N], dp[N], b[N];
int read(){
    char c; int x = 0, f = 1;
    c = getchar();
    while(c > '9' || c < '0'){ if(c == '-') f = -1; c = getchar();}
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x * f;
}
int Inc(int a, int b){
    return (a += b) >= Mod ? a - Mod : a;
}
int Mul(int a, int b){
    return 1ll * a * b % Mod;
}
int Qpow(int a, int b){
    int ans = 1;
    while(b){
        if(b & 1) ans = Mul(ans, a);
        a = Mul(a, a), b >>= 1;
    }
    return ans;
}
int Div(int a, int b){
    return Mul(a, Qpow(b, Mod - 2));
}
int main(){
    n = read(), k = read(), fac = 1;
    rep(i, 1, n) a[i] = read(), fac = 1ll * fac * i % Mod;
    dep(i, 1, n){
        int cnt = 0;
        for(int j = i; j <= n; j += i) if(a[j]) ++cnt;
        if(cnt & 1) a[i] = 1, ++limit;
        else a[i] = 0;
    }
    if(limit <= k) printf("%d", Mul(limit, fac));
    else{
        b[n] = 1;
        dep(i, k, n - 1) b[i] = Inc(Mul(Div(n - i, i), b[i + 1]), Div(n, i));
        dp[k] = k;
        rep(i, k + 1, limit) dp[i] = Inc(dp[i - 1], b[i]);
        printf("%d", Mul(dp[limit], fac));
    }
    return 0;
}