有nn个灯，状态为aia_i，按ii后ii的约数都改变状态，最后要使得所有灯都关上。

如果当前最优策略≤k\le k直接用最优策略，否则随机按一个灯。

求期望步数⋅n! mod 100003\cdot n! ~ \mathrm{mod}~100003

n≤100000n\le 100000

Luogu P3750

Solution

不难发现，最优策略就是贪心地从编号大的往编号小的走，碰到一个开着的灯就灭掉，正确性很显然

那么接下来考虑>k>k的时候如何dp

直接考虑设dp[i]dp[i]表示最优策略下还有ii步结束的期望步数

dp[i]={i(i<=k)1+indp[i−1]+n−indp[i+1](k<i≤n) dp[i]= \begin{cases} i &(i <= k)\\ 1 + \frac{i}{n}dp[i-1]+\frac{n-i}{n}dp[i+1]&(k<i\le n) \end{cases}

这是因为最优策略实际上与操作顺序无关，那么就有in\frac{i}{n}的概率改对，使得最小步数-1，n−in\frac{n-i}{n}的概率改错，使得最小步数+1

但是这个dp当k=0k=0的时候就有问题了

考虑改变dpdp数组的定义，维护差分，即dp[i]dp[i]表示最优策略下从ii步到i−1i-1步的期望步数

dp[i]={1(i<=k)in+n−in(dp[i+1]+dp[i]+1)(k<i≤n) dp[i]= \begin{cases} 1 &(i <= k)\\ \frac{i}{n}+\frac{n-i}{n}(dp[i+1]+dp[i]+1)&(k<i\le n) \end{cases}

化简一下下面这个式子就能直接算了

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#include <bits/stdc++.h>

#define x first
#define y second
#define y1 Y1
#define y2 Y2
#define mp make_pair
#define pb push_back

using namespace std;

typedef long long LL;
typedef pair<int, int> pii;

template <typename T> inline int Chkmax (T &a, T b) { return a < b ? a = b, 1 : 0; }
template <typename T> inline int Chkmin (T &a, T b) { return a > b ? a = b, 1 : 0; }

inline void proc_status()
{
	ifstream t ("/proc/self/status");
	cerr << string (istreambuf_iterator <char> (t), istreambuf_iterator <char> ()) <<endl;
}

template <typename T> T read ()
{
	T sum = 0, fl = 1; char ch = getchar();
	for (; !isdigit(ch); ch = getchar()) if (ch == '-') fl = -1;
	for (; isdigit(ch); ch = getchar()) sum = (sum << 3) + (sum << 1) + ch - '0';
	return sum * fl;
}

const int Maxn = 100000 + 100, Mod = 1e5 + 3;

int N, K, A[Maxn], Dp[Maxn];

inline int Pow (int a, int b)
{
	int ans = 1;
	while (b) { if (b & 1) ans = 1ll * ans * a % Mod; a = 1ll * a * a % Mod; b >>= 1; }
	return ans;
}

inline void Solve ()
{
	int cnt = 0;
	for (int i = N; i >= 1; --i)
	{
		if (!A[i]) continue;
		++cnt;
		for (int j = 1; j * j <= i; ++j)
		{
			if (i % j) continue;
			A[j] ^= 1;
			if (j * j != i) A[i / j] ^= 1;
		}
	}
	if (cnt <= K)
	{
		for (int i = 2; i <= N; ++i) cnt = 1ll * cnt * i % Mod;
		cout<<cnt<<endl;
		return ;
	}

	Dp[N] = 1;
	for (int i = N - 1; i >= 1; --i) Dp[i] = (1 + 1ll * (N - i) * (1 + Dp[i + 1]) % Mod * Pow(i, Mod - 2) % Mod) % Mod;

	int ans = K;
	for (int i = cnt; i > K; --i) (ans += Dp[i]) %= Mod;
	for (int i = 2; i <= N; ++i) ans = 1ll * ans * i % Mod;
	cout<<ans<<endl;
}

inline void Input ()
{
	N = read<int>(), K = read<int>();
	for (int i = 1; i <= N; ++i) A[i] = read<int>();
}

int main()
{
#ifndef ONLINE_JUDGE
	freopen("A.in", "r", stdin);
	freopen("A.out", "w", stdout);
#endif
	Input();
	Solve();
	return 0;
}