做法:模擬高斯消元

為什麼沒幾篇模擬高斯消元的題解啊？這不是最容易想到的嗎？

首先根據部分分提示，我們可以先想怎樣操作次數最少。比較明顯的是，我們每次選擇最後一個亮著的燈是最優的，因為選擇第 \(i\) 號燈能影響到的只有小於等於 \(i\) 的燈，而不影響後面的燈。因此，如果不去操作最後一個亮著的燈的話，以後去前面操作不會再讓其滅掉。

依照這種方法，我們可以得到一個最優的操作序列，並且可以對這個序列以任意順序進行操作。

然後再考慮隨機操作。我們發現隨便對一個燈進行操作，要麼恰好符合我們的期望，提前幫我們搞定一步；要麼操作了其它的點，我們的最少操作次數又多了一。

根據套路，我們設 \(f[i]\)

\[f[t] = t, t <= k \]

\[f[i] = 1 + \frac{i}{n}f[i-1] + \frac{n-i}{n}f[i+1],t>k \]

發現無法遞推，不過高斯消元足以通過前一半的資料。我們考慮對高斯消元進行優化。我們發現矩陣大概長這個樣子：

\[\begin{bmatrix} 1 & & & & & & & & & & 1\\ & 1 & & & & & & & & & 2\\ & & ... & & & & & & & & \\ & & & 1 & & & & & & & k\\ & & & -\frac{k+1}{n} & 1 & -\frac{n-k-1}{n} & & & & & 1\\ & & & & -\frac{k+2}{n} & 1 & -\frac{n-k-2}{n} & & & & 1\\ & & & & & -\frac{k+3}{n} & 1 & -\frac{n-k-3}{n} & & & 1\\ & & & & & & ... & ... & ... & & \\ & & & & & & & -\frac{n-1}{n} & 1 & -\frac{1}{n} & 1\\ & & & & & & & & -\frac{n}{n} & 1 &1 \end{bmatrix}\]

只有一條對角線，自然可以 \(O(n)\) 高斯消元求（避過那些那零減來減去的步驟就好)。不過空間開不下 \(n^2\) 陣列，需要手動模擬（詳見程式碼）。

複雜度：\(O(nlogn)\)（瓶頸在求逆元）

不過，這題比較噁心的一點是，模數只有 \(1e5\) 大小，很容易加加減減出現 \(0\)，使得我們無法求逆元，這時需要交換相鄰兩行，對其進行特殊處理（詳見程式碼）。如果模數開到 \(1e9\) 那麼大就基本沒有這個問題了。

\(Code:\)

#define N 101000
typedef long long ll;
template <typename T> inline void read(T &x) {
	x = 0; char c = getchar(); bool flag = false;
	while (!isdigit(c)) { if (c == '-')	flag = true; c = getchar(); }
	while (isdigit(c)) { x = (x << 1) + (x << 3) + (c ^ 48); c = getchar(); }
	if (flag)	x = -x;
}
using namespace std;
const int P = 100003;
inline ll quickpow(ll x, int k) {
	ll res = 1;
	while (k) {
		if (k & 1)	res = res * x % P;
		x = x * x % P;
		k >>= 1;
	}
	return res;
}
int n, K;
int state[N];
ll f[N], h[N][4];
bool tag[N];
inline void work() {
	for (register int i = 1; i <= K; ++i)	f[i] = i;
	h[K][1] = 1, h[K][3] = K;
	for (register int i = K + 1; i <= n; ++i) {
		h[i][0] = -i, h[i][1] = n, h[i][2] = i - n, h[i][3] = n;
	}
	for (register int i = K + 1; i <= n; ++i) {
		ll t = h[i][0];
		h[i][0] = 0; h[i][1] = (h[i][1] - t * h[i - 1][2]) % P;
		h[i][3] = (h[i][3] - t * h[i - 1][3]) % P;
		
		if (h[i][1]) {
			t = quickpow(h[i][1], P - 2);
			h[i][2] = h[i][2] * t % P; h[i][3] = h[i][3] * t % P;
			h[i][1] = 1;
		} else {//不存在逆元
			t = quickpow(h[i][2], P - 2);
			h[i][2] = 1, h[i][3] = h[i][3] * t % P;
			t = h[i + 1][1];
			h[i + 1][1] = 0, h[i + 1][3] = (h[i + 1][3] - t * h[i][3]) % P;//消去下面一行的中間的那個數
			tag[i + 1] = true;
			swap(h[i][2], h[i][1]);//“偏移”陣列
			swap(h[i], h[i + 1]);//交換兩行
		}
	}
	for (register int i = n - 1; i > K; --i) {
		if (!tag[i]) {
			ll t = h[i][2];
			h[i][2] = 0; h[i][3] = (h[i][3] - t * h[i + 1][3]) % P;
		} else {//對無逆元情況的特殊處理
			ll t = h[i - 1][2];
			h[i - 1][2] = 0, h[i - 1][3] = (h[i - 1][3] - t * h[i + 1][3]) % P;
			t = quickpow(h[i - 1][0], P - 2);
			h[i - 1][3] = h[i - 1][3] * t % P;
			swap(h[i - 1][1], h[i - 1][0]);
			--i;
		}
	}
	for (register int i = K + 1; i <= n; ++i)	f[i] = h[i][3];
}


vector<int> vec[N];
inline void get_ans() {
	for (register int i = 1; i <= n; ++i)
		for (register int j = i; j <= n; j += i)
			vec[j].push_back(i);//計算操作序列
	int nw = 0;
	for (register int i = n; i; --i) {
		if (state[i]) {
			for (register unsigned int j = 0; j < vec[i].size(); ++j) {
				state[vec[i][j]] ^= 1;//計算狀態的最少操作次數
			}
			++nw;
		}
	}
	ll res = 1;
	for (register int i = 1; i <= n; ++i)	res = res * i % P;
	printf("%lld\n", ((res * f[nw] % P) + P) % P);
}

int main() {
	read(n), read(K);
	for (register int i = 1; i <= n; ++i)	read(state[i]);
	work();
	get_ans();
	return 0;
}