题解 P3750 【[六省联考2017]分手是祝愿】

做法:模拟高斯消元

为什么没几篇模拟高斯消元的题解啊？这不是最容易想到的吗？

首先根据部分分提示，我们可以先想怎样操作次数最少。比较明显的是，我们每次选择最后一个亮着的灯是最优的，因为选择第 (i) 号灯能影响到的只有小于等于 (i) 的灯，而不影响后面的灯。因此，如果不去操作最后一个亮着的灯的话，以后去前面操作不会再让其灭掉。

依照这种方法，我们可以得到一个最优的操作序列，并且可以对这个序列以任意顺序进行操作。

然后再考虑随机操作。我们发现随便对一个灯进行操作，要么恰好符合我们的期望，提前帮我们搞定一步；要么操作了其它的点，我们的最少操作次数又多了一。

根据套路，我们设 (f[i]) 表示“最少需要 (i) 步”的状态，到结束的期望次数，并且有式子:

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

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

发现无法递推，不过高斯消元足以通过前一半的数据。我们考虑对高斯消元进行优化。我们发现矩阵大概长这个样子：

[egin{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
", ((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;
}