[SHOI2017]分手是祝愿

链接：https://ac.nowcoder.com/acm/problem/20437
来源：牛客网
题目：

B 君在玩一个游戏，这个游戏由 n 个灯和 n 个开关组成，给定这 n 个灯的初始状态，下标为从 1 到 n 的正整数。每个灯有两个状态亮和灭，我们用 1 来表示这个灯是亮的，用 0 表示这个灯是灭的，游戏 的目标是使所有灯都灭掉。

但是当操作第 i 个开关时，所有编号为 i 的约数（包括 1 和 i）的灯的状态都会被改变，即从亮变成灭，或者是从灭变成亮。B 君发现这个游戏很难，于是想到了这样的一个策略，每次等概率随机操作一个开关，直到所有灯都灭掉。这个策略需要的操作次数很多， B 君想到这样的一个优化。如果当前局面， 可以通过操作小于等于 k 个开关使所有灯都灭掉，那么他将不再随机，直接选择操作次数最小的操作方法（这个 策略显然小于等于 k 步）操作这些开关。

B 君想知道按照这个策略（也就是先随机操作，最后小于等于 k 步，使用操作次数最小的操作方法）的操作次数的期望。这个期望可能很大，但是 B 君发现这个期望乘以 n 的阶乘一定 是整数，所以他只需要知道这个整数对100003 取模之后的结果。

输入描述:

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

输出描述:

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

输入：

4 0
0 0 1 1

输出：

512

题解：

首先，我们应当明确三件事情：

1.如果直接考虑用最少的操作关掉所有灯，那么从最后面的灯开始关，如果灯亮着，关掉它，并翻转它所有因数编号的灯。一直遍历到第一盏灯为止。

2.假设关掉所有灯共需cnt次操作，那么任意交换这cnt次操作中的两次，不影响最终结果。

3.最少操作是唯一确定的。

假设f[i]为在最少操作次数为i的情况下，需要操作的平均次数。

由题意知，i<=k时，f[i]=i。

当i>k时，f[i]=(i/n)*f[i-1]+((n-i)/n)*f[i+1]+1（*）。

(i/n)*f[i-1]表示如果选中了n次操作中的其中i次，由条件1知其最少操作次数-1；

((n-i)/n)*f[i+1]表示如果选中了剩下的n-i次操作,那么你就得再花一次操作来抵消此次操作，因此最少操作次数 +1。

令d[i]=f[i]-f[i-1]，（*）式可化为d[i]=n/i+(n-i)/i*d[i+1]。

易知f[n]=f[n-1]+1，因此d[n]=1，可以通过递推方式将d[i]计算完毕，进而算出i>k时的所有f[i]。

由条件1计算得最少操作为cnt，那么最终答案为n!*f[cnt]%100003。

另：除法取模需借助逆元计算。

代码：

#include <bits/stdc++.h>
#define maxn 100010
#define ll long long
using namespace std;
const ll mod=100003;
ll n,k,f[maxn],s[maxn],d[maxn],opt[maxn];


ll ni[maxn],cnt=0;
void init(){
	ll i,j;
	memset(f,0,sizeof(f));
	memset(d,0,sizeof(d));
	memset(ni,0,sizeof(ni));
	ni[1]=1;
	for(i=2;i<mod;i++)ni[i]=-ni[mod%i]*(mod/i)%mod+mod;
	
	for(i=n;i>=1;i--){
		if(s[i]==1){
			cnt++;
			for(j=1;j<=sqrt(i+0.5);j++){
				if(i%j==0&&i>j){
					ll c=i/j;
					if(c>j){
						s[j]=s[j]^1;
						s[c]=s[c]^1;
					}
					else s[j]=s[j]^1;
				}
			}
		}
	}
}


int main(){
	int i;
	scanf("%lld%lld",&n,&k);
	for(i=1;i<=n;i++)scanf("%lld",&s[i]);
	init();
	d[n]=1;
	for(i=n-1;i>=1;i--){
		d[i]=((n-i)*ni[i]+mod)%mod*d[i+1]%mod+(n*ni[i]+mod)%mod;
		d[i]=d[i]%mod;
	}
	for(i=1;i<=k;i++)f[i]=i;
	for(i=k+1;i<=n;i++){
		f[i]=f[i-1]+d[i]+mod;
		f[i]=f[i]%mod;	
	}
	ll ans=1;
	for(i=1;i<=n;i++){
		ans*=i;
		ans+=mod;
		ans%=mod;
	}
	cout<<ans*f[cnt]%mod;
	return 0;
}