[SHOI2017] 分手是祝愿

Description

Zeit und Raum trennen dich und mich.
时空将你我分开。B 君在玩一个游戏，这个游戏由 n 个灯和 n 个开关组成，给定这 n 个灯的初始状态，下标为
从 1 到 n 的正整数。每个灯有两个状态亮和灭，我们用 1 来表示这个灯是亮的，用 0 表示这个灯是灭的，游戏
的目标是使所有灯都灭掉。但是当操作第 i 个开关时，所有编号为 i 的约数（包括 1 和 i）的灯的状态都会被
改变，即从亮变成灭，或者是从灭变成亮。B 君发现这个游戏很难，于是想到了这样的一个策略，每次等概率随机
操作一个开关，直到所有灯都灭掉。这个策略需要的操作次数很多， B 君想到这样的一个优化。如果当前局面，
可以通过操作小于等于 k 个开关使所有灯都灭掉，那么他将不再随机，直接选择操作次数最小的操作方法（这个
策略显然小于等于 k 步）操作这些开关。B 君想知道按照这个策略（也就是先随机操作，最后小于等于 k 步，使
用操作次数最小的操作方法）的操作次数的期望。这个期望可能很大，但是 B 君发现这个期望乘以 n 的阶乘一定
是整数，所以他只需要知道这个整数对 100003 取模之后的结果。

Input

第一行两个整数 n, k。
接下来一行 n 个整数，每个整数是 0 或者 1，其中第 i 个整数表示第 i 个灯的初始情况。
1 ≤ n ≤ 100000, 0 ≤ k ≤ n；

Output

输出一行，为操作次数的期望乘以 n 的阶乘对 100003 取模之后的结果。

Sample Input

4 0
0 0 1 1

Sample Output

512

Solution

看完题面就觉得这题非常的丧。
先考虑不随机按按钮的情况。
我们发现在最优情况下，我们按编号从大到小的顺序关灯，这样的步数是最小的，因为编号小的灯的按钮开关不会影响编号大的灯的开关，这样最坏情况下我们只要按 nnn 次就可以关闭所有灯。
然后我们发现， 50%" role="presentation" style="position: relative;">50%50%50\% 的数据 k=nk=nk = n，这不是送分吗，先打了再说。
就按编号从大到小 forforfor 一遍，如果这盏灯是开着的把所有编号为它的因数的灯都按一遍，最后统计就可以了。
如果我们求出的最优解 StepStepStep 大于 kkk 该怎么办呢。
我们设 f[i]" role="presentation" style="position: relative;">f[i]f[i]f[i] 为我们需要按 iii 次按钮时能结束的期望步数。
易得我们有 in" role="presentation" style="position: relative;">inin\frac{i}{n} 的概率按到我们需要的按钮 n−inn−in\frac{n - i}{n} 的概率按到我们不需要的按钮。
于是 f[i]=in×f[i−1]+n−in×f[i+1]+1f[i]=in×f[i−1]+n−in×f[i+1]+1f[i] = \frac{i}{n} \times f[i - 1] + \frac{n - i}{n} \times f[i + 1] + 1.
然而我们发现，我们最多有 n−1n−1n - 1 个方程而要求解 nnn 个变量，这显然是不可行的。
所以我们把 f" role="presentation" style="position: relative;">fff 数组差分一下，设 g[i]=f[i]−f[i−1]g[i]=f[i]−f[i−1]g[i] = f[i] - f[i - 1]。
整理得 g[i]=n+(n−i)×f[i+1]ig[i]=n+(n−i)×f[i+1]ig[i] = \frac{n + (n - i) \times f[i + 1]}{i}
最后答案为 k+∑Stepi=k+1g[i]k+∑i=k+1Stepg[i]k + \sum_{i = k + 1}^{Step} g[i]

Code

#include <algorithm>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <iostream>
#include <queue>
#include <set>
#include <stack>#define R register
#define ll long long
#define db double
#define sqr(_x) (_x) * (_x)
#define Cmax(_a, _b) ((_a) < (_b) ? (_a) = (_b), 1 : 0)
#define Cmin(_a, _b) ((_a) > (_b) ? (_a) = (_b), 1 : 0)
#define Max(_a, _b) ((_a) > (_b) ? (_a) : (_b))
#define Min(_a, _b) ((_a) < (_b) ? (_a) : (_b))
#define Abs(_x) (_x < 0 ? (-(_x)) : (_x))using namespace std;namespace Dntcry
{inline int read(){R int a = 0, b = 1; R char c = getchar();for(; c < '0' || c > '9'; c = getchar()) (c == '-') ? b = -1 : 0;for(; c >= '0' && c <= '9'; c = getchar()) a = (a << 1) + (a << 3) + c - '0';return a * b;}inline ll lread(){R ll a = 0, b = 1; R char c = getchar();for(; c < '0' || c > '9'; c = getchar()) (c == '-') ? b = -1 : 0;for(; c >= '0' && c <= '9'; c = getchar()) a = (a << 1) + (a << 3) + c - '0';return a * b;}const int Maxn = 100010, Mod = 100003;int n, k;int v[Maxn];ll f[Maxn], test, Ans, fac, inv[Maxn];int Main(){n = read(), k = read();fac = inv[0] = inv[1] = 1; for(R int i = 1; i <= n; i++) v[i] = read(), fac = 1ll * fac * i % Mod;for(R int i = n; i > 0; i--) if(v[i]) {R int lim = (int)sqrt(i);++test;for(R int j = 1; j <= lim; j++) if(i % j == 0) {v[j] ^= 1;if(j * j != i) v[i / j] ^= 1;}}if(test <= k) return !printf("%lld\n", test * fac % Mod);for(R int i = 2; i <= n; i++) inv[i] = (Mod - Mod / i) * inv[Mod % i] % Mod;f[n] = 1;for(R int i = n - 1; i >= k; i--) f[i] = 1ll * (n + 1ll * (n - i) * f[i + 1] % Mod) * inv[i] % Mod;for(R int i = test; i > k; i--) Ans = (Ans + f[i]) % Mod;Ans += k;printf("%lld\n", Ans * fac % Mod);return 0;}
}
int main()
{return Dntcry :: Main();
}